udf_OSGrid_WGS84.sql

create function dbo.udf_OSGrid_WGS84 (@east int, @north int)
returns geography
as
begin
	-- reference implementation: http://www.movable-type.co.uk/scripts/latlong-os-gridref.html

    if @east is null
        or @north is null
        return null

    declare @π float = pi()

    declare @a float = 6377563.396,
        @b float = 6356256.910,
        @f float = 1 / 299.3249646 					-- Airy 1830 major & minor semi-axes. and flattening
    declare @F0 float = 0.9996012717                -- NatGrid scale factor on central meridian
    declare @φ0 float = 49 * @π / 180,
        @λ0 float = -2 * @π / 180					-- NatGrid true origin is 49°N,2°W
    declare @N0 int = -100000,
        @E0 int = 400000							-- northing & easting of true origin, metres
    declare @e2 float = 1 - (@b * @b) / (@a * @a)	-- eccentricity squared
    declare @n float = (@a - @b) / (@a + @b)
    declare @n2 float = @n * @n
    declare @n3 float = @n * @n * @n
 
    declare @φ float = @φ0,
        @M float = 0

    while (@north - @N0 - @M >= 0.00001) -- ie until < 0.01mm
    begin
        set @φ = (@north - @N0 - @M) / (@a * @F0) + @φ
 
        declare @Ma float = (1 + @n + (5 / 4) * @n2 + (5 / 4) * @n3) * (@φ - @φ0)
        declare @Mb float = (3 * @n + 3 * @n * @n + (21 / 8) * @n3) * sin(@φ - @φ0) * cos(@φ + @φ0)
        declare @Mc float = ((15 / 8) * @n2 + (15 / 8) * @n3) * sin(2 * (@φ - @φ0)) * cos(2 * (@φ + @φ0))
        declare @Md float = (35 / 24) * @n3 * sin(3 * (@φ - @φ0)) * cos(3 * (@φ + @φ0))

		-- meridional arc
        set @M = @b * @F0 * (@Ma - @Mb + @Mc - @Md)                
    end 
 
    declare @cosφ float = cos(@φ),
        @sinφ float = sin(@φ)
    declare @nu float = @a * @F0 / sqrt(1 - @e2 * @sinφ * @sinφ)					-- nu = transverse radius of curvature
    declare @ρ float = @a * @F0 * (1 - @e2) / power(1 - @e2 * @sinφ * @sinφ, 1.5)  -- rho = meridional radius of curvature
    declare @η2 float = @nu / @ρ - 1
 
    declare @tanφ float = tan(@φ)
    declare @tan2φ float = @tanφ * @tanφ
    declare @tan4φ float = @tan2φ * @tan2φ
    declare @tan6φ float = @tan4φ * @tan2φ

    declare @secφ float = 1 / @cosφ
    declare @nu3 float = @nu * @nu * @nu
    declare @nu5 float = @nu3 * @nu * @nu
    declare @nu7 float = @nu5 * @nu * @nu

    declare @VII float = @tanφ / (2 * @ρ * @nu)
    declare @VIII float = @tanφ / (24 * @ρ * @nu3) * (5 + 3 * @tan2φ + @η2 - 9 * @tan2φ * @η2)
    declare @IX float = @tanφ / (720 * @ρ * @nu5) * (61 + 90 * @tan2φ + 45 * @tan4φ)
    declare @X float = @secφ / @nu
    declare @XI float = @secφ / (6 * @nu3) * (@nu / @ρ + 2 * @tan2φ)
    declare @XII float = @secφ / (120 * @nu5) * (5 + 28 * @tan2φ + 24 * @tan4φ)
    declare @XIIA float = @secφ / (5040 * @nu7) * (61 + 662 * @tan2φ + 1320 * @tan4φ + 720 * @tan6φ)
 
    declare @dE float = (@east - @E0)
    declare @dE2 float = @dE * @dE
    declare @dE3 float = @dE2 * @dE
    declare @dE4 float = @dE2 * @dE2
    declare @dE5 float = @dE3 * @dE2
    declare @dE6 float = @dE4 * @dE2
    declare @dE7 float = @dE5 * @dE2
    set @φ = @φ - @VII * @dE2 + @VIII * @dE4 - @IX * @dE6
    declare @λ float = @λ0 + @x * @dE - @XI * @dE3 + @XII * @dE5 - @XIIA * @dE7
 
	-- Convert polar geodetic to cartesian geocentric
    declare @c_sinφ float = sin(@φ),
        @c_cosφ float = cos(@φ)
    declare @c_sinλ float = sin(@λ),
        @c_cosλ float = cos(@λ)

    declare @eSq float = 2 * @f - @f * @f                      -- 1st eccentricity squared ≡ (a²-b²)/a²
    declare @ν float = @a / sqrt(1 - @eSq * @c_sinφ * @c_sinφ) -- radius of curvature in prime vertical

    declare @os_x float = @ν * @c_cosφ * @c_cosλ
    declare @os_y float = @ν * @c_cosφ * @c_sinλ
    declare @os_z float = @ν * (1 - @eSq) * @c_sinφ

	-- This is the lat/long in OSGB36 www.ordnancesurvey.co.uk/docs/support/guide-coordinate-systems-great-britain.pdf
	-- We need it in WGS84 www.icao.int/safety/pbn/documentation/eurocontrol/eurocontrol
    declare @tx float = 446.448,			  -- x-shift
        @ty float     = -125.157,			  -- y-shift
        @tz float     = 542.060,			  -- z-shift
        @s float      = -20.4894 / 1e6 + 1,     -- scale: normalise parts-per-million to (s+1)
        @rx float     = radians(0.1502 / 3600), -- x-rotation: normalise arcseconds to radians
        @ry float     = radians(0.2470 / 3600), -- y-rotation: normalise arcseconds to radians
        @rz float	  = radians(0.8421 / 3600)  -- z-rotation: normalise arcseconds to radians

    -- Apply Helmert transform
    declare @x2 float = @tx + @os_x * @s - @os_y * @rz + @os_z * @ry
    declare @y2 float = @ty + @os_x * @rz + @os_y * @s - @os_z * @rx
    declare @z2 float = @tz - @os_x * @ry + @os_y * @rx + @os_z * @s

	-- WGS84 major axis (a), minor axis (b), and flattening (f) for the ellipsoid
    declare @wgs84_a float = 6378137,
        @wgs84_b float = 6356752.314245,
        @wgs84_f float = 1 / 298.257223563

    declare @wgs84_e2 float = 2 * @wgs84_f - @wgs84_f * @wgs84_f   -- 1st eccentricity squared ≡ (a²-b²)/a²
    declare @ε2 float = @wgs84_e2 / (1 - @wgs84_e2) -- 2nd eccentricity squared ≡ (a²-b²)/b²
    declare @p float = sqrt(@x2 * @x2 + @y2 * @y2); -- distance from minor axis
    declare @R float = sqrt(@p * @p + @z2 * @z2); -- polar radius

    -- parametric latitude (Bowring eqn 17, replacing tanβ = z·a / p·b)
    declare @tanβ float = (@wgs84_b * @z2) / (@wgs84_a * @p) * (1 + @ε2 * @wgs84_b / @R)
    declare @sinβ float = @tanβ / sqrt(1 + @tanβ * @tanβ)
    declare @cosβ float = @sinβ / @tanβ

    -- geodetic latitude (Bowring eqn 18: tanφ = z+ε²bsin³β / p−e²cos³β)
    declare @wgs84_φ float = atn2(@z2 + @ε2 * @wgs84_b * @sinβ * @sinβ * @sinβ, @p - @wgs84_e2 * @wgs84_a * @cosβ * @cosβ * @cosβ)

    -- longitude
    declare @wgs84_λ float = atn2(@y2, @x2)


	-- Return 
    return geography::Point(degrees(@wgs84_φ), degrees(@wgs84_λ), 4326)
end