/** applet No. 1255
 *
 * template - Langevin dynamics 3D
 * - Lennard-Jones potential
 * - multi thread  LD3D :: display - asynchronous
 *
 * Created by Ikeuchi Mitsuru on February 02 2008.
 * Copyright (c) 2008 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2008.02.02  created
 *
 */

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

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

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

  Thread th = null; // for run()-paint() thread
  Thread thmd = null; // for LD3D() thread
  LJMolecule ma = new LJMolecule();
  LD3D md = new LD3D(ma);
  DisplayLD3D dsp = new DisplayLD3D();

  // for event
  Choice ch_mat1, ch_mat2, ch_tempControl, ch_view;
  Button bt_reset, bt_startStop, bt_step;
  Scrollbar sc_temp, sc_scale;

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

  int dgX,dgY, dgXb,dgYb;

  int dispMode = 3;
  int sleepTime = 100;
  int thCount = 0;


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

  public void init() {
    resize(630,460);
    setBackground(Color.white);

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

    ch_mat1 = new Choice();
    for (int i=0; i<ma.getKindMax(); i++) {
      ch_mat1.add(ma.getMatStr(i));
    }
    ch_mat1.addItemListener(this);
    ch_mat1.select(ma.getMatStr(1));

    ch_mat2 = new Choice();
    for (int i=0; i<ma.getKindMax(); i++) {
      ch_mat2.add(ma.getMatStr(i));
    }
    ch_mat2.addItemListener(this);
    ch_mat2.select(ma.getMatStr(0));

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

    ch_view = new Choice();
    ch_view.add("ball");
    ch_view.add("line"); 
    ch_view.add("ball+line");
    ch_view.add("real size");
    ch_view.add("temperature");
    ch_view.add("velocity");
    ch_view.add("vSpace");
    ch_view.add("line+rdf");
    ch_view.add("v+vdf");
    ch_view.addItemListener(this);
    ch_view.select("real size");

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

    bt_startStop= new Button("start/stop");
    bt_startStop.addActionListener(this); 

    bt_step = new Button("step");
    bt_step.addActionListener(this);

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

    sc_scale= new Scrollbar(Scrollbar.HORIZONTAL,100,10,50,210);
    sc_scale.addAdjustmentListener(this);

    addMouseListener(this);  
    addMouseMotionListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(2,6,5,0));
    //pnl.add(new Label("  "));
    pnl.add(bt_reset);
    pnl.add(bt_startStop);
    pnl.add(bt_step);
    pnl.add(new Label("  "));//pnl.add(ch_tempControl);
    pnl.add(new Label("  "));
    pnl.add(new Label("  "));

    pnl.add(ch_mat1);
    pnl.add(new Label("  "));
    pnl.add(new Label("  "));
    pnl.add(sc_temp);
    pnl.add(sc_scale);
    pnl.add(ch_view);
    add(pnl,"North");
  }

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

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


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

  public void itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_mat1){
      md.setMAT1( ch_mat1.getSelectedIndex() );
      thCount = 0;
    } else if (ev.getSource() == ch_mat2){
      md.setMAT2( ch_mat2.getSelectedIndex() );
      thCount = 0;
    } else if (ev.getSource() == ch_tempControl) {
      md.setTempMode( ch_tempControl.getSelectedIndex() );
    } 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);
      }
    } else if (ev.getSource() == bt_step){
      if (md.getStartSW()==0) {
        md.setStepSW();
      }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_temp) { 
      md.setContTemp( 1.0*(double)(sc_temp.getValue()) );
    } else if(ev.getSource() == sc_scale) { 
      dsp.setViewScale( 0.01*(double)(sc_scale.getValue()) );
    }
  }

  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;
    }
    dsp.setViewDelta((dgX-dgXb), (dgY-dgYb));
  }



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

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

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

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

  private void offPaint() {
    if (md != null) {
      dsp.plotGraph(gOff, md, dispMode, thCount);
    }
  }

}  



// ================================ display MD3D class ===========

class DisplayLD3D {

  Graphics gOff;
  LD3D md;

  private int Nmt;

  private double dt = md.dt;
  private double xMax = md.xMax;
  private double yMax = md.yMax;
  private double zMax = md.zMax;

  private double zDepth = xMax;
  private double dispWidth = 250.0;
  private double dispScale = (dispWidth/xMax);
  private double viewScale = 1.0;

  private double cx = 0.5*xMax;
  private double cy = 0.5*yMax;
  private double cz = 0.5*zMax;
  private double theta = -20.0*3.14/180.0;
  private double fai = 10.0*3.14/180.0;
  private double cosTh = Math.cos(theta);
  private double sinTh = Math.sin(theta);
  private double cosFi = Math.cos(fai);
  private double sinFi = Math.sin(fai);


  private int NN = 10000;
  private int srtz[][] = new int[2][NN];

