/** applet No. 1253
 *
 * template - Langevin dynamics 2D
 * - Lennard-Jones potential
 * - multi thread  LD2D :: display - asynchronous
 *
 * Created by Ikeuchi Mitsuru on January 26 2008.
 * Copyright (c) 2008 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2008.01.26  created
 *
 */

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

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

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

  Thread th = null; // for run()-paint() thread
  Thread thmd = null; // for LD2D() thread
  LJMolecule ma = new LJMolecule();
  LD2D md = new LD2D(ma);
  DrawGraph2D dg2d = new DrawGraph2D();

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

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

  private int dgX, dgY, dgXb, dgYb; // mouse

  private int sleepTime = 50;
  private double viewScale = 1.0;
  private int dispMode = 0;
  private 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("temp");
    ch_view.add("velocity");
    ch_view.add("vSpace");
    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); 

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

    sc_content= new Scrollbar(Scrollbar.HORIZONTAL,50,10,0,110);
    sc_content.addAdjustmentListener(this);

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

    sc_scale= new Scrollbar(Scrollbar.HORIZONTAL,100,10,50,510);
    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(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() );
    } else if (ev.getSource() == ch_mat2) {
      md.setMAT2( ch_mat2.getSelectedIndex() );
    } 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_content) { 
      md.setMat1content( sc_content.getValue() );
    } else if (ev.getSource() == sc_temp) { 
      md.setContTemp( 1.0*(double)(sc_temp.getValue()) );
    } else if (ev.getSource() == sc_scale) { 
      viewScale = 0.01*(double)(sc_scale.getValue());
      dg2d.setScale(viewScale);
    }
  }

  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;
    }
    dg2d.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) {
      dg2d.plotGraph2D(gOff, md, dispMode, thCount);
    }
  }

}



// =============================== draw graph 2D class ===========

class DrawGraph2D {

  private Graphics gOff;
  private LD2D md;

  private int Nmt = LD2D.Nmt;

  private double xMax = LD2D.xMax;
  private double yMax = LD2D.yMax;

  private double dispScale = 300.0/xMax;
  private double viewScale = 1.0;

  private int gx0 = 30;
  private int gy0 = 110;
  private int xImageLoc = gx0;
  private int yImageLoc = gy0;

  private int dragSW = 1;
  private double xShift = 0.0;
  private double yShift = 0.0;
  private double viewTheta = -15.0*3.14/180.0;
  private double viewFai = -72.0*3.14/180.0;
  private double cosTh = Math.cos(viewTheta); 
  private double sinTh = Math.sin(viewTheta);
  private double cosFi = Math.cos(viewFai);
  private double sinFi = Math.sin(viewFai);

  private int rdf[][] = new int[120][3];
  private double meanRDF[][] = new double [120][3];
  private int vdf[][] = new int [200][3];
  private double meanVDF[][] = new double [200][3];

  private int drawnCount = 0;


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

  DrawGraph2D() {
    
  }


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

