/* 
 *  Solves the Schroedinger equation for a wavepacket inside
 *  a square-well.
 *  Observe how the dymanical development of the wavefunction 
 *  displays the several square-well energy components of the wave.
 *  Sqwell.html: <applet code= Sqwell.class width=600 height=450></applet>
 *  by C.A. Bertulani, 08/15/2000
 */
import java.awt.*;
import java.applet.Applet;

public class Sqwell extends Applet implements Runnable {
  Thread runstring = null;   
  private Image offScreenImage;  // used for no flckering
  int ixmesh0, ixmeshN; 
  double y0, yN;
  double Stretchxx,Coefxx,Stretchyy,Coefyy,pi,k0,dx;
  int AppletHeight, AppletWidth; 
  double psr[][],psi[][],p2[];
  int jj, max;
  // Labels in applet
  String graphLabel1 = "Wavepacket in a square-well";
  String graphLabel2 = "(wait and see all kinds of wave phenomena)";
  // Font for caption
  Font smallArialFont = new Font("Arial", Font.PLAIN, 15);
  
  public void start()
  {
    if(runstring == null);
    {
     runstring = new Thread(this);
     runstring.start();
    }
  }

  public void stop()
  {
     if(runstring != null){
	   runstring.stop();	 
       runstring = null;
     }
  }

  public void run()
  {
    while(true){
      repaint();
      try{ Thread.sleep(10);
      }
      catch(InterruptedException e){ }
    }
  } 

  public void update(Graphics g)  //To avoid flicker
  {
    Graphics offScreenGraphics=offScreenImage.getGraphics();
	 //Erase the previous image
	Dimension d = size();
    offScreenGraphics.setColor(getBackground());
    offScreenGraphics.fillRect(0,0,AppletWidth,AppletHeight);
    offScreenGraphics.setColor(g.getColor());
	 //Paint the frame into the image
    paint(offScreenGraphics);
    g.drawImage(offScreenImage,0,0,this);

  }

  private int pixx(int ivalue)  { 
    // Given the  index of the mesh-xpoint return the pixel equivalent 
    int jvalue = (int)(Stretchxx*ivalue + Coefxx);
    return jvalue; 
  }  

  private int pixy(double yvalue)  { 
    // Given the y-value in desired scale, return the pixel equivalent 
	int jvalue = (int)(Stretchyy*yvalue + Coefyy);
    return jvalue; 
  }  

  void CoordUndo(double ixmesh0,double ixmeshN,double y0,double yN)
  {  // changes x-mesh-indexing and y-coordinate definition into
     // java applet coordinate system of pixels	
    double W = AppletWidth;
    double H = AppletHeight;
	// put graph within a margin from applet window
	double margin = 0.2;  // percent used in margin
    double xxmax = (1.0-margin)*W; // right x-pixel-limit  
    double xxmin = margin*W; // left x-pixel-limit
    double yymin = (1.0-margin)*H; // upper y-pixel-limit
    double yymax = margin*H; // lower y-pixel-limit
	// fit data within graph axis 
	// first transform mesh indices into x-window-pixel coordinates
    Stretchxx = (xxmax-xxmin)/(ixmeshN-ixmesh0);
    Coefxx = (xxmin*ixmeshN-xxmax*ixmesh0)/(ixmeshN-ixmesh0);
	// now transform graph y-values to y-window-pixel coordinates
    Stretchyy = (yymax-yymin)/(yN-y0);
    Coefyy = (yymin*yN-yymax*y0)/(yN-y0);
  }

  void axisg(int ix1,int ix2,double y1,double y2,Graphics g)
  {   // Plot the axis with desired cordinates
    int xpix1,xpix2,ypix1,ypix2;
    xpix1 = pixx(ix1); 
	xpix2 = pixx(ix2);
	ypix1 = pixy(y1);
	ypix2 = pixy(y2);
    g.drawLine(xpix1,ypix1,xpix2,ypix1);
    g.drawLine(xpix1,ypix1,xpix1,ypix2);
    g.drawLine(xpix2,ypix1,xpix2,ypix2);
 }

