/** applet No. 1023
 *
 * template - molecular dynamics 3D
 * - multi thread  MD :: display - synchronous
 *  ( MD free run and notice --> sycn display )
 *
 * Created by Ikeuchi Mitsuru on September 18 2006.
 * Copyright (c) 2006-2007 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2006.09.18  created
 * ver 0.0.2   2007.06.08  improved code
 *
 */

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

public class syncMTmd3d extends Applet
  implements MouseListener, MouseMotionListener,
    ItemListener, ActionListener, AdjustmentListener, Runnable {

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

  Thread th = null; // for run()-paint() thread
  MolecularDynamics3D md = null;
  Notice ntc = new Notice();

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

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

  int dispWidth = 630;
  int dispHeight = 420;

  int sleepTime = 50; // in ms
  int mdSleepTime = 10;

  int thCount = 0;

  int dispMode = 0; // 0-ball 1-near
  int tempMode = 2; // 0-adiabatic 1-scaling 2-wall cont
  double contTemp = 300.0;

  double dispScale = (250.0/6.0e-9); // 250(points)/6(nm) scale
  double viewScale = 1.0;

  int dgX, dgY, dgXb, dgYb;
  double theta = -20.0*3.14/180.0;
  double fai = 10.0*3.14/180.0;
  double dtheta = 0.0*3.14/180.0;
  double pai = 3.1415926536;

  // 3d display
  int NN = 2000;
  int srtz[][] = new int[2][NN];
  double px[] = new double[NN]; 
  double py[] = new double[NN]; 
  double pz[] = new double[NN]; 
  double wkx[] = new double[8]; 
  double wky[] = new double[8]; 
  double wkz[] = new double[8]; 
  double pwkx[] = new double[8]; 
  double pwky[] = new double[8]; 
  double pwkz[] = new double[8]; 
  int boxp[][] = {
   {0,1},{0,2},{0,4},{1,3},{1,5},{2,3},{2,6},{4,5},{4,6},{3,7},{5,7},{6,7}
  };


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

  public void init() {

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

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

    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("near");
    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_mdSleep = new Scrollbar(Scrollbar.HORIZONTAL,10,10,5,110);
    sc_mdSleep.addAdjustmentListener(this);

    sc_temp = new Scrollbar(Scrollbar.HORIZONTAL,300,10,10,610);
    sc_temp.addAdjustmentListener(this);

    // 3d control
    addMouseListener(this);  
    addMouseMotionListener(this);

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

  }

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

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


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

  public void itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_tempCont) {
      tempMode = ch_tempCont.getSelectedIndex();
      md.setTempMode(tempMode);
    } else if (ev.getSource() == ch_view) {
      dispMode = ch_view.getSelectedIndex();
    }
  }

  public void actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_reset){ 
      md.reset();
      thCount = 0;
    } else if (ev.getSource() == bt_startStop){
      if (md.getStartSW()==0) {
        md.setStartSW(1);
      } else {
        md.setStartSW(0);
      }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_mdSleep) { 
      mdSleepTime = sc_mdSleep.getValue();
      md.setSleepTime(mdSleepTime);
    } else if (ev.getSource() == sc_temp) { 
      contTemp = 1.0*(double)(sc_temp.getValue());
	  md.setTemperature(contTemp);
    }
  }

  // 3d control - mouse
  public void mousePressed(MouseEvent ev){ }
  public void mouseReleased(MouseEvent ev){
    dgXb = 0; dgYb = 0;
    dgX = 0; dgY = 0;
  }
  public void mouseClicked(MouseEvent ev){ }
  public void mouseEntered(MouseEvent ev){ } 
  public void mouseExited (MouseEvent ev){ }
  public void mouseMoved  (MouseEvent ev){ }
  public void mouseDragged(MouseEvent ev){
    dgXb = dgX; dgYb = dgY;
    dgX=ev.getX();
    dgY=ev.getY();
    if (dgXb==0 && dgYb==0) {
      dgXb = dgX; dgYb = dgY;
    }
    theta += 0.5*3.14/180.0*(dgX-dgXb);
    fai += 0.5*3.14/180.0*(dgY-dgYb);
  }



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

  public void run() {
    while (th != null) {
	  thCount += 1;
      offPaint();
      repaint();
      try {
        Thread.sleep(sleepTime);
      } catch (InterruptedException e) { }
      ntc.getup();
    }
  }

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

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


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

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

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

    gbx = 80;
    gby = 100;

    if (dispMode==0) {
      ball3dPlot(gbx,gby,1.0,0);

    } else if (dispMode==1) {
      ball3dPlot(gbx,gby,0.5,1);

    } else if (dispMode==2) {
      ;

    }

    gOff.setColor(Color.black);
    gOff.drawString("t="+(int)(md.getTime()*1.0e13+0.5)/10.0+" ps",10,40);
    if (md.getStartSW()==0) {
      gOff.drawString("stopped",630/6*1+10,40);
    } else {
      gOff.drawString("started",630/6*1+10,40);
    }
    gOff.drawString("md sleep="+mdSleepTime+" ",630/6*2+10,40);
    gOff.drawString("T cont.",dispWidth/6*3+10,40);
    gOff.drawString("Tc="+contTemp+" K",dispWidth/6*4+10,40);
    gOff.drawString("T="+(int)(md.temperature()*10+0.5)/10.0+" K",dispWidth/6*4+10,60);
    gOff.drawString("view",dispWidth/6*5+10,40);

    gOff.drawString("Box="+(int)(md.getxMax()*1.0e10+0.5)/10.0+
	  " x "+(int)(md.getyMax()*1.0e10+0.5)/10.0+
	  " x "+(int)(md.getzMax()*1.0e10+0.5)/10.0+" nm",dispWidth/6*4,100);

    gOff.setColor(Color.black);
    gOff.drawString("loop count",420,120);
    gOff.drawString("thread md ="+md.getcount()+"",420,140);
    gOff.drawString("thread disp ="+thCount+"",420,160);
    gOff.drawString("disp sleep = 50 ms",420,180);
  }


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

  void ball3dPlot(int gbx, int gby, double size, int mode) {
    int nn,i,j,jj,gx,gy,gz,g2x,g2y;
    double zm,sc,sz,sgm;

    setWaku();
    rotate(theta,fai);
    zSort();

    drawWaku(gbx,gby,0);

    zm = md.getzMax();
    sgm = md.getsgm();
    sc = dispScale*viewScale;

    nn = md.getNmt();
    for (i=0; i<nn; i++) {
      j = srtz[0][i];
      sz = sc*size*sgm;
      gx = (int)(sc*px[j])+gbx;
      gy = (int)(sc*py[j])+gby;
      gz = (int)(sz*(0.3*pz[j]/zm+0.7));

      if (mode==1) {
        for(jj=j+1;jj<nn;jj++) {
          if (md.distance(j,jj)<1.3*sgm) {
            g2x = (int)(sc*px[jj])+gbx;
            g2y = (int)(sc*py[jj])+gby;
            gOff.drawLine(gx,gy, g2x, g2y);
          }
        }
      }
      gOff.setColor(Color.getHSBColor(0.15f,0.8f,(float)(0.4*pz[j]/zm+0.4)));
      gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
    }
    drawWaku(gbx,gby,1);
  }

  //

  void drawWaku(int gbx, int gby, int imode) {
    int i,gx,gy,g2x,g2y,farPoint;
    double sc,tmp;

    sc = dispScale*viewScale;
    farPoint = findMin();

    for (i=0; i<12; i++) {
      gx = (int)(sc*pwkx[boxp[i][0]])+gbx;
      gy = (int)(sc*pwky[boxp[i][0]])+gby;
      g2x = (int)(sc*pwkx[boxp[i][1]])+gbx;
      g2y = (int)(sc*pwky[boxp[i][1]])+gby;
      if (imode==0) {
        if (boxp[i][0]==farPoint || boxp[i][1]==farPoint) { 
          gOff.setColor(Color.lightGray);
          gOff.drawLine(gx,gy, g2x, g2y);
        }
      } else if (imode==1) {
        if (boxp[i][0]!=farPoint && boxp[i][1]!=farPoint) { 
          gOff.setColor(Color.gray);
          if (tempMode==2) {
            tmp = md.temperature();
            if (tmp>contTemp) {
              gOff.setColor(Color.black); 
            } else {
              gOff.setColor(Color.red);
            }
          }
          gOff.drawLine(gx,gy, g2x, g2y);
        }
      }
    }
  }

  int findMin() {
    int i,im;
    double m;

    im = 0; m = pwkz[im];
    for (i=0; i<8; i++) {
      if (pwkz[i]<m) {
        m = pwkz[i];
        im = i;
      }
    }
    return(im);
  }

  //

  void setWaku() {
    int i,ix,iy,iz;
    double sgm,xm,ym,zm;

    i=0;
    sgm = md.getsgm();
    xm = md.getxMax(); ym = md.getyMax(); zm = md.getzMax();

    for (ix=0; ix<2; ix++) {
      for (iy=0; iy<2; iy++) {
        for (iz=0; iz<2; iz++) {
          wkx[i] = -0.5*sgm+(xm+sgm)*(double)(ix);
          wky[i] = -0.5*sgm+(ym+sgm)*(double)(iy);
          wkz[i] = -0.5*sgm+(zm+sgm)*(double)(iz);
          i += 1;
        }
      }
    }
  }

  //

  void rotate(double th,double fi) {
    int i,nn;
    double cosTh,sinTh,cosFi,sinFi,cx,cy,cz;
    double x,y,z,mx,my,mz;

    cosTh = Math.cos(th); sinTh = Math.sin(th);
    cosFi = Math.cos(fi); sinFi = Math.sin(fi);

    cx = 0.5*md.getxMax();
    cy = 0.5*md.getyMax();
    cz = 0.5*md.getzMax();

    nn = md.getNmt();
	for (i=0; i<nn; i++) {
	  my = md.getyy(i); mz = md.getzz(i);
      py[i] = cosFi*(my-cy)+sinFi*(mz-cz) + cy; 
      pz[i] = -sinFi*(my-cy)+cosFi*(mz-cz) + cz;
    } 
    for (i=0; i<nn; i++) {
	  mx = md.getxx(i);
      px[i] = cosTh*(mx-cx)+sinTh*(pz[i]-cz) + cx; 
      pz[i] = -sinTh*(mx-cx)+cosTh*(pz[i]-cz) + cz;
    } 

    for (i=0; i<8; i++) {
      pwky[i] = cosFi*(wky[i]-cy)+sinFi*(wkz[i]-cz) + cy; 
      pwkz[i] = -sinFi*(wky[i]-cy)+cosFi*(wkz[i]-cz) + cz;
    } 
    for (i=0; i<8; i++) {
      pwkx[i] = cosTh*(wkx[i]-cx)+sinTh*(pwkz[i]-cz) + cx; 
      pwkz[i] = -sinTh*(wkx[i]-cx)+cosTh*(pwkz[i]-cz) + cz;
    }  
  }

  //

  void zSort() {
    int nn,i,j,w;
    double zm;

    zm = md.getzMax();
    nn = md.getNmt();
    for (i=0; i<nn; i++) {
      srtz[0][i] = i;
      srtz[1][i] = (int)(32767.0*pz[i]/zm);
    }
    qSort(0,nn-1);
  }

  void qSort(int i,int j) {
    int k,p;

    if(i==j) return;
    p = i+1;
    while (p<=j && srtz[1][i]==srtz[1][p]) p++;
    if (p<=j) {
      if (srtz[1][i]>=srtz[1][p]) p = i;
      k = qSortPartition(i,j,srtz[1][p]);
      qSort(i,k-1);
      qSort(k,j);
    }
  }

  int qSortPartition(int i,int j, int x) {
    int l,r,w;

    l = i; r = j;
    while (l<=r) {
      while (l<=j && srtz[1][l]<x) l++;
      while (r>=i && srtz[1][r]>=x) r--;
      if (l>r) break;
      w = srtz[0][l]; srtz[0][l] = srtz[0][r]; srtz[0][r] = w;
      w = srtz[1][l]; srtz[1][l] = srtz[1][r]; srtz[1][r] = w;
      l++; r--;
    }
    return ( l );
  }

}



