/*
 * qm5.java - revised 3 Mar 08 - width 385, height 281
 * @author jackord@kw.igs.net
 * quantum harmonic oscillator
 *   Pulse A from the quantum string placed between the ground state
 *     turning points of the quantum oscillator
 *   direct solution of the time-dependent Schrodinger equation or
 *     expansion in terms of stationary states (a 54-term series)
 * classical mass/spring system motion given by Feynman algorithm
 */

import java.applet.Applet;
import java.awt.*;
import java.awt.event.ActionListener;
import java.awt.event.ActionEvent;

public class qm5 extends Applet
                 implements ActionListener {
  int kk=0; int first=0; int n=24; int n2=576;
  int nq=576; int np=384; int nc=288; int ns=54;
  double enf=0; double ens=0;
  String b1s="Schr Eqn"; Button b1=new Button(b1s);
  String b2s="Series"; Button b2=new Button(b2s);
  Image bim;
  Graphics bgr;
  double [][] y=new double[ns+1][np+1];                // Global arrays
  double [] ar=new double[ns+1];
  double [] ai=new double[ns+1];
  int [] xt=new int[121];
  int [] yt=new int[121];

  public void init() {
    setBackground(new Color(211, 211, 211));
    add(b1); add(b2);
    b1.addActionListener(this);
    b2.addActionListener(this);
  }

  public void paint(Graphics g) {
    int nf, c, d, yy, yen, xtp, sf, sy, xc, y1, y2, x2;
    long tt=System.currentTimeMillis(); int del=16;
    double [] yr=new double[nq+1];
    double [] yrr=new double[nq+1];
    double [] dyr=new double[nq+1];
    double [] yi=new double[nq+1];
    double [] yii=new double[nq+1];
    double [] dyi=new double[nq+1];
    double dt, du, pi, phi, s1, s2, re, im, w, energy;
    double xs, vx;
    if (first==0) {                                    // Initialize...
      bim=createImage(385, 281);                       // ...Animation buffer
      bgr=bim.getGraphics();
      b1.setBounds(2, 2, 65, 20);
      b2.setBounds(2, 27, 65, 20);
      first=1;
      bgr.setFont(new Font("TimesRoman", Font.PLAIN, 14));
      states();                                        // ...y[][]
      for (int i=0; i<=60; i=i+1) {                    // ...Potential function
        yt[60+i]=240-4*i; yt[60-i]=yt[60+i];
        xt[60+i]=np/2+(int)(n*Math.sqrt(4.*i/15.)+.5);
        xt[60-i]=np-xt[60+i];
      }
    }
    nf=320; c=160; sf=4000; sy=500; xc=192;         
    d=(nq-np)/2;
    yen=0; xtp=0; 
    pi=Math.PI; dt=n2*pi/nf/c; du=pi/nf/c/n2; w=2*pi/nf; 
    g.setColor(Color.black);
    g.drawRect(0, 0, 384, 280);
    if (kk>0) {
      for (int i=nc-n; i<=nc+n; i=i+1) {               // Init pulse A...
        phi=(i-(nc-n))*pi/n; yr[i]=(1-Math.cos(phi))/Math.sqrt(3*n); 
        yi[i]=-yr[i]*Math.sin(phi);
        yr[i]=yr[i]*Math.cos(phi);
      }
      for (int i=nc-n; i<=nc+n; i=i+1) {               // ...and dy's
        dyr[i]=dt*(yi[i-1]+yi[i+1]-2*yi[i])-du*(i-nc)*(i-nc)*yi[i];
        dyi[i]=-dt*(yr[i-1]+yr[i+1]-2*yr[i])+du*(i-nc)*(i-nc)*yr[i];
      }
      enf=0;                                           // ...enf
      for (int i=nc-n; i<=nc+n; i=i+1) {               
        s1=yr[i]*(-n2*(yr[i-1]+yr[i+1]-2*yr[i])+yr[i]*(i-nc)*(i-nc)/n2);
        s2=yi[i]*(-n2*(yi[i-1]+yi[i+1]-2*yi[i])+yi[i]*(i-nc)*(i-nc)/n2);
        enf=enf+s1+s2;
      }
      ens=0; 
      for (int k=1; k<=ns; k=k+2) {
        ar[k]=0; 
        for (int i=nc-n; i<=nc+n; i=i+1) {             // ...ar[]
          ar[k]=ar[k]+yr[i]*y[k][i-d];
        }
        ens=ens+ar[k]*ar[k]*(2*k-1);
      }
      for (int k=2; k<=ns; k=k+2) {                    // ...ai[]
        ai[k]=0;
        for (int i=nc-n; i<=nc+n; i=i+1) {
          ai[k]=ai[k]+yi[i]*y[k][i-d];
        }
        ens=ens+ai[k]*ai[k]*(2*k-1);                   // ...ens
      }
      if (kk==2) { energy=ens; } else { energy=enf; } 
      energy=(int)(1000*energy+.5)/1000.;
      xtp=(int)(24*Math.sqrt(energy)); xs=xc; vx=xtp;  // Initialize
      yen=240-(int)(15*energy+.5);                     // Plot parameters
      for (int jt=0; jt<=nf; jt=jt+1) {                // Main loop
        if (jt==0) {
          for (int jj=1; jj<=2; jj=jj+1) {             // yr and yi
            bgr.clearRect(0, 0, 384, 280);           
            bgr.setColor(Color.black);
            bgr.drawRect(0, 0, 384, 280); bgr.drawLine(0, 240, 384,240);
            if (kk==2) {                               // Series yr and yi
              for (int i=0; i<=np; i=i+1) {
                y[0][i]=0;
                for (int k=1; k<=ns; k=k+1) {
                  if (jj==1) { s1=ar[k]; } else { s1=ai[k]; }
                  y[0][i]=y[0][i]+s1*y[k][i];
                }
              }
              bgr.setColor(Color.magenta);             // Components
              for (int k=1; k<=ns; k=k+1) {
                s1=0;
                for (int i=1; i<=np; i=i+1) {
                  if (jj==1) { s2=ar[k]*y[k][i]; } else { s2=ai[k]*y[k][i]; }
                  bgr.drawLine(i-1, 100-(int)(sy*s1), i, 100-(int)(sy*s2)); s1=s2;
                }
              }
            }
            else {                                     // Initial yr and yi
              for (int i=0; i<=np; i=i+1) {
                if (jj==1) { s1=yr[i+d]; } else { s1= yi[i+d]; }
                y[0][i]=s1;
              }
            }
            bgr.setColor(Color.blue);
            for (int i=1; i<=np; i=i+1) {
              bgr.drawLine(i-1, 100-(int)(sy*y[0][i-1]), i, 100-(int)(sy*y[0][i]));
            }
            if (jj==1) {
              bgr.drawString("REAL COMPONENT", 130, 20);
            }
            else {
              bgr.drawString("IMAGINARY COMPONENT", 120, 20);
            }
            g.drawImage(bim, 0, 0, null);              // Show plot buffer
            tt=2000+System.currentTimeMillis();
            while (tt>System.currentTimeMillis()) { }
          }
        }
        bgr.clearRect(0, 0, 384, 280);
        bgr.setColor(Color.black);
        bgr.drawRect(0, 0, 384, 280); bgr.drawLine(0, 240, 384, 240);
        bgr.drawString("Frame "+jt+"/"+nf, 5, 70);
        bgr.drawString("E/gs = "+energy, 285, yen+6);
        bgr.setColor(Color.blue);
        for (int i=1; i<=np-1; i=i+1) {                
          if (kk==2) {                                 // Series 
            re=0; im=0;
            for (int k=1; k<=ns; k=k+1) {
              re=re+y[k][i]*(ar[k]*Math.cos((k-1)*jt*w)-ai[k]*Math.sin((k-1)*jt*w));
              im=im+y[k][i]*(ar[k]*Math.sin((k-1)*jt*w)+ai[k]*Math.cos((k-1)*jt*w));    //no - need 4 terms
            }
            yy=240-(int)(sf*(re*re+im*im)+.5);
          }
          else {                                        // Schrodinger equation
            yy=240-(int)(sf*(yr[i+d]*yr[i+d]+yi[i+d]*yi[i+d])+.5);
          }
          if (yy<240) { bgr.drawLine(i, yy, i, 239); }
        }
        x2=(int)(xs)-10; bgr.fillRect(x2, 250, 20, 20);  // Mass...
        bgr.setColor(Color.black); y1=260;
        for (int i=1; i<=x2; i=i+1) {                  // ...on a spring
          y2=260-(int)(10*Math.sin(16*pi*i/x2));
          bgr.drawLine(i-1, y1, i, y2); y1=y2;
        }
        bgr.setColor(Color.red);
        bgr.drawPolyline(xt, yt, 121);                 // Draw potential
        bgr.drawLine(168, 225, 216, 225);              // Ground state energy
        bgr.drawLine(216, 225, 216, 239);              //   ...and turning
        bgr.drawLine(168, 225, 168, 239);              //     ...points
        bgr.setColor(Color.cyan);
        bgr.drawLine(192-xtp, yen, 192+xtp, yen);      // Draw energy
        bgr.drawLine(192+xtp, yen, 192+xtp, 239);      // Draw turning...
        bgr.drawLine(192-xtp, yen, 192-xtp, 239);      //   ...points
        g.drawImage(bim, 0, 0, null);                  // Show plot buffer
        if (jt==0) {
          tt=2000+System.currentTimeMillis();
          while (tt>System.currentTimeMillis()) { }
        }
        xs=xs+vx*w; vx=vx+(xc-xs)*w;                   // Update mass on spring
        if (kk==1) {
          for (int jj=1;jj<=c; jj=jj+1) {              // Accuracy loop
            for (int i=1; i<=nq-1; i=i+1) {            // Project ahead dy/2
              yrr[i]=yr[i]+dyr[i]/2; yii[i]=yi[i]+dyi[i]/2;
            }
            for (int i=1; i<=nq-1; i=i+1) {             
              dyr[i]=dt*(yii[i-1]+yii[i+1]-2*yii[i])-du*(i-nc)*(i-nc)*yii[i];
              yr[i]=yr[i]+dyr[i];
              dyi[i]=-dt*(yrr[i-1]+yrr[i+1]-2*yrr[i])+du*(i-nc)*(i-nc)*yrr[i];
              yi[i]=yi[i]+dyi[i];
            }
          }
        }
        while (tt>System.currentTimeMillis()) { }
        tt=tt+del;
      }
      kk=0; 
    }
  } 

  public void states() {                               // Find y[][]
    double e, z, v, a, s, x, dx, tp;
    int ss, jj, tst, ii, nn;
    ss=1; dx=0.1; nn=192;
    for (int k=1; k<=ns; k=k+1) {
      e=2*k-1; z=1+ss; y[k][nn]=z; v=1-ss;              
      s=z*z/2; a=-e*z/n2; v=v+a*dx/2;                  // Feynman half-step
      x=0; jj=0; ii=0; tst=0; tp=n*Math.sqrt(e);
      do {                                             // Feynman loop
        jj=jj+1; x=x+dx; z=z+v*dx; a=(x*x/n2-e)*z/n2; v=v+a*dx;
        if ((x>tp)&&(v*z>0)) { z=0; tst=1; }
        s=s+z*z;
        if (jj==10) {
          jj=0; ii=ii+1; if (ii<=nn) { y[k][nn+ii]=z; }
        }
      } while (tst==0);
      s=Math.sqrt(2*s*dx);          
      for (int i=0; i<=nn; i=i+1) {                    // Normalize
        y[k][nn+i]=y[k][nn+i]/s; y[k][nn-i]=ss*y[k][nn+i];
      }
      ss=-ss;
    }
  }

  public void actionPerformed(ActionEvent e) {         // Buttons
    String tst;
    tst=e.getActionCommand();
    if (b1s.equals(tst)) { kk=1; }
    if (b2s.equals(tst)) { kk=2; }
    repaint();
  }

}