 public void Equation(Graphics g, int jj)
 {
   int pixx1,pixx2,pixy1,pixy2;
   double dt,x,fact2;
   int i,n;
   dt=dx*dx;   // a convenient choice for time step dt
   if(jj==1){
     // initial conditions
     x=0.0;
     for(i=0;i<max+1;i++){
       p2[i]=0.0;
     }
     for(i=0;i<max+1;i++){ // a gaussian wavepacket
       //real part of the initial wave function at t=0
       fact2=2.0*(x-5.0);
       psr[i][0] = (Math.cos(k0*x))*Math.exp(-0.5*fact2*fact2);
       //imaginary part of the initial wave function at t=0
       psi[i][0] = (Math.sin(k0*x))*Math.exp(-0.5*fact2*fact2);
       x = x+dx;
     }
   }//if ij
   //propagate solution through time. 
   for(n=0; n<jj+1; n++){
     // real part of the wave function and probability
     for( i=1;i<max;i++){ // Euler's method. You can improve that.
       psr[i][1]=psr[i][0] - dt*(psi[i+1][0] + psi[i-1][0]
                 -2.0*psi[i][0])/(2.0*dx*dx);
       p2[i]=psr[i][0]*psr[i][1]+psi[i][0]*psi[i][0];
     }
     //imaginary part of the wave function
     for(i=1;i<max;i++){
       psi[i][1] = psi[i][0]+dt*(psr[i+1][1] + psr[i-1][1]
                     -2.0*psr[i][1])/(2.0*dx*dx);
     }
     //only at certain time instants we want to know
     // the probability
     if(n==jj){
       for(i=1;i<max+1;i++){
         if(i>1){
            // transform coordinates to pixel coordinates			 
			pixx1=pixx(i-1);
			pixx2=pixx(i);
			pixy1=pixy(p2[i-1]);
			pixy2=pixy(p2[i]);
			g.setColor(Color.red);
			// draw probabilities
            g.drawLine(pixx1,pixy1,pixx2,pixy2);
          }
       }
     }
     // new iterations are now the old ones
     for(i=1;i<max+1;i++){
       psi[i][0]=psi[i][1];
       psr[i][0]=psr[i][1];
     }
  }
 }

  public void init()
 {
//   double ixmesh0, ixmeshN, y0, yN;   // plot limits in desired units
   psr = new double[751][2];  // real part of wavefunction
   psi = new double[751][2];  // imag part of wavefunction
   // create an image with the same size of the applet area
   // How big is the applet? 
   Dimension d = size(); 
   AppletHeight = d.height; 
   AppletWidth  = d.width; 
   offScreenImage = createImage(AppletWidth,AppletHeight);
   p2 = new double[751]; // probability
   jj = 1; // sets initial condition
   pi = 3.14159265358979323846;
   max = 750; // index of last point in the x mesh
   k0 = 8.0*pi;  // particle momentum
   dx = 0.02;  // mesh step size
   // limits on x-mesh and on desired y-values
   ixmesh0 = 0; 
   ixmeshN = max+1; 
   y0 = 0.0; 
   yN = 1.0; 
   CoordUndo(ixmesh0,ixmeshN,y0,yN); // transform pixels to coordinates
 }

 public void paint(Graphics g){
    axisg(ixmesh0,ixmeshN,y0,yN,g); // plot axis
    // find position for caption in pixels			 
    int pixxL =pixx(ixmesh0);
	int pixyL = AppletHeight - AppletHeight/10;
	g.setColor(Color.blue);
	g.setFont(smallArialFont);
    g.drawString(graphLabel1,pixxL, pixyL);  // draw caption
    g.drawString(graphLabel2,pixxL, pixyL+pixyL/20);  
    Equation(g,jj);  // iterate and plot wave
    jj = 30;
 }
}