// ======== notice ( MD free run and notice --> sycn display ) ===========

class Notice { 
   private int flag = 1;

   public synchronized void attention() { 
      notifyAll(); 
      flag = 0; 
   }

   public synchronized void getup() {
      while (flag==1) { 
         try { 
            wait(); 
         } catch (InterruptedException e) { } 
      } 
      flag = 1; 
   } 
}



// ======================= molecular dynamics 3D class ===========
 
class MolecularDynamics3D extends Thread {

  int NN = 2000;
  double xx[] = new double[NN];
  double yy[] = new double[NN];
  double zz[] = new double[NN];
  double vx[] = new double[NN];
  double vy[] = new double[NN];
  double vz[] = new double[NN];
  double ffx[] = new double[NN];
  double ffy[] = new double[NN];
  double ffz[] = new double[NN];

  double t = 0.0;
  double dt = 10.0e-15;
  double xMax = 6.0e-9;
  double yMax = 6.0e-9;
  double zMax = 6.0e-9;

  int sleepTime = 10; // in ms
  int waitSleepTime = 50;
  int Nmt = 0;
  int resetSW = 0; // 1-reset
  int startSW = 1; // 0-stop 1-start
  int tempMode = 2; // 0-adiabatic 1-scaling 2-wall cont
  double contTemp = 300.0; // in K

