#include "gwin.user.h"
#include <math.h>
double myatan2();
float t0[3],t1[3],t2[3],ti0[3],ti1[3],ti2[3];
double pi,pi1, pi2, pi3;
float xmerc();
float lats[500],longs[500];
int ortho;
int merc;
int setmercator(),setortho(),plot_grid(),interrupt_display();
int new_center(),redisplay(),setortho2();
float x,y;
int event;
char key;
int interrupt = 0;
int display_underway = 0;
float lat_cen,long_cen;
int ortho2;

main()
{
float xlat,xlong;
char line[30];

   pi1 = 3.1415926535897943/180.;
   pi2 = 180./3.1415926535897943;
   pi3 = 3.1415926535897943/2.0;

   ortho = 1;
   ortho2 = 0;
   merc = 0;

   printf("\n\nEnter lat_cen, long_cen:  ");
   scanf("%f %f",&lat_cen,&long_cen);

   ustart("high2",0.,640.,0.,400.);
   if(merc ) uwindo(-180.,180.,-180.,180.);
   if(ortho || ortho2) uwindo(-250.,250.,-160.,160.);

   setup_menus();

   make_rotation_matrices(lat_cen,long_cen);

   make_geo_grid();

   plot_grid();

   while(1 == 1) {
      ugrina(&x,&y,&event,&key);
      latlong(x,y,&xlat,&xlong);

      sprintf(line,"%6.2f  %6.2f",xlat,xlong);
      if (ortho || ortho2){
         uprint(100.,150.,line);
      }else{
         uprint(100.,170.,line);
      }

   }
}

latlong(x1,y1,xlat,xlong)
float x1,y1,*xlat,*xlong;
{
float xxx,yyy,zzz;
double xtemp,ytemp,ztemp;
float temp1,temp2;
float ph1,th1;
char line[30];

   if(ortho || ortho2){
      yyy = x1/150.;
      zzz = y1/150.;

      temp1 = yyy*yyy;
      temp2 = zzz*zzz;
      if(temp1 + temp2 > 1.0)return(1);

      xxx = sqrt((double)(1 - temp1 - temp2));

      xtemp = (double)(ti0[0] * xxx + ti0[1] * yyy + ti0[2] * zzz);
      ytemp = (double)(ti1[0] * xxx + ti1[1] * yyy + ti1[2] * zzz);
      ztemp = (double)(ti2[0] * xxx + ti2[1] * yyy + ti2[2] * zzz);

      if(fabs(xtemp) < 1e-15) xtemp = 1e-15;
      *xlong = pi2*myatan2(ytemp,xtemp);

      if(ztemp >  1.0) ztemp =  1.0;
      if(ztemp < -1.0) ztemp = -1.0;
      *xlat = pi2*(pi3 - acos(ztemp));

   }

   if(merc){

      ph1 = pi1 * y1;
      ph1 = 2.0 * atan(exp(ph1)) - pi3;
      ph1 = pi3 - ph1;

      th1 = pi1 * x1;

      xxx = sin((double)ph1) * cos((double)th1);
      yyy = sin((double)ph1) * sin((double)th1);
      zzz = cos((double)ph1);

      xtemp = (double)(ti0[0] * xxx + ti0[1] * yyy + ti0[2] * zzz);
      ytemp = (double)(ti1[0] * xxx + ti1[1] * yyy + ti1[2] * zzz);
      ztemp = (double)(ti2[0] * xxx + ti2[1] * yyy + ti2[2] * zzz);

      if(fabs(xtemp) < 1e-15) xtemp = 1e-15;
      *xlong = pi2*myatan2(ytemp,xtemp);

      if(ztemp >  1.0) ztemp =  1.0;
      if(ztemp < -1.0) ztemp = -1.0;
      *xlat = pi2*(pi3 - acos(ztemp));


   }
}

