/* plasma - particle-mesh method 2D */
/* coded by Ikeuchi Mitsuru  */
/* ver 0.0.1   2005.10.23    */

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

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

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

  Button bt_reset, bt_startStop;
  Choice ch_view;
  Scrollbar sc_scale;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 30;

  int resetSW = 0;
  int started = 1;
  int dispMode = 5;
  double viewScale = 1.0;


  int dispWidth = 630;
  int dispHeight = 320;

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

  double viewTheta = -15.0*3.14/180.0;
  double viewFai = -72.0*3.14/180.0;
  double dtheta = 0.0*3.14/180.0;
  double pai = 3.1415926536;
  double cosTh = Math.cos(viewTheta); 
  double sinTh = Math.sin(viewTheta);
  double cosFi = Math.cos(viewFai);
  double sinFi = Math.sin(viewFai);

  double t = 0.0;
  double dt = 1.0e-9;

  double dx = 1.0e-3;
  double dy = 1.0e-3;

  int Nsx = 40;
  int Nsy = 40;
  double xMax = Nsx*dx;
  double yMax = Nsy*dy;

  double dispSpan = xMax;
  double dispScale = (250.0/dispSpan);
  double ffn = 1000.0;

  int NNp = Nsx*Nsy*100;
  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];
  double ffx[] = new double [NNp];
  double ffy[] = new double [NNp];

  int section[][][] = new int [Nsx][Nsy][500];
  double charge[][] = new double [Nsx][Nsy];
  double field[][][] = new double [Nsx][Nsy][2];
  double potential[][] = new double [Nsx][Nsy];

  int xpts[] = new int[100];
  int ypts[] = new int[100];

  double amas = 1.66e-27;
  double emas = 9.11e-31;
  double eee = 1.602e-19;
  double particleData[][] = {
    /* mass(kg),charge,d(A), HSBcolor */
    {  1.0 *emas, -eee, 2.0 , 0.66 }, /*  0 electron */
    {  1.0 *amas, +eee, 1.06, 0.90 }, /*  1 H   */
    {  4.0 *amas, +eee, 2.33, 0.05 }, /*  2 He  */
    { 20.18*amas, +eee, 2.77, 0.10 }, /*  3 Ne  */
    { 39.95*amas, +eee, 4.17, 0.15 }, /*  4 Ar  */
    { 83.8 *amas, +eee, 4.76, 0.20 }, /*  5 Kr  */
    {131.3 *amas, +eee, 5.74, 0.25 }, /*  6 Xe  */
    {1.0e20*amas, +eee, 4.00, 0.01 }  /* 7 fixed charge */
  };

