/* plate flow - periodic express Morse 3D */
/* coded by Ikeuchi Mitsuru    */
/* ver 0.0.1   2005.04.23      */

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

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

/*-------------------------------------- define gloval -----*/

  Button bt_go;
  Choice ch_mat,ch_tcon,ch_view;
  Scrollbar sc_flow,sc_temp,sc_scale;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 50;
  double sgm = 3.4e-10;        /* L-J sigma for Ar (m)   */
  double dispScale = 5.0/sgm;

  double xMax = 180.0*sgm;      /* x-Box size in (m)      */
  double yMax = 60.0*sgm;      /* y-Box size in (m)      */
  double zMax = 60.0*sgm;      /* z-Box size in (m)      */


  int Nmt = 60*20*20+1600;          /* number of molecules    */
  int NN = Nmt + 1;
  double t = 0.0;              /* time (s)               */
  double dt = 5.0*1.0e-15;     /* time division (s)      */

  int Nsx = 45;
  int Nsy = 15;
  int Nsz = 15;

  int tempMode = 0;
  double contTemp = 300.0;
  double flowRate = 0.8;
  double scale = 0.6;
  int dispMode = 5;
  int started = 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[][] = {
  /*         mass         charge   D(eV)     A(A^-1)    r0(A)   */
  /* 0 W */ {183.85 *AU,  0.0*EE, 0.9906*EE, 1.4116/AA, 3.032*AA},
  /* 1 Mo*/ { 95.94 *AU,  0.0*EE, 0.8032*EE, 1.5079/AA, 2.976*AA},
  /* 2 Cr*/ { 51.996*AU,  0.0*EE, 0.4414*EE, 1.5721/AA, 2.754*AA},
  /* 3 Fe*/ { 55.847*AU,  0.0*EE, 0.4174*EE, 1.3885/AA, 2.845*AA},
  /* 4 Ni*/ { 58.71 *AU,  0.0*EE, 0.4205*EE, 1.4199/AA, 2.780*AA},
  /* 5 Al*/ { 26.98 *AU,  0.0*EE, 0.2703*EE, 1.1646/AA, 3.253*AA},
  /* 6 Pb*/ {207.19 *AU,  0.0*EE, 0.2348*EE, 1.1836/AA, 3.733*AA},
  /* 7 Cu*/ { 63.54 *AU,  0.0*EE, 0.3429*EE, 1.3588/AA, 2.866*AA},
  /* 8 Ag*/ {107.87 *AU,  0.0*EE, 0.3323*EE, 1.3690/AA, 3.115*AA},
  /* 9 Ca*/ { 40.08 *AU,  0.0*EE, 0.1623*EE, 0.8054/AA, 4.569*AA},
  /*10 Sr*/ { 87.62 *AU,  0.0*EE, 0.1513*EE, 0.7878/AA, 4.988*AA},
  /*11 Ba*/ {137.34 *AU,  0.0*EE, 0.1416*EE, 0.6570/AA, 5.373*AA},
  /*12 Na*/ { 22.99 *AU,  0.0*EE, 0.0633*EE, 0.5900/AA, 5.336*AA},
  /*13 K */ { 39.102*AU,  0.0*EE, 0.0542*EE, 0.4977/AA, 6.369*AA},
  /*14 Rb*/ { 85.47 *AU,  0.0*EE, 0.0464*EE, 0.4298/AA, 7.207*AA},
  /*15 Cs*/ {132.905*AU,  0.0*EE, 0.0449*EE, 0.4157/AA, 7.557*AA},
  /*16 Ne*/ { 20.183*AU,  0.0*EE, 0.0031*EE, 1.6500/AA, 3.076*AA},
  /*17 Ar*/ { 39.948*AU,  0.0*EE, 0.0104*EE, 1.3400/AA, 3.816*AA},
  /*18 Kr*/ { 83.80 *AU,  0.0*EE, 0.0141*EE, 1.2500/AA, 4.097*AA},
  /*19 Xe*/ {131.30 *AU,  0.0*EE, 0.0200*EE, 1.2400/AA, 4.467*AA},
  /*20 Hg*/ {200.59 *AU,  0.0*EE, 0.0734*EE, 1.4900/AA, 3.255*AA},
  /*21 fixed*/ { 1.0e40,  0.0*EE, 0.0104*EE, 1.3400/AA, 3.816*AA}
  };
  int NmatKind = 22;

  int MAT0 = 0;
  int MAT1 = 17;
  int MAT2 = 4;
  int MAT3 = 9;
  int MAT4 = 21;
  
  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][250];   /* registration */
  int section[][][][] = new int[Nsx][Nsy][Nsz][100];

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

  int kindT[] = new int[NmatKind];
  double forceTable[][][] = new double[5][5][10200];

  double hh = 1.0e-13;

  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 Nfx = 30;
  int Nfy = 10;
  int Nfz = 10;
  double flowx[][][] = new double[Nfx+1][Nfy+1][Nfz+1];
  double flowy[][][] = new double[Nfx+1][Nfy+1][Nfz+1];
  double flowz[][][] = new double[Nfx+1][Nfy+1][Nfz+1];
  double flown[][][] = new double[Nfx+1][Nfy+1][Nfz+1];
  double flowt[][][] = new double[Nfx+1][Nfy+1][Nfz+1];

  double meanFlowx[][][] = new double[Nfx+1][Nfy+1][Nfz+1];
  double meanFlowy[][][] = new double[Nfx+1][Nfy+1][Nfz+1];
  double meanFlowz[][][] = new double[Nfx+1][Nfy+1][Nfz+1];
  double meanFlown[][][] = new double[Nfx+1][Nfy+1][Nfz+1];
  double meanFlowt[][][] = new double[Nfx+1][Nfy+1][Nfz+1];

