Content deleted Content added
use precomputed ecc |
both code paths |
||
Line 271:
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=helmert.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
end
| |||