/** applet No. t06
 *
 * molecular dynamics 3D
 * Ar atom - Lennard-Jones 12-6 potential
 * - MD3D class, Display3D class
 *
 * Created by Ikeuchi Mitsuru on July 28 2007.
 * Copyright (c) 2007-2008 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2007.07.28  created
 * ver 0.0.2   2008.06.22  added setGUIPannel()
 *
 */

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

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

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

	Thread th = null; // for run()-paint() thread
	MD3D md = new MD3D();
	Display3D dsp = new Display3D(md);

	// for event
	private Choice ch_tempCont, ch_view;
	private Button bt_reset, bt_startStop, bt_step;
	private Scrollbar sc_temp;
	private int dgX, dgY, dgXb, dgYb; // mouse control

	// for off-paint buffer
	private Image imgOff;
	private Graphics gOff;

	private int sleepTime = 30;
	private int thCount = 0;


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

	public void init() {
		Dimension dim;

		resize(dsp.getDispWidth(),dsp.getDispHeight());
		setBackground(Color.white);

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

		setGUIPannel();
	}

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

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




	// ===================================== GUI control ===========

	private void setGUIPannel() {

		setChoice();
		setButton();
		setScrollbar();
		setMouse();

		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(new Label("  "));
		pnl.add(ch_tempCont);
		pnl.add(sc_temp);
		pnl.add(ch_view);
		add(pnl,"North");
	}


	// choice
	private void setChoice() {
		ch_tempCont = new Choice();
		ch_tempCont.add("adiabatic");
		ch_tempCont.add("T scaling"); 
		ch_tempCont.add("wall cont.");
		ch_tempCont.addItemListener(this);
		ch_tempCont.select("wall cont."); 

		ch_view = new Choice();
		ch_view.add("ball");
		ch_view.add("line");
		ch_view.addItemListener(this);
		ch_view.select("ball");
	}

	public void itemStateChanged(ItemEvent ev){
		if (ev.getSource() == ch_tempCont) {
			md.setTempMode( ch_tempCont.getSelectedIndex() );
		} else if (ev.getSource() == ch_view) {
			dsp.setDispMode( ch_view.getSelectedIndex() );
		}
	}


	// button
	private void setButton() {
		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); 
	}

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


	// scrollbar
	private void setScrollbar() {
		sc_temp = new Scrollbar(Scrollbar.HORIZONTAL,300,10,10,610);
		sc_temp.addAdjustmentListener(this);
	}

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


	// mouse
	private void setMouse() {
		addMouseListener(this);  
		addMouseMotionListener(this);
	}

	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) {
			md.timeEvolution();
			dsp.offPaint();
			repaint();
			try {
				Thread.sleep(sleepTime);
			} catch (InterruptedException e) { }
		}
	}

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

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

}




// ==================================== display 3D class ===========

class Display3D {

	private MD3D md;
	private Graphics gOff;
	private int dispWidth = 630;
	private int dispHeight = 460;

	private int dispMode = 0;

	private double theta = -20.0*3.14/180.0; // rotation angle theta
	private double fai = 10.0*3.14/180.0; // rotation angle fai

	private double dispScale = (250.0/(6.0e-9)); // 250dots/6nm
	private double viewScale = 1.0;

	private int NN = MD3D.getNmt();
	private int srtz[][] = new int[2][NN];
	private double px[] = new double[NN]; 
	private double py[] = new double[NN]; 
	private double pz[] = new double[NN]; 
	private double wkx[] = new double[8]; 
	private double wky[] = new double[8]; 
	private double wkz[] = new double[8]; 
	private double pwkx[] = new double[8]; 
	private double pwky[] = new double[8]; 
	private double pwkz[] = new double[8]; 
	private int boxp[][] = {
		{0,1},{0,2},{0,4},{1,3},{1,5},{2,3},{2,6},{4,5},{4,6},{3,7},{5,7},{6,7}
	};


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

	Display3D(MD3D amd) {
		md = amd;
	}


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

	public void offPaint() {
		int gbx,gby;
		double sc;

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

		gbx = 80;
		gby = 120;

		if (dispMode==0) {
			ball3dPlot(gbx,gby,1.0,0);

		} else if (dispMode==1) {
			ball3dPlot(gbx,gby,0.5,1);

		} else if (dispMode==2) {
			;

		}

		gOff.setColor(Color.blue);
		gOff.drawString("t="+(int)(md.getTime()*1.0e13+0.5)/10.0+" ps",10,40);
		gOff.setColor(Color.black);
		if (md.getStartSW()==0) {
			gOff.drawString("stopped",630/6*1+10,40);
		} else {
			gOff.drawString("started",630/6*1+10,40);
		}
		gOff.drawString("T cont.",630/6*3+10,40);
		gOff.drawString("Tc="+md.getContTemp()+" K",630/6*4+10,40);

		gOff.setColor(Color.blue);
		gOff.drawString("T="+(int)(md.getTemperature()*10+0.5)/10.0+" K",630/6*4+10,60);
		gOff.setColor(Color.black);
		gOff.drawString("view",630/6*5+10,40);

		gOff.drawString("Box="+(int)(md.getXMax()*1.0e10+0.5)/10.0+
			" x "+(int)(md.getYMax()*1.0e10+0.5)/10.0+
			" x "+(int)(md.getZMax()*1.0e10+0.5)/10.0+" nm",10,60);

	}


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

