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

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

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 es;		/* eccentricity squared			*/
  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
alberinvint(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;		/* sine and cos values			*/
double con;			/* temporary variable			*/
double temp;			/* temporary variable			*/
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 Standard Parallels on opposite sides of equator"
	  ,"alber-invinit");
   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 inverse equations--mapping x,y to lat/long
  -------------------------------------------------------------------*/
long
alberinv(double x, double y, double *lon, double *lat)
         			/* (O) X projection coordinate 	*/
         			/* (O) Y projection coordinate 	*/
            			/* (I) Longitude 		*/
            			/* (I) Latitude 		*/
{
double qsfnz(double eccent, double sinphi, double cosphi);			/* function to calculate q	*/
double rh1;			/* height above ellipsoid	*/
double qs;			/* function q			*/
double con;			/* temporary sign value		*/
double theta;			/* angle			*/
double adjust_lon(double x);    /* Function to adjust longitude to -180 - 180   */
double phi1z (double eccent, double qs, long int *flag);		/* function to calculate LAT phi	*/
long   flag;			/* error flag;			*/


flag = 0;
x -= false_easting;
y = rh - y + false_northing;;
if (ns0 >= 0)
   {
   rh1 = sqrt(x * x + y * y);
   con = 1.0;
   }
else
   {
   rh1 = -sqrt(x * x + y * y);
   con = -1.0;
   }
theta = 0.0;
if (rh1 != 0.0)
   theta = atan2(con * x, con * y);
con = rh1 * ns0 / r_major;
qs = (c - con * con) / ns0;
if (e3 >= 1e-10)
   {
   con = 1 - .5 * (1.0 - es) * log((1.0 - e3) / 1.0 + e3)/e3;
   if (fabs(fabs(con) - fabs(qs)) > .0000000001 )
      {
      *lat = phi1z(e3,qs,&flag);
      if (flag != 0)
         return(flag);
      }
   else
      {
      if (qs >= 0)
         *lat = .5 * PI;
      else
         *lat = -.5 * PI;
      }
   }
else
   {
   *lat = phi1z(e3,qs,&flag);
   if (flag != 0)
      return(flag);
   }

*lon = adjust_lon(theta/ns0 + lon_center);

return(OK);
}