/* gas expansion - smoothed particle hydrodynamics 1D */
/* coded by Ikeuchi Mitsuru   */
/* ver 0.0.1   2005.05.12     */
/* ver 0.0.2   2005.07.07 bug of dwwvdr(r,h) fixed */

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

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

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

  Button bt_reset, bt_start, bt_stop;
  Scrollbar sc_rate;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 30;

  int incRate = 10;
  int started = 1;
  int resetSW = 0;
  double kappa = 0.02;


  double xMax = 10.0;

  double t = 0.0;
  double dt = 0.01;

  int NNp = 100;
  double xx[] = new double [NNp];
  double vx[] = new double [NNp];
  double ax[] = new double [NNp];

  double hh[] = new double [NNp];
  double mass[] = new double [NNp];
  double density[] = new double [NNp];
  double pressure[] = new double [NNp];
  double energy[] = new double [NNp];
  double power[] = new double [NNp];

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

  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");
    bt_reset.addActionListener(this);

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

    bt_stop = new Button("stop");
    bt_stop.addActionListener(this);
    
    sc_rate= new Scrollbar(Scrollbar.HORIZONTAL,10,10,1,110);
    sc_rate.addAdjustmentListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
    pnl.add(bt_reset);
    pnl.add(bt_start);
    pnl.add(bt_stop);
    pnl.add(sc_rate);
    pnl.add(new Label("  "));
    pnl.add(new Label("  "));
    add(pnl,"North");

    setInitialCondition();
  }

  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;
      resetSW = 1;
      started = 0;
    } else if(ev.getSource() == bt_start){ 
      started = 1;
    } else if(ev.getSource() == bt_stop){
      started = 0;
    }
		}

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if(ev.getSource() == sc_rate) { 
      incRate = sc_rate.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 gb;

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

    gOff.setColor(Color.black);
    gOff.drawRect(30,80,400,200);

    gOff.setColor(Color.blue);
    gOff.drawString("density ni",50,100);
    dotPlot(density, 8.0);
    gOff.setColor(Color.lightGray);
    plotOval();

    gOff.setColor(Color.black);
    gOff.drawString("smoothed particle hydrodynamics - 1D",60,60);

    gOff.drawString("t="+(int)(t*10.0+0.5)/10.0+" ",630/6*0+10,40);
    gOff.drawString("rate="+incRate+" ",630/6*3+10,40);

    gb = 450;
    gOff.setColor(Color.black);
    gOff.drawString("dt = 0.01",gb,100);
    gOff.drawString("ni = sum(mj*Wij,j)",gb,120);
    gOff.drawString("dvi/dt",gb,140);
    gOff.drawString(" = sum(mj*aij*grad(Wij),j)",gb,155);
    gOff.drawString("    (aij = Pj/nj^2+Pi/ni^2)",gb,170);
    gOff.drawString("Pi = const*ni",gb,190);
    gOff.drawString("dxi/dt = vi",gb,210);
    gOff.drawString("xMax=10, h = 0.5",gb,230);
  }

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

  void dotPlot(double f[], double mag) {
    int i,gx,gy;

    gx = 30; gy = 80 + 200;
    for (i=0; i<NNp; i++) {
      gOff.fillOval(gx+(int)(xx[i]*40.0+0.5)-2,gy-(int)(mag*f[i])-2, 4, 4);

    }
  }

  void plotOval() {
    int i,gx,gy,ix;

    gx = 30; gy = 80 + 200;
    for (i=0; i<NNp; i++) {
	  ix = gx+(int)(xx[i]*40.0+0.5);
      gOff.drawOval(gx+(int)(xx[i]*40.0+0.5)-10,gy-10, 20, 20);
    }
  }


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

  void setInitialCondition() {
    int i,j;
    double s;

    t = 0.0;

    for (i=0; i<NNp; i++) {
	  mass[i] = 1.0;
	  xx[i] = 0.05*i;
	  vx[i] = 0.0;
	  ax[i] = 0.0;
	  hh[i] = 0.5;
    }

    for (i=0; i<NNp; i++) {
      energy[i] = 1.0;
	  power[i] = 0.0;
    }

    for (i=0; i<NNp; i++) {
	  s = 0.0;
      for (j=0; j<NNp; j++) {
        s += mass[j]*ww(i,j);
	  }
	  density[i] = s;
    }

  }

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

  void timeEvolution() {
    int i;

    if (resetSW==1) {
      setInitialCondition();
      resetSW = 0;
    }
    if (started==1) {
      for (i=0; i<incRate; i++) {
        timeStep();
      }
    }
  }

  void timeStep() {
    int i,j;
	double dtv2,s;

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

    for (i=0; i<NNp; i++) {
	  s = 0.0;
      for (j=0; j<NNp; j++) {
        s += mass[j]*ww(i,j);
	  }
	  density[i] = s;
    }

    for (i=0; i<NNp; i++) {
	  pressure[i] = 0.1*density[i];
    }

    for (i=0; i<NNp; i++) {
      vx[i] += dtv2*ax[i];
      xx[i] += vx[i]*dt;
    }
    forceCalc();
    for (i=0; i<NNp; i++) {
      vx[i] += dtv2*ax[i];
      energy[i] += power[i]*dt;
    }
	
    for (i=0; i<NNp; i++) {
      if (xx[i] < 0.0) { 
        xx[i] = 0.0; vx[i] = -vx[i];
      }
      if (xx[i] > xMax) {
        xx[i] = xMax; vx[i] = -vx[i];
      }
    }
  }

  void forceCalc() {
    int i,j;
    double fi,wi,ai,aj;

    for (i=0; i<NNp; i++) {
	  ax[i] = 0.0; power[i] = 0.0;
	  ai = pressure[i]/(density[i]*density[i]);
	  for (j=0; j<NNp; j++) {
	    aj = mass[j]*(ai+pressure[j]/(density[j]*density[j]))*dww(i,j);
		ax[i] += -aj;
		power[i] += 0.5*aj*(vx[i]-vx[j]);
	  }
    }
  }

  double distance(double x1, double x2) {
    double ret;

    ret = 0.0;
	if (x1>x2) {
	  ret = x1-x2;
	} else {
	  ret = x2-x1;
	}
    return ( ret );
  }

  double ww(int i, int j) {

    return ( kernel(xx[i]-xx[j], (hh[i]+hh[j])/2.0) );
  }

  double dww(int i, int j) {

    return ( dwwvdr(xx[i]-xx[j], (hh[i]+hh[j])/2.0) );
  }

  double d2ww(int i, int j) {

    return ( d2wwvdr2(xx[i]-xx[j], (hh[i]+hh[j])/2.0) );
  }

  double kernel(double r, double h) {
    double a,q,ret;

    ret = 0.0;
	q = r/h;
	a = 2.0/3.0/h;
	if (q<-2.0) {
      ret = 0.0;
	} else if (q<-1.0) {
      ret = a*0.25*(2.0+q)*(2.0+q)*(2.0+q);
    } else if (q<0.0) {
      ret = a*(1.0-1.5*q*q-0.75*q*q*q);
    } else if (q<1.0) {
      ret = a*(1.0-1.5*q*q+0.75*q*q*q);
    } else if (q<=2.0) {
      ret = a*0.25*(2.0-q)*(2.0-q)*(2.0-q);
	} else {
      ret = 0.0;
	}
    return ( ret );
  }

  double dwwvdr(double r, double h) {
    double a,q,ret;

    ret = 0.0;
	q = r/h;
	a = 2.0/3.0/h/h;
	if (q<-2.0) {
      ret = 0.0;
	} else if (q<-1.0) {
      ret = a*0.75*(2.0+q)*(2.0+q);
    } else if (q<0.0) {
      ret = a*(-3.0*q-2.25*q*q);
    } else if (q<1.0) {
      ret = a*(-3.0*q+2.25*q*q);
    } else if (q<=2.0) {
      ret = -a*0.75*(2.0-q)*(2.0-q);
	} else {
      ret = 0.0;
	}
    return ( ret );
  }

  double d2wwvdr2(double r, double h) {
    double a,q,ret;

    ret = 0.0;
	q = r/h;
	a = 2.0/3.0/h/h/h;
	if (q<-2.0) {
      ret = 0.0;
	} else if (q<-1.0) {
      ret = a*1.5*(2.0+q);
    } else if (q<0.0) {
      ret = a*(-3.0-4.5*q);
    } else if (q<1.0) {
      ret = a*(-3.0+4.5*q);
    } else if (q<=2.0) {
      ret = a*1.5*(2.0-q);
	} else {
      ret = 0.0;
	}
    return ( ret );
  }

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

}