/** 
 * @fileoverview This class contains all necessary functions to transform coordinates
 *	within a certain datum or between different datums
 *
 * @author Adriŕ Mercader adria.mercader@geodata.es
 * @version 1.0 
 */
 
 
//************************************************************************
// Public methods
//************************************************************************

/**
 * Transform geodesic to utm coordinates between differents datums
 * @param {double} lambda Longitude
 * @param {double} phi Latitude
 * @param {string} elipsoide_in Origin elipsoid (datum)
 * @param {string} elipsoide_out Destination elipsoid (datum)
 * @param {integer} parameters Transformation parameters code
 * @param {double} h Elipsoidal altitude (Optional, view docs)
 * @returns array with utm coordinates and the zone(x,y,zone)
*/
wcts.prototype.ged2utm_datum = function(lambda,phi,elipsoide_in,elipsoide_out,parameters,h) {
	
	//Deafult values
	if (!elipsoide_in) elipsoide_in = "WGS84";
	if (!elipsoide_out) elipsoide_out = "International 1924";
	if (!parameters) parameters = 1;
	if (!h) h = 0;

	var gec_datum1 = this._ged2gec(lambda,phi,elipsoide_in,h);
	var gec_datum2 = this._trans_gec(gec_datum1[0],gec_datum1[1],gec_datum1[2],parameters);
	var ged_datum2 = this._gec2ged(gec_datum2[0],gec_datum2[1],gec_datum2[2],elipsoide_out);
	var utm_datum2 = this.ged2utm(ged_datum2[0],ged_datum2[1],elipsoide_out);
	
	var results = new Array(utm_datum2[0],utm_datum2[1],utm_datum2[2]);
	return results;
}

/**
 * Transform utm to geodesic coordinates between differents datums
 * @param {double} x UTM X
 * @param {double} y UTM Y
 * @param {string} elipsoide_in Origin elipsoid (datum)
 * @param {string} elipsoide_out Destination elipsoid (datum)
 * @param {integer} parameters Transformation parameters code
 * @param {integer} zone UTM Zone
 * @returns array with geodesic coordinates (lambda,fi)
*/
wcts.prototype.utm2ged_datum = function (x,y,elipsoide_in,elipsoide_out,parameters,zone) {
	
	//Deafult values
	if (!elipsoide_in) elipsoide_in = "International 1924";
	if (!elipsoide_out) elipsoide_out = "WGS84";
	if (!parameters) parameters = 2;
	if (!zone) zone = 31;	
	
	var ged_datum1 = this.utm2ged(x,y,zone,elipsoide_in);
	var gec_datum1 = this._ged2gec(ged_datum1[0],ged_datum1[1],elipsoide_in);
	var gec_datum2 = this._trans_gec(gec_datum1[0],gec_datum1[1],gec_datum1[2],parameters);
	var ged_datum2 = this._gec2ged(gec_datum2[0],gec_datum2[1],gec_datum2[2],elipsoide_out);
	
	results = new Array(ged_datum2[0],ged_datum2[1]);
	return results;
	
}

/**
 * Transform geodesic to utm coordinates within the same datum
 * @param {double} lambda Longitude
 * @param {double} phi Latitude
 * @param {string} elipsoide Elipsoid (datum)
 * @returns array with utm coordinates and the zone(x,y,zone)
*/
wcts.prototype.ged2utm = function (lambda,phi,elipsoide) {

	//Deafult values
	if (!elipsoide) elipsoide = "WGS84";

	//Retrieve datum parameters
	var datum = this._datums[elipsoide];
	
	//Calculate zone
	var zone = Math.floor((lambda/6)+31);
	
	//Calculem el meridiŕ central del fus
	var meridia = (6 * zone) - 183;
	
	//Passem a radians les coordenades
	lambda = lambda * (this.pi /180);
	phi = phi * (this.pi /180);
	
	//Calculate angular distance between point longitude and xone central meridian
	var dif_lambda = lambda - (_deg2rad(meridia));
	
	//Calculate conversion parameters
	var A = Math.cos(phi)*Math.sin(dif_lambda);
	var xi = (1/2)*Math.log((1+A)/(1-A));
	var eta = Math.atan((Math.tan(phi))/(Math.cos(dif_lambda)))-phi;
	var ni = (datum["c"] /Math.pow(1+datum["e_2"]*Math.pow(Math.cos(phi),2),(1/2))  )*0.9996;
	var Zeta = (datum["e_2"]/2)*Math.pow(xi,2)*Math.pow(Math.cos(phi),2);
	var a1 = Math.sin(2*phi);
	var a2 = +a1 * (Math.pow(Math.cos(phi),2)); 
	var j2 = phi + (a1/2);
	var j4 = ((3*j2)+a2)/4;
	var j6 = (5*j4+a2*(Math.pow(Math.cos(phi),2)))/3;
	var alfa = (3/4)*datum["e_2"];
	var beta = (5/3)*Math.pow(alfa,2);
	var gamma = (35/27)*Math.pow(alfa,3);
	var B_fi = 0.9996*datum["c"]*(phi-(alfa*j2)+(beta*j4)-(gamma*j6));
	
	//Calculate coordinates and zone
	this.x = xi*ni*(1+Zeta/3)+500000;
	this.y = eta*ni*(1+Zeta)+B_fi;
	if (phi < 0) this.y = this.y + 10000000;
	this.zone = zone;
	
	var results = new Array(this.x,this.y,this.zone);
	return results;
}

