/** applet No. 0080
 *
 * time reverse - molecular dynamics 2D
 * - Lennard-Jones potential - Ar
 *
 * Created by Ikeuchi Mitsuru on September 09 2002.
 * Copyright (c) 2002 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2002.09.09  created
 * ver 0.0.2   2006.11.23  bug fixed
 *
 */

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

public class timeReverse extends Applet
  implements ItemListener, ActionListener, AdjustmentListener, Runnable {

/* -------------------------------------- set global ------ */

  Choice ch_tempCont, ch_view;
  Button bt_reset, bt_start, bt_stop, bt_reverse;
  Scrollbar sc_temp;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int dispWidth = 630;
  int dispHeight = 320;

  int sleepTime = 30; 
  double sgm = 3.40e-10;        // sigma of L-J potential for Ar (m)
  double eps = 0.0104*1.602e-19;// epsilon of L-J potential for Ar (J) 
  double mas = 39.948*1.661e-27;// mass of L-J potential for Ar (kg)

  double xMax = 6.0e-9;
  double yMax = 6.0e-9;

  int Nmt = 18;
  double t = 0.0;
  double dt = 5.0*1.0e-15;

  int resetSW = 0;
  int started = 1;
  int dispMode = 0;

  int tempMode = 0;             // 0-adiabatic 1-scaling 2-wall cont
  double contTemp = 200.0;      // control temperature 
  double boxDispWidth = 250.0;
  double dispScale = (boxDispWidth/xMax);
  double viewScale = 1.0;

  int NN = Nmt;
  double xx[] = new double[NN];
  double yy[] = new double[NN];
  double vx[] = new double[NN];
  double vy[] = new double[NN];
  double ffx[] = new double[NN];
  double ffy[] = new double[NN];

  int MM = 200;
  int mmp = 0;
  int mtt[] = new int[MM];
  double mxx[][] = new double[NN][MM];
  double myy[][] = new double[NN][MM];

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

  public void init() {

    resize(dispWidth,dispHeight);
    setBackground(Color.white);

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

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

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

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

    bt_start= new Button("start");
    bt_start.addActionListener(this); 

    bt_stop= new Button("stop");
    bt_stop.addActionListener(this); 

    bt_reverse= new Button("reverse");
    bt_reverse.addActionListener(this); 

    sc_temp = new Scrollbar(Scrollbar.HORIZONTAL,300,10,10,610);
    sc_temp.addAdjustmentListener(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_start);
    pnl.add(bt_stop);
    pnl.add(bt_reverse);
    pnl.add(new Label("  "));
    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 = null;
  }

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

  public void actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_reset){ 
      resetSW = 1;
	  if (dt<0) dt = -dt;
    } else if (ev.getSource() == bt_start){
      if (started==0) {
        started = 1;
        if (dt<0) dt = -dt;
      }
    } else if (ev.getSource() == bt_stop){
      started = 0;
    } else if (ev.getSource() == bt_reverse){
      if (started==0) {
        started = -1;
        if (dt>0) dt = -dt;
      }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if(ev.getSource() == sc_temp) { 
      contTemp = 1.0*(double)(sc_temp.getValue());
    }
  }

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

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

    gbx = 30;
    gby = 60;
    sc = viewScale*dispScale;

    if (dispMode==0) {
      ballPlot(gbx,gby,sc);

    } else if (dispMode==1) {
      velocityPlot(gbx,gby,sc, 0.05);

    } else if (dispMode==2) {
      tracePlot(gbx,gby,sc);

    } else if (dispMode==3) {
      ;

    }

    gOff.setColor(Color.black);
    gOff.drawString("t="+(int)(t*1.0e14+0.5)/100.0+" ps",10,40);
    if (started==0) {
      gOff.drawString("stop",dispWidth/6*1+10,40);
    } else if (started==1) {
      gOff.drawString("start",dispWidth/6*1+10,40);
    } else if (started==-1) {
      gOff.drawString("reverse",dispWidth/6*1+10,40);
    }
    gOff.drawString("T="+(int)(temperature()*10+0.5)/10.0+" K",dispWidth/6*4+10,40);
    gOff.drawString("view",dispWidth/6*5+10,40);

    gOff.drawString("Box="+(int)(xMax*1.0e10+0.5)/10.0+" x "+(int)(yMax*1.0e10+0.5)/10.0+" nm",500,80);
  }

