Module:Ordnance Survey coordinates/sandbox: Difference between revisions

Revision as of 06:08, 16 December 2023

Latest revision as of 14:18, 4 November 2024

(55 intermediate revisions by the same user not shown)

Line 45:

local abs = math.abs

local floor = math.floor

local ceil = math.ceil

local sin = math.sin

local cos = math.cos

Line 51 ⟶ 52:

local deg2rad = math.rad

local rad2deg = math.deg

local dr = math.deg(1.0)

local wgs84radius = 6378137.0

local wgs84ecc = 0.00669438003551279089034150031998869922791

--[[ DATUM SHIFT USING HELMERT TRANSFORMATION

local scaleGB = 0.9996012717

*

local eastingGB = 400000.0

* height above the ellipsoid is ignored

local northingGB = -100000.0

* latitude and longitude must be given in decimal degrees

local lat0GB = 49.0

*

local lon0GB = -2.0

* no transformation if abs(latitude)>89 or abs(longitude)>89

local Airyradius = 6377563.396

* (lifting this restriction requires some more lines of code)

local Airyecc = 0.006670539761597529073698869358812557054558

*

* Written in 2009 by user Telford at de.wikipedia

* Based on math described in "A Guide to Coordinate Systems in Great Britain"

* by the UK Ordnance Survey

* URL: https://www.ordnancesurvey.co.uk/documents/resources/guide-coordinate-systems-great-britain.pdf

]]

-- Datum parameters

local data = mw.loadData('Module:Ordnance Survey coordinates/data')

local OSGBglobe = data.OSGBglobe

local IEglobe = data.IEglobe

local WGSglobe = data.WGSglobe

local WGS2OSGBparam = data.WGS2OSGBparam

local OSGB2WGSparam = data.OSGB2WGSparam

local IE2WGSparam = data.IE2WGSparam

local function HelmertDatumShift ( latitude , longitude, param )

-- latitude and longitude are in degrees

-- return if latitude or longitude out of range

if abs(latitude) > 89. or abs(longitude) > 89. then

return {lat=latitude, lon=longitude}

end

param = param or WGS2OSGBparam

-- parameters for ellipsoids

-- Annex A.1, Ordnance Survey document --

local a1 = param.semimajor_src

local b1 = param.semiminor_src

local e1 = param.ecc_src

local a2 = param.semimajor_dst

local b2 = param.semiminor_dst

local e2 = param.ecc_dst

-- convert latitude and longitude to cartesian coordinates

-- math in Annex B.1, Ordnance Survey document

local phi = deg2rad ( latitude )

local lambda = deg2rad ( longitude )

local cosphi = cos ( phi )

local sinphi = sin ( phi )

local coslambda = cos ( lambda )

local sinlambda = sin ( lambda )

local ny = a1 / sqrt ( 1. - e1*sinphi*sinphi )

local x1 = ny * cosphi * coslambda

local y1 = ny * cosphi * sinlambda

local z1 = ny * ( 1. - e1 ) * sinphi

-- the helmert transformation proper

-- Equation 3, Section 6.2, Ordnance Survey document

local x2 = x1 + param.tx + ( param.s0 * x1 - param.rz * y1 + param.ry * z1 )

local y2 = y1 + param.ty + ( param.rz * x1 + param.s0 * y1 - param.rx * z1 )

local z2 = z1 + param.tz + ( -param.ry *x1 + param.rx * y1 + param.s0 * z1 )

-- convert cartesian coordinates to latitude and longitude

-- math in Annex B.2, of Ordnance Survey document

lambda = atan ( y2 / x2 )

local p2 = sqrt ( x2*x2 + y2*y2 )

phi = atan ( z2 / p2*(1.-e2) )

local n0 = 101

local phi_alt

repeat

phi_alt = phi

ny = a2 / sqrt ( 1. - e2 * sin(phi) * sin(phi) )