/*----------------------------- applet control functions -----*/

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

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

    bt_go= new Button("go");
    bt_go.addActionListener(this);

    ch_mat = new Choice();
/*
    ch_mat.addItem("W");  ch_mat.addItem("Mo"); ch_mat.addItem("Cr");
    ch_mat.addItem("Fe"); ch_mat.addItem("Ni"); ch_mat.addItem("Al");
    ch_mat.addItem("Pb"); ch_mat.addItem("Cu"); ch_mat.addItem("Ag"); 
    ch_mat.addItem("Ca"); ch_mat.addItem("Sr"); ch_mat.addItem("Ba");
    ch_mat.addItem("Na"); ch_mat.addItem("K");  ch_mat.addItem("Rb");
    ch_mat.addItem("Cs"); ch_mat.addItem("Ne"); 
*/
    ch_mat.addItem("Ar"); ch_mat.addItem("Kr"); ch_mat.addItem("Xe");
/*
    ch_mat.addItem("Hg");
*/
    ch_mat.addItemListener(this);
    ch_mat.select("Ar");

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

    ch_view = new Choice();
    ch_view.addItem("ball"); ch_view.addItem("line"); 
    ch_view.addItem("both"); ch_view.addItem("real");
    ch_view.addItem("temp"); ch_view.addItem("flow");
    ch_view.addItem("coarse T"); ch_view.addItem("coarse N");
    ch_view.addItemListener(this);
    ch_view.select("flow");

    sc_flow= new Scrollbar(Scrollbar.HORIZONTAL,80,10,1,110);
    sc_flow.addAdjustmentListener(this);

    sc_temp= new Scrollbar(Scrollbar.HORIZONTAL,30,10,1,210);
    sc_temp.addAdjustmentListener(this);

    sc_scale= new Scrollbar(Scrollbar.HORIZONTAL,60,10,20,210);
    sc_scale.addAdjustmentListener(this);

    addMouseListener(this);  
    addMouseMotionListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
    pnl.add(ch_mat);
    pnl.add(ch_tcon);
    pnl.add(ch_view);
/*
    pnl.add(new Label("  "));
*/
    pnl.add(sc_flow);
    pnl.add(sc_temp);
    pnl.add(sc_scale);
    add(pnl,"North");

    setInitialPosition();
    setWaku();
  }

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

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

  public void actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_go){
      if (started==0) {
        t = 0.0;
        gogogo();
        started = 1;
      }
    }
  }

  public void itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_mat){
      MAT1 = 17+ch_mat.getSelectedIndex();
      t = 0.0;
	  started = 0;
      setInitialPosition();
    } else if (ev.getSource() == ch_tcon) {
      tempMode = ch_tcon.getSelectedIndex();
    } else if (ev.getSource() == ch_view){
      dispMode = ch_view.getSelectedIndex();
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_flow) { 
      flowRate = 0.01*(double)(sc_flow.getValue());
      vAjustment();
    } else if (ev.getSource() == sc_temp) { 
      contTemp = 10.0*(double)(sc_temp.getValue());
      vAjustment();
    } else if(ev.getSource() == sc_scale) { 
      scale = 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;
    }
    theta += 0.5*3.14/180.0*(dgX-dgXb);
    fai += 0.5*3.14/180.0*(dgY-dgYb);
  }

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

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

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


