/****************************************************
 *
 *  Javascript classes for Geospatial coordinates and Coordinate Systems
 *
 *
 *  <c) 2007 Martin Davis
 *
 ****************************************************/

//===============================
// DMS class
//===============================

function DegMinSec()
{
  this.deg = 0.0;
  this.min = 0.0;
  this.sec = 0.0;
}

DegMinSec.prototype.setDegMinSec = function(d, m, s)
{
  this.deg = d;
  this.min = m;
  this.sec = s;
}


DegMinSec.prototype.getDecDegree = function()
{
  dd = Math.abs(this.deg) + this.min / 60 + this.sec / 3600;
  if (this.deg < 0) {
    dd = -dd;
  }
  return dd;
}

DegMinSec.prototype.setDecDegree = function(inputDD)
{
  dd = inputDD;
  if (inputDD < 0) {
    dd = -dd;
  }
   
  this.deg = Math.floor(dd);
  minDec = 60 * (dd - this.deg);
  this.min = Math.floor(minDec);
  secDec = 60 * (minDec - this.min);
  this.sec = secDec;

  if (inputDD < 0) {
    this.deg = -this.deg;
  }
}

//===============================
// Degrees/Radians conversion functions
//===============================


function degreesToRad(deg)
{
  return  deg / 180.0 * Math.PI;
}

function radToDegrees(deg)
{
  return  deg * 180.0 / Math.PI;
}


//===============================
// Point class
//===============================

function Point()
{
  this.x = 0.0;
  this.y = 0.0;
  this.zone = 0.0;
}

//===============================
// Speroid class
//===============================

function Spheroid(ia, ib, irf)
{
  this.a = ia;
  this.b = ib;
  this.rf = irf;

  if (this.b > 1.0) {
    this.f = 1.0 - (this.b / this.a);
    this.rf = 1.0 / this.f;
  }
  else {
    this.f = 1.0 / this.rf;
    this.b = this.a - this.a * this.f;
  }
  this.es = this.f + this.f - this.f * this.f;
  this.e = Math.sqrt(this.es);
}
	
Spheroid.prototype.getA = function()
{
return this.a; 
}

Spheroid.prototype.getB = function()
{
return this.b; 
}

Spheroid.prototype.getE = function()
{
return this.e; 
}

Spheroid.prototype.meridianRadiusOfCurvature = function(latitude)
{
    er = 1.0 - this.es * Math.sin(latitude) * Math.sin(latitude);
    el = Math.pow(er, 1.5);
    M0 = (a * (1.0 - this.es)) / el;
    return M0;
}

Spheroid.prototype.primeVerticalRadiusOfCurvature = function(latitude) 
{
    T1 = a * a;
    T2 = T1 * Math.cos(latitude) * Math.cos(latitude);
    T3 = b * b * Math.sin(latitude) * Math.sin(latitude);
    N0 = T1 / Math.sqrt(T2 + T3);
    return N0;
}

//===============================
// MeridianArcLength class
//===============================

function MeridianArcLength() 
{
  this.s = 0;
  this.a0 = 0;
  this.a2 = 0;
  this.a4 = 0;
  this.a6 = 0;
  this.a8 = 0;
}


MeridianArcLength.prototype.compute = function(spheroid, lat, diff) 
{
//  Returns the meridian arc length given the latitude
    a = spheroid.getA();
    e = spheroid.getE();
    e2 = e * e;
    e4 = e2 * e2;
    e6 = e4 * e2;
    e8 = e4 * e4;
    this.a0 = 1.0 - e2 / 4.0 - 3.0 * e4 / 64.0 - 5.0 * e6 / 256.0 - 175.0 * e8 / 16384.0;
    this.a2 = 3.0 / 8.0 * (e2 + e4 / 4.0 + 15.0 * e6 / 128.0 - 455.0 * e8 / 4096.0);
    this.a4 = 15.0 / 256.0 * (e4 + 3.0 * e6 / 4.0 - 77.0 * e8 / 128.0);
    this.a6 = 35.0 / 3072.0 * (e6 - 41.0 * e8 / 32.0);
    this.a8 = -315.0 * e8 / 131072.0;
    if (diff == 0) {
      this.s = a * (this.a0 * lat - this.a2 * Math.sin(2.0 * lat) + this.a4 * Math.sin(4.0 * lat)
           - this.a6 * Math.sin(6.0 * lat) + this.a8 * Math.sin(8.0 * lat));
    }
    else {
      this.s = this.a0 * lat - 2.0 * this.a2 * Math.cos(2.0 * lat) + 4.0 * this.a4 * Math.cos(4.0 * lat)
           - 6.0 * this.a6 * Math.cos(6.0 * lat) + 8.0 * this.a8 * Math.cos(8.0 * lat);
    }
}

