/* double-buffered X-windows/WIN32 packet animator       */
/* compile with;                                         */
/* gcc g2_packet.c -L/usr/X11R6/lib -lX11 -lm -lgd -lg2  */
/* or use "makemake.sh" to generate a makefile           */
/* requires g2                                           */
/* 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], potX[200], potY[200];
_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 < 5 || argc > 5)
    {
      printf ("\n\tThe program needs five command line parameters\n");
      printf ("\t\tusage; ./packet k0 sigma 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 2.0 3.0  1 1\n");
      exit (1);
    }
  
  k0 = (double) atof (argv[1]);
  sigma = (double) atof (argv[2]);
  type = (int) atoi (argv[3]);
  mode = (int) atoi (argv[4]);

  

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


  if (type < 0 || type > 1)
    {
      printf ("Invalid wavefunction display type, try again\n");
      exit (1);
    }
 
  if (mode < 0 || mode > 1)
    {
      printf ("Invalid wavefunction mode, 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.12/k0;
	     iexpt[n][T]=cexp(-0.5*I*k*k*t/M);
	  }
      }

      for(n=0;n<=200;n++){
	  k=k0+(double)(n-100)*dk;
	  
	  /* 0<=x<=20 */
	  for(p=0;p<1000;p++){
	     x = (double) p *dx;
	     iexpx[n][p]=(tanh(x-L)-I*k)*cexp(I*k*x)/(1.0+I*k);
	     
	  }
	  
      }

      for(n=0;n<200;n++){
	  potX[n]=(double)n*5.0*dx;
	  potY[n]=-4.0/(cosh(potX[n]-L)*cosh(potX[n]-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 (reflectionless potential)",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.12/k0;
      /* if we draw one mode, do it slowly */
      if(mode==0){
	  if(WINTYPE==0) sleep(1);
	  else
	      usleep(5);
      }

      
/* 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);
      for(n=0;n<199;n++)
	  g2_line(dev1, 15.0*potX[n], 200.0+20.0*potY[n], 15.0*potX[n+1], 200.0+20.0*potY[n+1]); 
      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);
      if(WINTYPE==0) sleep(1);
      g2_pen(dev1, white);
      g2_filled_rectangle(dev1, 0,0,600,400);
      
    }				/* end of time loop */
 


}