/*
 * rolling cone - periodic dynamics 3D
 *
 * coded by Ikeuchi Mitsuru
 * ver 0.0.1   2005.12.10  created
 * ver 0.0.2   2006.01.29  change sort method (buble -> quick)
 *
 */

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

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

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

  Choice ch_mat,ch_view;
  Button bt_reset, bt_startStop;
  Scrollbar sc_gravity, sc_spin, sc_scale;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 50;

  double t = 0.0;
  double dt = 2.5e-5;
  double xMax = 20.0*0.01+0.004;
  double yMax = 20.0*0.01;
  double zMax = 20.0*0.01+0.004;
  double zDepth = xMax;

  int Nmt = 10000;
  int NN = Nmt + 1;

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

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

  int resetSW = 0;
  int started = 1;
  int dispMode = 2;
  double gravity = 9.8;
  double spin = 2.0;
  int MAT1 = 2;

  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 = 0.01;
  double rc2 = rc*rc;
  double hh = rc/1000.0;


  double AU = 1.0e-3;
  double EE = 0.01;
  double AA = 0.001;
  double material[][] = {
  /*       0:mass(AU) 1:charge 2:D(EE) 3:A(/AA) 4:r0(AA) 5:color */
  /*fixed*/ { 10.0 *AU, 0.0*EE, 0.1*EE, 2.0/AA, 8.0*AA, 0.85 },
  /* 1 A */ {  5.0 *AU, 0.0*EE, 1.0*EE, 2.0/AA, 8.0*AA, 0.05 },
  /* 2 B */ {  5.0 *AU, 0.0*EE, 5.0*EE, 2.0/AA, 8.0*AA, 0.15 },
  /* 3 C */ { 10.0 *AU, 0.0*EE, 8.0*EE, 2.0/AA, 8.0*AA, 0.25 }
  };

  String materialStr[] = {
    "fixed" ,"A" ,"B","C"
  };
  int NmatKind = 4;
  
  int kind[] = new int[NN];
  double xx[] = new double[NN];
  double yy[] = new double[NN];
  double zz[] = new double[NN];
  double vx[] = new double[NN];
  double vy[] = new double[NN];
  double vz[] = new double[NN];
  double ffx[] = new double[NN];
  double ffy[] = new double[NN];
  double ffz[] = new double[NN];

  int bond[][] = new int[NN][20];
  int reg[][] = new int[NN][100];
  int section[][][][] = new int[Nsx][Nsy][Nsz][30];

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

  double forceTable[][][] = new double[NmatKind][NmatKind][1020];
  double dispColor[] = new double[NmatKind];

  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 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("fixed");
    ch_mat.addItem("A");
    ch_mat.addItem("B");
    ch_mat.addItem("C");
    ch_mat.addItemListener(this);
    ch_mat.select("B");

    ch_view = new Choice();
    ch_view.addItem("ball");
    ch_view.addItem("bond"); 
    ch_view.addItem("ball+bond");
    ch_view.addItem("velocity");
    ch_view.addItemListener(this);
    ch_view.select("ball+bond");

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

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

    sc_gravity= new Scrollbar(Scrollbar.HORIZONTAL,98,10,0,510);
    sc_gravity.addAdjustmentListener(this);

    sc_spin= new Scrollbar(Scrollbar.HORIZONTAL,200,10,0,310);
    sc_spin.addAdjustmentListener(this);

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

    addMouseListener(this);  
    addMouseMotionListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
    /* pnl.add(new Label("  ")); */
    pnl.add(bt_reset);
    pnl.add(bt_startStop);
    pnl.add(sc_spin);
    pnl.add(sc_gravity);
    pnl.add(sc_scale);
    pnl.add(ch_view);
    add(pnl,"North");

    setInitialCondition();
  }

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

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_gravity) { 
      gravity = 0.1*(double)(sc_gravity.getValue());
    } else if (ev.getSource() == sc_spin) { 
      spin = 0.01*(double)(sc_spin.getValue());
    } else if (ev.getSource() == sc_scale) { 
      viewScale = 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 {
        timeEvolution();
        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 gbx,gby,gx;

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

    rotate(theta,fai);
    zSort();

    gbx = 70;
    gby = 90;


    if (dispMode==0) {
       ballPlot(gbx,gby, 1.2, 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) {
      velocity3DPlot(gbx,gby, 20.0,0);

    } else if (dispMode==4) {
      ;

    } else if (dispMode==5) {
      ;

    }

/*
    potentialPlot(MAT1, MAT1);
*/
    gbx = dim.width/6;
    gOff.setColor(Color.black);
    gOff.drawString("t="+(int)(t*10000.0+0.5)/10.0+" ms",gbx*0+10,40);
    if (started==1) {
      gOff.drawString("started",630/6*1+10,40);
    } else {
      gOff.drawString("stop",630/6*1+10,40);
    }
    gOff.setColor(Color.black);
    gOff.drawString("spin0="+(int)(spin*100.0+0.5)/100.0+" m/s",gbx*2+10,40);
    gOff.drawString("gravity="+(int)(gravity/9.8*100.0+0.5)/100.0+" g",gbx*3+10,40);
    gOff.drawString("scale="+(int)(viewScale*100.0+0.5)+"%",gbx*4+10,40);
    gOff.drawString("view",gbx*5+10,40);

    gx = 450;
    gOff.setColor(Color.black);
    gOff.drawString("Box="+(int)(xMax*1000.0+0.5)/10.0+" x "+(int)(yMax*1000.0+0.5)/10.0+" x "+(int)(zMax*1000.0+0.5)/10.0+" cm",gx,60);

    gOff.drawString("regMax="+regMax()+"",530,80);
    gOff.drawString("secMax="+secMax()+"",530,100);

  }

/*-------------------------------  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,r,r0,col;

    drawWaku(gbx,gby,0);

    scl = dispScale*viewScale;
    sz = size*viewScale/AA;

    for (i=0; i<Nmt; i++) {
      j = srtz[0][i];
      gx = (int)(scl*px[j])+gbx;
      gy = (int)(scl*py[j])+gby;
      gr = (int)(sz*material[kind[j]][4]*(0.3*(pz[j]/zDepth)+0.7));
      if (gr<2) gr = 2;
      if (imode==1 || imode==2) {
        for (k=1; k<reg[j][0]; k++) {
          jj = reg[j][k];
          r = distance(j,jj);
          r0 = 0.5*(material[kind[j]][4]+material[kind[jj]][4]);
          if ( r<1.2*r0 ) {
            g2x = (int)(scl*px[jj])+gbx;
            g2y = (int)(scl*py[jj])+gby;
            if (bondQ(j,jj)==-1) {
              col = 0.33-50.0*(r/r0-1.0);if (col<0.01) col = 0.01;
              if (col>0.66) col = 0.66;
              gOff.setColor(Color.getHSBColor((float)(col),0.5f,0.8f));
            } else {
              gOff.setColor(Color.gray);
            }
            gOff.drawLine(gx,gy, g2x, g2y);
          }
        }
      }
      if (imode==0 || imode==2) {
        gOff.setColor(Color.getHSBColor(
          (float)material[kind[j]][5],0.8f,(float)(0.4*pz[j]/zDepth+0.3)));
        gOff.fillOval(gx-gr/2,gy-gr/2, gr, gr);
      }
    }
    drawWaku(gbx,gby,1);
  }

  void velocity3DPlot(int gbx, int gby, double mag, int imode) {
    int i,j,gx,gy,g2x,g2y,gr;
    double sc,sz,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;
	sz = 0.3*viewScale/AA;

    for (i=0; i<Nmt; i++) {
      j = srtz[0][i];
      gx = (int)(sc*px[j])+gbx;
      gy = (int)(sc*py[j])+gby;
      gr = (int)(sz*material[kind[j]][4]*(0.3*(pz[j]/zDepth)+0.7));

      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*material[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;
      if (vx[j]>0) {
        gOff.setColor(Color.blue);
      } else {
        gOff.setColor(Color.red);
      }
      gOff.drawLine(gx,gy, g2x,g2y);
	  
      gOff.setColor(Color.getHSBColor(
        (float)material[kind[j]][5],0.8f,(float)(0.4*pz[j]/zDepth+0.3)));
      gOff.fillOval(gx-gr/2,gy-gr/2, gr, gr);
    }

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

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

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

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

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

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

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

/*-------------------- potential / plot  ----------*/

  void potentialPlot(int ikind, int jkind) {
    int i,gx,gy;
    double v,a,d,r0,r;

    d = Math.sqrt(material[ikind][2]*material[jkind][2]);
    a = 0.5*(material[ikind][3]+material[jkind][3]);
    r0 = 0.5*(material[ikind][4]+material[jkind][4]);

    gx = 400; gy = 300;
    gOff.setColor(Color.black);
    gOff.drawLine(gx,gy,gx+100,gy);

    for (i=0; i<100; i++) {
      r = i*0.0001;
      v = 1000.0*(d*Math.exp(-2.0*a*(r-r0))-2.0*d*Math.exp(-a*(r-r0)));
      if (v>=0.0) {
        if(v<200.0) {
          gOff.setColor(Color.green);
          gOff.drawLine(gx+i,gy-(int)(v),gx+i,gy);
        }
      } else {
        gOff.setColor(Color.pink);
        gOff.drawLine(gx+i,gy,gx+i,gy-(int)(v));
      }
    }
  }

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

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

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

  void setInitialCondition() {
    int i; 
    double s;

    setForceTable();
    setWaku();

    t = 0.0;

    i = 0;
    i = setCone(i, 0.5*xMax, 0.65*yMax, 0.5*zMax,2);
    i = setFloor(i, 2, 3);
    Nmt = i;

    setBond();
    zSort();
  }

  int setFloor(int ii, int mat, int layers) {
    int i,ix,iy,iz,nx,nz;
    double a,ax,az;

    i = ii;

    a = material[mat][4]*0.707;
    nx = 2*((int)(xMax/2.0/a+0.5)); ax = xMax/nx;
    nz = 2*((int)(zMax/2.0/a+0.5)); az = zMax/nz;

    for (ix=0; ix<nx; ix++) {
      for (iy=0; iy<layers; iy++) {
        for (iz=0; iz<nz; iz++) {
          if ((ix+iy+iz)%2==1) {
              kind[i] = mat;
              xx[i] = 0.5*a+ax*(ix);
              yy[i] = yMax-a*(iy+0.5);
              zz[i] = 0.5*a+az*(iz);
		  
              vx[i] = 0.01*(Math.random()-0.5);
              vy[i] = 0.01*(Math.random()-0.5);
              vz[i] = 0.01*(Math.random()-0.5);
              ffx[i] = 0.0;
              ffy[i] = 0.0;
              ffz[i] = 0.0;
              i += 1;
          }
        }
      }
    }
    return ( i );
  }

  int setBall(int ii, double xpos, double ypos, double zpos,int mat,int iradius) {
    int i,ix,iy,iz,ir2;
    double a,r;

    i = ii;

    a = material[mat][4]*0.707;
	
    for (ix=-iradius; ix<iradius; ix++) {
      for (iy=-iradius; iy<iradius; iy++) {
        for (iz=-iradius; iz<iradius; iz++) {
          if ((ix+iy+iz+100)%2==1) {
            ir2 = ix*ix+iy*iy+iz*iz;
            if (ir2<iradius*iradius) {
              r = Math.sqrt((double)ir2);
              kind[i] = mat;
              xx[i] = xpos+a*ix;
              yy[i] = ypos+a*iy;
              zz[i] = zpos+a*iz;
		  
              vx[i] = -spin*iy/iradius+1.0;
              vy[i] = spin*ix/iradius;
              vz[i] = 0.01*(Math.random()-0.5);
              ffx[i] = 0.0;
              ffy[i] = 0.0;
              ffz[i] = 0.0;
              i += 1;
            }
          }
        }
      }
    }
    return ( i );
  }

  int setCone(int ii, double xpos, double ypos, double zpos,int mat) {
    int i,ix,iy,iz,ir2;
    double a,r;

    i = ii;

    a = material[mat][4]*0.707;
	
    for (ix=-8; ix<8; ix++) {
      for (iy=-8; iy<8; iy++) {
        for (iz=-8; iz<8; iz++) {
          if ((ix+iy+iz+100)%2==1) {
            ir2 = ix*ix+iy*iy;
            if (ir2<(8.0-0.4*(iz+8))*(8.0-0.4*(iz+8))) {
              r = Math.sqrt((double)ir2);
              kind[i] = mat;
              xx[i] = xpos+a*ix;
              yy[i] = ypos+a*iy;
              zz[i] = zpos+a*iz;
		  
              vx[i] = -spin*iy/8+1.0;
              vy[i] = spin*ix/8;
              vz[i] = 0.01*(Math.random()-0.5);
              ffx[i] = 0.0;
              ffy[i] = 0.0;
              ffz[i] = 0.0;
              i += 1;
            }
          }
        }
      }
    }
    return ( i );
  }

/*-------  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(material[i][2]*material[j][2]);
        a = 0.5*(material[i][3]+material[j][3]);
        r0 = 0.5*(material[i][4]+material[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 );
  }

/*------------  bond manager  ---------*/

  void setBond() {
    int i,j,ip;
    double r,r0;

    registration();

    for (i=0; i<Nmt; i++) {
      for (ip=0; ip<20; ip++) {
        bond[i][ip] = -1;
      }
    }

    for (i=0; i<Nmt; i++) {
      for(ip=1;ip<reg[i][0];ip++) {
        j = reg[i][ip];
        r0 = 0.5*(material[kind[i]][4]+material[kind[j]][4]);
        r = distancePreiodic(i,j);
        if (r<1.2*r0) bonding(i,j);
      }
    }
  }

  int bonding(int ii, int jj) {
    int i,ip;

    for (ip=0; ip<20; ip++) {
      if (bond[ii][ip]==-1) {
        bond[ii][ip] = jj;
        return(-1);
      }
    }
    return(0);
  }

  int unBonding(int ii, int jj) {
    int i,ip;

    for (ip=0; ip<20; ip++) {
      if (bond[ii][ip]==jj) {
        bond[ii][ip] = -1;
        return(-1);
      }
    }
    return(0);
  }

  int bondQ(int ii, int jj) {
    int i,ip;

    for (ip=0; ip<20; ip++) {
      if (bond[ii][ip]==jj) {
        return(-1);
      }
    }
    return(0);
  }

  double distance(int i, int j) {
    double xij,yij,zij;

    xij = xx[i]-xx[j]; yij = yy[i]-yy[j]; zij = zz[i]-zz[j];
    return(Math.sqrt(xij*xij+yij*yij+zij*zij));
  }

  double distancePreiodic(int i, int j) {
    double xi,xj,yi,yj,zi,zj;

    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 (zi-zj>0.5*zMax) { zj += zMax; };
    if (zi-zj<-0.5*zMax) { zj -= zMax; }; 
    return(Math.sqrt((xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)));
  }

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

  void timeEvolution() {
    int i;

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

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

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

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

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

    rr = 0.9;
    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] = 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] += 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] = gravity*material[kind[i]][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 (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,i,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;
        } else {
          if (bondQ(i,j)==-1) {
            unBonding(i,j);
          }
        }
      }
    }
  }

  double force(double r, int ii, int jj) {
    int i,ki,kj;
    double a,ret;

    ki = kind[ii]; kj = kind[jj];
    i = (int)(r/hh);
    a = r - i*hh;
    ret = ((hh-a)*forceTable[ki][kj][i] + a*forceTable[ki][kj][i+1])/hh;

    if (bondQ(ii,jj)==-1) return ( ret );
    return( 0.1*ret );
  }

/*-------------------------------  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 + 5.0*incRate*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); 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+(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 (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;
    }
  }

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

}  

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