make_rotation_matrices(lat_cen,long_cen)
float lat_cen,long_cen;
{

/* Rotation matrices to rotate lat/long to new lat/long */
/* so that (lat_cen,long_cen) becomes (0,0)             */

float sin_theta, cos_theta, theta_rot, sin_phi, cos_phi, phi_rot;

theta_rot = long_cen * pi1;
phi_rot   = lat_cen  * pi1;

sin_theta = sin(theta_rot);
cos_theta = cos(theta_rot);
sin_phi   = sin(phi_rot);
cos_phi   = cos(phi_rot);

t0[0] = cos_theta * cos_phi;
t0[1] = sin_theta * cos_phi;
t0[2] = sin_phi;
t1[0] = -sin_theta;
t1[1] = cos_theta;
t1[2] = 0.0;
t2[0] = -cos_theta * sin_phi;
t2[1] = -sin_theta * sin_phi;
t2[2] = cos_phi;

ti0[0] = cos_theta * cos_phi;
ti0[1] = -sin_theta;
ti0[2] = -cos_theta * sin_phi;
ti1[0] = sin_theta * cos_phi;
ti1[1] = cos_theta;
ti1[2] = -sin_theta * sin_phi;
ti2[0] = sin_phi;
ti2[1] = 0.0;
ti2[2] = cos_phi;

}

mercator(long_in,lat_in,x_out,y_out)
float lat_in,long_in,*y_out,*x_out;
{
float xla,xlo,x[3],xr[3],sth,cth,sph,cph;
double th1,ph1;

   xla = lat_in;
   xlo = long_in;
   th1 = pi1 * xlo;
   ph1 = pi3 - xla * pi1;

   sth = sin(th1);
   cth = cos(th1);
   sph = sin(ph1);
   cph = cos(ph1);

   x[0] = sph * cth;
   x[1] = sph * sth;
   x[2] = cph;

   xr[0] = t0[0] * x[0] + t0[1] * x[1] + t0[2] * x[2];
   xr[1] = t1[0] * x[0] + t1[1] * x[1] + t1[2] * x[2];
   xr[2] = t2[0] * x[0] + t2[1] * x[1] + t2[2] * x[2];

   *y_out = xmerc( pi2 * asin(xr[2]));
   *x_out = pi2 * myatan2(xr[1],xr[0]);
}

float xmerc(xlat)
float(xlat);
{
double a1;

   if (fabs(xlat) < 89.9999){
      if(xlat >= 0.0){
         a1 = pi1 * (45.0 + .5 * xlat);
         return((float) pi2 * log(tan(a1)));
      } else {
         a1 = -pi1 * (-45.0 + .5 * xlat);
         return((float) -pi2 * log(tan(a1)));
      }
   } else {

      if(xlat < 0) return(-750.);
      return(750.);
   }
}

orthographic(long_in,lat_in,x_out,y_out,frontside)
float lat_in,long_in,*y_out,*x_out;
int *frontside;
{
float xla,xlo,x[3],xr[3],sth,cth,sph,cph;
double th1,ph1;

   xla = lat_in;
   xlo = long_in;
   th1 = pi1 * xlo;
   ph1 = pi3 - xla * pi1;

   sth = sin(th1);
   cth = cos(th1);
   sph = sin(ph1);
   cph = cos(ph1);

   x[0] = sph * cth;
   x[1] = sph * sth;
   x[2] = cph;

   xr[0] = t0[0] * x[0] + t0[1] * x[1] + t0[2] * x[2];
   xr[1] = t1[0] * x[0] + t1[1] * x[1] + t1[2] * x[2];
   xr[2] = t2[0] * x[0] + t2[1] * x[1] + t2[2] * x[2];

   *y_out = 150.0 * xr[2];
   *x_out = 150.0 * xr[1];
   *frontside = 0;
   if(xr[0]>0.0)*frontside = 1;

}

make_geo_grid()
{
int i;
float x;

   i = 0;

   for(x = -90.0;x < 91.0 ;x += 10.0){
      lats[i++] = x;
   }
   lats[0] = -89.999;
   lats[i-1] = 89.999;
   i = 0;

   for(x = -180.0; x < 181.0; x += 10.0){
      longs[i++] = x;
   }
   longs[0] = -179.999;
   longs[i-1] = 179.999;
}

plot_grid()
{
   if(merc){
      mercplot();
      return(0);
   }
   if(ortho){
      orthoplot();
      return(0);
   }
   if(ortho2){
      ortho2plot();
      return(0);
   }
}

