/** applet No. 019
 *
 * template - Morse potential - molecular dynamics 2D
 *
 * Created by Ikeuchi Mitsuru on July 20 2002.
 * Copyright (c) 2002-2007 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2002.07.20  created
 * ver 0.0.2   2002.08.09  Morse potential (Di+Dj)/2 -> sqrt(Di*Dj)
 * ver 0.0.3   2002.12.11  add Hg
 * ver 0.0.4   2003.03.25  display nois reduced
 * ver 0.0.5   2003.07.17  improved in memory saving
 * ver 0.0.6   2003.11.15  add temp. control mode
 * ver 0.0.7   2004.02.04  improved pannel layout
 * ver 0.0.8   2004.06.02  improved verlet code
 * ver 0.0.9   2007.05.04  improved code
 * ver 0.1.0   2007.05.12  added setBoundary()
 * ver 0.1.1   2007.05.18  improved code
 * ver 0.1.2   2007.06.26  improved code
 *
 */

import java.awt.*;
import java.awt.event.*;
import java.applet.*;

public class Morse2D extends Applet
  implements ItemListener, ActionListener, AdjustmentListener, Runnable {

  // ------------------------------------ preset field -----------

  Thread th = null; // for run()-paint() thread

  // for event
  Choice ch_mat1, ch_mat2, ch_tempCont, ch_view;
  Button bt_reset, bt_startStop;
  Scrollbar sc_temp;

  // for off-paint buffer
  Dimension dim;
  Image imgOff;
  Graphics gOff;

  // MD2D
  double t = 0.0; // time (s)
  double dt = 5.0*1.0e-15; // time division (s)
  double xMax = 6.0e-9; // x-Box size (m)
  double yMax = 6.0e-9; // y-Box size (m)

  int tempMode = 2; // 0-adiabatic 1-scaling 2-wall cont
  double contTemp = 300.0; // control temperature (K)
  int resetSW = 0;
  int startSW = 1;

  int Nmt = 36; // number of molecules

  double AU = 1.661e-27;
  double EE = 1.602e-19;
  double AA = 1.0e-10;
  double Morse[][] = {
  // 0:mass		1:charge  2:D(eV)   3:A(A^-1)  4:r0(A)  5:color
    { 183.85 *AU, 0.0*EE, 0.9906*EE, 1.4116/AA, 3.032*AA, 0.30 }, //  0 W
    {  95.94 *AU, 0.0*EE, 0.8032*EE, 1.5079/AA, 2.976*AA, 0.33 }, //  1 Mo
    {  51.996*AU, 0.0*EE, 0.4414*EE, 1.5721/AA, 2.754*AA, 0.36 }, //  2 Cr
    {  55.847*AU, 0.0*EE, 0.4174*EE, 1.3885/AA, 2.845*AA, 0.39 }, //  3 Fe
    {  58.71 *AU, 0.0*EE, 0.4205*EE, 1.4199/AA, 2.780*AA, 0.42 }, //  4 Ni
    {  26.98 *AU, 0.0*EE, 0.2703*EE, 1.1646/AA, 3.253*AA, 0.45 }, //  5 Al
    { 207.19 *AU, 0.0*EE, 0.2348*EE, 1.1836/AA, 3.733*AA, 0.48 }, //  6 Pb
    {  63.54 *AU, 0.0*EE, 0.3429*EE, 1.3588/AA, 2.866*AA, 0.20 }, //  7 Cu
    { 107.87 *AU, 0.0*EE, 0.3323*EE, 1.3690/AA, 3.115*AA, 0.24 }, //  8 Ag
    {  40.08 *AU, 0.0*EE, 0.1623*EE, 0.8054/AA, 4.569*AA, 0.13 }, //  9 Ca
    {  87.62 *AU, 0.0*EE, 0.1513*EE, 0.7878/AA, 4.988*AA, 0.16 }, // 10 Sr
    { 137.34 *AU, 0.0*EE, 0.1416*EE, 0.6570/AA, 5.373*AA, 0.19 }, // 11 Ba
    {  22.99 *AU, 0.0*EE, 0.0633*EE, 0.5900/AA, 5.336*AA, 0.01 }, // 12 Na
    {  39.102*AU, 0.0*EE, 0.0542*EE, 0.4977/AA, 6.369*AA, 0.04 }, // 13 K
    {  85.47 *AU, 0.0*EE, 0.0464*EE, 0.4298/AA, 7.207*AA, 0.07 }, // 14 Rb
    { 132.905*AU, 0.0*EE, 0.0449*EE, 0.4157/AA, 7.557*AA, 0.10 }, // 15 Cs
    {  20.183*AU, 0.0*EE, 0.0031*EE, 1.6500/AA, 3.076*AA, 0.56 }, // 16 Ne
    {  39.948*AU, 0.0*EE, 0.0104*EE, 1.3400/AA, 3.816*AA, 0.60 }, // 17 Ar
    {  83.80 *AU, 0.0*EE, 0.0141*EE, 1.2500/AA, 4.097*AA, 0.64 }, // 18 Kr
    { 131.30 *AU, 0.0*EE, 0.0200*EE, 1.2400/AA, 4.467*AA, 0.68 }, // 19 Xe
    { 200.59 *AU, 0.0*EE, 0.0734*EE, 1.4900/AA, 3.255*AA, 0.72 }  // 20 Hg
  };

  String material[] = {
    "W" ,"Mo","Cr","Fe","Ni","Al","Pb","Cu","Ag","Ca",
    "Sr","Ba","Na","K" ,"Rb","Cs","Ne","Ar","Kr","Xe","Hg"
  };
  int NmatKind = 21;
  double MorseD[][] = new double[NmatKind][NmatKind];
  double MorseA[][] = new double[NmatKind][NmatKind];
  double Morser0[][] = new double[NmatKind][NmatKind];
  int MAT1 = 5; // Al
  int MAT2 = 0;

  int NN = Nmt;
  int kind[] = new int[NN]; // i-th particle kind
  double xx[] = new double[NN]; // i-th x-position
  double yy[] = new double[NN]; // i-th y-position
  double vx[] = new double[NN]; // i-th x-velocity
  double vy[] = new double[NN]; // i-th y-velocity
  double ffx[] = new double[NN]; // i-th x-force
  double ffy[] = new double[NN]; // i-th y-force

  // display
  int sleepTime = 30; // sleep time (ms)
  int dispWidth = 630;
  int dispHeight = 460;
  int dispMode = 0;
  double boxDispWidth = 320.0;
  double dispScale = (boxDispWidth/xMax);
  double viewScale = 1.0;


  // ------------------------------ applet main thread -----------

  public void init() {

    resize(dispWidth,dispHeight);
    setBackground(Color.white);

    dim = getSize();
    imgOff = createImage(dim.width,dim.height);
    gOff = imgOff.getGraphics();

    ch_mat1 = new Choice();
    for (int i=0; i<NmatKind; i++) {
      ch_mat1.add(material[i]);
    }
    ch_mat1.addItemListener(this);
    ch_mat1.select(material[MAT1]);

    ch_mat2 = new Choice();
    for (int i=0; i<NmatKind; i++) {
      ch_mat2.add(material[i]);
    }
    ch_mat2.addItemListener(this);
    ch_mat2.select(material[MAT2]);

    ch_tempCont = new Choice();
    ch_tempCont.add("adiabatic");
    ch_tempCont.add("T scaling"); 
    ch_tempCont.add("wall cont.");
    ch_tempCont.addItemListener(this);
    ch_tempCont.select("wall cont."); 

    ch_view = new Choice();
    ch_view.add("ball");
    ch_view.add("line");
    ch_view.add("velocity");
    ch_view.addItemListener(this);
    ch_view.select("ball");

    bt_reset= new Button("reset");
    bt_reset.addActionListener(this);

    bt_startStop= new Button("start/stop");
    bt_startStop.addActionListener(this); 
    
    sc_temp = new Scrollbar(Scrollbar.HORIZONTAL,300,10,10,1010);
    sc_temp.addAdjustmentListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
    // pnl.add(new Label("  "));
    pnl.add(ch_mat1); //pnl.add(bt_reset);
    pnl.add(bt_startStop);
    pnl.add(new Label("  "));
    pnl.add(ch_tempCont);
    pnl.add(sc_temp);
    pnl.add(ch_view);
    add(pnl,"North");

    setMorse();
    setInitialCondition();
  }

  public void start() {
    if (th == null) {
      th = new Thread(this);
      th.start();
    }
  }

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


  // ---------------------------------- event listener -----------

  public void itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_mat1) {
      MAT1 = ch_mat1.getSelectedIndex();
	  resetSW = 1;
    } else if (ev.getSource() == ch_mat2) {
      MAT2 = ch_mat2.getSelectedIndex();
	  resetSW = 1;
    } else if (ev.getSource() == ch_tempCont) {
      tempMode = ch_tempCont.getSelectedIndex();
    } else if (ev.getSource() == ch_view) {
      dispMode = ch_view.getSelectedIndex();
    }
  }

  public void actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_reset){ 
      resetSW = 1;
    } else if (ev.getSource() == bt_startStop){
      if (startSW==0) {
        startSW = 1;
      } else {
        startSW = 0;
      }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if(ev.getSource() == sc_temp) { 
      contTemp = 1.0*(double)(sc_temp.getValue());
      vAjustment();
    }
  }



  // ========================= run() - paint() loop ==============

  public void run() {
    while (th != null) {
      timeEvolution();
      offPaint();
      repaint();
      try {
        Thread.sleep(sleepTime);
      } catch (InterruptedException e) { }
    }
  }

  public void update(Graphics g) {
    paint(g);
  }

  public void paint(Graphics g) {
    g.drawImage(imgOff,0,0,this);
  }


  // --------------------------------------- off-paint -----------

  private void offPaint() {
    int gbx,gby;
    double sc;

    gOff.setColor(Color.white);
    gOff.fillRect(0,0,dim.width,dim.height);

    gbx = 30;
    gby = 90;
    sc = viewScale*dispScale;

    if (dispMode==0) {
      ballPlot(gbx,gby,sc);

    } else if (dispMode==1) {
      linePlot(gbx,gby,sc);

    } else if (dispMode==2) {
      velocityPlot(gbx,gby,sc, 0.05);

    } else if (dispMode==3) {
      ;

    }

    gOff.setColor(Color.blue);
    gOff.drawString("t="+(int)(t*1.0e13+0.5)/10.0+" ps",10,40);
    gOff.setColor(Color.black);
    if (startSW==0) {
      gOff.drawString("stopped",dispWidth/6*1+10,40);
    } else {
      gOff.drawString("started",dispWidth/6*1+10,40);
    }
    gOff.drawString("T cont.",dispWidth/6*3+10,40);
    gOff.drawString("Tc="+contTemp+" K",dispWidth/6*4+10,40);
	
    gOff.setColor(Color.blue);
    gOff.drawString("T="+(int)(temperature()*10+0.5)/10.0+" K",dispWidth/6*4+10,60);

    gOff.setColor(Color.black);
    gOff.drawString("view",dispWidth/6*5+10,40);

    gOff.drawString("Box="+(int)(xMax*1.0e10+0.5)/10.0+" x "+(int)(yMax*1.0e10+0.5)/10.0+" nm",10,60);
  }


  // ------------------------------------ plot methods -----------

  private void ballPlot(int gbx, int gby, double sc) {
    int i, ix, iy, ir, ballSize;

    drawBox(gbx, gby, sc);

    for (i=0; i<Nmt; i++) {
      ix = (int)(sc*xx[i])+gbx;
      iy = (int)(sc*yy[i])+gby;
	  ir = (int)(1.0*sc*Morse[kind[i]][4]);
      gOff.setColor(Color.getHSBColor((float)(Morse[kind[i]][5]),0.8f,0.9f ));
      gOff.fillOval(ix-ir/2,iy-ir/2, ir, ir);
    }
  }

  private void linePlot(int gbx, int gby, double sc) {
    int i,ix,iy,ir, j,jx,jy, ballSize;
    double r2,rc,rc2,col;

    drawBox(gbx, gby, sc);

    for (i=0; i<Nmt; i++) {
      ix = (int)(sc*xx[i])+gbx;
      iy = (int)(sc*yy[i])+gby;
      ir = (int)(0.4*sc*Morse[kind[i]][4]);
      for (j=0; j<Nmt; j++) {
        r2 = rsquare(i,j);
        rc = 0.5*(Morse[kind[i]][4]+Morse[kind[j]][4]);
        rc2 = 1.1*1.1*rc*rc;
        if (r2<=rc2) {
          jx = (int)(sc*xx[j])+gbx;
          jy = (int)(sc*yy[j])+gby;
          col = 0.33-9.0*(Math.sqrt(r2)/rc-1.0);
          if (col<0.01) col = 0.01;
          if (col>0.66) col = 0.66;
          gOff.setColor(Color.getHSBColor((float)col,0.8f,0.8f ));
          gOff.drawLine(ix,iy, jx,jy);
        }
      }
      gOff.setColor(Color.getHSBColor((float)(Morse[kind[i]][5]),0.9f,0.6f ));
      gOff.fillOval(ix-ir/2,iy-ir/2, ir, ir);
    }
  }

  private double rsquare(int i, int j) {
    double x,y;

    x = xx[i] - xx[j];
    y = yy[i] - yy[j];
    return x*x+y*y;
  }

  private void velocityPlot(int gbx, int gby, double sc, double mag) {
    int i, ix,iy,ir, ix2,iy2, ballSize;

    drawBox(gbx, gby, sc);

    gOff.setColor(Color.blue);
    for (i=0; i<Nmt; i++) {
      ix = (int)(sc*xx[i])+gbx;
      iy = (int)(sc*yy[i])+gby;
      ir = (int)(1.0*sc*Morse[kind[i]][4]);
      ix2 = ix+(int)(mag*vx[i]);
      iy2 = iy+(int)(mag*vy[i]);
      gOff.setColor(Color.getHSBColor((float)(Morse[kind[i]][5]),0.2f,0.9f ));
      gOff.fillOval(ix-ir/2,iy-ir/2, ir, ir);
      gOff.setColor(Color.getHSBColor((float)(Morse[kind[i]][5]),0.9f,0.3f ));
      gOff.drawLine(ix,iy, ix2, iy2);
    }
  }

  private void drawBox(int gbx, int gby, double sc) {
    int ix, iy, gball;
    double tmp;

    gball = (int)(3.4e-10*sc)/2;
    tmp = temperature();
    gOff.setColor(Color.gray); 
    if (tempMode==2) {
      gOff.setColor(Color.black); 
      if (tmp<contTemp) { gOff.setColor(Color.red); };
    }
    ix = (int)(sc*xMax);
    iy = (int)(sc*yMax);
    gOff.drawRect(gbx-gball,gby-gball,ix+2*gball,iy+2*gball);
  }



  // ======================== Molecular dymamics 2D ==============

  //  set initial condition

  private void setInitialCondition() {
    int i; 
    double s;

    t = 0.0;

    s = xMax/6.0;
    for (i=0; i<Nmt; i++) {
      kind[i] = MAT1;
      xx[i] = s*(double)(i%6)+0.5*s;
      yy[i] = s*(double)(i/6)+0.5*s;
      vx[i] = 500.0*(Math.random()-0.5);
      vy[i] = 500.0*(Math.random()-0.5);
      ffx[i] = 0.0;
      ffy[i] = 0.0;
    }
    vAjustment();
  }

  private void setMorse(){
    int i,j; 

    for (i=0; i<NmatKind; i++) {
      for (j=0; j<NmatKind; j++) {
        MorseD[i][j] = Math.sqrt(Morse[i][2]*Morse[j][2]);
        MorseA[i][j] = 0.5*(Morse[i][3]+Morse[j][3]);
        Morser0[i][j] = 0.5*(Morse[i][4]+Morse[j][4]);
      }
    }
  }


  // ---------------------------------- time evolution -----------

  private void timeEvolution() {
    int i;

    if (resetSW==1) {
      setInitialCondition();
      resetSW = 0;
    }

    if (startSW==1) {
      if (tempMode==1) {
        vAjustment();
      }
      for (i=0; i<20; i++) {
        timeStep();
      }
    }
  }

  private void timeStep() {
    int i; 
    double dtv2, a;

    dtv2 = dt/2.0;
    t = t + dt;

    for (i=0; i<Nmt; i++) {
      a = dtv2/Morse[kind[i]][0];
      vx[i] += a*ffx[i];
      vy[i] += a*ffy[i];
      xx[i] += vx[i]*dt;
      yy[i] += vy[i]*dt;
    }
    forceCalc();
    for (i=0; i<Nmt; i++) {
      a = dtv2/Morse[kind[i]][0];
      vx[i] += a*ffx[i];
      vy[i] += a*ffy[i];
    }

    setBoundary();
  }

  private void forceCalc() {
    int i,j;
    double xij,yij,rij,f,fxij,fyij;

    for(i=0;i<Nmt;i++) {
      ffx[i] = 0.0;
      ffy[i] = 0.0;
    }

    for(i=0;i<Nmt;i++) {
      for(j=i+1;j<Nmt;j++) {
        xij = xx[i]-xx[j];
        yij = yy[i]-yy[j];
        rij = Math.sqrt(xij*xij+yij*yij);
        f = force(rij,kind[i],kind[j]);
        fxij = f*xij/rij;
        fyij = f*yij/rij;

        ffx[i] += fxij;
        ffy[i] += fyij;

        ffx[j] += -fxij;
        ffy[j] += -fyij;
      }
    }
  }

  private double force(double r,int ki, int kj) {
    double d,a,r0,y;

    d = MorseD[ki][kj];
    a = MorseA[ki][kj];
    r0 = Morser0[ki][kj];
    y = Math.exp(-a*(r-r0));
    return ( 2.0*d*a*y*(y-1.0) );
  }

  private void setBoundary() {
    int i; 
    double rr, tmp;

    rr = 1.0;
    if (tempMode==2) {
      tmp = temperature();
      rr = Math.sqrt(contTemp/tmp);
      if (rr<0.5) { rr = 0.5; } else if (rr>1.2) { rr=1.2; };
    }

    for (i=0; i<Nmt; i++) {
      if (xx[i] < 0.0) { 
        xx[i] = 0.0; vx[i] = -rr*vx[i]; vy[i] = rr*vy[i]; 
      }
      if (xx[i] > xMax) {
        xx[i] = xMax; vx[i] = -rr*vx[i]; vy[i] = rr*vy[i]; 
      }
      if (yy[i] < 0.0) { 
        yy[i] = 0.0; vx[i] = rr*vx[i]; vy[i] = -rr*vy[i]; 
      }
      if (yy[i] > yMax) {
        yy[i] = yMax; vx[i] = rr*vx[i]; vy[i] = -rr*vy[i];
      }
    }
  }


  // --------------------------------------- utilities -----------

  private double temperature() {
    int i;
    double s;

    s = 0.0; 
    for (i=0; i<Nmt; i++) {
      s += 0.5*Morse[kind[i]][0]*(vx[i]*vx[i] + vy[i]*vy[i]);
    }
    return ( s/((double)Nmt*1.38e-23) );
  }

  //

  private void vAjustment() {
    int i;
    double tmp, r;

    tmp = temperature();
    r = Math.sqrt(contTemp/tmp);
    for (i=0; i<Nmt; i++) {
      vx[i] = r*vx[i];
      vy[i] = r*vy[i];
    }
  }

}