/**
 * Transform utm to geodesic coordinates within the same datum
 * @param {double} x UTM X
 * @param {double} y UTM Y
 * @param {integer} zone UTM Zone
 * @param {string} elipsoide Elipsoid (datum)
 * @returns array with geodesic coordinates and the zone(lambda,phi)
*/
wcts.prototype.utm2ged = function(x,y,zone,elipsoide) {
	
	//Deafult values
	if (!zone) zone = 31;
	if (!elipsoide) elipsoide = "WGS84";

	//Retrieve datum parameters
	var datum = this._datums[elipsoide];
	
	//Calculate zone central meridian
	var meridia = (6 * zone) - 183;
	
	//Calculate necessary parameters
	var phi_ = y / (6366197.724*0.9996);
	var ni = (datum["c"] /Math.pow(1+datum["e_2"]*Math.pow(Math.cos(phi_),2),(1/2))  )*0.9996;
	var a =(x-500000)/ni;
	var a1 = Math.sin(2*phi_);
	var a2 = a1 * (Math.pow(Math.cos(phi_),2)); 	
	var j2 = phi_ + (a1/2);
	var j4 = ((3*j2)+a2)/4;
	var j6 = (5*j4+a2*(Math.pow(Math.cos(phi_),2)))/3;
	var alfa = (3/4)*datum["e_2"];
	var beta = (5/3)*Math.pow(alfa,2);
	var gamma = (35/27)*Math.pow(alfa,3);
	var B_fi = 0.9996*datum["c"]*(phi_-(alfa*j2)+(beta*j4)-(gamma*j6));
	var b = (y-B_fi)/ni;
	var Zeta = ((datum["e_2"]*Math.pow(a,2))/2)*Math.pow((Math.cos(phi_)),2);
	var xi = a*(1-(Zeta/3));
	var eta = b*(1-Zeta)+phi_;
	var sin_h_xi = (Math.exp(xi)-Math.exp(xi*-1))/2;
	var dif_lambda = Math.atan(sin_h_xi/Math.cos(eta));
	var tau = Math.atan(Math.cos(dif_lambda)*Math.tan(eta));
	
	var phi_rad = phi_ + (1 + datum["e_2"]*(Math.pow(Math.cos(phi_),2)) - (3/2)*datum["e_2"]*Math.sin(phi_)*Math.cos(phi_)*(tau-phi_))*(tau-phi_);

	//Calculate coordinates
	this.lambda = (dif_lambda/this.pi)*180+meridia;
	this.phi = (phi_rad/this.pi)*180;
	
	var results = new Array(this.lambda,this.phi);
	return results;

}

/**
 * Transform geodesic coordinates in decimal format to Degrees/Minutes/Seconds format (ş ' '')
 * @param {double} value geodesic coordinate in decimal format
 * @returns array with coordiante en DMS format (hemisphere,degree,minutes,seconds)
*/
wcts.prototype.deg2dms = function(value) {
  
  var dms = new Array(3);
  
  dms[0] = ( value < 0 ) ? -1 : 1;
  value = Math.abs( value );
  dms[1] = parseInt(value);
  value = 60 * (value - dms[1]);
  dms[2] = parseInt(value);
  value = 60 * (value - dms[2]);
  dms[3] = Math.round(100*value)/100;

  return dms;
    
}

	
//************************************************************************
// Private methods
//************************************************************************