	private void ball3dPlot(int gbx, int gby, double size, int mode) {
		int i,j,jj,gx,gy,gz,g2x,g2y, Nmt;
		double sc,sz;

		setWaku();
		rotate(theta,fai);
		zSort();
		Nmt = md.getNmt();

		drawWaku(gbx,gby,0);

		sc = dispScale*viewScale;

		for (i=0; i<Nmt; i++) {
			j = srtz[0][i];
			sz = sc*size*md.getDia();
			gx = (int)(sc*px[j])+gbx;
			gy = (int)(sc*py[j])+gby;
			gz = (int)(sz*(0.3*pz[j]/md.getZMax()+0.7));

			if (mode==1) {
				for(jj=j+1;jj<Nmt;jj++) {
					if (distance(j,jj)<1.3*md.getDia()) {
						g2x = (int)(sc*px[jj])+gbx;
						g2y = (int)(sc*py[jj])+gby;
						gOff.drawLine(gx,gy, g2x, g2y);
					}
				}
			}
			gOff.setColor(Color.getHSBColor(0.15f,0.8f,(float)(0.4*pz[j]/md.getZMax()+0.4)));
			gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
		}
		drawWaku(gbx,gby,1);
	}

	private double distance(int i,int j) {
		double x,y,z;

		x = md.getxx(i)-md.getxx(j);
		y = md.getyy(i)-md.getyy(j);
		z = md.getzz(i)-md.getzz(j);

		return(Math.sqrt(x*x+y*y+z*z));
	}

	//  draw waku

	private void drawWaku(int gbx, int gby, int imode) {
		int i,gx,gy,g2x,g2y,farPoint;
		double sc,tmp;

		sc = dispScale*viewScale;
		farPoint = findMin();

		for (i=0; i<12; i++) {
			gx = (int)(sc*pwkx[boxp[i][0]])+gbx;
			gy = (int)(sc*pwky[boxp[i][0]])+gby;
			g2x = (int)(sc*pwkx[boxp[i][1]])+gbx;
			g2y = (int)(sc*pwky[boxp[i][1]])+gby;
			if (imode==0) {
				if (boxp[i][0]==farPoint || boxp[i][1]==farPoint) { 
					gOff.setColor(Color.lightGray);
					gOff.drawLine(gx,gy, g2x, g2y);
				}
			} else if (imode==1) {
				if (boxp[i][0]!=farPoint && boxp[i][1]!=farPoint) { 
					gOff.setColor(Color.gray);
					if (md.getTempMode()==2) {
						tmp = md.getTemperature();
						if (tmp>md.getContTemp()) {
							gOff.setColor(Color.black); 
						} else {
							gOff.setColor(Color.red);
						}
					}
					gOff.drawLine(gx,gy, g2x, g2y);
				}
			}
		}
	}

	private int findMin() {
		int i,im;
		double m;

		im = 0; m = pwkz[im];
		for (i=0; i<8; i++) {
			if (pwkz[i]<m) {
				m = pwkz[i];
				im = i;
			}
		}
		return im;
	}

	// set waku

	private void setWaku() {
		int i,ix,iy,iz;
		double sgm;

		i=0;
		sgm = 4.0e-10;

		for (ix=0; ix<2; ix++) {
			for (iy=0; iy<2; iy++) {
				for (iz=0; iz<2; iz++) {
					wkx[i] = -0.5*sgm+(md.getXMax()+sgm)*(double)(ix);
					wky[i] = -0.5*sgm+(md.getYMax()+sgm)*(double)(iy);
					wkz[i] = -0.5*sgm+(md.getZMax()+sgm)*(double)(iz);
					i += 1;
				}
			}
		}
	}

	// rotate