mercplot()
{
int i,j;
float x,y,xold;
int event,key;

   display_underway = 1;
   interrupt = 0;

   uerase();

   upset("colo",1.0);

   for(i=0;i<19;i++){
      mercator(longs[0],lats[i],&x,&y);
      umove(x,y);
      xold = x;
      for(j=1;j<37;j++){
         ugrina(&x,&y,&event,&key);
         if(interrupt) goto interrupted;
         mercator(longs[j],lats[i],&x,&y);
         if(fabs(xold-x) < 100.0)udraw(x,y);
         umove(x,y);
         xold = x;
      }
   }
   for(j=0;j<37;j++){
      mercator(longs[j],lats[0],&x,&y);
      umove(x,y);
      xold = x;
      for(i=1;i<19;i++){
         ugrina(&x,&y,&event,&key);
         if(interrupt) goto interrupted;
         mercator(longs[j],lats[i],&x,&y);
         if(fabs(xold-x) < 100.0)udraw(x,y);
         umove(x,y);
         xold = x;
      }
   }
interrupted:
   interrupt = 0;
   display_underway = 0;
}

orthoplot()
{
int i,j;
float x,y,xold,yold;
int event,key;
int frontside;

   display_underway = 1;
   interrupt = 0;

   uerase();

   upset("colo",1.0);

   for(i=0;i<19;i++){
      orthographic(longs[0],lats[i],&x,&y,&frontside);
      umove(x,y);
      xold = x;
      yold = y;
      for(j=1;j<37;j++){
         ugrina(&x,&y,&event,&key);
         if(interrupt) goto interrupted;
         orthographic(longs[j],lats[i],&x,&y,&frontside);
         if(fabs(xold-x) < 100.0 && frontside){
            umove(xold,yold);
            udraw(x,y);
         }else{
            umove(x,y);
         }
         xold = x;
         yold = y;
      }
   }
   for(j=0;j<37;j++){
      orthographic(longs[j],lats[0],&x,&y,&frontside);
      umove(x,y);
      xold = x;
      yold = y;
      for(i=1;i<19;i++){
         ugrina(&x,&y,&event,&key);
         if(interrupt) goto interrupted;
         orthographic(longs[j],lats[i],&x,&y,&frontside);
         if(fabs(xold-x) < 100.0 && frontside){
            umove(xold,yold);
            udraw(x,y);
         }else{
            umove(x,y);
         }
         xold = x;
         yold = y;
      }
   }
interrupted:
   interrupt = 0;
   display_underway = 0;
}

ortho2plot()
{
int i,j;
float x,y,xold,yold;
int event,key;
int frontside;

   frontside = 1;
   display_underway = 1;
   interrupt = 0;

   uerase();

/* plot backsides first */

   upset("colo",10.0);

   for(i=0;i<19;i++){
      orthographic(longs[0],lats[i],&x,&y,&frontside);
      umove(x,y);
      xold = x;
      yold = y;
      for(j=1;j<37;j++){
         ugrina(&x,&y,&event,&key);
         if(interrupt) goto interrupted;
         orthographic(longs[j],lats[i],&x,&y,&frontside);
         if(fabs(xold-x) < 100.0 && !frontside){
            umove(xold,yold);
            udraw(x,y);
         }else{
            umove(x,y);
         }
         xold = x;
         yold = y;
      }
   }
   for(j=0;j<37;j++){
      orthographic(longs[j],lats[0],&x,&y,&frontside);
      umove(x,y);
      xold = x;
      yold = y;
      for(i=1;i<19;i++){
         ugrina(&x,&y,&event,&key);
         if(interrupt) goto interrupted;
         orthographic(longs[j],lats[i],&x,&y,&frontside);
         if(fabs(xold-x) < 100.0 && !frontside){
            umove(xold,yold);
            udraw(x,y);
         }else{
            umove(x,y);
         }
         xold = x;
         yold = y;
      }
   }

/* now plot front sides */

   upset("colo",1.0);

   for(i=0;i<19;i++){
      orthographic(longs[0],lats[i],&x,&y,&frontside);
      umove(x,y);
      xold = x;
      yold = y;
      for(j=1;j<37;j++){
         ugrina(&x,&y,&event,&key);
         if(interrupt) goto interrupted;
         orthographic(longs[j],lats[i],&x,&y,&frontside);
         if(fabs(xold-x) < 100.0 && frontside){
            umove(xold,yold);
            udraw(x,y);
         }else{
            umove(x,y);
         }
         xold = x;
         yold = y;
      }
   }
   for(j=0;j<37;j++){
      orthographic(longs[j],lats[0],&x,&y,&frontside);
      umove(x,y);
      xold = x;
      yold = y;
      for(i=1;i<19;i++){
         ugrina(&x,&y,&event,&key);
         if(interrupt) goto interrupted;
         orthographic(longs[j],lats[i],&x,&y,&frontside);
         if(fabs(xold-x) < 100.0 && frontside){
            umove(xold,yold);
            udraw(x,y);
         }else{
            umove(x,y);
         }
         xold = x;
         yold = y;
      }
   }
interrupted:
   interrupt = 0;
   display_underway = 0;
}