  private double px[] = new double[NN]; 
  private double py[] = new double[NN]; 
  private double pz[] = new double[NN];

  private double pxx[] = new double[NN]; 
  private double pyy[] = new double[NN]; 
  private double pzz[] = new double[NN]; 

  private double wkx[] = new double[8]; 
  private double wky[] = new double[8]; 
  private double wkz[] = new double[8]; 

  private double pwkx[] = new double[8]; 
  private double pwky[] = new double[8]; 
  private double pwkz[] = new double[8]; 

  private 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}
  };

  private int rdf[] = new int[120];
  int vdf[][] = new int[4][220];
  double meanvdf[][] = new double[4][220];
  double ffvMaxwell[] = new double[220];

  private int drawnCount = 0;


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

  DisplayLD3D() {
    setWaku();
  }


  // ----------------------------------- class methods -----------

  public void plotGraph(Graphics gg, LD3D md3d, int dispMode, int thCount) {
    int mdCount, gbx, gby;
    double tmp;

    gOff = gg;
    md = md3d;

    if (md==null) return;
    mdCount = md.getCount();
    if (mdCount==drawnCount) return;
    drawnCount = mdCount;

    Nmt = md.getNmt();

    gOff.setColor(Color.white);
    gOff.fillRect(0,0,630,460);

    gbx = 70;
    gby = 130;

    if (dispMode==0) { // ball
      ballPlot(gbx,gby, 0.7, 0);

    } else if (dispMode==1) { // line
      ballPlot(gbx,gby, 0.5, 1);

    } else if (dispMode==2) { // ball + line
      ballPlot(gbx,gby, 0.5, 2);

    } else if (dispMode==3) { // real
      ballPlot(gbx,gby, 1.0, 0);

    } else if (dispMode==4) { // temperature
      ballPlot(gbx,gby, 0.7, 3);

    } else if (dispMode==5) { // velocity
      velocityPlot(gbx, gby, 10.0);

    } else if (dispMode==6) { // velocity
      vSpacePlot(gbx, gby, 1000.0);

    } else if (dispMode==7) {
      ballPlot(gbx,gby, 0.5, 1);
      rdfPlot(400,180, 0.5);

    } else if (dispMode==8) {
      velocityPlot(gbx, gby, 10.0);
      vdfPlot(400,180, 1.0);
    }


    tmp = md.temperature();
    gOff.setColor(Color.blue);
    gOff.drawString("t="+(int)(md.getTime()*1.0e13+0.5)/10.0+" ps",630/6*0+10,70);
    gOff.setColor(Color.black);

    gOff.drawString("Tc="+md.getContTemp()+"K",630/6*3+10,70);
    gOff.drawString("scale="+(int)(viewScale*100.0+0.499)+"%",630/6*4+10,70);
    gOff.drawString("view",630/6*5+10,70);

    gOff.setColor(Color.blue);
    gOff.drawString("T="+(int)(tmp*10+0.5)/10.0+" K",630/6*3+10,90);
    gOff.drawString("L="+(int)(xMax*1.0e10+0.5)/10.0+" nm",630/6*4+10,90);

    gOff.setColor(Color.black);
    gOff.drawString("regMax="+md.regMax()+"",630/6*5+10,130);
    gOff.drawString("secMax="+md.secMax()+"",630/6*5+10,150);

    gOff.setColor(Color.black);
    gOff.drawString("thread MD = "+md.getCount()+" ",20,440);
    gOff.drawString("thread disp = "+thCount+" ",220,440);
  }


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

  private void ballPlot(int gbx, int gby, double size, int imode) {
    int i,j,k,jj,gx,gy,g2x,g2y,gr;
	double sz,sc,tm;

    setPos();
    rotate();
    zSort();
    drawWaku(gbx,gby,0);

    sc = dispScale*viewScale;
    sz = size*sc;
    for (i=0; i<Nmt; i++) {
      j = srtz[0][i];
      gx = (int)(sc*px[j])+gbx;
      gy = (int)(sc*py[j])+gby;
      gr = (int)(sz*md.getDia(j)*(0.2*(pz[j]/zDepth)+0.8));
      if (gr<2) gr = 2;
      if (imode==1 || imode==2) {
        for (k=1; k<md.getnnReg(j); k++) {
          jj = md.getReg(j,k);
          if ( distance(j,jj)<0.6*(md.getDia(j)+md.getDia(jj)) ) {
            g2x = (int)(sc*px[jj])+gbx;
            g2y = (int)(sc*py[jj])+gby;
            gOff.setColor(Color.gray);
            gOff.drawLine(gx,gy, g2x, g2y);
          }
        }
      }
      if (imode==0 || imode==2) {
        gOff.setColor(Color.getHSBColor(
          (float)md.getColor(j),0.8f,(float)(0.2*pz[j]/zDepth+0.5)));
        gOff.fillOval(gx-gr/2,gy-gr/2, gr, gr);
      }
      if (imode==3) {
        tm = md.getTemp(j)/1000.0;
        if (tm>0.99) tm = 0.99;
        gOff.setColor(Color.getHSBColor(
          (float)md.getColor(j),0.8f,(float)(0.8*tm+0.2)));
        gOff.fillOval(gx-gr/2,gy-gr/2, gr, gr);
      }

    }
    drawWaku(gbx,gby,1);
  }

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

    x = px[i]-px[j];
    y = py[i]-py[j];
    z = pz[i]-pz[j];
    return Math.sqrt(x*x+y*y+z*z);
  }


  private void velocityPlot(int gbx, int gby, double mag) {
    int i,j,gx,gy,g2x,g2y;
    double sc,m, p2x,p2y,p2z,rotp2x,rotp2y,rotp2z;

    setPos();
    rotate();
    zSort();
    drawWaku(gbx,gby,0);

    sc = dispScale*viewScale;

    for (i=0; i<Nmt; i++) {
      j = srtz[0][i];
      gx = (int)(sc*px[j])+gbx;
      gy = (int)(sc*py[j])+gby;

      p2x = pxx[j] + mag*md.getvx(j)*20.0*dt;
      p2y = pyy[j] + mag*md.getvy(j)*20.0*dt;
      p2z = pzz[j] + mag*md.getvz(j)*20.0*dt;
      rotp2y = cosFi*(p2y-cy)+sinFi*(p2z-cz) + cy; 
      rotp2z = -sinFi*(p2y-cy)+cosFi*(p2z-cz) + cz;
      rotp2x = cosTh*(p2x-cx)+sinTh*(rotp2z-cz) + cx; 
      rotp2z = -sinTh*(p2x-cx)+cosTh*(rotp2z-cz) + cz;
      g2x = (int)(sc*rotp2x)+gbx;
      g2y = (int)(sc*rotp2y)+gby;

      gOff.setColor(Color.getHSBColor(
        (float)md.getColor(j),0.8f,(float)(0.2*pz[j]/zDepth+0.5)));
      gOff.drawLine(gx,gy, g2x,g2y);
    }

    drawWaku(gbx,gby,1);
  }

  private void vSpacePlot(int gbx, int gby, double mag) {
    int i,j,gx,gy,gr;
    double sc;

    setVelocity(mag);
    rotate();
    zSort();
    drawWaku(gbx,gby,0);

    sc = dispScale*viewScale;

    for (i=0; i<Nmt; i++) {
      j = srtz[0][i];
      gx = (int)(sc*px[j])+gbx;
      gy = (int)(sc*py[j])+gby;
      gr = (int)(5.0*viewScale);
      if (pxx[j]>=xMax/2.0) {
        gOff.setColor(Color.getHSBColor((float)md.getColor(j),0.8f,(float)(0.2*pz[j]/zDepth+0.5)));
      } else {
        gOff.setColor(Color.getHSBColor((float)md.getColor(j),0.4f,(float)(0.2*pz[j]/zDepth+0.5)));
	  }
      gOff.drawOval(gx-gr/2,gy-gr/2, gr, gr);
    }

    drawWaku(gbx,gby,1);
  }

  // plot utility

  private 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 (md.getTempMode()==2) {
            tmp = md.temperature();
            if (tmp>md.getContTemp()) {
              gOff.setColor(Color.black); 
            } else {
              gOff.setColor(Color.red);
            }
          }
          gOff.drawLine(gx,gy, g2x, g2y);
        }
      }
    }
  }

  private 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);
  }

  // setWaku

  private void setWaku() {
    int i,ix,iy,iz;
    double sgm;

    sgm = 4.0e-10;
    i=0;
    for (ix=0; ix<2; ix++) {
      for (iy=0; iy<2; iy++) {
        for (iz=0; iz<2; iz++) {
          wkx[i] = -0.5*sgm+(xMax+sgm)*(double)(ix);
          wky[i] = -0.5*sgm+(yMax+sgm)*(double)(iy);
          wkz[i] = -0.5*sgm+(zMax+sgm)*(double)(iz);
          i += 1;
        }
      }
    }
  }

  // sort z

  private void zSort() {
    int i,j,w;

    for (i=0; i<Nmt; i++) {
      srtz[0][i] = i;
      srtz[1][i] = (int)(32767.0*pz[i]/zMax);
    }
    qSort(0,Nmt-1);    
  }

  private 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);
    }
  }

  private 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 );
  }

  // rotate

  private void rotate() {
    int i;

    for (i=0; i<Nmt; i++) {
      py[i] = cosFi*(pyy[i]-cy)+sinFi*(pzz[i]-cz) + cy; 
      pz[i] = -sinFi*(pyy[i]-cy)+cosFi*(pzz[i]-cz) + cz;
      px[i] = cosTh*(pxx[i]-cx)+sinTh*(pz[i]-cz) + cx; 
      pz[i] = -sinTh*(pxx[i]-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;
      pwkx[i] = cosTh*(wkx[i]-cx)+sinTh*(pwkz[i]-cz) + cx; 
      pwkz[i] = -sinTh*(wkx[i]-cx)+cosTh*(pwkz[i]-cz) + cz;
    }  
  }


  private void setPos() {
    int i;

    for (i=0; i<Nmt; i++) {
      pxx[i] = md.getxx(i);
      pyy[i] = md.getyy(i);
      pzz[i] = md.getzz(i);
    } 
  }

  private void setVelocity(double mag) {
    int i;

    for (i=0; i<Nmt; i++) {
      pxx[i] = xMax/2.0 + mag*md.getvx(i)*dt;
      pyy[i] = yMax/2.0 + mag*md.getvy(i)*dt;
      pzz[i] = zMax/2.0 + mag*md.getvz(i)*dt;
    } 
  }

  // radial distribution set / plot

  private void rdfPlot(int gx, int gy, double mag) {
    int i,iy,hoganWidth,hoganHight,xDiv,yDiv;

    setRadialDistribution();

    hoganWidth = 200; hoganHight = 200;
    xDiv = hoganWidth/5; yDiv = hoganHight/2;

    gOff.setColor(Color.cyan);
    for (i=0; i<100; i++) {
      iy = (int)(mag*rdf[i]*(5000.0/(i*i+0.25)));
      if (iy>0) {
        gOff.fillRect(gx+i*2,gy+hoganHight-iy,2,iy);
      }
    }

    drawHogan(gx, gy, hoganWidth, hoganHight, xDiv, yDiv);
    gOff.setColor(Color.black);
    gOff.drawString("0",gx-5,gy+hoganHight+20);
    gOff.drawString("1 (nm)",gx+hoganWidth-20,gy+hoganHight+20);
    gOff.drawString("radial distribution",gx+40,gy+hoganHight+40);
  }

  private void drawHogan(int gx, int gy, int hoganWidth, int hoganHight, int xDiv, int yDiv) {
    int i;

    gOff.setColor(Color.gray);
    for (i=0; i<hoganWidth; i+=xDiv) {
      gOff.drawLine(gx+i,gy, gx+i, gy+hoganHight);
    }
    for (i=0; i<hoganHight; i+=yDiv) {
      gOff.drawLine(gx,gy+hoganHight-i, gx+hoganWidth, gy+hoganHight-i);
    }
    gOff.setColor(Color.black);
    gOff.drawRect(gx,gy,hoganWidth,hoganHight);
  }

  private void setRadialDistribution() {
    int i,j,k,ir;
    double rij;

    for (i=0; i<120; i++) {
      rdf[i]=0;
    }

    for (i=0; i<Nmt; i++) {
      for(k=1;k<md.getnnReg(i);k++) {
        j = md.getReg(i,k);
        rij = distance(i,j);
        if (rij<10.0e-10) {
          ir = (int)(rij*1.0e11+0.49999);
          rdf[ir] +=1;
        }
      }
    }
  }

  // verocity distribution function

  private void vdfPlot(int gx, int gy, double mag) {
    int i,iy,iy2,hoganWidth,hoganHight,mat;
    double ascale;

    mat = md.getMAT1();
    setVelocityDistribution(mat);

    ascale = 8000.0*mag;
    hoganWidth = 200; hoganHight = 200;

    // count
    gOff.setColor(Color.lightGray);
    for (i=0; i<200; i++) {
      iy = (int)(ascale/Nmt*vdf[0][i]);
      if (iy>0) {
        gOff.drawLine(gx+i,gy+200-iy,gx+i,gy+200);
      }
    }

    // Maxwell

      setMaxwell( md.getContTemp(), mat );
      gOff.setColor(Color.getHSBColor(0.01f,0.8f,0.6f));
      for (i=0; i<200-1; i++) {
        iy = (int)(ascale*ffvMaxwell[i]);
        iy2 = (int)(ascale*ffvMaxwell[i+1]);
        if (iy>0) {
          gOff.drawLine(gx+i,gy+200-iy,gx+i,gy+200-iy2);
        }
      }

      // MAT1
      gOff.setColor(Color.getHSBColor(md.getColorOf(mat),0.8f,0.6f));
      for (i=0; i<200-1; i++) {
        iy = (int)(ascale*meanvdf[1][i]);
        iy2 = (int)(ascale*meanvdf[1][i+1]);
        if (iy>0) {
          gOff.drawLine(gx+i,gy+200-iy,gx+i,gy+200-iy2);
        }
      }
    
      // MAT2
      //gOff.setColor(Color.getHSBColor(md.getColorOfKind(MAT2),0.8f,0.6f));
      //for (i=0; i<200-1; i++) {
      //  iy = (int)(ascale*meanvdf[2][i]);
      //  iy2 = (int)(ascale*meanvdf[2][i+1]);
      //  if (iy>0) {
      //    gOff.drawLine(gx+i,gy+200-iy,gx+i,gy+200-iy2);
      //  }
      //}


    // grid
    gOff.setColor(Color.gray);
    for (i=0; i<200; i+=50) {
      gOff.drawLine(gx+i,gy, gx+i, gy+200);
    }
    for (i=0; i<200; i+=100) {
      gOff.drawLine(gx,gy+200-i, gx+200, gy+200-i);
    }
    gOff.setColor(Color.black);
    gOff.drawRect(gx,gy,200,200);

    gOff.setColor(Color.black);
    gOff.drawString("0",gx-5,gy+hoganHight+20);
    gOff.drawString("2000 m/s",gx+hoganWidth-30,gy+hoganHight+20);
    gOff.drawString("velocity distribution",gx+30,gy+hoganHight+40);

    gOff.setColor(Color.lightGray);
    gOff.drawString("count",gx+120,gy+20);
    if (mat<21) {
      gOff.setColor(Color.getHSBColor(md.getColorOf(mat),0.8f,0.6f));
      gOff.drawString("mean distr.",gx+120,gy+40);
      //gOff.setColor(Color.getHSBColor(md.getColorOf(mat2),0.8f,0.6f));
      //gOff.drawString("mean distr.",gx+120,gy+60);
      gOff.setColor(Color.getHSBColor(0.01f,0.8f,0.6f));
      gOff.drawString("Maxwell(Tc)",gx+120,gy+80);
    }
  }

  private void setVelocityDistribution(int MAT1) {
    int i,iv,knd;
    double vv;

    for (i=0; i<220; i++) {
      vdf[0][i]=0; vdf[1][i]=0; vdf[2][i]=0;
    }

    for (i=0; i<Nmt; i++) {
      knd = md.getKind(i);
      vv = Math.sqrt(md.getv2(i));
      if (vv<2000.0) {
        iv = (int)(vv*0.1);
        vdf[0][iv] +=1;
        if (knd==MAT1) vdf[1][iv] +=1;
        //if (knd==MAT2) vdf[2][iv] +=1;
      }
    }

    for (i=0; i<220; i++) {
      meanvdf[0][i] = 0.99*meanvdf[0][i]+0.01*(vdf[0][i]/(Nmt*1.0));
      meanvdf[1][i] = 0.99*meanvdf[1][i]+0.01*(vdf[1][i]/(Nmt*1.0));
      meanvdf[2][i] = 0.99*meanvdf[2][i]+0.01*(vdf[2][i]/(Nmt*1.0));
    }

  }

  private void setMaxwell(double Tk, int mat) {
    int iv;
    double m,kT,a,v;

    m = md.getMassOf(mat);
    kT = 1.38e-23*Tk;
    a = m/kT*Math.sqrt(2.0*m/(3.14159*kT));

    for (iv=0; iv<200; iv++) {
      v = (double)(iv*10);
      ffvMaxwell[iv] = a*v*v*Math.exp(-0.5*m*v*v/kT)*10.0;
    }
  }


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

  public void setViewDelta(double xDisplacement, double yDisplacement) {
    theta += 0.5*3.14/180.0*xDisplacement;
    fai += 0.5*3.14/180.0*yDisplacement;
    cosTh = Math.cos(theta);
    sinTh = Math.sin(theta);
    cosFi = Math.cos(fai);
    sinFi = Math.sin(fai);
  }

  public void setViewScale(double s) { viewScale = s; }
  public double getViewScale() { return viewScale; }
}



