/*
 * 
 * Code is derived from the following sources
 * 
// Copyright 2005, 2006 - Morten Nielsen (www.iter.dk)
//
// This file is part of SharpMap.
// SharpMap is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
// 
// SharpMap is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public License
// along with SharpMap; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 

// SOURCECODE IS MODIFIED FROM ANOTHER WORK AND IS ORIGINALLY BASED ON GeoTools.NET:
/*
 *  Copyright (C) 2002 Urban Science Applications, Inc. 
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2.1 of the License, or (at your option) any later version.
 *
 *  This library is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */ 

// todo: delete this file?
package usda.weru.util.gis;

import java.awt.geom.Point2D;

/**
 *
 * @author wjr
 */
public class AlbersProjection {

	/**
	 *
	 */
	protected boolean _isInverse = false;

	/**
	 *
	 */
	protected boolean _isSpherical = false;

	/**
	 *
	 */
	protected double _e;

	/**
	 *
	 */
	protected double _es;

	/**
	 *
	 */
	protected double _semiMajor;

	/**
	 *
	 */
	protected double _semiMinor;
        double _falseEasting;
		double _falseNorthing;
		double e;				//eccentricity
		double e_sq = 0;
		double ro0;
		double lon_center;		//center longitude   
		

    private static final double epsilon = 1e-6;
    
    // default values
    double centralMeridian = -96.0;
    double standardParallel_1 = 29.5;
    double standardParallel_2 = 45.5;
    double lattitudeOfOrigin = 23.0;
    double semiMajor = 6378137.0;                       // equatorial radius
    double flattening = 298.257222101;
    double semiMinor = semiMajor * (1 - 1 / flattening);  // polar radius
    boolean isSpherical = (semiMajor == semiMinor);
    double falseEasting = 0;
    double falseNorthing = 0;
    
//    boolean secant = false;
    double lambda_0 = Double.NaN;
    double phi_1 = Double.NaN;
    double phi_2 = Double.NaN;
    double phi_0 = Double.NaN;
    double n = Double.NaN;
    double C = Double.NaN;
    double rho_0 = Double.NaN;
    double excentricitySquared = Double.NaN;
    double excentricity = Double.NaN;
    double ec = Double.NaN;

	/**
	 *
	 */
	public AlbersProjection() {
        init();
    }

	/**
	 *
	 * @param degrees
	 * @return
	 */
	public static double Degrees2Radians(double degrees) {
        return degrees * PI / 180;
    }

	/**
	 *
	 * @param radians
	 * @return
	 */
	public static double Radians2Degrees(double radians) {
        return radians * 180 / PI;
    }
    
    private void init() {
			this._isSpherical = (_semiMajor == _semiMinor);
			this._es = 1.0 - (_semiMinor * _semiMinor ) / ( _semiMajor * _semiMajor);
			this._e  = Math.sqrt(_es);
			lon_center = Degrees2Radians(centralMeridian);
			double lat0 = Degrees2Radians(lattitudeOfOrigin);
			double lat1 = Degrees2Radians(standardParallel_1);
			double lat2 = Degrees2Radians(standardParallel_2);
			this._falseEasting = Degrees2Radians(falseEasting);
			this._falseNorthing = Degrees2Radians(falseNorthing);
            this._semiMajor = semiMajor;
            this._semiMinor = semiMinor;

			if (Math.abs(lat1 + lat2) < epsilon)
				throw new IllegalArgumentException("Equal latitudes for standard parallels on opposite sides of Equator.");

			e_sq = 1.0 - Math.pow(this._semiMinor / this._semiMajor, 2);
			e = Math.sqrt(e_sq); //Eccentricity

			double alpha1 = alpha(lat1);
			double alpha2 = alpha(lat2);

			double m1 = Math.cos(lat1) / Math.sqrt(1 - e_sq * Math.pow(Math.sin(lat1), 2));
			double m2 = Math.cos(lat2) / Math.sqrt(1 - e_sq * Math.pow(Math.sin(lat2), 2));

			n = (Math.pow(m1, 2) - Math.pow(m2, 2)) / (alpha2 - alpha1);
			C = Math.pow(m1, 2) + (n * alpha1);

			ro0 = Ro(alpha(lat0));
    }
    


		// defines some usefull constants that are used in the projection routines
		/// <summary>
		/// PI
		/// </summary>

	/**
	 *
	 */
			protected static final double PI = Math.PI;
		/// <summary>
		/// Half of PI
		/// </summary>

	/**
	 *
	 */
			protected static final double HALF_PI = (PI*0.5);
		/// <summary>
		/// PI * 2
		/// </summary>

	/**
	 *
	 */
			protected static final double TWO_PI = (PI*2.0);
		/// <summary>
		/// EPSLN
		/// </summary>

	/**
	 *
	 */
			protected static final double EPSLN = 1.0e-10;
		/// <summary>
		/// S2R
		/// </summary>