//===============================
// Albers class
//===============================

function Albers(sph)
{
  this.sph = sph;
}

Albers.prototype.setParameters = function(centralMeridian,
      firstStandardParallel,
      secondStandardParallel,
      latitudeOfProjection,
      falseEasting,
      falseNorthing)
{
  this.L0 = centralMeridian * Math.PI / 180.0;
  this.phi1 = firstStandardParallel * Math.PI / 180.0;
  this.phi2 = secondStandardParallel * Math.PI / 180.0;
  this.phi0 = latitudeOfProjection * Math.PI / 180.0;
  this.X0 = falseEasting;
  this.Y0 = falseNorthing;

  m1 = this.albersM(this.phi1);
  m2 = this.albersM(this.phi2);
  q1 = this.albersQ(this.phi1);
  q2 = this.albersQ(this.phi2);
  q0 = this.albersQ(this.phi0);

  this.A_n = (m1 * m1 - m2 * m2) / (q2 - q1);
  this.A_C = m1 * m1 + this.A_n * q1;
  a = this.sph.getA();
  this.A_p0 = (a * Math.sqrt(this.A_C - this.A_n * q0)) / this.A_n;

  this.ptGeoRad = new Point();
}

/*
 * Computes the forward transformation into Albers.
 * 
 * Parameters:
 *   ptGeo - input geographic (lon, lat) coordinate, in degrees
 *   ptXY - output Albers coordinate, in metres
 * Returns:
 *   void
 */
Albers.prototype.forward = function(ptGeo, ptXY) 
{
  latRad = ptGeo.y / 180.0 * Math.PI;
  lonRad = ptGeo.x / 180.0 * Math.PI;
    a = this.sph.getA();
    que = this.albersQ(latRad);
    theta = this.A_n * (lonRad - this.L0);
    pee = (a * Math.sqrt(this.A_C - this.A_n * que)) / this.A_n;
    ptXY.x = pee * Math.sin(theta) + this.X0;
    ptXY.y = this.A_p0 - pee * Math.cos(theta) + this.Y0;
  }

/*
 * Computes the inverse transformation into Geographics.
 * 
 * Parameters:
 *   ptXY - input Albers coordinate, in metres
 *   ptGeo - output geographic (lon, lat) coordinate, in degrees
 * Returns:
 *   void
 */
Albers.prototype.inverse = function(ptXY, ptGeo) 
{
    a = this.sph.getA();
    e = this.sph.getE();
    es = e * e;
    x = ptXY.x - this.X0;
    y = ptXY.y - this.Y0;
    theta = Math.atan2(x, this.A_p0 - y);
    pee = Math.sqrt(x * x + Math.pow((this.A_p0 - y), 2.0));
    que = (this.A_C - (pee * pee * this.A_n * this.A_n) / (a * a)) / this.A_n;
    lon = this.L0 + theta / this.A_n;

    li = Math.asin(que / 2.0);
    delta = 10e010;
    do {
      j1 = Math.pow((1.0 - es * Math.pow(Math.sin(li), 2.0)), 2.0) / (2.0 * Math.cos(li));
      k1 = que / (1.0 - es);
      k2 = Math.sin(li) / (1.0 - es * Math.pow(Math.sin(li), 2.0));
      k3 = (1.0 / (2.0 * e)) * Math.log((1.0 - e * Math.sin(li)) / (1.0 + e * Math.sin(li)));
      lip1 = li + j1 * (k1 - k2 + k3);
      delta = Math.abs(lip1 - li);
      li = lip1;
    } while (delta > 1.0e-012);
    lat = li;

  ptGeo.x = radToDegrees(lon);
  ptGeo.y = radToDegrees(lat);

  }

Albers.prototype.albersQ = function(lat) 
{
    e = this.sph.getE();
    q = (1.0 - e * e) * (Math.sin(lat) / (1.0 - e * e * Math.pow(Math.sin(lat), 2.0))
         - (1.0 / (2.0 * e)) * Math.log((1.0 - e * Math.sin(lat)) / (1.0 + e * Math.sin(lat))));
    return q;
}

Albers.prototype.albersM = function(lat) 
{
    e = this.sph.getE();
    m = Math.cos(lat) / Math.sqrt(1.0 - e * e * Math.pow(Math.sin(lat), 2.0));
    return m;
}


//===============================
// Transverse Mercator class
//===============================

function TransverseMercator(sph)
{
  this.sph = sph;
}


  /**
   *@param  centralMeridian  in degrees
   */