/* ----------------------------- offPaint -------------------- */

  void offPaint() {
    int i,j,k,jj,farPoint;
    int gx,gy,gz,g2x,g2y,gbx,gby,gbz,gb;
    double sc,sz,tmp,tm,clr;

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

    rotate(theta,fai);
    zSort();
    tmp = temperature();
    gbx = 40;
    gby = 100;
    sc = dispScale*scale;
    sz = 3.6*scale/AA;

    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 (boxp[i][0]==farPoint || boxp[i][1]==farPoint) { 
        gOff.setColor(Color.lightGray);
        gOff.drawLine(gx,gy, g2x, g2y);
      } 
    }


    if (dispMode==0) {

      ballPlot(gbx,gby,sc);

    } else if (dispMode==1) {

      linePlot(gbx,gby,sc);

    } else if (dispMode==2) {

      bothPlot(gbx,gby,sc);

    } else if (dispMode==3) {

      realPlot(gbx,gby,sc);

    } else if (dispMode==4) {

      tempPlot(gbx,gby,sc);

    } else if (dispMode==5) {
	
      flowPlot(gbx, gby, sc);
	  
    } else if (dispMode==6) {
	
      coarseTempPlot(gbx, gby, sc);
	  
    } else if (dispMode==7) {
	
      coarseDensPlot(gbx, gby, sc);
	  
    }


    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 (boxp[i][0]!=farPoint && boxp[i][1]!=farPoint) { 
        gOff.setColor(Color.gray);
        if (tempMode==2) {
          if (tmp>contTemp) {
            gOff.setColor(Color.black); 
          } else {
            gOff.setColor(Color.red);
          }
        }
        gOff.drawLine(gx,gy, g2x, g2y);
      }
    }


    gbx = dim.width/6;
    gOff.setColor(Color.black);
    gOff.drawString("atom",gbx*0+10,40);
    gOff.drawString("temp control",gbx*1+10,40);
    gOff.drawString("view",gbx*2+10,40);
    gOff.drawString("flow="+(int)(flowRate*100.0+0.5)+"%",gbx*3+10,40);
    gOff.drawString("cont T="+contTemp+"K",gbx*4+10,40);
    gOff.drawString("scale="+(int)(scale*100.0+0.499)+"%",gbx*5+10,40);

    gOff.setColor(Color.blue);
    gOff.drawString("t="+(int)(t*1.0e13+0.5)/10.0+" ps",gbx*0+10,60);
    gOff.drawString("T="+(int)(tmp*10+0.5)/10.0+" K",gbx*4+10,60);
    gOff.drawString("L="+(int)(xMax*1.0e10+0.5)/10.0+" nm",gbx*5+10,60);
