/* double-buffered X-windows packet animator                 */
/* compile with;                                             */
/* gcc newwavepacket.c -L/usr/X11R6/lib -lplot -lXaw -lXmu -lXt -lSM -lICE -lXext -lX11 -lm                               */
/* requires GNU plotutils                                    */
/* creates X display in real time or GIF89a graphics file    */
/* J.R.Schmidt (2006)                                        */
/* uses the C99 complex number specification */
#include<stdio.h>
#include<stdlib.h>
#include<features.h>
#include<complex.h>
#include<math.h>
#include<plot.h>
#define PI 3.1415927
double dx, x, k0, sigma, t,firstX, firstY, secondX, secondY,norm,alpha, ratio;

_Complex double phase, psi, factor, factor2;

int n, m, time, T, N, type, displ;
int pl_handle;
FILE *fptr;
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 type time display\n\n");
      
      printf ("\t\tk0 is a double precision packet central wavevector\n");
      printf ("\t\tsigma is the width in Fourier space, double precision\n");
      printf ("\t\ttype is an integer, 0 to show psi*psi, 1 to show Re psi\n");
      printf ("\t\ttime is an integer number of time frames\n");
      printf ("\t\tdisplay=0 for X (real time), 1 for GIF89a\n");
      printf ("\t\texample; ./packet 2.0 0.5 1 100 0\n");
      exit (1);
    }
  dx = 0.04;
  ratio=0.1;  /* hbar over m */
  k0 = (double) atof (argv[1]);
  alpha = (double) atof (argv[2]);
  type = (int) atoi (argv[3]);
  time = (int) atoi (argv[4]);
  displ = (int) atoi (argv[5]);

  if (alpha <= 0.0)
    {
      printf ("Bad parameters, try again\n");
      exit (1);
    }
  if (type < 0 || type > 1)
    {
      printf ("Invalid wavefunction display type, try again\n");
      exit (1);
    }
  if (displ < 0 || displ > 1)
    {
      printf ("Invalid plotter display type, try again\n");
      exit (1);
    }


  if (displ == 0)
    {
      pl_parampl ("BITMAPSIZE", "400x250");
      pl_parampl ("BG_COLOR", "white");
      pl_parampl ("USE_DOUBLE_BUFFERING", "yes");
      pl_handle = pl_newpl ("X", stdin, fptr, stderr);
    }
  else
    {
      fptr = fopen ("packet.gif", "w+");
      pl_parampl ("BITMAPSIZE", "250x150");
      pl_parampl ("TRANSPARENT_COLOR", "white");
      pl_parampl ("GIF_ITERATIONS", "100");
      pl_parampl ("GIF_DELAY", "3");
      pl_handle = pl_newpl ("gif", stdin, fptr, stderr);
    }

  pl_selectpl (pl_handle);
  if (pl_openpl () < 0)		/* open Plotter */
    {
      fprintf (stderr, "Couldn't open Plotter\n");
      return 1;
    }

  if (type == 0)
    pl_fspace (-20.0, -10.0, 20.0, 10.0);
  else
    pl_fspace (-20.0, -10.0, 20.0, 10.0);
  pl_flinewidth (0.01);

  for (T = 0; T <= time; T++)
    {
      t = (double) T *0.1;

      pl_erase ();
/* draw the well in red */

      pl_pencolorname ("red");
      if (type == 0)
	{
	  pl_fline (-20.0, 0.0, 20.0, 0.0);
	  pl_fline (0.0, -10.0, 0.0, 10.0);
	}
      else
	{
	  pl_fline (-20.0, 0.0, 20.0, 0.0);
	  pl_fline (0.0, -10.0, 0.0, 10.0);
	}
      pl_fontname ("HersheySerif");
      pl_ffontsize (2.0);
      pl_pencolorname ("black");
      if (type == 0)
	{
	  pl_fmove (-8.9, 7.3);
	  pl_label ("\\*Q\\sp*\\ep(x,t)\\*Q(x,t)");
	}
      else
	{
	  pl_fmove (-8.9, 5.3);
	  pl_label ("Re\\*Q(x,t)");
	}


      firstX = -20.0;
      firstY = 0.0;
      norm=sqrt(sqrt(2.0*alpha/PI));
      
      factor=alpha+0.5*t*ratio*I;
      
      for (n = 0; n <= 1000; n++)
	{
	   
	  x = -20.0+(double) n *dx;
	  
	  factor2=(2.0*alpha*k0+(x+8.0)*I);
	  factor2=factor2*factor2/(4.0*factor)-alpha*k0*k0;
	  psi=norm*csqrt(PI/factor)*cexp(factor2);
	  secondX = x;
	  if (type == 0)
	      secondY = csqrt(psi*conj(psi));
	  else
	      secondY = creal(psi);

	  pl_fline (firstX, firstY, secondX, secondY);
	  firstX = secondX;
	  firstY = secondY;
	}
    }				/* end of time loop */
  pl_closepl ();
  pl_selectpl (0);
  pl_deletepl (pl_handle);
  if (displ == 1)
    close (fptr);

}