/** applet No. 1252
 *
 * p-n diode - Monte-Carlo simulation 3D
 * - multi thread  MCS3D :: display - asynchronous
 *
 * Created by Ikeuchi Mitsuru on January 20 2008.
 * Copyright (c) 2008 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2008.01.20  created
 *
 */

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

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

  // ------------------------------------ preset field -----------

  Thread th = null; // for run()-paint() thread
  Thread thmcs = null; // for MCS3D thread
  MCS3D mcs = new MCS3D(); // create MCS3D instance

  // for event
  Choice ch_dir, ch_view;
  Button bt_reset, bt_startStop, bt_step;
  Scrollbar sc_gate, sc_depth;

  // for off-paint buffer
  Dimension dim;
  Image imgOff;
  Graphics gOff;

  // display
  private int sleepTime = 30; // sleep time (ms)
  private int dispWidth = 630;
  private int dispHeight = 460;

  private int dispMode = 0;
  private double boxDispWidth = 320.0;

  private double dispScale = 1.0;
  private double viewScale = 1.0;
  private double viewDepth = 0.5;


  private double xShift = 0.0;
  private double yShift = 0.0;
  private int gx0 = 60;
  private int gy0 = 90;
  private int xImageLoc = gx0;
  private int yImageLoc = gy0;

  private double viewTheta = -15.0*3.14/180.0;
  private double viewFai = -60.0*3.14/180.0;
  private double cosTh = Math.cos(viewTheta); 
  private double sinTh = Math.sin(viewTheta);
  private double cosFi = Math.cos(viewFai);
  private double sinFi = Math.sin(viewFai);

  private int xpts[] = new int[400];
  private int ypts[] = new int[400];
  private double fpts[] = new double[400];
  int dgX, dgY, dgXb,dgYb; // mouse

  private int thCount = 0;


  // ------------------------------ applet main thread -----------

  public void init() {

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

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

    ch_dir = new Choice();
    ch_dir.add("positive");
    ch_dir.add("0"); 
    ch_dir.add("negative");
    ch_dir.addItemListener(this);
    ch_dir.select("0"); 

    ch_view = new Choice();
    ch_view.add("3D-plot");
    ch_view.add("z-depth");
    ch_view.add("grid pot");
    ch_view.addItemListener(this);
    ch_view.select("3D plot");

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

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

    bt_step = new Button("step");
    bt_step.addActionListener(this); 

    sc_gate = new Scrollbar(Scrollbar.HORIZONTAL,200,10,0,510);
    sc_gate.addAdjustmentListener(this);

    sc_depth = new Scrollbar(Scrollbar.HORIZONTAL,50,10,0,110);
    sc_depth.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(bt_step);
    pnl.add(ch_dir);
    pnl.add(sc_depth);
    pnl.add(ch_view);
    add(pnl,"North");
  }

  public void start() {
    if (thmcs == null) {
      thmcs = new Thread(mcs);
      thmcs.start();
    }
    if (th == null) {
      th = new Thread(this);
      th.start();
    }
  }

  public void stop() {
    if (th != null) th = null;
    if (thmcs != null) thmcs = null;
  }


  // ---------------------------------- event listener -----------

  public void itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_dir) {
      mcs.setcdir( 1 - ch_dir.getSelectedIndex() );
    } else if (ev.getSource() == ch_view) {
      dispMode = ch_view.getSelectedIndex();
    }
  }

  public void actionPerformed(ActionEvent ev){
    if (ev.getSource() == bt_reset){ 
      mcs.reset();
      thCount = 0;
    } else if (ev.getSource() == bt_startStop){
      if (mcs.getStartSW()==0) {
        mcs.setStartSW(1);
      } else {
        mcs.setStartSW(0);
      }
    } else if (ev.getSource() == bt_step){
      if (mcs.getStartSW()==0) {
        mcs.setStepSW();
      }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    double vs;

    if (ev.getSource() == sc_depth) {
      viewDepth =  0.01*(double)(sc_depth.getValue());
    } else if (ev.getSource() == sc_gate) {
      mcs.setGatePotential( 0.01*(double)(sc_gate.getValue()) - 2.0);
    }
  }

  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 (dispMode==0 || dispMode==1 || dispMode==2) {
      viewTheta += 0.5*3.14159/180.0*(dgX-dgXb);
      viewFai += 0.5*3.14159/180.0*(dgY-dgYb);
      cosTh = Math.cos(viewTheta); 
      sinTh = Math.sin(viewTheta);
      cosFi = Math.cos(viewFai);
      sinFi = Math.sin(viewFai);
    } else {
      xShift += (1.0/viewScale/dispScale)*(dgX-dgXb);
      yShift += (1.0/viewScale/dispScale)*(dgY-dgYb);
      xImageLoc = gx0+(int)(dispScale*viewScale*xShift);
      yImageLoc = gy0+(int)(dispScale*viewScale*yShift);
    }
  }



  // ========================= run() - paint() loop ==============

  public void run() {
    while (th != null) {
      try {
        Thread.sleep(sleepTime);
      } catch (InterruptedException e) { }
      thCount += 1;
      offPaint();
      repaint();
    }
  }

  public void update(Graphics g) {
    paint(g);
  }

  public void paint(Graphics g) {
    g.drawImage(imgOff,0,0,this);
  }


  // --------------------------------------- off-paint -----------

  private void offPaint() {

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

    if (dispMode==0) {
      walker3DPlot(60,90,1.0);

    } else if (dispMode==1) {
      walker3DPlot(60,90,viewDepth);

    } else if (dispMode==2) {
      gridPlot(60,90,viewDepth);

    } else if (dispMode==3) {
      ;

    } else if (dispMode==4) {
      ;

    }

    gOff.setColor(Color.blue);
    gOff.drawString("t="+mcs.getTime()+" ",10,40);
    gOff.setColor(Color.black);
    if (mcs.getStartSW()==0) {
      gOff.drawString("stopped",dispWidth/6*1+10,40);
    } else {
      gOff.drawString("started",dispWidth/6*1+10,40);
    }
    gOff.drawString("current",dispWidth/6*3+10,40);
    gOff.drawString("z="+(int)(viewDepth*100.0)+"%",dispWidth/6*4+10,40);
    gOff.drawString("view",dispWidth/6*5+10,40);

    gOff.setColor(Color.cyan);
    gOff.drawString("Nelec="+mcs.getNparticles(-1)+" ",dispWidth/6*0+10,60);
    gOff.setColor(Color.magenta);
    gOff.drawString("Nhole="+mcs.getNparticles(1)+" ",dispWidth/6*1+10,60);
    gOff.setColor(Color.black);
    gOff.drawString("I="+(int)(mcs.getMeanCurrent()*100.0)/100.0+"",dispWidth/6*3+10,60);

    gOff.setColor(Color.black);
    gOff.drawString("count MD="+mcs.getCount()+" ",400,440);
    gOff.drawString("disp="+thCount+" ",530,440);

    gOff.drawString("Vnn="+(float)mcs.potentialOfN()+" ",30,440);
    gOff.drawString("Vpp="+(float)mcs.potentialOfP()+" ",230,440);
  }


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

  private void walker3DPlot(int gbx, int gby, double depth) {
    int i, j, k, knd, ix, iy, ir, NNx, NNy,NNz, nz;
    float col;
    double sc,x,y,z,cxx,cyy,czz, ppx,ppy,ppz;



    NNx = mcs.getNNx();
    NNy = mcs.getNNy();
    NNz = mcs.getNNz();
	nz = (int)(depth*NNz);
    sc = (480.0/NNx);
	cxx = 0.5*sc*NNx; cyy = 0.5*sc*NNy; czz = 0.5*sc*NNz;

    ir = (int)sc; if (ir<1) ir = 1;

    for (i=0; i<NNx; i++) {
	  x = sc*i;
      for (j=0; j<NNy; j++) {
	    y = sc*j;
	    for (k=0; k<nz; k++) {
		  z = sc*k;
          knd = mcs.getField(i,j,k);
          if (knd==0 || knd==3) continue;
		  
          ppy = cosFi*(y-cyy)+sinFi*(z-czz) + cyy; 
          ppz = -sinFi*(y-cyy)+cosFi*(z-czz) + czz;
          ppx = cosTh*(x-cxx)+sinTh*(ppz-czz) + cxx; 
          ppz = -sinTh*(x-cxx)+cosTh*(ppz-czz) + czz;

          ix = (int)(ppx)+gbx;
          iy = (int)(ppy)+gby;
          if (knd==1) {
            gOff.setColor(Color.magenta);
          } else if (knd==-1) {
            gOff.setColor(Color.cyan);
          } else if (knd==2) {
            gOff.setColor(Color.red);
          } else if (knd==-2) {
            gOff.setColor(Color.blue);
          }
          gOff.fillRect(ix,iy, ir, ir);
		}
      }
    }
  }

  //

  private void gridPlot(int gbx, int gby, double depth) {
    int cx,cy,cz,Nsx,Nsy,Nsz;
    int i,j,k,ic;
    double sc,mag,z,px,py,pz,f;

    Nsx = mcs.getNNx()/4;
    Nsy = mcs.getNNy()/4;
    Nsz = mcs.getNNz()/4;
	k = (int)(depth*(Nsz-1));
    cx = Nsx/2;
    cy = Nsy/2;
    cz = Nsz/2;
    sc = (480.0/Nsx);
    ic = 1;
    mag = 10.0;

    for(j=0; j<Nsy; j+=ic) {
      for (i=0; i<Nsx; i++) {
        f = mcs.getPotential(i,j,k);
        z = mag*f+k;
        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;
        fpts[i] = f;
      }
      for (i=0; i<Nsx-1; i++) {
        setLineColor(fpts[i],fpts[i+1]);
        gOff.drawLine(xpts[i],ypts[i],xpts[i+ic],ypts[i+ic]);
      }
    }

    for(i=0; i<Nsx; i+=ic) {
      for (j=0; j<Nsy; j++) { 
        f = mcs.getPotential(i,j,k);
        z = mag*f+k;
        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;
        fpts[j] = f;
      }
      for (j=0; j<Nsy-1; j++) {
        setLineColor(fpts[j],fpts[j+1]);
        gOff.drawLine(xpts[j],ypts[j],xpts[j+ic],ypts[j+ic]);
      }
    }

  }

  private void setLineColor(double f1, double f2) {
    float col,s,d;
    double f;

    s = 0.8f;
    if (f1*f2<0.0) {
      gOff.setColor(Color.getHSBColor(0.33f,s,0.1f));
      return;
    }
    f = 0.5*(f1+f2);
    col = 0.50f; if (f<0.0) col = 0.05f;
    if (f>0.0) {
      d = (float)(0.5*f+0.1); if (d>=0.99f) d=0.99f;
    } else {
      f = -f;
      d = (float)(0.5*f+0.1); if (d>=0.99f) d=0.99f;
    }
    gOff.setColor(Color.getHSBColor(col,s,d));
  }

}



