/* Rotate and display a rotating globe correctly from any angle */ 
/* array begins at 1 and ends at N */
/* Uses all three Euler angles. Pretty cool  */
/* plane of projection is XY */
/* set up so that latitude-dependent rotation rate is possible */
/* J.R.Schmidt                             */
/* gcc -DXWIN -O3 -o globe g2_globe.c -L/usr/X11R6/lib -lX11 -lm -lgd -lg2 */
/* or on WIN32 platforms */
/* gcc -static -O2 -DPFD_DOUBLEBUFFER=1 -DTBSTYLE_EX_DOUBLEBUFFER=1 -DUNIX=1 -DX_DISPLAY_MISSING=1 -DDO_WIN32=1 -DHAVE_LIMITS_H=1 -DWIN32 g2_globe.c -I/usr/local/cross-tools/include/ -I/usr/local/cross-tools/i386-mingw32msvc/include/ -L/usr/local/cross-tools/lib/ -L/usr/local/cross-tools/i386-mingw32msvc/lib/ -lg2 -lm -lkernel32 -lgdi32 -lmsvcrt -luser32 -lcomdlg32 -o g2_globe.exe */
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<g2.h>
#ifdef WIN32
#include<g2_win32.h>
#define WINTYPE 0
#endif
#ifdef XWIN
#include<g2_X11.h>
#define WINTYPE 1
#endif
#define N 1024    /* this is number of plaquettes to be drawn */
#define Np 32 /* number of azimuthal divisions */
#define Nt 33 /* number of latitude divisions */
#define PI 3.1415926

struct plaquette {
    /* Each little "square" on the globe is called a "plaquette" */
    /* and these are given attributes of color, orientation, distance */
    /* from the viewer vantage point, and so forth. Additional attributes */
    /* are easy to construct, such as tangential velocity for a rotating */
    /* globe, radial velocity, temperature, and so forth. There is a routine */
    /* "getplaquetteattribs" where they are loaded */

  double vertex[4][3];  /* x and y coords of basepoint and other 3 points */
  double center[3];     /* x and y coordinates of its center */
  double norml[3];     /* normal vector to plaquette face */

  double latitude;
  double face;
  double vdist;         /* distance from vantage point */
  int clr;             /* grayscale shade             */
  double velocity[3]; /* velocity vector (rotating globe) */
} pl[N+2],tempPL;


double dist(double x1, double y1, double z1, double x2, double y2, double z2);
void getplaquettecenter( struct plaquette *Plaq);
void getplaquetteattribs( struct plaquette *Plaq, double u, double v, double w);
void getplaquettelatitude( struct plaquette *Plaq, double AXIS_x, double AXIS_y, double AXIS_z);
void getplaquettevelocity( struct plaquette *Plaq, double AXIS_x, double AXIS_y, double AXIS_z, double OMEGA);
/* use this function to control oblateness */
double function(double x, double y);
int mod(int A, int B);

double dphi,dtheta,Stheta,Sphi,Srad, vantx,vanty,vantz,R, omega, t, dt;
/* components of rotation matrix */
double a11,a12,a13,a21,a22,a23,a31,a32,a33;
double pts[8];
double Z[Np*Nt+2],X[Np*Nt+2],Y[Np*Nt+2];
double AXIS[3]; /* rotation axis */
int n,m,CL,i,j,l,ir,k,u,u1,u2,u3,TIM;
double dimensionx,vantagex, vantagey,vantagez;/* position of a vantage point */ 
double theta,phi,psi,tmpx,tmpy,tmpz;
int dev,dev1,white, black, yellow, red, color;