  public void plotGraph2D(Graphics gg, LD2D md2d, int dispMode, int thCount) {
    int gbh,mdCount;

    gOff = gg;
    md = md2d;

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

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

    if (dispMode==0) { // ball
      drawBox();
      ballPlot(0.6, 0.0);
      rdfPlot(400, 180, 1.0);

    } else if (dispMode==1) { // line
      drawBox();
      linePlot();
      rdfPlot(400, 180, 1.0);

    } else if (dispMode==2) { // ball + line
      drawBox();
      ballPlot(0.6, 0.0);
      linePlot();
      rdfPlot(400, 180, 1.0);

    } else if (dispMode==3) { // temp
      drawBox();
      ballPlot(0.6, 2000.0);
      rdfPlot(400, 180, 1.0);

    } else if (dispMode==4) { // velocity
      drawBox();
      velocityPlot(500.0);
      vdfPlot(400, 180, 2.0);

    } else if (dispMode==5) { // v-space
      vSpacePlot(320, 3000.0);
      vdfPlot(400, 180, 2.0);

    } else if (dispMode==6) {
      ;

    }

    gbh = 630/6;

    gOff.setColor(Color.black);
    gOff.drawString("Tcont="+md.getContTemp()+"K",gbh*3+10,70);
    gOff.drawString("scale="+(int)(viewScale*100.0+0.499)+"%",gbh*4+10,70);
    gOff.drawString("view",gbh*5+10,70);

    gOff.setColor(Color.blue);
    gOff.drawString("t="+(int)(md.getTime()*1.0e14+0.5)/100.0+" ps",gbh*0+10,90);
    gOff.drawString("T="+(int)(md.getTemperature()*10+0.5)/10.0+" K",gbh*3+10,90);
    gOff.drawString("Lx="+(int)(md.xMax*1.0e10+0.5)/10.0+" nm",gbh*4+10,90);

    gOff.setColor(Color.blue);
    gOff.drawString("regMax="+md.getRegMax()+"",540,90);

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


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

  private void ballPlot(double dia, double tempSpan) {
    int i,ix,iy,ir;
    double sc,tm;

    sc = viewScale*dispScale;
    for (i=0; i<Nmt; i++) {
      ix = xImageLoc + (int)(sc*md.getxx(i));
      iy = yImageLoc + (int)(sc*md.getyy(i));
      ir = (int)(sc*dia*md.getDia(i));
      tm = 0.95;
      if (tempSpan>0.0) {
        tm = md.getTemp(i);
        tm = 0.8*tm/tempSpan+0.199;
        if (tm>0.99) tm = 0.99;
      }
      gOff.setColor(Color.getHSBColor((float)(md.getColor(i)),0.95f,(float)tm));
      gOff.fillOval(ix-ir/2,iy-ir/2, ir, ir);
    }
  }

  private void linePlot() {
    int i,j,k,ix,iy,ix2,iy2;
    double sc,r2,dij;

    sc = viewScale*dispScale;
    for (i=0; i<Nmt; i++) {
      ix = xImageLoc + (int)(sc*md.getxx(i));
      iy = yImageLoc + (int)(sc*md.getyy(i));

      for(k=1;k<md.getnnReg(i);k++) {
        j = md.getReg(i, k);
        r2 = (md.getxx(i)-md.getxx(j))*(md.getxx(i)-md.getxx(j))+(md.getyy(i)-md.getyy(j))*(md.getyy(i)-md.getyy(j));
	dij = 1.1*(md.getDia(i)+md.getDia(j))/2.0;
        if (r2<dij*dij && j>i) {
          ix2 = xImageLoc + (int)(sc*md.getxx(j));
          iy2 = yImageLoc + (int)(sc*md.getyy(j));

          gOff.setColor(Color.gray);
          gOff.drawLine(ix,iy, ix2,iy2);
        }
      }
    }
  }

  private void velocityPlot(double vSpan) {
    int i,ix,iy,ix2,iy2;
    double sc;

    sc = viewScale*dispScale;
    for (i=0; i<Nmt; i++) {
      ix = xImageLoc + (int)(sc*md.getxx(i));
      iy = yImageLoc + (int)(sc*md.getyy(i));

      ix2 = ix+(int)(viewScale*10.0*md.getvx(i)/vSpan);
      iy2 = iy+(int)(viewScale*10.0*md.getvy(i)/vSpan);
      gOff.setColor(Color.getHSBColor((float)(md.getColor(i)),0.8f,0.8f));
      gOff.drawLine(ix,iy, ix2,iy2);
    }
  }

  //

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

    sc = viewScale*dispScale;
    gball = (int)(3.0*viewScale);

    gOff.setColor(Color.gray); 
    if (md.getTempMode()==2) {
      gOff.setColor(Color.black); 
      if (md.getTemperature()<md.getContTemp()) { 
        gOff.setColor(Color.red);
      }
    }
    ix = (int)(sc*xMax);
    iy = (int)(sc*yMax);
    gOff.drawRect(xImageLoc-gball,yImageLoc-gball,ix+2*gball,iy+2*gball);
  }

  //

  private void vSpacePlot(int plotWidth, double scale) {
    int i,ix,iy,ir,ix0,iy0;

    ix0 = xImageLoc + plotWidth/2;
    iy0 = yImageLoc + plotWidth/2;
    gOff.setColor(Color.black);
    gOff.drawLine(ix0, yImageLoc, ix0, yImageLoc+plotWidth);
    gOff.drawLine(xImageLoc, iy0, xImageLoc+plotWidth, iy0);

    for (i=0; i<Nmt; i++) {
      ix = (int)(ix0 + 0.5*plotWidth*md.getvx(i)/scale );
      iy = (int)(iy0 + 0.5*plotWidth*md.getvy(i)/scale );
      ir = (int)(0.6*md.getDia(i)/1.0e-10+0.5);

      gOff.setColor(Color.getHSBColor(md.getColor(i),0.4f,0.9f));
      gOff.fillOval(ix-ir/2,iy-ir/2, ir, ir);

    }
  }


  // redial distribution

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

    setRadialDistribution();

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

