/* wave - particle method 1D  */
/* coded by Ikeuchi Mitsuru   */
/* ver 0.0.1   2005.07.23     */

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

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

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

  Button bt_reset, bt_start, bt_stop;
  Choice ch_boundary;
  Scrollbar sc_temp0, sc_temp1;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 30;

  int incRate = 1;
  int started = 1;
  int resetSW = 0;
  int boundaryMode = 0;
  double kappa = 10000.0;

  double temp0 = 200.0;
  double temp1 = 0.0;


  double xMax = 10.0;

  double t = 0.0;
  double dt = 0.01;

  int tBefore = 0;
  int tNow = 1;
  int tNext = 2;

  int NNp = 200;
  double xx[] = new double [NNp];
  double vx[] = new double [NNp];

  double hh[] = new double [NNp];
  double mass[] = new double [NNp];
  double density[] = new double [NNp];

  double aaa[][] = new double [3][NNp];
  double d2aaa[] = new double [NNp];
  double dmaaa[] = 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);

    ch_boundary = new Choice();
    ch_boundary.addItem("free");
	ch_boundary.addItem("fixed"); 
    ch_boundary.addItemListener(this);
    ch_boundary.select("free");

    sc_temp0= new Scrollbar(Scrollbar.HORIZONTAL,200,10,0,210);
    sc_temp0.addAdjustmentListener(this);

    sc_temp1= new Scrollbar(Scrollbar.HORIZONTAL,0,10,0,210);
    sc_temp1.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(ch_boundary);
    pnl.add(sc_temp0);
    pnl.add(sc_temp1);
	/*
    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){ 
      resetSW = 1;
      started = 0;
    } else if(ev.getSource() == bt_start){ 
      started = 1;
    } else if(ev.getSource() == bt_stop){
      started = 0;
    }
  }

  public void itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_boundary) {
      boundaryMode = ch_boundary.getSelectedIndex();
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_temp0) { 
      temp0 = 1.0*sc_temp0.getValue();
    } else if (ev.getSource() == sc_temp1) { 
      temp1 = 1.0*sc_temp1.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.orange);
    densPlot();
	gOff.drawString("density",350,100);

    gOff.setColor(Color.blue);
    tempPlot();
    gOff.drawString("Y",250,100);


    gOff.setColor(Color.black);
    gOff.drawString("t="+(int)(t*100.0+0.5)/100.0+" ",630/6*0+10,40);
    gOff.drawString("boundary",630/6*3+10,40);
    gOff.drawString("T0="+(int)(temp0+0.5)+" ",630/6*4+10,40);
    gOff.drawString("T1="+(int)(temp1+0.5)+" ",630/6*5+10,40);

    gb = 440;
    gOff.setColor(Color.black);
    gOff.drawString("dt = 0.01",gb,100);
    gOff.drawString("ni = sum(mj*Wij)",gb,120);
    gOff.drawString("d^2 Y/dt^2 = K d^2 Y/dx^2",gb,140);
    gOff.drawString("Yi=Yi+sum(mj/nj*(Yj-Yi)*Wij)",gb,160);
    gOff.drawString("d2Yi=sum(mj/nj*(Yj-Yi)*d2Wij)",gb,180);
    gOff.drawString("d^2 Y/dt^2 = ",gb,200);
	gOff.drawString("  {Y(t+k)-2Y(t+k)+Y(t-k)}/k^2",gb,215);
    gOff.drawString("xMax=10, h = 0.5",gb,235);
  }

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

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

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

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

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


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

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

    t = 0.0;
	tBefore = 0;
    tNow = 1;
    tNext = 2;

    for (i=0; i<NNp; i++) {
	  mass[i] = 1.0;
	  xx[i] = xMax*Math.random();
	  hh[i] = 0.5;
    }

    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++) {
	  aaa[tBefore][i] = 0.0;
      aaa[tNow][i] = 0.0;
      aaa[tNext][i] = 0.0;
    }
    for (i=0; i<NNp; i++) {
	  if (xx[i]>0.2*xMax && xx[i]<0.5*xMax){
	    aaa[tBefore][i] = 200.0;
		aaa[tNow][i] = 200.0;
		aaa[tNext][i] = 200.0;
	  }
    }

  }

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

    t = t + dt;
    tBefore += 1; if (tBefore>2) tBefore = 0;
	tNow += 1; if (tNow>2) tNow = 0;
	tNext += 1; if (tNext>2) tNext = 0;

    for (i=0; i<NNp; i++) {
	  d2aaa[i] = 0.0;
	  for (j=0; j<NNp; j++) {
	    d2aaa[i] += mass[i]/density[j]*(aaa[tNow][j]-aaa[tNow][i])*d2ww(i,j);
	  }
    }
    for (i=0; i<NNp; i++) {
	  aaa[tNext][i] = 2.0*aaa[tNow][i] - aaa[tBefore][i]+ dt*dt*kappa*d2aaa[i];
    }

    for (i=0; i<NNp; i++) {
	  dmaaa[i] = 0.0;
	  for (j=0; j<NNp; j++) {
	    dmaaa[i] += mass[j]/density[j]*(aaa[tNext][j]-aaa[tNext][i])*ww(i,j);
	  }
    }
    for (i=0; i<NNp; i++) {
	  aaa[tNext][i] += dmaaa[i];
    }

    if (boundaryMode==1) {
      for (i=0; i<NNp; i++) {
	    if (xx[i]<0.5) {
	      aaa[tNext][i] = temp0;
	    }
	    if (xx[i]>xMax-0.5) {
	      aaa[tNext][i] = temp1;
	    } 
      }
    }

  }


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

}