/*
	original code by 
	Paul Hastings paul@sustainableGIS.com
	
	29-Aug-2005:
	ported to ActionScript by
	Tobias Houfek ibot@snafu.de, zillion@virtual-jabog32.de
*/


// some static vars/required values
iLONGITUDE = 1-1; // cf arrays start at 1... 
iLATITUDE = 2-1;  // <tobi> :) Flash Arrays start 0
iEASTING = 1-1;
iNORTHING = 2-1;


//WGS84
semiMajorAxis = 6378137.0;
inverseFlattening = 298.257223563;

a = semiMajorAxis;
f = 1.00 / inverseFlattening;

// UTM Parameters
longitude1 = -177; // Longitude of initial central meridian (Zone one) = 177 degrees west longitude
zoneWidth = 6; // Zone width = 6 degrees

centralScaleFactor = 0.9996;
k0 = centralScaleFactor;
falseEasting = 500000; // False easting = 500,000m
falseNorthing = 0; // False northing in the northern hemisphere, southern = 10,000,000m 

e2 = (2 * f) - (f * f); // Eccentricity
e = Math.sqrt(e2);
  
longitude0West = longitude1 - (zoneWidth * 1.5);
longitude0 = longitude0West + (zoneWidth * 0.5);

projectionName="UTM";
ellipsoidName="WGS84";


	errMessages = new Object();
	errMessages.NOT_ENOUGH_POINTS_MESSAGE = "ProjectionException: Less than 2 coordinates supplied";
	errMessages.INVALID_PROJECTION_NAME_MESAGE = "ProjectionException: unknown projection";
	errMessages.UTM_ZONE_NOT_DEFINED_MESSAGE = "ProjectionException: UTM zone not defined";
	errMessages.LONGITUDE_OUT_OF_RANGE_MESSAGE = "ProjectionException: longitude out of range";
	errMessages.LATITUDE_OUT_OF_RANGE_MESSAGE = "ProjectionException: latitude out of range";


function pow(a,b){
	return a^b;
}//pow

function sec(value) {
    return 1.0 / Math.cos(value);
} //sec

function toDegrees(radians){
		return (radians / Math.PI) * 180.0;
} //toDegrees

function toRadians(degrees) {
  	return (degrees / 180.0) * Math.PI;
} // toRadians

function getProjectionName(){
	var name="";
	if (UTMzone == undefined) 
		trace(errMessages.UTM_ZONE_NOT_DEFINED_MESSAGE);
	else {
		name=projectionName+" (Zone"+utmZone+")";
		return name;
	}
}//getProjectionName

function getEllipsoidName(){
	return ellipsoidName;
}//getEllipsoidName

function setEllipsoid(ellipsoid){

	return EllipsoidName;
}//setEllipsoid

function getUTMzone(){
	return UTMzone;
} //getUTMzone

// used for geographic to UTM, calculates somethings about UTM
function setUTMZone(longitude, latitude) {
	// make these known to app
	var thisZone=0;
	var hemisphere="N"; // default hemisphere
	if (Math.abs(longitude) > 180)
		trace(errMessages.LONGITUDE_OUT_OF_RANGE_MESSAGE+": "+longitude);
	if (Math.abs(latitude) > 90)
		trace(errMessages.LATITUDE_OUT_OF_RANGE_MESSAGE+": "+latitude);
	thisZone = Math.round((longitude+180.0)/zoneWidth+1.0);
	centralMeridian = thisZone * zoneWidth + longitude0;
	if (latitude < 0) {
		hemisphere = "S";
		falseNorthing = 10000000;
	}	
	utmZone = thisZone + hemisphere;
  //initErrorMessages();
} //setUTMzone

// used for UTM to geographic, need to know somethings about UTM
function initUTM(zone) {
	// make these known to app
	var hemisphere="N"; // default hemispehere
	var thisZone=0;
	if (zone.substr(-1,1) == "S") {
		hemisphere="S";
		falseNorthing = 10000000;
	}	
	
	//------------------------Hope it works out------
	//<tobi> OMG !! :|         original:
	//>>>>>thisZone=rereplace(zone,"N|S","","ALL");<<<<<<<<<
	//Flashs string methods are really limited
	if(zone.indexOf("N") == true){
				var zoneI=thisZone.indexOf("N");
  		}else if(zone.indexOf("S") == true){
	  			var zoneI=thisZone.indexOf("S");
  	}
		thisZone = zone.slice(0,zoneI)
	//-----------------------------------------------
	
	centralMeridian = thisZone * zoneWidth + longitude0;
	utmZone=thisZone & hemisphere;
  //initErrorMessages();
} //getUTMzone