phi = atan ( (z2 + e2 * ny * sin(phi)) / p2 )

n0 = n0 - 1

until abs ( phi - phi_alt ) <= 0.000000000001 or n0 <= 0

return {lat=rad2deg(phi),lon=rad2deg(lambda)}

end

-- LAT LONG TO OS GRID LIBRARY RESUMES HERE

local function northeast(lett,num,shift)

Line 96 ⟶ 179:

end

local function GBEN2LLEN2LL(e,n,datum)

local a = datum.semimajor*datum.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 lambda0ukn1 =lat0GB datum.n1

local n12 = n1*n1

-- scale factor @ central meridian

local F0ukn13 =scaleGB n12*n1

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

-- 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=Airyradius*F0uk

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

local buk=6356256.91*F0uk

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

local n1uk=0.00167322025032508731869331280635710896296

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

local e2uk=Airyecc

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+n1ukn1+1.25*pow(n1uk,2.0)n12+1.25*pow(n1uk,3.0)n13)*k3

j4=(3.0*n1ukn1+3.0*pow(n1uk,2.0)n12+2.625*pow(n1uk,3.0)n13)*sin(k3)*cos(k4)

j5=(1.875*pow(n1uk,2.0)n12+1.875*pow(n1uk,3.0)n13)*sin(2.0*k3)*cos(2.0*k4)

j6=35.0/24.0*pow(n1uk,3.0)n13*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 vx =auk/sqrt( 1.0-e2ukdatum.ecc*pow(sin(k),2.0))

local rv=v*(1.0-e2uk)a/sqrt(1.0-e2uk*pow(sin(k),2.0)x)

local r=v*(1.0-datum.ecc)/x

local h2=v/r-1.0

local y1=e-E0ukdatum.e0

j3local tank = tan(k)/(2.0*r*v)

local tank2 = tank*tank

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)

local tank4 = tank2*tank2

j5=tan(k)/(720.0*r*pow(v,5.0))*(61.0+90.0*pow(tan(k),2.0)+45.0*pow(tan(k),4.0))

local tank6 = tank2*tank4

local k9=k-y1*y1*j3+pow(y1,4.0)*j4-pow(y1,6.0)*j5

j6local cosk =1.0/( cos(k)*v)

local yv = y1/v -- dimensionless quantity in series expansion

local j7=1.0/(cos(k)*6.0*pow(v,3.0))*(v/r+2.0*pow(tan(k),2.0))

local yv2 = yv*yv

local j8=1.0/(cos(k)*120.0*pow(v,5.0))*(5.0+28.0*pow(tan(k),2.0)+24.0*pow(tan(k),4.0))

local j9yv3 =1.0/(cos(k) yv*5040.0*pow(v,7.0))yv2

local yv5 = yv3*yv2

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 yv7 = yv5*yv2

local long=phi0uk+dr*(y1*j6-y1*y1*y1*j7+pow(y1,5.0)*j8-pow(y1,7.0)*j9)

j3=tank/(2.0*r)

local lat=k9*dr

j4=tank/(24.0*r)*(5.0+3.0*tank2+h2-9.0*tank2*h2)

local v=Airyradius/sqrt(1.0-e2uk*pow(sin(k),2.0))

j5=tank/(720.0*r)*(61.0+90.0*tank2+45.0*tank4)

local cartxa=v*cos(k9)*cos(long/dr)

local cartyak9=vk-y1*cos(k9)yv*sin(long/drj3+yv3*j4-yv5*j5)

local cartzaj6=(1.0-e2uk)*v*sin/(k9cosk)

local j7=1.0/(cosk*6.0)*(v/r+2.0*tank2)

-- Helmert-Transformation from OSGB36 to WGS84 map date

local rotxj8=-01.15020/3600(cosk*120.0/dr)*(5.0+28.0*tank2+24.0*tank4)

local rotyj9=-1.0.2470/3600(cosk*5040.0/dr)

