I know this is and old post but thank you for sharing this code. I found it useful.
I converted your code/formulas above into T-sql functions, one for lat and one for long. Here is the latitude one in case anyone is interested for Florida State Plane East 0901 EPSG 2236:
SET ANSI_NULLS ON
GO
SET QUOTED_IDENTIFIER ON
GO
Alter FUNCTION [dbo].[Langd_Convert_Feet_to_DecimalDegrees_Latitude]
(
-- Add the parameters for the function here
@y Numeric(38,8), -- latitude
@x Numeric(38,8) -- longitude
)
RETURNS Numeric(38,8)
AS
BEGIN
-- Constants below from NAD 1983 State Plane Florida East FIPS 0901(Feet).
-- e_sqr is for NAD83 or GRS80 (http:en.wikipedia.org/wiki/NAD83#North_American_Datum_of_1983)
/* From Esri WKT https://epsg.io/2236
PROJCS["NAD83_Florida_East_ftUS",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID
["GRS_1980",6378137,298.257222101]],
PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",24.33333333333333],PARAMETER["central_meridian",-81],
PARAMETER["scale_factor",0.999941177],PARAMETER["false_easting",656166.667],PARAMETER["false_northing",0],
UNIT["Foot_US",0.30480060960121924]]
*/
Declare @false_easting Numeric(38,8) = 656166.667 -- original formula was (656166.6666666665 + (2/3)); -- false easting
Declare @false_northing Numeric(38,8) = 0 -- original formula was (0 + (1/3)) -- Should be 0. this is a 'fudge' factor.
Declare @a Numeric(38,8) = 6378137.0 --Equatorial Radius (meters) For GRS80
Declare @e_sqr Numeric(38,8) = 0.006694300229 --Eccentricity of ellipsoid (GRS80)
Declare @lat_o Numeric(38,8) = 24.33333333333333 --Origin Latitude FL state plane
Declare @long_o Numeric(38,8) = 81.0 -- original formula was (81.0 + ((1/3) * 0.01)) -- Central Meridan FL state plane, no negative sign for this formula
Declare @k_o Numeric(38,8) = 0.999941177 --0.9999411764705882 -- Scale factor
Declare @pointx_ft Numeric(38,8) = @x -- number ft
Set @pointx_ft = @x - @false_easting -- false easting accounted for.
Declare @pointy_ft Numeric(38,8) = @y --number ft
Set @pointy_ft = @y - @false_northing -- false northing accounted for.
Declare @e_sqr_2 Numeric(38,8) = power(@e_sqr, 2) -- need this to calculate M
Declare @e_sqr_3 Numeric(38,8) = power(@e_sqr, 3) -- need this to calculate M
Declare @radianlat_o Numeric(38,8) = @lat_o/(180/PI()) --degrees converted to radians 57.2957
Declare @radianLong_o Numeric(38,8) = @long_o/(180/PI()) --degrees converted to radians
Declare @pointx_m Numeric(38,8) = @pointx_ft/3.2808 --@pointx_ft/3.2808399
Declare @pointy_m Numeric(38,8) = @pointy_ft/3.2808 --@pointy_ft/3.2808399
Declare @Mpart1_o Numeric(38,8) = ((1-(@e_sqr/4)-((3*@e_sqr_2)/64)-((5*@e_sqr_3)/256))*@radianlat_o)
Declare @Mpart2_o Numeric(38,8) = ((((3*@e_sqr)/8)+((3*@e_sqr_2)/32)+((45*@e_sqr_3)/1024))*sin(2*@radianlat_o))
Declare @Mpart3_o Numeric(38,8) = ((((15*@e_sqr_2)/256) + ((45*@e_sqr_3)/1024))*4*sin(4*@radianlat_o) )
Declare @Mpart4_o Numeric(38,8) = ((35*@e_sqr_3)/3072)*4*sin(6*@radianlat_o)
Declare @M_o Numeric(38,8) = @a*(@Mpart1_o - @Mpart2_o + @Mpart3_o - @Mpart4_o)
Declare @e_prim_sqr Numeric(38,8) = @e_sqr / (1 - @e_sqr)
Declare @M Numeric(38,8) = @M_o + (@pointy_m / @k_o)
Declare @e1 Numeric(38,8) = (1-sqrt((1 - @e_sqr))) / (1+sqrt((1 - @e_sqr))) --note: constant only for @e_sqr of 0.00676866
Declare @u Numeric(38,8) = @M / ( @a*( 1-(@e_sqr/4)- (3*(power(@e_sqr,2)/64) - (5*(power(@e_sqr,3)/256) ))) ) --radians
Declare @lat1_part1 Numeric(38,8) = ((3 * (@e1/2)) - (27 * (power(@e1,3)/32) )) * sin(2*@u)
Declare @lat1_part2 Numeric(38,8) = ((21 * (power(@e1,2)/16)) - (55 * (power(@e1,4)/32) )) * sin(4*@u)
Declare @lat1_part3 Numeric(38,8) = ((151 * (power(@e1,3)/96)) * sin(6*@u) )
Declare @lat1_part4 Numeric(38,8) = ((1097 * (power(@e1,4)/512)) * sin(8*@u) )
Declare @lat1 Numeric(38,8) = @u + @lat1_part1 + @lat1_part2 + @lat1_part3 + @lat1_part4
Set @lat1 = @lat1 * (180/PI())
Declare @radian_lat1 Numeric(38,8) = @lat1/(180/PI()) --degrees converted to radians
Declare @C1 Numeric(38,8) = @e_prim_sqr * (0.5*(1 + cos(2*@radian_lat1)))
Declare @T1 Numeric(38,8) = (1-cos(2*@radian_lat1))/(1+cos(2*@radian_lat1))
Declare @N1 Numeric(38,8) = @a/(sqrt(1-(@e_sqr*0.5*(1-cos(2*@radian_lat1)))))
Declare @R1 Numeric(38,8) = (@a*(1-@e_sqr))/ sqrt(power((1-@e_sqr * (0.5*(1 - cos(2*@radian_lat1))) ),3) )
Declare @D Numeric(38,8) = @pointx_m/(@N1*@k_o)
Declare @Lat_final_part1 Numeric(38,8) = @lat1
Declare @Lat_final_part2 Numeric(38,8) = (@N1 * (tan(@radian_lat1) / @R1))
Declare @Lat_final_part3 Numeric(38,8) = power(@D,2)/2 - ((( 5 + (3*@T1) + (10*@C1) - (4* power(@C1,2)) - (9 * @e_prim_sqr)) * power(@D,4))/24)
Declare @Lat_final_part4 Numeric(38,8) = ((61 + (90 * @T1) + (298 * @C1) + (45 * power(@T1,2)) - (252 * @e_prim_sqr) - (3 * power(@C1,2)))* power(@D,6))/720
Declare @Lat_final Numeric(38,8) = @Lat_final_part1 - ((@Lat_final_part2 *(@Lat_final_part3 + @Lat_final_part4)) * (180/PI()))
Declare @Long_final_part1 Numeric(38,8) = @long_o
Declare @Long_final_part2 Numeric(38,8) = @D --D
Declare @Long_final_part3 Numeric(38,8) = (((1 + (2*@T1) + @C1))* (power(@D,3))) /6
Declare @Long_final_part4 Numeric(38,8) = 5 - (2 * @C1) + (28*@T1) - (3 * power(@C1,2)) + (8 * @e_prim_sqr) + (24 * power(@T1,2))
Declare @Long_final_part4a Numeric(38,8) = @Long_final_part4 * (power(@D,5) )
Declare @Long_final_part4b Numeric(38,8) = @Long_final_part4a /120
Declare @Long_final Numeric(38,8) = (-@Long_final_part1 + (((@Long_final_part2 - (@Long_final_part3 + @Long_final_part4b))) /cos(@radian_lat1)) * 180/PI())
RETURN @Lat_final
END