function getEasting() {
	var easting=0.0;
	term1 = (omega[3] / 6.0) * c[3] * (psi[2] - t[3]);
	term2 = (omega[5] / 120.0) * c[5] * ( 4 * psi[4] * (1.0 - 6.0 * t[3]) + psi[3] * (1.0 + 8.0 * t[3]) - psi[2] * 2.0 * t[3] + t[5] );
	term3 = (omega[7] / 5040.0) * c[7] * (61.0 - 479.0 * t[5] - t[7]);
	eDash = (k0 * nu * omega[2] * c[2]) * (1 + term1 + term2 + term3);
	easting = eDash + falseEasting;
    return easting;
} //geteasting

function getNorthing() {
	var northing =0.0;
	term1 = (omega[3] / 2.0) * nu * s * c[2];
	term2 = (omega[5] / 24.0) * nu * s * c[4] * (4.0 * psi[3] + psi[2] - t[3]);
	term3 = (omega[7] / 720.0) * nu * s * c[6] * ( 8.0 * psi[5] * (11.0 - 24.0 * t[3]) - 28.0 * psi[4] * (1.0 - 6.0 * t[3]) + psi[3] * (1.0 - 32.0 * t[3]) - psi[2] * (2.0 * t[3]) + t[5] );
	term4 = (omega[9] / 40320.0) * nu * s * c[8] * (1385.0 - 3111.0 * t[3] + 543.0 * t[5] - t[7]);
	nDash = k0 * (m + term1 + term2 + term3 + term4);
	northing = nDash + falseNorthing;
	return northing;
}//getnorthing

function getLatitude() {
	var latitude=0.0;
	term1 = (tDash[2] / (k0 * rhoDash)) * (eDash * x[2] / 2.0);
	term2 = (tDash[2] / (k0 * rhoDash)) * (eDash * x[4] / 24.0) * ( -4.0 * psiDash[3] + 9.0 * psiDash[2] * (1.0 - tDash[3]) + 12.0 * tDash[3] );
	term3 = (tDash[2] / (k0 * rhoDash)) * (eDash * x[6] / 720.0) * ( 8.0 * psiDash[5] * (11.0 - 24.0 * tDash[3]) - 12.0 * psiDash[4] * (21.0 - 71.0 * tDash[3]) + 15.0 * psiDash[3] * (15.0 - 98.0 * tDash[3] + 15.0 * tDash[5]) + 180.0 * psiDash[2] * (5.0 * tDash[3] - 3.0 * tDash[5]) + 360.0 * tDash[5] );
	term4 = (tDash[2] / (k0 * rhoDash)) * (eDash * x[8] / 40320.0) * (1385.0 + 3633.0 * tDash[3] + 4095.0 * tDash[5] + 1575.0 * tDash[7]);
	latitude = phiDash - term1 + term2 - term3 + term4;
	return toDegrees(latitude);
} //getLatitude
   
function getLongitude() {
	var longitude=0.0;
	var omegaL=0.0; // "global" omega is array
	if (centralMeridian == undefined) 
		trace(errMessages.UTM_ZONE_NOT_DEFINED_MESSAGE);
	term1 = x[2] * secDash;
	term2 = (x[4] / 6.0) * secDash * (psiDash[2] + 2.0 * tDash[3]);
	term3 = (x[6] / 120.0) * secDash * ( -4.0 * psiDash[4] * (1.0 - 6.0  * tDash[3]) + psiDash[3] * (9.0 - 68.0 * tDash[3]) + 72.0 * psiDash[2] * tDash[3] + 24.0 * tDash[5] );
	term4 = (x[8] / 5040.0) * secDash * (61.0 + 662.0 * tDash[3] + 1320.0 * tDash[5] + 720.0 * tDash[7]);
	omegaL = term1 - term2 + term3 - term4;
	longitude = toRadians(centralMeridian) + omegaL;
	return toDegrees(longitude);
} //getLongitude