j9=j9*(61.0+662.0*tank2+1320.0*tank4+720.0*tank6)

local rotz=-0.8421/3600.0/dr

local long=datum.lon0+rad2deg(yv*j6-yv3*j7+yv5*j8-yv7*j9)

local scale=-20.4894/1000000.0

local lat=rad2deg(k9)

local cartxb=446.448+(1.0+scale)*cartxa+rotz*cartya-roty*cartza

return {lat=lat,lon=long}

local cartyb=-125.157-rotz*cartxa+(1.0+scale)*cartya+rotx*cartza

end

local cartzb=542.06+roty*cartxa-rotx*cartya+(1.0+scale)*cartza

-- convert Cartesian to long/lat

local function GBEN2LL(e,n)

local awgs84=wgs84radius

local bwgs84latlong =6356752.3141 EN2LL(e,n,OSGBglobe)

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

local e2wgs84=wgs84ecc

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

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

Line 194 ⟶ 253:

local function IrishEN2LL(e,n)

local latlong = EN2LL(e,n,IEglobe)

-- True Origin is 8 deg W

local helmert = HelmertDatumShift ( latlong.lat, latlong.lon, IE2WGSparam)

local phi0ir=-8.0

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

-- True Origin is 53.5 deg N

local lambda0ir=53.5

-- scale factor @ central meridian

local F0ir=1.000035

-- True origin in 200 km E of false origin

local E0ir=200000.0

--True origin is 250km N of false origin

local N0ir=250000.0

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

local air=6377340.189*F0ir

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

local bir=6356034.447*F0ir

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

local n1ir=0.001673220384152058651484728058385228837777

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

local e2ir=0.006670540293336110419293763349975612794125

local k=(n-N0ir)/air+lambda0ir/dr

local nextcounter=0

local j3,j4,j5,j6,m

repeat

nextcounter=nextcounter+1

local k3=k-lambda0ir/dr

local k4=k+lambda0ir/dr

j3=(1.0+n1ir+1.25*pow(n1ir,2.0)+1.25*pow(n1ir,3.0))*k3

j4=(3.0*n1ir+3.0*pow(n1ir,2.0)+2.625*pow(n1ir,3.0))*sin(k3)*cos(k4)

j5=(1.875*pow(n1ir,2.0)+1.875*pow(n1ir,3.0))*sin(2.0*k3)*cos(2.0*k4)

j6=35.0/24.0*pow(n1ir,3.0)*sin(3.0*k3)*cos(3.0*k4)

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

k=k+(n-N0ir-m)/air

until abs(n-N0ir-m)<=0.000000000001 or nextcounter>=10000

local v=air/sqrt(1.0-e2ir*pow(sin(k),2.0))

local r=v*(1.0-e2ir)/(1.0-e2ir*pow(sin(k),2.0))

local h2=v/r-1.0

local y1=e-E0ir

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)

j5=tan(k)/(720.0*r*pow(v,5.0))*(61.0+90.0*pow(tan(k),2.0)+45.0*pow(tan(k),4.0))

local k9=k-y1*y1*j3+pow(y1,4.0)*j4-pow(y1,6.0)*j5

j6=1.0/(cos(k)*v)

local j7=1.0/(cos(k)*6.0*pow(v,3.0))*(v/r+2.0*pow(tan(k),2.0))

local j8=1.0/(cos(k)*120.0*pow(v,5.0))*(5.0+28.0*pow(tan(k),2.0)+24.0*pow(tan(k),4.0))

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=phi0ir+dr*(y1*j6-y1*y1*y1*j7+pow(y1,5.0)*j8-pow(y1,7.0)*j9)

local lat=k9*dr

-- convert long/lat to Cartesian coordinates

v=6377340.189/sqrt(1.0-e2ir*pow(sin(k),2.0))

local cartxa=v*cos(k9)*cos(long/dr)

