/* step barier                                                      */
/* requires libplot, compile with                                   */
/* gcc barrier.c -L/usr/X11R6/lib -lplot -lXaw -lXmu -lXt -lSM  -lICE -lXext -lX11 -lm                                          */
/* creates a file called "packet.gif"    J.R Schmidt (1999)         */
#include<stdio.h>
#include<strings.h>
#include<stdlib.h>
#include<math.h>
#include<plot.h>

#define N 101

double ReR[N],ImR[N],ReT[N], ImT[N],ReA[N],ImA[N],ReB[N],ImB[N];
double Amp[N], Den[N];

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;


void build_coefficients(double *ReRptr, double *ImRptr,double *ReTptr, double *ImTptr,double *ReAptr, double *ImAptr,double *ReBptr, double *ImBptr, double *Ampptr, double Length, double potential, double delta_k);

FILE *fptr;
char *name;
int dec,sign, num_char;
char title[30];
char *clock;

main(){
float in2;
double *ReRptr,*ImRptr,*ReTptr,*ImTptr,*ReAptr,*ImAptr,*ReBptr,*ImBptr,*Ampptr;
L=1.0;
dk=0.01;
dx=0.01;
dt=0.15;
 
/* get data from user */
  printf ("\t\tQUANTUM PACKET STEP POTENTIAL SIMULATION\n\n");
  printf ("\tEnter a potential value V_0 between 0.0 and 10.0\n");
  scanf ("%f", &in2);
  V = (double) in2;
  printf ("\tYou chose a potential height of V=%f\n", in2);
  printf ("\tEnter a title string, 30 characters or less, NO WHITESPACE!\n");
  scanf ("%s", title);
  printf ("\tCreating animated GIF with 125 timesteps\n\n");

  for(n=0;n<101;n++){
ReR[n]=0.0;
ImR[n]=0.0;
ReT[n]=0.0;
ImT[n]=0.0;
ReA[n]=0.0;
ImA[n]=0.0;
ReB[n]=0.0;
ImB[n]=0.0;
Amp[n]=0.0;}
ReRptr=&ReR[0];
ImRptr=&ImR[0];
ReTptr=&ReT[0];
ImTptr=&ImT[0];
ReAptr=&ReA[0];
ImAptr=&ImA[0];
ReBptr=&ReB[0];
ImBptr=&ImB[0];
Ampptr=&Amp[0];
build_coefficients(ReRptr,ImRptr,ReTptr,ImTptr,ReAptr,ImAptr,ReBptr,ImBptr,Ampptr,L, V, dk);

 
  pl_parampl ("USE_DOUBLE_BUFFERING", "yes");
  pl_handle = pl_newpl ("X", stdin, fptr, stderr);
  pl_selectpl (pl_handle);
  if (pl_openpl () < 0)         /* open Plotter */
    {
      fprintf (stderr, "Couldn't open Plotter\n");
      return 1;
    }
  pl_fspace (-24.0, -4.0 * V, 24.0, 8.0 * V);
 

 for(time=1;time<125;time++){
pl_erase ();
/* 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 (-12.0, 0.0, 12.0, 0.0);
      pl_fline (0.0, -0.4 * V, 0.0, 3.0 * V);
/* draw a label string */
      pl_fmove (-9.0, -0.3*V);
      pl_ffontsize (1.0);
      pl_fontname ("HersheySerif");
      pl_label (title);

t=-6.0+(double)time*dt;
     firstX = -12.0;
      firstY = 0.0;
      /* get the time for display */
clock=fcvt((double)time,num_char,&dec,&sign);
    pl_fmove (9.0, -0.3*V);   
    pl_label (clock);



 for(l=-1200;l<=-100;l++){
x=(double)l*dx;
Rsum=0.0;
Imsum=0.0;

 for(m=0;m<N;m++){
k=(double)(m-50)*dk+2.5;
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[m]*Rleft-ImR[m]*Imleft;
Impsi=Imright+ReR[m]*Imleft+ImR[m]*Rleft;
Rsum=Rsum+dk*Repsi*Amp[m];
Imsum=Imsum+dk*Impsi*Amp[m];
}
 psi2=Rsum*Rsum+Imsum*Imsum;
          secondX = x;
          secondY = psi2;
          pl_fline (firstX, firstY, secondX, secondY);
          firstX = secondX;
          firstY = secondY;

 }

for(l=-99;l<=99;l++){
x=(double)l*dx;
Rsum=0.0;
Imsum=0.0;
for(m=0;m<N;m++){
k=(double)(m-50)*dk+2.5;
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[m]-Imright*ImA[m]+Rleft*ReB[m]-Imleft*ImB[m];
 Impsi=Imright*ReA[m]+Rright*ImA[m]+Imleft*ReB[m]+Rleft*ImB[m];
Rsum=Rsum+dk*Repsi*Amp[m];
Imsum=Imsum+dk*Impsi*Amp[m];
}
psi2=Rsum*Rsum+Imsum*Imsum;
          secondX = x;
          secondY = psi2;
          pl_fline (firstX, firstY, secondX, secondY);
          firstX = secondX;
          firstY = secondY;
 }

for(l=100;l<=1200;l++){
x=(double)l*dx;
Rsum=0.0;
Imsum=0.0;
for(m=0;m<N;m++){
k=(double)(m-50)*dk+2.5;
arg4=k*x-2.0*k*L-k*k*t/2.0;
Rright=cos(arg4);
Imright=sin(arg4);
Repsi=ReT[m]*Rright-ImT[m]*Imright;
Impsi=ReT[m]*Imright+ImT[m]*Rright;
Rsum=Rsum+dk*Repsi*Amp[m];
Imsum=Imsum+dk*Impsi*Amp[m];
}
psi2=Rsum*Rsum+Imsum*Imsum;
          secondX = x;
          secondY = psi2;
          pl_fline (firstX, firstY, secondX, secondY);
          firstX = secondX;
          firstY = secondY;
 }
 }     /* end of time loop */
 pl_closepl ();
  pl_selectpl (0);
  pl_deletepl (pl_handle);
 
  printf("\tDone. The file is called packet.gif. View it with Netscape\n\n");
return;
}


