/* plasma - direct simulation Monte Carlo method for plasmas 1D */
/* coded by Ikeuchi Mitsuru   */
/* ver 0.0.1   2005.10.08     */

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

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

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

  Button bt_reset, bt_startStop;
  Choice ch_kind, ch_kind2, ch_col, ch_view;
  Scrollbar sc_temp, sc_temp2, sc_adsorption, sc_scale;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 50;

  double t = 0.0;
  double dt = 2.0e-10;/* s */

  double dx = 1.0e-3;/* m */
  double dy = 1.0e-3;/* m */
  double vvc = dx*dy;
  int Nsx = 100;
  int Nsy = 10;
  double xMax = Nsx*dx;
  double yMax = Nsy*dy;

  double ffn = 5.0e4;
  double sigmaMax = 2.0*4.19e-10;
  double vvmax2 = 2000.0;

  int dispMode = 12;
  int resetSW = 0;
  int started = 1;
  int incRate = 1;
  double dispScale = (400.0/xMax);
  double scale = 1.0;
  double adsorptionEnergy = 0.1; /* eV */
  int collisionMode = 1;
  int kindOfGas = 1001; /* H+ */
  int kindOfGas2 = 1000; /* electron */
  double setTemp = 10000.0;
  double setTemp2 = 10000.0;

  int dispWidth = 630;
  int dispHeight = 320;

  int dgX, dgY, dgXb, dgYb;
  double xShift = 0.0;
  double yShift = 0.0;

  int NNp = Nsx*Nsy*30;
  int kind[] = new int [NNp];
  double xx[] = new double [NNp];
  double yy[] = new double [NNp];
  double vx[] = new double [NNp];
  double vy[] = new double [NNp];

  int mark[] = new int [NNp];

  double gasData[][] = {
    /* mass(x1.66e-27 kg),  d(x1.0e-10 m), HSBcolor */
    { 0.000549,2.0,0.66 }, /*  0 dummy (electron) */
    {  1.0 , 1.06, 0.90 }, /*  1 H   */
    {  4.0 , 2.33, 0.05 }, /*  2 He  */
    { 20.18, 2.77, 0.10 }, /*  3 Ne  */
    { 39.95, 4.17, 0.15 }, /*  4 Ar  */
    { 83.8 , 4.76, 0.20 }, /*  5 Kr  */
    {131.3 , 5.74, 0.25 }  /*  6 Xe  */
  };

  int section[][][] = new int [Nsx][Nsy][200];
  int cellAttribute[][] = new int [Nsx][Nsy];

  double density[][] = new double [Nsx][Nsy];
  double temperature[][] = new double [Nsx][Nsy];
  double pressure[][] = new double [Nsx][Nsy];

  double charge[] = new double [Nsx];
  double field[] = new double [Nsx];
  double potential[] = new double [Nsx];

  int vDistribution[][] = new int [200][3];
  int vSpace[][] = new int [250][250];

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

  public void init() {

    resize(dispWidth,dispHeight);
    setBackground(Color.white);

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

    bt_reset= new Button("reset");
    bt_reset.addActionListener(this);

    bt_startStop= new Button("start/stop");
    bt_startStop.addActionListener(this); 

    ch_kind = new Choice();
    ch_kind.addItem("H"); ch_kind.addItem("He"); ch_kind.addItem("Ne");
	ch_kind.addItem("Ar"); ch_kind.addItem("Kr"); ch_kind.addItem("Xe");
    ch_kind.addItemListener(this);
    ch_kind.select("H");

    ch_kind2 = new Choice();
    ch_kind2.addItem("H"); ch_kind2.addItem("He"); ch_kind2.addItem("Ne");
	ch_kind2.addItem("Ar"); ch_kind2.addItem("Kr"); ch_kind2.addItem("Xe");
    ch_kind2.addItemListener(this);
    ch_kind2.select("He");

    ch_col = new Choice();
    ch_col.addItem("null");
    ch_col.addItem("rigid sphere");
    ch_col.addItemListener(this);
    ch_col.select("rigid sphere");

    ch_view = new Choice();
    ch_view.addItem("particle");
    ch_view.addItem("collision");
    ch_view.addItem("velocity");
    ch_view.addItem("sample 1000");
    ch_view.addItem("veloc. 1000");
    ch_view.addItem("density");
    ch_view.addItem("mass flow");
    ch_view.addItem("v space");
    ch_view.addItem("temperature");
    ch_view.addItem("pressure");

    ch_view.addItem("charge");
    ch_view.addItem("eField");
	ch_view.addItem("potential");

    ch_view.addItemListener(this);
    ch_view.select("potential");

    sc_temp= new Scrollbar(Scrollbar.HORIZONTAL,10000,10,100,20010);
    sc_temp.addAdjustmentListener(this);

    sc_temp2= new Scrollbar(Scrollbar.HORIZONTAL,10000,10,100,20010);
    sc_temp2.addAdjustmentListener(this);

    sc_adsorption= new Scrollbar(Scrollbar.HORIZONTAL,100,10,0,160);
    sc_adsorption.addAdjustmentListener(this);

    sc_scale= new Scrollbar(Scrollbar.HORIZONTAL,100,10,50,1010);
    sc_scale.addAdjustmentListener(this);

    addMouseListener(this);  
    addMouseMotionListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
/*
    pnl.add(new Label("  "));
*/
    pnl.add(ch_kind);
    pnl.add(bt_startStop);
    pnl.add(sc_temp);
    pnl.add(sc_temp2);
    pnl.add(sc_scale);
    pnl.add(ch_view);
    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 itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_kind) {
      kindOfGas = 1001 + ch_kind.getSelectedIndex();
	  resetSW = 1;
    } else if (ev.getSource() == ch_kind2) {
      kindOfGas2 = 1 + ch_kind2.getSelectedIndex();
	  resetSW = 1;
    } else if (ev.getSource() == ch_view) {
      dispMode = ch_view.getSelectedIndex();
    } else if (ev.getSource() == ch_col) {
      collisionMode = ch_col.getSelectedIndex();
    }
  }

  public void actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_reset){ 
      resetSW = 1;
    } else if (ev.getSource() == bt_startStop){
	  if (started==0) {
        started = 1;
	  } else {
        started = 0;
	  }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_temp) { 
      setTemp = 1.0*sc_temp.getValue();
      setTemperature(setTemp,kindOfGas);
    } else if (ev.getSource() == sc_temp2) { 
      setTemp2 = 1.0*sc_temp2.getValue();
      setTemperature(setTemp2,kindOfGas2);
    } else if (ev.getSource() == sc_adsorption) { 
      adsorptionEnergy = 0.001*sc_adsorption.getValue();
    } else if (ev.getSource() == sc_scale) { 
      scale = 0.01*sc_scale.getValue();
    }
  }

  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;
    }
    xShift += (1.0/scale/dispScale)*(dgX-dgXb);
    yShift += (1.0/scale/dispScale)*(dgY-dgYb);
  }

  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 gx, gbx, gby;
    double sc;

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

    sc = scale*dispScale;
    gbx = 15+(int)(sc*xShift);
    gby = 80+(int)(sc*yShift);

    gOff.setColor(Color.black);
	gOff.drawRect(gbx,gby,(int)(sc*xMax+0.5),(int)(sc*yMax+0.5));
    gOff.setColor(Color.gray);
    gOff.drawLine(gbx,gby,gbx+(int)(sc*xMax+0.5),gby);
    gOff.drawLine(gbx,gby+(int)(sc*yMax+0.5),gbx+(int)(sc*xMax+0.5),gby+(int)(sc*yMax+0.5));

    if (dispMode==0) {
	  drawBlock(gbx,gby,sc);
      ovalPlot(gbx,gby,sc,NNp);

	} else if (dispMode==1) {
	  drawBlock(gbx,gby,sc);
	  markPlot(gbx,gby,sc);

	} else if (dispMode==2) {
	  drawBlock(gbx,gby,sc);
      velosityPlot(gbx,gby,sc,NNp,20.0);

	} else if (dispMode==3) {
	  drawBlock(gbx,gby,sc);
      ovalPlot(gbx,gby,sc,1000);

	} else if (dispMode==4) {
	  drawBlock(gbx,gby,sc);
      velosityPlot(gbx,gby,sc,1000,20.0);

	} else if (dispMode==5) {
      densityPlot(gbx,gby,sc);
      drawBlock(gbx,gby,sc);

	} else if (dispMode==6) {
	  drawBlock(gbx,gby,sc);
      massflowPlot(gbx,gby,sc,5.0e3);

	} else if (dispMode==7) {
	  vSpacePlot(gbx,gby,0.1);

	} else if (dispMode==8) {
      setField();
	  fieldPlot(temperature,gbx,gby,sc,0.0000005);
      drawBlock(gbx,gby,sc);
	
	} else if (dispMode==9) {
      setField();
      fieldPlot(pressure,gbx,gby,sc,2.0e4);
	  drawBlock(gbx,gby,sc);

	} else if (dispMode==10) {
	  drawBlock(gbx,gby,sc);
      ovalPlot(gbx,gby,sc,1000);
      field1DPlot(charge, 15, 120, 0.5/1.6e-19);

	} else if (dispMode==11) {
	  drawBlock(gbx,gby,sc);
      ovalPlot(gbx,gby,sc,1000);
      field1DPlot(field, 15, 120, 0.5);

	} else if (dispMode==12) {
	  drawBlock(gbx,gby,sc);
      ovalPlot(gbx,gby,sc,1000);
      field1DPlot(potential, 15, 120, 100.0);

	} else if (dispMode==13) {
	  ;

	} else if (dispMode==14) {
	  ;

	}

    vDistributionPlot(420,110,0.2);


    gOff.setColor(Color.black);
    gOff.drawString("t="+(int)(t*1.0e9*10.0+0.5)/10.0+" ns",630/6*0+10,40);

	if (started==1) {
      gOff.drawString("started",630/6*1+10,40);
	} else {
      gOff.drawString("stopped",630/6*1+10,40);
	}

    gOff.drawString("setT="+(int)(setTemp+0.5)+" K",630/6*2+10,40);
    gOff.drawString("setTe="+(int)(setTemp2+0.5)+" K",630/6*3+10,40);

    gOff.drawString("scale="+(int)(scale*100.0)/100.0+" ",630/6*4+10,40);
    gOff.drawString("view",630/6*5+10,40);

    gOff.drawString("max Nc="+maxSection()+"",630/6*0+10,60);
    gOff.drawString("T="+(int)(temperatureOf(kindOfGas)*10.0+0.5)/10.0+" K",630/6*2+10,60);
    gOff.drawString("Te="+(int)(temperatureOf(kindOfGas2)*10.0+0.5)/10.0+" K",630/6*3+10,60);
    gOff.drawString("Lx="+(float)(xMax*1000.0)+"(mm), Ly="+(float)(yMax*1000.0)+"(mm)",630/6*4+10,60);


    gOff.setColor(Color.black);
    gx = 450;
    gOff.drawString("N="+NNp+", FN="+(float)ffn+"",gx,80);
    gOff.drawString("sigmaMax="+(float)(sigmaMax)+"(m^2)",gx,100);
    gOff.drawString("vmax="+(float)(vvmax2/2.0)+"(m/s)",gx,120);
    gOff.drawString("Etot="+(float)totalKineticEnergy()+"(J)",gx,140);

  }

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

  void ovalPlot(int gx, int gy, double sc,int Nsample) {
    int i,ir;

    for (i=0; i<Nsample; i++) {
	  ir = (int)(0.5*diameterOf(kind[i])*scale+2.0);
	  if (kind[i]>0) {
	    if (kind[i]>1000) {
		  gOff.setColor(Color.getHSBColor((float)colorOf(kind[i]),0.9f,0.9f));
          gOff.drawOval(gx+(int)(xx[i]*sc+0.5)-ir/2,gy+(int)(yy[i]*sc+0.5)-ir/2, ir, ir);
		} else {
	      gOff.setColor(Color.getHSBColor((float)colorOf(kind[i]),0.9f,0.9f));
          gOff.fillOval(gx+(int)(xx[i]*sc+0.5)-ir/2,gy+(int)(yy[i]*sc+0.5)-ir/2, ir, ir);
		}
      } else if (kind[i]<0) {
	    gOff.setColor(Color.getHSBColor((float)colorOf(kind[i]),0.9f,0.4f));
        gOff.drawOval(gx+(int)(xx[i]*sc+0.5)-ir/2,gy+(int)(yy[i]*sc+0.5)-ir/2, ir, ir);
	  }
    }
  }

  void markPlot(int gx, int gy, double sc) {
    int i,ir;

    ir = 2;
    for (i=0; i<NNp; i++) {
	  if (mark[i]==0) {
	    gOff.setColor(Color.green);
	  } else {
        gOff.setColor(Color.red);
	  }
      gOff.drawOval(gx+(int)(xx[i]*sc+0.5)-ir/2,gy+(int)(yy[i]*sc+0.5)-ir/2, ir, ir);
    }
  }

  void velosityPlot(int gx, int gy, double sc, int Nsample, double mag) {
    int i,ix,iy,ix2,iy2;

    for (i=0; i<Nsample; i++) {
      if (kind[i]>0) {
        ix = gx+(int)(xx[i]*sc+0.5);
        iy = gy+(int)(yy[i]*sc+0.5);
        if (vx[i]>0.0) {
          gOff.setColor(Color.getHSBColor(0.66f,0.8f,0.8f));
        } else {
          gOff.setColor(Color.getHSBColor(0.01f,0.8f,0.8f));
        }

        ix2 = ix+(int)(vx[i]*dt*sc*mag+0.5);
        iy2 = iy+(int)(vy[i]*dt*sc*mag+0.5);
        gOff.drawLine(ix,iy,ix2,iy2);
      }
    }
  }

  void densityPlot(int gx, int gy, double sc) {
    int ix,iy,iss;
    double ss,d;
	
    ss = sc*xMax/(double)Nsx;
    iss = (int)(ss+0.999);
    for (ix=0; ix<Nsx; ix++) {
      for (iy=0; iy<Nsy; iy++) {
        d = 0.015*section[ix][iy][0]; if (d>0.999) d = 0.999;
        gOff.setColor(Color.getHSBColor(0.33f,0.8f,(float)(d)));
        gOff.fillRect(gx+(int)(ss*ix+0.5),gy+(int)(ss*iy+0.5), iss, iss);
      }
    }
  }

  void massflowPlot(int gx, int gy, double sc,double mag) {
    int ix,iy,ip,iss,ii,igx,igy,igx2,igy2;
    double ss,m,mvx,mvy;
	
    ss = sc*xMax/(double)Nsx;
    iss = (int)(ss+0.999);

    for (ix=0; ix<Nsx; ix++) {
      for (iy=0; iy<Nsy; iy++) {
        mvx = 0.0; mvy = 0.0;
        for (ip=1; ip<=section[ix][iy][0]; ip++) {
          ii = section[ix][iy][ip];
          m = massOf(kind[ii])/(29.0*1.66e-27);
          mvx += m*vx[ii]; mvy += m*vy[ii];
        }
        igx = gx+(int)(ss*(ix+0.5)+0.5);
        igy = gy+(int)(ss*(iy+0.5)+0.5);

        if (mvx>0.0) {
          gOff.setColor(Color.getHSBColor(0.66f,0.8f,0.8f));
        } else {
          gOff.setColor(Color.getHSBColor(0.01f,0.8f,0.8f));
        }
        igx2 = igx+(int)(mvx*dt*sc*mag+0.5);
        igy2 = igy+(int)(mvy*dt*sc*mag+0.5);
        gOff.drawLine(igx,igy,igx2,igy2);
      }
    }
  }

  void fieldPlot(double fld[][], int gx, int gy, double sc, double mag) {
    int ix,iy,iss;
    double ss,f;
	
    ss = sc*xMax/(double)Nsx;
    iss = (int)(ss+0.999);
    for (ix=0; ix<Nsx; ix++) {
      for (iy=0; iy<Nsy; iy++) {
        f = fld[ix][iy]*mag; if (f>0.999) f = 0.999;
        gOff.setColor(Color.getHSBColor((float)(0.66-0.66*f),0.9f,0.9f));
        gOff.fillRect(gx+(int)(ss*ix+0.5),gy+(int)(ss*iy+0.5), iss, iss);
      }
    }
  }

  void vSpacePlot(int gx, int gy, double mag) {
    int i,j,ip;
    double v,f;

    for(i=0;i<250;i++) {
      for(j=0;j<250;j++) {
	    vSpace[i][j] = 0;
	  }
    }

    for(ip=0;ip<NNp;ip++) {
      if (kind[ip]>0) {
        i = (int)((vx[ip]+1250000.0)/10000.0+0.5);
        if (i<0) { 
          i = 0;
        } else if (i>249) {
          i = 249;
        }
        j = (int)((vy[ip]+1250000.0)/10000.0+0.5);
        if (j<0) { 
          j = 0;
        } else if (j>249) {
          j = 249;
        }

        vSpace[i][j] += 1;
	  }
    }

    gOff.setColor(Color.getHSBColor(0.66f,0.9f,0.9f));
    gOff.fillRect(gx,gy,250,250);

    for (i=0; i<250; i++) {
      for (j=0; j<250; j++) {
        f = vSpace[i][j]*mag; if (f>0.999) f = 0.999;
        if (f>0.0) {
          gOff.setColor(Color.getHSBColor((float)(0.66-0.66*f),0.9f,0.9f));
          gOff.fillRect(gx+i,gy+j,1,1);
        }
      }
    }
    gOff.setColor(Color.yellow);
    gOff.drawLine(gx,gy+125,gx+250,gy+125);
    gOff.drawLine(gx+125,gy,gx+125,gy+250);
  }

  void field1DPlot(double ff[], int gx, int gy, double yMag) {
    int i,n;
    double col;

    for (i=0; i<Nsx; i++) {
	  n = (int)(yMag*ff[i]+100.0);
	  if (n>100) {
	    col = 0.01;
	    gOff.setColor(Color.getHSBColor((float)col,0.7f,0.8f));
	    gOff.drawRect(gx+i*4,gy+100-(n-100),4,(n-100));
      } else {
	    col = 0.66;
	    gOff.setColor(Color.getHSBColor((float)col,0.7f,0.8f));
	    gOff.drawRect(gx+i*4,gy+100,4,-(n-100));
	  }
    }
  }

  void vDistributionPlot(int gx, int gy, double yMag) {
    int i,n1,n2;;
    double vdiv,c1,c2,cm;

    vdiv = 200.0;
	if (kindOfGas%1000>2) vdiv = 100.0;
	if (kindOfGas%1000>4) vdiv = 50.0;

    setVelocityDistribution(vdiv,10000.0);

    c1 = colorOf(kindOfGas);
    c2 = colorOf(kindOfGas2);
    cm = 0.5*(c1+c2);
    for (i=0; i<200; i++) {
	  n1 = (int)(yMag*vDistribution[i][1]);
	  n2 = (int)(yMag*vDistribution[i][2]);
	  if (n1>=n2) {
	    gOff.setColor(Color.getHSBColor((float)cm,0.9f,0.6f));
        gOff.drawLine(gx+i,gy+200,gx+i,gy+200-n2);
	    gOff.setColor(Color.getHSBColor((float)c1,0.9f,0.9f));
        gOff.drawLine(gx+i,gy+200-n2-1,gx+i,gy+200-n1);
	  } else {
	    gOff.setColor(Color.getHSBColor((float)cm,0.9f,0.6f));
        gOff.drawLine(gx+i,gy+200,gx+i,gy+200-n1);
	    gOff.setColor(Color.getHSBColor((float)c2,0.9f,0.9f));
        gOff.drawLine(gx+i,gy+200-n1-1,gx+i,gy+200-n2);
	  }
    }
	gOff.setColor(Color.getHSBColor((float)c1,0.9f,0.9f));
    gOff.drawString(" "+(int)(vdiv*10.0+0.5)/10.0+" m/s/div",gx+100,gy+80);
    gOff.setColor(Color.getHSBColor((float)c2,0.9f,0.9f));
    gOff.drawString("10000 m/s/div",gx+100,gy+100);
  }

  void drawBlock(int gx, int gy, double sc) {
    int ix,iy,iss,ca;
    double ss;

    ss = sc*xMax/(double)Nsx;
    iss = (int)(ss+0.999);
	
    for (ix=0; ix<Nsx; ix++) {
      for (iy=0; iy<Nsy; iy++) {
	    ca = cellAttribute[ix][iy];
        if (ca>0) {
		  if (ca==1) {
		    gOff.setColor(Color.getHSBColor(0.85f,0.9f,0.9f));
		  } else if (ca==2) {
		    gOff.setColor(Color.getHSBColor(0.85f,0.9f,0.2f));
		  } else if (ca==3) {
		    gOff.setColor(Color.getHSBColor(0.85f,0.9f,0.5f));
		  }
          gOff.fillRect(gx+(int)(ss*ix+0.5),gy+(int)(ss*iy+0.5), iss, iss);
        }
      }
    }
  }

