/** applet No. 0318
 *
 * moon travel - faster Morse MD
 *
 * Created by Ikeuchi Mitsuru on July 06 2003.
 * Copyright (c) 2003 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2003.07.06  created
 * ver 0.0.2   2006.12.11  bug fixed
 *
 */

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

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

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

  Thread th = null;

  Button bt_reset, bt_start, bt_stop;
  Choice ch_view;
  Scrollbar sc_scale, sc_temp;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 50;
  int atomSizeX = 25; 
  int atomSizeY = 25; 
  int dispR = 10;           /* display molecular radious */
  int paintSizeX = atomSizeX*dispR; 
  int paintSizeY = atomSizeY*dispR;
  int wakuX = paintSizeX + dispR;
  int wakuY = paintSizeY + dispR;

  int Nmt = 1161;             /* number of molecules    */
  int NN = Nmt + 1;
  double t = 0.0;              /* time (s)               */
  double dt = 5.0*1.0e-15;     /* time division (s)      */
  double sgm = 3.40e-10;       /* L-J sigma for Ar (m)   */

  int resetSW = 0;
  int startSW = 1;
  int goSW = 0;
  int count = 0;
  double contTemp = 4000.0;
  double velocity = 1000.0;
  double magn = 0.5;
  int dispMode = 0;

  double rc = 1.0e-9;
  double rc2 = rc*rc;

  double rreg = 1.5e-9;
  double rreg2 = rreg*rreg;

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

  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},
  /*21 bld*/{ 39.948*AU,  0.0*EE, 0.0104*EE, 1.3400/AA, 3.816*AA},
  /*22 fix*/{ 1.0      ,  0.0*EE, 0.0104*EE, 1.3400/AA, 3.816*AA}
  };
  int NmatKind = 23;
  double MorseD[][] = new double[NmatKind][NmatKind];
  double MorseA[][] = new double[NmatKind][NmatKind];
  double Morser0[][] = new double[NmatKind][NmatKind];
  int MAT0 = 17;
  int MAT1 = 0;
  int MAT2 = 0;
  int MAT3 = 0;
  int MAT4 = 22;
  int MAT5 = 21; /* do not change */
  
  int kind[] = new int[NN];       /* i-th molec kind */
  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 */

  int kindT[] = new int[NmatKind];
  double forceTable[][][] = new double[6][6][12000];

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

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

  public void init() {
    resize(630,320);
    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);

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

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

    sc_scale = new Scrollbar(Scrollbar.HORIZONTAL,50,10,20,210);
    sc_scale.addAdjustmentListener(this);

    sc_temp = new Scrollbar(Scrollbar.HORIZONTAL,40,10,10,110);
    sc_temp.addAdjustmentListener(this);

    ch_view = new Choice();
    ch_view.addItem("ball");
    ch_view.addItem("color"); 
    ch_view.addItem("line");
    ch_view.addItem("both");
    ch_view.addItemListener(this);
    ch_view.select("ball");

    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_start);
    pnl.add(bt_stop);
    pnl.add(sc_temp);
    pnl.add(sc_scale);
    pnl.add(ch_view);
    add(pnl,"North");

    setMorse();
    setInitialPosition();

    velocity = Math.sqrt(3.0*1.38e-23*contTemp/Morse[MAT0][0]);
  }

  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_view){
      dispMode = ch_view.getSelectedIndex();
    }
  }

  public void actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_reset){
      resetSW = 1;
	  goSW = 0;
    } else if(ev.getSource() == bt_start){
	  startSW = 1; 
      if (goSW==0) {
        excite();
        goSW = 1;
      }
    } else if(ev.getSource() == bt_stop) {
      startSW = 0;
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_scale) { 
      magn = 0.01*(double)(sc_scale.getValue());
    } else if (ev.getSource() == sc_temp) { 
      contTemp = 100.0*(double)(sc_temp.getValue());
      velocity = Math.sqrt(3.0*1.38e-23*contTemp/Morse[MAT0][0]);
    }
  }

  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/magn)*(dgX-dgXb);
    yShift += (1.0/magn)*(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,j,k,ix,iy,ir,ix2,iy2;
    int gbx,gby;
    double tm,r2,rd2;

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

    gbx = 30+(int)(magn*xShift);
    gby = 120+(int)(magn*yShift);

    if (dispMode==0) {

      for (i=0; i<Nmt; i++) {
        ix = gbx + (int)(magn*2.5*xx[i]/AA);
        iy = gby + (int)(magn*2.5*yy[i]/AA);
        ir = 1 + (int)(magn*1.5*(Morse[kind[i]][4])/AA);

        gOff.setColor(Color.getHSBColor(
          (float)(ballColor(kind[i])),0.99f,0.99f));
        gOff.fillOval(ix-ir/2,iy-ir/2, ir, ir);
      }

    } else if (dispMode==1) {

      for (i=0; i<Nmt; i++) {
        ix = gbx + (int)(magn*2.5*xx[i]/AA);
        iy = gby + (int)(magn*2.5*yy[i]/AA);
        ir = 1 + (int)(magn*1.5*(Morse[kind[i]][4])/AA);

        tm = 0.5*Morse[kind[i]][0]*(vx[i]*vx[i]+vy[i]*vy[i])/1.38e-23;
        gOff.setColor(Color.getHSBColor(
          (float)(ballColor(kind[i])),0.99f,(float)(0.3+0.6*tm/2000.0) ));
        gOff.fillOval(ix-ir/2,iy-ir/2, ir, ir);
      }

    } else if (dispMode==2) {

      for (i=0; i<Nmt-1; i++) {
        rd2 = 1.2*1.2*Morse[kind[i]][4]*Morse[kind[i]][4];
        ix = gbx + (int)(magn*2.5*xx[i]/AA);
        iy = gby + (int)(magn*2.5*yy[i]/AA);
        ir = 1 + (int)(magn*1.5*(Morse[kind[i]][4])/AA);

        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<rd2 && j>i) {
            ix2 = gbx + (int)(magn*2.5*xx[j]/AA);
            iy2 = gby + (int)(magn*2.5*yy[j]/AA);

            gOff.setColor(Color.getHSBColor( (float)(1.0-r2/rd2),0.9f,0.9f ));
            gOff.drawLine(ix,iy, ix2,iy2);
          }
        }
      }

    } else if (dispMode==3) {

      for (i=0; i<Nmt-1; i++) {
        rd2 = 1.2*1.2*Morse[kind[i]][4]*Morse[kind[i]][4];
        ix = gbx + (int)(magn*2.5*xx[i]/AA);
        iy = gby + (int)(magn*2.5*yy[i]/AA);
        ir = 1 + (int)(magn*1.5*(Morse[kind[i]][4])/AA);

        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<rd2 && j>i) {
            ix2 = gbx + (int)(magn*2.5*xx[j]/AA);
            iy2 = gby + (int)(magn*2.5*yy[j]/AA);

            gOff.setColor(Color.getHSBColor( (float)(1.0-r2/rd2),0.9f,0.9f ));
            gOff.drawLine(ix,iy, ix2,iy2);
          }
        }
      }

      for (i=0; i<Nmt; i++) {
        ix = gbx + (int)(magn*2.5*xx[i]/AA);
        iy = gby + (int)(magn*2.5*yy[i]/AA);
        ir = 1 + (int)(magn*1.0*(Morse[kind[i]][4])/AA);

        gOff.setColor(Color.getHSBColor(
          (float)(ballColor(kind[i])),0.99f,0.99f));
        gOff.fillOval(ix-ir/2,iy-ir/2, ir, ir);
      }

    }


    gOff.setColor(Color.black);
    gOff.drawString("T fuel="+(int)(contTemp)+"K",630/6*3+20,40);
    gOff.drawString("scale="+(int)(magn*100.0+0.5)+"%",630/6*4+20,40);
    gOff.drawString("view"+(int)(magn*100.0+0.5)+"%",630/6*5+20,40);
    gOff.setColor(Color.blue);
    gOff.drawString("t="+Math.floor(t*1.0e15+0.5)+"(fs)",10,40);

    gOff.drawString("speed="+(int)(xspeed()*10.0+0.5)/10.0+" m/s",450,60);
  }

  double ballColor(int knd) {
    double c;

    c = 0.0;
    if (knd==MAT0) {
      c = 0.01;
    } else if (knd==MAT1) {
      c = 0.15;
    } else if (knd==MAT2) {
      c = 0.33;
    } else if (knd==MAT3) {
      c = 0.50;
    } else if (knd==MAT4) {
      c = 0.66;
    } else if (knd==MAT5) {
      c = 0.80;
    }
    return (c);
  }