void build_coefficients(double *ReRptr, double *ImRptr,double *ReTptr, double *ImTptr,double *ReAptr, double *ImAptr,double *ReBptr, double *ImBptr,double *Ampptr, double Length, double potential, double delta_k)
  {
int m;
double k,lam,lam2,k2,eplus, eminus, s2lam,c2lam,fac;
double ratio, rat2, prod, diff, diff2, sum, sum2, Den;
double normalization;

normalization=pow(N/3.1415962,0.25);

 for(m=0;m<N;m++){
k=(double)(m-50)*delta_k+2.5;
if (k*k<2.0*potential){
lam=sqrt(2.0*potential-k*k);
k2=k*k;
lam2=2.0*potential-k2;
eplus=exp(lam*Length);
eminus=exp(-lam*Length);
s2lam=sinh(2.0*Length*lam);
c2lam=cosh(2.0*Length*lam);
prod=2.0*k*lam;
diff2=k2-lam2;
sum2=2.0*potential;
 Den=prod*prod+sum2*sum2*s2lam*s2lam;
 *(ReRptr+m)=sum2*diff2*s2lam*s2lam/Den;
 *(ImRptr+m)=-prod*sum2*s2lam*c2lam/Den;
 *(ReTptr+m)=prod*prod*c2lam/Den;
 *(ImTptr+m)=prod*diff2*s2lam/Den;
 *(ReBptr+m)=eplus*k*(k*diff2*s2lam+prod*lam*c2lam)/Den;
 *(ImBptr+m)=eplus*k*(lam*diff2*s2lam-prod*k*c2lam)/Den;
 *(ReAptr+m)=-eminus*k*(k*diff2*s2lam-prod*lam*c2lam)/Den;
 *(ImAptr+m)=eminus*k*(lam*diff2*s2lam+prod*k*c2lam)/Den;
}

else{
lam=sqrt(k*k-2.0*potential);
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));
*(ReBptr+m)=2.0*(ratio-1.0)*4.0*ratio*(cos(fac))/Den;
*(ImBptr+m)=2.0*(ratio-1.0)*2.0*sum2*(sin(fac))/Den;
*(ReAptr+m)=2.0*(ratio+1.0)*4.0*ratio*(cos(fac))/Den;
*(ImAptr+m)=2.0*(ratio+1.0)*2.0*sum2*(sin(fac))/Den;
*(ReRptr+m)=4.0*(1.0-rat2)*sum2*(sin(fac))*(sin(fac))/Den;
*(ImRptr+m)=-8.0*ratio*(1-rat2)*(sin(fac))*(cos(fac))/Den;
*(ReTptr+m)=4.0*ratio*(4.0*ratio*cos(fac))/Den;
*(ImTptr+m)=4.0*ratio*(2.0*sum2*sin(fac))/Den;
}

*(Ampptr+m)=normalization*exp(-4.0*(k-2.5)*(k-2.5));
 }
  }