/* curve ball - periodic express Morse 2D */
/* coded by Ikeuchi Mitsuru   */
/* ver 0.0.1   2005.04.02     */

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

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

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

  Choice ch_mat,ch_tcon,ch_view;
  Scrollbar 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 = 6.0/sgm;

  double xMax = 300.0*sgm;      /* x-Box size in (m)      */
  double yMax = 200.0*sgm;      /* y-Box size in (m)      */


  int Nmt = 150*100+4000;             /* 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 = 60; /* xMax/5.0  */
  int Nsy = 40; /* yMax/5.0  */

  int Ncx = 30;
  int Ncy = 20;

  int tempMode = 0;
  double contTemp = 300.0;
  double scale = 0.27;
  int dispMode = 5;
  int resetFlag = 0;
  int resetCount = 0;

  int dgX, dgY, dgXb, dgYb;
  double xShift = 0.0;
  double yShift = 0.0;

  double rc = 2.6*sgm;
  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 = 17;
  int MAT1 = 2;
  int MAT2 = 3;
  int MAT3 = 4;
  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 vx[] = new double[NN];   /* i-th x-velocity */
  double vy[] = new double[NN];   /* i-th y-velocity */
  double ffx[] = new double[NN];  /* i-th x-force    */
  double ffy[] = new double[NN];  /* i-th y-force    */

  int reg[][] = new int[NN][200];   /* registration */
  int section[][][] = new int[Nsx][Nsy][100];

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

  double flowx[][] = new double[Ncx+1][Ncy+1];
  double flowy[][] = new double[Ncx+1][Ncy+1];
  double flown[][] = new double[Ncx+1][Ncy+1];
  double flowt[][] = new double[Ncx+1][Ncy+1];

  double meanFlowx[][] = new double[Ncx+1][Ncy+1];
  double meanFlowy[][] = new double[Ncx+1][Ncy+1];
  double meanFlown[][] = new double[Ncx+1][Ncy+1];
  double meanFlowt[][] = new double[Ncx+1][Ncy+1];

  double traceBallPos[][] = new double[127][2];
  int traceCount = 0;

  double hh = 1.0e-13;

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

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

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

    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_tcon = new Choice();
    ch_tcon.addItem("adiabatic"); ch_tcon.addItem("t scaling"); 
    ch_tcon.addItemListener(this);

    ch_view = new Choice();
    ch_view.addItem("ball"); ch_view.addItem("line"); 
    ch_view.addItem("both"); ch_view.addItem("temp");
    ch_view.addItem("dens"); ch_view.addItem("flow");ch_view.addItem("mean T");
    ch_view.addItemListener(this);
    ch_view.select("flow");

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

    sc_scale= new Scrollbar(Scrollbar.HORIZONTAL,27,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(new Label("  "));
    pnl.add(ch_tcon);
    pnl.add(ch_view);
    pnl.add(new Label("  "));
    pnl.add(sc_temp);
    pnl.add(sc_scale);
    add(pnl,"North");

    setInitialPosition();
  }

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

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

  public void itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_mat) {
      MAT1 = ch_mat.getSelectedIndex();
      t = 0.0;
	  resetFlag = 1;
      resetCount = 2;
    } 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_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;
    }
    xShift += (1.0/scale/dispScale)*(dgX-dgXb);
    yShift += (1.0/scale/dispScale)*(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,ix,iy;
    int gball,gbx,gby,gb;
    double sc,sz,tmp;

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

    setFlow();
    sc = scale*dispScale;
    sz = 0.8*scale*dispScale;
    gbx = 30+(int)(sc*xShift);
    gby = 65+(int)(sc*yShift);
    gball = (int)(6.0*scale);

    tmp = temperature();


    if (dispMode==0) {
      dispBalls(gbx,gby);

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

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

    } else if (dispMode==3) {
      dispTempBalls(gbx,gby);

    } else if (dispMode==4) {
      dispDensity(gbx,gby);

    } else if (dispMode==5) {
      dispFlow(gbx,gby);

    } else if (dispMode==6) {
      dispMeanTemp(gbx,gby);

    }


    gOff.setColor(Color.gray);
    ix = (int)(sc*xMax);
    iy = (int)(sc*yMax);
    gOff.fillRect(gbx-gball/2,gby-gball/2,gball,iy+gball);
    gOff.fillRect(gbx-gball/2,gby-gball/2,ix+gball,gball);
    gOff.fillRect(gbx+ix-gball/2+1,gby-gball/2,gball,iy+gball);
    gOff.fillRect(gbx-gball/2,gby+iy-gball/2+1,ix+gball,gball);

    gOff.drawLine(gbx-gball/2,gby+iy/2, gbx+ix-gball/2,gby+iy/2);
    gOff.drawLine(gbx+ix/3,gby-gball/2, gbx+ix/3,gby+iy+gball/2);

    gOff.setColor(Color.black);
    for (i=0; i<100; i++) {
      ix = gbx+(int)(sc*traceBallPos[i][0]);
      iy = gby+(int)(sc*traceBallPos[i][1]);
      gOff.fillOval(ix-3,iy-3, 6, 6);
    }

    gbx = dim.width/6;
    gOff.setColor(Color.black);
    gOff.drawString("temp control",gbx*1+10,40);
    gOff.drawString("view",gbx*2+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*3+10,40);
    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.drawString("ball position",530,100);
    gOff.drawString("X="+(int)(ballPosX()*1.0e11+0.5)/100.0+" nm",530,120);
    gOff.drawString("Y="+(int)(ballPosY()*1.0e11+0.5)/100.0+" nm",530,140);
  }

  int regMax() {
    int i,m;

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

/*--------------------------------- display methods  ----------*/

  void dispBalls(int gx, int gy) {
    int i,ir,ix,iy;
    double sc,sz;

    sc = scale*dispScale;
    sz = 0.8*scale*dispScale;
    for (i=0; i<Nmt; i++) {
      ix = (int)(sc*xx[i])+gx;
      iy = (int)(sc*yy[i])+gy;
      ir = (int)(sz*Morse[kind[i]][4] );
      gOff.setColor(Color.getHSBColor((float)(colorOf(i)),0.9f,0.9f));
      gOff.fillOval(ix-ir/2,iy-ir/2, ir, ir);
    }
  }

  void dispLines(int gx, int gy) {
    int i,j,k,ir,ix,iy,ix2,iy2;
    double sc,sz,r2;

    sc = scale*dispScale;
    sz = 0.8*scale*dispScale;
    for (i=0; i<Nmt; i++) {
      ix = (int)(sc*xx[i])+gx;
      iy = (int)(sc*yy[i])+gy;

      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]);
        if (r2<0.26*(Morse[kind[i]][4]+Morse[kind[j]][4])*(Morse[kind[i]][4]+Morse[kind[j]][4]) && j>i) {
          ix2 = (int)(sc*xx[j])+gx;
          iy2 = (int)(sc*yy[j])+gy;
          gOff.setColor(Color.gray);
          gOff.drawLine(ix,iy, ix2,iy2);
        }
      }
    }
  }

  void dispLinesAndBalls(int gx, int gy) {
    int i,j,k,ir,ix,iy,ix2,iy2;
    double sc,sz,r2;

    sc = scale*dispScale;
    sz = 0.8*scale*dispScale;
    for (i=0; i<Nmt; i++) {
      ix = (int)(sc*xx[i])+gx;
      iy = (int)(sc*yy[i])+gy;
      ir = (int)(0.6*sz*Morse[kind[i]][4] );

      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]);
        if (r2<0.26*(Morse[kind[i]][4]+Morse[kind[j]][4])*(Morse[kind[i]][4]+Morse[kind[j]][4]) && j>i) {
          ix2 = (int)(sc*xx[j])+gx;
          iy2 = (int)(sc*yy[j])+gy;
          gOff.setColor(Color.gray);
          gOff.drawLine(ix,iy, ix2,iy2);
        }
      }
      gOff.setColor(Color.getHSBColor((float)(colorOf(i)),0.95f,0.95f));
      gOff.fillOval(ix-ir/2,iy-ir/2, ir, ir);
    }
  }

  void dispTempBalls(int gx, int gy) {
    int i,ir,ix,iy;
    double sc,sz,tm;

    sc = scale*dispScale;
    sz = 0.8*scale*dispScale;
    for (i=0; i<Nmt; i++) {
      ix = (int)(sc*xx[i])+gx;
      iy = (int)(sc*yy[i])+gy;
      ir = (int)(sz*Morse[kind[i]][4] );
	  
      tm = 0.5*Morse[kind[i]][0]*(vx[i]*vx[i]+vy[i]*vy[i])/1.38e-23;
      tm = tm/3000.0; if (tm>0.7) tm = 0.7;
      gOff.setColor(Color.getHSBColor((float)(colorOf(i)),0.95f,(float)(tm+0.29)));
      gOff.fillOval(ix-ir/2,iy-ir/2, ir, ir);
    }
  }

  void dispDensity(int gx, int gy) {
    int i,j,ix,iy;
    double sc,sz,dn;

    sc = scale*dispScale;
    sz = 0.8*scale*dispScale;
    for (i=0; i<Ncx; i++) {
      for (j=0; j<Ncy; j++) {
        ix = (int)(sc*(i)*xMax/Ncx)+gx;
        iy = (int)(sc*(j)*yMax/Ncy)+gy;
        dn = 0.02*flown[i][j]; if (dn>0.99) dn = 0.99;
        gOff.setColor(Color.getHSBColor(0.33f,(float)dn,0.9f));
        gOff.fillRect(ix,iy,(int)(sc*xMax/Ncx)+1,(int)(sc*yMax/Ncy)+1);
      }
    }
    setBalls(gx,gy);
  }

  void dispFlow(int gx, int gy) {
    int i,j,ix,iy,ix2,iy2;
    double sc,sz;

    sc = scale*dispScale;
    sz = 0.8*scale*dispScale;
    for (i=0; i<Ncx; i++) {
      for (j=0; j<Ncy; j++) {
        ix = (int)(sc*(i+0.5)*xMax/Ncx)+gx;
        iy = (int)(sc*(j+0.5)*yMax/Ncy)+gy;
        ix2 = ix + (int)(0.01*meanFlowx[i][j]);
        iy2 = iy + (int)(0.01*meanFlowy[i][j]);
        if (meanFlowx[i][j]>0) {
          gOff.setColor(Color.blue);
        } else {
          gOff.setColor(Color.red);
        }
        gOff.drawLine(ix,iy, ix2,iy2);
      }
    }
    setBalls(gx,gy);
  }

  void dispMeanTemp(int gx, int gy) {
    int i,j,ix,iy;
    double sc,sz,tm;

    sc = scale*dispScale;
    sz = 0.8*scale*dispScale;
    for (i=0; i<Ncx; i++) {
      for (j=0; j<Ncy; j++) {
        ix = (int)(sc*(i)*xMax/Ncx)+gx;
        iy = (int)(sc*(j)*yMax/Ncy)+gy;
        tm = 0.0;
        if (flown[i][j]>0.0) tm = flowt[i][j]/flown[i][j]/3000.0;
        if (tm>0.66) tm = 0.66;
        gOff.setColor(Color.getHSBColor((float)(0.67-tm),0.8f,0.8f));
        gOff.fillRect(ix,iy,(int)(sc*xMax/Ncx)+1,(int)(sc*yMax/Ncy)+1);
      }
    }
    setBalls(gx,gy);
  }

  void setBalls(int gx, int gy) {
    int i, ir,ix,iy;
    double sc,sz;

    sc = scale*dispScale;
    sz = 0.8*scale*dispScale;
    for (i=0; i<Nmt; i++) {
      if (kind[i]!=MAT0) {
        ix = (int)(sc*xx[i])+gx;
        iy = (int)(sc*yy[i])+gy;
        ir = (int)(sz*Morse[kind[i]][4] );
        gOff.setColor(Color.getHSBColor((float)(colorOf(i)),0.9f,0.9f));
        gOff.fillOval(ix-ir/2,iy-ir/2, ir, ir);
      }
    }
  }

  double colorOf(int i) {
    double col;

    col = 0.0;
    if (kind[i]==MAT0) {
      col = 0.33;
    } else if (kind[i]==MAT1) {
      col = 0.01;
    }  else if (kind[i]==MAT2) {
      col = 0.50;
    }  else if (kind[i]==MAT3) {
      col = 0.85;
    } else if (kind[i]==MAT4) {
      col = 0.66;
    }
    return ( col );
  }


