/* many particles - short range Morse 2D */
/* coded by Ikeuchi Mitsuru              */
/* ver 0.0.1   2004.12.26                */

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

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

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

  Button bt_go;
  Choice cm1,ctc,cvw;
  Scrollbar scf,scr,scs;
  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 = 8.0/sgm;

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


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


  int tempMode = 0;
  double contTemp = 300.0;
  double scale = 0.12;
  int dispMode = 4;
  double setTemp = 10000.0;
  int started = 0;

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

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

  int MAT1 = 17;

  int kind[] = new int[NN];
  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    */
  double fbx[] = new double[NN];  /* i-th x-force before dt */
  double fby[] = new double[NN];  /* i-th y-force before dt */

  int reg[][] = new int[NN][30];   /* registration */

  double forceTable[] = new double[20000];

  double flowx[][] = new double[31][11];
  double flowy[][] = new double[31][11];
  double flown[][] = new double[31][11];
  double flowt[][] = new double[31][11];

  double meanFlowx[][] = new double[31][11];
  double meanFlowy[][] = new double[31][11];
  double meanFlown[][] = new double[31][11];
  double meanFlowt[][] = new double[31][11];

  double hh = 1.0e-13;
  double mas = 39.948*1.661e-27;
  double r0 = 3.032*AA;

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

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

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

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

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

    cvw = new Choice();
    cvw.addItem("ball"); cvw.addItem("line"); 
    cvw.addItem("both"); cvw.addItem("temp");
	cvw.addItem("coarse n"); cvw.addItem("coarse T");
    cvw.addItemListener(this);
    cvw.select("coarse n");

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

    scf= new Scrollbar(Scrollbar.HORIZONTAL,60,10,0,110);
    scf.addAdjustmentListener(this);

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

    scs= new Scrollbar(Scrollbar.HORIZONTAL,12,10,10,210);
    scs.addAdjustmentListener(this);

    addMouseListener(this);  
    addMouseMotionListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
    pnl.add(cm1);
    pnl.add(ctc);
    pnl.add(cvw);
    pnl.add(bt_go);
    pnl.add(scr);
    pnl.add(scs);
    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 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() == cm1) {
      MAT1 = 17+cm1.getSelectedIndex();
      t = 0.0;
      started = 0;
      setInitialPosition();
    } else if (ev.getSource() == ctc) {
      tempMode = ctc.getSelectedIndex();
    } else if (ev.getSource() == cvw) {
      dispMode = cvw.getSelectedIndex();
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == scr) { 
      contTemp = 10.0*(double)(scr.getValue());
      vAjustment();
    } else if (ev.getSource() == scs) { 
      scale = 0.01*(double)(scs.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,ir,j,k,ix2,iy2;
    int gball,gbx,gby,gb;
    double sc,sz,tmp,tm,cd,r2,r2c;

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

    setFlow();
    gball = (int)(scale*10);
    sc = scale*dispScale;
    sz = 0.9*scale*dispScale;
    gbx = 30+(int)(sc*xShift);
    gby = 80+(int)(sc*yShift);

    tmp = temperature();
    gOff.setColor(Color.gray); 
    if (tempMode==2) {
      gOff.setColor(Color.black); 
      if (tmp<contTemp) { gOff.setColor(Color.red); };
    }
    ix = (int)(sc*xMax);
    iy = (int)(sc*yMax);
    gOff.drawRect(gbx-gball/2,gby-gball/2,ix+gball,iy+gball);

    if (dispMode==0) {

      ir = (int)(sz*Morse[MAT1][4] );
      for (i=0; i<Nmt; i++) {
        ix = (int)(sc*xx[i])+gbx;
        iy = (int)(sc*yy[i])+gby;
        gOff.setColor(Color.green);
        if (kind[i]>0) {
          gOff.setColor(Color.green);
        } else {
          gOff.setColor(Color.gray);
        }
        gOff.fillOval(ix-ir/2,iy-ir/2, ir, ir);
      }

    } else if (dispMode==1) {

      r2c = 1.69*r0*r0;
      gOff.setColor(Color.gray);
      for (i=0; i<Nmt; i++) {
        ix = (int)(sc*xx[i])+gbx;
        iy = (int)(sc*yy[i])+gby;

        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<r2c && j>i) {
            ix2 = (int)(sc*xx[j])+gbx;
            iy2 = (int)(sc*yy[j])+gby;

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

    } else if (dispMode==2) {

      r2c = 1.69*r0*r0;
      ir = (int)(0.7*sz*Morse[MAT1][4] );
      for (i=0; i<Nmt; i++) {
        ix = (int)(sc*xx[i])+gbx;
        iy = (int)(sc*yy[i])+gby;

        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<r2c && j>i) {
            ix2 = (int)(sc*xx[j])+gbx;
            iy2 = (int)(sc*yy[j])+gby;

            gOff.setColor(Color.gray);
            gOff.drawLine(ix,iy, ix2,iy2);
          }
        }
        if (kind[i]>0) {
          gOff.setColor(Color.green);
        } else {
          gOff.setColor(Color.gray);
        }
        gOff.fillOval(ix-ir/2,iy-ir/2, ir, ir);
      }

    } else if (dispMode==3) {

      ir = (int)(sz*Morse[MAT1][4] );
      for (i=0; i<Nmt; i++) {
        ix = (int)(sc*xx[i])+gbx;
        iy = (int)(sc*yy[i])+gby;
        tm = 0.5*mas*(vx[i]*vx[i]+vy[i]*vy[i])/1.38e-23;
        tm = tm/3000.0+0.3;
        if (tm>0.99) {
          tm = 0.99;
        }
        if (kind[i]>0) {
          gOff.setColor(Color.getHSBColor(0.33f,0.95f,(float)(tm)));
        } else {
          gOff.setColor(Color.gray);
        }
        gOff.fillOval(ix-ir/2,iy-ir/2, ir, ir);
      }

    } else if (dispMode==4) {

      for (i=0; i<30; i++) {
        for (j=0; j<10; j++) {
          ix = (int)(sc*(i)*xMax/30.0)+gbx;
          iy = (int)(sc*(j)*yMax/10.0)+gby;
		  cd = 0.004*flown[i][j]; if (cd>0.99) cd = 0.99;
		  gOff.setColor(Color.getHSBColor(0.33f,(float)cd,0.9f));
          gOff.fillRect(ix,iy,(int)(sc*xMax/30.0)+1,(int)(sc*yMax/10.0)+1);
        }
      }
      ir = (int)(sz*Morse[MAT1][4] );
      for (i=0; i<Nmt; i++) {
        if (kind[i]<0) {
          ix = (int)(sc*xx[i])+gbx;
          iy = (int)(sc*yy[i])+gby;
          gOff.setColor(Color.gray);
          gOff.fillOval(ix-ir/2,iy-ir/2, ir, ir);
        }
      }

    } else if (dispMode==5) {

      for (i=0; i<30; i++) {
        for (j=0; j<10; j++) {
          ix = (int)(sc*(i)*xMax/30.0)+gbx;
          iy = (int)(sc*(j)*yMax/10.0)+gby;
          tm = flowt[i][j];
          if (flown[i][j]>0.0) {
            tm = tm/flown[i][j]/1500.0;
          } else {
            tm = 0.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/30.0)+1,(int)(sc*yMax/10.0)+1);
        }
      }
      ir = (int)(sz*Morse[MAT1][4] );
      for (i=0; i<Nmt; i++) {
        if (kind[i]<0) {
          ix = (int)(sc*xx[i])+gbx;
          iy = (int)(sc*yy[i])+gby;
          gOff.setColor(Color.gray);
          gOff.fillOval(ix-ir/2,iy-ir/2, ir, ir);
        }
      }

    }

    for (i=0; i<30; i++) {
      for (j=0; j<10; j++) {
        ix = (int)(sc*(i+0.5)*xMax/30.0)+gbx;
        iy = (int)(sc*(j+0.5)*yMax/10.0)+gby;
        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);
      }
    }

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

    gOff.setColor(Color.blue);
    gOff.drawString("t="+Math.floor(t*1.0e13+0.5)/10.0+"(ps)",630/6*3+10,60);
    gOff.drawString("T="+temperature()+" K",630/6*4+10,60);
    gOff.drawString("Lx="+(int)(xMax*1.0e11+0.5)/100.0+" nm",630/6*5+10,60);
