/* Draw  plane polarized light figures for illustration    */
/* gcc lightwave3.c -L/usr/X11R6/lib -lplot -lXaw -lXmu -lXt -lSM  -lICE -lXext -lX11 -lm                                                */
/* Jeffrey R. Schmidt */
/* X output  */
#include<stdio.h>
#include<strings.h>
#include<stdlib.h>
#include<math.h>
#include<plot.h>
#define N 200                              
#define PI 3.1415927

double EX[N],EY[N],EZ[N], EX0[N],EY0[N], EZ0[N]; /* base and tips of arrows */
double BX[N],BY[N],BZ[N],BX0[N],BY0[N],BZ0[N]; /* base and tips of arrows */
double dx,omega,E,size,theta,phi,psi,tmpx, tmpy, tmpz,one[3],X[3],Y[3],Z[3];
double a11,a12,a13,a21,a22,a23,a31,a32,a33,dt,t;
int n,m,k,time;
int pl_handle;
char buffer[6],*labels;
FILE *fptr;
char *name;
int dec,sign, num_char;
char *prefix;
void draw_arrow(double x0, double y0, double x, double y);
void draw_arrow2(double x0, double y0, double x, double y);
main (int argc, char *argv[])
{
 char *base_ptr=".pnm";
 char name[10], sbuff[8], *which; 
    srand(19);
    size=6.0;
    dx=size/100.0;
    omega=0.7*PI;
    E=size/4.0;
num_char=0; /* increase this if you expect very many frames */
 strcpy(name,base_ptr); /* place file extension in name variable */ 
    theta=0.8;
    phi=0.3;
    psi=0.1;
    dt=0.07123;



    /* pl_parampl ("USE_DOUBLE_BUFFERING", "yes"); */
    pl_parampl ("BITMAPSIZE", "400x400");
       pl_parampl ("BG_COLOR", "white");
       pl_parampl ("TRANSPARENT_COLOR", "white");
       pl_parampl ("GIF_ITERATIONS", "100"); 
       pl_parampl ("GIF_DELAY", "5"); 
  pl_handle = pl_newpl ("GIF", stdin, stdout, stderr);
  pl_selectpl (pl_handle);
  if (pl_openpl () < 0)         /* open Plotter */
    {
      fprintf (stderr, "Couldn't open Plotter\n");
      return 1;
    }
  pl_fspace (-size, -size, size, size);
  pl_flinewidth (0.005);
 for(time=1;time<100;time++){
t=(double)time*dt;
  pl_erase(); 

 for(n=0;n<N;n++){
     EX0[n]=-0.8*size+(double)n*dx;
     BX0[n]=-0.8*size+(double)n*dx;
     EY0[n]=0.0;
     EZ0[n]=0.0;
     BZ0[n]=0.0;
     BY0[n]=0.0;
     EX[n]=-0.8*size+(double)n*dx;
     EY[n]=E*cos(omega*(double)n*dx-t);
     EZ[n]=-E*sin(omega*(double)n*dx-t);
     BX[n]=-0.8*size+(double)n*dx;
     BZ[n]=-E*sin(omega*(double)n*dx-t);
     BY[n]=-E*cos(omega*(double)n*dx-t);
 }
 one[0]=0.0;
 one[1]=0.0;
 one[2]=0.0;
 X[0]=0.6*size;
 X[1]=0.0;
 X[2]=0.0;
 Y[1]=0.5*size;
 Y[0]=0.0;
 Y[2]=0.0;
 Z[2]=0.5*size;
 Z[1]=0.0;
 Z[0]=0.0;

 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);

tmpx=a11*X[0]+a12*X[1]+a13*X[2];
tmpy=a21*X[0]+a22*X[1]+a23*X[2];
tmpz=a31*X[0]+a32*X[1]+a33*X[2];
X[0]=tmpx;
X[1]=tmpy;
X[2]=tmpz;
tmpx=a11*Y[0]+a12*Y[1]+a13*Y[2];
tmpy=a21*Y[0]+a22*Y[1]+a23*Y[2];
tmpz=a31*Y[0]+a32*Y[1]+a33*Y[2];
Y[0]=tmpx;
Y[1]=tmpy;
Y[2]=tmpz;
tmpx=a11*Z[0]+a12*Z[1]+a13*Z[2];
tmpy=a21*Z[0]+a22*Z[1]+a23*Z[2];
tmpz=a31*Z[0]+a32*Z[1]+a33*Z[2];
Z[0]=tmpx;
Z[1]=tmpy;
Z[2]=tmpz;

/* now rotate the arrows using phi and theta */
  for(n=0;n<N;n++){
tmpx=a11*EX0[n]+a12*EY0[n]+a13*EZ0[n];
tmpy=a21*EX0[n]+a22*EY0[n]+a23*EZ0[n];
tmpz=a31*EX0[n]+a32*EY0[n]+a33*EZ0[n];
EX0[n]=tmpx;
EY0[n]=tmpy;
EZ0[n]=tmpz;
tmpx=a11*EX[n]+a12*EY[n]+a13*EZ[n];
tmpy=a21*EX[n]+a22*EY[n]+a23*EZ[n];
tmpz=a31*EX[n]+a32*EY[n]+a33*EZ[n];
EX[n]=tmpx;
EY[n]=tmpy;
EZ[n]=tmpz;
tmpx=a11*BX0[n]+a12*BY0[n]+a13*BZ0[n];
tmpy=a21*BX0[n]+a22*BY0[n]+a23*BZ0[n];
tmpz=a31*BX0[n]+a32*BY0[n]+a33*BZ0[n];
BX0[n]=tmpx;
BY0[n]=tmpy;
BZ0[n]=tmpz;
tmpx=a11*BX[n]+a12*BY[n]+a13*BZ[n];
tmpy=a21*BX[n]+a22*BY[n]+a23*BZ[n];
tmpz=a31*BX[n]+a32*BY[n]+a33*BZ[n];
BX[n]=tmpx;
BY[n]=tmpy;
BZ[n]=tmpz;
  }

pl_pencolorname ("black");
for(n=0;n<=120;n++){
    draw_arrow(EX0[n],EY0[n],EX[n],EY[n]);
}
pl_ffontname("HersheySerif");
 draw_arrow2(one[0],one[1],X[0],X[1]);
 draw_arrow2(one[0],one[1],Y[0],Y[1]);
 draw_arrow2(one[0],one[1],Z[0],Z[1]);
 pl_fmove(1.1*X[0],1.1*X[1]);
 pl_label("x");
 pl_fmove(1.1*Y[0],1.1*Y[1]);
 pl_label("y");
 pl_fmove(1.1*Z[0],1.1*Z[1]);
 pl_label("z");

pl_pencolorname ("red");
for(n=0;n<=120;n++){
    draw_arrow(BX0[n],BY0[n],BX[n],BY[n]);
}
 }

  pl_closepl ();
  pl_selectpl (0);
  pl_deletepl (pl_handle);

 
}



void draw_arrow(double x0, double y0, double x, double y){
    double tip1x, tip1y, tip2x, tip2y;
    tip1x=x0+0.85*(x-x0)+0.075*(y-y0);
    tip1y=y0+0.85*(y-y0)-0.075*(x-x0);
    tip2x=x0+0.85*(x-x0)-0.075*(y-y0);
    tip2y=y0+0.85*(y-y0)+0.075*(x-x0);
    pl_fline(x0, y0 ,x, y);
    pl_fline(x, y ,tip1x, tip1y);
    pl_fline(x, y ,tip2x, tip2y);
 
}
void draw_arrow2(double x0, double y0, double x, double y){
    double tip1x, tip1y, tip2x, tip2y;
    tip1x=x0+0.9*(x-x0)+0.06*(y-y0);
    tip1y=y0+0.9*(y-y0)-0.06*(x-x0);
    tip2x=x0+0.9*(x-x0)-0.06*(y-y0);
    tip2y=y0+0.9*(y-y0)+0.06*(x-x0);
    pl_fline(x0, y0 ,x, y);
    pl_fline(x, y ,tip1x, tip1y);
    pl_fline(x, y ,tip2x, tip2y);
 
}