/* double-buffered packet animator  for g2-lib   */
#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 k1, k2, k3, C,M, t,firstX, firstY, secondX, secondY, x, dx;
double Amp[201];
_Complex double psi;
int n, m, p, time, T, N, type, mode;
int dev,dev1, white, black,red;

main (int argc, char *argv[])
{
  if (argc < 2 || argc > 2)
    {
      printf ("\n\tThe program needs a parameter\n");
      printf ("\t\tusage; ./a.out mode\n\n");           
      printf ("\t\tmode=22, 23, 24, 25\n");
      exit (1);
    }

  mode = (int) atoi (argv[1]);
  if (mode < 22 || mode > 25)
    {
      printf ("try again\n");
      exit (1);
    }

    time = 1000; /* number of time steps */
    k1=2.0*PI;
    k2=2.0*PI+0.5;
    k3=2.0*PI-0.5;
    C=1.0;
    dx = 0.04;
    M=1.0;  /* mass/hbar */
      	      
      dev=g2_open_vd();
      if(WINTYPE==1)
	  dev1=g2_open_X11(600,400);
      if(WINTYPE==0)	 
	  dev1=g2_open_win32(600,400,"Quantum 441 Assignment 4",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;    
	      g2_pen(dev1, red);
	      g2_line(dev1, 0,200, 600, 200);
	      g2_line(dev1,1,200-250,1, 200+250); 
	      g2_pen(dev1, black);
	      firstX = 0.0;
	      firstY = 0.0;
	      for (n = 0; n < 1000; n++)
		  {
		      x = (double) n *dx; 
		      /* assemble the packet */
		      if(mode==22)
			  psi=0.707*(cexp(I*k1*x-I*k1*C*t)+cexp(I*k2*x-I*k2*C*t));
		      if(mode==23)
			  psi=0.707*(cexp(I*k1*x-0.5*I*k1*k1*t/M)+cexp(I*k2*x-0.5*I*k2*k2*t/M)); 
		      if(mode==24)
			  psi=(0.5*cexp(I*k3*x-0.5*I*k3*k3*t/M)+0.5*cexp(I*k2*x-0.5*I*k2*k2*t/M)+0.707*cexp(I*k1*x-0.5*I*k1*k1*t/M)); 

		      if(mode==25)
			  psi=(0.5*cexp(I*k3*x-I*k3*C*t)+0.5*cexp(I*k2*x-I*k2*C*t)+0.707*cexp(I*k1*x-I*k1*C*t)); 
		      secondX = x;
		      secondY = creal(psi);
		      g2_line(dev1,15*firstX, 200+20*firstY, 15*secondX, 200+20*secondY);
		      firstX = secondX;
		      firstY = secondY;
		  }
	      if(mode==22)g2_string(dev1,200, 100,"non-dispersive sum of 2 waves");
	      if(mode==23)g2_string(dev1,200, 100,"dispersive sum of 2 waves");
	      if(mode==24)g2_string(dev1,200, 100,"dispersive sum of 3 waves");
	      if(mode==25)g2_string(dev1,200, 100,"non-dispersive sum of 3 waves");

	      g2_flush(dev1);
	      if(WINTYPE==0){sleep(1);}
	      if(WINTYPE==1){usleep(5000);}
	      g2_pen(dev1, white);
	      g2_filled_rectangle(dev1, 0,0,600,400);      
	  }  /* end of time loop */
}