/* barrierwave.c; produces PNM graphics files for mpeg_encode       */
/* to create mpeg movie of a wave colliding with a square barrier   */
/* requires libplot, compile with                                   */
/* gcc barrierwave2.c -L/usr/X11R6/lib -lplot -lXaw -lXmu -lXt -lSM */
/*  -lICE -lXext -lX11 -lm                                          */
/* run in an otherwise empty directory. it will create any number   */
/* of PNM graphics, labeled by consecutive integers                 */
/* J. R. Schmidt (1998)                                             */
#include<stdio.h>
#include<strings.h>
#include<stdlib.h>
#include<math.h>
#include<plot.h>

double ReR,ImR,ReT, ImT,ReA,ImA,ReB,ImB;
double  Den;
int m;
double k,lam,lam2,k2,eplus, eminus, s2lam,c2lam,fac;
double ratio, rat2, prod, diff, diff2, sum, sum2, Den;
double dx, x,dk,k,L,V,dt,t;
double Rsum, Imsum, arg1,arg2,Rright,Imright,Rleft,Imleft,Repsi,Impsi,psi2;
double lambda, elamx, arg3, arg4;
int n,m,l,time;
int pl_handle;
double firstX, firstY, secondX, secondY;




FILE *fptr;
char *name;
int dec,sign, num_char;
char *prefix;