/**
 * Transform geodesic to geocentric coordinates within the same datum
 * @param {double} lambda Longitude
 * @param {double} phi Latitude
 * @param {double} h Elipsoidal altitude (Optional, view docs)
 * @param {string} elipsoide Elipsoid (datum)
 * @returns array with geocentric coordinates (x,y,z)
*/
function _ged2gec(lambda,phi,elipsoide,h){
	//Deafult values
	if (!elipsoide) elipsoide = "WGS84";
	if (!h) h = 0;
	
	//Retrieve datum parameters
	var datum = this._datums[elipsoide];
	
	//Degree to radian transformation
	lambda = _deg2rad(lambda);
	phi = _deg2rad(phi);
	
	//Calculate parameter
	var n =(Math.pow(datum["a"],2))/Math.sqrt(Math.pow(datum["a"],2)* Math.pow(Math.cos(phi),2)  +  Math.pow(datum["b"],2)*Math.pow(Math.sin(phi),2)  );
	
	//Calculate coordinates
	this.x_ecef = +(n+h)*Math.cos(phi)*Math.cos(lambda);
	this.y_ecef = +(n+h)*Math.cos(phi)*Math.sin(lambda);
	this.z_ecef =((Math.pow(datum["b"],2)/Math.pow(datum["a"],2))*n+h)*Math.sin(phi);

	var results = new Array(this.x_ecef,this.y_ecef,this.z_ecef);
	return results;
}

/**
 * Transforms geocentric coordinates between datums
 * @param {double} x_ecef X
 * @param {double} y_ecef Y
 * @param {double} z_ecef Z
 * @param {integer} parameters Transformation parameter code (view trans_params)
 * @returns array with geocentric coordinates (x,y,z)
*/

function _trans_gec(x_ecef,y_ecef,z_ecef,parameters) {
	//Retrieve transformation parameters
	var params = this._trans_params[parameters];
	
	//Transform scale factor
	var S = params["S"]/ 1000000;
	
	//Calculate coordinates transformed to the destination datum
	this.x_trans = +params["AX"]+((1+S)*(x_ecef+(y_ecef*_deg2rad(params["RZ"]/60/60))-(z_ecef*_deg2rad(params["RY"]/60/60))));
	this.y_trans = +params["AY"]+((1+S)*(-(_deg2rad(params["RZ"]/60/60)*x_ecef)+y_ecef+(z_ecef*_deg2rad(params["RX"]/60/60))));
	this.z_trans = +params["AZ"]+((1+S)*((x_ecef*_deg2rad(params["RY"]/60/60))-(y_ecef*_deg2rad(params["RX"]/60/60))+z_ecef));

	var results = array(this.x_trans,this.y_trans,this.z_trans);
	return results;
}

/**
 * Transforms geocentric to geodesic coordinates within the same datum
 * @param {double} x_ecef X
 * @param {double} y_ecef Y
 * @param {double} z_ecef Z
 * @param {string} elipsoide Elipsoid (datum)
 * @returns array with geodesic coordinates (lambda,phi,h)
*/
function _gec2ged(x_ecef,y_ecef,z_ecef,elipsoide) {
	
	//Deafult values
	if (!elipsoide) elipsoide = "International 1924";

	//Retrieve datum parameters
	var datum = this._datums[elipsoide];
	
	//Calculate necessary parameters
	var e2 = (Math.pow(datum["a"],2)-Math.pow(datum["b"],2))/(Math.pow(datum["a"],2));
	var e_2 = (Math.pow(datum["a"],2)-Math.pow(datum["b"],2))/(Math.pow(datum["b"],2));
	var p = +Math.sqrt(Math.pow(x_ecef,2)+Math.pow(y_ecef,2));
	var theta = Math.atan((z_ecef*datum["a"])/(p*datum["b"]));
	var phi_rad = +Math.atan((+z_ecef+e_2*datum["b"]*Math.pow(Math.sin(theta),3))/(p-e2*datum["a"]*Math.pow(Math.cos(theta),3)));
	var lambda_rad =+Math.atan(y_ecef/x_ecef);
	var n = +Math.pow(datum["a"],2)/(Math.sqrt(Math.pow(datum["a"],2)*Math.pow(Math.cos(phi_rad),2)+ Math.pow(datum["b"],2)*Math.pow(Math.sin(phi_rad),2)));
	var h = +(p/Math.cos(phi_rad))-n;
	
	//Calculate coordinates (decimal geodesic)
	this.lambda =_rad2deg(lambda_rad);
	this.phi =_rad2deg(phi_rad);
	this.h = h;

	var results = new Array(this.lambda,this.phi,this.h);
	return results;

}