local cartya=v*cos(k9)*sin(long/dr)

local cartza=(1.0-e2ir)*v*sin(k9)

-- Helmert-Transformation from Ireland 1965 to WGS84 map date

local rotx=1.042/3600.0/dr

local roty=0.214/3600.0/dr

local rotz=0.631/3600.0/dr

local scale=8.15/1000000.0

local cartxb=482.53+(1.0+scale)*cartxa+rotz*cartya-roty*cartza

local cartyb=-130.596-rotz*cartxa+(1.0+scale)*cartya+rotx*cartza

local cartzb=564.557+roty*cartxa-rotx*cartya+(1.0+scale)*cartza

-- convert Cartesian to long/lat

local awgs84=wgs84radius

local bwgs84=6356752.3141

local e2wgs84=wgs84ecc

local lambdaradwgs84=atan(cartyb/cartxb)

long=lambdaradwgs84*dr

local pxy=sqrt(pow(cartxb,2.0)+pow(cartyb,2.0))

local phinewwgs84=atan(cartzb/pxy/(1.0-e2wgs84))

local phiradwgs84

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>=10000

lat=phinewwgs84*dr

return {region="IE",lat=lat,long=long}

end

Line 296 ⟶ 280:

local function NGR2LL(ngr)

local result = {}

local _ = 0

ngr, _ = mw.ustring.gsub(mw.ustring.upper(ngr),"[%s%p]","")

local first, last, lett, num = mw.ustring.find(ngr,"^([A-Z]+)(%d+)$")

Line 318 ⟶ 303:

local function trim(s)

local _ = 0

s, _ = mw.ustring.gsub(s,"^%s+","")

s, _ = mw.ustring.gsub(s,"%s+$","")

Line 339 ⟶ 325:

end

-- generate URL of geohack map

function oscoord.main(frame)

local function geohack(main_args, other_args, LL)

local args = getArgs(frame)

-- create geohack link. Example:

local input = args[1]

-- https://geohack.toolforge.org/geohack.php?pagename=Mount_Whitney&params=36.578580925_N_118.29199495_W_type:mountain_region:US-CA_scale:100000_source:NGS

local url = 'https://geohack.toolforge.org/geohack.php?'

local input = main_args[1]

local namearg = main_args["name"]

local current_page = mw.title.getCurrentTitle()

local pagename = mw.uri.encode( current_page.prefixedText, 'WIKI' )

if not empty(pagename) then

url = url..'pagename='..pagename..'&'

end

url = url..'params='..mw.ustring.format('%.6f',LL.lat)..'_N_'

if LL.long < 0 then

url = url..mw.ustring.format('%.6f',-LL.long)..'_W'

else

url = url..mw.ustring.format('%.6f',LL.long)..'_E'

end

for _, a in ipairs(other_args) do

url = url..'_'..a

end

if not mw.ustring.find(input,"region") and LL.region then

url = url..'_region:'..LL.region

end

if not mw.ustring.find(input,"scale") and

not mw.ustring.find(input,"type") and

not mw.ustring.find(input,"dim") and LL.prec then

url = url..'_dim:'..floor(50*LL.prec+0.5)..'m'

end

if not empty(namearg) then

url = url .. "&title=" .. mw.uri.encode(namearg)

end

return url

end

-- function to generate direct link to OS map

local function directLink(main_args, other_args, LL)

-- create link to Bing server for OS maps. Example:

-- https://www.bing.com/maps/?mkt=en-gb&v=2&cp=56.796026%7E-5.01307&lvl=16.0&sp=Point.56.796029_-5.004711_Ben+Nevis&sty=s&style=s

local current_page = mw.title.getCurrentTitle()

local namearg = mw.uri.encode( main_args["name"] or current_page.prefixedText, 'QUERY' )

local args = {}

for _, a in ipairs(other_args) do

local splitOut = mw.text.split(a, ':', true)

