(:
Port of the Java code at http://www.webcoding.co.uk/?page=UKGridToLatLon
into XQuery
    
Converts a UK Grid reference easting and northing into a decimal lattitude and longitude

Author: Adam Retter <adam dot retter at devon dot gov dot uk>

:)

declare namespace math="http://exist-db.org/xquery/math";


declare function local:UKGridToLatLong($easting as xs:double, $northing as xs:double) as xs:double+
{
	let 
		$deg2rad := math:pi() div 180,
		$rad2deg := 180.0 div math:pi(),
		$pi := math:pi(),

		$wgs_axis := 6378137,
		$wgs_eccent := 0.00669438,
		$latorig := 49 * $deg2rad,
		$lonorig := -2 * $deg2rad,
		$scale := 0.9996012717,
		$falseeast := 400000,
		$falsenorth := -100000,
		$axis := 6377563.396,
		$eccent := 0.00667054,

		$ep := $easting,
		$np := $northing,

		$p := $eccent div 8,
		$a := $axis * (1 - (2 * $p) - (3 * $p * $p) - (10 * $p * $p * $p)),
		$b := $axis * ((6 * $p) + (12 * $p * $p) + (45 * $p * $p * $p)) div 2,
		$c := $axis * ((15 * $p * $p) + (90 * $p * $p * $p)) div 4,
		$mo := $a * $latorig - $b * math:sin(2 * $latorig) + $c * math:sin(4 * $latorig),
		$mp := $mo + (($np - $falsenorth) div $scale),
		$phidash := $mp div $a,
		$phif := $phidash + (($b * math:sin(2 * $phidash)) - ($c * math:sin(4 * $phidash))) div ($a - (2 * $b * math:cos(2 * $phidash))),

		$uf := $axis div math:sqrt(1 - ($eccent * (math:sin($phif) * math:sin($phif)))),
		$h := ($ep - $falseeast) div ($scale * $uf),
		$nsqd := $eccent * (math:cos($phif) * math:cos($phif)) div (1 - $eccent),
		$tsqd := math:power(math:sin($phif) div math:cos($phif), 2),
		$lambdap := $lonorig + ( 1 div math:cos($phif)) * (($h - (($h * $h * $h) div 6) * (1 + (2 * $tsqd) + $nsqd))),
		$phip := $phif - ((1 + $nsqd) * (math:sin($phif) div math:cos($phif)) * ((($h * $h) div 2) - (($h * $h * $h * $h) div 24) * (5 + 3 * $tsqd))),

		$result1 := local:ll_to_cartesian($phip, $lambdap, $axis, $eccent, 0),

		$x := $result1[1],
		$y := $result1[2],
		$z := $result1[3],

		$result2 := local:cartesian_to_ll($x, $y, $z, $wgs_eccent, $wgs_axis),

		$lat := $result2[1],
		$lon := $result2[2]

	return
	
		($lat * $rad2deg, $lon * $rad2deg)
};


declare function local:cartesian_to_ll_lat_loop($errvalue as xs:double, $lat as xs:double, $z as xs:double, $eccentsq as xs:double, $v as xs:double, $p as xs:double) as xs:double
{
	if($errvalue > 0.001)then
	(
		let $lat0 := math:atan(($z + $eccentsq * $v * math:sin($lat)) div $p),
		$errvalue := math:abs($lat0 - $lat) return

		local:cartesian_to_ll_lat_loop($errvalue, $lat0, $z, $eccentsq, $v, $p)
	)
	else
	(
		$lat
	)
};


declare function local:cartesian_to_ll($x as xs:double, $y as xs:double, $z as xs:double, $eccentsq as xs:double, $a as xs:double) as xs:double+
{
	let 
		$lon := math:atan($y div $x),
		$p := math:sqrt(($x * $x) + ($y * $y)),
		$lat := math:atan($z div ($p * (1 - $eccentsq))),
		$v := $a div (math:sqrt(1 - $eccentsq * (math:sin($lat) * math:sin($lat)))),
		$errvalue := 1.0,

		$lat0 := 0,

		$finalLat := local:cartesian_to_ll_lat_loop($errvalue, $lat, $z, $eccentsq, $v, $p)

	return

		($finalLat, $lon)
};

declare function local:ll_to_cartesian($lat as xs:double, $lon as xs:double, $a as xs:double, $eccentsq as xs:double, $h as xs:double) as xs:double+
{
	let
		$v := $a div (math:sqrt(1 - $eccentsq * (math:power(math:sin($lat), 2)))),
		$x := ($v + $h) * math:cos($lat) * math:cos($lon),
		$y := ($v + $h) * math:cos($lat) * math:sin($lon),
		$z := ((1 - $eccentsq) * $v + $h) * math:sin($lat)
	return

		($x, $y, $z)
};

local:UKGridToLatLong(309950, 093393)