/**
 * Converts degrees to radians
 * @param {double} value in degrees
 * @returns double in radians
*/
function _deg2rad(value) {
	return value * (Math.PI /180);
}

/**
 * Converts radians to degrees
 * @param {double} value in radians
 * @returns double in degrees
*/
function _rad2deg(value) {
	return value * (180 / Math.PI);
}


//************************************************************************
// Constructor
//************************************************************************
/**
 * Constructs a new wcts object.
 * @class
 * This is the main object   
 * @constructor
 * @return A new wcts object
 */
function wcts(){
  //constants
  this.pi = Math.PI;
  
  /**
  * Datum parameters
 	* a: major semiaxis, b: minor semiaxis, e: excentricity, e_: 2nd excentricity, e_2: square of 2nd excentricity, c: polar radii of corbature
	* Note: ED50 =+- International 1909/1924
  * @type array
  * @private
  */
  this._datums = new Array();
		this._datums["Airy 1830"] = {"a" : 6377563.396,"b" : 6356256.91, "e" : 0.081673372, "e_" : 0.081947146, "e_2" : 0.006715335, "c" : 6398941.302}; 
		this._datums["Airy Modified 1965"] = {"a" : 6377340.189,"b" : 6356034.4479, "e" : 0.081673374, "e_" : 0.081947147, "e_2" : 0.006715335, "c" : 6398717.348};
		this._datums["Bessel 1841"] = {"a" : 6377397.155000,"b" : 6356078.962840, "e" : 0.081696831, "e_" : 0.081970841, "e_2" : 0.006719219, "c" : 6398786.848};
		this._datums["Clarke 1866"] = {"a" : 6378206.400000,"b" : 6356583.800000, "e" : 0.082271854, "e_" : 0.082551711, "e_2" : 0.006814785, "c" : 6399902.552};
		this._datums["Clarke 1880"] = {"a" : 6378249.145000,"b" : 6356514.869550, "e" : 0.0824834, "e_" : 0.082765428, "e_2" : 0.006850116, "c" : 6400057.735};
		this._datums["Fischer 1960"] = {"a" : 6378166.000000,"b" : 6356784.280000, "e" : 0.081813341, "e_" : 0.082088529, "e_2" : 0.006738527, "c" : 6399619.64};
		this._datums["Fischer 1968"] = {"a" : 6378150.000000,"b" : 6356768.330000, "e" : 0.081813348, "e_" : 0.082088536, "e_2" : 0.006738528, "c" : 6399603.59};
		this._datums["GRS80"] = {"a" : 6378137.000000,"b" : 6356752.314140, "e" : 0.081819191, "e_" : 0.082094438, "e_2" : 0.006739497, "c" : 6399593.626};
		this._datums["Hayford 1909"] = {"a" : 6378388.000000,"b" : 6356911.946130, "e" : 0.08199189, "e_" : 0.08226889, "e_2" : 0.00676817, "c" : 6399936.608};
		this._datums["Helmert 1906"] = {"a" : 6378200.000000,"b" : 6356818.170000, "e" : 0.081813333, "e_" : 0.082088521, "e_2" : 0.006738525, "c" : 6399653.75};
		this._datums["Hough 1960"] = {"a" : 6378270.000000,"b" : 6356794.343479, "e" : 0.08199189, "e_" : 0.08226889, "e_2" : 0.00676817, "c" : 6399818.209};
		this._datums["International 1909"] = {"a" :	6378388.000000,"b" : 6356911.946130, "e" : 0.08199189, "e_" : 0.08226889, "e_2" : 0.00676817, "c" : 6399936.608};
		this._datums["International 1924"] = {"a" :	6378388.000000,"b" : 6356911.946130, "e" : 0.08199189, "e_" : 0.08226889, "e_2" : 0.00676817, "c" : 6399936.608};
		this._datums["Krasovsky 1940"] = {"a" : 6378245.000000,"b" : 6356863.018800, "e" : 0.081813334, "e_" : 0.082088522, "e_2" : 0.006738525, "c" : 6399698.902};
		this._datums["Mercury 1960"] = {"a" : 6378166.000000,"b" : 6356784.283666, "e" : 0.081813334, "e_" : 0.082088522, "e_2" : 0.006738525, "c" : 6399619.636};
		this._datums["Mercury Modificado 1968"] = {"a" : 6378150.000000,"b" : 6356768.337303, "e" : 0.081813334, "e_" : 0.082088522, "e_2" : 0.006738525, "c" : 6399603.582};
		this._datums["Nuevo International 1967"] = {"a" : 6378157.500000,"b" : 6356772.200000, "e" : 0.081820233, "e_" : 0.08209549, "e_2" : 0.00673967, "c" : 6399614.744};
		this._datums["Sudamericano 1969"] = {"a" : 6378160.000000,"b" : 6356774.720000, "e" : 0.081820178, "e_" : 0.082095436, "e_2" : 0.006739661, "c" : 6399617.224};
		this._datums["Walbeck 1817"] = {"a" : 6376896.000000,"b" : 6355834.846700, "e" : 0.081206823, "e_" : 0.081475916, "e_2" : 0.006638325, "c" : 6398026.943};
		this._datums["WGS66"] = {"a" : 6378145.000000,"b" : 6356759.769356, "e" : 0.08182018, "e_" : 0.082095437, "e_2" : 0.006739661, "c" : 6399602.174};
		this._datums["WGS72"] = {"a" : 6378135.000000,"b" : 6356750.519915, "e" : 0.081818811, "e_" : 0.082094054, "e_2" : 0.006739434, "c" : 6399591.419};
		this._datums["WGS84"] = {"a" : 6378137.000000,"b" : 6356752.314, "e" : 0.081819191, "e_" : 0.082094438, "e_2" : 0.006739497, "c" : 6399593.626};
		
  /**
  * Transformation parameters
  * -Parameters codes
  *					1: Espanya (Peninsula Ibčrica Principal) de ETRS89 (~WGS84) a ED50	
  *					2: Espanya (Peninsula Ibčrica Principal) de ED50 a ETRS89 (~WGS84)	
  *					3: Espanya (Nordest Península Ibčrica) de ETRS89 (~WGS84) a ED50	
  *					4: Espanya (Nordest Península Ibčrica) de ED50 a ETRS89 (~WGS84)	
  *					5: Espanya (Illes Balears) de ETRS89 (~WGS84) a ED50	
  *					6: Espanya (Illes Balears) de ED50 a ETRS89 (~WGS84)	
  *	- AX,AY,AZ: distance increments RX,RY,RZ: rotations S: scale factor
  * @type array
  * @private
  */
	this._trans_params = new Array();
		this._trans_params[1] = {"AX" : 131.0320,"AY" : 100.2510,"AZ" : 163.3540,"RX" : -1.2438,"RY" : -0.0195,"RZ" : -1.1436,"S" : -9.3900};
		this._trans_params[2] = {"AX" : -131.0320,"AY" : -100.2510,"AZ" : -163.3540,"RX" : 1.2438,"RY" : 0.0195,"RZ" : 1.1436,"S" : 9.3900};
		this._trans_params[3] = {"AX" : 178.3830,"AY" : 83.1720,"AZ" : 221.2930,"RX" : 0.5401,"RY" : -0.5319,"RZ" : -0.1263,"S" : -21.2000};
		this._trans_params[4] = {"AX" : -178.3830,"AY" : -83.1720,"AZ" : -221.2930,"RX" : -0.5401,"RY" : 0.5319,"RZ" : 0.1263,"S" : 21.2000};
		this._trans_params[5] = {"AX" : 181.4609,"AY" : 90.2931,"AZ" : 187.1902,"RX" : 0.1435,"RY" : 0.4922,"RZ" : -0.3935,"S" : -17.5700};
		this._trans_params[6] = {"AX" : -181.4609,"AY" : -90.2931,"AZ" : -187.1902,"RX" : -0.1435,"RY" : -0.4922,"RZ" : 0.3935,"S" : 17.5700};

}