// =================== Monte-Carlo simulation 3D class ===========

class MCS3D implements Runnable {

  private int t = 0; // system time
  private int dt = 1; // time division
  // private int hh = 1; // space division v0=hh/dt

  private int NBlock = 4;
  private int Nsx = 48;
  private int Nsy = 32;
  private int Nsz = 32;
  private int NNx = Nsx*NBlock;
  private int NNy = Nsy*NBlock;
  private int NNz = Nsy*NBlock;

  private double alpha = 1.0; // vd/v0 = a = alpha E ; drift velocity
  private int cdir = 0; // current direction { -1, 0, 1 }
  private double pExtinction = 0.5; // electron-hole collision-extinction probability
  private double gatePotential = 0.0;
  private double iMean = 0.0; // mean current (through current source)
  private double issMean = 0.0;
  private double iddMean = 0.0;

  // point difinition
  private int p_x1 = 4*8;
  private int p_x2 = p_x1 + 4*16;
  private int p_x3 = p_x2 + 4*16;
  private int p_y1 = 4*8;
  private int p_y2 = p_y1 + 4*16;
  private int p_z1 = 4*8;
  private int p_z2 = p_z1 + 4*16;

  // region difinition .... xPos, yPos, zPos, xSize, ySize, zSize
  private int region_world[] = { 0, 0, 0, NNx, NNy, NNz };
  private int region_n[] = { p_x1, p_y1, p_z1, p_x2-p_x1, p_y2-p_y1, p_z2-p_z1 }; // n-type
  private int region_p[] = { p_x2, p_y1, p_z1, p_x3-p_x2, p_y2-p_y1, p_z2-p_z1 }; // p-type
  private int region_nTerm[] = { p_x1, p_y1+12, p_z1+8, 12, p_y2-p_y1-16, p_z2-p_z1-16 }; // n terminal
  private int region_pTerm[] = { p_x3-12, p_y1+8, p_z1+8, 12, p_y2-p_y1-16, p_z2-p_z1-16 }; // p terminal

