// ------------------------------------------------------------------------
// GeoPosition - hold positions in osgr, wgs84
// ------------------------------------------------------------------------
// .Based on code by Steve Loughran.
// Here : http://www.koders.com/cpp/fid179B65F9F342183C947E111C2AA0B6E6B3F3829A.aspx
// .Original OSGB to WGS84 c++ code written by  Chuck Gantz- chuck.gantz@globalstar.com
// Here: http://www.gpsy.com/gpsinfo/geotoutm/
// .This version Alistair Rutherford, July 2005. arutherford@netthreads.co.uk
// Here : www.gtraffic.info
// Here: www.netthreads.co.uk
// ------------------------------------------------------------------------
// Version - Date - Comment
// 1.0 - 23/07/05 - Initial version.
// 1.1 - 24/07/05 - Updated to return results in array
// ------------------------------------------------------------------------

function GeoPosition()
{
    // Math Constants
    this.PI = 3.14159265;
    this.FOURTHPI = this.PI / 4;
    this.DEG2RAD = this.PI / 180;
    this.RAD2DEG = 180.0 / this.PI;
    
    // Variables
    this.longitude = 0.0;
    this.latitude = 0.0;
    this.refEasting = 0.0;
    this.refNorthing = 0.0;
    
    this.osgbEasting = 0;
    this.osgbNorthing = 0;
    this.osgbGridSquare = "";
}

// ---------------------------------------------------------------------
// osgbToLL - Convert the OSGR36 grid ref value to WGS84 lat and long
//------------------------------------------------------------------------
function GeoPosition_osgbToLL(osgbZone, osbgEasting, osgbNorthing)
{
    // Converts OSGB coords to lat/long.  Equations from USGS Bulletin 1532 
    // East Longitudes are positive, West longitudes are negative. 
    // North latitudes are positive, South latitudes are negative
    // Lat and Long are in decimal degrees. 
    
    var k0 = 0.9996012717;
    var a = 0.0;
    var eccPrimeSquared = 0.0;
    var N1 = 0.0;
    var T1 = 0.0;
    var C1 = 0.0;
    var R1 = 0.0;
    var D = 0.0;
    var M = 0.0;
    var LongOrigin = -2.0;
    var LatOrigin = 49.0;
    
    var LatOriginRad = LatOrigin * this.DEG2RAD;
    
    var mu = 0.0;
    var phi1 = 0.0;
    var phi1Rad = 0.0;
    var x = 0.0;
    var y = 0.0;
    
    var majoraxis = a = 6377563.396; // Airy
    var minoraxis = 6356256.91; // Airy
    
    var eccSquared = (majoraxis * majoraxis - minoraxis * minoraxis) / (majoraxis * majoraxis);
    var e1 = (1.0-Math.sqrt(1.0-eccSquared))/(1.0+Math.sqrt(1.0-eccSquared));
    
    // Only calculate M0 once since it is based on the origin of the OSGB projection, which is fixed
    var M0 = a*((1 - eccSquared/4 - 3*eccSquared*eccSquared/64 - 5*eccSquared*eccSquared*eccSquared/256)*LatOriginRad - (3*eccSquared/8 + 3*eccSquared*eccSquared/32 + 45*eccSquared*eccSquared*eccSquared/1024)*Math.sin(2*LatOriginRad) + (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*Math.sin(4*LatOriginRad) - (35*eccSquared*eccSquared*eccSquared/3072)*Math.sin(6*LatOriginRad));
    
    // Calculate refEasting and refNorthing
    this.osgbSquareToRefCoords(osgbZone);
    
    x = osbgEasting - 400000.0 + this.refEasting; // remove 400,000 meter false easing for longitude
    y = osgbNorthing + 100000.0 + this.refNorthing; // remove 100,000 meter false easing for longitude
    
    eccPrimeSquared = (eccSquared)/(1-eccSquared);
    
    M = M0 + y / k0;
    mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/64-5*eccSquared*eccSquared*eccSquared/256));
    
    phi1Rad = mu + (3*e1/2-27*e1*e1*e1/32)*Math.sin(2*mu) + (21*e1*e1/16-55*e1*e1*e1*e1/32)*Math.sin(4*mu) +(151*e1*e1*e1/96)*Math.sin(6*mu);
            
    phi1 = phi1Rad*this.RAD2DEG;
    
    N1 = a/Math.sqrt(1-eccSquared*Math.sin(phi1Rad)*Math.sin(phi1Rad));
    T1 = Math.tan(phi1Rad)*Math.tan(phi1Rad);
    C1 = eccPrimeSquared*Math.cos(phi1Rad)*Math.cos(phi1Rad);
    R1 = a*(1-eccSquared)/Math.pow(1-eccSquared*Math.sin(phi1Rad)*Math.sin(phi1Rad), 1.5);
    D = x/(N1*k0);
    
    this.latitude = phi1Rad - (N1*Math.tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*eccPrimeSquared)*D*D*D*D/24 +(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C1*C1)*D*D*D*D*D*D/720);
    this.latitude = this.latitude * this.RAD2DEG;
    
    this.longitude = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccPrimeSquared+24*T1*T1) *D*D*D*D*D/120)/Math.cos(phi1Rad);
    this.longitude = LongOrigin + this.longitude * this.RAD2DEG;
    
    return new Array(this.latitude, this.longitude);
}