/*-------------------------------------  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*(Morse[kind[i]][0])*(vx[i]*vx[i] + vy[i]*vy[i]);
    }
    return ((double)Math.floor(s/((double)Nmt*1.38e-23)*10.0)/10.0 );
  }

  double xspeed() {
    int i,n;
    double s;

    s = 0.0; 
    n = 0;
    for (i=0; i<Nmt; i++) {
      if (kind[i]==MAT5) {
        s += vx[i];
        n += 1;
      }
    }
    return ( s/(double)n );
  }


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

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

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

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

  void setMorse(){
    int i,j; 

    for (i=0; i<NmatKind; i++) {
      for (j=0; j<NmatKind; j++) {
        MorseD[i][j] = Math.sqrt(Morse[i][2]*Morse[j][2]);
        MorseA[i][j] = 0.5*(Morse[i][3]+Morse[j][3]);
        Morser0[i][j] = 0.5*(Morse[i][4]+Morse[j][4]);
      }
    }
  }

  void setForceTable(){
    int i,j,ir,k0,k1; 
    double d,a,r0,y,r,sc;
    int kkk[] = {MAT0,MAT1,MAT2,MAT3,MAT4,MAT5};

    for (i=0; i<6; i++) {
      k0=kkk[i];
      for (j=0; j<6; j++) {
        k1=kkk[j];

        d = Math.sqrt(Morse[k0][2]*Morse[k1][2]);
        a = 0.5*(Morse[k0][3]+Morse[k1][3]);
        r0 = 0.5*(Morse[k0][4]+Morse[k1][4]);

        for (ir=100; ir<10000; ir++) {
          r = (double)ir*1.0e-13;
          y = Math.exp(-a*(r-r0));
          forceTable[i][j][ir]=2.0*d*a*y*(y-1.0);
        }
        for (ir=0; ir<100; ir++) {
          forceTable[i][j][ir]=forceTable[i][j][100];
        }

        for (ir=6000; ir<10000; ir++) {
          sc = 1.0-((double)(ir-6000))/4000.0;
          forceTable[i][j][ir] = sc*forceTable[i][j][ir];
        }
      }
    }

    i = 6-1;
    j = 6-1;
    d = 8.0*EE;
    a = Morse[MAT5][3];
    r0 = Morse[MAT5][4];

        for (ir=100; ir<10000; ir++) {
          r = (double)ir*1.0e-13;
          y = Math.exp(-a*(r-r0));
          forceTable[i][j][ir]=2.0*d*a*y*(y-1.0);
        }
        for (ir=0; ir<100; ir++) {
          forceTable[i][j][ir]=forceTable[i][j][100];
        }

        for (ir=6000; ir<10000; ir++) {
          sc = 1.0-((double)(ir-6000))/4000.0;
          forceTable[i][j][ir] = sc*forceTable[i][j][ir];
        }

  }

  void setInitialPosition(){
    int i,ix,iy; 
    double a,yd,x0,y0;

    kindT[MAT0] = 0; kindT[MAT1] = 1; kindT[MAT2] = 2;
    kindT[MAT3] = 3; kindT[MAT4] = 4; kindT[MAT5] = 5;
    setForceTable();

    a = Morse[MAT0][4]*1.0;
    x0 = 2.0*sgm;
    y0 = 2.0*sgm;

    i = 0;
    for (ix=0; ix<50; ix++) {
      for (iy=0; iy<20; iy++) {

        if (ix<30) {
          kind[i] = MAT0;
        } else {
          kind[i] = MAT5;
        }
        xx[i] = a*(ix+0.5*(iy%2))+x0;
        yy[i] = a*0.866*iy+y0;

        vx[i] = 200.0*Math.random()-100.0;
        vy[i] = 200.0*Math.random()-100.0;

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

        i += 1;

      }
    }

    for (i=1000; i<1000+70; i++) {
        kind[i] = MAT4;
        xx[i] = a*(i-1000)+x0;
        yy[i] = y0-2.0*sgm;

        vx[i] = 0.0;
        vy[i] = 0.0;

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

    for (i=1070; i<1070+70; i++) {
        kind[i] = MAT4;
        xx[i] = a*(i-1070)+x0;
        yy[i] = y0-2.0*sgm+22.0*sgm;

        vx[i] = 0.0;
        vy[i] = 0.0;

        ffx[i] = 0.0;
        ffy[i] = 0.0;
        fbx[i] = 0.0;
        fby[i] = 0.0;
    }
    for (i=1140; i<1140+21; i++) {
        kind[i] = MAT4;
        xx[i] = x0-1.0*sgm;
        yy[i] = y0-2.0*sgm+a*(i-1140);

        vx[i] = 0.0;
        vy[i] = 0.0;

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

  }


  void excite() {
    int i;

    for (i=0; i<Nmt; i++) {
      if (kind[i]==MAT0) {
        vx[i] = 0.707*velocity*(2.0*Math.random()-1.0);
        vy[i] = 0.707*velocity*(2.0*Math.random()-1.0);
      }
    }
  }

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

  void calcposition(){
    int i;

    if (resetSW==1) {
      setInitialPosition();
	  resetSW = 0;
	}
    if (startSW==1) {
      registration();
      for (i=0; i<10; i++) {
        timeEvolution();
      }
    }
  }

  void timeEvolution(){
    int i;

    t = t + dt;

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

  }

  void forceCalc() {
    int i,j,k;
    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];
        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,kind[i],kind[j]);
          fxij = f*(xx[i]-xx[j])/rij;
          fyij = f*(yy[i]-yy[j])/rij;
        } else {
          fxij = 0.0;
          fyij = 0.0;
        }
        ffx[i] +=  fxij;
        ffy[i] +=  fyij;
        ffx[j] += -fxij;
        ffy[j] += -fyij;
      }
    }

  }

  double force(double r, int knd0, int knd1) {
    return ( forceTable[kindT[knd0]][kindT[knd1]][(int)(r*1.0e13)] );
  }

}  

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