/** applet No. 1074
 *
 * template - dynamics 2D
 * - multi thread  D2D :: display - asynchronous
 *
 * Created by Ikeuchi Mitsuru on December 10 2006.
 * Copyright (c) 2006-2007 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2006.12.10  created
 * ver 0.0.2   2006.12.13  improved code
 * ver 0.0.3   2007.06.08  improved code
 *
 */

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

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

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

  Thread th = null; // for run()-paint() thread
  Dynamics2D d2d = null;
  DisplayD2D dsp = new DisplayD2D();

  // for event
  Choice ch_mat,ch_view;
  Button bt_reset, bt_startStop, bt_step;
  Scrollbar sc_gravity, sc_spin, sc_scale;

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

  int dgX, dgY, dgXb, dgYb;

  int sleepTime = 30;
  int dispMode = 2;
  int thCount = 0;


  // ------------------------------ applet main thread -----------

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

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

    ch_mat = new Choice();
    ch_mat.add("fixed");
    ch_mat.add("A");
    ch_mat.add("B");
    ch_mat.add("C");
    ch_mat.addItemListener(this);
    ch_mat.select("B");

    ch_view = new Choice();
    ch_view.add("ball");
    ch_view.add("bond"); 
    ch_view.add("ball+bond");
    ch_view.add("velocity");
    ch_view.add("force");
    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);

    bt_step= new Button("step");
    bt_step.addActionListener(this);

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

    sc_spin= new Scrollbar(Scrollbar.HORIZONTAL,300,10,0,410);
    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");

  }

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

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


  // ---------------------------------- event listener -----------

  public void itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_view){
      dispMode = ch_view.getSelectedIndex();
    }
  }

  public void actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_reset){ 
      d2d.reset();
      thCount = 0;
    } else if (ev.getSource() == bt_startStop){
      if (d2d.getStartSW()==0) {
        d2d.setStartSW(1);
      } else {
        d2d.setStartSW(0);
      }
    } else if (ev.getSource() == bt_step){
      if (d2d.getStartSW()==0) {
        d2d.setStepSW();
      }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_gravity) { 
      d2d.setGravity( 0.1*(double)(sc_gravity.getValue()) );
    } else if (ev.getSource() == sc_spin) { 
      d2d.setSpin( 0.01*(double)(sc_spin.getValue()) );
    } else if (ev.getSource() == sc_scale) { 
      dsp.setScale( 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;
    }
    dsp.setViewDelta( (dgX-dgXb), (dgY-dgYb) );
  }



  // ========================= run() - paint() loop ==============

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

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

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

  private void offPaint() {
    if (d2d != null) {
      dsp.plotGraph(gOff, d2d, dispMode, thCount);
    }
  }

}  



// ================================= display D2D class ===========

class DisplayD2D {

  Graphics gOff;
  Dynamics2D d2d;

  private double xMax = Dynamics2D.xMax;
  private double yMax = Dynamics2D.yMax;

  private int Nmt = 0;

  private double dispWidth = 400.0;
  private double dispScale = (dispWidth/xMax);
  private double viewScale = 1.0;

  private double xShift = 0.0;
  private double yShift = 0.0;

  private int drawnCount = 0;


  // ------------------------------- class constructor -----------

  DisplayD2D() {
  }


  // ----------------------------------- class methods -----------

