/** applet No. 1014
 *
 * sectiom profile - molecular dynamics 3D
 *
 * Created by Ikeuchi Mitsuru on September 09 2006.
 * Copyright (c) 2006-2007 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2006.09.09  created
 * ver 0.0.2   2006.09.10  view (added section thickness) improved
 * ver 0.0.3   2006.09.13  added NN selectiom (8^3 - 16^3)
 * ver 0.0.4   2006.09.15  bug fixed
 * ver 0.0.5   2006.09.16  added neighbor statistics
 * ver 0.0.6   2006.11.26  added N=343 selection
 * ver 0.0.7   2007.06.08  improved code
 *
 */

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

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

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

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

  // for event
  Choice ch_mat,ch_nCube,ch_tempControl,ch_view;
  Button bt_reset, bt_startStop, bt_step;
  Scrollbar sc_temp,sc_scale,sc_zSection,sc_zThick;

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

  int sleepTime = 30;

  double t = 0.0; // time (s)
  double dt = 10.0*1.0e-15; // time division (s)
  double xMax = 10.0*1.0e-9; // x-Box size in (m)
  double yMax = 10.0*1.0e-9; // y-Box size in (m)
  double zMax = 10.0*1.0e-9; // z-Box size in (m)
  double zDepth = 10.0*1.0e-9;

  int nCube = 8;
  int Nmt = nCube*nCube*nCube;

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

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

  int resetSW = 0;
  int started = 1;
  int stepSW = 0;
  int dispMode = 3;
  int tempMode = 2;
  double contTemp = 300.0;
  double zSection = 1.0;
  double zThick = 1.0;

  int dgX, dgY, dgXb,dgYb;
  double cx = 0.5*xMax;
  double cy = 0.5*yMax;
  double cz = 0.5*zMax;
  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;

  double rc = 1.0e-9;
  double rc2 = rc*rc;

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

  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;

  int MAT1 = 17;
  int NN = 5000;
  int kind[] = new int[NN]; // i-th molec kind
  double xx[] = new double[NN]; // i-th x-position
  double yy[] = new double[NN]; // i-th y-position
  double zz[] = new double[NN]; // i-th z-position
  double vx[] = new double[NN]; // i-th x-velocity
  double vy[] = new double[NN]; // i-th y-velocity
  double vz[] = new double[NN]; // i-th z-velocity
  double ffx[] = new double[NN]; // i-th x-force
  double ffy[] = new double[NN]; // i-th y-force
  double ffz[] = new double[NN]; // i-th z-force

  int reg[][] = new int[NN][750]; // registration
  int section[][][][] = new int[Nsx][Nsy][Nsz][300];

  double hh = 1.0e-12;
  double forceTable[][][] = new double[NmatKind][NmatKind][1020];

  int srtz[][] = new int[2][NN]; // z-sort work space

  double dispColor[] = new double[NmatKind];

  int neighbor[] = new int[NN];
  double neighborColor[] = {
    0.65, 0.60, 0.55, 0.50, 0.45,
    0.40, 0.35, 0.30, 0.25, 0.20,
    0.15, 0.10, 0.05, 0.05, 0.05,
    0.01, 0.01, 0.01, 0.01, 0.01
  };
  int neighborDistance[][] = new int[16][30];

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

  int rdf[] = new int[120];
  int vdf[] = new int[220];

  double virial[][] = new double[3][256];
  double memTemp[] = new double[256];
  int virialp = 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_mat = new Choice();
