/* File: convert_albers.c

   Simple program to convert an x/y location, in meters, to a lat/lon
   location for NLCD 1992/2001 data that give the corner-point coordinates
   in x/y. The transformation assumes the x/y location is with respect
   to an Albers Equal Area projection using the NAD83 datum.

   When run, the program takes, as command-line arguments, the x and y
   locations; the corresponding latitude and longitude are printed to
   standard output.

   References:
   Snyder, John P., "Map Projections -- A Working Manual", 
   USGS Professional Paper 1395, United State Government Printing Office, 
   Washington D.C., 1987.


   Michael G. Duda, NCAR/MMM
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define RE_NAD83 6378137.0  
#define E_NAD83  0.0818187034  /* Eccentricity */
#define D2R      0.017453293   /* PI/180 */

int main(int argc, char ** argv)
{
   double x, y, lat, lon;
   double truelat1, truelat2, stand_lon, lat1, lon1;
   double rho, rho0, theta, m1, m2, sinphi, C, nc, q1, q2, q0, q, beta, a, b, c;

   /* True latitudes and standard longitude used for NLCD */
   truelat1  =  29.5;
   truelat2  =  45.5;
   stand_lon = -96.0;
 
   /* Lat and lon of origin */
   lat1      =  23.0;
   lon1      = -96.0;

   if (argc != 3)
   {
      fprintf(stderr,"Usage: convert_albers x y\n\n");
      return 1;
   }

   x = atof(argv[1]);
   y = atof(argv[2]);

   m1 = cos(truelat1*D2R)/sqrt(1.0-pow((E_NAD83*sin(truelat1*D2R)),2.0));
   m2 = cos(truelat2*D2R)/sqrt(1.0-pow((E_NAD83*sin(truelat2*D2R)),2.0));

   sinphi = sin(truelat1*D2R);
   q1 = (1.0-pow(E_NAD83,2.0))*((sinphi/(1.0-pow((E_NAD83*sinphi),2.0))) - 1.0/(2.0*E_NAD83) * log((1.0-E_NAD83*sinphi)/(1.0+E_NAD83*sinphi)));

   sinphi = sin(truelat2*D2R);
   q2 = (1.0-pow(E_NAD83,2.0))*((sinphi/(1.0-pow((E_NAD83*sinphi),2.0))) - 1.0/(2.0*E_NAD83) * log((1.0-E_NAD83*sinphi)/(1.0+E_NAD83*sinphi)));

   nc = (pow(m1,2.0) - pow(m2,2.0)) / (q2 - q1);

   C = pow(m1,2.0) + nc*q1;

   sinphi = sin(lat1*D2R);
   q0 = (1.0-pow(E_NAD83,2.0))*((sinphi/(1.0-pow((E_NAD83*sinphi),2.0))) - 1.0/(2.0*E_NAD83) * log((1.0-E_NAD83*sinphi)/(1.0+E_NAD83*sinphi)));

   rho0 = RE_NAD83 * sqrt(C - nc * q0) / nc;

   rho = sqrt(pow(x,2.0) + pow((rho0 - y),2.0));

   q = (C - pow((rho*nc/RE_NAD83),2.0)) / nc;

   beta = asin(q/(1.0 - log((1.0-E_NAD83)/(1.0+E_NAD83))*(1.0-pow(E_NAD83,2.0))/(2.0*E_NAD83)));

   a = 1./3.*pow(E_NAD83,2.) + 31./180.*pow(E_NAD83,4.) + 517./5040.*pow(E_NAD83,6.);
   b =                         23./360.*pow(E_NAD83,4.) + 251./3780.*pow(E_NAD83,6.);
   c =                                                    761./45360.*pow(E_NAD83,6.);

   theta = atan2(x, (rho0 - y));

   lat = (beta + a*sin(2.0*beta) + b*sin(4.0*beta) + c*sin(6.0*beta))/D2R;
   lon = stand_lon + (theta/D2R)/nc;

   printf("known_lat = %f\n",lat);
   printf("known_lon = %f\n",lon);
   
   return 0;
}