setup_menus()
{
   uamenu(1,0,0,"PROJECTION      ",' ',0,(USHORT)MIDRAWN|MENUENABLED,7);
   uamenu(1,1,0,"   ORTHOGRAPHIC ",'O',6,(USHORT)MIDRAWN|ITEMTEXT|HIGHCOMP
          |COMMSEQ|ITEMENABLED|CHECKIT,setortho);
   uamenu(1,2,0,"   MERCATOR     ",'M',5,(USHORT)MIDRAWN|ITEMTEXT|HIGHCOMP
          |COMMSEQ|ITEMENABLED,setmercator);
   uamenu(1,3,0,"   ORTHO2       ",'2',3,(USHORT)MIDRAWN|ITEMTEXT|HIGHCOMP
          |COMMSEQ|ITEMENABLED,setortho2);

   uamenu(2,0,0,"DISPLAY    ",' ',0,(USHORT)MIDRAWN|MENUENABLED,7);
   uamenu(2,1,0,"REDISPLAY  ",'R',0,(USHORT)MIDRAWN|ITEMTEXT|HIGHCOMP
          |COMMSEQ|ITEMENABLED,redisplay);
   uamenu(2,2,0,"INTERRUPT  ",'I',0,(USHORT)MIDRAWN|ITEMTEXT|HIGHCOMP
          |COMMSEQ|ITEMENABLED,interrupt_display);
   uamenu(2,3,0,"NEW CENTER ",'C',0,(USHORT)MIDRAWN|ITEMTEXT|HIGHCOMP
          |COMMSEQ|ITEMENABLED,new_center);

}

redisplay()
{
   plot_grid();
   interrupt=1;  /* cancel previous invocation if any */
}

setortho()
{
   uwindo(-250.,250.,-160.,160.);
   ortho = 1;
   merc = 0;
   ortho2 = 0;
   uamenu(1,1,0,"   ORTHOGRAPHIC ",'O',6,(USHORT)MIDRAWN|ITEMTEXT|HIGHCOMP
          |COMMSEQ|ITEMENABLED|CHECKIT|CHECKED,setortho);
   redisplay();
}

setmercator()
{
   uwindo(-180.,180.,-180.,180.);
   merc = 1;
   ortho = 0;
   ortho2 = 0;
   uamenu(1,2,0,"   MERCATOR     ",'M',5,(USHORT)MIDRAWN|ITEMTEXT|HIGHCOMP
          |COMMSEQ|ITEMENABLED|CHECKIT|CHECKED,setmercator);
   redisplay();
}

setortho2()
{
   uwindo(-250.,250.,-160.,160.);
   ortho2 = 1;
   ortho = 0;
   merc = 0;
   uamenu(1,3,0,"   ORTHO2       ",'2',3,(USHORT)MIDRAWN|ITEMTEXT|HIGHCOMP
          |COMMSEQ|ITEMENABLED|CHECKIT|CHECKED,setortho2);
   redisplay();
}

interrupt_display()
{
    interrupt = 1;
}

new_center()
{
char coords[255];

   uerase();
   uwindo(0.,100.,0.,100.);
   ugetstring(10.,75.,70.,"Enter latitude and longitude:  ",coords);

   while (sscanf(coords,"%f %f",&lat_cen,&long_cen) < 2){
      ugetstring(10.,75.,70.,
         "Error, try again.  Enter latitude and longitude:  ",coords);
   }
   make_rotation_matrices(lat_cen,long_cen);

   if(merc)setmercator();
   if(ortho)setortho();
   if(ortho2)setortho2();
}

double myatan2(y,x)
double x,y;
{
double temp;
double pi = 3.1415926535897943;
double piov2 = 0.5 * 3.1415926535897943;

   if(fabs(x) > 0.0){
      temp = atan((double)y/x);
   }else{
      if(y < 0.0){
         temp = -piov2;
      }else{
         temp =  piov2;
         if(fabs(x) < 1e-15 && fabs(y) <1e-15) temp = 0.0;
      }
   }

   if(x < 0.0 ){
      if(y >= 0.0) {
         temp =   pi + temp;
      }else{
         temp =  -pi + temp;
      }
   }
   return(temp);
}