args[splitOut[1]] = splitOut[2]

end

if not args.scale and not args.type and not args.dim then

args.dim = LL.prec and tostring(floor(50*LL.prec+0.5))..'m'

end

args.viewport_cm = 10

local zoom = require('Module:Infobox dim')._zoom

local lvl = zoom(args) or 12

local url = mw.ustring.format('https://www.bing.com/maps/?mkt=en-gb&v=2&cp=%.6f~%.6f&lvl=%d&sp=Point.%.6f_%.6f',LL.lat,LL.long,lvl,LL.lat,LL.long)

if not empty(namearg) then

url = url..'_'..namearg

end

url = url..'&sty=s&style=s'

return url

end

function oscoord._main(main_args)

local input = main_args[1]

if empty(input) then

return warning(nil)

end

local linktitle = argsmain_args[2]

local nameargrawurl = argsyesno(main_args["namerawurl"])

local rawurldirect = yesno(argsmain_args["rawurldirect"])

local args = split(input,'_')

local LL

local restargs = 1

local current_page = mw.title.getCurrentTitle()

local pagename = mw.uri.encode( current_page.prefixedText, 'WIKI' )

if #args >= 2 and alldigits(args[2]) then

if mw.ustring.sub(args[1],1,1) == 'i' then

Line 380 ⟶ 423:

if not empty(LL.err) then

return linktitle ..warning(LL.err)

end

-- https://geohack.toolforge.org/geohack.php?pagename=Mount_Whitney&params=36.578580925_N_118.29199495_W_type:mountain_region:US-CA_scale:100000_source:NGS

local url = ''

if not rawurl then

url = url..'['

end

url = url..'https://geohack.toolforge.org/geohack.php?'

if not empty(pagename) then

url = url..'pagename='..pagename..'&'

end

LL.lat = LL.lat or 0

LL.long = LL.long or 0

urlLL.lat = url..'params='..ceil(LL.lat..'_N_'*1000000)/1000000

if LL.long <= 0 thenceil(LL.long*1000000)/1000000

local other_args = {}

url = url..(-LL.long)..'_W'

for i = restargs, #args do

table.insert(other_args, args[i])

end

local url

if not direct then

url = geohack(main_args, other_args, LL)

else

url = directLink(main_args, other_args, LL)

url = url..LL.long..'_E'

end

for i = restargs,#args do

url = url..'_'..args[i]

end

if not mw.ustring.find(input,"region") and LL.region then

url = url..'_region:'..LL.region

end

if not mw.ustring.find(input,"scale") and

not mw.ustring.find(input,"type") and

not mw.ustring.find(input,"dim") and LL.prec then

url = url..'_dim:'..floor(50*LL.prec+0.5)..'m'

end

if not empty(namearg) then

url = url .. "&title=" .. mw.uri.encode(namearg)

end

if not rawurl then

url = '['..url..' '..linktitle..']'

end

return url

end

function oscoord.oscoord_oscoord(frameargs)

local output = '<span class="plainlinks nourlexpansion" style="white-space: nowrap">' .. oscoord.main_main(frameargs) .. '</span>'

if namespace == 0 then

output = output .. '[[Category:Articles with OS grid coordinates]]'

end

return output

end

function oscoord.main(frame)

local args = getArgs(frame)

return oscoord._main(args)

end

function oscoord.oscoord(frame)

local args = getArgs(frame)

return oscoord._oscoord(args)

end

--[[

****** LAT/LONG CONVERSION TO OS GRID REF FUNCTIONS *****

* Convert latitude longitude geographical coordinates to

* Uses the WGS-84 ellipsoid and only works for GB grid ref

* Transverse Mercator coordinates.

*

* Uses the WGS-84 ellipsoid by default

*

* See also:

Line 458 ⟶ 492:

* along with this program; if not, write to the Free Software

* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

Converted to Lua by User:Hike395 on 2023-12-15

]]