	private void rotate(double th,double fi) {
		int i,Nmt;
		double cx,cy,cz;
		double x,y,z,cosTh,sinTh,cosFi,sinFi;

		cx = 0.5*md.getXMax(); // x-center of rotation
		cy = 0.5*md.getYMax(); // y-center of rotation
		cz = 0.5*md.getZMax(); // z-center of rotation

		Nmt = md.getNmt();
		cosTh = Math.cos(th); sinTh = Math.sin(th);
		cosFi = Math.cos(fi); sinFi = Math.sin(fi);

		for (i=0; i<Nmt; i++) {
			y = md.getyy(i); z = md.getzz(i);
			py[i] = cosFi*(y-cy)+sinFi*(z-cz) + cy; 
			pz[i] = -sinFi*(y-cy)+cosFi*(z-cz) + cz;
		} 
		for (i=0; i<Nmt; i++) {
			x = md.getxx(i);
			px[i] = cosTh*(x-cx)+sinTh*(pz[i]-cz) + cx; 
			pz[i] = -sinTh*(x-cx)+cosTh*(pz[i]-cz) + cz;
		} 

		for (i=0; i<8; i++) {
			pwky[i] = cosFi*(wky[i]-cy)+sinFi*(wkz[i]-cz) + cy; 
			pwkz[i] = -sinFi*(wky[i]-cy)+cosFi*(wkz[i]-cz) + cz;
		} 
		for (i=0; i<8; i++) {
			pwkx[i] = cosTh*(wkx[i]-cx)+sinTh*(pwkz[i]-cz) + cx; 
			pwkz[i] = -sinTh*(wkx[i]-cx)+cosTh*(pwkz[i]-cz) + cz;
		}  
	}

	// sort z

	private void zSort() {
		int i,j,w,Nmt;

		Nmt = md.getNmt();

		for (i=0; i<Nmt; i++) {
			srtz[0][i] = i;
			srtz[1][i] = (int)(32767.0*pz[i]/md.getZMax());
		}
		qSort(0,Nmt-1);    
	}

	private void qSort(int i,int j) {
		int k,p;

		if(i==j) return;
		p = i+1;
		while (p<=j && srtz[1][i]==srtz[1][p]) p++;
		if (p<=j) {
			if (srtz[1][i]>=srtz[1][p]) p = i;
			k = qSortPartition(i,j,srtz[1][p]);
			qSort(i,k-1);
			qSort(k,j);
		}
	}

	private int qSortPartition(int i,int j, int x) {
		int l,r,w;

		l = i; r = j;
		while (l<=r) {
			while (l<=j && srtz[1][l]<x) l++;
			while (r>=i && srtz[1][r]>=x) r--;
			if (l>r) break;
			w = srtz[0][l]; srtz[0][l] = srtz[0][r]; srtz[0][r] = w;
			w = srtz[1][l]; srtz[1][l] = srtz[1][r]; srtz[1][r] = w;
			l++; r--;
		}
		return l;
	}


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

	public void setGraphics(Graphics gg) {
		gOff = gg;
	}

	public void setViewDelta(double xDisplacement, double yDisplacement) {

		theta += 0.5*3.14/180.0*xDisplacement;
		fai += 0.5*3.14/180.0*yDisplacement;
	}

	public void setDispMode(int m) { dispMode = m; }

	public int getDispWidth() { return dispWidth; }
	public int getDispHeight() { return dispHeight; }
}




// ========================= molecular dynamics 3D class ===========

class MD3D {

	private double t = 0.0; // time (s)
	private double dt = 10.0*1.0e-15; // time division (s)
	private double xMax = 6.0e-9; // x-Box size (m)
	private double yMax = 6.0e-9; // y-Box size (m)
	private double zMax = 6.0e-9; // z-Box size (m)

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

	private int tempMode = 2; // 0-adiabatic 1-scaling 2-wall cont
	private double contTemp = 300.0; // control temperature (K)

	private static int Nmt = 125; // number of molecules
	private int NN = Nmt; // dimension of array
	private double xx[] = new double[NN]; // i-th molecules x-position
	private double yy[] = new double[NN]; // i-th molecules y-position
	private double zz[] = new double[NN]; // i-th molecules z-position
	private double vx[] = new double[NN]; // i-th molecules x-velocity
	private double vy[] = new double[NN]; // i-th molecules y-velocity
	private double vz[] = new double[NN]; // i-th molecules z-velocity
	private double ffx[] = new double[NN]; // i-th molecules x-force
	private double ffy[] = new double[NN]; // i-th molecules y-force
	private double ffz[] = new double[NN]; // i-th molecules z-force

	private int resetSW = 0; // 1-reset : reset and set resetSW = 0
	private int startSW = 1; // 0-stop 1-start
	private int stepSW = 0;  // 0-stop 1-step
	private int count = 0;


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

	MD3D() {
		setInitialCondition();
	}

	//

	private void setInitialCondition() {
		int i,ix,iy,iz; 
		double s;

		t = 0.0;
		count = 0;

		s = xMax/5.0;
		i = 0;
		for (ix=0; ix<5; ix++) {
			for (iy=0; iy<5; iy++) {
				for (iz=0; iz<5; iz++) {
					xx[i] = s*(ix+0.5);
					yy[i] = s*(iy+0.5);
					zz[i] = s*(iz+0.5);
					vx[i] = 500.0*(Math.random()-0.5);
					vy[i] = 500.0*(Math.random()-0.5);
					vz[i] = 500.0*(Math.random()-0.5);
					ffx[i] = 0.0;
					ffy[i] = 0.0;
					ffz[i] = 0.0;
					i += 1;
				}
			}
		}
		vAjustment(contTemp);
	}


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