/*
    gOff.setColor(Color.black);
    gOff.drawString("regMax="+regMax()+"",gbx*5+10,100);
    gOff.drawString("secMax="+secMax()+"",gbx*5+10,120);
*/
  }

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

  double ballColor(int knd) {
    double c;

    c = 0.0;
    if (knd==MAT0) {
      c = 0.01;
    } else if (knd==MAT1) {
      c = 0.33;
    } else if (knd==MAT2) {
      c = 0.48;
    } else if (knd==MAT3) {
      c = 0.66;
    } else if (knd==MAT4) {
      c = 0.15;
    }
    return ( c );
  }

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

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

  void ballPlot(int gbx, int gby, double scl) {
    int i,j,gx,gy,gz;
	double sz,clr;
	
      sz = 3.6*scale/AA;
      for (i=0; i<Nmt; i++) {
        j = srtz[0][i];
        gx = (int)(scl*px[j])+gbx;
        gy = (int)(scl*py[j])+gby;
        gz = (int)(0.5*sz*Morse[kind[j]][4]*(0.2*(pz[j]/zMax)+0.7));
        clr = ballColor(kind[j]);
        gOff.setColor(Color.getHSBColor(
          (float)clr,0.8f,(float)(0.2*pz[j]/zMax+0.5)));
        gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
	}
  }

  void linePlot(int gbx, int gby, double scl) {
    int i,j,k,jj,gx,gy,g2x,g2y;
	double sz;
	
    sz = 3.6*scale/AA;

      for (i=0; i<Nmt; i++) {
        j = srtz[0][i];
        gx = (int)(scl*px[j])+gbx;
        gy = (int)(scl*py[j])+gby;
        for (k=1; k<reg[j][0]; k++) {
          jj = reg[j][k];
          if ( distance(j,jj)<0.6*(Morse[kind[j]][4]+Morse[kind[jj]][4]) ) {
              g2x = (int)(scl*px[jj])+gbx;
              g2y = (int)(scl*py[jj])+gby;
              gOff.setColor(Color.gray);
              gOff.drawLine(gx,gy, g2x, g2y);
          }
        }
      }
  }

  void bothPlot(int gbx, int gby, double scl) {
    int i,j,k,jj,gx,gy,g2x,g2y,gz;
	double sz,clr;
	
    sz = 3.6*scale/AA;

      for (i=0; i<Nmt; i++) {
        j = srtz[0][i];
        gx = (int)(scl*px[j])+gbx;
        gy = (int)(scl*py[j])+gby;
        gz = (int)(0.5*sz*Morse[kind[j]][4]*(0.2*(pz[j]/zMax)+0.7));
        for (k=1; k<reg[j][0]; k++) {
          jj = reg[j][k];
          if ( distance(j,jj)<0.6*(Morse[kind[j]][4]+Morse[kind[jj]][4]) ) {
              g2x = (int)(scl*px[jj])+gbx;
              g2y = (int)(scl*py[jj])+gby;
              gOff.setColor(Color.gray);
              gOff.drawLine(gx,gy, g2x, g2y);
          }
        }
		clr = ballColor(kind[j]);
		gOff.setColor(Color.getHSBColor(
          (float)clr,0.8f,(float)(0.2*pz[j]/zMax+0.5)));
        gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
      }
  }

  void realPlot(int gbx, int gby, double scl) {
    int i,j,gx,gy,gz;
    double sz,clr;
	
      sz = 3.6*scale/AA;
      for (i=0; i<Nmt; i++) {
        j = srtz[0][i];
        gx = (int)(scl*px[j])+gbx;
        gy = (int)(scl*py[j])+gby;
        gz = (int)(0.9*sz*Morse[kind[j]][4]*(0.2*(pz[j]/zMax)+0.7));
        clr = ballColor(kind[j]);
        gOff.setColor(Color.getHSBColor(
          (float)clr,0.8f,(float)(0.2*pz[j]/zMax+0.5)));
        gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
	  }
  }

  void tempPlot(int gbx, int gby, double scl) {
    int i,j,gx,gy,gz;
    double sz,tm,clr;
	
      sz = 3.6*scale/AA;
      for (i=0; i<Nmt; i++) {
        j = srtz[0][i];
        gx = (int)(scl*px[j])+gbx;
        gy = (int)(scl*py[j])+gby;
        gz = (int)(0.5*sz*Morse[kind[j]][4]*(0.2*(pz[j]/zMax)+0.7));
        clr = ballColor(kind[j]);
        tm = Morse[kind[j]][0]*(vx[j]*vx[j]+vy[j]*vy[j]+vz[j]*vz[j])/(3.0*1.38e-23)/6000.0;
        if(tm>0.66) {
          tm = 0.66;
        }
        gOff.setColor(Color.getHSBColor(
          (float)(0.66-tm),0.9f,0.9f));
        gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
	}
  }

  void flowPlot(int gbx, int gby, double sc) {
    int i,j,ix,iy,iz,gx,gy,gz,g2x,g2y;
    double cosTh,sinTh,cosFi,sinFi;
    double x,y,z,x2,y2,z2,ppx,ppy,ppz,ppx2,ppy2,ppz2,sz,clr;	

    cosTh = Math.cos(theta); sinTh = Math.sin(theta);
    cosFi = Math.cos(fai); sinFi = Math.sin(fai);
    sz = 3.6*scale/AA;

      for (i=0; i<Nmt; i++) {
        j = srtz[0][i];
		if (kind[j]==MAT4) {
        gx = (int)(sc*px[j])+gbx;
        gy = (int)(sc*py[j])+gby;
        gz = (int)(0.5*sz*Morse[kind[j]][4]*(0.2*(pz[j]/zMax)+0.7));
        clr = ballColor(kind[j]);
        gOff.setColor(Color.getHSBColor(
          (float)clr,0.8f,(float)(0.2*pz[j]/zMax+0.5)));
        gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
		}
	  }

      for (ix=0; ix<Nfx; ix++) {
	    x = (ix+0.5)*xMax/Nfx;
	    for (iy=0; iy<Nfy; iy++) {
		  y = (iy+0.5)*yMax/Nfy;
	      for (iz=0; iz<Nfy; iz++) {
		    z = (iz+0.5)*zMax/Nfz;
			
            ppy = cosFi*(y-cy)+sinFi*(z-cz) + cy; 
            ppz = -sinFi*(y-cy)+cosFi*(z-cz) + cz;
            ppx = cosTh*(x-cx)+sinTh*(ppz-cz) + cx; 
            ppz = -sinTh*(x-cx)+cosTh*(ppz-cz) + cz;
            gx = (int)(sc*ppx)+gbx;
            gy = (int)(sc*ppy)+gby;

            x2 = x + 300.0*meanFlowx[ix][iy][iz]*dt;
            y2 = y + 300.0*meanFlowy[ix][iy][iz]*dt;
            z2 = z + 300.0*meanFlowz[ix][iy][iz]*dt;
            ppy2 = cosFi*(y2-cy)+sinFi*(z2-cz) + cy; 
            ppz2 = -sinFi*(y2-cy)+cosFi*(z2-cz) + cz;
            ppx2 = cosTh*(x2-cx)+sinTh*(ppz2-cz) + cx; 
            ppz2 = -sinTh*(x2-cx)+cosTh*(ppz2-cz) + cz;
            g2x = (int)(sc*ppx2)+gbx;
            g2y = (int)(sc*ppy2)+gby;
            if (meanFlowx[ix][iy][iz]>0) {
              gOff.setColor(Color.blue);
            } else {
              gOff.setColor(Color.red);
            }
            gOff.drawLine(gx,gy, g2x, g2y);
		}
	  }
    }
  }

  void coarseTempPlot(int gbx, int gby, double sc) {
    int i,j,ix,iy,iz,gx,gy,gz;
    double cosTh,sinTh,cosFi,sinFi;
    double x,y,z,x2,y2,z2,ppx,ppy,ppz,sz,tm,clr;	

	cosTh = Math.cos(theta); sinTh = Math.sin(theta);
    cosFi = Math.cos(fai); sinFi = Math.sin(fai);
    sz = 3.6*scale/AA;

      for (i=0; i<Nmt; i++) {
        j = srtz[0][i];
		if (kind[j]==MAT4) {
        gx = (int)(sc*px[j])+gbx;
        gy = (int)(sc*py[j])+gby;
        gz = (int)(0.5*sz*Morse[kind[j]][4]*(0.2*(pz[j]/zMax)+0.7));
        clr = ballColor(kind[j]);
        gOff.setColor(Color.getHSBColor(
          (float)clr,0.8f,(float)(0.2*pz[j]/zMax+0.5)));
        gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
		}
	  }

      gz = 6;
      for (ix=0; ix<Nfx; ix++) {
	    x = (ix+0.5)*xMax/Nfx;
	    for (iy=0; iy<Nfy; iy++) {
		  y = (iy+0.5)*yMax/Nfy;
	      for (iz=0; iz<Nfy; iz++) {
		    z = (iz+0.5)*zMax/Nfz;
			
            ppy = cosFi*(y-cy)+sinFi*(z-cz) + cy; 
            ppz = -sinFi*(y-cy)+cosFi*(z-cz) + cz;
            ppx = cosTh*(x-cx)+sinTh*(ppz-cz) + cx; 
            ppz = -sinTh*(x-cx)+cosTh*(ppz-cz) + cz;
            gx = (int)(sc*ppx)+gbx;
            gy = (int)(sc*ppy)+gby;

            tm = meanFlowt[ix][iy][iz]/meanFlown[ix][iy][iz]/2000;
			if (tm>0.66) tm = 0.66;
            gOff.setColor(Color.getHSBColor( (float)(0.66-tm), 0.9f, 0.9f ));
            gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
		}
	  }
    }
  }

  void coarseDensPlot(int gbx, int gby, double sc) {
    int i,j,ix,iy,iz,gx,gy,gz;
    double cosTh,sinTh,cosFi,sinFi;
    double x,y,z,x2,y2,z2,ppx,ppy,ppz,sz,dens,clr;	

    cosTh = Math.cos(theta); sinTh = Math.sin(theta);
    cosFi = Math.cos(fai); sinFi = Math.sin(fai);
    sz = 3.6*scale/AA;

    for (i=0; i<Nmt; i++) {
      j = srtz[0][i];
      if (kind[j]==MAT4) {
        gx = (int)(sc*px[j])+gbx;
        gy = (int)(sc*py[j])+gby;
        gz = (int)(0.5*sz*Morse[kind[j]][4]*(0.2*(pz[j]/zMax)+0.7));
        clr = ballColor(kind[j]);
        gOff.setColor(Color.getHSBColor(
          (float)clr,0.8f,(float)(0.2*pz[j]/zMax+0.5)));
        gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
      }
    }

    gz = 6;
    for (ix=0; ix<Nfx; ix++) {
      x = (ix+0.5)*xMax/Nfx;
      for (iy=0; iy<Nfy; iy++) {
        y = (iy+0.5)*yMax/Nfy;
	       for (iz=0; iz<Nfy; iz++) {
		        z = (iz+0.5)*zMax/Nfz;
			
          ppy = cosFi*(y-cy)+sinFi*(z-cz) + cy; 
          ppz = -sinFi*(y-cy)+cosFi*(z-cz) + cz;
          ppx = cosTh*(x-cx)+sinTh*(ppz-cz) + cx; 
          ppz = -sinTh*(x-cx)+cosTh*(ppz-cz) + cz;
          gx = (int)(sc*ppx)+gbx;
          gy = (int)(sc*ppy)+gby;

          dens = 0.03*meanFlown[ix][iy][iz];
          if (dens>0.99) dens = 0.99;
          gOff.setColor(Color.getHSBColor( (float)(0.66-0.66*dens), 0.9f, 0.9f ));
          gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
        }
      }
    }
  }