    gOff.setColor(Color.lightGray);
    for (i=10; i<100; i++) {
      n1 = (int)(mag*meanRDF[i][0]*(1600.0/(i*i)));
      if (n1>0) {
        gOff.fillRect(gx+i*2,gy+hoganHight-n1,2,n1);
      }
    }

    gOff.setColor(Color.getHSBColor(md.getColorOf(md.getMAT1()),0.9f,0.6f));
    for (i=10; i<100-1; i++) {
      n1 = (int)(mag*meanRDF[i][1]*(1600.0/(i*i)));
      n2 = (int)(mag*meanRDF[i+1][1]*(1600.0/(i*i)));
      gOff.drawLine(gx+i*2,gy+hoganHight-n1,gx+i*2+2,gy+hoganHight-n2);
    }

    gOff.setColor(Color.getHSBColor(md.getColorOf(md.getMAT2()),0.9f,0.6f));
    for (i=10; i<100-1; i++) {
      n1 = (int)(mag*meanRDF[i][2]*(1600.0/(i*i)));
      n2 = (int)(mag*meanRDF[i+1][2]*(1600.0/(i*i)));
      gOff.drawLine(gx+i*2,gy+hoganHight-n1,gx+i*2+2,gy+hoganHight-n2);
    }

    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 (n/r^2)",gx+20,gy+hoganHight+40);
  }

  private void setRadialDistribution() {
    int i,j,k,ir,m1,m2,ki,kj;
    double xij,yij,rij;

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

    m1 = md.getMAT1(); m2 = md.getMAT2();
    for (i=0; i<Nmt; i++) {
      for(k=1;k<md.getnnReg(i);k++) {
        j = md.getReg(i, k);
        xij = md.getxx(i)-md.getxx(j); yij = md.getyy(i)-md.getyy(j);
        rij = Math.sqrt(xij*xij+yij*yij);
        if (rij<10.0e-10) {
          ki = md.getKind(i); kj = md.getKind(j);
          ir = (int)(rij*1.0e11+0.5); if (ir>119) ir=119;
          rdf[ir][0] +=1;
          if (ki==m1 && kj==m1) {
            rdf[ir][1] += 1;
          } else if (ki==m2 && kj==m2) {
            rdf[ir][2] += 1;
          }
        }
      }
    }

    for (i=0; i<120; i++) {
      meanRDF[i][0] = 0.9*meanRDF[i][0]+0.1*rdf[i][0];
      meanRDF[i][1] = 0.9*meanRDF[i][1]+0.1*rdf[i][1];
      meanRDF[i][2] = 0.9*meanRDF[i][2]+0.1*rdf[i][2];
    }
  }

  // velocity distribution

  private void vdfPlot(int gx, int gy, double mag) {
    int i,n1,n2;
    double yMag;

    setVelocityDistribution();

    yMag = mag*0.005*Nmt;
    gOff.setColor(Color.lightGray);
    for (i=0; i<200-1; i++) {
      n1  = (int)(yMag*meanVDF[i][0]);
      gOff.drawLine(gx+i,gy+200-n1,gx+i,gy+200);
    }

    gOff.setColor(Color.getHSBColor(md.getColorOf(md.getMAT1()),0.9f,0.6f));
    for (i=0; i<200-2; i++) {
      n1 = (int)(yMag*meanVDF[i][1]);
      n2 = (int)(yMag*meanVDF[i+1][1]);
      gOff.drawLine(gx+i,gy+200-n1,gx+i+1,gy+200-n2);
    }

    gOff.setColor(Color.getHSBColor(md.getColorOf(md.getMAT2()),0.9f,0.6f));
    for (i=0; i<200-2; i++) {
      n1 = (int)(yMag*meanVDF[i][2]);
      n2 = (int)(yMag*meanVDF[i+1][2]);
      gOff.drawLine(gx+i,gy+200-n1,gx+i+1,gy+200-n2);
    }

    drawHogan(gx, gy, 200, 200, 50, 100);
    gOff.setColor(Color.black);
    gOff.drawString("0",gx-5,gy+220);
    gOff.drawString("1000",gx+100-15,gy+220);
    gOff.drawString("(m/s)",gx+200-30,gy+220);
    gOff.drawString("velocity distribution",gx+20,gy+235);
  }

  private void setVelocityDistribution() {
    int i,iv,m1,m2,knd;
    double vx,vy,v;

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

    m1 = md.getMAT1(); m2 = md.getMAT2();
    for(i=0;i<Nmt;i++) {
      knd = md.getKind(i);
      vx = md.getvx(i); vy = md.getvy(i);
      v = Math.sqrt(vx*vx+vy*vy);
      iv = (int)(v/10.0); if (iv>199) iv = 199;
      vdf[iv][0] += 1;
      if (knd==m1) {
        vdf[iv][1] += 1;
      } else if (knd==m2) {
        vdf[iv][2] += 1;
      }
    }

    for(i=0;i<200;i++) {
      meanVDF[i][0] = 0.98*meanVDF[i][0]+0.02*vdf[i][0];
      meanVDF[i][1] = 0.98*meanVDF[i][1]+0.02*vdf[i][1];
      meanVDF[i][2] = 0.98*meanVDF[i][2]+0.02*vdf[i][2];
    }
  }

  //

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

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

  public void setViewDelta(double xDisplacement, double yDisplacement) {

    if (dragSW==1) {
      xShift += (1.0/viewScale/dispScale)*xDisplacement;
      yShift += (1.0/viewScale/dispScale)*yDisplacement;
      xImageLoc = gx0+(int)(dispScale*viewScale*xShift);
      yImageLoc = gy0+(int)(dispScale*viewScale*yShift);
    } else if (dragSW==2) {
      viewTheta += 0.5*3.14/180.0*xDisplacement;
      viewFai += 0.5*3.14/180.0*yDisplacement;
      cosTh = Math.cos(viewTheta); 
      sinTh = Math.sin(viewTheta);
      cosFi = Math.cos(viewFai);
      sinFi = Math.sin(viewFai);
    }
  }
  public void setDragSW(int sw) { dragSW = sw; }
  public void setScale(double sc) { viewScale = sc; }

}