--[[

* datum shift using Helmert transformation

* parameters for transformation WGS84 --> OSGB36 hard coded

* possibility to implement other parameters should be evident

*

* height above the ellipsoid is ignored

* latitude and longitude must be given in decimal degrees

*

* no transformation if abs(latitude)>89 or abs(longitude)>89

* (lifting this restriction requires some more lines of code)

*

* Written in 2009 by user Telford at de.wikipedia

* Based on math described in "A Guide to Coordinate Systems in Great Britain"

* by the UK Ordnance Survey

* URL: https://www.ordnancesurvey.co.uk/documents/resources/guide-coordinate-systems-great-britain.pdf

]]

local function HelmertDatumShift ( latitude , longitude )

-- return if latitude or longitude out of range

if abs(latitude) > 89. or abs(longitude) > 89. then

return {lat=latitude, lon=longitude}

end

-- parameters for ellipsoids

-- Annex A.1, Ordnance Survey document --

local a1 = wgs84radius

local b1 = 6356752.3141

local a2 = Airyradius

local b2 = 6356256.910

-- parameters for transformation WGS84 --> OSGB36

-- Table 4, Section 6.6, Ordnance Survey document

local tx = -446.448

local ty = 125.157

local tz = -542.060

local s0 = 1.0000204894

local rx = deg2rad ( -.1502/3600. )

local ry = deg2rad ( -.2470/3600. )

local rz = deg2rad ( -.8241/3600. )

-- convert latitude and longitude to cartesian coordinates

-- math in Annex B.1, Ordnance Survey document

local phi = deg2rad ( latitude )

local lambda = deg2rad ( longitude )

local cosphi = cos ( phi )

local sinphi = sin ( phi )

local coslambda = cos ( lambda )

local sinlambda = sin ( lambda )

local aq = a1 * a1

local e1 = ( aq - b1 * b1 ) / aq

local ny = a1 / sqrt ( 1. - e1*sinphi*sinphi )

local x1 = ny * cosphi * coslambda

local y1 = ny * cosphi * sinlambda

local z1 = ny * ( 1. - e1 ) * sinphi

-- the helmert transformation proper

-- Equation 3, Section 6.2, Ordnance Survey document

local x2 = tx + ( s0 * x1 - rz * y1 + ry * z1 )

local y2 = ty + ( rz * x1 + s0 * y1 - rx * z1 )

local z2 = tz + ( -ry *x1 + rx * y1 + s0 * z1 )

-- convert cartesian coordinates to latitude and longitude

-- math in Annex B.2, of Ordnance Survey document

local aq = a2 * a2

local e2 = ( aq - b2 * b2 ) / aq

lambda = atan ( y2 , x2 )

local p2 = sqrt ( x2*x2 + y2*y2 )

phi = atan ( z2 , p2*(1.-e2) )

local limit = 1. / b2

local n0 = 11

repeat

local phi_alt = phi

ny = a2 / sqrt ( 1. - e2 * sin(phi) * sin(phi) )

phi = atan ( z2 + e2 * ny * sin(phi) , p2 )

n0 = n0 - 1

until abs ( phi - phi_alt ) <= limit or n0 <= 0

return {lat=rad2deg(phi),lon=rad2deg(lambda)}

end

local function find_M ( lat_rad )

local e = AiryeccOSGBglobe.ecc

local e2 = e*e

local e3 = e*e*e

return AiryradiusOSGBglobe.semimajor * ( ( 1 - e/4 - 3 * e * ee2/64

- 5 * e * e * ee3/256

) * lat_rad

- ( 3 * e/8 + 3 * e * ee2/32

+ 45 * e * e * ee3/1024

) * sin(2 * lat_rad)

+ ( 15 * e * ee2/256 +

45 * e * e * ee3/1024

) * sin(4 * lat_rad)

- ( 35 * e * e * ee3/3072

) * sin(6 * lat_rad) )

