/* evacuation and adsorption */
/* coded by Ikeuchi Mitsuru  */
/* ver 0.0.1   2002.07.28    */
/* ver 0.0.2   2003.08.09  flicker noise reduced */
/* ver 0.0.3   2004.06.26  improved pannel layout */

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

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

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

  Button bt_open;
  Scrollbar sad, scr;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 50;
  int atomSizeX = 40; 
  int atomSizeY = 20;  
  int dispR = 10;         /* display molecular radious */
  int paintSizeX = atomSizeX*dispR; 
  int paintSizeY = atomSizeY*dispR;
  int paintSizeX2 = 10*dispR; 
  int paintSizeY2 = 4*dispR;
  int wakuX = paintSizeX + dispR;
  int wakuY = paintSizeY + dispR;

  int Nmt = 10*10;                  /* number of molecules    */
  int NN = Nmt + 1;
  double t = 0.0;                /* time (s)               */
  double dt = 20.0*1.0e-15;      /* time division (s)      */
  double sgm = 3.40e-10;         /* L-J sigma for Ar (m)   */
  double eps = 1.67e-21;         /* L-J epsilon for Ar (J) */
  double mas = 39.95*1.67e-27;   /* mass of Ar (kg)        */

  double adsrbQ = 0.1;
  double contTemp = 300.0;
  double load = 300.0*1.0e-12; /* pN/atom */
  int gate = 1;

  double xMax = atomSizeX*sgm; /* x-Box size in (m)      */
  double yMax = atomSizeY*sgm; /* y-Box size in (m)      */
  double zMax = 0.0;           /* z-Box size in (m)      */

  int xBoxMin1 = 0;
  int xBoxMax1 = 20;
  int yBoxMin1 = 0;
  int yBoxMax1 = 20;

  int xBoxMin2 = xBoxMax1;
  int xBoxMax2 = xBoxMax1+20;
  int yBoxMin2 = 8;
  int yBoxMax2 = 12;

  int xDispMin1 = xBoxMin1*dispR;
  int xDispMax1 = xBoxMax1*dispR;
  int yDispMin1 = yBoxMin1*dispR;
  int yDispMax1 = yBoxMax1*dispR;

  int xDispMin2 = xBoxMin2*dispR;
  int xDispMax2 = xBoxMax2*dispR;
  int yDispMin2 = yBoxMin2*dispR;
  int yDispMax2 = yBoxMax2*dispR;

  double xMin1 = xBoxMin1*sgm;
  double xMax1 = xBoxMax1*sgm; 
  double yMin1 = yBoxMin1*sgm;
  double yMax1 = yBoxMax1*sgm;

  double xMin2 = xBoxMin2*sgm;
  double xMax2 = xBoxMax2*sgm; 
  double yMin2 = yBoxMin2*sgm;
  double yMax2 = yBoxMax2*sgm;


  int box[] = new int[NN];     /* i-th box number */  
  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 platef = 0.0;
  double platex = xMin2+2.0*sgm;
  double platev = 0.0; 
  double platem = 1.0*mas; 

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

  public void init() {
    resize(600,280);
    setBackground(Color.white);

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

    bt_open = new Button("open");
    bt_open.addActionListener(this); 

    sad= new Scrollbar(Scrollbar.HORIZONTAL,10,10,0,60);
    sad.addAdjustmentListener(this);
    scr= new Scrollbar(Scrollbar.HORIZONTAL,30,10,10,70);
    scr.addAdjustmentListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,4,5,0));
    pnl.add(new Label("  "));
    pnl.add(bt_open);
    pnl.add(sad);
    pnl.add(scr);
    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_open) {
       if (gate==1) {
         t = 0.0;
         gate = 0;
       }
     } 
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if(ev.getSource() == sad) { 
      adsrbQ = 0.01*(double)(sad.getValue());
    } else if(ev.getSource() == scr) { 
      contTemp = 10.0*(double)(scr.getValue());
    }
  }

  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;
    double tmp;

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

    gOff.setColor(Color.black); 
    tmp = temperature();
    if (tmp<contTemp) { gOff.setColor(Color.red); };
    gOff.drawRect(xDispMin1+30,yDispMin1+30,xDispMax1+dispR,yDispMax1+dispR);
    gOff.setColor(Color.black);
    gOff.drawRect(xDispMin2+dispR+30,yDispMin2+30,
               xDispMax2-xDispMax1,yDispMax2-yDispMin2+dispR);
    gOff.setColor(Color.white); 
    gOff.drawRect(xDispMin2+dispR+30,yDispMin2+30,1,yDispMax2-yDispMin2+dispR);
    gOff.drawRect(xDispMax2+dispR+30,yDispMin2+30,1,yDispMax2-yDispMin2+dispR);
    
    if (gate==1) {
      gOff.setColor(Color.green);
      ix = (int)((double)dispR*platex/sgm);
      gOff.drawRect(ix+dispR+30,yDispMin2+30,2,yDispMax2-yDispMin2+dispR);
    }

    gOff.setColor(Color.blue);
    for (i=0; i<Nmt; i++) {
      ix = (int)(paintSizeX*xx[i]/xMax)+30;
      iy = (int)(paintSizeY*yy[i]/yMax)+30;
      gOff.fillOval(ix,iy, 10, 10);
    }

    gOff.setColor(Color.black);
    gOff.drawString("Q="+adsrbQ+"eV",600/4*2+10,40);
    gOff.drawString("cont T="+contTemp+"K",600/4*3+10,40);
    gOff.setColor(Color.blue);
    gOff.drawString("t="+Math.floor(t*1.0e12)+"(ps)",600/4*2+10,60);
    gOff.drawString("T="+temperature()+" K",600/4*3+10,60);
    gOff.drawString("N="+kosu()+" ",30+80,wakuY+50);
  }