  private int NNp = 100000; // number of particles
  private int NN = 100000;
  private byte kind[] = new byte[NN]; // particle kind -1-electron, 0-dead, 1-hole
  private int px[] = new int[NN]; // particle x-pos
  private int py[] = new int[NN]; // particle y-pos
  private int pz[] = new int[NN]; // particle y-pos

  private byte field[][][] = new byte[NNx][NNy][NNz]; // -2-acceptor,-1-electron, 0-space, 1-hole, 2-donor, 3-wall
  private float potential[][][] = new float[Nsx][Nsy][Nsz]; // potential
  private float charge[][][] = new float[Nsx][Nsy][Nsz]; // charge

  private int sleepTime = 1;
  private int stopSleepTime = 100;
  private int resetSW = 0;
  private int startSW = 1;
  private int stepSW = 0;
  private int statMode = 2;
  private int count = 0;


  // ------------------------------- class constructor -----------

  MCS3D() {
    setInitialCondition();
  }


  // -------------------------------------- thread run -----------

  public void run() {
    while (true) {
      count += 1;
      timeEvolution();
      if (sleepTime>0) {
        try {
          Thread.sleep(sleepTime);
        } catch (InterruptedException e) { }
      }
    }
  }


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

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

    t = 0;
    count = 0;

    fillField(3, region_world); //clear (outer space = 3)