end

Line 588 ⟶ 530:

local lat_rad = deg2rad( latitude )

local e = AiryeccOSGBglobe.ecc

local e_prime_sq = e / (1-e)

local v = AiryradiusOSGBglobe.semimajor / sqrt(1 - e * sin(lat_rad)*sin(lat_rad))

local Ttank = pow( tan(lat_rad), 2)

local T = tank*tank

local T2 = T*T

local C = e_prime_sq * pow( cos(lat_rad), 2)

local A = deg2rad( longitude2 -longitude_originOSGBglobe.lon0 ) * cos(lat_rad)

local A2 = A*A

local A3 = A2*A

local A4 = A2*A2

local A5 = A3*A2

local A6 = A3*A3

local M = find_M( lat_rad )

local M0 = 0.0

if latitude_origin ~= 0 then

M0 = find_M(deg2rad( latitude_originOSGBglobe.lat0 ))

end

local northing = northingGBOSGBglobe.n0 + scaleGBOSGBglobe.scale *

( (M - M0) + v*tan(lat_rad) *

( A*AA2/2

+ (5 - T + 9*C + 4*C*C) * pow(A,4)A4/24

+ (61 - 58*T + T*T T2

+ 600*C - 330*e_prime_sq) * pow(A,6)A6/720 ))

local easting = eastingGBOSGBglobe.e0 + scaleGBOSGBglobe.scale * v *

( A

+ (1-T+C)*pow(A,3)A3/6

+ (5 - 18*T + pow(T,2)T2 + 72*C

- 58 * e_prime_sq)*pow(A,5)A5/120 )

return {northing=northing,easting=easting}

Line 621 ⟶ 570:

* OSBG36 Easting and Northing

]]

local function LatLon2OSGB36( latlon, prec )

local tm = LatLonOrigin2TM(latlon.lat, latlon.lon)

if not tm then return 'TM failure'..tostring(latlon.lat)..' '..tostring(latlon.lon) end

if not tonumber(prec) then prec = 5 end

prec = floor(prec+0.5)

if prec > 5 then prec = 5 end

if prec < 1 then prec = 1 end

-- fix by Roger W Haworth

Line 633 ⟶ 586:

local letters = "ABCDEFGHJKLMNOPQRSTUVWXYZ"

local c1indx1 = mw.ustring.sub(letters,18-floor(grid_y/5)*5+floor(grid_x/5),1)

local c2indx2 = mw.ustring.sub(letters,21-(grid_y%5)*5+grid_x%5,1)

local c1 = mw.ustring.sub(letters,indx1,indx1)

local c2 = mw.ustring.sub(letters,indx2,indx2)

local eeasting = mw.ustring.format('%05d', tm.easting%100000)

local nnorthing = mw.ustring.format('%05d', tm.northing%100000)

local grid = pow(10.0,5.0-prec)

easting = floor(easting/grid)

northing = floor(northing/grid)

local fmt = '%0'..prec..'d'

local e = mw.ustring.format(fmt, easting)

local n = mw.ustring.format(fmt, northing)

return c1..c2..e..n

Line 643 ⟶ 604:

end

function oscoord._WGS2OSGB(lat,lon,prec)

return LatLon2OSGB36(HelmertDatumShift(lat,lon,WGS2OSGBparam),prec)

end

function oscoord.WGS2OSGB(frame)

local args = getArgs(frame)

return args[1] and args[2] and oscoord._WGS2OSGB(args[1],args[2],args[3]) or ''

end

function oscoord.LL2OS(frame)

local args = getArgs(frame)

if not args[1] or not args[2] then return '' end

local gridRef = oscoord._WGS2OSGB(args[1],args[2],args.prec)

if not gridRef or #gridRef == 0 then return '' end

if args[3] then

gridRef = gridRef..'_'..args[3]

end

return oscoord._oscoord({gridRef,args.linktitle,name=args.name})

end