TransverseMercator.prototype.setParameters = function(centralMeridian) 
{
    this.L0 = degreesToRad(centralMeridian);
}


/*
 * Computes the inverse transformation into Geographics.
 * 
 * Parameters:
 *   ptXY - input TM coordinate, in metres
 *   ptGeo - output geographic (lon, lat) coordinate, in degrees
 * Returns:
 *   void
 */
TransverseMercator.prototype.inverse = function(ptXY, ptGeo) 
{
    L1 = this.footPointLatitude(p.y);
    a = this.sph.getA();
    b = this.sph.getB();
    ep2 = (a * a - b * b) / (b * b);
    
    // N1 = the radius of curvature of the spheroid in the prime vertical plane
    // at the foot point latitude
    N1 = this.sph.primeVerticalRadiusOfCurvature(L1);
    // M1 = meridian radius of curvature at the foot point latitude
    M1 = this.sph.meridianRadiusOfCurvature(L1);
    n12 = ep2 * Math.pow(Math.cos(L1), 2.0);
    n1 = Math.sqrt(n12);
    n14 = n12 * n12;
    n16 = n14 * n12;
    n18 = n14 * n14;

    t1 = Math.tan(L1);
    t12 = t1 * t1;
    t14 = t12 * t12;
    t16 = t14 * t12;

    u0 = t1 * Math.pow(p.x, 2.0) / (2.0 * M1 * N1);
    u1 = t1 * Math.pow(p.x, 4.0) / (24.0 * M1 * Math.pow(N1, 3.0));
    u2 = t1 * Math.pow(p.x, 6.0) / (720.0 * M1 * Math.pow(N1, 5.0));
    u3 = t1 * Math.pow(p.x, 8.0) / (40320.0 * M1 * Math.pow(N1, 7.0));
    v1 = 5.0 + 3.0 * t12 + n12 - 4.0 * n14 - 9.0 * n12 * t12;
    v2 = 61.0 - 90.0 * t12 + 46.0 * n12 + 45.0 * t14 - 252.0 * t12 * n12 - 3.0 * n14
         + 100.0 * n16 - 66.0 * t12 * n14 - 90.0 * t14 * n12 + 88.0 * n18 + 225.0 * t14 * n14
         + 84.0 * t12 * n16 - 192.0 * t12 * n18;
    v3 = 1385.0 + 3633.0 * t12 + 4095.0 * t14 + 1575.0 * t16;

    lat = L1 - u0 + u1 * v1 - u2 * v2 + u3 * v3;

    XdN1 = p.x / N1;
    u0 = XdN1;
    u1 = Math.pow(XdN1, 3.0) / 6.0;
    u2 = Math.pow(XdN1, 5.0) / 120.0;
    u3 = Math.pow(XdN1, 7.0) / 5040.0;
    v1 = 1.0 + 2.0 * t12 + n12;
    v2 = 5.0 + 6.0 * n12 + 28.0 * t12 - 3.0 * n14 + 8.0 * t12 * n12 + 24.0 * t14 - 4.0 * n16 + 4.0 * t12 * n14 + 24.0 * t12 * n16;
    v3 = 61.0 + 662.0 * t12 + 1320.0 * t14 + 720.0 * t16;
    lon = 1.0 / Math.cos(L1) * (u0 - u1 * v1 + u2 * v2 - u3 * v3) + this.L0;
    
    ptGeo.x = radToDegrees(lon);
    ptGeo.y = radToDegrees(lat);

}

/*
 * Computes the forward transformation into Transverse Mercator.
 * 
 * Parameters:
 *   ptGeo - input geographic (lon, lat) coordinate, in degrees
 *   ptXY - output Albers coordinate, in metres
 * Returns:
 *   void
 */
