/* materials - dynamics 2D      */
/* coded by Ikeuchi Mitsuru     */
/* ver 0.0.1   2004.01.03       */

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

public class materialsD2D extends Applet
  implements ActionListener, AdjustmentListener, Runnable {

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

  Button bt_reset;
  Scrollbar scr,scf,sce;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 50;
  double um = 1.0;             /* unit of length (m)     */
  double dispScale = 10.0/um;

  double xMax = 40.0*um;       /* x-Box size in (m)      */
  double yMax = 20.0*um;       /* y-Box size in (m)      */

  int Nmt = 600;               /* number of molecules    */
  int NN = Nmt + 1;
  double t = 0.0;               /* time (s)               */
  double dt = 0.00025;           /* time division (s)      */

  double scale = 1.0;
  double gravity = 1.0;

  double rc = 1.2*um;
  double rc2 = rc*rc;

  double material[][] = {
   /* mm         ra          kc         elong  */
   { 4000.0,    1.0*um,    100.0e6,      0.07   },
   { 8000.0,    1.0*um,     60.0e6,      0.20   },
   { 2000.0,    1.0*um,     40.0e6,      0.15   }
  };
  int NNm = 3;


  double ra[][] = new double[NNm][NNm];
  double ra1[][] = new double[NNm][NNm];
  double ra2[][] = new double[NNm][NNm];
  double ram[][] = new double[NNm][NNm];
  double fa[][] = new double[NNm][NNm];
  double kc[][] = new double[NNm][NNm];

  double rcmp[][] = new double[NNm][NNm];
  double rten[][] = new double[NNm][NNm];


  int kind[] = new int[NN];       /* i-th molec kind */
  int bond[][] = new int[NN][NN]; /* i-th bond       */
  double mm[] = new double[NN];   /* i-th mass       */
  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][100];   /* registration */

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

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

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

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

    scr= new Scrollbar(Scrollbar.HORIZONTAL,10,10,0,110);
    add(new Label("g=")); add(scr);
    scr.addAdjustmentListener(this);

    scf= new Scrollbar(Scrollbar.HORIZONTAL,100,10,10,210);
    add(new Label("kConstant=")); add(scf);
    scf.addAdjustmentListener(this);

    sce= new Scrollbar(Scrollbar.HORIZONTAL,7,10,5,60);
    add(new Label("elongation=")); add(sce);
    sce.addAdjustmentListener(this);

    setSpringConst();
    setInitialPosition();
    repaint();
  }

  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_reset) {
      t = 0.0;

      setInitialPosition();

      repaint();
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if(ev.getSource() == scr) { 
      gravity = (double)(scr.getValue())/10.0;
      repaint();
    } else if(ev.getSource() == scf) { 
      material[0][2] = (1.0e6)*(double)(scf.getValue());
      setSpringConst();
      repaint();
    } else if(ev.getSource() == sce) { 
      material[0][3] = 0.01*(double)(sce.getValue());
      setSpringConst();
      repaint();
    }
  }

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

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

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

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

  void offPaint() {
    int i,j,ix,iy,ir,ixx,iyy,ki,kj;
    int gbx,gby,gb;
    double sc,d;

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

    gbx = 30;
    gby = 35;
    sc = scale*dispScale;
    ir = 8;

    ix = (int)(sc*xMax);
    iy = (int)(sc*yMax);
    gOff.setColor(Color.black);
    gOff.drawRect(gbx-ir/2,gby-ir/2,ix+ir,iy+ir);

    gOff.setColor(Color.lightGray);
    for (i=0; i<Nmt; i++) {
      ki = kind[i];
      ix = (int)(sc*xx[i])+gbx;
      iy = (int)(sc*yy[i])+gby;
      ir = 8;
      if (ki==0) {
        gOff.setColor(Color.getHSBColor( 0.66f,
          (float)(0.1+material[ki][2]/2.5e8),0.7f ));
      } else if (ki==1) {
        gOff.setColor(Color.getHSBColor( 0.50f,
          (float)(0.1+material[ki][2]/2.5e8),0.7f ));
      } else if (ki==2) {
        gOff.setColor(Color.getHSBColor( 0.33f,
          (float)(0.1+material[ki][2]/2.5e8),0.7f ));
      }
      gOff.fillOval(ix-ir/2,iy-ir/2, ir, ir);
    }

    for (i=0; i<Nmt; i++) {
      ki = kind[i];
      ix = (int)(sc*xx[i])+gbx;
      iy = (int)(sc*yy[i])+gby;
      for (j=i+1; j<Nmt; j++) {
        kj = kind[j];
        if (bond[i][j]>0) {
          d = Math.sqrt((xx[i]-xx[j])*(xx[i]-xx[j])+(yy[i]-yy[j])*(yy[i]-yy[j]));
          if (d<rcmp[ki][kj]) {
            gOff.setColor(Color.yellow);
          } else if (d<rten[ki][kj]) {
            gOff.setColor(Color.gray);
          } else if (d<ra1[ki][kj]) {
            gOff.setColor(Color.blue);
          } else {
            gOff.setColor(Color.red);
          }
          ixx = (int)(sc*xx[j])+gbx;
          iyy = (int)(sc*yy[j])+gby;
          gOff.drawLine(ix,iy,ixx,iyy);
        }
      }
    }

    gb = 460;
    gOff.setColor(Color.black);
    gOff.drawString("g="+gravity+"(g)",gb,50);
    gOff.drawString("Lx="+(int)(xMax+0.5)+"(m)",gb,70);
    gOff.drawString("Ly="+(int)(yMax+0.5)+"(m)",gb,90);
    gOff.setColor(Color.blue);
    gOff.drawString("t="+(int)(t*100.0)/100.0+"(s)",gb,130);

    gOff.setColor(Color.getHSBColor( 0.66f,
          (float)(0.1+material[0][2]/2.5e8),0.7f ));
    gOff.drawString("k="+(int)(material[0][2]/1.0e5+0.5)/10.0+"(MN/m)",gb,170);
    gOff.drawString("e="+(int)(material[0][3]*100+0.5)+"(%)",gb+100,170);

    gOff.setColor(Color.getHSBColor( 0.50f,
          (float)(0.1+material[1][2]/2.5e8),0.7f ));
    gOff.drawString("k="+(int)(material[1][2]/1.0e5+0.5)/10.0+"(MN/m)",gb,190);
    gOff.drawString("e="+(int)(material[1][3]*100+0.5)+"(%)",gb+100,190);

    gOff.setColor(Color.getHSBColor( 0.33f,
          (float)(0.1+material[2][2]/2.5e8),0.7f ));
    gOff.drawString("k="+(int)(material[2][2]/1.0e5+0.5)/10.0+"(MN/m)",gb,210);
    gOff.drawString("e="+(int)(material[2][3]*100+0.5)+"(%)",gb+100,210);

  }


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

  void setSpringConst() {
    int i,j;
    double rai,raj,ei,ej,rmax;

    rmax = 0.0;
    for (i=0; i<NNm; i++) {
      for (j=0; j<NNm; j++) {
        rai = material[i][1];
        raj = material[j][1];
        ei = material[i][3];
        ej = material[j][3];
        ra[i][j] = 0.5*(rai+raj);
        ra1[i][j] = 0.5*((1.0+0.5*ei)*rai+(1.0+0.5*ej)*raj);
        ra2[i][j] = 0.5*((1.0+ei)*rai+(1.0+ej)*raj);
        ram[i][j] = 0.5*(ra1[i][j]+ra2[i][j]);
        fa[i][j] = 2.0/((ra2[i][j]-ra1[i][j])*(ra2[i][j]-ra1[i][j]));
        rcmp[i][j] = (1.0-0.05*(ei+ej))*ra[i][j];
        rten[i][j] = (1.0+0.05*(ei+ej))*ra[i][j];
        if (ra2[i][j]>rmax) {
          rmax = ra2[i][j];
        }
        kc[i][j] = Math.sqrt(material[i][2]*material[j][2]);
      }
    }
    rc = 1.05*rmax;
    rc2 = rc*rc;

  }

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

  void setInitialPosition(){
    int i,ix,iy,mt; 
    double a,x,y;

    a = ra[kind[0]][kind[0]];

    i = 0;
    for (ix=0; ix<40; ix++) {
      for (iy=0; iy<40; iy++) {
        x = a*(ix+0.5*(iy%2));
        y = a*0.866*iy;

        mt = -1;
        if (x>5 && x<35 && y>6 && y<10) {
          mt = 0;
        } else if (x>10 && x<15 && y>10 && y<15) {
          mt = 1;
        } else if (x>25 && x<30 && y>10 && y<15) {
          mt = 1;
        }
        if (y>15.2 && y<20) {
          mt = 2;
        }

        if (mt>=0) {
          kind[i] = mt;
          xx[i] = x;
          yy[i] = y;
          mm[i] = material[kind[i]][0];
          vx[i] = 0.2*Math.random()-0.1;
          vy[i] = 0.2*Math.random()-0.1;

          ffx[i] = 0.0;
          ffy[i] = 0.0;
          fbx[i] = 0.0;
          fby[i] = 0.0;

          i += 1;
        }
      }
    }
    Nmt = i;
    bonding();

  }


  void bonding() {
    int i,j;
    double r;

    for (i=0; i<Nmt; i++) {
      for (j=i+1; j<Nmt; j++) {
        r = Math.sqrt((xx[i]-xx[j])*(xx[i]-xx[j])+(yy[i]-yy[j])*(yy[i]-yy[j]));
        if (r<ra2[kind[i]][kind[j]]) {
          bond[i][j]=1;
        }
      }
    }
  }

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

  void registration() {
    int i,j,k;
    double r2,rreg,rreg2;

    rreg = rc+1.0*20.0*dt;
    rreg2 = rreg*rreg;

    for(i=0;i<Nmt;i++) {
      k = 1;
      for(j=i+1;j<Nmt;j++) {
        r2=(xx[i]-xx[j])*(xx[i]-xx[j])+(yy[i]-yy[j])*(yy[i]-yy[j]);
        if (r2<rreg2) {
          reg[i][k]=j;
          k += 1;
        }
      }
      reg[i][0]=k;
    }
  }



