// ********************************************************************************************
// Projection algorithms were used from:
// Map Projections: A Working Manual
// John P. Snyder
// USGS Professional Paper 1395
// 1987
// pages 97-103
// 
// *** All projections are ellipsoid-based ***
//
//Implemented by Dylan Beaudette Summer 2005
//University of California at Davis
//Davis, Ca 95616
// ********************************************************************************************

/**
Projection algorithms were used from: Map Projections A Working Manual. John P. Snyder. USGS Professional Paper 1395, 1987. pages 97-103
\file aea.js
\author Dylan Beaudette
\date Summer 2005
\brief Forward and reverse elipsoid based projections to Albers Equal Area Conic

*/


// note all JS trig functions require arguments to be in RADIANS!
/// constant for conversion to radians from degrees
var degrees2radians = (Math.PI / 180.0);

/// constant for conversion back to degrees
var radians2degrees = (180.0 / Math.PI);

// ********************************************************************************************
// Ellipsoid parameters for NAD83 -- GRS80 ellipsoid
// from Elements of Cartography 5th Ed. Robinson et. al. 1995.
// ********************************************************************************************

///major axis
var a = 6378137 ;
///minor axis
var b = 6356752.3 ;
///eccentricity squared
var e2 = (a * a - b * b) / (a * a);
///eccentricity 
var e = Math.sqrt(e2) ;

// ********************************************************************************************
// Map constants: used in both forward and inverse projections
// ********************************************************************************************

// latitude and longitude of origin : these are the parameters for our project
// these values corrospond to phi_0 and lambda_0 in USGS paper 1395
	
	///proj4-style lat of origin (40 deg N)
	var lat_0 = 40.0 ;
	///proj4-style lon of origin (96 deg W)
	var lon_0 = -96.0 ;
	///convert to radians
	var phi_0 = lat_0 * degrees2radians ;
	///convert to radians
	var lamda_0 = lon_0 * degrees2radians;
	
	
// std parallels : phi_1 and phi_2 in USGS paper 1395
	
	///proj4-style std parallel 1 (20 deg N)
	var lat_1 = 20.0 ;
	///proj4-style std parallel 2 (60 deg N)
	var lat_2 = 60.0 ;
	///convert to radians
	var phi_1 = lat_1 * degrees2radians ;
	///convert to radians
	var phi_2 = lat_2 * degrees2radians ;


// m-terms: m_1 and m_2 --  Equation 14-15
	var frac_sin_sq_1 = 1 - (e2 * Math.pow(Math.sin(phi_1), 2)) ;
	var frac_sin_sq_2 = 1 - (e2 * Math.pow(Math.sin(phi_2), 2)) ;
	
	var m_1 = Math.cos(phi_1 )/ Math.sqrt(frac_sin_sq_1 )  ;
	var m_2 = Math.cos(phi_2 )/ Math.sqrt(frac_sin_sq_2 )  ;


// q-terms: q_0 q_1 q_2 -- Equation 3-12
	var log_term_0 = Math.log( (1 - e * Math.sin(phi_0) ) / (1 + e * Math.sin(phi_0) ) );
	var log_term_1 = Math.log( (1 - e * Math.sin(phi_1) ) / (1 + e * Math.sin(phi_1) ) );
	var log_term_2 = Math.log( (1 - e * Math.sin(phi_2) ) / (1 + e * Math.sin(phi_2) ) );
	
	var q_0 =  (1 - e2) * ( Math.sin(phi_0)/(1 - e2 * Math.pow(Math.sin(phi_0), 2)) - (1/(2 * e)) *  log_term_0) ;
	var q_1 =  (1 - e2) * ( Math.sin(phi_1)/(1 - e2 * Math.pow(Math.sin(phi_1), 2)) - (1/(2 * e)) *  log_term_1) ;
	var q_2 =  (1 - e2) * ( Math.sin(phi_2)/(1 - e2 * Math.pow(Math.sin(phi_2), 2)) - (1/(2 * e)) *  log_term_2) ;

// n, C, and rho terms:
	/// Equation 14-14 
	var n =  (Math.pow(m_1, 2) - Math.pow(m_2, 2) ) / (q_2 - q_1);
	/// Equation 14-13
	var C =  Math.pow(m_1, 2) + n * q_1;
	/// Equation 14-12a
	var rho_0 = a * Math.sqrt(C - n * q_0) / n ;

	