// ---------------------------------------------------------------------
// osgbSquareToRefCoords - Calculate the os grid ref square easting and
// northing offsets.
// ---------------------------------------------------------------------
function GeoPosition_osgbSquareToRefCoords(osgbGridSquare)
{
    var pos = 0;
    var x_multiplier = 0;
    var y_multiplier = 0;
    
    var gridSquare = "VWXYZQRSTULMNOPFGHJKABCDE";
    
    // find 500km offset
    var ch = osgbGridSquare[0];
    switch (ch)
    {
        case 'S':
            x_multiplier = 0;
            y_multiplier = 0;
        break;
        
        case 'T':
            x_multiplier = 1;
            y_multiplier = 0;
        break;

        case 'N':
            x_multiplier = 0;
            y_multiplier = 1;
        break;
    
        case 'O':
            x_multiplier = 1;
            y_multiplier = 1;
        break;

        case 'H':
            x_multiplier = 0;
            y_multiplier = 2;
        break;

        case 'J':
            x_multiplier = 1;
            y_multiplier = 2;
        break;
    }
    
    this.refEasting = x_multiplier * 500000;
    this.refNorthing = y_multiplier * 500000;
    
    // find 100km offset and add to 500km offset to get coordinate of 
    // square point is in
    var index = gridSquare.indexOf(osgbGridSquare[1])
    
    this.refEasting += ((index % 5) * 100000);
    this.refNorthing += (trunc(index / 5) * 100000);
}
 
