context: // Copyright 2005 - 2009 - Morten Nielsen (www.sharpgis.net)
  //
  // This file is part of ProjNet.
  // ProjNet 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.
  // 
  // ProjNet 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 ProjNet; 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
   *
   */
  
  using System;
  using System.Collections.Generic;
  using ProjNet.CoordinateSystems;
  using ProjNet.CoordinateSystems.Transformations;
  
  namespace ProjNet.CoordinateSystems.Projections
  {
  	/// <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>
  	internal class AlbersProjection : MapProjection
  	{
  		double _falseEasting;
  		double _falseNorthing;
  		double C;				//constant c 
  		double e;				//eccentricity
  		double e_sq = 0;
  		double ro0;
  		double n;
  		double lon_center;		//center longitude   
  		
  		#region Constructors
  
  		/// <summary>
  		/// Creates an instance of an Albers projection object.
  		/// </summary>
  		/// <param name="parameters">List of parameters to initialize the projection.</param>
  		/// <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_false_origin</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_false_origin</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>latitude_of_1st_standard_parallel</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>latitude_of_2nd_standard_parallel</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>easting_at_false_origin</term><description>The easting value assigned to the false origin.</description></item>
  		/// <item><term>northing_at_false_origin</term><description>The northing value assigned to the false origin.</description></item>
  		/// </list>
  		/// </remarks>
  		public AlbersProjection(List<ProjectionParameter> parameters)
  			: this(parameters, false)
  		{
  		}
  		
  		/// <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>
  		public AlbersProjection(List<ProjectionParameter> parameters, bool isInverse)
  			: base(parameters, isInverse)
  		{
  			this.Name = "Albers_Conic_Equal_Area";			
  
  			//Retrieve parameters
  			ProjectionParameter longitude_of_center = GetParameter("longitude_of_center");
  			ProjectionParameter latitude_of_center = GetParameter("latitude_of_center");
  			ProjectionParameter standard_parallel_1 = GetParameter("standard_parallel_1");
  			ProjectionParameter standard_parallel_2 = GetParameter("standard_parallel_2");
  			ProjectionParameter false_easting = GetParameter("false_easting");
  			ProjectionParameter false_northing = GetParameter("false_northing");
  			//Check for missing parameters
  			if (longitude_of_center == null)
  			{
  				longitude_of_center = GetParameter("central_meridian"); //Allow for altenative name
  				if (longitude_of_center == null)
  					throw new ArgumentException("Missing projection parameter 'longitude_of_center'");
  			}
  			if (latitude_of_center == null)
  			{
  				latitude_of_center = GetParameter("latitude_of_origin"); //Allow for altenative name
  				if (latitude_of_center == null)
  					throw new ArgumentException("Missing projection parameter 'latitude_of_center'");
  			}
  			if (standard_parallel_1 == null)
  				throw new ArgumentException("Missing projection parameter 'standard_parallel_1'");
  			if (standard_parallel_2 == null)
  				throw new ArgumentException("Missing projection parameter 'standard_parallel_2'");
  			if (false_easting == null)
  				throw new ArgumentException("Missing projection parameter 'false_easting'");
  			if (false_northing == null)
  				throw new ArgumentException("Missing projection parameter 'false_northing'");
  
  			lon_center = Degrees2Radians(longitude_of_center.Value);
  			double lat0 = Degrees2Radians(latitude_of_center.Value);
  			double lat1 = Degrees2Radians(standard_parallel_1.Value);
  			double lat2 = Degrees2Radians(standard_parallel_2.Value);
  			this._falseEasting = false_easting.Value * _metersPerUnit;
  			this._falseNorthing = false_northing.Value * _metersPerUnit;
  
  			if (Math.Abs(lat1 + lat2) < double.Epsilon)
  				throw new ArgumentException("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));
  			/*
  			double sin_p0 = Math.Sin(lat0);
  			double cos_p0 = Math.Cos(lat0);
  			double q0 = qsfnz(e, sin_p0, cos_p0);
  
  			double sin_p1 = Math.Sin(lat1);
  			double cos_p1 = Math.Cos(lat1);
  			double m1 = msfnz(e,sin_p1,cos_p1);
  			double q1 = qsfnz(e,sin_p1,cos_p1);
  
  
  			double sin_p2 = Math.Sin(lat2);
  			double cos_p2 = Math.Cos(lat2);
  			double m2 = msfnz(e,sin_p2,cos_p2);
  			double q2 = qsfnz(e,sin_p2,cos_p2);
  
  			if (Math.Abs(lat1 - lat2) > EPSLN)
  				ns0 = (m1 * m1 - m2 * m2)/ (q2 - q1);
  			else
  				ns0 = sin_p1;
  			C = m1 * m1 + ns0 * q1;
  			rh = this._semiMajor * Math.Sqrt(C - ns0 * q0)/ns0;
  			*/
  		}
  		#endregion
  
  		#region Public methods
  		/// <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>
          public override double[] DegreesToMeters(double[] lonlat)
  		{
  			double dLongitude = Degrees2Radians(lonlat[0]);
  			double dLatitude = Degrees2Radians(lonlat[1]);
  
  
gt: 			double a = alpha(dLatitude);