/*-------------------------- statistics -------- */

  double totalKineticEnergy() {
    int i;
    double tke;

    tke = 0.0;
    for(i=0;i<NNp;i++) {
      if (kind[i]>0) {
        tke += 0.5*massOf(kind[i])*(vx[i]*vx[i]+vy[i]*vy[i]);
	  }
    }
	return( ffn*tke );
  }

  double temperatureOf(int knd) {
    int i,n;
    double tke;

    n=0; tke = 0.0;
    for(i=0;i<NNp;i++) {
      if (kind[i]==knd) {
	    n += 1;
        tke += (vx[i]*vx[i]+vy[i]*vy[i]);
	  }
    }
	return( (0.5*massOf(knd)*tke/n)/1.38e-23 );
  }

  int numberOfPartickes(int knd) {
    int i,n;

    n = 0;
    for(i=0;i<NNp;i++) {
      if (kind[i]==knd) {
        n += 1;
	  }
    }
	return( n );
  }

  void setVelocityDistribution(double vscale1,double vscale2) {
    int i,iv;
    double v;

    for(i=0;i<200;i++) {
      vDistribution[i][1] = 0;
      vDistribution[i][2] = 0;
    }


    for(i=0;i<NNp;i++) {
      if (kind[i]>0) {
        v = Math.sqrt(vx[i]*vx[i]+vy[i]*vy[i]);
		if (kind[i]==kindOfGas) {
		  iv = (int)(v/vscale1); if (iv>199) iv = 199;
          vDistribution[iv][1] += 1;
		} else if (kind[i]==kindOfGas2) {
          iv = (int)(v/vscale2); if (iv>199) iv = 199;
		  vDistribution[iv][2] += 1;
		}
	  }
    }
  }

  void setField() {
    int i,j,k,ip,n;
    double ke,vmx,vmy;

    for(i=0;i<Nsx;i++) {
      for(j=0;j<Nsy;j++) {
        n = section[i][j][0];
        density[i][j] = n*ffn/vvc;

        temperature[i][j] = 0.0;
        if (n>10) {
		  vmx = 0.0; vmy = 0.0;
          for(ip=1;ip<=section[i][j][0];ip++) {
            k = section[i][j][ip];
			vmx += vx[k]; vmy += vy[k];
          }
		  vmx = vmx/n;
		  vmy = vmy/n;

          ke = 0.0; 
          for(ip=1;ip<=section[i][j][0];ip++) {
            k = section[i][j][ip];
            ke += 0.5*massOf(kind[k])*((vx[k]-vmx)*(vx[k]-vmx)+(vy[k]-vmy)*(vy[k]-vmy));
          }
          temperature[i][j] = ke/n/(1.0*1.38e-23);
        }

        pressure[i][j] =  density[i][j]*1.38e-23*temperature[i][j];
      }
    }
  }


