/* *******************************************************************
 *                                                                   *
 * Copyright (c) L-DGO/MIT/JGOFS                                     *
 *                                                                   *
 *                                                                   *
 * File             : remap.c                                        *
 *                                                                   *
 * Purpose          :                                                *         
 *                                                                   *
 * Version Number   : 1.1                                            *
 *                                                                   *
 * Revision History :                                                *
 *                                                                   *
 * Date                              Developer                       *
 * ----                              ---------                       *
 *                                                                   *
 * Sat Oct 17  1992 10.00           Glenn Flierl                     *
 *                                                                   *   
 *                                                                   *
 ******************************************************************* */

/*
 main entry points are
procedure getinfo;
  reads map.var and sets up constants
  useful global variables are  maptype,lat0,lat1,long0,mag
procedure remap(maptype:integer;lambda,phi:real; var x,y:real);
  transforms lat,lon (phi,lambda) to x,y
  either -1.5<=x<=1.5 or -1<=y<=1 depending on projection
*/
#include <math.h>
#include <stdio.h>
#include <string.h>
#ifdef PC
#include <stdlib.h>
#else
double atof();
#endif

#define sqrt2 1.4142135623731
#define pi 3.1415926535898
#define halfpi 1.5707963267949
#define twopi 6.2831853071796
#define radian 1.7453292519943e-2
#define radiandiv100 1.7453292519943e-4
#define numproj 12

float notvisible= -10.0;   /* flag for point visibility */
float maxlat= 1.397;   /* ~80 degrees. */
float lambmap[22] ={
1,0.9981044,0.9924095,0.9829041,0.9695689,0.952375,0.9312824,
0.9062371,0.8771676,0.8439795,0.8065483,0.764708,0.7182353,
0.6668243,0.6100447,0.5472709,0.4775446,0.3992848,0.3095366,
0.2011957,0.0,0.0};

struct mapitem{
  char mapname[40];
  int (*mapcommand)();
  };

float lat0,long0,lat1,mag,dely;
int my_quad,oldquad;
char proj[40];
char mapfile[80];

float arctanh(x)
float x;
{
float a,t;
t=fabs(x);
if (t>= 1) a=1000;
else a = 0.5 * log((1.0 + t)/(1.0 - t));
if(x<0)return -a;
return a;
}  /* ArcTanH. */

float meridian(lambda, lambda0)
float lambda,lambda0;
/* returns difference between current longitude and map center. */
{
float dellam;
dellam = lambda - lambda0;
if (dellam < -pi) dellam = dellam + twopi;
else if (dellam >  pi) dellam = dellam - twopi;
return dellam;
}   /* meridian. */

lambert(lambda,phi,x,y)
float lambda,phi,*x,*y;
{
int k;
*x=fabs(phi)*40.0/pi;
k = *x;
*x = *x-k;
*y = lambmap[k]*(1.0- *x)+ *x * lambmap[k+1];
*x = *y * lambda*1.5/pi;
*y = sqrt(1.0 - *y * *y)/0.614022403*1.5/pi;
if (phi<0.0) *y = - *y;
}

mercator(lambda, phi, x, y)
float lambda,phi,*x,*y;
{
if (fabs(phi) < maxlat){
  *x = lambda/pi*1.5;
  *y = arctanh(sin(phi))/pi*1.5;
  }
else
  {*x = notvisible;*y=0;};
}   /* mercator. */

equicyl(lambda, phi, x, y )
float lambda,phi,*x,*y;
{
*x = lambda * cos(lat1)/pi*1.5;
*y = phi/pi*1.5;
}   /* equicyl. */

sinusoidal(lambda, phi, x, y )
float lambda,phi,*x,*y;
{
*x = cos(phi) * lambda/pi*1.5;
*y = phi/pi*1.5;
}  /* sinusoidal. */

hammer(lambda, phi, x, y)
float lambda,phi,*x,*y;
{
float k, cosphi, halflambda;
halflambda = 0.5*lambda;
cosphi=cos(phi);
k = 0.5 / sqrt(1 +cosphi * cos(halflambda));
*x = 2 * k * cosphi * (sin(halflambda))*1.5;
*y = k * sin(phi)*1.5;
}  /* hammer. */

orthographic(lambda, phi,x, y)
float lambda,phi,*x,*y;
{
float cosl, sinphi, cosphi;
cosphi =cos(phi);    sinphi= sin(phi);
cosl   =cos(lambda)*cosphi;
if (cosl >= 0){
  *x = cosphi * sin(lambda);
  *y = sinphi;
  }
else 
  { *x=notvisible;*y=0;};
}  /* orthographic. */

gnomic(lambda,phi, x,y)
float lambda,phi,*x,*y;
{
float cl;
cl=cos(lambda);
if (cl<=0 || fabs(phi)>maxlat){
  *x=notvisible;
  *y=0;
  }
else
  {
  *x=sin(lambda)/cl/3.0;
  *y=sin(phi)/cos(phi)/cl/3.0;
  };
}