/*----------------------------------  findMin()  ----------*/

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

/*-------------------------------------  sort z  ----------*/

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

    for (i=0; i<Nmt; i++) {
      j = srtz[0][i];
      srtz[1][i] = (int)(32767.0*pz[j]/zMax);
    } 
    for (i=0; i<100; i++) {
      for (j=0; j<Nmt-1; j++) {
        if (srtz[1][j]>srtz[1][j+1]) {
          w = srtz[0][j]; srtz[0][j]=srtz[0][j+1]; srtz[0][j+1]=w;
          w = srtz[1][j]; srtz[1][j]=srtz[1][j+1]; srtz[1][j+1]=w;
        }
      }
    }    
  }

  void zSortInit() {
    int i,j,w;

    for (i=0; i<Nmt; i++) {
      srtz[0][i] = i;
      srtz[1][i] = (int)(32767.0*pz[i]/zMax);
    } 
    for (i=0; i<Nmt; i++) {
      for (j=0; j<Nmt-1; j++) {
        if (srtz[1][j]>srtz[1][j+1]) {
          w = srtz[0][j]; srtz[0][j]=srtz[0][j+1]; srtz[0][j+1]=w;
          w = srtz[1][j]; srtz[1][j]=srtz[1][j+1]; srtz[1][j+1]=w;
        }
      }
    }    
  }