/*-------------------------- set initial condition --------- */
/*  cellAttribute -  0:space, 1:wall, 2:absorb, 3:adsorb     */
/*  kind - 1-999:molecule, 1000:electron, 1001-1999: ionized */
/*    0:killed, kind<0:adsorped particle (-kind means above) */

  void setInitialCondition() {
    int i,j,ix,iy;

    t = 0.0;

    for(i=0;i<Nsx;i++) {
      for(j=0;j<Nsy;j++) {
        cellAttribute[i][j] = 0;
		/*
        if (i==Nsx*8/10+1 && j>5) cellAttribute[i][j] = 1;
		if (i==Nsx*8/10+10 && j<5) cellAttribute[i][j] = 2;
        if (i==Nsx*8/10+15 && j>5) cellAttribute[i][j] = 3;
		*/
      }
    }

    for (i=0; i<NNp; i++) {
	  if (i%2==0) {
        kind[i] = kindOfGas;
        xx[i] = xMax*(Math.random()*0.999);
        yy[i] = yMax*Math.random();
	  } else {
	    kind[i] = kindOfGas2;
	    xx[i] = xMax*(Math.random()*0.999);
	    yy[i] = yMax*Math.random();
	  }
      vx[i] = 1000.0*(Math.random()-0.5);
      vy[i] = 1000.0*(Math.random()-0.5);
    }

    setSigmaMax();
    setTemperature(setTemp,kindOfGas);
    setTemperature(setTemp2,kindOfGas2);
  }

  void setSigmaMax() {
    int i;
    double s,sMax;

    sMax = 0.0;
    for(i=0;i<NNp;i++) {
	  if (kind[i]>0) {
        s = sigmaOf(kind[i]);
        if (s>sMax) sMax = s;
	  }
    }
    sigmaMax = 2.0*sMax;
  }

  void setTemperature(double temp0, int knd) {
    int i,n;
    double m,tke,tmp,r;

    m = massOf(knd);
    n = 0; tke = 0.0;
    for(i=0;i<NNp;i++) {
      if (kind[i]==knd) {
	    n += 1;
        tke += 0.5*m*(vx[i]*vx[i]+vy[i]*vy[i]);
	  }
    }
	tmp = tke/n/1.38e-23;
	r = Math.sqrt(temp0/(tmp+1.0e-20));
    for(i=0;i<NNp;i++) {
      if (kind[i]==knd) {
	    vx[i] = r*vx[i];
	    vy[i] = r*vy[i];
	  }
    }
  }

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

  void timeEvolution() {
    int i;

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


  void timeStep() {
    int i;

    t = t + dt;

	setElectricField();
    movement();
	dividingSection();

    for(i=0;i<NNp;i++) {
      mark[i] = 0;
    }
	if (collisionMode>0) {
      interaction();
    }
  }


/*--------------------------  set electric field   -------------*/

  void setElectricField() {
    int i,j,k,ip;
    double a,ch;

    for(i=0;i<Nsx;i++) {
      charge[i] = 0.0;
      for(j=0;j<Nsy;j++) {
        for(ip=1;ip<=section[i][j][0];ip++) {
          k = section[i][j][ip];
          charge[i] += chargeOf(kind[k]);
        }
      }
    }

    poisson(100);

	field[0] = -(potential[1]-potential[0])/dx;
    for (i=1; i<Nsx-1; i++) {
      field[i] = -(potential[i+1]-potential[i-1])/(2.0*dx);
    }
	field[Nsx-1] = -(potential[Nsx-1]-potential[Nsx-2])/dx;
  }

  void poisson(int iterMax) {
    int i,j,iter;
    double h2;

    h2 = dx*dx;
	potential[0] = 0.0; potential[Nsx-1] = 0.0;
    for (iter=1; iter<iterMax; iter++) {
      for (i=1; i<Nsx-1; i++) {
        potential[i] = potential[i]+0.66*(potential[i+1]+potential[i-1]-2.0*potential[i]+h2*charge[i]*ffn/(Nsy*dy*dx)/(8.854e-12));
      }
    }
  }

/*-------------------------------  movement   -------------*/

  void movement() {
    int i,ix,ca;
    double rr;

    rr = 1.0;

    for (i=0; i<NNp; i++) {
	  if (kind[i]>0) {
	    ix = (int)(Nsx*xx[i]/xMax); if (ix>=Nsx) ix = Nsx-1;
        vx[i] += field[ix]*chargeOf(kind[i])/massOf(kind[i])*dt;
        xx[i] += vx[i]*dt;
        yy[i] += vy[i]*dt;

        ca = cellAttributeAt(xx[i],yy[i]);
        if (ca>0) {
          if (ca==1) {
            wallCell(i);
		  } else if (ca==2) {
            absorbCell(i);
		  } else if (ca==3) {
            adsorbCell(i);
		  }
        }
	  } else if (kind[i]<0) {
	    if (Math.random()<1.0*Math.exp(-adsorptionEnergy/0.025)) {
		  kind[i] = -kind[i];
		  vx[i] = -vx[i]; vy[i] = -vy[i];
		}
	  }
    }

    rr = 1.0;
    for (i=0; i<NNp; i++) {
	  if (kind[i]>0) {
        if (xx[i] < 0.0) { 
          xx[i] = 0.0; vx[i] = -rr*vx[i]; vy[i] = rr*vy[i]; 
        }
        if (xx[i] > xMax) {
          xx[i] = xMax; vx[i] = -rr*vx[i]; vy[i] = rr*vy[i]; 
        }
        if (yy[i] < 0.0) { 
          yy[i] = yy[i] + yMax; 
        }
        if (yy[i] > yMax) {
          yy[i] = yy[i] - yMax;
        }
      }
    }
  }

  int cellAttributeAt(double x, double y) {
    int ix,iy;

    ix = (int)(Nsx*x/xMax); if (ix>=Nsx) ix = Nsx-1;
    iy = (int)(Nsy*y/yMax); if (iy>=Nsy) iy = Nsy-1;
    return( cellAttribute[ix][iy] );
  }

  void wallCell(int i) {
    int at10, at01;

    xx[i] -= vx[i]*dt;
    yy[i] -= vy[i]*dt;

    at10 = cellAttributeAt(xx[i]+vx[i]*dt, yy[i]);
    at01 = cellAttributeAt(xx[i], yy[i]+vy[i]*dt);
    if (at10==1) vx[i] = -vx[i];
    if (at01==1) vy[i] = -vy[i];
    if (at10==0 && at01==0) {
      vx[i] = -vx[i]; vy[i] = -vy[i];
	}
  }

  void absorbCell(int i) {
    kind[i] = 0;
  }

  void adsorbCell(int i) {
    kind[i] = -kind[i];
    xx[i] -= vx[i]*dt;
    yy[i] -= vy[i]*dt;
  }

/*-------------------------------  dividing section   -------------*/

  void dividingSection() {
    int i,j,ip,iq;

    for(i=0;i<Nsx;i++) {
      for(j=0;j<Nsy;j++) {
        section[i][j][0] = 0;
      }
    }

    for(ip=0;ip<NNp;ip++) {
	  if (kind[ip]>0) {
        i = (int)(Nsx*xx[ip]/xMax); if (i>=Nsx) i = Nsx-1;
        j = (int)(Nsy*yy[ip]/yMax); if (j>=Nsy) j = Nsy-1;
        iq = section[i][j][0]+1;
        section[i][j][0] = iq;
        section[i][j][iq] = ip;
      }
    }
  }

  int maxSection() {
    int i,j,m;

    m = 0;
    for(i=0;i<Nsx;i++) {
      for(j=0;j<Nsy;j++) {
        if (section[i][j][0]>m) m = section[i][j][0];
      }
    }
    return ( m );
  }

/*-------------------------------  interaction   -------------*/

  void interaction() {
    int ic,jc,ii,ir,i,j,nn,immc,pi,pj;
    double vmax,mmc;

    vvmax2 = 2.0*maxSpeed();

    for(ic=0;ic<Nsx;ic++) {
      for(jc=0;jc<Nsy;jc++) {
        nn = section[ic][jc][0];
        if (nn>1) {
          mmc = 0.5*nn*nn*ffn*sigmaMax*vvmax2*dt/vvc;
		  immc = (int)mmc;
          for(ii=0;ii<immc;ii++) {
		    collisionInTheCell(ic, jc);
          }
		  if ( Math.random()<(mmc-immc) ) {
		    collisionInTheCell(ic, jc);
		  }
        }
      }
    }
  }

  double maxSpeed() {
    int i;
    double v2,vm2;

    vm2 = 0.0;
    for(i=0;i<NNp;i++) {
	  if (kind[i]>0) {
        v2 = vx[i]*vx[i]+vy[i]*vy[i];
        if (v2>vm2) vm2 = v2;
	  }
    }
    return( Math.sqrt(vm2) );
  }

  void collisionInTheCell(int ic, int jc) {
    int i,j,pi,pj,nnc;

    nnc = section[ic][jc][0];
    do {
      i = (int)(nnc*Math.random()); if (i>=nnc) i = nnc-1;
      j = (int)(nnc*Math.random()); if (j>=nnc) j = nnc-1;
    } while (i==j);
    pi = section[ic][jc][i+1]; pj = section[ic][jc][j+1];
    collision(pi,pj);
  }


  void collision(int i, int j) {
    double vrel,a,b,vcmx,vcmy,th,costh,sinth;

    vrel = Math.sqrt((vx[i]-vx[j])*(vx[i]-vx[j])+(vy[i]-vy[j])*(vy[i]-vy[j]));
    if ( (vrel/vvmax2)> Math.random() && (clossSection(kind[i], kind[j], vrel)/sigmaMax) > Math.random() ) {
      a = massOf(kind[i])/( massOf(kind[i])+massOf(kind[j]) ); b = 1.0 - a;
      vcmx = a*vx[i] + b*vx[j]; vcmy = a*vy[i] + b*vy[j];
      th = 2.0*3.14159265*Math.random();
      costh = Math.cos(th); sinth = Math.sin(th);
      vx[i] = vcmx + vrel*costh*b;
      vx[j] = vcmx - vrel*costh*a;
      vy[i] = vcmy + vrel*sinth*b;
      vy[j] = vcmy - vrel*sinth*a;
	  mark[i] = 1; mark[j] = 1;
    }
  }

  double clossSection(int kindi, int kindj, double vrel) {

    if (kindi!=1000 && kindj!=1000) {
      return( gasData[kindi%1000][1]+gasData[kindj%1000][1] );
	} else if (kindi==1000 && kindj==1000) {
      return( 2.0+2.0 );
	}

    if (kindi==1000) {
      return( gasData[kindj%1000][1]+0.5 );
	} else {
      return( gasData[kindi%1000][1]+0.5 );
	}
  }

/*-------------------- property of particle kind  -------------*/

  double sigmaOf(int kind) {

    if (kind!=1000) {
      return( gasData[kind%1000][1]*1.0e-10 );
	}
    
	return( 3.0e-10 );
  }

  double massOf(int kind) {

    if (kind!=1000) {
      return( gasData[kind%1000][0]*1.66e-27 );
	}

    return( 9.11e-31 );
  }

  double chargeOf(int kind) {

    if (kind<1000) {
      return( 0.0 );
	}
    if (kind==1000) {
      return( -1.602e-19 );
	} else {
      return( 1.602e-19 );
	}
  }

  double diameterOf(int kind) {

    if (kind!=1000) {
      return( gasData[kind%1000][1] );
	}
    
	return( 0.1 );
  }

  double colorOf(int kind) {

    if (kind!=1000) {
      return( gasData[kind%1000][2] );
	}
    
	return( 0.66 );
  }

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

}