/*-------------------------------------  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];
    }
    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]);
    }
    return ((double)Math.floor(s/((double)Nmt*1.38e-23)*10.0)/10.0 );
  }

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

    for (ix=0; ix<Ncx; ix++) {
      for (iy=0; iy<Ncy; iy++) {
        meanFlowx[ix][iy] = 0.0;
        meanFlowy[ix][iy] = 0.0;
        meanFlown[ix][iy] = 0.0;
        meanFlowt[ix][iy] = 0.0;
      }
    }
  }

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

    for (ix=0; ix<Ncx; ix++) {
      for (iy=0; iy<Ncy; iy++) {
        flowx[ix][iy] = 0.0;
        flowy[ix][iy] = 0.0;
        flown[ix][iy] = 0;
        flowt[ix][iy] = 0.0;
      }
    }

    for (i=0; i<Nmt; i++) {
      if (kind[i]==MAT0) {
        ix = (int)(xx[i]*Ncx/xMax);
        iy = (int)(yy[i]*Ncy/yMax);
        flowx[ix][iy] += vx[i];
        flowy[ix][iy] += vy[i];
        flown[ix][iy] += 1.0;
        flowt[ix][iy] += Morse[kind[i]][0]*(vx[i]*vx[i]+vy[i]*vy[i])/(2.0*1.38e-23);
      }
    }

    for (ix=0; ix<Ncx; ix++) {
      for (iy=0; iy<Ncy; iy++) {
        meanFlowx[ix][iy] = 0.99*meanFlowx[ix][iy]+0.01*flowx[ix][iy];
        meanFlowy[ix][iy] = 0.99*meanFlowy[ix][iy]+0.01*flowy[ix][iy];
        meanFlown[ix][iy] = 0.99*meanFlown[ix][iy]+0.01*flown[ix][iy];
        meanFlowt[ix][iy] = 0.99*meanFlowt[ix][iy]+0.01*flowt[ix][iy];
      }
    }
  }

  double ballPosX() {
    int i,n;
    double s;

    n = 0; s = 0.0;
    for (i=0; i<Nmt; i++) {
      if (kind[i]==MAT1) {
        n += 1;
        s += xx[i];
      }
    }
    if (n<1.0) n=1;
    return(s/n);
  }

  double ballPosY() {
    int i,n;
    double s;

    n = 0; s = 0.0;
    for (i=0; i<Nmt; i++) {
      if (kind[i]==MAT1) {
        n += 1;
        s += yy[i];
      }
    }
    if (n<1.0) n=1;
    return(s/n);
  }

  void traceBallInit(double posx, double posy) {
    int i;

    for (i=0; i<100; i++) {
      traceBallPos[i][0] = posx;
      traceBallPos[i][1] = posy;
    }
  }

  void traceBall(double posx, double posy) {
    int i;
	
    for (i=99; i>=1; i--) {
      traceBallPos[i][0] = traceBallPos[i-1][0];
      traceBallPos[i][1] = traceBallPos[i-1][1];
    }
    traceBallPos[0][0] = posx;
    traceBallPos[0][1] = posy;
  }

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

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

  void setForceTable(){
    int i,j,ir,k0,k1; 
    double d,a,r0,y,r,sc;
    int kkk[] = {0,1,2,3,4};

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

    kkk[0]=MAT0; kkk[1]=MAT1;kkk[2]=MAT2;kkk[3]=MAT3;kkk[4]=MAT4;
    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; 
    double s;

    setForceTable();

    i = 0;
    i= setBall(i,100.0*sgm,100.0*sgm, MAT1);
/*
    i = setLine(i,50.0*sgm,180.0*sgm,50.0*sgm,220.0*sgm,MAT4);
    i = setLine(i,150.0*sgm,180.0*sgm,150.0*sgm,220.0*sgm,MAT4);
    i = setLine(i,250.0*sgm,180.0*sgm,250.0*sgm,220.0*sgm,MAT4);
*/
    i = fillGas(i,MAT0);
    Nmt = i - 1;

    vAjustment();
    traceBallInit(ballPosX(), ballPosY());
  }


  int setBall(int icount,double x0,double y0,int knd) {
    int ic,ix,iy,lx,ly;
    double ax,ay,x,y,ra,ra2,rb,rb2,r2;
	
    lx = 50; ly = 80;
	
    ic = icount;

    ra = 0.5*(Morse[MAT4][4]+Morse[knd][4]);
    ra2 = ra*ra;
    rb = 25.0*Morse[knd][4];
    rb2 = rb*rb;
    ax = 1.0*Morse[knd][4]; ay = 0.866*ax;
    for (ix=0; ix<lx; ix++) {
      for (iy=0; iy<ly; iy++) {
        x = x0 + ax*(ix-lx/2)+0.5*ax*(iy%2);
        y = y0 + ay*(iy-ly/2);
        r2 = (x-x0)*(x-x0)+(y-y0)*(y-y0);
        if( r2<rb2 ) {
          xx[ic] = x;
          yy[ic] = y;
          kind[ic] = knd;
          vx[ic] = 400.0*(y-y0)/rb+100.0;
          vy[ic] = 400.0*(-(x-x0)/rb);
          ffx[ic] = 0.0;
          ffy[ic] = 0.0;
          ic += 1;
        }
      }
    }
    return(ic);
  }

  int setLine(int icount,double x0,double y0,double x1,double y1,int knd) {
    int ic,i,n;
    double len,x,y;
	
    ic = icount;
    len = Math.sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0));
    n = (int)(len/(1.2*sgm)+1.0);
    for (i=0; i<=n; i++) {
      x = x0 + i*(x1-x0)/n;
      y = y0 + i*(y1-y0)/n;
      xx[ic] = x;
      yy[ic] = y;
      kind[ic] = knd;
      vx[ic] = 0.0;
      vy[ic] = 0.0;
      ffx[ic] = 0.0;
      ffy[ic] = 0.0;
      ic += 1;
    }
    return(ic);
  }

  int fillGas(int icount,int knd) {
    int ic,ix,iy;
    double s,x,y,s2;

    ic = icount;
    s = xMax/150.0; s2 = 2.0*sgm;
    for (ix=0; ix<150; ix++) {
      for (iy=0; iy<100; iy++) {
        x = s*(double)(ix)+0.5*s;
        y = s*(double)(iy)+0.5*s;
        if ( objectNotFound(icount,x,y,s2)==-1 ) {
          kind[ic] = knd;
          xx[ic] = x;
          yy[ic] = y;
          vx[ic] = 500.0*Math.random()-250.0;
          vy[ic] = 500.0*Math.random()-250.0;
          ffx[ic] = 0.0;
          ffy[ic] = 0.0;
          ic += 1;
        }
      }
    }
    return(ic);
  }

  int objectNotFound(int icount,double x0,double y0,double r0) {
    int i,rt;
    double x,y;
	
    rt = -1;
    for (i=0; i<icount; i++) {
      x = xx[i]; y=yy[i];
      if (x>x0-r0 && x<x0+r0 && y>y0-r0 && y<y0+r0) {
        rt=0;
        break;
      }
    }
    return(rt);
  }

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

  void registration() {
    int i,j,ii,jj,ip,jp,kp,i0,i1,j0,j1,iq;
    double r2,rreg,rreg2;
    double xij,yij;
    double x,y;

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

      for(i=i0;i<=i1;i++) {
	    ii = (i+Nsx)%Nsx;
        for(j=j0;j<=j1;j++) {
          jj = (j+Nsy)%Nsy;
          for(iq=1;iq<=section[ii][jj][0];iq++) {
            jp = section[ii][jj][iq];
            if (jp>ip) {
              xij = xx[ip] - xx[jp]; yij = yy[ip] - yy[jp];
              if (xij>0.5*xMax) xij -= xMax;
              if (xij<-0.5*xMax) xij += xMax;  
              if (yij>0.5*yMax) yij -= yMax;
              if (yij<-0.5*yMax) yij += yMax; 
              r2=xij*xij+yij*yij;
              if (r2<rreg2) {
                reg[ip][kp]=jp;
                kp += 1;
              }
            }
          }
        }
      }
      reg[ip][0]=kp;

    }
  }

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


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

  void calcposition(){
    int i;

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

    if (resetFlag==1) {
      setInitialPosition();
      resetCount -= 1;
      if (resetCount<=0) resetFlag = 0;
    }

    traceCount += 1;
    if (traceCount>=100) {
      traceCount = 0;
      traceBall(ballPosX(), ballPosY());
    }

    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];
      xx[i] += vx[i]*dt;
      yy[i] += vy[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];
    }

    for (i=0; i<Nmt; i++) {
      if (xx[i] < 0.0) { 
        xx[i] += xMax; 
      }
      if (xx[i] > xMax) {
        xx[i] -= xMax; 
      }
      if (yy[i] < 0.0) { 
        yy[i] += yMax; 
      }
      if (yy[i] > yMax) {
        yy[i] -= yMax;
      }
    }

  }

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

    for(i=0;i<Nmt;i++) {
      ffx[i] = 0.0;
      ffy[i] = 0.0;
    }

    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];
        if (xij>0.5*xMax) xij -= xMax;
        if (xij<-0.5*xMax) xij += xMax;  
        if (yij>0.5*yMax) yij -= yMax;
        if (yij<-0.5*yMax) yij += yMax; 
        r2=xij*xij+yij*yij;
        if (r2<rc2) {
          rij = Math.sqrt(r2);
          f = 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;
        }
      }
    }
  }

  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 applet -----*/