/*******************************************************************************
NAME                     ALBERS CONICAL EQUAL AREA 

PURPOSE:	Transforms input longitude and latitude to Easting and Northing
		for the Albers Conical Equal Area projection.  The longitude
		and latitude must be in radians.  The Easting and Northing
		values will be returned in meters.

PROGRAMMER              DATE
----------              ----
T. Mittan,       	Feb, 1992

ALGORITHM REFERENCES

1.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
    Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
    State Government Printing Office, Washington D.C., 1987.

2.  Snyder, John P. and Voxland, Philip M., "An Album of Map Projections",
    U.S. Geological Survey Professional Paper 1453 , United State Government
    Printing Office, Washington D.C., 1989.
*******************************************************************************/
#include "gctpc/cproj.h"
#include "gctpc/report.h"

/* Variables common to all subroutines in this code file
  -----------------------------------------------------*/
  static double r_major;	/* major axis				*/
  static double r_minor;	/* minor axis				*/
  static double c;		/* constant c				*/
  static double e3;		/* eccentricity 			*/
  static double rh;		/* heigth above elipsoid		*/
  static double ns0;		/* ratio between meridians		*/
  static double lon_center;	/* center longitude			*/
  static double false_easting;	/* x offset in meters			*/
  static double false_northing;	/* y offset in meters			*/

/* Initialize the Albers projection
  -------------------------------*/
long
alberforint(double r_maj, double r_min, double lat1, double lat2, double lon0, double lat0, double false_east, double false_north)

             			/* major axis				*/
             			/* minor axis				*/
            			/* first standard parallel		*/
            			/* second standard parallel		*/
            			/* center longitude			*/
            			/* center lattitude			*/
                  		/* x offset in meters			*/
                   		/* y offset in meters			*/
{
double msfnz(double eccent, double sinphi, double cosphi);			/* function to compute small m		*/
double qsfnz(double eccent, double sinphi, double cosphi);			/* function to compute small q		*/
double sin_po,cos_po;		/* sin and cos values			*/
double con;			/* temporary variable			*/
double es,temp;			/* eccentricity squared and temp var	*/
double ms1;			/* small m 1				*/
double ms2;			/* small m 2				*/
double qs0;			/* small q 0				*/
double qs1;			/* small q 1				*/
double qs2;			/* small q 2				*/

false_easting = false_east;
false_northing = false_north;
lon_center = lon0;
if (fabs(lat1 + lat2) < EPSLN)
   {
   p_error("Equal latitudes for St. Parallels on opposite sides of equator",
	  "alber-forinit");
   return(31);
   }
r_major = r_maj;
r_minor = r_min;
temp = r_minor / r_major;
es = 1.0 - SQUARE(temp);
e3 = sqrt(es);

sincos(lat1, &sin_po, &cos_po);
con = sin_po;

ms1 = msfnz(e3,sin_po,cos_po);
qs1 = qsfnz(e3,sin_po,cos_po);

sincos(lat2,&sin_po,&cos_po);

ms2 = msfnz(e3,sin_po,cos_po);
qs2 = qsfnz(e3,sin_po,cos_po);

sincos(lat0,&sin_po,&cos_po);

qs0 = qsfnz(e3,sin_po,cos_po);

if (fabs(lat1 - lat2) > EPSLN)
   ns0 = (ms1 * ms1 - ms2 *ms2)/ (qs2 - qs1);
else
   ns0 = con;
c = ms1 * ms1 + ns0 * qs1;
rh = r_major * sqrt(c - ns0 * qs0)/ns0;

/* Report parameters to the user
  -----------------------------*/
ptitle("ALBERS CONICAL EQUAL-AREA");
radius2(r_major, r_minor);
stanparl(lat1,lat2);
cenlonmer(lon_center);
origin(lat0);
offsetp(false_easting,false_northing);

return(OK);
}

/* Albers Conical Equal Area forward equations--mapping lat,long to x,y
  -------------------------------------------------------------------*/
long
alberfor(double lon, double lat, double *x, double *y)
           			/* (I) Longitude 		*/
           			/* (I) Latitude 		*/
          			/* (O) X projection coordinate 	*/
          			/* (O) Y projection coordinate 	*/
{
double sin_phi,cos_phi;		/* sine and cos values		*/
double qs;			/* small q			*/
double theta;			/* angle			*/ 
double rh1;			/* height above ellipsoid	*/
double qsfnz(double eccent, double sinphi, double cosphi);			/* function to calculate q	*/
double adjust_lon(double x);    /* Function to adjust longitude to -180 - 180   */

sincos(lat,&sin_phi,&cos_phi);
qs = qsfnz(e3,sin_phi,cos_phi);
rh1 = r_major * sqrt(c - ns0 * qs)/ns0;
theta = ns0 * adjust_lon(lon - lon_center); 
*x = rh1 * sin(theta) + false_easting;
*y = rh - rh1 * cos(theta) + false_northing;

return(OK);
}