/*-------------------------------------  rotate  ----------*/

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

/*-------------------------------------  setWaku  ---------*/

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

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

/*-------------------------------------  statictics  ----------*/


  double rmsV() {
    int i;
    double s;

    s = 0.0; 
    for (i=0; i<Nmt; i++) {
      s += vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i];
    }
    return ((double)Math.floor(Math.sqrt(s/(double)Nmt)*10.0)/10.0);
  }

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

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

    for (ix=0; ix<Nfx; ix++) {
      for (iy=0; iy<Nfy; iy++) {
        for (iz=0; iz<Nfz; iz++) {
          meanFlowx[ix][iy][iz] = 0.0;
          meanFlowy[ix][iy][iz] = 0.0;
          meanFlowz[ix][iy][iz] = 0.0;
          meanFlown[ix][iy][iz] = 0.0;
          meanFlowt[ix][iy][iz] = 0.0;
        }
      }
    }
  }

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

    r = 0.01;
    for (ix=0; ix<Nfx; ix++) {
      for (iy=0; iy<Nfy; iy++) {
        for (iz=0; iz<Nfz; iz++) {
          flowx[ix][iy][iz] = 0.0;
          flowy[ix][iy][iz] = 0.0;
          flowz[ix][iy][iz] = 0.0;
          flown[ix][iy][iz] = 0.0;
          flowt[ix][iy][iz] = 0.0;
        }
      }
    }

    for (i=0; i<Nmt; i++) {
      if (kind[i]>0) {
        ix = (int)(xx[i]*Nfx/xMax);
        iy = (int)(yy[i]*Nfy/yMax);
		iz = (int)(zz[i]*Nfz/zMax);
        flowx[ix][iy][iz] += vx[i];
        flowy[ix][iy][iz] += vy[i];
        flowz[ix][iy][iz] += vz[i];
        flown[ix][iy][iz] += 1.0;
        flowt[ix][iy][iz] += (Morse[kind[i]][0])*(vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i])/(3.0*1.38e-23);
      }
    }

    for (ix=0; ix<Nfx; ix++) {
      for (iy=0; iy<Nfy; iy++) {
        for (iz=0; iz<Nfz; iz++) {
          meanFlowx[ix][iy][iz] = (1.0-r)*meanFlowx[ix][iy][iz]+r*flowx[ix][iy][iz];
          meanFlowy[ix][iy][iz] = (1.0-r)*meanFlowy[ix][iy][iz]+r*flowy[ix][iy][iz];
          meanFlowz[ix][iy][iz] = (1.0-r)*meanFlowz[ix][iy][iz]+r*flowz[ix][iy][iz];
          meanFlown[ix][iy][iz] = (1.0-r)*meanFlown[ix][iy][iz]+r*flown[ix][iy][iz];
          meanFlowt[ix][iy][iz] = (1.0-r)*meanFlowt[ix][iy][iz]+r*flowt[ix][iy][iz];
        }
      }
    }
  }

