/* diffusion mesh-free differnce method */
/* coded by Ikeuchi Mitsuru   */
/* ver 0.0.1   2005.05.07     */

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

public class diffusionMF1D 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;
  double kappa = 0.01;

  double t = 0.0;
  double dt = 0.1;

  double re0=0.2;
  int NNp = 100;
  double xx[] = new double [NNp];
  double density[] = new double [NNp];
  double lambda[] = new double [NNp];
  double temp[] = new double [NNp]; 
  double wrk[] = new double [NNp];

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

  public void init() {

    resize(630,300);
    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;
      setInitialCondition();
      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);
    dotPlot(density, 10.0);

    gOff.setColor(Color.green);
    dotPlot(lambda, 10000.0);

    gOff.setColor(Color.red);
    dotPlot(temp, 2.0);

    gOff.setColor(Color.black);
    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 = 460;
    gOff.setColor(Color.red);
    gOff.drawString("temperature",gb,100);
    gOff.setColor(Color.blue);
    gOff.drawString("density",gb,120);
    gOff.setColor(Color.green);
    gOff.drawString("lambda",gb,140);
	
    gOff.setColor(Color.black);
	gOff.drawString("dT/dt = D d2T/dx2",gb,200);
    gOff.drawString("dt = 0.1",gb,220);
    gOff.drawString("D = 0.01",gb,240);
    gOff.drawString("L=5.0 (100 points)",gb,260);
  }

/*------------------------- 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]*80.0+0.5)-2,gy-(int)(mag*f[i])-2, 4, 4);

    }
  }

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

  void setInitialCondition() {
    int i;

    t = 0.0;

    for (i=0; i<NNp; i++) {
      xx[i] = 0.05*i;
      if (i<NNp/2) {
        temp[i] = 100.0;
	  } else {
        temp[i] = 0.0;
	  }
    }

	setDensity();
	setLambda();
  }

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

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

  void setLambda() {
    int i,j;
    double s,sw,wij;

    for (i=0; i<NNp; i++) {
	  s = 0.0; sw = 0.0;
	  for (j=0; j<NNp; j++) {
        if (j!=i) {
		  wij = ww(i,j);
	      s += (xx[j]-xx[i])*(xx[j]-xx[i])*wij;
	      sw += wij;
	    }
	  }
	  lambda[i] = s/sw;
    }
  }

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

  void timeEvolution() {
    int i;

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

  void timeStep() {
    int i;

    t = t + dt;

    setLaplacianOf(temp,wrk);

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

  void setLaplacianOf(double aa[], double wk[]) {
    int i,j;
    double s;

    for (i=0; i<NNp; i++) {
	  s = 0.0;
	  for (j=0; j<NNp; j++) {
        if (j!=i) {
	      s += (aa[j]-aa[i])*ww(i,j);
	    }
	  }
	  wk[i] = (2.0/lambda[i]/density[i])*s;
    }
  }

  double ww(int i, int j) {
    double rij;

    rij = xx[i]-xx[j]; if (rij<0.0) rij = -rij;
    return ( weight(rij, re0) );
  }

  double weight(double r, double re) {
    double ret;

    ret = 0.0;
	if (r<re) {
      ret = re/r - 1.0;
	}
    return ( ret );
  }

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

}