	/**
	 *
	 */
			protected static final double S2R = 4.848136811095359e-6;
		/// <summary>
		/// MAX_VAL
		/// </summary>

	/**
	 *
	 */
			protected static final double MAX_VAL = 4;
		/// <summary>
		/// prjMAXLONG
		/// </summary>

	/**
	 *
	 */
			protected static final double prjMAXLONG = 2147483647;
		/// <summary>
		/// DBLLONG
		/// </summary>

	/**
	 *
	 */
			protected static final double DBLLONG = 4.61168601e18;


		/// <summary>
		/// Function to compute the constant small m which is the radius of
		/// a parallel of latitude, phi, divided by the semimajor axis.
		/// </summary>

	/**
	 *
	 * @param eccent
	 * @param sinphi
	 * @param cosphi
	 * @return
	 */
			protected static double msfnz (double eccent, double sinphi, double cosphi)
		{
			double con;

			con = eccent * sinphi;
			return((cosphi / (Math.sqrt(1.0 - con * con))));
		}
		
		/// <summary>
		/// Function to compute constant small q which is the radius of a 
		/// parallel of latitude, phi, divided by the semimajor axis. 
		/// </summary>

	/**
	 *
	 * @param eccent
	 * @param sinphi
	 * @param cosphi
	 * @return
	 */
			protected static double qsfnz (double eccent, double sinphi, double cosphi)
		{
			double con;

			if (eccent > 1.0e-7)
			{
				con = eccent * sinphi;
				return (( 1.0- eccent * eccent) * (sinphi /(1.0 - con * con) - (.5/eccent)*
					Math.log((1.0 - con)/(1.0 + con))));
			}
			else
				return(2.0 * sinphi);
		}

		/// <summary>
		/// Function to calculate the sine and cosine in one call.  Some computer
		/// systems have implemented this function, resulting in a faster implementation
		/// than calling each function separately.  It is provided here for those
		/// computer systems which don`t implement this function
		/// </summary>
//		protected static void sincos(double val, out double sin_val, out double cos_val) 
//
//		{ 
//			sin_val = Math.sin(val); 
//			cos_val = Math.cos(val);
//		}
		/// <summary>
		/// Function to compute the constant small t for use in the forward
		/// computations in the Lambert Conformal Conic and the Polar
		/// Stereographic projections.
		/// </summary>

	/**
	 *
	 * @param eccent
	 * @param phi
	 * @param sinphi
	 * @return
	 */
	
		protected static double tsfnz(	double eccent,
			double phi,
			double sinphi
			)
		{
			double con;
			double com;

			con = eccent * sinphi;
			com = .5 * eccent; 
			con = Math.pow(((1.0 - con) / (1.0 + con)),com);
			return (Math.tan(.5 * (HALF_PI - phi))/con);
		}
		

 


//		/// <summary>
//		/// Converts a longitude value in degrees to radians.
//		/// </summary>
//		/// <param name="x">The value in degrees to convert to radians.</param>
//		/// <param name="edge">If true, -180 and +180 are valid, otherwise they are considered out of range.</param>
//		/// <returns></returns>
//		static protected double LongitudeToRadians( double x, bool edge) 
//		{
//			if (edge ? (x>=-180 && x<=180) : (x>-180 && x<180))
//			{
//				return  Degrees2Radians(x);
//			}
//			throw new ArgumentOutOfRangeException("x",x," not a valid longitude in degrees.");
//		}
//
//  
//		/// <summary>
//		/// Converts a latitude value in degrees to radians.
//		/// </summary>
//		/// <param name="y">The value in degrees to to radians.</param>
//		/// <param name="edge">If true, -90 and +90 are valid, otherwise they are considered out of range.</param>
//		/// <returns></returns>
//		static protected double LatitudeToRadians(double y, bool edge)
//		{
//			if (edge ? (y>=-90 && y<=90) : (y>-90 && y<90))
//			{
//				return  Degrees2Radians(y);
//			}
//			throw new ArgumentOutOfRangeException("x",y," not a valid latitude in degrees.");
//		}
//		#endregion
//	}
//}