/*-------------------------------------  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,n;
    double s;

    n = 0; s = 0.0; 
    for (i=0; i<Nmt; i++) {
      if (box[i]==0) { 
        n += 1;
        s += vx[i]*vx[i] + vy[i]*vy[i];
      }
    }
    if (n==0) { 
      return (0.0); 
    } else {
      return ((double)Math.floor(mas*s/((double)n*2.0*1.38e-23)*10.0)/10.0 );
    }
  }

  int kosu() {
    int i,n;

    n = 0; 
    for (i=0; i<Nmt; i++) {
      if (box[i]==0) { 
        n += 1;
      }
    }
    return (n); 
  }

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

  void setInitialPosition(){
    int i,ix,iy; 

    for (i=0; i<Nmt; i++) { 
      xx[i] = 2.0*sgm*(double)(i/10)+1.0*sgm;
      yy[i] = 2.0*sgm*(double)(i%10)+1.0*sgm;
      vx[i] = 500.0*Math.random()-250.0;
      vy[i] = 500.0*Math.random()-250.0;
      ffx[i] = 0.0;
      ffy[i] = 0.0;
    }
  }

  void calcposition(){
    int i;

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

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

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

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

    platef = 0.0;
    for (i=0; i<Nmt; i++) {
      if (box[i]==0) { 
        if (xx[i] < xMin1) { 
          xx[i] = xMin1; vx[i] = -rr*vx[i]; vy[i] = rr*vy[i]; 
        }
        if (xx[i] > xMax1) {
          if (yy[i]>yMin2 && yy[i]<yMax2 ) {
            box[i]=1;
          } else {
            xx[i] = xMax1; vx[i] = -rr*vx[i]; vy[i] = rr*vy[i];
          } 
        }
        if (yy[i] < yMin1) { 
          yy[i] = yMin1; vx[i] = rr*vx[i]; vy[i] = -rr*vy[i]; 
        }
        if (yy[i] > yMax1) {
          yy[i] = yMax1; vx[i] = rr*vx[i]; vy[i] = -rr*vy[i];
        }
      } else if (box[i]==1) { 
        if (xx[i] < xMin2) { 
          if (yy[i]>yMin2 && yy[i]<yMax2 ) {
            box[i]=0;
          } else {
            xx[i] = xMin2; vx[i] = -vx[i]; vy[i] = vy[i];
          } 
        }
        if (xx[i] > platex && gate==1) {
          platef += 2.0*mas*vx[i]/dt;
          xx[i] = platex; vx[i] = -vx[i]; vy[i] = vy[i];
        }
        if (xx[i] > xMax2) {
          box[i] = 2;
          /* xx[i] = xMax2; vx[i] = -vx[i]; vy[i] = vy[i]; */
        }
        if (yy[i] < yMin2) { 
          yy[i] = yMin2; vx[i] = vx[i]; vy[i] = -vy[i]; 
        }
        if (yy[i] > yMax2) {
          yy[i] = yMax2; vx[i] = vx[i]; vy[i] = -vy[i];
        }
      }
    } 

    platev += 0.5*dt*(platef-load)/platem;
    platex += platev*dt;
    if (platex>xMax2) { 
      platex=xMax2; platev=0.0; 
    } else if (platex<xMin2+sgm) { 
      platex=xMin2+sgm; platev=0.0;
    }    
  }

  void forceCalc() {
    int i,j;
    double 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(j=i+1;j<Nmt;j++) {
        rij = Math.sqrt((xx[i]-xx[j])*(xx[i]-xx[j])+(yy[i]-yy[j])*(yy[i]-yy[j]));
        f = force(rij);
        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;
      }
    }
    for(i=0;i<Nmt;i++) {
      if (box[i]==0) {
        if (xx[i]<1.0*sgm) {
          ffx[i] -= adsrbQ*1.6e-19/sgm;
        }
        if (yy[i]<1.0*sgm) {
          ffy[i] -= adsrbQ*1.6e-19/sgm;
        }
        if (yy[i]>yMax1-sgm) {
          ffy[i] += adsrbQ*1.6e-19/sgm;
        }
      }
    }
  }

  double force(double r) {
    double ri,r3,r6;

    ri = sgm/r;
    r3 = ri*ri*ri;
    r6 = r3*r3;
    return (24.0*eps*r6*(2.0*r6-1.0)/r);
  }

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

}  

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