/** Poisson solver (SOR method)
 *
 *  -div(grad V)) + mue dV/dt = rho/eps
 *  -(V(i+1)+V(i-1)-2V(i))/dx^2 + mue (V(t+dt)-V(t))/dt = rho/eps
 *  mue (V(t+dt)-V(t))/dt = (V(i+1)+V(i-1)-2V(i))/dx^2 + rho/eps
 *  V(t+dt) = V(t) + (dt/mue) (V(i+1)+V(i-1)-2V(i))/dx^2 + rho/eps
 *  V(t+dt) = V(t) + (1/mue)(dt/dx^2) { V(i+1)+V(i-1)-2V(i) + rho/eps*dx^2 }
 *
 *
 *  cell size; Lc = 0.001 (m)
 *  mean free path; lambda ~ 0.002 (m)
 *  No. of representative particle in the cell; Nc ~ 30 (particles)
 *  maximum relative speed; Vmax ~ 10^6 (m/s) for electron
 *  time interval; dt = 0.2*Lc/Vmax = 2.0e-10 (s)
 *  maximum cross section of the collision; sgmMax = 2D = 2*4.2e-10 (m)
 *  number density; n = Nc*FN/Lc^2, lambda = 1/(n*sgmMax)
 *    n = 1/(lambda*sgmMax) = 1/(0.002*8.4e-10) = 6.0e11 (particles/m^2)
 *  represent numbers of the real molecule; FN
 *    FN = n*Lc^2/Nc = 6e11*1e-6/30 = 2e4 (particles) --> 5.0e4 (sgm<sgmMax,V<Vmax)
 *  the number of collisions in the cell; Mc
 *    Mc = 0.5*Nc*Nc*FN*sgmMax*Vmax*dt/Lc^2 ~ 4(times)
 */