// ====================== Lennard-Jones molecule class ===========

class LJMolecule {

  public static final double AU = 1.661e-27;
  public static final double EE = 1.602e-19;
  public static final double AA = 1.0e-10;
  public static final double KB = 1.38e-23;
  public static final double FS = 1.0e-15;

  public static final double Molec[][] = {
    //  0-mass   1-charge 2-E(in kB) 3-sigma(A) 4-dt(s) 5-color
    {  32.000*AU,  0.0*EE, 507.0*KB, 3.585*AA, 20.0*FS, 0.70 }, //0 CH3OH
    {  46.000*AU,  0.0*EE, 415.0*KB, 4.370*AA, 20.0*FS, 0.80 }, //1 C2H5OH
    {  60.000*AU,  0.0*EE, 577.0*KB, 4.549*AA, 20.0*FS, 0.90 }  //2 n-C3H7OH
  };

  private static final String matStr[] = {
    "CH3OH","C2H5OH","C3H7OH"
  };
  private static final int kindMax = 3;

  private double hh = 1.0e-12;
  public static double rCutoff = 1.0e-9;
  private double forceTable[][][] = new double[kindMax][kindMax][1020];


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

  LJMolecule() {
    setForceTable();
  }

  //

  private void setForceTable() {
    int i,j,ir; 
    double eps,sgm,r,ri,r3,r6;

    for (i=0; i<kindMax; i++) {
      for (j=0; j<kindMax; j++) {
        eps = Math.sqrt(Molec[i][2]*Molec[j][2]);
        sgm = 0.5*(Molec[i][3]+Molec[j][3]);
        for (ir=10; ir<1000; ir++) {
          r = (double)ir*hh;
          ri = sgm/r;
          r3 = ri*ri*ri;
          r6 = r3*r3;
          forceTable[i][j][ir] = cutoff(r)*(24.0*eps*r6*(2.0*r6-1.0)/r);
        }
        for (ir=0; ir<10; ir++) {
          forceTable[i][j][ir]=forceTable[i][j][10];
        }
      }
    }
  }