/*-------------------------------------  v ajustment  ---------*/

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

/*-------------------------------------  gogogo  ---------*/

  void gogogo() {
    int i;
    double x,y,z,x0,y0,z0,rr;

    x0 = 0.17*xMax; y0 = 0.5*yMax; z0 = 0.5*zMax;
    rr = 0.35*yMax;
    for (i=0; i<Nmt; i++) {
      x = xx[i] - x0;
      y = yy[i] - y0;
	  z = zz[i] - z0;
      if (x*x+y*y+z*z<rr*rr) {
        vx[i] = 4000.0*Math.random()-2000.0;
        vy[i] = 4000.0*Math.random()-2000.0;
		vz[i] = 4000.0*Math.random()-2000.0;
      }
    }
  }


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

  void setForceTable(){
    int i,j,ir,k0,k1; 
    double d,a,r0,y,r,sc;
    int kkk[] = {MAT0,MAT1,MAT2,MAT3,MAT4};

    kindT[MAT0] = 0; kindT[MAT1] = 1; kindT[MAT2] = 2;
    kindT[MAT3] = 3; kindT[MAT4] = 4;

    for (i=0; i<5; i++) {
      k0=kkk[i];
      for (j=0; j<5; j++) {
        k1=kkk[j];

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

        for (ir=100; ir<10000; 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<100; ir++) {
          forceTable[i][j][ir]=forceTable[i][j][100];
        }

      }
    }

  }

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

/*-----------------------------  set initial position   ---------*/

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

    setForceTable();

    s = xMax/60.0;
	i=0;
	i = setPlate(i, 0.3*xMax, 0.2*yMax);
	i = fillGas(i, 0.3*xMax);
    Nmt = i - 1;

    setFlowInit();
    vAjustment();
    zSortInit();
  }

  int setPlate(int iCount, double x0, double r0) {
    int i,iy,iz;
    double a,b,r2,x,y,z;

    i = iCount;
    a = 1.2*sgm;
	b = 0.866*a;
	r2 = r0*r0;
    x = x0;

    for (iy=0; iy<40; iy++) {
      for (iz=0; iz<40; iz++) {
	    y = 0.5*yMax + a*(iy-20.0)+0.5*a*(iz%2);
	    z = 0.5*zMax + b*(iz-20.0);
        if ( (y-0.5*yMax)*(y-0.5*yMax)+(z-0.5*zMax)*(z-0.5*zMax)<r2 ) {
          kind[i] = MAT4;
          xx[i] = x;
          yy[i] = y;
          zz[i] = z;
          vx[i] = 0.0;
          vy[i] = 0.0;
          vz[i] = 0.0;
          ffx[i] = 0.0;
          ffy[i] = 0.0;
          ffz[i] = 0.0;
		  i += 1;
		}
      }
    }
	return(i);
  }

  int fillGas(int iCount, double x0) {
    int i,ix,iy,iz;
    double s,x,y,z;

    i = iCount;
    s = xMax/60.0;

    for (ix=0; ix<60; ix++) {
	  x = s*(double)(ix)+0.5*s;
      for (iy=0; iy<20; iy++) {
	    y = s*(double)(iy)+0.5*s;
        for (iz=0; iz<20; iz++) {
	      z = s*(double)(iz)+0.5*s;
          if ( x< x0-sgm || x>x0+sgm) {
          kind[i] = MAT1;
          xx[i] = x;
          yy[i] = y;
          zz[i] = z;
          vx[i] = 800.0*Math.random()-400.0+100.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;
		  i += 1;
		  }
		}
      }
    }
	return(i);
  }

