Module:Ordnance Survey coordinates/sandbox: Difference between revisions

Revision as of 19:14, 16 December 2023

Revision as of 19:34, 16 December 2023

Line 53:

local dr = math.deg(1.0)

-- Transverse Mercator parameters

local scaleGB = 0.9996012717

local eastingGB = 400000.0

local northingGB = -100000.0

local lat0GB = 49.0

local lon0GB = -2.0

Line 80 ⟶ 74:

semimajor = 6377563.396,

semiminor = 6356256.910,

ecc = 0.006670539761597529073698869358812557054558,

n1 = local n1uk=0.00167322025032508731869331280635710896296,

local scaleGB scale = 0.9996012717,

local eastingGB e0 = 400000.0,

local northingGB n0 = -100000.0,

lat0 = 49.0,

lon0 = -2.0

}

Line 243:

end

local function GBEN2LLEN2LL(e,n,datum)

local auka = OSGBglobedatum.semimajor*F0ukdatum.scale

-- True Origin is 2 deg W

local phi0ukb =lon0GB datum.semiminor*datum.scale

local lat0rad = deg2rad(datum.lat0)

-- True Origin is 49 deg N

local lambda0ukk=lat0GB(n-datum.n0)/a+lat0rad

-- scale factor @ central meridian

local F0uk=scaleGB

-- True origin in 400 km E of false origin

local E0uk=eastingGB

--True origin is 100 km S of false origin

local N0uk=northingGB

-- semi-major axis (in line to equator) 0.996012717 is yer scale @ central meridian

local auk= OSGBglobe.semimajor*F0uk

--semi-minor axis (in line to poles)

local buk=OSGBglobe.semiminor*F0uk

-- flatness=a1-b1/(a1+b1)

local n1uk=0.00167322025032508731869331280635710896296

-- first eccentricity squared=2*f-f^2where f=(a1-b1)/a1

local e2uk=OSGBglobe.ecc

local k=(n-N0uk)/auk+lambda0uk/dr

local nextcounter=0

local j3, j4, j5, j6, m

repeat

nextcounter=nextcounter+1

local k3=k-lambda0uk/drlat0rad

local k4=k+lambda0uk/drlat0rad

j3=(1.0+n1ukdatum.n1+1.25*pow(n1ukdatum.n1,2.0)+1.25*pow(n1ukdatum.n1,3.0))*k3

j4=(3.0*n1ukdatum.n1+3.0*pow(n1ukdatum.n1,2.0)+2.625*pow(n1ukdatum.n1,3.0))*sin(k3)*cos(k4)

j5=(1.875*pow(n1ukdatum.n1,2.0)+1.875*pow(n1ukdatum.n1,3.0))*sin(2.0*k3)*cos(2.0*k4)

j6=35.0/24.0*pow(n1ukdatum.n1,3.0)*sin(3.0*k3)*cos(3.0*k4)

m=bukb*(j3-j4+j5-j6)

k=k+(n-N0ukdatum.n0-m)/auka

until abs(n-N0ukdatum.n0-m)<=0.000000001 or nextcounter>=100

local v=auka/sqrt(1.0-e2ukdatum.ecc*pow(sin(k),2.0))

local r=v*(1.0-e2ukdatum.ecc)/(1.0-e2ukdatum.ecc*pow(sin(k),2.0))

local h2=v/r-1.0

local y1=e-E0ukdatum.e0

j3=tan(k)/(2.0*r*v)

j4=tan(k)/(24.0*r*pow(v,3.0))*(5.0+3.0*pow(tan(k),2.0)+h2-9.0*pow(tan(k),2.0)*h2)

Line 289 ⟶ 274:

local j9=1.0/(cos(k)*5040.0*pow(v,7.0))

local j9=j9*(61.0+662.0*pow(tan(k),2.0)+1320.0*pow(tan(k),4.0)+720.0*pow(tan(k),6.0))

local long=phi0ukdatum.lon0+dr*rad2deg(y1*j6-y1*y1*y1*j7+pow(y1,5.0)*j8-pow(y1,7.0)*j9)

local lat=rad2deg(k9*dr)

return {lat=lat,lon=long}

local helmert = HelmertDatumShift ( lat, long, OSGB2WGSparam)

end

local function GBEN2LL(e,n)

local e2uklatlong = EN2LL(e,n,OSGBglobe.ecc)

local helmert = HelmertDatumShift ( latlong.lat, longlatlong.lon, OSGB2WGSparam)

return {region="GB",lat=helmert.lat,long=helmert.lon}

end