  private double cutoff(double r) {
    double rt;

    rt = 0.0;
    if (r<0.8*rCutoff) {
      rt = 1.0;
    } else if (r>=0.8*rCutoff && r<rCutoff) {
      rt = 0.5+0.5*Math.cos(3.1415926536*(r-0.8*rCutoff)/(0.2*rCutoff));
    } else {
      rt = 0.0;
    }
    return( rt );
  }


  // ----------------------------------- class methods -----------

  public double force(double r, int ki, int kj) {
    int i;
    double a;

    i = (int)(r/hh);
    a = r - i*hh;
    return ( ((hh-a)*forceTable[ki][kj][i] + a*forceTable[ki][kj][i+1])/hh );
  }


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

  public int getKindMax() { return kindMax; }
  public String getMatStr(int knd) { return matStr[knd]; }

  public double getMassOf(int knd) { return Molec[knd][0]; }
  public double getDiaOf(int knd) { return 1.12246*Molec[knd][3]; }
  public double getColorOf(int knd) { return Molec[knd][5]; }
  public double getrCutoff() { return rCutoff; }
}



// ================== Langevin dynamics - LJ potential ===========

class LD3D implements Runnable {

  private LJMolecule ma;

  private double t = 0.0;
  public static double dt = 5.0*1.0e-15;
  public static double xMax = 20.0*1.0e-9;
  public static double yMax = 20.0*1.0e-9;
  public static double zMax = 20.0*1.0e-9;