function meridianDistance(phi) {
    var A0 = 1.0 - (e2 / 4.0) - (3.0 * e2*e2 / 64.0) - (5.0 * e2*e2*e2 / 256.0);
    var A2 = (3.0/8.0) * (e2 + e2*e2 / 4.0 + 15.0 * e2*e2*e2 / 128.0);
    var A4 = (15.0/256.0) * (e2*e2 + 3.0 * e2*e2*e2 / 4.0);
    var A6 = 35.0 * e2*e2*e2 / 3072.0;
	var mM = 0.0;
	mM = a * ( (A0 * phi) - A2 * Math.sin(2.0 * phi) + A4 * Math.sin(4.0 * phi) - A6 * Math.sin(6.0 * phi) );
    return mM;
} //meridianDistance

function footPointLatitude(mM) {
    var n = f / (2.0 - f);
    var n2 = n*n;
    var n3 = n2*n;
    var n4 = n3*n;
    var phiDash=0.0;
    var G = a * (1.0 - n) * (1.0 - n2) * (1.0 + (9.0/4.0) * n2 + (225.0/64.0) * n4) * (Math.PI/180.0);
    var sigma = (mM * Math.PI) / (180.0 * G);
    phiDash = sigma + ((3.0 * n / 2.0) - (27.0 * n3 / 32.0)) * Math.sin(2.0 * sigma) + ((21.0 * n2 / 16.0) - (55.0 * n4 / 32.0)) * Math.sin(4.0 * sigma) + (151.0 * n3 / 96.0) * Math.sin(6.0 * sigma) + (1097.0 * n4 / 512.0) * Math.sin(8.0 * sigma);
    return phiDash;
} //footPointLatitude

function RoC(s,returnWhat) {
	var rho = a * (1.0 - e2) / pow(1.0 - (e2 * s), 1.5);
	var nu = a / pow((1.0 - e2 * s), 0.5);
	var psi = nu / rho;

    if (returnWhat == "rho") return rho;
	if (returnWhat == "Nu") return Nu;
	if (returnWhat == "Psi") return psi;
} //RoC 

function GeographicalToGrid(latitude, longitude) {
	var latRadians = toRadians(latitude);
	if (centralMeridian == undefined) 
		trace(errMessages.UTM_ZONE_NOT_DEFINED_MESSAGE);
	 // calculate cos^n(latitude) for n in (1..7)
	
	c = new Array;
	c[1] = 1.0;
	c[2] = Math.cos(latRadians);
	c[3] = c[2] * c[2];
	c[4] = c[3] * c[2];
	c[5] = c[4] * c[2];
	c[6] = c[5] * c[2];
	c[7] = c[6] * c[2];
	c[8] = c[7] * c[2];
	 
	t = new Array;
	t[1] = 1.0;
	t[2] = Math.tan(latRadians);
	t[3] = t[2] * t[2];
	t[4] = t[3] * t[2];
	t[5] = t[4] * t[2];
	t[6] = t[5] * t[2];
	t[7] = t[6] * t[2];
	
	omega = new Array;
	omega[1] = 1.0;
	omega[2] = toRadians(longitude - centralMeridian);
	omega[3] = omega[2] * omega[2];      
	omega[4] = omega[3] * omega[2];
	omega[5] = omega[4] * omega[2];
	omega[6] = omega[5] * omega[2];      
	omega[7] = omega[6] * omega[2];      
	omega[8] = omega[7] * omega[2];      
	omega[9] = omega[8] * omega[2];      
	 
	 s = Math.sin(latRadians);
	
	 nu = RoC(s*s,"Nu");
	 
	 psi = new Array;
	 psi[1] = 1.0;
	 psi[2] = RoC(s*s,"Psi");
	 psi[3] = psi[2] * psi[2];
	 psi[4] = psi[3] * psi[2];
	 psi[5] = psi[4] * psi[2];
	 
	 m = meridianDistance(latRadians);
} //GeographicalToGrid