main(int argc, char *argv[]){
dtheta=PI/(double)Nt;
dphi=2.0*PI/(double)Np;
vantagex=0.0;
vantagey=0.0;
vantagez=5.0;
 if(argc < 4 || argc > 4){
     printf("./globe theta phi psi\n");
     printf("All three Euler angles in radians. Suggestion; 1.9 0.0 0.0\n");
     exit(1);}
theta=atof(argv[1]);
phi=atof(argv[2]);
psi=atof(argv[3]);

 dt=0.05;
 omega=0.2;   /* axial rotation rate (uniform) */
 AXIS[0]=0.0, AXIS[1]=0.0, AXIS[3]=1.0;

  /* create a spherical grid and get the scale */ 
 X[0]=0.0;Y[0]=0.0;Z[0]=function(0.0,0.0); /* north pole */
  dimensionx=function(0.0,0.0);
  for(n=1;n<=Nt;n++){
    for(m=0;m<=Np;m++){
      Sphi=(double)m*dphi;
      Stheta=(double)n*dtheta;
      /* now remap the surface */
      R=function(Sphi,Stheta);
      if(R>dimensionx) dimensionx=1.2*R;
      X[Np*(n-1)+m+1]=R*sin((double)n*dtheta)*cos((double)m*dphi);
      Y[Np*(n-1)+m+1]=R*sin((double)n*dtheta)*sin((double)m*dphi);
      Z[Np*(n-1)+m+1]=R*cos((double)n*dtheta);}
  }

      dev=g2_open_vd();
      if(WINTYPE==1)
	  dev1=g2_open_X11(600,600);
      if(WINTYPE==0)
	  dev1=g2_open_win32(600,600,"Globe",0);

      g2_attach(dev,dev1);
      g2_set_auto_flush(dev1,0);
      white=g2_ink(dev1, 1.0, 1.0, 1.0);
      black=g2_ink(dev1, 0.0, 0.0, 0.0);	  
      red=g2_ink(dev1, 0.8, 0.0, 0.0);


/* when drawing the plaquettes, fill them in with pl[n].clr */
/* construct axis of rotation */
AXIS[0]=sin(psi)*sin(theta);
AXIS[1]=cos(psi)*sin(theta);
AXIS[2]=cos(theta);

 for(TIM=0;TIM<500;TIM++){       /* time animation loop */

      t=(double)TIM*dt;
      phi=phi+omega*dt;           /* delete this line for static plot */
 

for(n=1;n<Nt;n++){
  for(m=0;m<Np;m++){
u=Np*(n-1)+m+1;
u1=Np*(n-1)+mod(m+1,Np)+1;
u2=Np*n+mod(m+1,Np)+1;
u3=Np*n+m+1;
 pl[u].vertex[0][0]=X[u];
 pl[u].vertex[0][1]=Y[u];
 pl[u].vertex[0][2]=Z[u];
 pl[u].vertex[1][0]=X[u1];
 pl[u].vertex[1][1]=Y[u1];
 pl[u].vertex[1][2]=Z[u1];
 pl[u].vertex[2][0]=X[u2];
 pl[u].vertex[2][1]=Y[u2];
 pl[u].vertex[2][2]=Z[u2];
 pl[u].vertex[3][0]=X[u3];
 pl[u].vertex[3][1]=Y[u3];
 pl[u].vertex[3][2]=Z[u3];
 /* comment out these two lines if you do not want diff. rotation rate */
 getplaquettecenter(&pl[u]);
 getplaquettelatitude(&pl[u],AXIS[0],AXIS[1],AXIS[2]);
  }}

/* create rotation matrix */
 
 a11=cos(psi)*cos(phi)-cos(theta)*sin(phi)*sin(psi);
 a12=cos(psi)*sin(phi)+cos(theta)*cos(phi)*sin(psi);
 a13=sin(psi)*sin(theta);
 a21=-sin(psi)*cos(phi)-cos(theta)*sin(phi)*cos(psi);
 a22=-sin(psi)*sin(phi)+cos(theta)*cos(phi)*cos(psi);
 a23=cos(psi)*sin(theta);
 a31=sin(theta)*sin(phi);
 a32=-sin(theta)*cos(phi);
 a33=cos(theta);
/* rotate the plaquette points */
 for(n=1;n<=N;n++){
   for(m=0;m<4;m++){
tmpx=a11*pl[n].vertex[m][0]+a12*pl[n].vertex[m][1]+a13*pl[n].vertex[m][2];
tmpy=a21*pl[n].vertex[m][0]+a22*pl[n].vertex[m][1]+a23*pl[n].vertex[m][2];
tmpz=a31*pl[n].vertex[m][0]+a32*pl[n].vertex[m][1]+a33*pl[n].vertex[m][2];
pl[n].vertex[m][0]=tmpx;
pl[n].vertex[m][1]=tmpy;
pl[n].vertex[m][2]=tmpz;
}}

 /* now compute positions of centers and normals */
 for(n=1;n<=N;n++){
 pl[n].vdist=dist(vantagex,vantagey,vantagez,pl[n].center[0],pl[n].center[1],pl[n].center[2]);
 /* the (1.0, 0.3, 0.2) is the position of a light source */
 /* illuminating the globe. You may not want it */
 getplaquettecenter(&pl[n]);
 getplaquetteattribs(&pl[n],1.0,0.3,0.2);


 }

 /* add a "sunspot" or two*/
 pl[511].clr=30000;
 pl[456].clr=30000;
 pl[512].clr=0;
 getplaquettevelocity(&pl[512], AXIS[0], AXIS[1], AXIS[2], omega);
 /* print out radial velocity of sunspot */
 /* printf("%f\t%f\n", t, pl[512].velocity[2]); */

 for(n=0;n<N;n++){

/* color the plaquettes */
     color=g2_ink(dev1, (double)pl[n].clr/(double)65535,(double)pl[n].clr/(double)65535,(double)pl[n].clr/(double)65535); 
     // g2_pen(dev1, color); 
 g2_pen(dev1, black); 
/* only draw plaquettes with normals pointing towards viewer */
 if(pl[n].norml[2]>=0.0){
     pts[0]=300.0+500.0*pl[n].vertex[0][0];
     pts[1]=300.0+500.0*pl[n].vertex[0][1];
     pts[2]=300.0+500.0*pl[n].vertex[1][0];
     pts[3]=300.0+500.0*pl[n].vertex[1][1];

     pts[4]=300.0+500.0*pl[n].vertex[2][0];
     pts[5]=300.0+500.0*pl[n].vertex[2][1];
     pts[6]=300.0+500.0*pl[n].vertex[3][0];
     pts[7]=300.0+500.0*pl[n].vertex[3][1];

     g2_polygon(dev1, 4,pts);
     // g2_filled_polygon(dev1, 4,pts);

     // pl_colorname("black");

   


     // pl_fline(pl[n].vertex[2][0],pl[n].vertex[2][1],pl[n].vertex[3][0],pl[n].vertex[3][1]);
     //pl_fline(pl[n].vertex[0][0],pl[n].vertex[0][1],pl[n].vertex[1][0],pl[n].vertex[1][1]);
     // pl_fline(pl[n].vertex[3][0],pl[n].vertex[3][1],pl[n].vertex[0][0],pl[n].vertex[0][1]);
     //pl_fline(pl[n].vertex[1][0],pl[n].vertex[1][1],pl[n].vertex[2][0],pl[n].vertex[2][1]);
 }
 }
 g2_flush(dev1);
      
      g2_pen(dev1, white);
      g2_filled_rectangle(dev1, 0,0,600,600);
  } /*end of TIM animation loop */


 g2_close(dev);
}


