/* gcc surfaceplotter.c -L/usr/X11R6/lib -lplot -lXaw -lXmu -lXt -lSM     */
/*  -lICE -lXext -lX11 -lm                                                */
#include<stdio.h>
#include<strings.h>
#include<stdlib.h>
#include<math.h>
#include<plot.h>

#define N 31 /* should be odd */
/* creates a N X N grid, with middle being x-y axes */
/* rotates figure by two Euler angles for different */
/* camera angles                                    */ 
#define PI 3.1415927

double Z[N*N],X[N*N],Y[N*N];
double theta,phi,tmpx,tmpy,tmpz,dx,dy,dt,x,y,dimensionx, dimensiony,t;
double axes[4][2];
double function(double x, double y, double t);
 int mod(int A, int B);

int n,m,time,Tmax;
int pl_handle;
FILE *fptr;
main (int argc, char *argv[])
{

  dx = 0.2;
  dy = 0.2;
  dt = 0.05;
  dimensionx=(double)N*dx;
  dimensiony=(double)N*dy;
  theta=(double)atof(argv[1]);
  phi=(double)atof(argv[2]);
  Tmax=atoi(argv[3]);
  fptr=fopen("3d.gif","w+");
 


  pl_parampl ("BG_COLOR", "white");
  pl_parampl ("BITMAPSIZE", "250x250");
  pl_parampl ("TRANSPARENT_COLOR", "white");
  pl_parampl ("GIF_ITERATIONS", "100");
  pl_parampl ("GIF_DELAY", "3");
  pl_handle = pl_newpl ("GIF", stdin, fptr, stderr);
  pl_selectpl (pl_handle);
  if (pl_openpl () < 0)         /* open Plotter */
    {
      fprintf (stderr, "Couldn't open Plotter\n");
      return 1;
    }
  pl_fspace (-0.71*dimensionx, -0.71*dimensiony,0.71*dimensionx, 0.71*dimensiony);
  pl_flinewidth (0.01);

  for(time=1;time<Tmax;time++){
t=(double)time*dt;

 axes[0][0]=-(double)(N+1)*dx/2; /* points to draw x,y axes */
  axes[0][1]=0.0;
  axes[1][0]=(double)(N+1)*dx/2;
  axes[1][1]=0.0;
  axes[2][1]=-(double)(N+1)*dy/2;
  axes[2][0]=0.0;
  axes[3][1]=(double)(N+1)*dy/2;
  axes[3][0]=0.0;

  for(n=0;n<N*N;n++){
    X[n]=(double)(mod(n,N))*dx;
    Y[n]=(double)(mod((n-X[n])/N,N))*dy;
    /* now apply the offset */
    X[n]=X[n]-(double)(N-1)*dx/2.0;
    Y[n]=Y[n]-(double)(N-1)*dy/2.0;
    Z[n]=function(X[n],Y[n],t);
}

/* now rotate the plot using phi and theta */
  for(n=0;n<N*N;n++){
tmpx=cos(phi)*X[n]+sin(phi)*Y[n];
tmpy=-cos(theta)*sin(phi)*X[n]+cos(theta)*cos(phi)*Y[n]+sin(theta)*Z[n];
X[n]=tmpx;
 Y[n]=tmpy;
 /*printf("%f\t%f\n",X[n],Y[n]);*/ }
  /* rotate axis lines too */
  for(n=0;n<4;n++){
tmpx=cos(phi)*axes[n][0]+sin(phi)*axes[n][1];
tmpy=-cos(theta)*sin(phi)*axes[n][0]+cos(theta)*cos(phi)*axes[n][1];
axes[n][0]=tmpx;
 axes[n][1]=tmpy;}


pl_erase();
 pl_pencolorname ("black");
  /* now build up the grid */
  for(n=0;n<N*(N-1);n++){
    if(mod(n,N)!=N-1){
pl_fline(X[n],Y[n],X[n+1],Y[n+1]);
 pl_fline(X[n],Y[n],X[n+N],Y[n+N]);}
    } 
  for(n=N*(N-1);n<N*N-1;n++){
    pl_fline(X[n],Y[n],X[n+1],Y[n+1]); }

  for(n=1;n<N-1;n++){
    pl_fline(X[n*N],Y[n*N],X[(n+1)*N],Y[(n+1)*N]); }

  for(n=0;n<N-1;n++){
    pl_fline(X[n*N+N-1],Y[n*N+N-1],X[(n+1)*N+N-1],Y[(n+1)*N+N-1]); }
  /* draw the axes in red */
 pl_pencolorname ("red");
pl_fline(axes[0][0],axes[0][1],axes[1][0],axes[1][1]);
pl_fline(axes[2][0],axes[2][1],axes[3][0],axes[3][1]);
  /* put in some labels */
pl_fmove(axes[1][0],axes[1][1]);
pl_ffontname("HersheySerif");
pl_ffontsize(0.35);
pl_label("x");
pl_fmove(axes[3][0],axes[3][1]);
pl_label("y"); 
pl_fmove(-0.69*dimensionx, -0.50*dimensiony);
pl_ffontsize(0.45);
 pl_label("\\*Q(x,t)\\*Q\\sp*\\ep(x,t), \\*w\\sbx\\eb=2\\*p, \\*w\\sby\\eb=\\*p  ");        
  } /* end of time loop */         
  pl_closepl ();
  pl_selectpl (0);
  pl_deletepl (pl_handle);
fclose(fptr);

}

double function(double x, double y, double t){
double phi1,phi2,omega1,omega2,r,arg1,arg2;
phi1=0.0;
phi2=1.5;
omega1=6.28;
omega2=3.14;
r=0.5;
arg1=x-1.414*r*cos(omega1*t-phi1);
arg1=arg1*arg1;
arg2=y-1.414*r*cos(omega2*t-phi2);
arg2=arg2*arg2;
return (2.0*exp(-arg1-arg2));
}


 int mod(int A, int B){
  if(A>=0)
return(A-B*(A/B));
  else{
       do{A=A+B;}while(A<0);
       return(A);}
  }