	public void timeEvolution() {
		int i;

		if (resetSW==1) {
			setInitialCondition();
			resetSW = 0;
		}
	
		if (startSW==1) {
			timeStep20();

		} else {
			if (stepSW==1) {
				timeStep20();
				stepSW = 0;
			}
		}
	}

	private void timeStep20() {
		int i;

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

	private void timeStep() {
		int i; 
		double dtv2, a;

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

		for (i=0; i<Nmt; i++) {
			a = dtv2/mas;
			vx[i] += a*ffx[i];
			vy[i] += a*ffy[i];
			vz[i] += a*ffz[i];
			xx[i] += vx[i]*dt;
			yy[i] += vy[i]*dt;
			zz[i] += vz[i]*dt;
		}
		forceCalc();
		for (i=0; i<Nmt; i++) {
			a = dtv2/mas;
			vx[i] += a*ffx[i];
			vy[i] += a*ffy[i];
			vz[i] += a*ffz[i];
		}

		setBoundary();
	}

	private void forceCalc() {
		int i,j;
		double xij,yij,zij,rij,f,fxij,fyij,fzij;

		for(i=0;i<Nmt;i++) {
			ffx[i] = 0.0; ffy[i] = 0.0; ffz[i] = 0.0;
		}

		for(i=0;i<Nmt;i++) {
			for(j=i+1;j<Nmt;j++) {
				xij = xx[i]-xx[j]; yij = yy[i]-yy[j]; zij = zz[i]-zz[j];
				rij = Math.sqrt(xij*xij + yij*yij + zij*zij);
				f = force(rij);
				fxij = f*xij/rij;
				fyij = f*yij/rij;
				fzij = f*zij/rij;

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

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

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

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

		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++) {
			if (xx[i] < 0.0) { 
				xx[i] = 0.0; vx[i] = -rr*vx[i]; vy[i] = rr*vy[i]; vz[i] = rr*vz[i];
			}
			if (xx[i] > xMax) {
				xx[i] = xMax; vx[i] = -rr*vx[i]; vy[i] = rr*vy[i]; vz[i] = rr*vz[i]; 
			}
			if (yy[i] < 0.0) { 
				yy[i] = 0.0; vx[i] = rr*vx[i]; vy[i] = -rr*vy[i]; vz[i] = rr*vz[i]; 
			}
			if (yy[i] > yMax) {
				yy[i] = yMax; vx[i] = rr*vx[i]; vy[i] = -rr*vy[i]; vz[i] = rr*vz[i];
			}
			if (zz[i] < 0.0) { 
				zz[i] = 0.0; vx[i] = rr*vx[i]; vy[i] = rr*vy[i]; vz[i] = -rr*vz[i];
			}
			if (zz[i] > zMax) {
				zz[i] = zMax; vx[i] = rr*vx[i]; vy[i] = rr*vy[i]; vz[i] = -rr*vz[i];
			}
		}
	}


	// --------------------------------------- utilities -----------

	private double temperature() {
		int i;
		double s,ke;

		s = 0.0; 
		for (i=0; i<Nmt; i++) {
			s += (vx[i]*vx[i] + vy[i]*vy[i]+ vz[i]*vz[i]);
		}
		ke = 0.5*mas*s/((double)Nmt);
		return ( ke/(1.5*1.38e-23) );
	}

	private void vAjustment(double temp) {
		int i;
		double tp, r;

		tp = temperature();
		r = Math.sqrt(temp/tp);
		for (i=0; i<Nmt; i++) {
			vx[i] = r*vx[i];
			vy[i] = r*vy[i];
			vz[i] = r*vz[i];
		}
	}


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

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

	public void setTempMode(int m) { tempMode = m; }
	public int getTempMode() { return tempMode; }
	public void setContTemp(double tmp) { contTemp = tmp; }
	public double getContTemp() { return contTemp; }

	// get data
	public static int getNmt() { return Nmt; }
	public double getXMax() { return xMax; }
	public double getYMax() { return yMax; }
	public double getZMax() { return zMax; }
	public double getDia() { return sgm; }

	public double getTime() { return t; }
	public double getxx(int i) { return xx[i]; }
	public double getyy(int i) { return yy[i]; }
	public double getzz(int i) { return zz[i]; }
	public double getvx(int i) { return vx[i]; }
	public double getvy(int i) { return vy[i]; }
	public double getvz(int i) { return vz[i]; }
	public double getTemperature() { return temperature(); }

}