double dist(double x1, double y1, double z1, double x2, double y2, double z2){
  /* computes squared point to point distance */
return( (x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2) );
}

void getplaquetteattribs( struct plaquette *Plaq, double u, double v, double w){
double del1x,del1y,del1z,del2x,del2y,del2z,dotprod,nlength,viewlength;
del1x=Plaq->vertex[1][0]-Plaq->vertex[0][0];
del1y=Plaq->vertex[1][1]-Plaq->vertex[0][1];
del1z=Plaq->vertex[1][2]-Plaq->vertex[0][2];
del2x=Plaq->vertex[3][0]-Plaq->vertex[0][0];
del2y=Plaq->vertex[3][1]-Plaq->vertex[0][1];
del2z=Plaq->vertex[3][2]-Plaq->vertex[0][2];
Plaq->norml[0]=del1y*del2z-del2y*del1z;
Plaq->norml[1]=del1z*del2x-del2z*del1x;
Plaq->norml[2]=del1x*del2y-del2x*del1y;
nlength=(Plaq->norml[0])*(Plaq->norml[0])+(Plaq->norml[1])*(Plaq->norml[1])+(Plaq->norml[2])*(Plaq->norml[2]);
nlength=sqrt(nlength);
viewlength=sqrt(u*u+v*v+w*w);
if(nlength!=0.0)
Plaq->face=((Plaq->norml[0])*u+(Plaq->norml[1])*v+(Plaq->norml[2])*w)/(nlength*viewlength);
else
Plaq->face=0.0;
/* this coloration scheme is a "soft" shadowing suitable for illumination */
/* thats what (u,v,w) are, the location of a light source */
/* so we dot the normal into this vector to determine shading */
/* from some chosen position of a light source */
Plaq->clr=(int)(65535.0*(1.0+tanh(2.0*Plaq->face))/2.0); 
/* All white */
// Plaq->clr=65535;

}

void getplaquettevelocity( struct plaquette *Plaq, double AXIS_x, double AXIS_y, double AXIS_z, double OMEGA){
    Plaq->velocity[0]=OMEGA*(AXIS_y*Plaq->center[2]-AXIS_z*Plaq->center[1]);
    Plaq->velocity[1]=-OMEGA*(AXIS_x*Plaq->center[2]-AXIS_z*Plaq->center[0]);
    Plaq->velocity[2]=OMEGA*(AXIS_x*Plaq->center[1]-AXIS_y*Plaq->center[0]);
}


void getplaquettelatitude( struct plaquette *Plaq,double AXIS_x, double AXIS_y, double AXIS_z){
    double lat, len;
    
    len=(Plaq->center[0])*(Plaq->center[0])+(Plaq->center[1])*(Plaq->center[1])+(Plaq->center[2])*(Plaq->center[2]);
    lat=Plaq->center[0]*AXIS_x+Plaq->center[1]*AXIS_y+Plaq->center[2]*AXIS_z;
    lat=acos(lat/sqrt(len));
    Plaq->latitude=lat;

}


void getplaquettecenter( struct plaquette *Plaq){
Plaq->center[0]=0.25*(Plaq->vertex[0][0]+Plaq->vertex[1][0]+ Plaq->vertex[2][0]+Plaq->vertex[3][0]);
Plaq->center[1]=0.25*(Plaq->vertex[0][1]+Plaq->vertex[1][1]+ Plaq->vertex[2][1]+Plaq->vertex[3][1]);
Plaq->center[2]=0.25*(Plaq->vertex[0][2]+Plaq->vertex[1][2]+ Plaq->vertex[2][2]+Plaq->vertex[3][2]);
}




double function(double x, double y){
double f;
/* Sphere */
 f=1.0; /* a globe, Y_{lm}=Y_{00} */
/* you clould make an oblate sun with an admixture of Y_{lm}*/
 f=1.0+0.1*sin(y)*sin(y);
return(0.2*f);
}

 int mod(int A, int B){
  if(A>=B)
return(A-B);
  else
 return(A-B*(A/B));
  }