/************************************************************************/
/************************************************************************/
/* Module PROJ.C- handles coordinate conversion (map projection)        */
/* functions included are:                                              */
/*	PROJ (projection number)                                        */
/*		called from main program                                */
/*		to choose a projection                                  */
/*	PASPEC (longitude,latitude,first parallel, second parallel)     */
/*		called from main program                                */
/*		to set the aspect of the projection.                    */
/*	PMAPLM (west,east,south,north)                                  */
/*		called from main program                                */
/*		to set the geographic area                              */
/*	PDRNGE (xmin,xmax,ymin,ymax,pixasp)                                    */
/*		called from main program                                */
/*		to set the equivalent plotting coordinate range.        */
/*	PSET(void)                                                      */
/*		called from MAP and UNMAP once after any of the         */
/*		above initialization calls                              */
/*		to compute the necessary scale and translation factors  */
/*	MAP(longitude,latitude,x,y)                                     */
/*		may be called to convert from                           */
/*		geographic coordinates to plotting coordinates          */
/*	UNMAP(longitude,latitude,x,y)                                   */
/*		may be called to convert from                           */
/*		plotting coordinates to geographic coordinates          */
/*                                                                      */
/************************************************************************/

#include <math.h>

#define PI 3.141593
#define DEGRAD PI/180.
#define RADDEG 180./PI

static int iproj;
static int aspect=0;
static double plon,plat,plat1,plat2;
static double rlonw,rlone,rlats,rlatn;
static double sx,sy,tx,ty;
static double dxmin,dxmax,dymin,dymax,pixasp;
static double n,F,rho0,B,k;
static int debug=0;
static int setflg=1;
static double R=1.;
static double k0=1.;
/**********************************************************************/
/* PROJ- sets the map projection number                               */
/* Argument- *dum1 is a number denoting the projection, 1=linear,     */
/*   2=Mercator, 3=Transverse Mercator, 4=Stereographic,              */
/*   5=Lambert Conformal Conic                                        */

int proj(int *dum1)
{
  iproj=*dum1; /* save the projection number */
  if(debug) printf("PROJ: %d \n",iproj);
  setflg=1;    /* setup will be performed on first call to MAP or UNMAP */
  aspect=0;    /* PSET will select a reasonable aspect */
  return iproj;
}

/********************************************************************/
/* PASPEC- sets the aspect and standard parallels for projection    */
/* Arguments- *dum1 is longitude, *dum2 is latitude (used in 3,4,5) */
/*   *dum3 and *dum4 are standard parallels for 5                   */

void paspec(float *dum1,float *dum2, float *dum3, float *dum4)
{
  plon=*dum1*DEGRAD;
  plat=*dum2*DEGRAD;
  plat1=*dum3*DEGRAD;
  plat2=*dum4*DEGRAD;
  aspect=1;
  setflg=1;
  return;
}

/********************************************************************/
/* PDRNGE sets the device coordinate range for the map area         */
/* Arguments- left,right,bottom,top plotting coordinate limits      */

void pdrnge(float *dum1, float *dum2, float *dum3, float *dum4,
  float *dum5 )
{
  dxmin=*dum1;
  dxmax=*dum2;
  dymin=*dum3;
  dymax=*dum4;
  pixasp=*dum5;
  setflg=1;
  if(debug) printf("PDRNGE: %lf %lf %lf %lf %lf\n",
    dxmin,dxmax,dymin,dymax,pixasp);
}

/********************************************************************/
/* PMAPLM- sets the geographic coordinate range for the map area    */
/* Arguments- West, East, South, and North geographic bounds        */

void pmaplm(float *dum1,float *dum2,float *dum3, float *dum4)
{
  rlonw=*dum1;
  rlone=*dum2;
  rlats=*dum3;
  rlatn=*dum4;
  setflg=1;
  return;
}
/***********************************************************************/
/*PSET- sets up scale and translation factors to fit geographic space  */
/*through the requested projection into the device space.              */
/*It is called by either MAP or UNMAP when the setflg variable is set  */
/*Method- Project each of NCHECK * NCHECK geographic points evenly     */
/*  spaced over the map area to determine the projected map extent.    */
/*  This method handles odd shaped projections very well.              */
/*  Scale and translation factors are computed to spread the projection */
/*  coordinates over the plotting system coordinates.                   */

