/* double-buffered X-windows packet animator                 */
/* compile with;                                             */
/* gcc g2_packet.c -L/usr/X11R6/lib -lX11 -lm -lgd -lg2  */
/* requires GNU plotutils                                    */
/* creates X display in real time or GIF89a graphics file    */
/* J.R.Schmidt (2000/2006)                                        */
#include<stdio.h>
#include<stdlib.h>
#include<complex.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 PI 3.1415927
double V0, dx, x, k,dk,k0, sigma, t,firstX, firstY, secondX, secondY,norm,sigma,M,factor,L;
double Amp[201];
_Complex double phase, psi, exponential, kp, iexpt[201][1000],iexpx[201][1000];

int n, m, p, time, T, N, type, displ,mode;
int dev,dev1, white, black,red;

double modulus(__complex__ double X);
main (int argc, char *argv[])
{
  if (argc < 6 || argc > 6)
    {
      printf ("\n\tThe program needs five command line parameters\n");
      printf ("\t\tusage; ./packet k0 sigma V0 type mode\n\n");     
      printf ("\t\tk0 is a double precision packet central wavevector\n");
      printf ("\t\tsigma is the width in Fourier space, double precision (>1) \n");
      printf ("\t\ttype is an integer, 0 to show psi*psi, 1 to show Re psi\n");      printf ("\t\t mode =0 for single k-vector, mode=1 for a packet. \n");
      printf ("\t\texample; ./packet 6.0 2.0 4.0 1 1\n");
      exit (1);
    }
  
  k0 = (double) atof (argv[1]);
  sigma = (double) atof (argv[2]);
  V0=(double)atof(argv[3]);
  type = (int) atoi (argv[4]);
  mode = (int) atoi (argv[5]);

  

  if (sigma <= 0.85 || sigma > 3.0)
    {
      printf ("Bad sigma, try again\n");
      exit (1);
    }
 if (V0 <= 0.0 || V0 > 20.0)
    {
      printf ("Bad V0, try again\n");
      exit (1);
    }

  if (type < 0 || type > 1)
    {
      printf ("Invalid wavefunction display type, try again\n");
      exit (1);
    }
 
 

    /* most calculations can be taken out of the time/position loops */

  /* nymber of time steps */
    time = 1000;
    /* step starts at x=L */
    L=20.0;
 
    dx = 0.04;
    M=2.0;  /* mass/hbar */
      
      dk=2.0*(1.0/sigma)/100.0;
      factor=1.414*sigma/sqrt(PI);
      /* sigma only appears squared after this point */
      sigma=sigma*sigma;
  
      for(n=0;n<=200;n++){
	  k=(double)(n-100)*dk;
	  Amp[n]=exp(-sigma*k*k/2.0);
      }

      for(n=0;n<=200;n++){
	  k=k0+(double)(n-100)*dk;
	  for(T=0;T<time;T++){
	     t = (double) T *0.02;
	     iexpt[n][T]=cexp(-0.5*I*k*k*t/M);
	  }
      }

      for(n=0;n<=200;n++){
	  k=k0+(double)(n-100)*dk;
	  if(0.5*k*k/M>=V0)
	      kp=sqrt(k*k-2.0*M*V0);
	  else
	      kp=I*sqrt(2.0*M*V0-k*k);
	  /* 0<=x<=20 */
	  for(p=0;p<=500;p++){
	     x = (double) p *dx;
	     iexpx[n][p]=cexp(I*k*(x))+((k-kp)/(k+kp))*cexp(-I*k*x+2.0*I*k*L);
	  }
	  /* 0<=x<=20 */
	  for(p=500;p<1000;p++){
	     x = (double) p *dx;
	     iexpx[n][p]=(2.0*k/(k+kp))*cexp(I*kp*(x)+I*k*L-I*kp*L);
	  }


      }
	      
      dev=g2_open_vd();
      if(WINTYPE==1)
	  dev1=g2_open_X11(600,600);
      if(WINTYPE==0)	 
	  dev1=g2_open_win32(600,400,"Quantum Packet Simulator",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);
  for (T = 0; T < time; T++)
    {
	
      t = (double) T *0.02;
      /* if we draw one mode, do it slowly */
      if(mode==0)sleep(10);

      
/* draw the well in red */
      
      
      g2_pen(dev1, red);
      g2_line(dev1, 0,200, 600, 200);
      g2_line(dev1,1,200-250,1, 200+250);
      g2_line(dev1,300, 200, 300, 200+20*0.5*V0); 
      g2_line(dev1,300, 200+20*0.5*V0,600, 200+20*0.5*V0); 
      g2_pen(dev1, black);


      firstX = 0.0;
      firstY = 0.0;
      

      for (n = 0; n < 1000; n++)
	{
	    x = (double) n *dx; 
	    /* assemble the packet */
	  psi=0.0;
	  if(mode==0){
	      psi=Amp[100]*iexpx[100][n]*iexpt[100][T];
	      
	  }
	  else{
	  for(p=-100;p<=100;p++){   
	      psi=psi+factor*Amp[p+100]*iexpx[p+100][n]*iexpt[p+100][T];
	  }

	  psi=psi*dk;
	  }

	  secondX = x;
	  if (type == 0)
	      secondY = csqrt(psi*conj(psi));
	  else
	      secondY = creal(psi);

	  g2_line(dev1,15*firstX, 200+20*firstY, 15*secondX, 200+20*secondY);
	  firstX = secondX;
	  firstY = secondY;
	}
      g2_flush(dev1);
      sleep(1);
      g2_pen(dev1, white);
      g2_filled_rectangle(dev1, 0,0,600,400);
      
    }				/* end of time loop */
 


}