  double sgm = 3.40e-10;        // sigma of L-J potential for Ar (m)
  double eps = 0.0104*1.602e-19;// epsilon of L-J potential for Ar (J)
  double mas = 39.948*1.661e-27;// mass of L-J potential for Ar (kg)

  Notice note;
  int mdCount = 0;


  // ------------------------------- class constructor -----------

  MolecularDynamics3D(Notice ntc) {
    note = ntc;
    setInitialCondition();
  }


  // -------------------------------------- thread run -----------

  public void run() {
    while (true) {
      mdCount += 1;
      timeEvolution();
      try {
        this.sleep(sleepTime);
      } catch (InterruptedException e) { }
      note.attention();
    }
  }


  // --------------------------- set initial condition -----------

  void setInitialCondition() {
    int i,ix,iy,iz; 
    double s;

    t = 0.0;
    mdCount = 0;

    s = xMax/5.0;
	i = 0;
    for (ix=0; ix<5; ix++) {
      for (iy=0; iy<5; iy++) {
        for (iz=0; iz<5; iz++) {
          xx[i] = s*(ix+0.5);
          yy[i] = s*(iy+0.5);
          zz[i] = s*(iz+0.5);
          vx[i] = 500.0*(Math.random()-0.5);
          vy[i] = 500.0*(Math.random()-0.5);
          vz[i] = 500.0*(Math.random()-0.5);
          ffx[i] = 0.0;
          ffy[i] = 0.0;
          ffz[i] = 0.0;
          i += 1;
        }
      }
    }
	Nmt = i;
    vAjustment(contTemp);
  }


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