  public static int Nmt = 16*16*8;

  private int Nsx = 20;
  private int Nsy = 20;
  private int Nsz = 20;

  private int MAT1 = 1; // C2H5OH
  private int MAT2 = 0; // CH3OH
  private int mat2content = 50; // 0-100

  private double rc = ma.rCutoff;
  private double rc2 = rc*rc;

  private int tempMode = 0;
  private double contTemp = 300.0;

  private double gamma = 1.0e-13;
  private double kBoltzmann = 1.38e-23;
  private double mH2O = 18*1.67e-27;
  private double alpha = 0.183*Math.sqrt(3.0*mH2O*kBoltzmann)/dt;

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

  private int reg[][] = new int[NN][300];
  private int section[][][][] = new int[Nsx][Nsy][Nsz][100];

  private double virial[][] = new double[3][256];
  private double memTemp[] = new double[256];
  private int virialp = 0;

  private int resetSW = 0;
  private int startSW = 1;
  private int stepSW = 0;
  private int sleepTime = 1;
  private int stopSleepTime = 100;
  private int count = 0;


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

  LD3D(LJMolecule ma3d) {
    ma = ma3d;
    setInitialCondition();
  }


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

  public void run() {
    while (true) {
      count += 1;
      timeEvolution();
      if (sleepTime>0) {
        try {
          Thread.sleep(sleepTime);
        } catch (InterruptedException e) { }
      }
    }
  }


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

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