main(){
 char *base_ptr=".pnm";
 char name[10], sbuff[8], *which;
L=2.0;
V=4.01;
k=2.8;
dx=0.01;
dt=0.01;
num_char=0; /* increase this if you expect very many frames */
 strcpy(name,base_ptr); /* place file extension in name variable */
/* create the trasmission/reflection coefficients */
if (k*k<2.0*V){
lam=sqrt(2.0*V-k*k);
k2=k*k;
lam2=2.0*V-k2;
eplus=exp(lam*L);
eminus=exp(-lam*L);
s2lam=sinh(2.0*L*lam);
c2lam=cosh(2.0*L*lam);
prod=2.0*k*lam;
diff2=k2-lam2;
sum2=2.0*V;
 Den=prod*prod+sum2*sum2*s2lam*s2lam;
 ReR=sum2*diff2*s2lam*s2lam/Den;
 ImR=-prod*sum2*s2lam*c2lam/Den;
 ReT=prod*prod*c2lam/Den;
 ImT=prod*diff2*s2lam/Den;
 ReB=eplus*k*(k*diff2*s2lam+prod*lam*c2lam)/Den;
 ImB=eplus*k*(lam*diff2*s2lam-prod*k*c2lam)/Den;
 ReA=-eminus*k*(k*diff2*s2lam-prod*lam*c2lam)/Den;
 ImA=eminus*k*(lam*diff2*s2lam+prod*k*c2lam)/Den;
}

else{
lam=sqrt(k*k-2.0*V);
ratio=lam/k;
sum2=ratio*ratio+1.0;
diff2=ratio*ratio-1;
rat2=ratio*ratio;
fac=2.0*lam*L;

Den=16.0*rat2+4.0*(1.0-rat2)*(1.0-rat2)*(sin(fac))*(sin(fac));
ReB=2.0*(ratio-1.0)*4.0*ratio*(cos(fac))/Den;
ImB=2.0*(ratio-1.0)*2.0*sum2*(sin(fac))/Den;
ReA=2.0*(ratio+1.0)*4.0*ratio*(cos(fac))/Den;
ImA=2.0*(ratio+1.0)*2.0*sum2*(sin(fac))/Den;
ReR=4.0*(1.0-rat2)*sum2*(sin(fac))*(sin(fac))/Den;
ImR=-8.0*ratio*(1-rat2)*(sin(fac))*(cos(fac))/Den;
ReT=4.0*ratio*(4.0*ratio*cos(fac))/Den;
ImT=4.0*ratio*(2.0*sum2*sin(fac))/Den;
}



 for(time=1;time<500;time++){
t=-6.0+(double)time*dt;
firstX=-10.0;
firstY=0.0;
which=fcvt(time,num_char,&dec,&sign); /* which is the value of time var */
strncat(which,name,4);
fptr=fopen(which,"w+");
/* set up a plotter */
pl_parampl("BITMAPSIZE","300x300");
pl_handle=pl_newpl("pnm",stdin,fptr,stderr);
pl_selectpl(pl_handle);
if (pl_openpl () < 0)       /* open Plotter */
         {
           fprintf (stderr, "Couldn't open Plotter\n");
           return 1;
         }
pl_fspace (-10.0, -1.2*V, 10.0, 1.2*V);
pl_flinewidth(0.1);
/* draw the barrier in red */

pl_pencolorname("red");
pl_fline(-L,0.0,-L,V);
pl_fline(-L,V,L,V);
pl_fline(L,V,L,0.0);
pl_pencolorname("black");
/* draw in the axes */
pl_fline(-10.0,0.0, 10.0, 0.0);
pl_fline(0.0, -1.2*V, 0.0, 1.2*V);
/* draw a label string */
pl_fmove(-9.0, -V);
pl_ffontsize (1.0);
pl_fontname ("HersheySerif");
pl_label("Square Barrier, k=2.8, V=4.01");




 for(l=-1000;l<=-200;l++){
x=(double)l*dx;
arg1=k*x-k*k*t/2.0;
arg2=k*x+2.0*k*L+k*k*t/2.0;
Rright=cos(arg1);
Imright=sin(arg1);
Rleft=cos(arg2);
Imleft=-sin(arg2);
Repsi=Rright+ReR*Rleft-ImR*Imleft;
Impsi=Imright+ReR*Imleft+ImR*Rleft;
secondX=x;secondY=Repsi;
pl_fline(firstX,firstY,secondX,secondY);
/*fprintf(fptr,"%f\t%f\n",x,Repsi);*/
firstX=secondX;firstY=secondY;
 }

for(l=-199;l<=199;l++){
x=(double)l*dx;

if(k*k<2.0*V){
 lambda=sqrt(2.0*V-k*k);
 elamx=exp(lambda*x);
 arg3=k*L+k*k*t/2.0;
 Rright=elamx*cos(arg3);
 Imright=-elamx*sin(arg3);
 Rleft=(cos(arg3))/elamx;
 Imleft=(-sin(arg3))/elamx;}
else{
 lambda=sqrt(k*k-2.0*V);
 arg3=lambda*x-lambda*L-k*L-t*k*k/2.0;
 arg4=-(lambda*x-lambda*L+k*L+t*k*k/2.0);
 Rright=cos(arg3);
 Imright=sin(arg3);
 Rleft=cos(arg4);
 Imleft=sin(arg4);
}
 Repsi=Rright*ReA-Imright*ImA+Rleft*ReB-Imleft*ImB;
 Impsi=Imright*ReA+Rright*ImA+Imleft*ReB+Rleft*ImB;

psi2=Rsum*Rsum+Imsum*Imsum;
secondX=x;secondY=Repsi;
pl_fline(firstX,firstY,secondX,secondY);
/*fprintf(fptr,"%f\t%f\n",x,Repsi);*/
firstX=secondX;firstY=secondY;
 }

for(l=200;l<=1000;l++){
x=(double)l*dx;

arg4=k*x-2.0*k*L-k*k*t/2.0;
Rright=cos(arg4);
Imright=sin(arg4);
Repsi=ReT*Rright-ImT*Imright;
Impsi=ReT*Imright+ImT*Rright;

psi2=Rsum*Rsum+Imsum*Imsum;
secondX=x;secondY=Repsi;
pl_fline(firstX,firstY,secondX,secondY);
/*fprintf(fptr,"%f\t%f\n",x,Repsi);*/
firstX=secondX;firstY=secondY;
 }
pl_closepl ();             
pl_selectpl (0);          
pl_deletepl (pl_handle);       
fclose(fptr);
 }     /* end of time loop */
return;
}