Rewriting a Java Script Function

[Simon Knight]

Hi,
I realise that I'm well off topic but I suspect that many live coders have tried other languages so I'll ask. I am rewriting a function that converts a position in latitude and longitude into a UTM reference. It should just be straight forward maths but as I have never really used languages such as Java Script I think it worth checking my interpretation with experts! So here is the code.

Code: Select all

// adjust zone & central meridian for Svalbard
    if (zone==32 && latBand=='X' && this.lon<  9) { zone--; λ0 -= (6).toRadians(); }

I'm reading this as if variable zone is 32 AND latBand is X and lon is less than 9 then do the stuff inside the braces:
1) subtract 1 from zone
2) convert 6 to radians and subtract the radians value from λ0 and store the result back in λ0

How did I do?
Thanks
Simon K.

[capellan]

Hi Simon,

Are you using the JavaScript code of this webpage?
https://www.movable-type.co.uk/scripts/ ... -mgrs.html

Al

[Simon Knight]

Hi Al,

Yes. Although I am only converting the function that converts latitude and longitude to UTM coordinates:

/**
* Converts latitude/longitude to UTM coordinate.
*
* Implements Karney’s method, using Krüger series to order n^6, giving results accurate to 5nm for
* distances up to 3900km from the central meridian.
*
* @returns {Utm} UTM coordinate.
* @throws {Error} If point not valid, if point outside latitude range.
*
* @example
* var latlong = new LatLon(48.8582, 2.2945);
* var utmCoord = latlong.toUtm(); // utmCoord.toString(): '31 N 448252 5411933'
*/
LatLon.prototype.toUtm = function() {
if (isNaN(this.lat) || isNaN(this.lon)) throw new Error('Invalid point');
if (!(-80<=this.lat && this.lat<=84)) throw new Error('Outside UTM limits');

var falseEasting = 500e3, falseNorthing = 10000e3;

var zone = Math.floor((this.lon+180)/6) + 1; // longitudinal zone
var λ0 = ((zone-1)*6 - 180 + 3).toRadians(); // longitude of central meridian

// ---- handle Norway/Svalbard exceptions
// grid zones are 8° tall; 0°N is offset 10 into latitude bands array
var mgrsLatBands = 'CDEFGHJKLMNPQRSTUVWXX'; // X is repeated for 80-84°N
var latBand = mgrsLatBands.charAt(Math.floor(this.lat/8+10));
// adjust zone & central meridian for Norway
if (zone==31 && latBand=='V' && this.lon>= 3) { zone++; λ0 += (6).toRadians(); }
// adjust zone & central meridian for Svalbard
if (zone==32 && latBand=='X' && this.lon< 9) { zone--; λ0 -= (6).toRadians(); }
if (zone==32 && latBand=='X' && this.lon>= 9) { zone++; λ0 += (6).toRadians(); }
if (zone==34 && latBand=='X' && this.lon< 21) { zone--; λ0 -= (6).toRadians(); }
if (zone==34 && latBand=='X' && this.lon>=21) { zone++; λ0 += (6).toRadians(); }
if (zone==36 && latBand=='X' && this.lon< 33) { zone--; λ0 -= (6).toRadians(); }
if (zone==36 && latBand=='X' && this.lon>=33) { zone++; λ0 += (6).toRadians(); }

var φ = this.lat.toRadians(); // latitude ± from equator
var λ = this.lon.toRadians() - λ0; // longitude ± from central meridian

var a = this.datum.ellipsoid.a, f = this.datum.ellipsoid.f;
// WGS 84: a = 6378137, b = 6356752.314245, f = 1/298.257223563;

var k0 = 0.9996; // UTM scale on the central meridian

// ---- easting, northing: Karney 2011 Eq 7-14, 29, 35:

var e = Math.sqrt(f*(2-f)); // eccentricity
var n = f / (2 - f); // 3rd flattening
var n2 = n*n, n3 = n*n2, n4 = n*n3, n5 = n*n4, n6 = n*n5; // TODO: compare Horner-form accuracy?

var cosλ = Math.cos(λ), sinλ = Math.sin(λ), tanλ = Math.tan(λ);

var τ = Math.tan(φ); // τ ≡ tanφ, τʹ ≡ tanφʹ; prime (ʹ) indicates angles on the conformal sphere
var σ = Math.sinh(e*Math.atanh(e*τ/Math.sqrt(1+τ*τ)));

var τʹ = τ*Math.sqrt(1+σ*σ) - σ*Math.sqrt(1+τ*τ);

var ξʹ = Math.atan2(τʹ, cosλ);
var ηʹ = Math.asinh(sinλ / Math.sqrt(τʹ*τʹ + cosλ*cosλ));

var A = a/(1+n) * (1 + 1/4*n2 + 1/64*n4 + 1/256*n6); // 2πA is the circumference of a meridian

var α = [ null, // note α is one-based array (6th order Krüger expressions)
1/2*n - 2/3*n2 + 5/16*n3 + 41/180*n4 - 127/288*n5 + 7891/37800*n6,
13/48*n2 - 3/5*n3 + 557/1440*n4 + 281/630*n5 - 1983433/1935360*n6,
61/240*n3 - 103/140*n4 + 15061/26880*n5 + 167603/181440*n6,
49561/161280*n4 - 179/168*n5 + 6601661/7257600*n6,
34729/80640*n5 - 3418889/1995840*n6,
212378941/319334400*n6 ];

var ξ = ξʹ;
for (var j=1; j<=6; j++) ξ += α[j] * Math.sin(2*j*ξʹ) * Math.cosh(2*j*ηʹ);

var η = ηʹ;
for (var j=1; j<=6; j++) η += α[j] * Math.cos(2*j*ξʹ) * Math.sinh(2*j*ηʹ);

var x = k0 * A * η;
var y = k0 * A * ξ;

// ---- convergence: Karney 2011 Eq 23, 24

var pʹ = 1;
for (var j=1; j<=6; j++) pʹ += 2*j*α[j] * Math.cos(2*j*ξʹ) * Math.cosh(2*j*ηʹ);
var qʹ = 0;
for (var j=1; j<=6; j++) qʹ += 2*j*α[j] * Math.sin(2*j*ξʹ) * Math.sinh(2*j*ηʹ);

var γʹ = Math.atan(τʹ / Math.sqrt(1+τʹ*τʹ)*tanλ);
var γʺ = Math.atan2(qʹ, pʹ);

var γ = γʹ + γʺ;

// ---- scale: Karney 2011 Eq 25

var sinφ = Math.sin(φ);
var kʹ = Math.sqrt(1 - e*e*sinφ*sinφ) * Math.sqrt(1 + τ*τ) / Math.sqrt(τʹ*τʹ + cosλ*cosλ);
var kʺ = A / a * Math.sqrt(pʹ*pʹ + qʹ*qʹ);

var k = k0 * kʹ * kʺ;

// ------------

// shift x/y to false origins
x = x + falseEasting; // make x relative to false easting
if (y < 0) y = y + falseNorthing; // make y in southern hemisphere relative to false northing

// round to reasonable precision
x = Number(x.toFixed(6)); // nm precision
y = Number(y.toFixed(6)); // nm precision
var convergence = Number(γ.toDegrees().toFixed(9));
var scale = Number(k.toFixed(12));

var h = this.lat>=0 ? 'N' : 'S'; // hemisphere

return new Utm(zone, h, x, y, this.datum, convergence, scale);
};