    i = 0;
    i = setNtype(i, region_n); // n-type
    i = setPtype(i, region_p); // p-type
    NNp = i;

    initPotential();
  }

  private int setNtype(int ip, int region[]) { // n-type semiconductor
    int ii;

    fillField(0, region); // inner space = 0

    ii = ip;
    ii = setBlock(ii,-1, 1, 1, 1, region); // electron
    ii = setBlock(ii, 2, 2, 2, 2, region); // donor
    return ii;
  }

  private int setPtype(int ip, int region[]) { // p-type semiconductor
    int ii;

    fillField(0, region); // inner space = 0

    ii = ip;
    ii = setBlock(ii, 1, 1, 2, 1, region); // hole
    ii = setBlock(ii,-2, 2, 1, 2, region); // acceptor
    return ii;
  }

  private int setItype(int ip, int region[]) { // intrinsic semiconductor

    fillField(0, region); // inner space = 0

    return ip;
  }

  private void fillField(int f, int region[]) {
    int ix,iy,iz, xPos,yPos,zPos, xSize,ySize,zSize;

    xPos = region[0]; yPos = region[1]; zPos = region[2];
	xSize = region[3]; ySize = region[4]; zSize = region[5];

    for (ix=xPos; ix<xPos+xSize; ix++) {
      for (iy=yPos; iy<yPos+ySize; iy++) {
	    for (iz=zPos; iz<zPos+zSize; iz++) {
          field[ix][iy][iz] = (byte)f;
		}
      }
    }
  }

  private int setBlock(int ii, int knd, int sx, int sy, int sz, int region[]) {
    int i,j,k, nx,ny,nz, ip;
    int ix,iy,iz, xPos,yPos,zPos, xSize,ySize,zSize;

    xPos = region[0]; yPos = region[1]; zPos = region[2];
	xSize = region[3]; ySize = region[4]; zSize = region[5];

    ip = ii;
    nx = xSize/4;
    ny = ySize/4;
    nz = zSize/4;

    for (i=0; i<nx; i++) {
      for (j=0; j<ny; j++) {
	    for (k=0; k<nz; k++) {
          ix = i*4+sx+xPos;
          iy = j*4+sy+yPos;
		  iz = k*4+sz+zPos;
          field[ix][iy][iz] = (byte)knd;
          if (knd==-1 || knd==1) {
            kind[ip] = (byte)knd;
            px[ip] = ix;
            py[ip] = iy;
		    pz[ip] = iz;
            ip += 1;
		  }
        }
      }
    }
    return ip;
  }

  private void initPotential() {
    int i,j,k;

    for (i=0; i<Nsx; i++) {
      for (j=0; j<Nsy; j++) {
	    for (k=0; k<Nsz; k++) {
          potential[i][j][k] = 0.0f;
          charge[i][j][k] = 0.0f;
		}
      }
    }
    setCharge();
    poisson(10);
  }


  // ---------------------------------- time evolution -----------

  private void timeEvolution() {

    if (resetSW==1) {
      setInitialCondition();
      resetSW = 0;
    }
	
    if (startSW==1) {
      timeStep();

    } else {
      if (stepSW==1) {
        timeStep();
        stepSW = 0;
      }
      try {
        Thread.sleep(stopSleepTime);
      } catch (InterruptedException e) { }
    }
  }

  private void timeStep() {
    int i,k, ir,ix,iy,iz; 
    double r, a, q;

    t = t + dt;

    setCharge();
    poisson(1);

    for (i=0; i<NNp; i++) {
      k = kind[i];
      if (k==0) continue;
      r = Math.random()*6.0;
      if (r<2.0) { // x-direction
        a = alpha*k*exField(px[i],py[i],pz[i]);
        if (r>=1+a) {
          ix = px[i]+1;
        } else {
          ix = px[i]-1;
        }
        iy = py[i];
		iz = pz[i];
      } else if (r<4.0) { // y-direction
        r = r - 2.0;
        ix = px[i];
        a = alpha*k*eyField(px[i],py[i],pz[i]);
        if (r>=1+a) {
          iy = py[i]+1;
        } else {
          iy = py[i]-1;
        }
		iz = pz[i];
      } else { // z-direction
        r = r - 4.0;
        ix = px[i];
		iy = py[i];
        a = alpha*k*ezField(px[i],py[i],pz[i]);
        if (r>=1+a) {
          iz = pz[i]+1;
        } else {
          iz = pz[i]-1;
        }	
      }
      move(i,ix,iy,iz);
    }

    q = 0;
    if (cdir>=0 || potentialOf(region_nTerm)<1.0) {
      q = currentSource(0.5*cdir, region_pTerm, region_nTerm);
    }
    iMean = 0.001*q + 0.999*iMean;
  }

  private int move(int i, int ix, int iy, int iz) {
    int ib,jb,kb,irem;

    if (ix<0 || ix>=NNx || iy<0 || iy>=NNy || iz<0 || iz>=NNz) return 0;
    if (field[ix][iy][iz]!=0) {
      ib = px[i]; jb = py[i]; kb = pz[i];
      if (field[ix][iy][iz]*field[ib][jb][kb]!=-1) return 0;
      // electron - hole extinction
      if (Math.random()>pExtinction) return 0;
      removeParticle(i);
      irem = searchParticle(ix,iy,iz);
      if (irem>=0) removeParticle(irem);

      return 0;
    }
    ib = px[i]; jb = py[i]; kb = pz[i];
    field[ib][jb][kb] = 0;
    field[ix][iy][iz] = kind[i];
    px[i] = ix;
    py[i] = iy;
	pz[i] = iz;
    return -1;
  }

  // set charge

  private void setCharge() {
    int i,j,k, ii,jj,kk, knd,cg;

    for (i=1; i<Nsx-1; i++) {
      for (j=1; j<Nsy-1; j++) {
	    for (k=1; k<Nsz-1; k++) {
          cg = 0;
          for (ii = i*NBlock;ii<i*NBlock+4; ii++) {
            for (jj = j*NBlock;jj<j*NBlock+4; jj++) {
			  for (kk = k*NBlock;kk<k*NBlock+4; kk++) {
                knd = field[ii][jj][kk];
                if (knd==1 || knd==2) { // hole or donor
                  cg += 1;
                } else if (knd==-1 || knd==-2) { // electron or acceptor
                  cg -= 1;
				}
              }
            }
          }
          charge[i][j][k] = cg;
		}
      }
    }
  }

  // Poisson

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

    for (iter=1; iter<=iterMax; iter++) {
      for (i=1; i<Nsx-1; i++) {
        for(j=1; j<Nsy-1; j++) {
		  for(k=1; k<Nsz-1; k++) {
		    potential[i][j][k] += 0.20f*( potential[i+1][j][k]+potential[i-1][j][k]
			  +potential[i][j+1][k]+potential[i][j-1][k]+potential[i][j][k+1]+potential[i][j][k-1]
			  -6.0f*potential[i][j][k]+0.1f*charge[i][j][k] );
		  }
        }
      }
    }
  }

  private float exField(int ix, int iy, int iz) {
    int iix,iiy,iiz;

    iix = ix/NBlock; iiy = iy/NBlock; iiz = iz/NBlock;
    if (iix==0 || iix==Nsx-1) return 0.0f;
    return 0.5f*(potential[iix+1][iiy][iiz]-potential[iix-1][iiy][iiz]);
  }

  private float eyField(int ix, int iy, int iz) {
    int iix,iiy,iiz;

    iix = ix/NBlock; iiy = iy/NBlock; iiz = iz/NBlock;
    if (iiy==0 || iiy==Nsy-1) return 0.0f;
    return 0.5f*(potential[iix][iiy+1][iiz]-potential[iix][iiy-1][iiz]);
  }

  private float ezField(int ix, int iy, int iz) {
    int iix,iiy,iiz;

    iix = ix/NBlock; iiy = iy/NBlock; iiz = iz/NBlock;
    if (iiz==0 || iiz==Nsz-1) return 0.0f;
    return 0.5f*(potential[iix][iiy][iiz+1]-potential[iix][iiy][iiz-1]);
  }

  // set potential
 
  private int setPotentialN(double pot, int region[]) { // return charge changed
    int ip,q;
    double p;

    p = potentialOf(region);
    q = 0;
    if (p>1.02*pot) {
      ip = putParticle(-1, region);
      if (ip>=0) q = -1;
    } else if (p<0.98*pot) {
      ip = removableParticleIndex(-1, region);
      if (ip>=0) {
        removeParticle(ip);
        q = 1;
      }
    }
    return q;
  }
 
  private int setPotentialP(double pot, int region[]) { // return charge changed
    int ip,q;
    double p;

    p = potentialOf(region);
    q = 0;
    if (p>1.02*pot) {
      ip = removableParticleIndex( 1, region);
      if (ip>=0) {
        removeParticle(ip);
        q = -1;
      }
    } else if (p<0.98*pot) {
      ip = putParticle( 1, region);
      if (ip>=0) q = 1;
    }
    return q;
  }
 
  // current source

  private int currentSource(double pp, int source[], int sink[]) { // return charge changed
    int rn,rp,q;

    q = 0;

    if (pp>=0.0 && Math.random()<pp) {
      putParticle(-1, sink); // add electron into sink(n-type)
      putParticle( 1, source); // add hole into source(p-type)
      q = 1;
    } else if (pp<0.0 && Math.random()<-pp ) {
      rn = removableParticleIndex(-1, sink); // remove electron from sink(n-type)
      rp = removableParticleIndex( 1, source); // remove hole from source(p-type)
      if (rn>=0 && rp>=0) {
        removeParticle(rn);
        removeParticle(rp);
        q = -1;
      }
    }
	
    return q;
  }

  //

  private int putParticle(int knd, int region[]) { // return particle index
    int n,nxy,i,ir;
    int ix,iy,iz, xPos,yPos,zPos, xSize,ySize,zSize;

    xPos = region[0]; yPos = region[1]; zPos = region[2];
	xSize = region[3]; ySize = region[4]; zSize = region[5];

    n = xSize*ySize*zSize;
    nxy = xSize*ySize;
    if (n<=0) return -1;
    if (knd!=-1 && knd!=1) return -1;
    ir = (int)(n*Math.random());
    for (i=0; i<n; i++) {
      ix = (((ir+i)%n)%nxy)%xSize+xPos;
      iy = (((ir+i)%n)%nxy)/xSize+yPos;
	  iz = ((ir+i)%n)/nxy+zPos;
      if (field[ix][iy][iz]==0) {
        return creatParticle(knd, ix, iy, iz);
      }
    }
    return -1;
  }

  private int removableParticleIndex(int knd, int region[]) { // return particle index
    int n,nxy, i,ir;
    int ix,iy,iz, xPos,yPos,zPos, xSize,ySize,zSize;

    xPos = region[0]; yPos = region[1]; zPos = region[2];
	xSize = region[3]; ySize = region[4]; zSize = region[5];

    n = xSize*ySize*zSize;
    nxy = xSize*ySize;
    if (n<=0) return -1;
    if (knd!=-1 && knd!=1) return -1;
    ir = (int)(n*Math.random());
    for (i=0; i<n; i++) {
      ix = (((ir+i)%n)%nxy)%xSize+xPos;
      iy = (((ir+i)%n)%nxy)/xSize+yPos;
	  iz = ((ir+i)%n)/nxy+zPos;
      if (field[ix][iy][iz]==knd) {
        return searchParticle(ix,iy,iz);
      }
    }
    return -1;
  }

  // creat/remove particle

  private int creatParticle(int knd, int ix, int iy, int iz) { // return particle index
    int idx;

    if (ix<0 || ix>=NNx || iy<0 || iy>=NNy || iz<0 || iz>=NNz) return -1;
    if (knd!=1 && knd!=-1) return -1;
    if (field[ix][iy][iz]!=0) return -1;
    idx = searchIndex(); if (idx==-1) return -1;
	
    field[ix][iy][iz] = (byte)knd;
    kind[idx] = (byte)knd;
    px[idx] = ix;
    py[idx] = iy;
    pz[idx] = iz;
    return idx;
  }

  private int searchIndex() {
    int i;

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

  //

  private void removeParticle(int i) {
    int ii,jj,kk;

    ii = px[i]; jj = py[i]; kk = pz[i];
    field[ii][jj][kk] = 0;
    kind[i] = 0;
  }

  //

  private int searchParticle(int ix, int iy, int iz) {
    int i;

    for (i=0; i<NNp; i++) {
      if (kind[i]!=0 && px[i]==ix && py[i]==iy && pz[i]==iz) return i;
    }
    return -1;
  }

  //

  private double potentialOf(int region[]) {
    int xPos,yPos,zPos, xSize,ySize,zSize;

    xPos = region[0]; yPos = region[1]; zPos = region[2];
	xSize = region[3]; ySize = region[4]; zSize = region[5];
    return potential[(xPos+xSize/2)/4][(yPos+ySize/2)/4][(zPos+zSize/2)/4];
  }

  public double potentialOfN() {
    return potentialOf(region_nTerm);
  }

  public double potentialOfP() {
    return potentialOf(region_pTerm);
  }

  // ----------------------------------------- utility -----------

  public int getNparticles(int knd) {
    int i,n;

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

  }


  // ------------------------------------- I/O methods -----------

  // control
  public void reset() { resetSW = 1; count = 0; }
  public void setStartSW(int sw) { startSW = sw; }
  public int getStartSW() { return startSW; }
  public void setStepSW() { stepSW = 1; }
  public int getCount() { return count; }

  public void setcdir(int d) { cdir = d; }
  public int getcdir() { return cdir; }
  public void setGatePotential(double gv) { gatePotential = gv; }
  public double getGatePotential() { return gatePotential; }

  // get data
  public int getNNp() { return NNp; }
  public int getNNx() { return NNx; }
  public int getNNy() { return NNy; }
  public int getNNz() { return NNz; }

  public int getTime() { return t; }
  public int getKind(int i) { return kind[i]; }
  public int getPx(int i) { return px[i]; }
  public int getPy(int i) { return py[i]; }
  public int getField(int ix, int iy, int iz) { return field[ix][iy][iz]; }
  public double getPotential(int ix, int iy, int iz) { return potential[ix][iy][iz]; }
  public double getMeanCurrent() { return 100.0*iMean; }
}