/* ---------------------------------------------- */

  void calcposition(){
    int i;

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

      timeEvolution();

    }
  }


  void timeEvolution() {
    int i; 
    double rr;

    t = t + dt;

    rr = 0.95;

    for (i=0; i<Nmt; i++) {
      xx[i] += vx[i]*dt+0.5*dt*dt*ffx[i]/mm[i];
      yy[i] += vy[i]*dt+0.5*dt*dt*ffy[i]/mm[i];
      fbx[i] = ffx[i];
      fby[i] = ffy[i];
    }

    forceCalc();

    for (i=0; i<Nmt; i++) {
      vx[i] += 0.5*dt*(ffx[i]+fbx[i])/mm[i];
      vy[i] += 0.5*dt*(ffy[i]+fby[i])/mm[i];
    }

    for (i=0; i<Nmt; i++) {
      if (xx[i] < 0.0) { 
        xx[i] = 0.0; vx[i] = -rr*vx[i]; vy[i] = rr*vy[i]; 
      }
      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 r2,rij,f,fxij,fyij;

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

    for(i=0;i<Nmt;i++) {
      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<rc2) {
          rij = Math.sqrt(r2);
          f = force(rij,i,j);
          fxij = f*(xx[i]-xx[j])/rij;
          fyij = f*(yy[i]-yy[j])/rij;

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

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


  double force(double r,int ii, int jj) {
    int ki,kj;
    double f;

    ki = kind[ii];
    kj = kind[jj];
    if (r>ra2[ki][kj]) {
      bond[ii][jj]=0;
    } 

    f = 0.0;

    if (bond[ii][jj]==0) {
      if (r<ra[ki][kj]) {
        f = kc[ki][kj]*(r-ra[ki][kj]);
      }
    } else {
      if (r<=ra1[ki][kj]) {
        f = kc[ki][kj]*(r-ra[ki][kj]);
      } else if (r>=ra2[ki][kj]) {
        f = 0.0;
      } else if (ra1[ki][kj]<=r && r<ram[ki][kj]) {
        f = kc[ki][kj]*(r-ra[ki][kj])*(1.0-fa[ki][kj]*(r-ra1[ki][kj])*(r-ra1[ki][kj]));
      } else if (ram[ki][kj]<=r && r<ra2[ki][kj]) {
        f = kc[ki][kj]*(r-ra[ki][kj])*fa[ki][kj]*(r-ra2[ki][kj])*(r-ra2[ki][kj]);
      }
    }

    return ( -f );

  }


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

}  

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