areacyl(lambda,phi, x,y)
float lambda,phi,*x,*y;
{
*x=lambda/pi*1.5;
*y=sin(phi)/pi*1.5;
}

peters(lambda,phi, x,y)
float lambda,phi,*x,*y;
{
*x=lambda/pi*1.5;
*y=sin(phi)/pi*3.0;
}

stereographic(lambda,phi, x,y)
float lambda,phi,*x,*y;
{
float ct,st,cr,sr;
ct=cos(phi);cr=cos(lambda);
if (ct*cr<-0.99){
  *x=notvisible;*y=0;
  }
else
  {
  *x=2*ct*sin(lambda)/(1+ct*cr)/5.0;
  *y=2*sin(phi)/(1+ct*cr)/5.0;
  };
}

equiplanar(lambda,phi,x,y)
float lambda,phi,*x,*y;
{
float ct,cph;
ct=cos(phi);cph=ct*cos(lambda);
cph=acos(cph);
if (cph<0) cph=pi+cph;
*x=ct*sin(lambda);
*y=sin(phi);
if (*x==0.0 && *y==0.0) cph=0.0; else cph=cph/pi/sqrt(*x * *x+ *y * *y);
*x=cph* *x;*y=cph* *y;
}

areaplanar(lambda,phi, x,y)
float lambda,phi,*x,*y;
{
float ct,cph;
ct=cos(phi);cph=ct*cos(lambda);
*x=ct*sin(lambda);
*y=sin(phi);
if (*x==0.0 && *y==0.0) cph=0.0; else cph=sqrt((2-2*cph)/(*x * *x+ *y
   * *y))/2.0;
*x=cph* *x;*y=cph* *y;
}

struct mapitem maplist[numproj]=
{"mercator",mercator,
 "lambert",lambert,
 "equicylinder",equicyl,
 "sinusoidal",sinusoidal,
 "hammer",hammer,
 "orthographic",orthographic,
 "gnomic",gnomic,
 "areacylinder",areacyl,
 "stereographic",stereographic,
 "equiplanar",equiplanar,
 "areaplanar",areaplanar,
 "peters",peters,
};
int (*maptype)();

remap(lambda,phi, x,y)
float lambda,phi,*x,*y;
{
float st0,st1,ct0,ct1,sl,cl,a1,b1,c1;
/*write(phi,' ',lambda,' ');*/
lambda=meridian(radian*lambda,long0);
phi=phi*radian;
if (lat0 != 0){
  if (fabs(lat0)==90){
    st0=lat0/fabs(lat0);ct0=0;
    }
  else
    {
    st0=sin(lat0);ct0=cos(lat0);
    };
  st1=sin(phi);ct1=cos(phi);
  cl=cos(lambda);sl=sin(lambda);
  a1=st0*st1+ct0*ct1*cl;
  b1=ct0*st1-ct1*st0*cl;
  c1=ct1*sl;
  if (phi != 0.0 || fabs(lat0)!=90)
    {
    phi=asin(b1);
    lambda=asin(c1/sqrt(fabs(1.0-b1*b1+1.0e-10)));
    }
  else
    {
/*    write(phi,' ',lambda,' ');*/
    if (lat0==90) phi=fabs(lambda)-halfpi;else phi=halfpi-fabs(lambda);
    if (lambda>0) lambda=halfpi;else lambda = -halfpi;
    a1=0.0;
/*    writeln(phi,' ',lambda);*/
    };
  if (a1<0)
    {
    if (c1<0) lambda= -pi-lambda;else lambda=pi-lambda;
    };
  };
  if (lambda>halfpi) my_quad=1;
  else if (lambda < -halfpi) my_quad = -1;
  else my_quad=0;
  if (my_quad*oldquad == -1) newquad();
  oldquad=my_quad;

/*writeln(phi,' ',lambda);*/
(*maptype)(lambda,phi,x,y);
*x = *x * mag;
*y = *y * mag - dely;
}

getinfo()
{
float xt,yt;
int i;
char s[81],*t,*v;
FILE *fil;

maptype=maplist[0].mapcommand;
if((fil=fopen("map.var","r")) == NULL){
  printf("&x error - Missing map.var\n");
  exit(1);
};
while (fgets(s,81,fil) != NULL) 
  if (s[0] !='#'){
    t=strtok(s,"=");
    v=strtok(NULL," #\n");
    if(!strcmp(t,"projection")){
      for(i=0;i<numproj;i++) if(!strcmp(maplist[i].mapname,v)){
        maptype=maplist[i].mapcommand;
	break;
	};
      }
    else if (!strcmp(t,"cent_lat"))lat0=atof(v);
    else if (!strcmp(t,"cent_lon"))long0=atof(v);
    else if (!strcmp(t,"shift_lat"))lat1=atof(v);
    else if (!strcmp(t,"mag"))mag=atof(v);
  };
fclose(fil);
if (fabs(lat0)<90) lat0 *= radian;
long0 *= radian;
dely=0.0;
if (lat1 != 0.0){remap(long0,lat1,&xt,&yt);dely=yt;};
}