// ---------------------------------------------------------------------
// llToOSGB - convert WGS84 to OSGB
// --------------------------------------------------------------------- 
function GeoPosition_llToOSGB(lat, lng)
{
    // Converts lat/long to OSGB coords.  Equations from USGS Bulletin 1532 
    // East Longitudes are positive, West longitudes are negative. 
    // North latitudes are positive, South latitudes are negative
    // Lat and Long are in decimal degrees
    // Written by Chuck Gantz- chuck.gantz@globalstar.com
    
    var a;
    var eccSquared;
    var k0 = 0.9996012717;
    
    var LongOrigin = -2;
    var LongOriginRad = LongOrigin * this.DEG2RAD;
    var LatOrigin = 49;
    var LatOriginRad = LatOrigin * this.DEG2RAD;
    var eccPrimeSquared;
    var N, T, C, A, M;
    
    var LatRad =lat*this.DEG2RAD;
    var LongRad = lng*this.DEG2RAD;
    var easting, northing;
    
    var majoraxis = a = 6377563.396;//Airy
    var minoraxis = 6356256.91;//Airy
    
    eccSquared = (majoraxis * majoraxis - minoraxis * minoraxis) / (majoraxis * majoraxis);
    
    // only calculate M0 once since it is based on the origin 
    // of the OSGB projection, which is fixed
    var M0 = a*((1 - eccSquared/4 - 3*eccSquared*eccSquared/64 - 5*eccSquared*eccSquared*eccSquared/256)*LatOriginRad - (3*eccSquared/8+ 3*eccSquared*eccSquared/32+ 45*eccSquared*eccSquared*eccSquared/1024)*Math.sin(2*LatOriginRad) + (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*Math.sin(4*LatOriginRad) - (35*eccSquared*eccSquared*eccSquared/3072)*Math.sin(6*LatOriginRad));
    
    var eccPrimeSquared = (eccSquared)/(1-eccSquared);
    
    N = a/Math.sqrt(1-eccSquared*Math.sin(LatRad)*Math.sin(LatRad));
    T = Math.tan(LatRad)*Math.tan(LatRad);
    C = eccPrimeSquared*Math.cos(LatRad)*Math.cos(LatRad);
    A = Math.cos(LatRad)*(LongRad-LongOriginRad);
    
    M = a*((1 - eccSquared/4 - 3*eccSquared*eccSquared/64 - 5*eccSquared*eccSquared*eccSquared/256)*LatRad - (3*eccSquared/8 + 3*eccSquared*eccSquared/32 + 45*eccSquared*eccSquared*eccSquared/1024)*Math.sin(2*LatRad) + (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*Math.sin(4*LatRad) - (35*eccSquared*eccSquared*eccSquared/3072)*Math.sin(6*LatRad));
    
    easting = (k0*N*(A+(1-T+C)*A*A*A/6 + (5-18*T+T*T+72*C-58*eccPrimeSquared)*A*A*A*A*A/120));
    easting += 400000.0; //false easting
    
    northing = (k0*(M-M0+N*Math.tan(LatRad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24 + (61-58*T+T*T+600*C-330*eccPrimeSquared)*A*A*A*A*A*A/720)));
    northing -= 100000.0;//false northing
    
    this.coordsToOSGBSquare(easting, northing);
    
    return new Array(this.osgbGridSquare, this.osgbEasting, this.osgbNorthing);
}

function GeoCode_coordsToOSGBSquare(easting, northing)
{
    var GridSquare = "VWXYZQRSTULMNOPFGHJKABCDE";
    var posx; //positions in grid
    var posy; 
    
    this.osgbEasting = (easting + 0.5); //round to nearest int
    this.osgbNorthing = (northing + 0.5); //round to nearest int
    
    //find correct 500km square
    posx = this.osgbEasting / 500000;
    posy = this.osgbNorthing / 500000;
    this.osgbGridSquare += GridSquare[trunc(posx) + trunc(posy) * 5 + 7];
    
    //find correct 100km square
    posx = this.osgbEasting % 500000; //remove 500 km square
    posy = this.osgbNorthing % 500000; //remove 500 km square
    posx = posx / 100000; //find 100 km square
    posy = posy / 100000; //find 100 km square
    this.osgbGridSquare += GridSquare[trunc(posx) + trunc(posy) * 5];
    
    //remainder is northing and easting
    this.osgbNorthing = this.osgbNorthing % 500000; 
    this.osgbNorthing = trunc(this.osgbNorthing % 100000);
    
    this.osgbEasting = this.osgbEasting % 500000;
    this.osgbEasting = trunc(this.osgbEasting % 100000);
}
 
function trunc(val)
{
    if (val < 0)
        return -(Math.floor(Math.abs(val)));
    else
        return Math.floor(Math.abs(val));
} 
		
GeoPosition.prototype.osgbToLL = GeoPosition_osgbToLL;
GeoPosition.prototype.osgbSquareToRefCoords = GeoPosition_osgbSquareToRefCoords;
GeoPosition.prototype.llToOSGB = GeoPosition_llToOSGB;
GeoPosition.prototype.coordsToOSGBSquare = GeoCode_coordsToOSGBSquare;