/*
    gOff.drawString("regMax="+regMax()+" ",630/6*5+10,80);
*/
  }

/*-------------------------------------  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*mas*(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<31; ix++) {
      for (iy=0; iy<11; 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<31; ix++) {
      for (iy=0; iy<11; iy++) {
        flowx[ix][iy] = 0.0;
        flowy[ix][iy] = 0.0;
        flown[ix][iy] = 0.0;
        flowt[ix][iy] = 0.0;
      }
    }

    for (i=0; i<Nmt; i++) {
      if (kind[i]>0) {
        ix = (int)(xx[i]*30.0/xMax);
        iy = (int)(yy[i]*10.0/yMax);
        flowx[ix][iy] += vx[i];
        flowy[ix][iy] += vy[i];
        flown[ix][iy] += 1.0;
        flowt[ix][iy] += mas*(vx[i]*vx[i]+vy[i]*vy[i])/(3.0*1.38e-23);
      }
    }

    for (ix=0; ix<31; ix++) {
      for (iy=0; iy<11; 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];
      }
    }
  }

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

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

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

    x0 = 100*sgm; y0 = yMax/2.0;
    rr = 50*sgm;
    for (i=0; i<Nmt; i++) {
      x = xx[i] - x0;
      y = yy[i] - y0;
      if (x*x+y*y<rr*rr) {
        vx[i] = 4000.0*Math.random()-2000.0;
        vy[i] = 4000.0*Math.random()-2000.0;
      }
    }
  }

/*-------------------------------  set force table   ---------*/

  void setForceTable(int knd) {
    int ir; 
    double d,a,r0,y,r;

    d = Morse[knd][2];
    a = Morse[knd][3];
    r0 = Morse[knd][4];

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

    for (ir=0; ir<100; ir++) {
      forceTable[ir]=forceTable[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; 
    double s,x,y;

    mas = Morse[MAT1][0];
    r0 = Morse[MAT1][4];
    rc = 1.6*r0;
    rc2 = rc*rc;
    setForceTable(MAT1);

    s = xMax/300.0;
    i = 0;
    for (ix=0; ix<300; ix++) {
      for (iy=0; iy<100; iy++) {
        kind[i] = 1;
        x = s*(double)(ix)+0.5*s;
        y = s*(double)(iy)+0.5*s;
 
        xx[i] = x;
        yy[i] = y;
        vx[i] = 500.0*Math.random()-250.0;
        vy[i] = 500.0*Math.random()-250.0;
        ffx[i] = 0.0;
        ffy[i] = 0.0;
        fbx[i] = 0.0;
        fby[i] = 0.0;
        i += 1;

      }
    }

    Nmt = i - 1;
    vAjustment();
    setFlowInit();
  }

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

  void registration() {
    int i,j,k;
    double xi,xj,yi,yj;
    double r2,rreg,rreg2;

    rreg = rc+2000.0*20.0*dt;
    rreg2 = rreg*rreg;
    for(i=0;i<Nmt;i++) {
      k = 1;
      for(j=i+1;j<Nmt;j++) {
        xi = xx[i]; xj = xx[j]; yi = yy[i]; yj = yy[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; };
*/
        r2=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj);
        if (r2<rreg2) {
          reg[i][k]=j;
          k += 1;
        }
      }
      reg[i][0]=k;
    }
  }

  int regMax() {
    int i,m;

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

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

  void calcposition(){
    int i;

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

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

  void timeEvolution(){
    int i; 
    double rr, tmp;

    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++) {
	  if (kind[i]>0) {
        xx[i] += vx[i]*dt+0.5*dt*dt*ffx[i]/mas;
        yy[i] += vy[i]*dt+0.5*dt*dt*ffy[i]/mas;
      }
      fbx[i] = ffx[i];
      fby[i] = ffy[i];
    }
    forceCalc();
    for (i=0; i<Nmt; i++) {
	  if (kind[i]>0) {
        vx[i] += 0.5*dt*(ffx[i]+fbx[i])/mas;
        vy[i] += 0.5*dt*(ffy[i]+fby[i])/mas;
      }
    }

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

  }

  void forceCalc() {
    int i,j,k;
    double xi,xj,yi,yj;
    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];
        xi = xx[i]; xj = xx[j]; yi = yy[i]; yj = yy[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; };
*/
        r2=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj);
        if (r2<rc2) {
          rij = Math.sqrt(r2);
          f = force(rij);
          fxij = f*(xi-xj)/rij;
          fyij = f*(yi-yj)/rij;

          ffx[i] +=  fxij;
          ffy[i] +=  fyij;

          ffx[j] += -fxij;
          ffy[j] += -fyij;
        }
      }
    }

  }

  double force(double r) {
    int i;
    double a;

    i = (int)(r/hh);
    a = r - i*hh;

    return ( ((hh-a)*forceTable[i] + a*forceTable[i+1])/hh );
  }

}  

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