  public void plotGraph(Graphics gg, Dynamics2D ad2d, int dispMode, int thCount) {
    int d2dCount, gx, gy, gbx;
    double sc;

    gOff = gg;
    d2d = ad2d;
    if (d2d==null) return;
    d2dCount = d2d.getCount();
    if (d2dCount==drawnCount) return;
    drawnCount = d2dCount;

    Nmt = d2d.getNmt();
    sc = dispScale*viewScale;
    gx = 30+(int)(sc*xShift);
    gy = 80+(int)(sc*yShift);

    gOff.setColor(Color.white);
    gOff.fillRect(0,0,630,420);


    drawWaku(gOff,gx,gy,sc, viewScale);
    if (dispMode==0) {
      ballPlot(gOff,gx,gy,sc, 0.8, 0.9);

    } else if (dispMode==1) {
      bondPlot(gOff,gx,gy,sc);

    } else if (dispMode==2) {
      ballPlot(gOff,gx,gy,sc, 0.4, 0.9);
      bondPlot(gOff,gx,gy,sc);

    } else if (dispMode==3) {
      ballPlot(gOff,gx,gy,sc, 0.4, 0.9);
      velocityPlot(gOff,gx,gy,sc,viewScale*20.0);

    } else if (dispMode==4) {
      ballPlot(gOff,gx,gy,sc, 0.4, 0.9);
      forcePlot(gOff,gx,gy,sc,viewScale*1.0);

    } else if (dispMode==5) {
      ;

    }

    gbx = 630/6;
    gOff.setColor(Color.black);
    gOff.drawString("t="+(int)(d2d.getTime()*10000.0+0.5)/10.0+" ms",gbx*0+10,40);
    if (d2d.getStartSW()==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)(d2d.getSpin()*100.0+0.5)/100.0+" m/s",gbx*2+10,40);
    gOff.drawString("gravity="+(int)(d2d.getGravity()/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);

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

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


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

  private void ballPlot(Graphics gOff, int gx, int gy, double sc, double dmag, double dens) {
    int i,ix,iy,ir;
    double tm;

    for (i=0; i<Nmt; i++) {
      ix = (int)(sc*d2d.getxx(i))+gx;
      iy = (int)(sc*d2d.getyy(i))+gy;
      ir = (int)(sc*dmag*d2d.getr0(i));
      gOff.setColor(Color.getHSBColor(d2d.getColorOf(i),(float)dens,0.8f));
      gOff.fillOval(ix-ir/2,iy-ir/2, ir, ir);
    }
  }

  private void bondPlot(Graphics gOff, int gx, int gy, double sc) {
    int i,j,k,ix,iy,ix2,iy2;
    double xi,xj,yi,yj,r0,r2,rbond2,col;

    for (i=0; i<Nmt; i++) {
      xi = d2d.getxx(i); yi = d2d.getyy(i);
      ix = (int)(sc*xi)+gx;
      iy = (int)(sc*yi)+gy;

      for(k=1; k<d2d.getReg(i,0); k++) {
        j = d2d.getReg(i,k);
		xi = d2d.getxx(i); yi = d2d.getyy(i);
        xj = d2d.getxx(j); yj = d2d.getyy(j);
        r2 = (xi-xj)*(xi-xj)+(yi-yj)*(yi-yj);
        r0 = 0.5*(d2d.getr0(i)+d2d.getr0(j));
        rbond2 = r2/(r0*r0);
        if (rbond2<2.0 && j>i) {
          if (d2d.bondQ(i,j)==-1) {
            col = 0.33-20.0*(rbond2-1.0);
            if (col<0.01) col = 0.01;
            if (col>0.99) col = 0.99;
            gOff.setColor(Color.getHSBColor((float)col,0.9f,0.9f));
            ix2 = (int)(sc*xj)+gx;
            iy2 = (int)(sc*yj)+gy;
            gOff.drawLine(ix,iy, ix2,iy2);
          } else if (rbond2<1.1) {
            gOff.setColor(Color.gray);
            ix2 = (int)(sc*xj)+gx;
            iy2 = (int)(sc*yj)+gy;
            gOff.drawLine(ix,iy, ix2,iy2);
          }
        }
      }
    }
  }

  private void velocityPlot(Graphics gOff, int gx, int gy, double sc, double vmag) {
    int i,ix,iy,ir,ix2,iy2;
    double vx;

    for (i=0; i<Nmt; i++) {
      ix = (int)(sc*d2d.getxx(i))+gx;
      iy = (int)(sc*d2d.getyy(i))+gy;
      vx = d2d.getvx(i);
      ix2 = ix+(int)(vmag*vx);
      iy2 = iy+(int)(vmag*d2d.getvy(i));
      if (vx>0) {
        gOff.setColor(Color.blue);
      } else {
        gOff.setColor(Color.red);
      }
      gOff.drawLine(ix,iy, ix2,iy2);
    }

  }

  private void forcePlot(Graphics gOff, int gx, int gy, double sc, double fmag) {
    int i,ix,iy,ir,ix2,iy2;
    double ffx;

    for (i=0; i<Nmt; i++) {
      ix = (int)(sc*d2d.getxx(i))+gx;
      iy = (int)(sc*d2d.getyy(i))+gy;
	  ffx = d2d.getffx(i);
      ix2 = ix+(int)(fmag*ffx);
      iy2 = iy+(int)(fmag*d2d.getffy(i));

      if (ffx>0) {
        gOff.setColor(Color.blue);
      } else {
        gOff.setColor(Color.red);
      }

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

  private void drawWaku(Graphics gOff, int gx, int gy, double sc, double viewScale) {
    int ix,iy,gball,idx,idy;

    gball = (int)(3.0*viewScale);
    gOff.setColor(Color.black); 
    ix = (int)(sc*xMax);
    iy = (int)(sc*yMax);
    idx = ix+2*gball;
    idy = iy+2*gball;
    gOff.drawRect(gx-gball,gy-gball,idx,idy);
  }


  // ------------------------------------- I/O methods -----------

  public void setViewDelta(double xDisplacement, double yDisplacement) {
    xShift += (1.0/viewScale/dispScale)*xDisplacement;
    yShift += (1.0/viewScale/dispScale)*yDisplacement;
  }

  public void setScale(double sc) { viewScale = sc; }
}



// ================================= dynamics 2D class ===========

class Dynamics2D extends Thread {

  private double t = 0.0;
  private double dt = 2.5e-5;
  public static final double xMax = 40.0*0.01; // m
  public static final double yMax = 30.0*0.01; // m

  private int Nmt = 10000;
  private int Nsx = 40;
  private int Nsy = 30;

  private double gravity = 9.8; // m/s2
  private double spin = 3.0;
  private int incRate = 20;
  private double vx0 = 4.0;

  private double rc = 0.01; // cutoff=1cm
  private double rc2 = rc*rc;
  private double hh = rc/1000.0; // delta x

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

  private String materialStr[] = {
    "fixed" ,"A" ,"B","C"
  };
  private int NmatKind = 4;

  private int NN = Nmt + 1;
  private int kind[] = new int[NN];
  private double xx[] = new double[NN];
  private double yy[] = new double[NN];
  private double vx[] = new double[NN];
  private double vy[] = new double[NN];
  private double ffx[] = new double[NN];
  private double ffy[] = new double[NN];

  private int bond[][] = new int[NN][10];
  private int reg[][] = new int[NN][100];
  private int section[][][] = new int[Nsx][Nsy][30];

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

  private int resetSW = 0;
  private int startSW = 1;
  private int stepSW = 0;
  private int goSW = 0;
  private int sleepTime = 1;
  private int stopSleepTime = 100;
  private int count = 0;


  // ------------------------------- class constructor -----------

  Dynamics2D() {
    setInitialCondition();
  }


  // -------------------------------------- thread run -----------

  public void run() {
    while (true) {
      count += 1;
      timeEvolution();
      if (sleepTime>0) {
        try {
          Thread.sleep(sleepTime);
        } catch (InterruptedException e) { }
      }
    }
  }


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

  private void setInitialCondition() {
    int i; 
    double s;

    setForceTable();

    t = 0.0;
    count = 0;
    goSW = 0;

    i = 0;
    i = setBall(i, 14.0*0.01, yMax-20.0*0.01,2);
    i = setFloor(i, 7, 1);
    Nmt = i;

    setBond();
  }

  private int setBall(int ii, double xpos, double ypos,int mat) {
    int i,ix,iy,iz,ir2;
    double a,r2,r02,x,y;

    i = ii;

    a = material[mat][4];
	r02 = 55.0*a*a;
    for (ix=-8; ix<=8; ix++) {
      for (iy=-8; iy<=8; iy++) {
        kind[i] = mat;
        x = a*ix+0.5*a*((iy+10)%2);
        y = 0.866*a*iy;
        r2 = x*x+y*y;
        if (r2<r02) {
          kind[i] = mat;
          xx[i] = xpos+x;
          yy[i] = ypos+y;
          vx[i] = -spin*y/0.08+1.0;
          vy[i] = spin*x/0.08;
          ffx[i] = 0.0;
          ffy[i] = 0.0;

          i += 1;
        }
      }
    }
    return ( i );
  }

  private int setFloor(int ii, int layer,int mat) {
    int i,ix,iy,iz,ir2,nx;
    double a,aa;

    i = ii;

    a = material[mat][4];
    nx = (int)(xMax/a);
    aa = xMax/nx;
    for (ix=0; ix<nx; ix++) {
      for (iy=0; iy<layer; iy++) {
        kind[i] = mat;
        xx[i] = 0.25*aa+aa*ix+0.5*aa*(iy%2);
        yy[i] = yMax-0.866*a*iy;
        vx[i] = 0.0;
        vy[i] = 0.0;
        ffx[i] = 0.0;
        ffy[i] = 0.0;

        i += 1;
      }
    }
    return ( i );
  }

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

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

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

  //

  void gogogo(int knd) {
    int i;

    if (goSW==1) {
      for(i=0;i<Nmt;i++) {
        if (kind[i]==knd) {
          vx[i] += vx0;
        }
      }
      goSW = 2;
    }
  }


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

  private void timeEvolution(){

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

    if (startSW==1) {
      timeStep();

    } else {
      if (stepSW==1) {
        timeStep();
        stepSW = 0;
      }
      try {
        Thread.sleep(stopSleepTime);
      } catch (InterruptedException e) { }
    }
  }

  private void timeStep() {
    int i;

    if (goSW==1) {
	  gogogo(2);
	}
    registration();
    for (i=0; i<incRate; i++) {
      timeStep1();
    }
  }

  private void timeStep1() {
    int i; 
    double dtv2, a2;

    t += dt;

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

    setBoundary();
  }

  private 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] = gravity*material[kind[i]][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];  
        r2 = (xi-xj)*(xi-xj)+(yi-yj)*(yi-yj);
        if (r2<rc2) {
          rij = Math.sqrt(r2);
          f = force(rij,i,j);
          fxij = f*(xi-xj)/rij;
          fyij = f*(yi-yj)/rij;

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

          ffx[j] += -fxij;
          ffy[j] += -fyij;
        } else {
          if (bondQ(i,j)==-1) unBonding(i,j);
        }
      }
    }
  }

  private 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) {
      // damp
      //ret += -0.5*(material[kind[ii]][5]+material[kind[ii]][5])*((xx[ii]-xx[jj])*(vx[ii]-vx[jj])+(yy[ii]-yy[jj])*(vy[ii]-vy[jj]));
      return ret;
    }
    return( 0.01*ret );
  }

  private void setBoundary() {
    int i; 
    double rr;

    rr = 0.9;
    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];
      }
    }
  }

  // registration

  private void registration() {
    int i,j,ii,jj,ip,jp,kp,i0,i1,j0,j1,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); if (i0<0) i0 = 0;
      i1 = (int)(Nsx*(xx[ip]+rreg)/xMax); if (i1>=Nsx) i1 = Nsx-1;
      j0 = (int)(Nsy*(yy[ip]-rreg)/yMax); if (j0<0) j0 = 0;
      j1 = (int)(Nsy*(yy[ip]+rreg)/yMax); if (j1>=Nsy) j1 = Nsy-1;

      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) {
              xi = xx[ip]; xj = xx[jp]; yi = yy[ip]; yj = yy[jp];
              r2=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj);
              if (r2<rreg2) {
                reg[ip][kp]=jp;
                kp += 1;
              }
            }
          }
        }
      }
      reg[ip][0]=kp;
    }

  }

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


  // ----------------------------------------- utility -----------

  // bond manager

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

    registration();

    for (i=0; i<Nmt; i++) {
      for (ip=0; ip<10; 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 = distance(i,j);
        if (r<1.2*r0) bonding(i,j);
      }
    }
  }

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

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

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

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

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

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

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

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

  //

  public int regMax() {
    int i,m;

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

  public int secMax() {
    int i,j,k,m;

    m = 0;
    for(i=0;i<Nsx;i++) {
      for(j=0;j<Nsy;j++) {
        if (section[i][j][0]>m) {
          m = section[i][j][0];
        }
      }
    }
    return(m);
  }


  // ------------------------------------- I/O methods -----------

  // control
  public void reset() { resetSW = 1; }
  public void setStartSW(int sw) { startSW = sw; }
  public int getStartSW() { return startSW; }
  public void setStepSW() { stepSW = 1; }
  public int getCount() { return count; }

  public void setGravity(double g) { gravity = g; }
  public double getGravity() { return gravity; }
  public void setSpin(double s) { spin = s; }
  public double getSpin() { return spin; }
  public void setvx0(double v) { vx0 = v; resetSW = 1; }
  public double getvx0() { return vx0; }
  public void setGo() { if (goSW==0) { goSW = 1; } }

  // get data
  public int getNmt() { return Nmt; }
  public double getTime() { return t; }

  public double getxx(int i) { return xx[i]; }
  public double getyy(int i) { return yy[i]; }
  public double getvx(int i) { return vx[i]; }
  public double getvy(int i) { return vy[i]; }
  public double getffx(int i) { return ffx[i]; }
  public double getffy(int i) { return ffy[i]; }
  public double getr0(int i) { return material[kind[i]][4]; }
  public float getColorOf(int i) { return (float)(material[kind[i]][6]); }

  public int getReg(int i,int k) { return reg[i][k]; }
}