|
local function HelmertDatumShift ( latitude , longitude, param )
-- latitude and longitude are in degrees
-- return if latitude or longitude out of range
local long=phi0uk+dr*(y1*j6-y1*y1*y1*j7+pow(y1,5.0)*j8-pow(y1,7.0)*j9)
local lat=k9*dr
local helmert = HelmertDatumShift ( lat, long, OSGB2WGSparam)
local v=OSGBglobe.semimajor/sqrt(1.0-e2uk*pow(sin(k),2.0))
return {region="GB",lat= helmert.lat,long= longhelmert.lon} ▼ local cartxa=v*cos(k9)*cos(long/dr)
local cartya=v*cos(k9)*sin(long/dr)
local cartza=(1.0-e2uk)*v*sin(k9)
-- Helmert-Transformation from OSGB36 to WGS84 map date
local rotx=-0.1502/3600.0/dr
local roty=-0.2470/3600.0/dr
local rotz=-0.8421/3600.0/dr
local scale=-20.4894/1000000.0
local cartxb=446.448+(1.0+scale)*cartxa+rotz*cartya-roty*cartza
local cartyb=-125.157-rotz*cartxa+(1.0+scale)*cartya+rotx*cartza
local cartzb=542.06+roty*cartxa-rotx*cartya+(1.0+scale)*cartza
-- convert Cartesian to long/lat
local awgs84=WGSglobe.semimajor
local bwgs84=WGSglobe.semiminor
local e2wgs84=WGSglobe.ecc
local lambdaradwgs84=atan(cartyb/cartxb)
long=lambdaradwgs84*dr
local pxy=sqrt(pow(cartxb,2.0)+pow(cartyb,2.0))
local phiradwgs84
local phinewwgs84=atan(cartzb/pxy/(1.0-e2wgs84))
nextcounter=0
repeat
phiradwgs84=phinewwgs84
nextcounter=nextcounter+1
v=awgs84/sqrt(1.0-e2wgs84*pow(sin(phiradwgs84),2.0))
phinewwgs84=atan((cartzb+e2wgs84*v*sin(phiradwgs84))/pxy)
until abs(phinewwgs84-phiradwgs84)<=0.000000000001 or nextcounter>=100
lat=phinewwgs84*dr
▲ return {region="GB",lat=lat,long=long}
end
| 