/* diffusion - 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 diffusionPM1D 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 = 10;
  int started = 1;
  int resetSW = 0;
  int tempMode = 0;
  double kappa = 0.5;

  double temp0 = 200.0;
  double temp1 = 0.0;


  double xMax = 10.0;

  double t = 0.0;
  double dt = 0.01;

  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 temp[] = new double [NNp];
  double d2temp[] = new double [NNp];
  double dmtemp[] = 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("insulate");
	ch_boundary.addItem("heatBath"); 
    ch_boundary.addItemListener(this);
    ch_boundary.select("insulate");

    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) {
      tempMode = 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("temperature",250,100);


    gOff.setColor(Color.black);
    gOff.drawString("t="+(int)(t*10.0+0.5)/10.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("dT/dt = d^2 T/dx^2",gb,140);
    gOff.drawString("Ti=Ti+sum(mj/nj*(Tj-Ti)*Wij)",gb,160);
    gOff.drawString("d2Ti=sum(mj/nj*(Tj-Ti)*d2Wij)",gb,180);
    gOff.drawString("xMax=10, h = 0.5",gb,200);
  }

/*------------------------- 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)(temp[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;

    for (i=0; i<NNp; i++) {
	  mass[i] = 1.0;
	  xx[i] = /*xMax/NNp*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++) {
	  if (xx[i]<0.5*xMax){
	    temp[i] =  200.0;
	  } else {
        temp[i] =  0.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;

    for (i=0; i<NNp; i++) {
	  d2temp[i] = 0.0;
	  for (j=0; j<NNp; j++) {
	    d2temp[i] += mass[i]/density[j]*(temp[j]-temp[i])*d2ww(i,j);
	  }
    }

    for (i=0; i<NNp; i++) {
	  temp[i] += dt*kappa*d2temp[i];
    }

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

    if (tempMode==1) {
      for (i=0; i<NNp; i++) {
	    if (xx[i]<0.5) {
	      temp[i] = temp0;
	    } else if (xx[i]>xMax-0.5) {
	      temp[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 -------------- */

}