/* double-buffered packet animator  for g2-lib   */
#include<stdio.h>
#include<stdlib.h>
#include<complex.h>
#include<curses.h>
#include<math.h>
#include<g2.h>
#define LINES 35
#define COLS 100
#ifdef WIN32
#include<g2_win32.h>
#define WINTYPE 0
#endif
#ifdef XWIN
#include<g2_X11.h>
#define WINTYPE 1
#endif
#define PI 3.1415927
_Complex double iexpt[201][1000],iexpx[201][1000],psi,R[201], A[201],Denom;
int dev,dev1, white, black,red,time,lines, cols;
double V0, a, M, x, dx, y, r, par, xmax, first_x,test,k,dk,k0,sigma,P,C;
double t,firstX, firstY,secondX,secondY,Amp[201],res1, res2, res3;
int n,i,j, MAX,which,T,Ti,mode,type,p,pa,dispersion,q;
char label[50], input[10], output[10];
double function(double param, double x);
double diff_function(double param, double x);


main (int argc, char *argv[])
{

 WINDOW *pop_ptr;
    pop_ptr=initscr();
    echo();
    resize_term(LINES, COLS);
 
    if(!has_colors()){
        endwin();
        exit(0);
    }
    if(start_color() != OK){
        endwin();
        exit(0);
    }

init_pair(1, COLOR_RED, COLOR_BLACK);
    init_pair(2, COLOR_YELLOW, COLOR_BLACK);
    init_pair(3, COLOR_BLUE, COLOR_BLACK);
    init_pair(0, COLOR_WHITE, COLOR_BLACK);
    init_pair(6, COLOR_BLUE, COLOR_WHITE);
    init_pair(5, COLOR_CYAN, COLOR_BLACK);

   move(1, 10);
    attrset(COLOR_PAIR(6));
    printw("Hard-core quantum well bound states");
    attrset(COLOR_PAIR(0));
    f10:   move(3,1);
    printw("Packet central wavenumber k0 = ");
    attron(A_UNDERLINE);
    attrset(COLOR_PAIR(1));
    getstr(input);
    attrset(COLOR_PAIR(0));
    attroff(A_UNDERLINE);
    k0=(double)atof(input);
    if(k0<0.1 || k0>10.0){
             pop_ptr=newwin(4,25,10, 30);
             box(pop_ptr,'|','-');
             mvwprintw(pop_ptr,2,2,"0.1<=k0<=10.0");
             wrefresh(pop_ptr);  /* draw */
             napms(1500);
             wclear(pop_ptr);    /* clear */
             wrefresh(pop_ptr);  /* redraw */
             delwin(pop_ptr);    /* delete */
             refresh();
        move(3,1);
        for(i=0;i<strlen(input);i++)
        delch();
        goto f10;
    }

    f20:   move(5,1);
    printw("Packet sigma (x-space) = ");
    attron(A_UNDERLINE);
    attrset(COLOR_PAIR(1));
    getstr(input);
    attrset(COLOR_PAIR(0));
    attroff(A_UNDERLINE);
    sigma=(double)atof(input);
    if(sigma<0.5 || sigma > 4.0){
             pop_ptr=newwin(4,25,10, 30);
             box(pop_ptr,'|','-');
             mvwprintw(pop_ptr,2,2,"0.5<=sigma<=4.0");
             wrefresh(pop_ptr);  /* draw */
             napms(1500);
             wclear(pop_ptr);    /* clear */
             wrefresh(pop_ptr);  /* redraw */
             delwin(pop_ptr);    /* delete */
             refresh();
        move(5,1);
        for(i=0;i<strlen(input);i++)
        delch();
        goto f20;
    }


   f30:   move(7,1);
    printw("Well depth V0 = ");
    attron(A_UNDERLINE);
    attrset(COLOR_PAIR(1));
    getstr(input);
    attrset(COLOR_PAIR(0));
    attroff(A_UNDERLINE);
    V0=(double)atof(input);
    if(V0<0.0 || V0 > 20.0){
             pop_ptr=newwin(4,25,10, 30);
             box(pop_ptr,'|','-');
             mvwprintw(pop_ptr,2,2,"0.0<=V0<=20.0");
             wrefresh(pop_ptr);  /* draw */
             napms(1500);
             wclear(pop_ptr);    /* clear */
             wrefresh(pop_ptr);  /* redraw */
             delwin(pop_ptr);    /* delete */
             refresh();
        move(7,1);
        for(i=0;i<strlen(input);i++)
        delch();
        goto f30;
    }

 f33:   move(7,31);
    printw("Wall at a = ");
    attron(A_UNDERLINE);
    attrset(COLOR_PAIR(1));
    getstr(input);
    attrset(COLOR_PAIR(0));
    attroff(A_UNDERLINE);
    a=(double)atof(input);
    if(a<1.0 || a > 10.0){
             pop_ptr=newwin(4,25,10, 30);
             box(pop_ptr,'|','-');
             mvwprintw(pop_ptr,2,2,"1<=a<=10");
             wrefresh(pop_ptr);  /* draw */
             napms(1500);
             wclear(pop_ptr);    /* clear */
             wrefresh(pop_ptr);  /* redraw */
             delwin(pop_ptr);    /* delete */
             refresh();
        move(7,31);
        for(i=0;i<strlen(input);i++)
        delch();
        goto f33;
    }




   f40:   move(9,1);
    printw("type 0=|psi|^2 or 1=Re psi; ");
    attron(A_UNDERLINE);
    attrset(COLOR_PAIR(1));
    getstr(input);
    attrset(COLOR_PAIR(0));
    attroff(A_UNDERLINE);
    type=atoi(input);
    if(type<0 || type>1){
             pop_ptr=newwin(4,25,10, 30);
             box(pop_ptr,'|','-');
             mvwprintw(pop_ptr,2,2,"Not a valid choice");
             wrefresh(pop_ptr);  /* draw */
             napms(1500);
             wclear(pop_ptr);    /* clear */
             wrefresh(pop_ptr);  /* redraw */
             delwin(pop_ptr);    /* delete */
             refresh();
        move(9,1);
        for(i=0;i<strlen(input);i++)
        delch();
        goto f40;
    }


 f50:   move(11,1);
    printw("mode 0=single k=k0 or 1=packet; ");
    attron(A_UNDERLINE);
    attrset(COLOR_PAIR(1));
    getstr(input);
    attrset(COLOR_PAIR(0));
    attroff(A_UNDERLINE);
    mode=atoi(input);
    if(mode<0 || mode>1){
             pop_ptr=newwin(4,25,10, 30);
             box(pop_ptr,'|','-');
             mvwprintw(pop_ptr,2,2,"Not a valid choice");
             wrefresh(pop_ptr);  /* draw */
             napms(1500);
             wclear(pop_ptr);    /* clear */
             wrefresh(pop_ptr);  /* redraw */
             delwin(pop_ptr);    /* delete */
             refresh();
        move(11,1);
        for(i=0;i<strlen(input);i++)
        delch();
        goto f50;
    }

f60:   move(13,1);
    printw("Dispersion; non-relativistic=0, relativistic=1; ");
    attron(A_UNDERLINE);
    attrset(COLOR_PAIR(1));
    getstr(input);
    attrset(COLOR_PAIR(0));
    attroff(A_UNDERLINE);
    dispersion=atoi(input);
    if(dispersion<0 || dispersion>1){
             pop_ptr=newwin(4,25,10, 30);
             box(pop_ptr,'|','-');
             mvwprintw(pop_ptr,2,2,"Not a valid choice");
             wrefresh(pop_ptr);  /* draw */
             napms(1500);
             wclear(pop_ptr);    /* clear */
             wrefresh(pop_ptr);  /* redraw */
             delwin(pop_ptr);    /* delete */
             refresh();
        move(11,1);
        for(i=0;i<strlen(input);i++)
        delch();
        goto f60;
    }


   
  
   
    
    dk=3.0*(1.0/sigma)/100.0;
    time=1000;
    dx = 0.04;
    pa=(int)floor(a/dx);
    a=(double)pa*dx;
     M=2.0;  /* mass/hbar */
     C=3.0;

move(15,1);
printw("Using mass M = %f, C = %f, hbar=1.0, 200 Fourier modes with",M,C);
move(16,1);
printw("%2.4f<=k<=%2.4f,  dk=%1.5f", k0-100.0*dk, k0+100.0*dk, dk);
move(17,1);
printw("Well depth %2.4f between 0.0, a=%2.2f.", V0, a);
move(18,1);
 if(dispersion==0)
     printw("Non-relativistic waves, E=k^2/2M");
 else
     printw("Ultra-relativistic waves, E=kC, C=%2.2f", C); 

 /* figure out where low-lying resonances are */
 q=0;
 do{
     q=q+1;
 }while(((double)q+0.5)*PI/a<sqrt(2.0*M*V0));

 res1=sqrt( (((double)q+0.5)*PI/a)*(((double)q+0.5)*PI/a)-2.0*M*V0);
 res2=sqrt( (((double)(q+1)+0.5)*PI/a)*(((double)(q+1)+0.5)*PI/a)-2.0*M*V0);
 res3=sqrt( (((double)(q+2)+0.5)*PI/a)*(((double)(q+2)+0.5)*PI/a)-2.0*M*V0); 
move(19,1);
 printw("Resonances excited at k0 = %f, %f, %f", res1, res2, res3); 
refresh();


     for(n=0;n<=200;n++){
	  k=k0+(double)(n-100)*dk;
	  P=sqrt(k*k+2.0*M*V0);
	  Amp[n]=(1.414*sigma/sqrt(PI))*exp(-sigma*sigma*(k-k0)*(k-k0)/2.0);
	  Denom= sin(P*a)+I*(P/k)*cos(P*a);	  
	  R[n]=cexp(-2.0*I*k*a)*(sin(P*a)-I*(P/k)*cos(P*a))/Denom;
	  A[n]=2.0*cexp(-I*k*a)/Denom;	  	 
      }
      

      for(n=0;n<=200;n++){
	  k=k0+(double)(n-100)*dk;
	  for(Ti=0;Ti<time;Ti++){
	      t = (double) (Ti-400) *0.1/(k0);
	     if(dispersion==0){
		 iexpt[n][Ti]=cexp(-0.5*I*k*k*t/M);}
	     if(dispersion==1){
		 iexpt[n][Ti]=cexp(-I*fabs(k)*C*t);}
	  }
      }

      for(n=0;n<=200;n++){
	  k=k0+(double)(n-100)*dk; 
	  P=sqrt(k*k+2.0*M*V0);
	  
	  /* 0<=x<=20 */
	  for(p=0;p<=pa;p++){
	     x = (double) p *dx;
	     iexpx[n][p]=Amp[n]*A[n]*sin(P*x);
	  }
	 
	 
	  /* 0<=x<=20 */
	  for(p=pa;p<1000;p++){
	      x = (double) p *dx;
	      iexpx[n][p]=Amp[n]*(cexp(-I*k*x)+R[n]*cexp(I*k*x));
	  }


      }
  	      
      dev=g2_open_vd();
      if(WINTYPE==1)
	  dev1=g2_open_X11(600,400);
      if(WINTYPE==0)	 
	  dev1=g2_open_win32(600,400,"Quantum  states",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 (Ti = 0; Ti < time; Ti++)
	  {	
	      t = (double) T *0.04;    
	      g2_pen(dev1, red);
	      g2_line(dev1, 10, 200-20.0*V0, 10, 350);
	      g2_line(dev1, 10+15.0*a,200, 600, 200);
	      g2_line(dev1, 10,200-20.0*V0,10+15.0*a,200-20.0*V0); 
	      g2_line(dev1, 10+15.0*a,200-20.0*V0, 10+15.0*a,200);
	      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=iexpx[100][n]*iexpt[100][Ti];
	      
	  }
	  else{
	  for(p=-100;p<=100;p++){   
	      psi=psi+iexpx[p+100][n]*iexpt[p+100][Ti];
	  }

	  psi=psi*dk;
	  }

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

	  g2_line(dev1,10.0+15.0*firstX, 210+20.0*firstY, 10.0+15.0*secondX, 210+20.0*secondY);
	  firstX = secondX;
	  firstY = secondY;
		  }



	      g2_flush(dev1);
	       if(WINTYPE==0){sleep(1);}
	       if(WINTYPE==1){usleep(100);}

	      g2_pen(dev1, white);
	      g2_filled_rectangle(dev1, 0,0,600,600);   
	      
	      }  /* end of time loop */
move(22, 1);
attrset(COLOR_PAIR(0));
printw("<Enter> to quit ");
attrset(COLOR_PAIR(2));
getstr(input);
   endwin();
   exit(0);
    
}


double function(double param, double x){
    return(1.0/tan(x)+sqrt(param/(x*x)-1.0));
}

double diff_function(double param, double x){
    return(-1.0/(sin(x)*sin(x))-(param/(x*sqrt(param-x*x))));
}