/**
\brief Forward AEA projection: latlon to AEA.
\param[in] latlon
 	an array: latlon[0] = x, latlon[1] = y, latlon[2]
\return aea
	an array: aea[0] = x, aea[1] = y units are in meters
\see shape.js for definition of point object 
 */
function geographic2aea(latlon)
	{
	//init a new array
	var aea = new Array();
	
	//extract our x and y latlon values and convert to radians
	var lamda = latlon[0] * degrees2radians ; 
	var phi = latlon[1]  * degrees2radians;
	
	
	// Calculate q
	// Equation 3-12 with phi = input latitude
	var log_term = Math.log( (1 - e * Math.sin(phi) ) / (1 + e * Math.sin(phi) ) );
	var q = (1 - e2) * ( Math.sin(phi)/(1 - e2 * Math.pow(Math.sin(phi), 2)) - (1/(2 * e)) *  log_term) ;
	
	// Equation 14-12 : Calculate rho
	var rho = a * Math.sqrt(C - (n*q)) / n ; 
	
	
	// Equation 14-4 : Calculate theta
	var theta = n*(lamda - lamda_0) ;
	
		
	// Equation 14-1
	var x = rho * Math.sin(theta) ;
	
	
	// Equation 14-2
	var y = rho_0 - rho * Math.cos(theta) ;
	
	//setup the point to return
	aea[0] = x ;
	aea[1] = y ;
	
	//return the latitude and longitude as as point class.
	return aea;	
	} //end function	

	

/**
\brief Inverse AEA projection: AEA to latlon.
\param[in] aea
 	an array: aea[0] = x, aea[1] = y
\return latlon
	an array: latlon[0] = x, latlon[1] = y, latlon[2] = error_factor
\see shape.js for definition of point object
\note According to the USGS paper, recursion should be used to calculate the exact latlon. None is used here as an approximation is good enough. 
 */
function aea2geographic(aea)
	{
	
	//init a tolerance for an acceptible error in phi
	var phi_tol = 0.001;
	//init a new array
	var latlon = new Array();
	
	//extract our x and y AEA map coordinates
	var x = aea[0] ; 
	var y = aea[1];
	
	//init some new variables specific to an inverse projection:
	// Equation 14-11
	var theta = Math.atan2(x,(rho_0 - y)) ;
	// Equation 14-10
	var rho = Math.sqrt( (Math.pow(x, 2) + Math.pow(rho_0 - y, 2) ) ) ;
	// Equation 14-19
	var q = (C - Math.pow(rho, 2) * Math.pow(n, 2) / Math.pow(a, 2) ) / n ;
	
	//calculate the latitude phi: eq. 3-16 USGS paper 1395
	var phi_guess = Math.asin(q/2); //estimate an initial phi
	var first_term = Math.pow(1 - e2 * Math.pow(Math.sin(phi_guess), 2), 2)/(2 * Math.cos(phi_guess)) ;
	var second_term = (q / (1 - e2) ) - ( (Math.sin(phi_guess)) / (1 - e2 * Math.pow(Math.sin(phi_guess),2)) ) ;
	var third_term = ( 1 / (2*e) ) * Math.log( (1 - e * Math.sin(phi_guess))/(1 + e * Math.sin(phi_guess)) ) ;
	var phi_initial = phi_guess + first_term * (second_term + third_term);
	
	// a bit of recursion to make inverse projection more accurate
	// right now it is non-functioning
	// points are off by about 36m to the SSW
	var phi_error1 = phi_initial - phi_guess;
	
// 	if(phi_error1 > phi_tol)
// 		{
// 		var phi_next = phi_initial + first_term * (second_term + third_term);
// 		}
	
	
	phi_final = phi_initial * radians2degrees; //convert to degrees
	
	//calculate the longitude: eq. 14-9 USGS paper 1395
	var lamda = lamda_0 + (theta/n);
	lamda = lamda * radians2degrees; //convert to degrees
	
	// setup the point to return
	// longitude:
	latlon[0] = lamda ;
	// latitude
	latlon[1] = phi_final ;
	// error in latitude
	latlon[2] = phi_error1 * radians2degrees;
	
	//return the latitude and longitude!
	return latlon;
	}