    /// <summary>
	///		Implements the Albers projection.
	/// </summary>
	/// <remarks>
	/// 	<para>Implements the Albers projection. The Albers projection is most commonly
	/// 	used to project the United States of America. It gives the northern
	/// 	border with Canada a curved appearance.</para>
	/// 	
	///		<para>The <a href="http://www.geog.mcgill.ca/courses/geo201/mapproj/naaeana.gif">Albers Equal Area</a>
	///		projection has the property that the area bounded
	///		by any pair of parallels and meridians is exactly reproduced between the 
	///		image of those parallels and meridians in the projected domain, that is,
	///		the projection preserves the correct area of the earth though distorts
	///		direction, distance and shape somewhat.</para>
	/// </remarks>


		
		/// <summary>
		/// Creates an instance of an Albers projection object.
		/// </summary>
		/// <remarks>
		/// <para>The parameters this projection expects are listed below.</para>
		/// <list type="table">
		/// <listheader><term>Items</term><description>Descriptions</description></listheader>
		/// <item><term>latitude_of_center</term><description>The latitude of the point which is not the natural origin and at which grid coordinate values false easting and false northing are defined.</description></item>
		/// <item><term>longitude_of_center</term><description>The longitude of the point which is not the natural origin and at which grid coordinate values false easting and false northing are defined.</description></item>
		/// <item><term>standard_parallel_1</term><description>For a conic projection with two standard parallels, this is the latitude of intersection of the cone with the ellipsoid that is nearest the pole.  Scale is true along this parallel.</description></item>
		/// <item><term>standard_parallel_2</term><description>For a conic projection with two standard parallels, this is the latitude of intersection of the cone with the ellipsoid that is furthest from the pole.  Scale is true along this parallel.</description></item>
		/// <item><term>false_easting</term><description>The easting value assigned to the false origin.</description></item>
		/// <item><term>false_northing</term><description>The northing value assigned to the false origin.</description></item>
		/// </list>
		/// </remarks>
		/// <param name="parameters">List of parameters to initialize the projection.</param>
		/// <param name="isInverse">Indicates whether the projection forward (meters to degrees or degrees to meters).</param>
		
        /// <summary>
		/// Converts coordinates in decimal degrees to projected meters.
		/// </summary>
		/// <param name="lonlat">The point in decimal degrees.</param>
		/// <returns>Point in projected meters</returns>
		
	/**
	 *
	 * @param lonlat
	 * @return
	 */
			
        public Point2D DegreesToMeters(Point2D lonlat)
		{
			double dLongitude = Degrees2Radians(lonlat.getX());
			double dLatitude = Degrees2Radians(lonlat.getY());

			double a = alpha(dLatitude);
			double ro = Ro(a);
			double theta = n * (dLongitude - lon_center);
			return new Point2D.Double(
				_falseEasting + ro * Math.sin(theta),
				_falseNorthing + ro0 - (ro * Math.cos(theta)));			
		}

		/// <summary>
		/// Converts coordinates in projected meters to decimal degrees.
		/// </summary>
		/// <param name="p">Point in meters</param>
		/// <returns>Transformed point in decimal degrees</returns>
 
	/**
	 *
	 * @param p
	 * @return
	 */
			public Point2D MetersToDegrees(Point2D p) 
		{
			double theta = Math.atan((p.getX() - _falseEasting) / (ro0 - (p.getY() - _falseNorthing)));
			double ro = Math.sqrt(Math.pow(p.getX() - _falseEasting, 2) + Math.pow(ro0 - (p.getY() - _falseNorthing), 2));
			double q = (C - Math.pow(ro, 2) * Math.pow(n, 2) / Math.pow(this._semiMajor, 2)) / n;
			double b = Math.sin(q / (1 - ((1 - e_sq) / (2 * e)) * Math.log((1 - e) / (1 + e))));

			double lat = Math.asin(q * 0.5);
			double preLat = Double.MAX_VALUE;
			int iterationCounter = 0;
			while (Math.abs(lat - preLat) > 0.000001)
			{
				preLat = lat;
				double sin = Math.sin(lat);
				double e2sin2 = e_sq * Math.pow(sin, 2);
				lat += (Math.pow(1 - e2sin2, 2) / (2 * Math.cos(lat))) * ((q / (1 - e_sq)) - sin / (1 - e2sin2) + 1 / (2 * e) * Math.log((1 - e * sin) / (1 + e * sin)));
				iterationCounter++;
				if (iterationCounter > 25)
					throw new IndexOutOfBoundsException("Transformation failed to converge in Albers backwards transformation");
			}
			double lon = lon_center + (theta / n);
			return new Point2D.Double(Radians2Degrees(lon), Radians2Degrees(lat));
		}
		//private double ToAuthalic(double lat)
		//{
		//    return Math.Atan(Q(lat) / Q(Math.PI * 0.5));
		//}
		//private double Q(double angle)
		//{
		//    double sin = Math.Sin(angle);
		//    double esin = e * sin;
		//    return Math.Abs(sin / (1 - Math.Pow(esin, 2)) - 0.5 * e) * Math.Log((1 - esin) / (1 + esin)));
		//}
		private double alpha(double lat)
		{
			double sin = Math.sin(lat);
			double sinsq = Math.pow(sin,2);
			return (1 - e_sq) * (((sin / (1 - e_sq * sinsq)) - 1/(2 * e) * Math.log((1 - e * sin) / (1 + e * sin))));
		}

		private double Ro(double a)
		{
			return this._semiMajor * Math.sqrt((C - n * a)) / n;
		}

	}