    count = 0;
    t = 0.0;

    s = xMax/16.0;
    for (i=0; i<Nmt; i++) {
      kind[i] = MAT1;
      xx[i] = s*(double)(i/256)+0.5*s;
      yy[i] = s*(double)((i%256)/16)+0.5*s;
      zz[i] = s*(double)((i%256)%16)+0.5*s;
      mr = Math.sqrt(ma.AU/ma.getMassOf(kind[i]));
      vx[i] = 800.0*mr*(Math.random()-0.5);
      vy[i] = 800.0*mr*(Math.random()-0.5);
      vz[i] = 800.0*mr*(Math.random()-0.5);
      ffx[i] = 0.0;
      ffy[i] = 0.0;
      ffz[i] = 0.0;
    }

    vAjustment();
    initVirial();
  }

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

  private void timeEvolution(){

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

    if (startSW==1) {
      timeStep();

    } else {
      if (stepSW==1) {
        timeStep();
        stepSW = 0;
      }
      try {
        Thread.sleep(stopSleepTime);
      } catch (InterruptedException e) { }
    }
  }

  private void timeStep() {
    int i;

    if (tempMode==1) {
      vAjustment();
    }
    registration();
    for (i=0; i<20; i++) {
      time1Step();
    }
    setVirial();
  }


  private void time1Step() {
    int i; 
    double dtv2, a2, rr, tmp;

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

    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++) {
      a2 = dtv2/ma.getMassOf(kind[i]);
      vx[i] += a2*ffx[i];
      vy[i] += a2*ffy[i];
      vz[i] += a2*ffz[i];
      xx[i] += vx[i]*dt;
      yy[i] += vy[i]*dt;
      zz[i] += vz[i]*dt;
    }
    forceCalc();
    for (i=0; i<Nmt; i++) {
      a2 = dtv2/ma.getMassOf(kind[i]);
      vx[i] += a2*ffx[i];
      vy[i] += a2*ffy[i];
      vz[i] += a2*ffz[i];
    }

    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];
      }
    }
  }

  private void forceCalc() {
    int i,j,k;
    double ffr,frz,the,a;
    double r2,rij,f,fxij,fyij,fzij;

    a = alpha*Math.sqrt(contTemp);
    for(i=0;i<Nmt;i++) {
      ffr = a*gaussianRandom();
      frz = 2.0*Math.random() - 1.0;
      the = 2.0*Math.PI*Math.random();
      ffx[i] = -gamma*vx[i] + ffr*Math.sqrt(1.0-frz*frz)*Math.cos(the);
      ffy[i] = -gamma*vy[i] + ffr*Math.sqrt(1.0-frz*frz)*Math.sin(the);
      ffz[i] = -gamma*vz[i] + ffr*frz;
    }

    for(i=0;i<Nmt;i++) {
      for(k=1;k<reg[i][0];k++) {
        j = reg[i][k];
        r2=(xx[i]-xx[j])*(xx[i]-xx[j])+(yy[i]-yy[j])*(yy[i]-yy[j])+(zz[i]-zz[j])*(zz[i]-zz[j]);
        if (r2<rc2) {
          rij = Math.sqrt(r2);
          f = ma.force(rij,kind[i],kind[j]);
          fxij = f*(xx[i]-xx[j])/rij;
          fyij = f*(yy[i]-yy[j])/rij;
          fzij = f*(zz[i]-zz[j])/rij;

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

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

  //registration

  private void registration() {
    int i,j,k,ip,jp,kp,i0,i1,j0,j1,k0,k1,iq;
    double r2,rreg,rreg2;
    double x,y,z;

    preRegistration();

    rreg = rc + 3000.0*20*dt;
    rreg2 = rreg*rreg;

    for(ip=0;ip<Nmt;ip++) {
      kp = 1;

      i0 = (int)(Nsx*(xx[ip]-rreg)/xMax); if (i0<0) i0 = 0;
      i1 = (int)(Nsx*(xx[ip]+rreg)/xMax); if (i1>=Nsx) i1 = Nsx-1;
      j0 = (int)(Nsy*(yy[ip]-rreg)/yMax); if (j0<0) j0 = 0;
      j1 = (int)(Nsy*(yy[ip]+rreg)/yMax); if (j1>=Nsy) j1 = Nsy-1;
      k0 = (int)(Nsz*(zz[ip]-rreg)/zMax); if (k0<0) k0 = 0;
      k1 = (int)(Nsz*(zz[ip]+rreg)/zMax); if (k1>=Nsz) k1 = Nsz-1;

      for(i=i0;i<=i1;i++) {
        for(j=j0;j<=j1;j++) {
          for(k=k0;k<=k1;k++) {
            for(iq=1;iq<=section[i][j][k][0];iq++) {
              jp = section[i][j][k][iq];
              if (jp>ip) {
                r2=(xx[ip]-xx[jp])*(xx[ip]-xx[jp])+(yy[ip]-yy[jp])*(yy[ip]-yy[jp])+(zz[ip]-zz[jp])*(zz[ip]-zz[jp]);
                if (r2<rreg2) {
                  reg[ip][kp]=jp;
                  kp += 1;
                }

              }
            }
          }
        }
      }
      reg[ip][0]=kp;

    }
  }

  private void preRegistration() {
    int i,j,k,ip,iq;

    for(i=0;i<Nsx;i++) {
      for(j=0;j<Nsy;j++) {
        for(k=0;k<Nsz;k++) {
          section[i][j][k][0] = 0;
        }
      }
    }

    for(ip=0;ip<Nmt;ip++) {
      i = (int)(Nsx*xx[ip]/xMax); if (i>=Nsx) i = Nsx-1;
      j = (int)(Nsy*yy[ip]/yMax); if (j>=Nsy) j = Nsy-1;
      k = (int)(Nsz*zz[ip]/zMax); if (k>=Nsz) k = Nsz-1;
      iq = section[i][j][k][0]+1;
      section[i][j][k][0] = iq;
      section[i][j][k][iq] = ip;
    }
  }


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

  private double gaussianRandom() { // Box-Muller method
    double rnd;
    rnd = Math.sqrt( -2.0*Math.log(Math.random()) );
    return rnd*Math.sin(2.0*Math.PI*Math.random());
  }

  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];
      vz[i] = r*vz[i];
    }
  }

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

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

  private double temp(int i) {
    double ke;

    ke = 0.5*ma.getMassOf(kind[i])*(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    return ( ke/(1.5*1.38e-23) );
  }



  public int regMax() {
    int i,m;

    m = 0;
    for(i=0;i<Nmt;i++) {
      if (reg[i][0]>m) {
        m = reg[i][0];
      }
    }
    return(m);
  }

  public int secMax() {
    int i,j,k,m;

    m = 0;
    for(i=0;i<Nsx;i++) {
      for(j=0;j<Nsy;j++) {
        for(k=0;k<Nsz;k++) {
          if (section[i][j][k][0]>m) {
            m = section[i][j][k][0];
          }
        }
      }
    }
    return(m);
  }

  // pressure

  private void initVirial() {
    int i;

    for (i=0; i<256; i++) {
      virial[0][i] = 0.0;
      virial[1][i] = 0.0;
      virial[2][i] = 0.0;
      memTemp[i] = 0.0;
    }
    virialp = 0;
  }

  private void setVirial() {
    int i;
    double xVir,yVir,zVir;

    xVir = 0.0; yVir = 0.0; zVir = 0.0;
    for (i=0; i<Nmt; i++) {
      xVir += ffx[i]*xx[i];
      yVir += ffy[i]*yy[i];
      zVir += ffz[i]*zz[i];
    }
    virial[0][virialp] = xVir;
    virial[1][virialp] = yVir;
    virial[2][virialp] = zVir;
    memTemp[virialp] = temperature();
    virialp = (virialp+1)%256;
  }

  private double meanTemp() {
    int i;
    double st,tm;

    st = 0;
    for (i=0; i<256; i++) {
      st += memTemp[i];
    }
    tm = st/256.0;
    return( tm );
  }

  public double meanPressure() {
    int i;
    double s,vir,tm,vol,p;

    s = 0.0;
    for (i=0; i<256; i++) {
      s += virial[0][i]+virial[1][i]+virial[2][i];
    }
    vir = s/256.0;
    tm = meanTemp();
    vol = xMax*yMax*zMax;
    p = (Nmt*1.38e-23*tm+vir/3.0)/vol;
    return( p );
  }

  public double pressure() {
    int i,im;
    double vir,tm,vol,p;

    im = (virialp-1+256)%256;
    vir = virial[0][im]+virial[1][im]+virial[2][im];
    tm = temperature();
    vol = xMax*yMax*zMax;
    p = (Nmt*1.38e-23*tm+vir/3.0)/vol;
    return( p );
  }


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

  // control
  public void reset() { resetSW = 1; }
  public void setStartSW(int sw) { startSW = sw; }
  public int getStartSW() { return startSW; }
  public void setStepSW() { stepSW = 1; }
  public int getCount() { return count; }

  public void setTempMode(int tmode) { tempMode = tmode; }
  public int getTempMode() { return tempMode; }
  public void setContTemp(double tmp) { contTemp = tmp; vAjustment(); }
  public double getContTemp() { return contTemp; }

  public void setMAT1(int m) { MAT1 = m; resetSW = 1; }
  public int getMAT1() { return MAT1; }
  public void setMAT2(int m) { MAT2 = m; resetSW = 1; }
  public int getMAT2() { return MAT2; }
  public void setMat2content(int m2c) { mat2content = m2c; resetSW = 1; }
  public int getMat2content() { return mat2content; }

  // get data
  public double getTime() { return t; }
  public int getNmt() { return Nmt; }
  public float getColorOf(int knd) { return (float)ma.getColorOf(knd); }
  public float getMassOf(int knd) { return (float)ma.getMassOf(knd); }

  public double getDia(int i) { return ma.getDiaOf(kind[i]); }
  public double getMass(int i) { return ma.getMassOf(kind[i]); }
  public float getColor(int i) { return (float)ma.getColorOf(kind[i]); }

  public int getKind(int i) { return kind[i]; }
  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 getv2(int i) { return vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i]; }
  public double getTemp(int i) { return temp(i); }

  public int getnnReg(int i) { return reg[i][0]; }
  public int getReg(int i, int k) { return reg[i][k]; }
}