/*
    ch_mat.add("W");  ch_mat.add("Mo"); ch_mat.add("Cr");
    ch_mat.add("Fe"); ch_mat.add("Ni"); ch_mat.add("Al");
    ch_mat.add("Pb"); ch_mat.add("Cu"); ch_mat.add("Ag"); 
    ch_mat.add("Ca"); ch_mat.add("Sr"); ch_mat.add("Ba");
    ch_mat.add("Na"); ch_mat.add("K");  ch_mat.add("Rb");
    ch_mat.add("Cs"); 
*/
    ch_mat.add("Ne"); 
    ch_mat.add("Ar");
    ch_mat.add("Kr");
    ch_mat.add("Xe");
    ch_mat.add("Hg");
    ch_mat.addItemListener(this);
    ch_mat.select("Ar");

    ch_nCube = new Choice();
    ch_nCube.add("N=343");
    ch_nCube.add("N=512");
    ch_nCube.add("N=729"); 
    ch_nCube.add("N=1000");
    ch_nCube.add("N=1331");
    ch_nCube.add("N=1728"); 
    ch_nCube.add("N=2197");
    ch_nCube.add("N=2744");
    ch_nCube.add("N=3375"); 
    ch_nCube.add("N=4096");
    ch_nCube.addItemListener(this);
    ch_nCube.select("N=512"); 

    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("momentum");
    ch_view.add("line+rdf");
    ch_view.add("v+vdf");
    ch_view.add("neighbor");
    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,10,610);
    sc_temp.addAdjustmentListener(this);

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

    sc_zSection = new Scrollbar(Scrollbar.HORIZONTAL,100,10,10,110);
    sc_zSection.addAdjustmentListener(this);

    sc_zThick = new Scrollbar(Scrollbar.HORIZONTAL,100,10,10,110);
    sc_zThick.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(ch_nCube);
    pnl.add(bt_reset);
    pnl.add(bt_startStop);
    pnl.add(bt_step);
    pnl.add(new Label(" section "));
    pnl.add(ch_view);

    pnl.add(ch_mat);
    pnl.add(ch_tempControl);
    pnl.add(sc_temp);
    pnl.add(sc_scale);
    pnl.add(sc_zSection);
    pnl.add(sc_zThick);
    add(pnl,"North");

    setForceTable();
    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_mat){
      MAT1 = 16 + ch_mat.getSelectedIndex();
	  setInitialCondition();
    } else if (ev.getSource() == ch_nCube) {
      nCube = 7 + ch_nCube.getSelectedIndex();
	  Nmt = nCube*nCube*nCube;
	  setInitialCondition();
    } else if (ev.getSource() == ch_tempControl) {
      tempMode = ch_tempControl.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 (started==0) {
        started = 1;
      } else {
        started = 0;
      }
    } else if (ev.getSource() == bt_step){
      if (started==0) stepSW = 1;
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_temp) { 
      contTemp = 1.0*(double)(sc_temp.getValue());
      vAjustment();
    } else if(ev.getSource() == sc_scale) { 
      viewScale = 0.01*(double)(sc_scale.getValue());
    } else if(ev.getSource() == sc_zSection) { 
      zSection = 0.01*(double)(sc_zSection.getValue());
    } else if(ev.getSource() == sc_zThick) { 
      zThick = 0.01*(double)(sc_zThick.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;
    }
    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) {
      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 -----------

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

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

    rotate(theta,fai);
    zSort();
    tmp = temperature();
    gbx = 70;
    gby = 140;


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

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

    } else if (dispMode==2) {
      ballPlot(gbx,gby, 0.5, 2);

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

    } else if (dispMode==4) {
      ballPlot(gbx,gby, 1.0, 3);

    } else if (dispMode==5) {
      velocity3DPlot(gbx,gby, 10.0,0);

    } else if (dispMode==6) {
      velocity3DPlot(gbx,gby, 10.0,1);

    } else if (dispMode==7) {
      ballPlot(gbx,gby, 0.5, 1);
      rdfPlot(400,140, 1.0*512.0/Nmt);

    } else if (dispMode==8) {
      velocity3DPlot(gbx,gby, 10.0,0);
      vdfPlot(400,140, 10.0*512.0/Nmt,MAT1);

    } else if (dispMode==9) {
      neighborPlot(gbx,gby, 0.7, 0);

    }


    gbx = dim.width/6;
    gOff.setColor(Color.black);
    gOff.drawString("t="+(int)(t*1.0e13+0.5)/10.0+" ps",gbx*0+10,70);
    gOff.drawString("temp control",gbx*1+10,70);

    gOff.setColor(Color.black);
    gOff.drawString("cont T="+contTemp+"K",gbx*2+10,70);
    gOff.drawString("scale="+(int)(viewScale*100.0+0.5)+"%",gbx*3+10,70);
    gOff.drawString("view(d="+(int)(zSection*100.0+0.5)/100.0+")",gbx*4+10,70);
    gOff.drawString("thick="+(int)(zThick*100.0+0.5)/100.0+" ",gbx*5+10,70);
	gOff.drawString("secton depth and thickness",gbx*4+0,90);

    gOff.setColor(Color.blue);
    gOff.drawString("T="+(int)(tmp*10+0.5)/10.0+" K",gbx*2+10,90);
    gOff.drawString("L="+(int)(xMax*1.0e10+0.5)/10.0+" nm",gbx*3+10,90);
		
    gOff.drawString("P="+(int)(pressure()/1.0e4+0.5)/100.0+"MPa",gbx*4+0,110);
    gOff.drawString("<P>="+(int)(meanPressure()/1.0e4+0.5)/100.0+"MPa",gbx*4+0,130);
    gOff.drawString("nkT="+(int)((Nmt/(xMax*yMax*zMax)*1.38e-23*tmp)/1.0e4+0.5)/100.0+"MPa",gbx*4+0,150);

    //gOff.setColor(Color.black);
    //gOff.drawString("regMax="+regMax()+"",gbx*5+10,100);
    //gOff.drawString("secMax="+secMax()+"",gbx*5+10,120);

  }


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

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

    drawWaku(gbx,gby,0);

    scl = dispScale*viewScale;
    sz = size*scl;
    for (i=0; i<Nmt; i++) {
      j = srtz[0][i];
	  if (zz[j]>zSection*zMax || zz[j]<(zSection-zThick)*zMax) continue;
      gx = (int)(scl*px[j])+gbx;
      gy = (int)(scl*py[j])+gby;
      gr = (int)(sz*Morse[kind[j]][4]*(0.2*(pz[j]/zDepth)+0.8));
	  if (gr<2) gr = 2;
	  if (imode==1 || imode==2) {
        for (k=1; k<reg[j][0]; k++) {
          jj = reg[j][k];
		  rrr = 2.0*distance(j,jj)/(Morse[kind[j]][4]+Morse[kind[jj]][4]);
          if ( rrr<1.2 ) {
            g2x = (int)(scl*px[jj])+gbx;
            g2y = (int)(scl*py[jj])+gby;
			col = 0.33*(1.21 - rrr)/0.2;
            gOff.setColor(Color.getHSBColor((float)col,0.8f,(float)(0.2*pz[j]/zDepth+0.7)));
            gOff.drawLine(gx,gy, g2x, g2y);
          }
        }
	  }
	  if (imode==0 || imode==2) {
        gOff.setColor(Color.getHSBColor(
          0.15f,0.8f,(float)(0.2*pz[j]/zDepth+0.5)));
        gOff.fillOval(gx-gr/2,gy-gr/2, gr, gr);
      }
      if (imode==3) {
	    tm = Morse[kind[j]][0]*(vx[j]*vx[j]+vy[j]*vy[j]+vz[j]*vz[j])/(3.0*1.38e-23)/600.0;
		if (tm>0.99) tm = 0.99;
        gOff.setColor(Color.getHSBColor(
          0.15f,0.8f,(float)(0.8*tm+0.2)));
        gOff.fillOval(gx-gr/2,gy-gr/2, gr, gr);
      }

    }
    drawWaku(gbx,gby,1);
  }

  double distance(int j, int jj) {
    double r2;

    r2 = (xx[j]-xx[jj])*(xx[j]-xx[jj])+(yy[j]-yy[jj])*(yy[j]-yy[jj])+(zz[j]-zz[jj])*(zz[j]-zz[jj]);
    return ( Math.sqrt(r2) );
  }

  void setNeighbor() {
    int i,j,k,nb1,nb2;
    double rrr;

    for (i=0; i<Nmt; i++) {
	  neighbor[i] = 0;
    }
    for (i=0; i<16; i++) {
      for (j=1; j<30; j++) {
        neighborDistance[i][j] = 0;
      }
    }
	
    for (i=0; i<Nmt; i++) {
      for (k=1; k<reg[i][0]; k++) {
        j = reg[i][k];
        rrr = 2.0*distance(i,j)/(Morse[kind[i]][4]+Morse[kind[j]][4]);
        if ( rrr<1.2 ) {
          neighbor[i] += 1;
          neighbor[j] += 1;
        }
      }
    }

    for (i=0; i<Nmt; i++) {
      nb1 = neighbor[i];
      for (k=1; k<reg[i][0]; k++) {
        j = reg[i][k];
        nb2 = neighbor[j];
        rrr = 2.0*distance(i,j)/(Morse[kind[i]][4]+Morse[kind[j]][4]);
        if ( rrr<1.2 ) {
          neighborDistance[nb1][(int)(rrr*20)] += 1;
          neighborDistance[nb2][(int)(rrr*20)] += 1;
        }
      }
    }
  }

  void neighborPlot(int gbx, int gby, double size, int imode) {
    int i,j,k,jj,gx,gy,g2x,g2y,gr;
	double sz,scl,tm,rrr,col;

    setNeighbor();
    drawWaku(gbx,gby,0);

    scl = dispScale*viewScale;
    sz = size*scl;
    for (i=0; i<Nmt; i++) {
      j = srtz[0][i];
	  if (zz[j]>zSection*zMax || zz[j]<(zSection-zThick)*zMax) continue;
      gx = (int)(scl*px[j])+gbx;
      gy = (int)(scl*py[j])+gby;
      gr = (int)(sz*Morse[kind[j]][4]*(0.2*(pz[j]/zDepth)+0.8));
	  if (gr<2) gr = 2;
	  if (imode==1 || imode==2) {
        for (k=1; k<reg[j][0]; k++) {
          jj = reg[j][k];
		  rrr = 2.0*distance(j,jj)/(Morse[kind[j]][4]+Morse[kind[jj]][4]);
          if ( rrr<1.2 ) {
            g2x = (int)(scl*px[jj])+gbx;
            g2y = (int)(scl*py[jj])+gby;
			col = 0.33*(1.21 - rrr)/0.2;
            gOff.setColor(Color.getHSBColor((float)col,0.8f,(float)(0.2*pz[j]/zDepth+0.7)));
            gOff.drawLine(gx,gy, g2x, g2y);
          }
        }
	  }
	  if (imode==0 || imode==2) {
        gOff.setColor(Color.getHSBColor(
          (float)neighborColor[neighbor[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 = Morse[kind[j]][0]*(vx[j]*vx[j]+vy[j]*vy[j]+vz[j]*vz[j])/(3.0*1.38e-23)/600.0;
		if (tm>0.99) tm = 0.99;
        gOff.setColor(Color.getHSBColor(
          0.15f,0.8f,(float)(0.8*tm+0.2)));
        gOff.fillOval(gx-gr/2,gy-gr/2, gr, gr);
      }

    }
    drawWaku(gbx,gby,1);
    neighborTable(400, 200);
  }

  void neighborTable(int gbx, int gby) {
    int i,j,n;

    gOff.setColor(Color.black);
    for (j=16; j<24; j+=2) {
      gOff.drawString(""+(int)(j*100.0/20.0)/100.0+"-",gbx+(j-14)*20,gby);
    }
    for (i=1; i<=12; i++) {
      gOff.drawString("("+i+")",gbx,gby+20*i);
      for (j=16; j<24; j+=2) {
        n = neighborDistance[i][j] + neighborDistance[i][j+1];
        gOff.drawString(""+n+"",gbx+(j-14)*20,gby+20*i);
      }
    }
  }

  void velocity3DPlot(int gbx, int gby, double mag, int imode) {
    int i,j,gx,gy,g2x,g2y;
    double sc,m,p2x,p2y,p2z,rotp2x,rotp2y,rotp2z;
    double cosTh,sinTh,cosFi,sinFi;

    drawWaku(gbx,gby,0);

    cosTh = Math.cos(theta); sinTh = Math.sin(theta);
    cosFi = Math.cos(fai); sinFi = Math.sin(fai);
    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;

      if (imode==0) {
        p2x = xx[j]+mag*vx[j]*20*dt;
        p2y = yy[j]+mag*vy[j]*20*dt;
        p2z = zz[j]+mag*vz[j]*20*dt;
      } else {
	    m = mag*20*dt*0.01*Morse[kind[j]][0]/AU;
        p2x = xx[j]+m*vx[j];
        p2y = yy[j]+m*vy[j];
        p2z = zz[j]+m*vz[j];
	  }
      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(
        0.15f,0.8f,(float)(0.2*pz[j]/zDepth+0.5)));
      gOff.drawLine(gx,gy, g2x,g2y);
    }

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

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

  //

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

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

  //

  void rotate(double th,double fi) {
    int i;
    double cosTh,sinTh,cosFi,sinFi;

    cosTh = Math.cos(th); sinTh = Math.sin(th);
    cosFi = Math.cos(fi); sinFi = Math.sin(fi);
    for (i=0; i<Nmt; i++) {
      py[i] = cosFi*(yy[i]-cy)+sinFi*(zz[i]-cz) + cy; 
      pz[i] = -sinFi*(yy[i]-cy)+cosFi*(zz[i]-cz) + cz;
    } 
    for (i=0; i<Nmt; i++) {
      px[i] = cosTh*(xx[i]-cx)+sinTh*(pz[i]-cz) + cx; 
      pz[i] = -sinTh*(xx[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;
    } 
    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 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]);
	  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);
  }

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

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

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

    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]);
        rij = Math.sqrt(r2);
        if (rij<10.0e-10) {
          ir = (int)(rij*1.0e11+0.49999);
          rdf[ir] +=1;
        }
      }
    }
  }

  //

  void vdfPlot(int gx, int gy, double mag, int knd) {
    int hoganWidth,hoganHight,xDiv,yDiv;

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

    kindPlot(gx,gy+hoganHight, mag, -1);
    kindPlot(gx,gy+hoganHight, mag, knd);
    drawHogan(gx, gy, hoganWidth, hoganHight, xDiv, yDiv);

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

  void kindPlot(int gx, int ghy, double mag, int knd) {
    int i,iy; 

    setVelocityDistribution(knd);

    gOff.setColor(Color.lightGray);
	if (knd>=0) {
      gOff.setColor(Color.getHSBColor((float)Morse[knd][5],0.8f,0.8f));
	}
    for (i=0; i<200; i++) {
      iy = (int)(mag*vdf[i]);
	  if (iy>0) {
	    gOff.drawLine(gx+i,ghy-iy,gx+i,ghy);
      }
    }
  }

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

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

    for (i=0; i<Nmt; i++) {
      v2 = vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i];
      vv = Math.sqrt(v2);
      if (vv<2000.0) {
        iv = (int)(vv*0.1);
        if (kind[i]==knd || knd==-1) {
          vdf[iv] +=1;
        }
      }
    }
  }



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

  //  set initial condition

  void setInitialCondition() {
    int i,nc,nc2; 
    double s;

    setWaku();

    t = 0.0;

    nc = nCube; nc2 = nCube*nCube;
    s = xMax/(1.0*nCube);
    for (i=0; i<Nmt; i++) {
      kind[i] = MAT1;
      xx[i] = s*(double)((i%nc2)%nc)+0.5*s;
      yy[i] = s*(double)((i%nc2)/nc)+0.5*s;
      zz[i] = s*(double)(i/nc2)+0.5*s;
      vx[i] = 800.0*Math.random()-400.0;
      vy[i] = 800.0*Math.random()-400.0;
      vz[i] = 800.0*Math.random()-400.0;
      ffx[i] = 0.0;
      ffy[i] = 0.0;
      ffz[i] = 0.0;
    }

    vAjustment();
    zSort();
    initVirial();
  }

  //  set forceTable[i][j][ir]

  void setForceTable(){
    int i,j,ir; 
    double d,a,r0,y,r,sc;

    for (i=0; i<NmatKind; i++) {
      for (j=0; j<NmatKind; j++) {

        d = Math.sqrt(Morse[i][2]*Morse[j][2]);
        a = 0.5*(Morse[i][3]+Morse[j][3]);
        r0 = 0.5*(Morse[i][4]+Morse[j][4]);

        for (ir=10; ir<1000; ir++) {
          r = (double)ir*hh;
          y = Math.exp(-a*(r-r0));
          forceTable[i][j][ir] = cutoff(r)*2.0*d*a*y*(y-1.0);
        }
        for (ir=0; ir<10; ir++) {
          forceTable[i][j][ir]=forceTable[i][j][10];
        }

      }
    }

  }

  double cutoff(double r) {
    double rt;

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


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

  void timeEvolution() {
    int i;

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

    if (started==1) {
      registration();
      for (i=0; i<20; i++) {
        timeStep();
      }
      setVirial();
	}

    if (started==0 && stepSW==1) {
      stepSW = 0;
      registration();
      for (i=0; i<20; i++) {
        timeStep();
      }
      setVirial();
	}

    if (tempMode==1) {
      vAjustment();
    }
  }

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

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

    for (i=0; i<Nmt; i++) {
      a2 = dtv2/Morse[kind[i]][0];
      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/Morse[kind[i]][0];
      vx[i] += a2*ffx[i];
      vy[i] += a2*ffy[i];
      vz[i] += a2*ffz[i];
    }

    setBoundary();
  }

  void forceCalc() {
    int i,j,k;
    double r2,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(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 = 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;
        }
      }
    }
  }

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

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

  // registration

  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;

    }
  }

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

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

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

  //

  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] + vz[i]*vz[i]);
    }
    return ((double)Math.floor(s/((double)Nmt*1.5*1.38e-23)*100.0)/100.0 );
  }

  int regMax() {
    int i,m;

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

  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

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

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

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

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

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

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

}