// ====================== 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;
  private 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 class ===========

class LD2D implements Runnable {

  private LJMolecule ma;

  private double t = 0.0;
  public static double dt = 5.0*1.0e-15;
  public static double xMax = 25.0*1.0e-9;
  public static double yMax = 25.0*1.0e-9;
  public static int Nmt = 40*20;

  private int Nsx = 50;
  private int Nsy = 50;

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

  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(2.0*mH2O*kBoltzmann)/dt;

  public final int NN = 10000;
  private int kind[] = new int[NN];
  private double xx[] = new double[NN];
  private double yy[] = new double[NN];
  private double vx[] = new double[NN];
  private double vy[] = new double[NN];
  private double ffx[] = new double[NN];
  private double ffy[] = new double[NN];

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

  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 -----------

  LD2D(LJMolecule ma2d) {
    ma = ma2d;
    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;

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

    vAjustment();
  }


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

  private void timeEvolution() {

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

    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<10; i++) {
      time1Step();
    }
  }

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

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

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

    setBoundary();
  }

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

    rc2 = ma.getrCutoff()*ma.getrCutoff();

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

    for(i=0;i<Nmt;i++) {
      for(k=1;k<reg[i][0];k++) {
        j = reg[i][k];
        xij = xx[i] - xx[j]; yij = yy[i] - yy[j];
        r2=xij*xij+yij*yij;
        if (r2<rc2) {
          rij = Math.sqrt(r2);
          f = ma.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 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];
      }
    }
  }

  //registration

  private void registration() {
    int i,j,ip,jp,kp,i0,i1,j0,j1,iq;
    double r2,rreg,rreg2;

    preRegistration();

    rreg = ma.getrCutoff() + 3000.0*10*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;

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

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

    for(i=0;i<Nsx;i++) {
      for(j=0;j<Nsy;j++) {
        section[i][j][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;
      iq = section[i][j][0]+1;
      section[i][j][0] = iq;
      section[i][j][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];
    }
  }

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

    s = 0.0; 
    for (i=0; i<Nmt; i++) {
      s += ma.getMassOf(kind[i])*(vx[i]*vx[i] + vy[i]*vy[i]);
    }
    return (0.5*s/((double)Nmt*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]);
    return (ke/1.38e-23);
  }

  private int regMax() {
    int i,m;

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


  // ------------------------------------- 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 setMat1content(int m2c) { mat2content = m2c; resetSW = 1; }
  public int getMat2content() { return mat2content; }

  // get data
  public double getTime() { return t; }
  public double getdt() { return dt; }
  public int getNmt() { return Nmt; }
  public float getColorOf(int knd) { return (float)ma.getColorOf(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 getvx(int i) { return vx[i]; }
  public double getvy(int i) { return vy[i]; }
  public double getTemp(int i) { return temp(i); }

  public double getTemperature() { return temperature(); }

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