void pset(void)
{
  #define LARGE 1000000.
  #define NCHECK 10.       /* gives 100 sample points */

  void map();
/* variables to receive the projection extent */
  double pxmax=-LARGE;
  double pxmin=LARGE;
  double pymax=-LARGE;
  double pymin=LARGE;
/* sample increments */
  double checkx,checky;
/* sample points */
  float rlon,rlat;
/* projected point */
  float x,y;

  setflg=0;

  if ( ! aspect )
  {
/* establish reasonable defualts for aspect and standard parallels */
    if(rlatn==90.) plat=90.*DEGRAD;
    else plat=(rlatn+rlats)/2.*DEGRAD;
    plon=(rlone+rlonw)/2.*DEGRAD;
    plat1=(rlats+(1./6.)*(rlatn-rlats))*DEGRAD;
    plat2=(rlats+(5./6.)*(rlatn-rlats))*DEGRAD;
    if(debug) printf("plat:%f plon:%f plat1:%f plat2:%f\n",
      plat,plon,plat1,plat2);
  }

  if(iproj==5)
  {
/*compute constants for Lambert Conformal Conic */
    n=log(cos(plat1)/cos(plat2))/
      log(tan(PI/4.+plat2/2.)/tan(PI/4.+plat1/2.));
    F=cos(plat1)*pow(tan(PI/4.+plat1/2.),n)/n;
    rho0=R*F/pow(tan(PI/4.+plat/2.),n);
    if(debug) printf("Lambert constants:\n:%lf F:%lf rho0:%lf\n",
      n,F,rho0);
  }

/* determine the project extent by sampling NCHECK*NCHECK points */
  checkx=(rlone-rlonw)/NCHECK;
  checky=(rlatn-rlats)/NCHECK;
  sx=1.;
  sy=1.;
  tx=ty=0.;                         /* reset translation factors */
  for (rlat=rlats;rlat<=rlatn;rlat=rlat+checky)
  {
    for(rlon=rlonw;rlon<=rlone;rlon=rlon+checkx)
    {
      map(&rlon,&rlat,&x,&y);      /* project the sample point */
      if(x>pxmax) pxmax=x;
      if(x<pxmin) pxmin=x;         /* testing for projection extent */
      if(y>pymax) pymax=y;
      if(y<pymin) pymin=y;
    }
  }
/* compute scales and translations to fit the projection extent to the
     plotting  coordinate range */
  sx=(dxmax-dxmin)/(pxmax-pxmin);
  sy=(dymax-dymin)/(pymax-pymin);
/* correct for nonsquare pixel */

/* select scale of smaller absolute magnitude to maintain the
map aspect ratio, the pixel aspect ratio must be considered */
  if(fabs(sx)<fabs(sy*pixasp))
  {
    if(sy<0.) sy=-fabs(sx/pixasp);
    else sy=fabs(sx/pixasp);
    if(debug) printf("sx is limiting, sx:%lf, sy:%lf\n",sx,sy);
  }
  else
  {
    if(sx<0.) sx=-fabs(sy*pixasp);
    else sx=fabs(sy*pixasp);
    if(debug) printf("sy is limiting, sx:%lf, sy:%lf\n",sx,sy);
  }
/*  sy=sy/1.2190; measured aspect ratio for EGA */
if(debug) printf("sx %lf sy %lf pixasp %lf\n",sx,sy,pixasp);
/* translation factors */
  tx=dxmin-(pxmin*sx);
  ty=dymin-(pymin*sy);
  setflg=0;
  if(debug) printf("PSET: %lf %lf %lf %lf %lf %lf %lf %lf\n",
    rlonw,rlone,rlats,rlatn,tx,ty,sx,sy);
  return;
}

/**********************************************************************/
/* MAP- computes plotting coordinates from geographic coordinates     */
/* Inputs- rlond,rlatd- geo cordinates in degrees                     */
/* Outputs- x,y- plotting coordinates                                 */

