/*
 * nm7.java - revised 22 Nov 07 - width 481, height 221
 * @author jackord@kw.igs.net
 * Pulse B incident on a variable thickness segment where the mass density
 *   increases by a factor 1.44 followed by a region where the mass density
 *   increases again by a factor of 1.44
 * the motion is generated either directly from Newton's Law
 *   or from a superposition of (non-orthogonal) normal mode displacements
 */

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

public class nm7 extends Applet
                 implements ActionListener, Runnable {
  int kk=0; int plt=0; int first=0; int dl=0;          // Declarations
  int n=480; int nm=1; int ns=120;
  double pi=Math.PI;
  double [] k=new double[ns+1];
  double [][] q=new double[ns+1][n+1];
  String b1s="Pulse B"; Button b1=new Button(b1s);
  String b2s="Motion"; Button b2=new Button(b2s);
  Choice ch=new Choice();
  Checkbox ch1=new Checkbox("Series");
  Image bim;
  Graphics bgr;
  Thread tmr;

  public void init() {
    setLayout(new FlowLayout(FlowLayout.LEFT));
    setBackground(new Color(211, 211, 211));
    ch1.setBackground(getBackground());
    ch.setBackground(getBackground());
    ch.addItem("D(film)=0"); ch.addItem("D(film)=5"); ch.addItem("D(film)=10");
    add(b1); b1.addActionListener(this);
    add(b2); b2.addActionListener(this);
    add(ch); add(ch1);
  }

  public void start() {
    if (tmr==null) {
      tmr=new Thread(this);
      tmr.setPriority(Thread.MIN_PRIORITY);
      tmr.start();
    }
  }

  public void stop() {
    tmr=null; plt=0;
  }

  public void run() {
    repaint();
  }

  double f(double k1) {                                // The upper right...
    double c1, c, c2, s1, s, s2, k, k2, a, b;          
    a=pi/n;
    k=1.2*k1; k2=1.2*k;
    c1=Math.cos(a*k1*n/2); s1=Math.sin(a*k1*n/2);
    c=Math.cos(a*k*dl); s=Math.sin(a*k*dl);
    c2=Math.cos(a*k2*(n/2-dl)); s2=Math.sin(a*k2*(n/2-dl));
    b=k*k*s*s1*s2-k1*c1*k*c*s2-k1*c1*k2*c2*s-k*c*k2*c2*s1;
    return b;                                          // ...matrix element
  }

  public void findk() {                                // Find the k-values
    int cnt=0;                                         //   for which the
    double a, b, c, e, x;                              //   matrix element = 0
    a=f(.5); e=0.00001; x=.6;
    do {                                               // Scan for zero crossings
      b=f(x);                                          //   in the range of interest
      if (a*b<0) { cnt=cnt+1; k[cnt]=x; }
      a=b; x=x+.05;
    } while (cnt<ns);
    for (int i=nm; i<=ns; i=i+1) {                     // Use a bisection algorithm
      a=k[i]-.1; b=k[i];                               //   to refine the values
      while ((b-a)>e) {
        c=(a+b)/2; if (f(c)*f(a)>0) { a=c; } else { b=c; }
      }
      k[i]=(a+b)/2;
    }
  }

  public void fey() {                                  // Find normal modes
    double z, v, a, c, e, s, dx;                       //   from nm to ns
    int ii, jj;
    c=pi/n; c=c*c; dx=0.05;
    for (int i=nm; i<=ns; i=i+1) {                      
      e=k[i]*k[i]; a=-e*c; s=0;
      z=0; v=1; ii=0; jj=0;
      for (int j=1; j<=9600; j=j+1) {
        z=z+v*dx;
        if (j==4800) { a=1.44*a; }
        if (j==(4800+20*dl)) { a=1.44*a; }
        v=v+a*z*dx; s=s+z*z; jj=jj+1;
        if (jj==20) {
          jj=0; ii=ii+1; q[i][ii]=z;
        }
      }
      s=Math.sqrt(s*dx);
      for (int j=1; j<=480; j=j+1) {                   // Normalize them
        q[i][j]=q[i][j]/s;
      }
    }
  }

  public void paint(Graphics g) {
    int n=480;                                         // Declarations
    int [] yy=new int[n+1];
    double [] y=new double[n+1];
    double [] v=new double[n+1];
    double [] b=new double[ns+1];
    double [] d=new double[ns+1];
    double [] bb=new double[ns+1];
    double [] dd=new double[ns+1];
    double [] w=new double[ns+1];
    double [][] o=new double[ns+1][ns+1];
    double m, e1, e2, s1, s2, z;
    double dt=0.05;
    int no, nw, nf, jt, del, delt;
    no=48; nw=96; nf=320; dl=0; e1=0; e2=0; del=24;    // Constants
    dl=5*ch.getSelectedIndex();
    g.setColor(Color.black);
    g.drawRect(0, 0, 480, 220);
    if (first==0) {
      bim=createImage(481, 221);                       // Animation buffer
      bgr=bim.getGraphics();
      first=1;
    }
    if (kk>0) {
      for (int i=no; i<=no+nw; i=i+1) {                // Init Pulse
        z=(i-no)*2*pi/nw;
        y[i]=40*(1-Math.cos(z));
        v[i]=-40*2*pi/nw*Math.sin(z);
        v[i]=v[i]*Math.cos(4*z)+y[i]*4*2*pi/nw*Math.sin(4*z);
        y[i]=y[i]*Math.cos(4*z);
        e1=e1+v[i]*v[i];                               // Incident energy (KE)
      }

      if (ch1.getState()==true) {                      // Series
        findk();                                       // Find k[i]
        fey();                                         // Find y[i][j]
        for (int i=nm; i<=ns; i=i+1) {                  
          w[i]=k[i]*pi/n;                              // No dispersion
          for (int j=1; j<=n-1; j=j+1) {               // Expansion coefficients
            bb[i]=bb[i]+y[j]*q[i][j];                  //   (first approximation)
            dd[i]=dd[i]+v[j]/w[i]*q[i][j];
          }
        }
        for (int i=nm; i<=ns-1; i=i+1) {               // Mode overlap integrals
          for (int j=i+1; j<=ns; j=j+1) {
            s1=0;
            for (int jj=1; jj<=n-1; jj=jj+1) {
              s1=s1+q[i][jj]*q[j][jj];
            }
            o[i][j]=s1; o[j][i]=s1;
          }
        }
        for (int i=nm; i<=ns; i=i+1) {                 // Expansion coefficients
          s1=0; s2=0;                                  //   (second approximation)
          for (int j=nm; j<=ns; j=j+1) {
            s1=s1+bb[j]*o[i][j];
            s2=s2+dd[j]*o[i][j];
          }
          b[i]=bb[i]-s1;
          d[i]=dd[i]-s2;
        }
        s1=0;
        for (int i=nm; i<=ns; i=i+1) {
          s1=s1+b[i]*q[i][no+nw/2];
        }
        for (int i=nm; i<=ns; i=i+1) {                 // Expansion coefficients
          b[i]=b[i]*80/s1; d[i]=d[i]*80/s1;            //   (final correction)
        }
        for (int j=1; j<=n-1; j=j+1) {
          y[j]=0;
          for (int i=nm; i<=ns; i=i+1) {
            y[j]=y[j]+b[i]*q[i][j];
          }
        }
      }
      else {
        for (int i=no; i<=no+nw; i=i+1) {              // Feynman half-step
          v[i]=v[i]+(y[i-1]+y[i+1]-2*y[i])*dt/2;
        }
      }
      jt=0; yy[0]=120; yy[n]=yy[0];
      long tt=System.currentTimeMillis();
      do {
        bgr.setColor(getBackground());                 // Clear plot buffer
        bgr.fillRect(0, 0, 480, 220);
        bgr.setColor(Color.black); bgr.drawRect(0, 0, 480, 220);
        bgr.drawString("Frame "+jt+"/"+nf, 380, 20);
        bgr.setColor(Color.blue);
        for (int i=1; i<=n; i=i+1) {                   // Draw blue string
          yy[i]=120-(int)(y[i]+.5);
          if (i==n/2) { bgr.setColor(Color.green); }
          if (i==n/2+dl) { bgr.setColor(Color.red); }
          bgr.drawLine(i-1, yy[i-1], i, yy[i]);
        }
        g.drawImage(bim, 0, 0, null);                  // Show plot buffer
        jt=jt+1;                                       // Step time
        if (ch1.getState()==true) {                    // Series
          for (int i=1; i<=n-1; i=i+1) {
            y[i]=0;
            for (int k=nm; k<=ns; k=k+1) {
              y[i]=y[i]+q[k][i]*(b[k]*Math.cos(w[k]*jt)+d[k]*Math.sin(w[k]*jt));
            }
          }
        }
        else {
          for (int j=1; j<=20; j=j+1) {
            for (int i=1; i<=n-1; i=i+1) {               // Step y[]
              y[i]=y[i]+v[i]*dt;
            }
            m=1;
            for (int i=1; i<=n-1; i=i+1) {               // Step v[] (F=ma)
              if (i==n/2) { m=1.44*m; }
              if (i==n/2+dl) { m=1.44*m; }
              v[i]=v[i]+(y[i-1]+y[i+1]-2*y[i])/m*dt;
            }
          }
        }
        if (tt>System.currentTimeMillis()) { delt=del; } else { delt=0; }
        try {
          Thread.sleep(delt);
        }
        catch (InterruptedException e) {
          stop();
        }
        tt=tt+del;
        if (jt==nf+1) {
          for (int i=1; i<=n/2-1; i=i+1) {
            if (ch1.getState()==true) {
              v[i]=0;
              for (int k=nm; k<=ns; k=k+1) {
                v[i]=v[i]+w[k]*q[k][i]*(d[k]*Math.cos(w[k]*nf)-b[k]*Math.sin(w[k]*nf));
              }
            }
            e2=e2+v[i]*v[i];                           // Reflected energy (KE)
          }
          e2=(int)(10000*e2/e1)/100.;
          g.drawString("Reflects "+e2+" % of incident energy", 150, 45);
          kk=0; plt=0;
        }
      } while (plt==1);
    }
    stop();
  } 

  public void actionPerformed(ActionEvent e) {         // Buttons
    String tst;
    tst=e.getActionCommand();
    if (b1s.equals(tst)) { kk=1; }                     // Pulse
    if (b2s.equals(tst)) { if (kk>0) { plt=1; } }      // Motion
    start();
  }
}