/* ----------------------------- 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_view = new Choice();
    ch_view.addItem("particle1000");
    ch_view.addItem("particle all");
    ch_view.addItem("rho(x,y)");
    ch_view.addItem("V(x,y)");
    ch_view.addItem("E(x,y)");
    ch_view.addItem("grid V");
    ch_view.addItem("grid rho");

    ch_view.addItemListener(this);
    ch_view.select("grid V");

    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(bt_reset);
    pnl.add(bt_startStop);
    pnl.add(new Label("  "));
    pnl.add(new Label("  "));
    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 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 itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_view) {
      dispMode = ch_view.getSelectedIndex();
      if (dispMode<=4) {
        dragSW=1;
      } else {
        dragSW=2;
      }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if(ev.getSource() == sc_scale) { 
      viewScale = 0.01*(double)(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;
    }
	
    if (dragSW==1) {
      xShift += (1.0/viewScale/dispScale)*(dgX-dgXb);
      yShift += (1.0/viewScale/dispScale)*(dgY-dgYb);
    } else if (dragSW==2) {
      viewTheta += 0.5*3.14/180.0*(dgX-dgXb);
      viewFai += 0.5*3.14/180.0*(dgY-dgYb);
      cosTh = Math.cos(viewTheta); 
      sinTh = Math.sin(viewTheta);
      cosFi = Math.cos(viewFai);
      sinFi = Math.sin(viewFai);
    }
  }


  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 = viewScale*dispScale;
    gbx = 30+(int)(sc*xShift);
    gby = 60+(int)(sc*yShift);

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

    } else if (dispMode==1) {
      ovalPlot(gbx,gby,sc,NNp);

    } else if (dispMode==2) {
      fieldPlot(charge,gbx,gby,sc,1.0e15);

    } else if (dispMode==3) {
      fieldPlot(potential,gbx,gby,sc,1000.0);

    } else if (dispMode==4) {
      vectorPlot(field,gbx,gby,sc,0.01);

    } else if (dispMode==5) {
      gridPlot(potential,1000.0);

    } else if (dispMode==6) {
      gridPlot(charge,1.0e14);

    }

    gbx = 420;
    dispParameter();
    gOff.drawString("Lx="+(float)(xMax*1000.0)+"(mm)",gbx,80);
    gOff.drawString("Ly="+(float)(yMax*1000.0)+"(mm)",gbx,100);
    gOff.drawString("Ti="+(int)(temperatureOf(1)*10.0+0.5)/10.0+" K",gbx,120);
    gOff.drawString("Te="+(int)(temperatureOf(0)*10.0+0.5)/10.0+" K",gbx,140);
    gOff.drawString("max Nc="+maxSection()+"",gbx,160);
    gOff.drawString("density="+(float)(NNp*ffn/(xMax*yMax))+"(/m2)",gbx,180);
  }


  void dispParameter() {
    gOff.setColor(Color.black);
	gOff.drawString("t="+(int)(t*1.0e11+0.5)/100.0+"(ns)",dispWidth/6*0+10,40);
    if (started==1) {
      gOff.drawString("started",dispWidth/6*1+10,40);
    } else {
      gOff.drawString("stopped",dispWidth/6*1+10,40);
    }
    gOff.drawString("scale="+(int)(viewScale*100.0+0.5)/100.0+"",dispWidth/6*4+10,40);
    gOff.drawString("view",dispWidth/6*5+10,40);
  }

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

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

    gOff.setColor(Color.black);
    gOff.drawRect(gx,gy,(int)(sc*xMax+0.5),(int)(sc*yMax+0.5));

    for (i=0; i<Nsample; i++) {
	  ir = (int)(0.5*particleData[kind[i]][2]+2.0);
      gOff.setColor(Color.getHSBColor((float)particleData[kind[i]][3],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);
    }
  }

  void fieldPlot(double ff[][], 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 = ff[ix][iy]*mag; if (f>0.999) f = 0.999;
		if (f<-0.999) f = -0.999;
        gOff.setColor(Color.getHSBColor((float)(0.33-0.33*f),0.9f,0.9f));
        gOff.fillRect(gx+(int)(ss*ix+0.5),gy+(int)(ss*iy+0.5), iss, iss);
      }
    }
  }

  void vectorPlot(double vf[][][], int gx, int gy, double sc, double mag) {
    int ix,iy,iss,igx,igy,igx2,igy2;
    double ss;

    gOff.setColor(Color.black);
    gOff.drawRect(gx,gy,(int)(sc*xMax+0.5),(int)(sc*yMax+0.5));

    ss = sc*xMax/(double)Nsx;
    iss = (int)(ss+0.999);

    for (ix=0; ix<Nsx; ix++) {
      for (iy=0; iy<Nsy; iy++) {
        igx = gx+(int)(ss*(ix+0.5)+0.5);
        igy = gy+(int)(ss*(iy+0.5)+0.5);

        igx2 = igx+(int)(vf[ix][iy][0]*sc*mag+0.5);
        igy2 = igy+(int)(vf[ix][iy][1]*sc*mag+0.5);
	
        if (vf[ix][iy][0]>0.0) {
          gOff.setColor(Color.getHSBColor(0.66f,0.8f,0.8f));
        } else {
          gOff.setColor(Color.getHSBColor(0.01f,0.8f,0.8f));
        }
        gOff.drawLine(igx,igy,igx2,igy2);
      }
    }
  }

/* -------------------------- grid plot -------------- */

  void gridPlot(double ff[][],double weight) {
    int i,j,k,gbx,gby;
    int cx,cy,cz;
    double z,f,px,py,pz,sc,col;

    cx = Nsx/2;
    cy = Nsy/2;
    cz = 0;
    sc = 8.0;

    gbx = 40+(int)(cx/sc);
    gby = 40+(int)(cy/sc);

    for(j=0; j<Nsy; j += 1) {
      for (i=0; i<Nsx; i++) { 
        f = weight*ff[i][j];
        z = f + cz;
        py = cosFi*(j-cy)+sinFi*(z-cz) + cy; 
        pz = -sinFi*(j-cy)+cosFi*(z-cz) + cz;
        px = cosTh*(i-cx)+sinTh*(pz-cz) + cx; 
        pz = -sinTh*(i-cx)+cosTh*(pz-cz) + cz;
        xpts[i] = (int)(sc*px)+gbx;
        ypts[i] = (int)(sc*py)+gby;
      }
      for (i=0; i<Nsx-1; i++) {
	    if (ff[i][j]+ff[i+1][j]>=0) { col=0.10; } else { col=0.66; };
        gOff.setColor(Color.getHSBColor((float)col,0.9f,0.8f));
        gOff.drawLine(xpts[i],ypts[i],xpts[i+1],ypts[i+1]);
      }
    }

    for(i=0; i<Nsx; i += 1) {
      for (j=0; j<Nsy; j++) {
        f = weight*ff[i][j];
        z = f + cz;
        py = cosFi*(j-cy)+sinFi*(z-cz) + cy; 
        pz = -sinFi*(j-cy)+cosFi*(z-cz) + cz;
        px = cosTh*(i-cx)+sinTh*(pz-cz) + cx; 
        pz = -sinTh*(i-cx)+cosTh*(pz-cz) + cz;
        xpts[j] = (int)(sc*px)+gbx;
        ypts[j] = (int)(sc*py)+gby;
      }
      for (j=0; j<Nsy-1; j++) { 
        if (ff[i][j]+ff[i][j+1]>=0) { col=0.10; } else { col=0.66; };
        gOff.setColor(Color.getHSBColor((float)col,0.9f,0.8f));
        gOff.drawLine(xpts[j],ypts[j],xpts[j+1],ypts[j+1]);
      }
    }
  }


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

  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 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 );
  }

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

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

    t = 0.0;

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

    clearPotential();
	setTemperature(10000, 0);
	setTemperature(10000, 1);
  }

  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() {

    if (resetSW==1) {
      setInitialCondition();
      started = 0;
      resetSW = 0;
    }

    if (started==1) {
      timeStep();
    }
  }

  void timeStep() {
    int i;

    t = t + dt;

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

  }