TransverseMercator.prototype.forward = function(ptGeo, ptXY) 
{
  lonRad = degreesToRad(ptGeo.x);
  latRad = degreesToRad(ptGeo.y);
  
    a = this.sph.getA();
    b = this.sph.getB();

    // ep2 = the second eccentricity squared.
    ep2 = (a * a - b * b) / (b * b);
    // N = the radius of curvature of the spheroid in the prime vertical plane
    N = this.sph.primeVerticalRadiusOfCurvature(latRad);

    n2 = ep2 * Math.pow(Math.cos(latRad), 2.0);
    n = Math.sqrt(n2);
    n4 = n2 * n2;
    n6 = n4 * n2;
    n8 = n4 * n4;

    t = Math.tan(latRad);
    t2 = t * t;
    t4 = t2 * t2;
    t6 = t4 * t2;
    S = new MeridianArcLength();
    S.compute(this.sph, latRad, 0);

    cosLat = Math.cos(latRad);
    sinLat = Math.sin(latRad);

    L = lonRad - this.L0;// 'L' for lambda (longitude) - must be in radians
    L2 = L * L;
    L3 = L2 * L;
    L4 = L2 * L2;
    L5 = L4 * L;
    L6 = L4 * L2;
    L7 = L5 * L2;
    L8 = L4 * L4;

    u0 = L * cosLat;
    u1 = L3 * Math.pow(cosLat, 3.0) / 6.0;
    u2 = L5 * Math.pow(cosLat, 5.0) / 120.0;
    u3 = L7 * Math.pow(cosLat, 7.0) / 5040.0;
    v1 = 1.0 - t2 + n2;
    v2 = 5.0 - 18.0 * t2 + t4 + 14.0 * n2 - 58.0 * t2 * n2 + 13.0 * n4 + 4.0 * n6 - 64.0 * n4 * t2 - 24.0 * n6 * t2;
    v3 = 61.0 - 479.0 * t2 + 179.0 * t4 - t6;
    x = u0 + u1 * v1 + u2 * v2 + u3 * v3;

    u0 = L2 / 2.0 * sinLat * cosLat;
    u1 = L4 / 24.0 * sinLat * Math.pow(cosLat, 3.0);
    u2 = L6 / 720.0 * sinLat * Math.pow(cosLat, 5.0);
    u3 = L8 / 40320.0 * sinLat * Math.pow(cosLat, 7.0);
    v1 = 5.0 - t2 + 9.0 * n2 + 4.0 * n4;
    v2 = 61.0 - 58.0 * t2 + t4 + 270.0 * n2 - 330.0 * t2 * n2 + 445.0 * n4 + 324.0 * n6 - 680.0 * n4 * t2
         + 88.0 * n8 - 600.0 * n6 * t2 - 192.0 * n8 * t2;
    v3 = 1385.0 - 311.0 * t2 + 543.0 * t4 - t6;
    y = S.s / N + u0 + u1 * v1 + u2 * v2 + u3 * v3;

    ptXY.x = N * x;
    ptXY.y = N * y;
}

TransverseMercator.prototype.footPointLatitude = function(y) 
{
// returns the footpoint Latitude given the y coordinate

    a = this.sph.getA();
    newlat = y / a;
    i = 0;
    S = new MeridianArcLength();
    do {
      Lat1 = newlat;
      i++;
      if (i == 100) {
          //Prevent infinite loop. 
          break; 
      }
      S.compute(this.sph, Lat1, 0);
      flat = S.s - y;
      dflat = a * (S.a0 - 2.0 * S.a2 * Math.cos(2.0 * Lat1) + 4.0 * S.a4 * Math.cos(4.0 * Lat1)
           - 6.0 * S.a6 * Math.cos(6.0 * Lat1) + 8.0 * S.a8 * Math.cos(8.0 * Lat1));
      newlat = Lat1 - flat / dflat;
      //Increased tolerance from 1E-16 to 1E-15. 1E-16 was causing an infinite loop.
    } while (Math.abs(newlat - Lat1) > 1.0e-015);
    Lat1 = newlat;
    return Lat1;
}

//===============================
// Universal Transverse Mercator class
//===============================

function UTM(sph)
{
  this.sph = sph;
  this.SCALE_FACTOR = 0.9996;
  this.FALSE_EASTING = 500000.0;
  this.FALSE_NORTHING = 0.0;

}

/**
 *@param  zone  the UTM zone
 */
UTM.prototype.setParameters = function(zone) 
{
  this.tm = new TransverseMercator(this.sph);
  this.tm.setParameters(this.centMer(zone));
}

UTM.prototype.forward = function(ptGeo, ptXY) 
{
  this.tm.forward(ptGeo, ptXY);
  ptXY.x = this.SCALE_FACTOR * ptXY.x + this.FALSE_EASTING;
  ptXY.y = this.SCALE_FACTOR * ptXY.y + this.FALSE_NORTHING;
}

UTM.prototype.inverse = function(ptXY, ptGeo) 
{
  p = new Point();	
  p.x = (ptXY.x - this.FALSE_EASTING) / this.SCALE_FACTOR;
  p.y = (ptXY.y - this.FALSE_NORTHING) / this.SCALE_FACTOR;
  this.tm.inverse(p, ptGeo);
}

UTM.prototype.centMer = function(zone) 
{
  return -180 - 3 + 6 * zone;
}

UTM.prototype.zone = function(lon) 
{
  return Math.floor((lon + 180) / 6) + 1;
}