  void timeEvolution() {
    int i;

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

    if (startSW==1) {

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

    } else {
      try {
        this.sleep(waitSleepTime);
      } catch (InterruptedException e) { }
    }
  }

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

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

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

    setBoundary();
  }

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

    for(i=0;i<Nmt;i++) {
      ffx[i] = 0.0; ffy[i] = 0.0; ffz[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]; zij = zz[i]-zz[j];
        rij = Math.sqrt(xij*xij + yij*yij + zij*zij);
        f = force(rij);
        fxij = f*xij/rij;
        fyij = f*yij/rij;
        fzij = f*zij/rij;

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

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

  double force(double r) {
    double ri,r3,r6;

    ri = sgm/r;
    r3 = ri*ri*ri;
    r6 = r3*r3;
    return (24.0*eps*r6*(2.0*r6-1.0)/r);
  }


  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]; vz[i] = rr*vz[i];
      }
      if (xx[i] > xMax) {
        xx[i] = xMax; vx[i] = -rr*vx[i]; vy[i] = rr*vy[i]; vz[i] = rr*vz[i]; 
      }
      if (yy[i] < 0.0) { 
        yy[i] = 0.0; vx[i] = rr*vx[i]; vy[i] = -rr*vy[i]; vz[i] = rr*vz[i]; 
      }
      if (yy[i] > yMax) {
        yy[i] = yMax; vx[i] = rr*vx[i]; vy[i] = -rr*vy[i]; vz[i] = rr*vz[i];
      }
      if (zz[i] < 0.0) { 
        zz[i] = 0.0; vx[i] = rr*vx[i]; vy[i] = rr*vy[i]; vz[i] = -rr*vz[i];
      }
      if (zz[i] > zMax) {
        zz[i] = zMax; vx[i] = rr*vx[i]; vy[i] = rr*vy[i]; vz[i] = -rr*vz[i];
      }
    }
  }


  // ----------------------------------------- utility -----------

  void vAjustment(double temp) {
    int i;
    double tp, r;

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

  double temperature() {
    int i;
    double s,ke;

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

  double distance(int i,int j) {
    double x,y,z;

    x = xx[i]-xx[j];
    y = yy[i]-yy[j];
    z = zz[i]-zz[j];

    return(Math.sqrt(x*x+y*y+z*z));
  }


  // ------------------------------------- I/O methods -----------

  // control
  public void reset() { resetSW = 1; }
  public void setStartSW(int sw) { startSW = sw; }
  public int getStartSW() { return startSW; }
  public int getNmt() { return Nmt; }
  public void setSleepTime(int mdSleepTime) { sleepTime = mdSleepTime; }
  public int getcount() { return mdCount; }

  public void setTempMode(int mode) { tempMode = mode; }
  public void setTemperature(double temp) { contTemp = temp; }

  // get data
  public double getTime() { return t; }
  public double getxMax() { return xMax; }
  public double getyMax() { return yMax; }
  public double getzMax() { return zMax; }

  public double getxx(int i) { return xx[i]; }
  public double getyy(int i) { return yy[i]; }
  public double getzz(int i) { return zz[i]; }
  public double getvx(int i) { return vx[i]; }
  public double getvy(int i) { return vy[i]; }
  public double getvz(int i) { return vz[i]; }

  public double getsgm() { return sgm; }

}