/*------------------------------------  plot method  ----------*/

  void ballPlot(int gbx, int gby, double sc) {
    int i, ix, iy, ballSize;

    ballSize = (int)(sgm*sc);
    drawBox(gbx, gby, sc, ballSize);

    gOff.setColor(Color.blue);
    for (i=0; i<Nmt; i++) {
      ix = (int)(sc*xx[i])+gbx;
      iy = (int)(sc*yy[i])+gby;
      gOff.fillOval(ix-ballSize/2,iy-ballSize/2, ballSize, ballSize);
    }
  }

  void velocityPlot(int gbx, int gby, double sc, double mag) {
    int i, ix, iy, ix2, iy2, ballSize;

    ballSize = (int)(sgm*sc);
    drawBox(gbx, gby, sc, ballSize);

    gOff.setColor(Color.blue);
    for (i=0; i<Nmt; i++) {
      ix = (int)(sc*xx[i])+gbx;
      iy = (int)(sc*yy[i])+gby;
      ix2 = ix+(int)(mag*vx[i]);
      iy2 = iy+(int)(mag*vy[i]);
      gOff.setColor(Color.getHSBColor(0.66f,0.2f,0.9f ));
      gOff.fillOval(ix-ballSize/2,iy-ballSize/2, ballSize, ballSize);
      gOff.setColor(Color.blue);
      gOff.drawLine(ix,iy, ix2, iy2);
    }
  }

  void tracePlot(int gbx, int gby, double sc) {
    int i, j, k, ix, iy, ix2, iy2, ballSize;

    ballSize = (int)(sgm*sc);
    drawBox(gbx, gby, sc, ballSize);

    for (i=0; i<Nmt; i++) {
      ix = (int)(sc*xx[i])+gbx;
      iy = (int)(sc*yy[i])+gby;
      gOff.setColor(Color.getHSBColor(0.66f,0.2f,0.9f ));
      gOff.fillOval(ix-ballSize/2,iy-ballSize/2, ballSize, ballSize);
    }

    for (i=0; i<Nmt; i++) {
	  for (j=MM-2; j>1; j--) {
	    k = (mmp-j+MM)%MM;
	    ix = (int)(sc*mxx[i][k])+gbx;
		iy = (int)(sc*myy[i][k])+gby;
		k = (mmp-j-1+MM)%MM;
        ix2 = (int)(sc*mxx[i][k])+gbx;
		iy2 = (int)(sc*myy[i][k])+gby;
		if (mtt[k]==1) {
		  gOff.setColor(Color.getHSBColor(0.66f,0.8f,0.6f ));
		} else if (mtt[k]==-1) {
		  gOff.setColor(Color.getHSBColor(0.05f,0.8f,0.6f ));
		} else {
		  gOff.setColor(Color.getHSBColor(0.66f,0.01f,0.9f ));
		}
		gOff.drawLine(ix,iy,ix2,iy2);
	  }
    }
  }

  void drawBox(int gbx, int gby, double sc, int ballSize) {
    int ix, iy, gball;
    double tmp;

    gball = ballSize/2;
    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,gby-gball,ix+2*gball,iy+2*gball);
  }

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

  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 ( s/((double)Nmt*1.38e-23) );
  }

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

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

  void setInitialCondition() {
    int i; 
    double s;

    t = 0.0;

    s = xMax/20.0;
    for (i=0; i<Nmt; i++) {
      xx[i] = s*3.0*(double)(i/6)+2.5*s;
      yy[i] = s*3.0*(double)(i%6)+2.5*s;
      vx[i] = 800.0*Math.random()-400.0;
      vy[i] = 800.0*Math.random()-400.0;
      ffx[i] = 0.0;
      ffy[i] = 0.0;
    }
    vAjustment();
	initTrace();
  }

  void initTrace() {
    int i,j;

    for (i=0; i<Nmt; i++) {
      for (j=0; j<MM; j++) {
	    mxx[i][j] = xx[i];
		myy[i][j] = yy[i];
		mtt[i] = 0;
	  }
    }
    mmp = 0;
  }

/* ----------------------------- timeEvolution -------------- */

  void timeEvolution() {
    int i;

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

    if (started!=0) {
      for (i=0; i<20; i++) {
        timeStep();
      }
	  if (tempMode==1) vAjustment();
	  memoryTrace();
    }
  }

  void memoryTrace() {
    int i;

    for (i=0; i<Nmt; i++) {
      mxx[i][mmp] = xx[i];
      myy[i][mmp] = yy[i];
	  mtt[mmp] = started;
	}
    mmp = (mmp+1)%MM;
  }

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

    dtv2 = dt/2.0;
    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++) {
      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;
    }

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

  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 applet -------------- */

}