This is my unfinished code. Its just a case of working through the routine above and wondering how code can be so difficult to read. Its all untested and I'm not certain my hyperbolic functions are correct.

Code: Select all

Function ConvertToUTM pLat, pLng
   
   # pLat pLon are numbers.  negative values indicate 
   -- 0-180W in Latitude or 
   -- 0-90S in Longitude 
   
   if pLat is not a number then return "bad value - latitude"
   if pLng is not a number then return "bad value - longitude"
   
   if pLat <= -80 OR pLat >= 84 then return "Outside UTM Limits - Latitude"
   if pLng < 180 OR pLat > 180 then return "Outside UTM Limits - Longitude"
   
   put 500000 into tFalseEasting 
   put 10000000 into tFalseNorthing
   
   
   put the floor of ((pLng+180)/6) + 1  into tZone-- longitude zone (sixty 6 degree zones)
   
   put ((tZone-1)*6)-180 + 3 into tLngCM // Longitude of centeral meridan of zone
   -- JavaScript code uses lambda0 for this variable
   
   //  ------ handle special cases Norway/Svalbard exceptions
   // grid zones are 8 degrees tall; 0 deg North is offset 10 into latitude bands array
   
   put "CDEFGHJKLMNPQRSTUVWXX" into tMgrsLatBands  -- X is repeated for posn's 80-84 degrees North
   put char (Floor((pLat/8) + 10 )) of tMgrsLatBands into tLatBand  -- single letter - N is equator to 8 deg north
   
   -- adjust tZone and central meridian for Norway special case
   
   if (tZone=31 AND tLatBand="V" AND pLon>=3) then 
      add 1 to tzone
      add 6 to tLngCM  -- still in degrees
   end if
   
   if (tZone=32 AND tLatBand="X" AND pLon<9) then 
      subtract 1 from tzone
      subtract 6 from tLngCM  -- still in degrees
   end if
   
   if (tZone=32 AND tLatBand="X" AND pLon>=9) then 
      add 1 to tzone
      add 6 to tLngCM  -- still in degrees
   end if
   
   if (tZone=34 AND tLatBand="X" AND pLon<21) then 
      subtract 1 from tzone
      subtract 6 from tLngCM  -- still in degrees
   end if
   
   if (tZone=34 AND tLatBand="X" AND pLon>=21) then 
      add 1 to tzone
      add 6 to tLngCM  -- still in degrees
   end if
   
   if (tZone=36 AND tLatBand="X" AND pLon<33) then 
      subtract 1 from tzone
      subtract 6 from tLngCM  -- still in degrees
   end if
   
   if (tZone=36 AND tLatBand="X" AND pLon>=33) then 
      add 1 to tzone
      add 6 to tLngCM  -- still in degrees
   end if
   
   // *****************************************
   
   put pLat*pi/180 into tPhi                -- Latitude in radians
   put (pLng-tLngCM)*pi/180 into tLambda  -- longitude ref the longitude of central meriadian in radians
   
   -- WGS84 datum conversion factors
   put  6378137 into a
   put 1/298.257223563 into f
   put 0.9996 into kZero -- UTM scale on the central meridian
   put sqrt(f*(2-f)) into e -- eccentricity
   put f / (2-f) into n -- 3rd flattening
   put n*n into n2
   put n*n2 into n3
   put n*n3 into n4
   put n*n4 into n5
   put n*n5 into n6
   
   put cos(tLambda) into cosLambda
   put sin(tLambda) into sinLambda
   put tan(tLambda) into tanLambda
   
   put tan(tPhi) into tTau
   put e*tTau/(sqrt(1+tTau*tTau))
   
   put sinh(e*atanh(e*tTau/sqrt(1+(tTau*tTau)))) into tDelta
   
   
end ConvertToUTM

Private Function aTanh pX
   --returns the inverse hyperbolic tangent
   put (ln((1+pX)/(1-pX)))/2 into tAnswer
   return tAnswer
   
end aTanh

Private Function Sinh pX
   put (exp(pX) - exp(-pX))/2 into tAnswer
   return tAnswer
end Sinh

Simon

[Simon Knight]

My first version will be based on the source Java Script function. I expect that the Livecode version will be slow because of all the calculations that are being conducted. However, my second version will move the setting of the variables that do not rely on the Lat Long passed to a set up function that stores to values to script local variables.