void map(float *rlond, float *rlatd, float *x,float *y)
{
/*projected coordinates */
  double px,py;
/* geo coordinates in radians */
  double rlonr,rlatr;
/* temporary variables */
  double rho,theta;
/* if any parameters have been changes, recompute scales and translations */
  if(setflg) pset();
  rlonr=*rlond*DEGRAD;
  rlatr=*rlatd*DEGRAD;

  switch(iproj)
  {
    case 1 :
    {
/* Linear or Cylindrical Equidistant */
      px=rlonr;
      py=rlatr;
      break;
    }
    case 2 :
    {
/* Mercator- ref 1, page 47 */
      px=R*(rlonr-plon);
      py=R*log(tan(PI/4.+rlatr/2.));
      break;
    }
    case 3 :
/* Transverse Mercator- Ref 1, page 67 */
    {
      B=cos(rlatr)*sin(rlonr-plon);
      px=0.5*R*log((1.+B)/(1.-B));
      py=R*(atan(tan(rlatr)/cos(rlonr-plon))-plat);
      break;
    }
    case 4 :
    {
/* Stereographic- Ref 1, page 158 */
      k=2./(1.+sin(plat)*sin(rlatr)+
	cos(plat)*cos(rlatr)*cos(rlonr-plon));
      px=R*k*cos(rlatr)*sin(rlonr-plon);
      py=R*k*(cos(plat)*sin(rlatr)-sin(plat)*cos(rlatr)*cos(rlonr-plon));
      break;
    }
    case 5 :
    {
/* Lambert Conformal Conic- Ref 1, page 105 */
      if((rlonr-plon)>PI) rlonr=rlonr-PI;
      if((rlonr-plon)<-PI) rlonr=rlonr+PI;
      theta=n*(rlonr-plon);
      if(*rlatd==90.) rho=0.;
      else rho=R*F/pow(tan(PI/4.+rlatr/2.),n);
      px=rho*sin(theta);
      py=rho0-rho*cos(theta);
      break;
    }
  }
/* compute plotting coordinates from projection coordinates */
  *x=px*sx+tx;
  *y=py*sy+ty;
  if(debug) printf("MAP: %f %f %f %f\n",*x,*y,*rlond,*rlatd);
  return;
}
/***********************************************************************/
/*UNMAP- computes geographic coordinates from plotting coordinates     */
/* Inputs- x,y in plotting coordinates                                 */
/* Outputs- rlond,rlatd - geo coordinates in degrees                   */
/* Returns 0 if point within geo area, else 1                          */

int unmap (float *rlond, float *rlatd,float *x, float *y)
{
/* projection coordinates */
  double px,py;
/* temporary variables */
  double c,rho,D,theta;
/* geo coords in radians */
  double rlonr,rlatr;
/* if any parameters have changes, compute new scale and translations */
  if(setflg) pset();
/* compute projection coords from plotting coords */
  px=(*x-tx)/sx;
  py=(*y-ty)/sy;
/* execute the inverse projection code */
  switch(iproj)
  {
    case 1 :
    {
/* Linear */
      rlonr=px;
      rlatr=py;
      break;
    }
    case 2 :
/* inverse of Mercator Ref 1, page 50 */
    {
      rlatr=PI/2.-2.*atan(pow(2.7182818,-py/R));
      rlonr=px/R+plon;
      break;
    }
    case 3 :
/* inverse of Transverse Mercator Ref 1, page 67 */
    {
      D=py/(R*k0)+plat;
      rlatr=asin(sin(D)/cosh(px/R*k0));
      rlonr=plon+atan2(sinh(px/R*k0),cos(D));
      break;
    }
    case 4 :
/* inverse of Stereographic Ref 1, page 159 */
    {
      rho=sqrt(px*px+py*py);
      c=2.*atan2(rho,(2*R*k0));
      if(rho==0.) rlatr=plon;
      else
	rlatr=asin(cos(c)*sin(plat)+(py*sin(c)*cos(plat)/rho));
      if(plat==90.)
	rlonr=plon+atan2(px,(-py));
      else
	if(plat==-90.)
	  rlonr=plon+atan2(px,py);
	else
	  rlonr=plon+atan2(px*sin(c),
	    (rho*cos(plat)*cos(c)-py*sin(plat)*sin(c)));
      break;
    }
    case 5 :
/* inverse of Lambert Conformal Conic Ref 1 page 105 */
    {
      rho=sqrt(px*px+(rho0-py)*(rho0-py));
      if(n<0.) rho=-rho;
      theta=atan2(px,(rho0-py));
      if(rho==0.)
      {
	if(n<0.) rlatr=-PI/2.;
	else rlatr=PI/2.;
      }
      else rlatr=2.*pow(atan2(R*F,rho),1./n)-PI/2;
      rlonr=theta/n+plon;
      break;
    }
  }
  *rlond=rlonr*RADDEG;
  *rlatd=rlatr*RADDEG;
  if(debug) printf("UNMAP: %f %f %f %f\n",*x,*y,*rlond,*rlatd);
  if (*rlond<rlonw || *rlond>rlone || *rlatd<rlats ||*rlatd>rlatn)
    return 1;
    else return 0;

}