/*-------------------------------  registration   -------------*/

  void registration() {
    int i,j,k,ii,jj,kk,ip,jp,kp,i0,i1,j0,j1,k0,k1,iq;
    double r2,rreg,rreg2;
    double xi,xj,yi,yj,zi,zj;

    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+(double)Nsx)-Nsx;
      i1 = (int)(Nsx*(xx[ip]+rreg)/xMax);
      j0 = (int)(Nsy*(yy[ip]-rreg)/yMax+(double)Nsy)-Nsy;
      j1 = (int)(Nsy*(yy[ip]+rreg)/yMax);
      k0 = (int)(Nsz*(zz[ip]-rreg)/zMax+(double)Nsz)-Nsz;
      k1 = (int)(Nsz*(zz[ip]+rreg)/zMax);

      for(i=i0;i<=i1;i++) {
	    ii = (i+Nsx)%Nsx;
        for(j=j0;j<=j1;j++) {
		  jj = (j+Nsy)%Nsy;
		  for(k=k0;k<=k1;k++) {
            kk = (k+Nsz)%Nsz;
	        for(iq=1;iq<=section[ii][jj][kk][0];iq++) {
              jp = section[ii][jj][kk][iq];
              if (jp>ip) {
                xi = xx[ip]; xj = xx[jp]; yi = yy[ip]; yj = yy[jp]; zi = zz[ip]; zj = zz[jp];
                if (xi-xj>0.5*xMax) xj += xMax;
                if (xi-xj<-0.5*xMax) xj -= xMax;
                if (yi-yj>0.5*yMax) yj += yMax;
                if (yi-yj<-0.5*yMax) yj -= yMax;
                if (zi-zj>0.5*zMax) zj += zMax;
                if (zi-zj<-0.5*zMax) zj -= zMax;
                r2=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj);
                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;
    }
  }

/*-------------------------------  molecules motion   ---------*/

  void calcposition(){
    int i;

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

    registration();

    for (i=0; i<20; i++) {
      timeEvolution();
    }
  }

  void timeEvolution(){
    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/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];
    }

    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) {
		if (Math.random()<flowRate) {
          xx[i] -= xMax;
        } else {
          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] += yMax; 
      }
      if (yy[i] > yMax) {
        yy[i] -= yMax;
      }
      if (zz[i] < 0.0) { 
        zz[i] += zMax; 
      }
      if (zz[i] > zMax) {
        zz[i] -= zMax; 
      }
    }
  }

  void forceCalc() {
    int i,j,k;
    double xi,xj,yi,yj,zi,zj;
    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];
        xi = xx[i]; xj = xx[j]; yi = yy[i]; yj = yy[j]; zi = zz[i]; zj = zz[j];
        if (xi-xj>0.5*xMax) xj += xMax;
        if (xi-xj<-0.5*xMax) xj -= xMax;
        if (yi-yj>0.5*yMax) yj += yMax;
        if (yi-yj<-0.5*yMax) yj -= yMax;
        if (zi-zj>0.5*zMax) zj += zMax;
        if (zi-zj<-0.5*zMax) zj -= zMax;
        r2 = (xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj);
        if (r2<rc2) {
          rij = Math.sqrt(r2);
          f = force(rij,kind[i],kind[j]);
          fxij = f*(xi-xj)/rij;
          fyij = f*(yi-yj)/rij;
          fzij = f*(zi-zj)/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[kindT[ki]][kindT[kj]][i] + 
      a*forceTable[kindT[ki]][kindT[kj]][i+1])/hh );
  }

/*----- end of molecules motion -----*/

}  

/*----- end of applet -----*/