function GridToGeographical(easting,northing)  {
	nDash = northing - falseNorthing;
    m = nDash / k0;
    phiDash = footPointLatitude(m);
	tDash = arrayNew(1);
	tDash[1] = 1.0;
    tDash[2] = Math.tan(phiDash);
    tDash[3] = tDash[2] * tDash[2];
    tDash[4] = tDash[3] * tDash[2];
    tDash[5] = tDash[4] * tDash[2];
    tDash[6] = tDash[5] * tDash[2];
    tDash[7] = tDash[6] * tDash[2];
	sDash = Math.sin(phiDash);
    rhoDash = RoC(sDash*sDash,"rho");
    nuDash = RoC(sDash*sDash,"Nu");
    psiDash = new Array;
    psiDash[1] = 1.0;
    psiDash[2] = RoC(sDash*sDash,"Psi");
    psiDash[3] = psiDash[2] * psiDash[2];
    psiDash[4] = psiDash[3] * psiDash[3];
    psiDash[5] = psiDash[4] * psiDash[2];
    eDash = easting - falseEasting;
	
	x = new Array;
    x[1] = 1.0;
    x[2] = eDash / (k0 * nuDash);
    x[3] = x[2] * x[2];
    x[4] = x[3] * x[2];
    x[5] = x[4] * x[2];
    x[6] = x[5] * x[3];
	x[7] = x[6] * x[2];
    x[8] = x[7] * x[2];
	secDash = sec(phiDash);
} // GridToGeographical

function getGridConvergence() {
	var gamma = 0.0;
	term1 = -omega[2] * s;
	term2 = (omega[4] / 3.0) * s * c[3] * (2 * psi[3] - psi[2]);
	term3 = (omega[6] / 15.0) * s * c[5] * (psi[5] * (11.0 - 24.0 * t[3]) - psi[4] * (11.0 - 36.0 * t[3]) + 2 * psi[3] * (1.0 - 7.0 * t[3]) + psi[2] * t[3]);
	term4 = (omega[8] / 315.0) * s * c[7] * (17.0 - 26.0 * t[3] + 2.0 * t[5]);
	gamma = term1 + term2 + term3 + term4;
	return gamma;
} //getGridConvergence

function getPointScaleFactorGeo() {
	// GeographicalToGrid
	var k=0.0;
	term1 = (omega[3] / 2.0) * psi[2] * c[3];
	term2 = (omega[5] / 24.0) * c[5] * ( 4 * psi[4] * (1.0 - 6.0 * t[3]) + psi[3] * (1.0 + 24.0 * t[3]) - 4.0 * psi[2] * t[3] );
	term3 = (omega[7] / 720.0) * c[7] * (61.0 - 148.0 * t[3] + 16.0 * t[5]);
	k = k0 + k0 * term1 + k0 * term2 + k0 * term3;
	return k;
} //getPointScaleFactorGeo

function getPointScaleFactorGrid() {
	// GridToGeographical
	var k=0.0;
	var x = eDash*eDash / (k0*k0 * rhoDash * nuDash);
	var x2 = x*x;
	var x3 = x2*x;
	term1 = x / 2.0;
	term2 = (x2 / 24.0) * ( 4.0 * psiDash[2] * (1.0 - 6.0 * tDash[3]) - 3.0 * (1.0 - 16.0 * tDash[3]) - 24.0 * tDash[3] / psiDash[2] );
	term3 = x3 / 720.0;
	k = k0 + k0 * term1 + k0 * term2 + k0 * term3;
	return k;
} //getPointScaleFactor

function projectionToGeographic(points){
	//var i = 1;
	for (var i=1; i<=points.length; i++) {
	if (points[i].length < 2) {
		trace(errMessages.NOT_ENOUGH_POINTS_MESSAGE+": "+i);
	}
	GridToGeographical(points[i][iEASTING], points[i][iNORTHING]);
	points[i][iLONGITUDE] = getLongitude();
	points[i][iLATITUDE] = getLatitude();
	}
	return points;
} //projectionToGeographic

function geographicToProjection(points){
	//var i = 1;
	for (var i=1; i<=points.length; i++) {
	if (points[i].length < 2) {
		trace(errMessages.NOT_ENOUGH_POINTS_MESSAGE+": "+i);
		}
		GeographicalToGrid(points[i][iLATITUDE], points[i][iLONGITUDE]);
		points[i][iEASTING] = getEasting();
		points[i][iNORTHING] = getNorthing();
	}
	return points;	
} //geographicToProjection