/*-------------------------------  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++) {

        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;
    }

  }

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

  void setElectricField(int iterMax) {
    int i,j;

    setCharge();
    poisson(iterMax);

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

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

  void setCharge() {
    int i,j,k,ip;

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

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

    h2 = dx*dx;
    for (iter=1; iter<iterMax; iter++) {
      for (i=1; i<Nsx-1; i++) {
        for(j=1;j<Nsy-1;j++) {
          potential[i][j] = potential[i][j]+0.33*(potential[i+1][j]+potential[i-1][j]+potential[i][j+1]+potential[i][j-1]-4.0*potential[i][j]+h2*charge[i][j]/(dy*dx)/(8.854e-12));
        }
      }
    }
  }


  void clearPotential() {
    int i,j;

    for (i=1; i<Nsx-1; i++) {
      for(j=1;j<Nsy-1;j++) {
          potential[i][j] = 0.0;
      }
    }
  }

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

  void movement() {
    int i,ix,iy;
    double dtv2,rr,q,m;

    dtv2 = dt/2.0;
    for (i=0; i<NNp; i++) {
      m = massOf(kind[i]);
	  q = chargeOf(kind[i]);
      vx[i] += dtv2*ffx[i]/m;
      vy[i] += dtv2*ffy[i]/m;
      xx[i] += vx[i]*dt;
      yy[i] += vy[i]*dt;
    }

    for (i=0; i<NNp; i++) {
	  q = chargeOf(kind[i]);
      ix = (int)(Nsx*xx[i]/xMax); if (ix>=Nsx) ix = Nsx-1;
      iy = (int)(Nsy*yy[i]/yMax); if (iy>=Nsy) iy = Nsy-1;
      ffx[i] = field[ix][iy][0]*q;
      ffy[i] = field[ix][iy][1]*q;

    }

    for (i=0; i<NNp; i++) {
	  m = massOf(kind[i]);
      vx[i] += dtv2*ffx[i]/m;
      vy[i] += dtv2*ffy[i]/m;
    }

    rr = 1.0;
    for (i=0; i<NNp; i++) {

        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] = 0.0; vx[i] = rr*vx[i]; vy[i] = -rr*vy[i]; 
        }
        if (yy[i] > yMax) {
          yy[i] = yMax; vx[i] = rr*vx[i]; vy[i] = -rr*vy[i];
        }

    }
  }

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

  double massOf(int kind) {

    return( ffn*particleData[kind][0] );
  }

  double chargeOf(int kind) {

    return( ffn*particleData[kind][1] );
  }

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

}

/** particle in cell (simple version)
 *  (perfect collision-free plasma)
 *
 *  super-particle (represent ffn particles)
 *
 *  1. dividing section (mesh)
 *    registration section[ix][iy][pointer]
 *
 *  2. set electric field 
 *    set charge[ix][iy]
 *    solve Poisson's equation  potential[ix][iy]
 *    electric field field[][] = -grad potential[][]
 *
 *  3. move particle(vericity Verlet algorithm)
 *   v = v + (f/m)*(dt/2)
 *   r = r + v*dt
 *   f = q E
 *   v = v + (f/m)*(dt/2)
 *
 */