/** applet No. 1251
 *
 * n-cannel MOS FET - Monte-Carlo simulation 2D
 * - multi thread  MCS2D :: display - asynchronous
 *
 * Created by Ikeuchi Mitsuru on January 19 2008.
 * Copyright (c) 2008 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2008.01.19  created
 *
 */

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

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

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

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

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

  // 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 = 1;
  private double boxDispWidth = 320.0;

  private double dispScale = 1.0;
  private double viewScale = 1.0;

  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("e+h");
    ch_view.add("e+h+pot");
    ch_view.add("grid pot");
    ch_view.addItemListener(this);
    ch_view.select("e+h+pot");

    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_scale = new Scrollbar(Scrollbar.HORIZONTAL,100,10,50,510);
    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(bt_step);
    pnl.add(sc_gate);
    pnl.add(sc_scale);
    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_scale) {
      vs = viewScale;
      viewScale =  0.01*(double)(sc_scale.getValue());
      xImageLoc = xImageLoc - (int)(0.5*mcs.getNNx()*(viewScale - vs));
      yImageLoc = yImageLoc - (int)(0.5*mcs.getNNy()*(viewScale - vs));
      xShift = (xImageLoc-gx0)/(dispScale*viewScale);
      yShift = (yImageLoc-gy0)/(dispScale*viewScale);
    } 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==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) {
      walkerPlot();
      linePotentialPlot();

    } else if (dispMode==1) {
      potentialPlot();
      walkerPlot();

    } else if (dispMode==2) {
      gridPlot(10.0);

    } else if (dispMode==3) {
      ;

    }

    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("Vgate="+(int)(mcs.getGatePotential()*100.0)/100.0+"",dispWidth/6*3+10,40);
    gOff.drawString("scale="+(int)(viewScale*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("n channel junction FET :  VSS = 0.0, VDD = 3.0",160,80);


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

  }


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

  private void walkerPlot() {
    int i, j, k, ix, iy, ir, NNx, NNy;
    float col;
    double sc;

    NNx = mcs.getNNx();
    NNy = mcs.getNNy();
    sc = 1.0*viewScale;
    ir = (int)sc; if (ir<1) ir = 1;

    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        k = mcs.getField(i,j);
        if (k==0 || k==3) continue;
        ix = (int)(sc*i)+xImageLoc;
        iy = (int)(sc*j)+yImageLoc;
        if (k==1) {
          gOff.setColor(Color.magenta);
        } else if (k==-1) {
          gOff.setColor(Color.cyan);
        } else if (k==2) {
          gOff.setColor(Color.red);
        } else if (k==-2) {
          gOff.setColor(Color.blue);
        }
        gOff.fillRect(ix,iy, ir, ir);
      }
    }
  }

  private void linePotentialPlot() {
    int i,j,Nsx,gby;

    Nsx = mcs.getNNx()/4;
    j = mcs.getNNy()/4/2;
    gby = (int)(viewScale*160.0)+yImageLoc;

    for(i=0; i<Nsx; i++) {
      xpts[i] = (int)(viewScale*i*4)+xImageLoc;
      ypts[i] = (int)(viewScale*(160.0-50.0*mcs.getPotential(i,j)))+yImageLoc;
      fpts[j] = mcs.getPotential(i,j);
    }
    gOff.setColor(Color.lightGray);
    gOff.drawLine(xImageLoc,gby,(int)(viewScale*Nsx*4-4)+xImageLoc,gby);
    for (i=0; i<Nsx-2; i++) {
      gOff.setColor(Color.black);
      gOff.drawLine(xpts[i],ypts[i],xpts[i+1],ypts[i+1]);
    }
  }

  private void potentialPlot() {
    int i,j,ix,iy,Nsx,Nsy;
    double p,col;

    Nsx = mcs.getNNx()/4;
    Nsy = mcs.getNNy()/4;

    for (i=0; i<Nsx; i++) {
      for (j=0; j<Nsy; j++) {
        p = mcs.getPotential(i,j);
        if (p>0.0) {
          col = 0.61;
        } else {
          col = 0.01; p = -p;
        }
        p = 0.5*p; if (p>1.0) p=1.0;
        gOff.setColor(Color.getHSBColor((float)col,0.5f,(float)p));
        ix = (int)(i*4*viewScale)+xImageLoc;
        iy = (int)(j*4*viewScale)+yImageLoc;
        gOff.fillRect(ix,iy, (int)(4*viewScale+1), (int)(4*viewScale+1));
      }
    }
  }

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

    Nsx = mcs.getNNx()/4;
    Nsy = mcs.getNNy()/4;
    cx = Nsx/2;
    cy = Nsy/2;
    cz = 0;
    sc = (480.0/Nsx);
    ic = 2;

    gbx = 60;
    gby = 90;

    for(j=0; j<Nsy; j+=ic) {
      for (i=0; i<Nsx; i++) {
        f = mcs.getPotential(i,j);
        z = mag*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;
        fpts[i] = f;
      }
      for (i=0; i<Nsx-2; 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);
        z = mag*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;
        fpts[j] = f;
      }
      for (j=0; j<Nsy-2; 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 2D class ===========

class MCS2D 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 = 120;
  private int Nsy = 80;
  private int NNx = Nsx*NBlock;
  private int NNy = 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.01; // 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 = 60;
  private int p_x2 = p_x1 + 120;
  private int p_x3 = p_x2 + 120;
  private int p_x4 = p_x3 + 120;
  private int p_y1 = 60;
  private int p_y2 = p_y1 + 76;
  private int p_y3 = p_y2 + 4;
  private int p_y4 = p_y3 + 120;

  // area difinition .... xPos, yPos, xSize, ySize
  private int area_world[] = { 0, 0, NNx, NNy };
  private int area_ns[] = { p_x1, p_y3, p_x2-p_x1, p_y4-p_y3 }; // source
  private int area_p[]  = { p_x2, p_y3, p_x3-p_x2, p_y4-p_y3 }; // p-type substrate
  private int area_nd[] = { p_x3, p_y3, p_x4-p_x3, p_y4-p_y3 }; // drain
  private int area_ng[] = { p_x2, p_y1, p_x3-p_x2, p_y2-p_y1 }; // gate
  private int area_nsTerm[] = { p_x1, p_y3+16, 8, p_y4-p_y3-32 }; // sorce terminal
  private int area_pTerm[] = { p_x2+16, p_y4-8, p_x3-p_x2-32, 8 }; // p-type terminal
  private int area_ndTerm[] = { p_x4-8, p_y3+16, 8, p_y4-p_y3-32 }; // drain terminal
  private int area_ngTerm[] = { p_x2+16, p_y1, p_x3-p_x2-32, 8 }; // gate terminal

  private int NNp = 20000; // number of particles
  private int NN = 20000;
  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 byte field[][] = new byte[NNx][NNy]; // -2-acceptor,-1-electron, 0-space, 1-hole, 2-donor, 3-wall
  private float potential[][] = new float[Nsx][Nsy]; // potential
  private float charge[][] = new float[Nsx][Nsy]; // 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 -----------

  MCS2D() {
    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, area_world); //clear (outer space = 3)

    i = 0;
    i = setNtype(i, area_ns); // source
    i = setPtype(i, area_p); // p-type substrate
    i = setNtype(i, area_nd); // drain
    i = setNtype(i, area_ng); // gate
    NNp = i;

    initPotential();
  }

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

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

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

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

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

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

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

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

    return ip;
  }

  private void fillField(int f, int area[]) {
    int ix,iy, xPos,yPos,xSize,ySize;

    xPos = area[0]; yPos = area[1]; xSize = area[2]; ySize = area[3];

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

  private int setBlock(int ii, int knd, int sx, int sy, int area[]) {
    int i,j, nx,ny,ip;
    int ix,iy, xPos,yPos,xSize,ySize;

    xPos = area[0]; yPos = area[1]; xSize = area[2]; ySize = area[3];

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

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

  private void initPotential() {
    int i,j;

    for (i=0; i<Nsx; i++) {
      for (j=0; j<Nsy; j++) {
        potential[i][j] = 0.0f;
        charge[i][j] = 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; 
    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()*4.0;
      if (r<2.0) { // x-direction
        a = alpha*k*exField(px[i],py[i]);
        if (r>=1+a) {
          ix = px[i]+1;
        } else {
          ix = px[i]-1;
        }
        iy = py[i];
      } else { // y-direction
        r = r - 2.0;
        ix = px[i];
        a = alpha*k*eyField(px[i],py[i]);
        if (r>=1+a) {
          iy = py[i]+1;
        } else {
          iy = py[i]-1;
        }	
      }
      move(i,ix,iy);
    }

    setPotentialN(gatePotential, area_ngTerm);
    setPotentialP(-1.0, area_pTerm);
    q = setPotentialN( 0.0, area_nsTerm);
	issMean = 0.001*q + 0.999*issMean;
    q = setPotentialN( 3.0, area_ndTerm);
    iddMean = 0.001*q + 0.999*iddMean;
	iMean = 0.001*(iddMean - issMean)+0.999*iMean;
  }

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

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

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

  // set charge

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

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

  // Poisson

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

    for (iter=1; iter<=iterMax; iter++) {
      for (i=1; i<Nsx-1; i++) {
        for(j=1;j<Nsy-1;j++) {
          potential[i][j] += 0.33f*(potential[i+1][j]+potential[i-1][j]+potential[i][j+1]+potential[i][j-1]-4.0f*potential[i][j]+0.0625f*charge[i][j]);
        }
      }
    }
  }

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

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

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

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

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

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

    p = potentialOf(area);
    q = 0;
    if (p>1.02*pot) {
      ip = removableParticleIndex( 1, area);
      if (ip>=0) {
        removeParticle(ip);
        q = -1;
      }
    } else if (p<0.98*pot) {
      ip = putParticle( 1, area);
      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 area[]) { // return particle index
    int n,i,ix,iy,ir;
    int xPos,yPos,xSize,ySize;

    xPos = area[0]; yPos = area[1]; xSize = area[2]; ySize = area[3];

    n = 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)%xSize+xPos;
      iy = ((ir+i)%n)/xSize+yPos;
      if (field[ix][iy]==0) {
        return creatParticle(knd, ix, iy);
      }
    }
    return -1;
  }

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

    xPos = area[0]; yPos = area[1]; xSize = area[2]; ySize = area[3];

    n = 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)%xSize+xPos;
      iy = ((ir+i)%n)/xSize+yPos;
      if (field[ix][iy]==knd) {
        return searchParticle(ix,iy);
      }
    }
    return -1;
  }

  // creat/remove particle

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

    if (ix<0 || ix>=NNx || iy<0 || iy>=NNy) return -1;
    if (knd!=1 && knd!=-1) return -1;
    if (field[ix][iy]!=0) return -1;
    idx = searchIndex(); if (idx==-1) return -1;
	
    field[ix][iy] = (byte)knd;
    kind[idx] = (byte)knd;
    px[idx] = ix;
    py[idx] = iy;
    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;

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

  //

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

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

  //

  private double potentialOf(int area[]) {
    int xPos,yPos,xSize,ySize;

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

  /*
  public double gatePotential() {
   return potentialOf(area_p1);
  }

  public double sourcePotential() {
   return potentialOf(area_ns);
  }

  public double drainPotential() {
   return potentialOf(area_nd);
  }
  */


  // ----------------------------------------- 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 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) { return field[ix][iy]; }
  public double getPotential(int ix, int iy) { return potential[ix][iy]; }
  public double getMeanCurrent() { return 100.0*iMean; }
}