/** applet No. 1237
 *
 * grain growth in poly-crystal tungsten (NTP)
 * - periodic molecular dynamics 2D
 * - selected Morse potential
 * - multi thread  SMD :: display - asynchronous
 *
 * Created by Ikeuchi Mitsuru on November 18 2007.
 * Copyright (c) 2007 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2007.11.18  created
 *
 */

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

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

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

  Thread th = null; // for run()-paint() thread
  Thread thmd = null; // PeriodicMD2D thread
  Particles ma = new Particles(0); // W
  PeriodicMD2D md = new PeriodicMD2D(ma);
  DrawGraph2D dg2d = new DrawGraph2D();

  // for event
  Choice ch_mat1, ch_mat2, ch_tempControl, ch_view;
  Button bt_reset, bt_startStop, bt_step, bt_restart;
  Scrollbar sc_content, sc_temp, sc_press, sc_xLen, sc_scale;

  TextArea ta;
  String CRLF = System.getProperty("line.separator");

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

  private int dgX, dgY, dgXb, dgYb; // mouse

  private int printCount = 0;
  private int sleepTime = 50;
  private double viewScale = 1.0;
  private int dispMode = 3;
  private int thCount = 0;


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

  public void init() {
    resize(630,460);
    setBackground(Color.white);

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

    ch_mat1 = new Choice();
    for (int i=0; i<ma.getKindMax(); i++) {
      ch_mat1.add(ma.getMatStr(i));
    }
    ch_mat1.addItemListener(this);
    ch_mat1.select(ma.getMatStr(4));

    ch_mat2 = new Choice();
    for (int i=0; i<ma.getKindMax(); i++) {
      ch_mat2.add(ma.getMatStr(i));
    }
    ch_mat2.addItemListener(this);
    ch_mat2.select(ma.getMatStr(0));

    ch_tempControl = new Choice();
    ch_tempControl.add("adiabatic");
    ch_tempControl.add("t scaling"); 
    //ch_tempControl.add("wall cont.");
    ch_tempControl.addItemListener(this);
    ch_tempControl.select("t scaling"); 

    ch_view = new Choice();
    ch_view.add("ball");
    ch_view.add("line");
    ch_view.add("ball+line");
    ch_view.add("grain");
    ch_view.add("ball+grain");
    ch_view.add("defect");
    ch_view.add("temp");
    ch_view.add("velocity");
    ch_view.add("vSpace");
    ch_view.addItemListener(this);
    ch_view.select("grain");

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

    bt_restart = new Button("restart");
    bt_restart.addActionListener(this);

    sc_content= new Scrollbar(Scrollbar.HORIZONTAL,50,10,0,110);
    sc_content.addAdjustmentListener(this);

    sc_temp= new Scrollbar(Scrollbar.HORIZONTAL,100,10,1,210);
    sc_temp.addAdjustmentListener(this);

    sc_press= new Scrollbar(Scrollbar.HORIZONTAL,200,10,0,410);
    sc_press.addAdjustmentListener(this);

    sc_xLen= new Scrollbar(Scrollbar.HORIZONTAL,300,10,280,320);
    sc_xLen.addAdjustmentListener(this);

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

    addMouseListener(this);  
    addMouseMotionListener(this);

    setLayout(null);
    ta = new TextArea(" ",5,20);
    add(ta);
    ta.setBounds(400,120,220,100);// posX, posY, sizeX, sizeY
    ta.setText("time(ps), defect, grain size"+CRLF);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(2,6,5,0));
    //pnl.add(new Label("  "));
    pnl.add(bt_reset);
    pnl.add(bt_startStop);
	pnl.add(bt_step);
    pnl.add(ch_tempControl);
    pnl.add(new Label("  "));
    pnl.add(ch_view);

    pnl.add(bt_restart);
    pnl.add(new Label("  "));
    pnl.add(new Label("  "));
    pnl.add(sc_temp);
    pnl.add(sc_press);
    pnl.add(sc_scale);
    add(pnl,"North");

  }

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

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


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

  public void itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_mat1) {
      md.setMAT1( ch_mat1.getSelectedIndex() );
    } else if (ev.getSource() == ch_mat2) {
      md.setMAT2( ch_mat2.getSelectedIndex() );
    } else if (ev.getSource() == ch_tempControl) {
      md.setTempMode( ch_tempControl.getSelectedIndex() );
    } else if (ev.getSource() == ch_view) {
      dispMode = ch_view.getSelectedIndex();
    }
  }

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

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_content) { 
      md.setMat1content( sc_content.getValue() );
    } else if (ev.getSource() == sc_temp) { 
      md.setContTemp( 10.0*(double)(sc_temp.getValue()) );
    } else if (ev.getSource() == sc_press) { 
      md.setContPress( 0.01*(double)(sc_press.getValue())-2.0 );
    } else if (ev.getSource() == sc_xLen) { 
      md.setxLength( 1.0e-10*(double)(sc_xLen.getValue()) );
    } else if (ev.getSource() == sc_scale) { 
      viewScale = 0.01*(double)(sc_scale.getValue());
      dg2d.setScale(viewScale);
    }
  }

  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 (dgX<400) dg2d.setViewDelta( (dgX-dgXb), (dgY-dgYb) );
  }



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

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

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

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

  private void offPaint() {
    if (md != null) {
      dg2d.plotGraph2D(gOff, md, dispMode, thCount);
    }
  }

  private void printTA() {
    double tp;

    tp = md.getTime()*1.0e12;
    if (tp>=printCount+1.0) {
      ta.append(""+(printCount+1)+", "+md.getNumberOfDefect()+", "+
	    (int)(md.getMeanGrainSize()*1.0e9*10.0+0.5)/10.0+CRLF);
      printCount = printCount+1;
    }
  }
}



// =============================== draw graph 2D class ===========

class DrawGraph2D {

  private Graphics gOff;
  private PeriodicMD2D md;

  private int Nmt;

  private double xMax;
  private double yMax;

  private double plotSize = 300.0;
  private double dispScale = plotSize/md.getxMax();
  private double viewScale = 1.0;

  private int gx0 = 30;
  private int gy0 = 110;
  private int xImageLoc = gx0;
  private int yImageLoc = gy0;

  private int dragSW = 1;
  private double xShift = 0.0;
  private double yShift = 0.0;
  private double viewTheta = -15.0*3.14/180.0;
  private double viewFai = -72.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 rdf[][] = new int[120][3];
  private double meanRDF[][] = new double [120][3];
  private int vdf[][] = new int [200][3];
  private double meanVDF[][] = new double [200][3];
  private int adf[] = new int [220];

  private int drawnCount = 0;


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

  DrawGraph2D() {
  }


  // ----------------------------------- class methods -----------

  public void plotGraph2D(Graphics gg, PeriodicMD2D md2d, int dispMode, int thCount) {
    int gbh,mdCount;

    gOff = gg;
    md = md2d;

    if (md2d==null) return;
    mdCount = md.getCount();
    if (mdCount==drawnCount) return;
    drawnCount = mdCount;

    Nmt = md.getNmt();
    xMax = md.getxMax(); yMax = md.getyMax();

    gOff.setColor(Color.white);
    gOff.fillRect(0,0,630,460);

    if (dispMode==0) { // ball
      ballPlot(0.6, 0.0);
      drawBox();
      rdfPlot(400, 220, 0.4);

    } else if (dispMode==1) { // line
      linePlot();
      drawBox();
      rdfPlot(400, 220, 0.4);

    } else if (dispMode==2) { // ball + line
      ballPlot(0.6, 0.0);
      linePlot();
      drawBox();
      rdfPlot(400, 220, 0.4);

    } else if (dispMode==3) { // grain
      drawBox();
      grainPlot();
      adfPlot(400, 220, 0.4);

    } else if (dispMode==4) { // ball + grain
      drawBox();
      ballPlot(0.6, 0.0);
      grainPlot();
	  adfPlot(400, 220, 0.4);

    } else if (dispMode==5) { // defect
      drawBox();
      nearPlot();


    } else if (dispMode==6) { // temp
      ballPlot(0.6, 2000.0);
      drawBox();
      rdfPlot(400, 220, 0.4);

    } else if (dispMode==7) { // velocity
      velocityPlot(500.0);
      drawBox();
      vdfPlot(400, 220, 0.01);

    } else if (dispMode==8) { // v-space
      vSpacePlot(320, 3000.0);
      vdfPlot(400, 220, 0.01);

    } else if (dispMode==9) {
      ;

    }

    gbh = 630/6;

    gOff.setColor(Color.black);
    gOff.drawString("poly-crystal W",gbh*0+10,70);
    gOff.drawString("Tcont="+md.getContTemp()+"K",gbh*3+10,70);
    gOff.drawString("Pcont="+(int)(md.getContPress()*100.0+0.5)/100.0+"",gbh*4+10,70);
    gOff.drawString("scale="+(int)(viewScale*100.0+0.499)+"%",gbh*5+10,70);

    gOff.setColor(Color.blue);
    gOff.drawString("box="+(int)(md.getxMax()*1.0e12+0.5)/1000.0+"x" 
	  +(int)(md.getyMax()*1.0e12+0.5)/1000.0+"nm",140,70);

    gOff.setColor(Color.blue);
    gOff.drawString("t="+(int)(md.getTime()*1.0e14+0.5)/100.0+" ps",gbh*0+10,90);
    gOff.drawString("N="+md.getNmt()+" ",gbh*2+10,90);
    gOff.drawString("T="+(int)(md.getTemperature()*10+0.5)/10.0+" K",gbh*3+10,90);
    gOff.drawString("P="+(int)(md.getPressure()*1000.0)/1000.0+"N/m",gbh*4+10,90);

    //gOff.drawString("Pxx="+(int)(md.getxxPress()*1000.0)/1000.0+"N/m",400,110);
    //gOff.drawString("Pyy="+(int)(md.getyyPress()*1000.0)/1000.0+"N/m",400,125);
    //gOff.drawString("Pxy="+(int)(md.getxyPress()*1000.0)/1000.0+"N/m",400,140);

    //gOff.setColor(Color.black);
    //gOff.drawString("thread MD = "+md.getCount()+" ",20,440);
    //gOff.drawString("thread disp = "+thCount+" ",220,440);

    gOff.setColor(Color.black);
    gOff.drawString("time(ps), defect, grain size(nm)",400,115);
  }


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

  private void ballPlot(double dia, double tempSpan) {
    int i,ix,iy,ir;
    double sc,tm;

    sc = viewScale*dispScale;
    for (i=0; i<Nmt; i++) {
      ix = xImageLoc + (int)(sc*md.getxx(i));
      iy = yImageLoc + (int)(sc*md.getyy(i));
      ir = (int)(sc*dia*md.getDia(i));
      tm = 0.95;
      if (tempSpan>0.0) {
        tm = md.getTemp(i);
        tm = 0.8*tm/tempSpan+0.199;
        if (tm>0.99) tm = 0.99;
      }
      gOff.setColor(Color.getHSBColor((float)(md.getColor(i)),0.95f,(float)tm));
      gOff.fillOval(ix-ir/2,iy-ir/2, ir, ir);
    }
  }

  private void grainPlot() {
    int i,j,k,ix,iy,ix2,iy2;
    double sc,x,y,r2,dij,gcol;

    for (i=0; i<200; i++) {
      adf[i] = 0;
    }

    sc = viewScale*dispScale;
    for (i=0; i<Nmt; i++) {
      ix = xImageLoc + (int)(sc*md.getxx(i));
      iy = yImageLoc + (int)(sc*md.getyy(i));

      for(k=1;k<md.getnnReg(i);k++) {
        j = md.getReg(i, k);
        x = md.getxx(i)-md.getxx(j); y = md.getyy(i)-md.getyy(j);
        r2 = x*x+y*y;
        dij = 1.1*(md.getDia(i)+md.getDia(j))/2.0;
        if (r2<dij*dij && j>i) {
          ix2 = xImageLoc + (int)(sc*md.getxx(j));
          iy2 = yImageLoc + (int)(sc*md.getyy(j));
          gcol = grainColor(x,y);
          adf[(int)(gcol*200.0)] += 1;
          gOff.setColor(Color.getHSBColor((float)gcol,0.9f,0.9f));
          gOff.drawLine(ix,iy, ix2,iy2);
        }
      }
    }
  }

  private double grainColor(double x, double y) {
    double th,col;

    th = 3.0*(Math.atan2(y,x)/Math.PI + 1.0);
    return (th-(int)th);
  }

  private void linePlot() {
    int i,j,k,ix,iy,ix2,iy2;
    double sc,r2,dij;

    sc = viewScale*dispScale;
    for (i=0; i<Nmt; i++) {
      ix = xImageLoc + (int)(sc*md.getxx(i));
      iy = yImageLoc + (int)(sc*md.getyy(i));

      for(k=1;k<md.getnnReg(i);k++) {
        j = md.getReg(i, k);
        r2 = (md.getxx(i)-md.getxx(j))*(md.getxx(i)-md.getxx(j))+(md.getyy(i)-md.getyy(j))*(md.getyy(i)-md.getyy(j));
        dij = 1.1*(md.getDia(i)+md.getDia(j))/2.0;
        if (r2<dij*dij && j>i) {
          ix2 = xImageLoc + (int)(sc*md.getxx(j));
          iy2 = yImageLoc + (int)(sc*md.getyy(j));

          gOff.setColor(Color.gray);
          gOff.drawLine(ix,iy, ix2,iy2);
        }
      }
    }
  }

  private void nearPlot() {
    int i,j,k,ix,iy,ix2,iy2, ir, n;
    double sc,x,y,r2,dij,gcol;

    sc = viewScale*dispScale;
    for (i=0; i<Nmt; i++) {
      ix = xImageLoc + (int)(sc*md.getxx(i));
      iy = yImageLoc + (int)(sc*md.getyy(i));

      ir = 0;
      if (md.getMeanNear(i)<5.5) {
        ir = 2;
      }
      gOff.setColor(Color.getHSBColor(0.66f,0.9f,0.8f));
      if (ir>0) gOff.fillOval(ix-ir/2,iy-ir/2, ir, ir);
    }

    gOff.setColor(Color.red);
    n = md.getNumberOfDefect();
    gOff.drawString("no. of defect="+n+" ",400,400);
    gOff.drawString("mean grain size="+(int)(md.getMeanGrainSize()*1.0e9*10.0+0.5)/10.0+" nm",400,420);
    //if (1.0*n>0.4*Nmt) initMeanNear();
  }

  //

  private void velocityPlot(double vSpan) {
    int i,ix,iy,ix2,iy2;
    double sc;

    sc = viewScale*dispScale;
    for (i=0; i<Nmt; i++) {
      ix = xImageLoc + (int)(sc*md.getxx(i));
      iy = yImageLoc + (int)(sc*md.getyy(i));

      ix2 = ix+(int)(viewScale*10.0*md.getvx(i)/vSpan);
      iy2 = iy+(int)(viewScale*10.0*md.getvy(i)/vSpan);
      gOff.setColor(Color.getHSBColor((float)(md.getColor(i)),0.8f,0.8f));
      gOff.drawLine(ix,iy, ix2,iy2);
    }
  }

  //

  private void drawBox() {
    int ix,iy,gball;
    double sc;

    sc = viewScale*dispScale;
    gball = (int)(2.0*viewScale);

    gOff.setColor(Color.gray); 
    ix = (int)(sc*xMax);
    iy = (int)(sc*yMax);
    gOff.drawRect(xImageLoc-gball,yImageLoc-gball,ix+2*gball,iy+2*gball);
  }

  //

  private void vSpacePlot(int plotWidth, double scale) {
    int i,ix,iy,ir,ix0,iy0;

    ix0 = xImageLoc + plotWidth/2;
    iy0 = yImageLoc + plotWidth/2;
    gOff.setColor(Color.black);
    gOff.drawLine(ix0, yImageLoc, ix0, yImageLoc+plotWidth);
    gOff.drawLine(xImageLoc, iy0, xImageLoc+plotWidth, iy0);

    for (i=0; i<Nmt; i++) {
      ix = (int)(ix0 + 0.5*plotWidth*md.getvx(i)/scale );
      iy = (int)(iy0 + 0.5*plotWidth*md.getvy(i)/scale );
      ir = (int)(0.6*md.getDia(i)/1.0e-10+0.5);

      gOff.setColor(Color.getHSBColor(md.getColor(i),0.4f,0.9f));
      gOff.fillOval(ix-ir/2,iy-ir/2, ir, ir);

    }
  }

  // redial distribution

  private void rdfPlot(int gx, int gy, double mag) {
    int i,n1,n2,hoganWidth,hoganHight,xDiv,yDiv;
    double a;

    setRadialDistribution();

    hoganWidth = 200; hoganHight = 200;
    xDiv = hoganWidth/5; yDiv = hoganHight/2;
    a = mag*(25000.0/Nmt);

    gOff.setColor(Color.lightGray);
    for (i=10; i<100; i++) {
      n1 = (int)(a*meanRDF[i][0]/i);
      if (n1>0) {
        gOff.fillRect(gx+i*2,gy+hoganHight-n1,2,n1);
      }
    }

    gOff.setColor(Color.getHSBColor(md.getColorOf(md.getMAT1()),0.9f,0.6f));
    for (i=10; i<100-1; i++) {
      n1 = (int)(a*meanRDF[i][1]/i);
      n2 = (int)(a*meanRDF[i+1][1]/i);
      gOff.drawLine(gx+i*2,gy+hoganHight-n1,gx+i*2+2,gy+hoganHight-n2);
    }

    gOff.setColor(Color.getHSBColor(md.getColorOf(md.getMAT2()),0.9f,0.6f));
    for (i=10; i<100-1; i++) {
      n1 = (int)(a*meanRDF[i][2]/i);
      n2 = (int)(a*meanRDF[i+1][2]/i);
      gOff.drawLine(gx+i*2,gy+hoganHight-n1,gx+i*2+2,gy+hoganHight-n2);
    }

    drawHogan(gx, gy, hoganWidth, hoganHight, xDiv, yDiv);
    gOff.setColor(Color.black);
    gOff.drawString("0",gx-5,gy+hoganHight+20);
    gOff.drawString("1 (nm)",gx+hoganWidth-20,gy+hoganHight+20);
    gOff.drawString("radial distribution (n/r)",gx+20,gy+hoganHight+40);
  }

  private void setRadialDistribution() {
    int i,j,k,ir,m1,m2,ki,kj;
    double xij,yij,rij;

    for (i=0; i<120; i++) {
      rdf[i][0]=0;
      rdf[i][1]=0;
      rdf[i][2]=0;
    }

    m1 = md.getMAT1(); m2 = md.getMAT2();
    for (i=0; i<Nmt; i++) {
      for(k=1;k<md.getnnReg(i);k++) {
        j = md.getReg(i, k);
        xij = md.getxx(i)-md.getxx(j); yij = md.getyy(i)-md.getyy(j);
        rij = Math.sqrt(xij*xij+yij*yij);
        if (rij<10.0e-10) {
          ki = md.getKind(i); kj = md.getKind(j);
          ir = (int)(rij*1.0e11+0.5); if (ir>119) ir=119;
          rdf[ir][0] +=1;
          if (ki==m1 && kj==m1) {
            rdf[ir][1] += 1;
          } else if (ki==m2 && kj==m2) {
            rdf[ir][2] += 1;
          }
        }
      }
    }

    for (i=0; i<120; i++) {
      meanRDF[i][0] = 0.9*meanRDF[i][0]+0.1*rdf[i][0];
      meanRDF[i][1] = 0.9*meanRDF[i][1]+0.1*rdf[i][1];
      meanRDF[i][2] = 0.9*meanRDF[i][2]+0.1*rdf[i][2];
    }
  }

  // velocity distribution

  private void vdfPlot(int gx, int gy, double mag) {
    int i,n1,n2;
    double yMag;

    setVelocityDistribution();

    yMag = mag*0.005*Nmt;
    gOff.setColor(Color.lightGray);
    for (i=0; i<200-1; i++) {
      n1  = (int)(yMag*meanVDF[i][0]);
      gOff.drawLine(gx+i,gy+200-n1,gx+i,gy+200);
    }

    gOff.setColor(Color.getHSBColor(md.getColorOf(md.getMAT1()),0.9f,0.6f));
    for (i=0; i<200-2; i++) {
      n1 = (int)(yMag*meanVDF[i][1]);
      n2 = (int)(yMag*meanVDF[i+1][1]);
      gOff.drawLine(gx+i,gy+200-n1,gx+i+1,gy+200-n2);
    }

    gOff.setColor(Color.getHSBColor(md.getColorOf(md.getMAT2()),0.9f,0.6f));
    for (i=0; i<200-2; i++) {
      n1 = (int)(yMag*meanVDF[i][2]);
      n2 = (int)(yMag*meanVDF[i+1][2]);
      gOff.drawLine(gx+i,gy+200-n1,gx+i+1,gy+200-n2);
    }

    drawHogan(gx, gy, 200, 200, 50, 100);
    gOff.setColor(Color.black);
    gOff.drawString("0",gx-5,gy+220);
    gOff.drawString("1000",gx+100-15,gy+220);
    gOff.drawString("(m/s)",gx+200-30,gy+220);
    gOff.drawString("velocity distribution",gx+20,gy+235);
  }

  private void setVelocityDistribution() {
    int i,iv,m1,m2,knd;
    double vx,vy,v;

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

    m1 = md.getMAT1(); m2 = md.getMAT2();
    for(i=0;i<Nmt;i++) {
      knd = md.getKind(i);
      vx = md.getvx(i); vy = md.getvy(i);
      v = Math.sqrt(vx*vx+vy*vy);
      iv = (int)(v/10.0); if (iv>199) iv = 199;
      vdf[iv][0] += 1;
      if (knd==m1) {
        vdf[iv][1] += 1;
      } else if (knd==m2) {
        vdf[iv][2] += 1;
      }
    }

    for(i=0;i<200;i++) {
      meanVDF[i][0] = 0.98*meanVDF[i][0]+0.02*vdf[i][0];
      meanVDF[i][1] = 0.98*meanVDF[i][1]+0.02*vdf[i][1];
      meanVDF[i][2] = 0.98*meanVDF[i][2]+0.02*vdf[i][2];
    }
  }

  // azimuth distribution

  private void adfPlot(int gx, int gy, double mag) {
    int i,n1,n2;
    double yMag;

    yMag = 10000.0/Nmt*mag;
    gOff.setColor(Color.lightGray);
    for (i=0; i<200; i++) {
      n1  = (int)(yMag*adf[i]);
      gOff.setColor(Color.getHSBColor((float)(i/200.0),0.9f,0.8f));
      gOff.drawLine(gx+i,gy+200-n1,gx+i,gy+200);
    }

    drawHogan(gx, gy, 200, 200, 50, 100);
    gOff.setColor(Color.black);
    gOff.drawString("0",gx-5,gy+220);
    gOff.drawString("PI/3",gx+200-15,gy+220);
    gOff.drawString("azimuth distribution",gx+20,gy+235);
  }

  //

  private void drawHogan(int gx, int gy, int hoganWidth, int hoganHight, int xDiv, int yDiv) {
    int i;

    gOff.setColor(Color.gray);
    for (i=0; i<hoganWidth; i+=xDiv) {
      gOff.drawLine(gx+i,gy, gx+i, gy+hoganHight);
    }
    for (i=0; i<hoganHight; i+=yDiv) {
      gOff.drawLine(gx,gy+hoganHight-i, gx+hoganWidth, gy+hoganHight-i);
    }
    gOff.setColor(Color.black);
    gOff.drawRect(gx,gy,hoganWidth,hoganHight);
  }


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

  public void setViewDelta(double xDisplacement, double yDisplacement) {

    if (dragSW==1) {
      xShift += (1.0/viewScale/dispScale)*xDisplacement;
      yShift += (1.0/viewScale/dispScale)*yDisplacement;
      xImageLoc = gx0+(int)(dispScale*viewScale*xShift);
      yImageLoc = gy0+(int)(dispScale*viewScale*yShift);
    } else if (dragSW==2) {
      viewTheta += 0.5*3.14/180.0*xDisplacement;
      viewFai += 0.5*3.14/180.0*yDisplacement;
      cosTh = Math.cos(viewTheta); 
      sinTh = Math.sin(viewTheta);
      cosFi = Math.cos(viewFai);
      sinFi = Math.sin(viewFai);
    }
  }
  public void setDragSW(int sw) { dragSW = sw; }
  public void setScale(double sc) { viewScale = sc; }

}



// =================================== Particles class ===========
//
// selected Morse atom
//

class Particles {

  private static double AU = 1.661e-27;
  private static double EE = 1.602e-19;
  private static double AA = 1.0e-10;

  private static double Morse[][] = {
  // 0:mass     1:charge  2:D(eV)   3:A(A^-1)  4:r0(A)   5:color 
    {183.85 *AU, 0.0*EE, 0.9906*EE, 1.4116/AA, 3.032*AA, 0.30 }, // 0 W
    { 95.94 *AU, 0.0*EE, 0.8032*EE, 1.5079/AA, 2.976*AA, 0.33 }, // 1 Mo
    { 51.996*AU, 0.0*EE, 0.4414*EE, 1.5721/AA, 2.754*AA, 0.36 }, // 2 Cr
    { 55.847*AU, 0.0*EE, 0.4174*EE, 1.3885/AA, 2.845*AA, 0.39 }, // 3 Fe
    { 58.71 *AU, 0.0*EE, 0.4205*EE, 1.4199/AA, 2.780*AA, 0.42 }, // 4 Ni
    { 26.98 *AU, 0.0*EE, 0.2703*EE, 1.1646/AA, 3.253*AA, 0.45 }, // 5 Al
    {207.19 *AU, 0.0*EE, 0.2348*EE, 1.1836/AA, 3.733*AA, 0.48 }, // 6 Pb
    { 63.54 *AU, 0.0*EE, 0.3429*EE, 1.3588/AA, 2.866*AA, 0.20 }, // 7 Cu
    {107.87 *AU, 0.0*EE, 0.3323*EE, 1.3690/AA, 3.115*AA, 0.24 }, // 8 Ag
    { 40.08 *AU, 0.0*EE, 0.1623*EE, 0.8054/AA, 4.569*AA, 0.13 }, // 9 Ca
    { 87.62 *AU, 0.0*EE, 0.1513*EE, 0.7878/AA, 4.988*AA, 0.16 }, // 10 Sr
    {137.34 *AU, 0.0*EE, 0.1416*EE, 0.6570/AA, 5.373*AA, 0.19 }, // 11 Ba
    { 22.99 *AU, 0.0*EE, 0.0633*EE, 0.5900/AA, 5.336*AA, 0.01 }, // 12 Na
    { 39.102*AU, 0.0*EE, 0.0542*EE, 0.4977/AA, 6.369*AA, 0.04 }, // 13 K
    { 85.47 *AU, 0.0*EE, 0.0464*EE, 0.4298/AA, 7.207*AA, 0.07 }, // 14 Rb
    {132.905*AU, 0.0*EE, 0.0449*EE, 0.4157/AA, 7.557*AA, 0.10 }, // 15 Cs
    { 20.183*AU, 0.0*EE, 0.0031*EE, 1.6500/AA, 3.076*AA, 0.56 }, // 16 Ne
    { 39.948*AU, 0.0*EE, 0.0104*EE, 1.3400/AA, 3.816*AA, 0.60 }, // 17 Ar
    { 83.80 *AU, 0.0*EE, 0.0141*EE, 1.2500/AA, 4.097*AA, 0.64 }, // 18 Kr
    {131.30 *AU, 0.0*EE, 0.0200*EE, 1.2400/AA, 4.467*AA, 0.68 }, // 19 Xe
    {200.59 *AU, 0.0*EE, 0.0734*EE, 1.4900/AA, 3.255*AA, 0.72 }  // 20 Hg
  };

  private static String matStr[] = {
    "W" ,"Mo","Cr","Fe","Ni","Al","Pb","Cu","Ag","Ca",
    "Sr","Ba","Na","K" ,"Rb","Cs","Ne","Ar","Kr","Xe","Hg"
  };
  private static int kindMax = 21;

  private static double hh = 1.0e-12;
  private static double rCutoff = 1.0e-9;
  private double forceTable[][][] = new double[kindMax][kindMax][];

  private static double dt = 5.0*1.0e-15;


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

  Particles() {
    int i,j;

    for (i=0; i<kindMax; i++) {
      for (j=0; j<kindMax; j++) {
        forceTable[i][j] = new double[1020];
      }
    }
    for (i=0; i<kindMax; i++) {
      for (j=0; j<kindMax; j++) {
        setForceTable(i, j);
      }
    }
  }

  Particles(int knd) {
    forceTable[knd][knd] = new double[1020];
    setForceTable(knd,knd);
  }

  Particles(int ki, int kj) {
    forceTable[ki][ki] = new double[1020];
    forceTable[kj][kj] = new double[1020];
    forceTable[ki][kj] = new double[1020];
    forceTable[kj][ki] = new double[1020];
    setForceTable(ki,ki);
    setForceTable(kj,kj);
    setForceTable(ki,kj);
    setForceTable(kj,ki);
  }

  //

  private void setForceTable(int ki, int kj) {
    int ir; 
    double d,a,r0,y,r,sc;

    d = Math.sqrt(Morse[ki][2]*Morse[kj][2]);
    a = 0.5*(Morse[ki][3]+Morse[kj][3]);
    r0 = 0.5*(Morse[ki][4]+Morse[kj][4]);
    for (ir=10; ir<1000; ir++) {
      r = (double)ir*hh;
      y = Math.exp(-a*(r-r0));
      forceTable[ki][kj][ir] = cutoff(r)*2.0*d*a*y*(y-1.0);
    }
    for (ir=0; ir<10; ir++) {
      forceTable[ki][kj][ir]=forceTable[ki][kj][10];
    }
  }

  private double cutoff(double r) {
    double rt;

    rt = 0.0;
    if (r<0.8*rCutoff) {
      rt = 1.0;
    } else if (r>=0.8*rCutoff && r<rCutoff) {
      rt = 0.5+0.5*Math.cos(3.1415926536*(r-0.8*rCutoff)/(0.2*rCutoff));
    } else {
      rt = 0.0;
    }
    return( rt );
  }


  // ----------------------------------- class methods -----------

  public double force(double r, int ki, int kj) {
    int i;
    double a;

    i = (int)(r/hh);
    a = r - i*hh;
    return ( ((hh-a)*forceTable[ki][kj][i] + a*forceTable[ki][kj][i+1])/hh );
  }


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

  public static int getKindMax() { return kindMax; }
  public static String getMatStr(int knd) { return matStr[knd]; }

  public static double getMassOf(int knd) { return Morse[knd][0]; }
  public static double getDiaOf(int knd) { return Morse[knd][4]; }
  public static double getColorOf(int knd) { return Morse[knd][5]; }
  public static double getrCutoff() { return rCutoff; }

  public static double getdt() { return dt; }
}



// ==================== periodic molecular dynamics 2D ===========

class PeriodicMD2D implements Runnable {

  private Particles particle;

  private double t = 0.0;
  private double dt = Particles.getdt();
  private static double xMax = 90.0*1.0e-9;
  private static double yMax = 90.0*1.0e-9;

  private int Nmt = 20000;

  private int Nsx = 180;
  private int Nsy = 180;

  private int MAT1 = 0; // W
  private int MAT2 = 0; // W
  private int mat2content = 50; // 0-100

  private int tempMode = 1;
  private double contTemp = 1000.0;
  private double contPress = 0.0;
  private double xVirial = 0.0;
  private double yVirial = 0.0;
  private double xyVirial = 0.0;
  private double xLength = 0.97*xMax;

  public final int NN = 100000;
  private int kind[] = new int[NN];
  private double xx[] = new double[NN];
  private double yy[] = new double[NN];
  private double vx[] = new double[NN];
  private double vy[] = new double[NN];
  private double ffx[] = new double[NN];
  private double ffy[] = new double[NN];

  private int reg[][] = new int[NN][100];
  private int section[][][] = new int[Nsx][Nsy][50];

  private int near[] = new int[NN];
  private float meanNear[] = new float[NN];

  private int nBlock = 18;
  private int azim[][] = new int[nBlock][nBlock];

  private int resetSW = 0;
  private int startSW = 1;
  private int stepSW = 1;
  private int restartSW = 0;
  private int scaleSW = 0;  
  private int sleepTime = 1;
  private int stopSleepTime = 100;
  private int count = 0;


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

  PeriodicMD2D(Particles p) {
    particle = p;
	setAzim();
    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; 
    double nm;

    count = 0;
    t = 0.0;

    nm = 1.0e-9;
    xMax = 90.0*nm; yMax = 90.0*nm;
    i = 0;

    for (ix=0; ix<nBlock; ix++) {
      for (iy=0; iy<nBlock; iy++) {
        i = setNanoBlock(i,MAT1, (0.1+5.0*ix)*nm, (0.1+5.0*iy)*nm,
             4.8*nm, 4.8*nm, azim[ix][iy]*Math.PI/15.0 );
      }
    }
    Nmt = i;

    scaleLength(xLength);
    initVirial();
    vAjustment();
    initMeanNear();
  }

  private void initVirial() {
    xVirial = 0.0; yVirial = 0.0; xyVirial = 0.0;
  }

  private void setAzim() {
    int ix,iy; 

    for (ix=0; ix<nBlock; ix++) {
      for (iy=0; iy<nBlock; iy++) {
        azim[ix][iy] = (int)(5.0*Math.random());
      }
    }
  }

  private int setNanoBlock(int ii, int knd, double x0, double y0, double xLen, double yLen, double th) {
    int ip,i,j,nx,ny;
    double a,b, len, x,y,xp,yp, sth,cth;

    ip = ii;
    a = 0.98*particle.getDiaOf(knd);
    b = 0.866025*a;
    len = xLen; if (len<yLen) len = yLen;
    len = 1.5*len;
    nx = (int)(len/b) + 1;
    ny = (int)(len/a) + 1;
    sth = Math.sin(th);
    cth = Math.cos(th);

    for (i=0; i<nx; i++) {
      x = b*i - len/2.0;
      for (j=0; j<ny; j++) {
        y = a*j - len/2.0; if (i%2==1) y += 0.5*a;
        xp = x0 + xLen/2.0 + cth*x - sth*y;
        yp = y0 + yLen/2.0 + sth*x + cth*y;
        if (xp>=x0 && xp<=x0+xLen && yp>=y0 && yp<=y0+yLen) {
          kind[ip] = knd;
          xx[ip] = xp;
          yy[ip] = yp;
          vx[ip] = 500.0*(Math.random()-0.5);
          vy[ip] = 500.0*(Math.random()-0.5);
          ffy[ip] = 0.0;
          ffy[ip] = 0.0;
          ip += 1;
        }
      }
    }
    return ip;
  }

  //

  private void initMeanNear() {
    int i;

    for (i=0; i<Nmt; i++) {
      meanNear[i] = 5.6f;
    }
  }


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

  private void timeEvolution(){

    if (resetSW==1) {
      setAzim();
      setInitialCondition();
      scaleLength(xLength);
      resetSW = 0;
    }

    if (restartSW==1) {
      setInitialCondition();
      scaleLength(xLength);
      restartSW = 0;
    }

    if (scaleSW==1) {
      //scaleLength(xLength);
      scaleSW = 0;
    }

    if (startSW==1) {
      timeStep();

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

  private void timeStep() {
    int i;

    if (tempMode==1) {
      vAjustment();
    }
    registration();
    for (i=0; i<10; i++) {
      time1Step();
      volControl();
    }
  }

  private void time1Step(){
    int i; 
    double dtv2, a2;

    dtv2 = dt/2.0;
    t = t + dt;

    for (i=0; i<Nmt; i++) {
      a2 = dtv2/particle.getMassOf(kind[i]);
      vx[i] += a2*ffx[i];
      vy[i] += a2*ffy[i];
      xx[i] += vx[i]*dt;
      yy[i] += vy[i]*dt;
    }
    forceCalc();
    for (i=0; i<Nmt; i++) {
      a2 = dtv2/particle.getMassOf(kind[i]);
      vx[i] += a2*ffx[i];
      vy[i] += a2*ffy[i];
    }

    setBoundary();
  }

  private void forceCalc() {
    int i,j,k;
    double xij,yij;
    double rc2,r2,rij,f,fxij,fyij,xvir,yvir,xyvir;

    rc2 = particle.getrCutoff()* particle.getrCutoff();

    for(i=0;i<Nmt;i++) {
      ffx[i] = 0.0;
      ffy[i] = 0.0;
    }
    xvir = 0.0; yvir = 0.0; xyvir = 0.0;

    for(i=0;i<Nmt;i++) {
      for(k=1;k<reg[i][0];k++) {
        j = reg[i][k];
        xij = xx[i] - xx[j]; yij = yy[i] - yy[j];
        if (xij>0.5*xMax) xij -= xMax;
        if (xij<-0.5*xMax) xij += xMax;  
        if (yij>0.5*yMax) yij -= yMax;
        if (yij<-0.5*yMax) yij += yMax; 
        r2=xij*xij+yij*yij;
        if (r2<rc2) {
          rij = Math.sqrt(r2);
          f = particle.force(rij,kind[i],kind[j]);
          fxij = f*xij/rij;
          fyij = f*yij/rij;

          ffx[i] += fxij;
          ffy[i] += fyij;
          ffx[j] -= fxij;
          ffy[j] -= fyij;

          xvir += fxij*xij;
          yvir += fyij*yij;
          xyvir +=  fxij*yij;
        }
      }
    }
    xVirial = 0.99*xVirial + 0.01*xvir;
    yVirial = 0.99*yVirial + 0.01*yvir;
    xyVirial = 0.99*xyVirial + 0.01*xyvir;
  }


  private void setBoundary() {
    int i; 
    double rr, tmp;

    rr = 1.0;
    if (tempMode==2) {
      tmp = temperature();
      rr = Math.sqrt(contTemp/tmp);
      if (rr<0.5) { rr = 0.5; } else if (rr>1.2) { rr=1.2; };
    }

    for (i=0; i<Nmt; i++) {
      if (xx[i] < 0.0) { 
        xx[i] += xMax; 
      }
      if (xx[i] > xMax) {
        xx[i] -= xMax; 
      }
      if (yy[i] < 0.0) { 
        yy[i] += yMax; 
      }
      if (yy[i] > yMax) {
        yy[i] -= yMax;
      }
    }
  }

  //registration and set meanNear

  private void registration() {
    int i,j,ii,jj,ip,jp,kp,i0,i1,j0,j1,iq;
    double xij,yij,r2,rreg,rreg2,dij;

    preRegistration();

    rreg = particle.getrCutoff() + 3000.0*10*dt;
    rreg2 = rreg*rreg;

    for (i=0; i<Nmt; i++) {
      near[i] = 0;
    }

    for(ip=0;ip<Nmt;ip++) {
      kp = 1;

      i0 = (int)(Nsx*(xx[ip]-rreg)/xMax+(double)Nsx)-Nsx;
      i1 = (int)(Nsx*(xx[ip]+rreg)/xMax);
      j0 = (int)(Nsy*(yy[ip]-rreg)/yMax+(double)Nsy)-Nsy;
      j1 = (int)(Nsy*(yy[ip]+rreg)/yMax);

      for(i=i0;i<=i1;i++) {
        ii = (i+Nsx)%Nsx;
        for(j=j0;j<=j1;j++) {
          jj = (j+Nsy)%Nsy;
          for(iq=1;iq<=section[ii][jj][0];iq++) {
            jp = section[ii][jj][iq];
            if (jp>ip) {
              xij = xx[ip] - xx[jp]; yij = yy[ip] - yy[jp];
              if (xij>0.5*xMax) xij -= xMax;
              if (xij<-0.5*xMax) xij += xMax;  
              if (yij>0.5*yMax) yij -= yMax;
              if (yij<-0.5*yMax) yij += yMax; 
              r2=xij*xij+yij*yij;
              if (r2<rreg2) {
                reg[ip][kp]=jp;
                kp += 1;
              }
              dij = 1.1*(particle.getDiaOf(kind[ip])+particle.getDiaOf(kind[jp]))/2.0;
              if (r2<dij*dij) {
                near[ip] +=1;
                near[jp] +=1;
              }
            }
          }
        }
      }
      reg[ip][0]=kp;
    }

    for (i=0; i<Nmt; i++) {
      meanNear[i] = 0.9f*meanNear[i]+0.1f*near[i];
    }
  }

  private void preRegistration() {
    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<Nmt;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;
    }
  }


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

  private void vAjustment() {
    int i;
    double tmp, r;

    tmp = temperature();
    r = Math.sqrt(contTemp/tmp);
    for (i=0; i<Nmt; i++) {
      vx[i] = r*vx[i];
      vy[i] = r*vy[i];
    }
  }

  private double temperature() {
    int i;
    double s;

    s = 0.0; 
    for (i=0; i<Nmt; i++) {
      s += particle.getMassOf(kind[i])*(vx[i]*vx[i] + vy[i]*vy[i]);
    }
    return (0.5*s/((double)Nmt*1.38e-23));
  }

  private double temp(int i) {
    double ke;

    ke = 0.5* particle.getMassOf(kind[i])*(vx[i]*vx[i] + vy[i]*vy[i]);
    return (ke/1.38e-23);
  }

  //

  private double pressure() {
    double vol;

    vol = xMax*yMax;
    return (Nmt*1.38e-23*temperature()+(xVirial+yVirial)/2.0)/vol;
  }

  private double xxPress() {
    int i;
    double nkt,vol;

    nkt= 0.0;
    for (i=0; i<Nmt; i++) {
      nkt += particle.getMassOf(kind[i])*vx[i]*vx[i];
    }

    vol = xMax*yMax;
    return( (nkt+xVirial)/vol );
  }

  private double yyPress() {
    int i;
    double nkt,vol;

    nkt= 0.0;
    for (i=0; i<Nmt; i++) {
      nkt += particle.getMassOf(kind[i])*vy[i]*vy[i];
    }

    vol = xMax*yMax;
    return( (nkt+yVirial)/vol );
  }

  private double xyPress() {
    int i;
    double nkt,vol;

    nkt= 0.0;
    for (i=0; i<Nmt; i++) {
      nkt += particle.getMassOf(kind[i])*vx[i]*vy[i];
    }

    vol = xMax*yMax;
    return( (nkt+xyVirial)/vol );
  }

  // vol control

  private void volControl() {
    double d,a,x,y,z;

    d = dt*20.0;
    a = 1.0e-4;

    x = (1.0+a*(xxPress()-contPress))*xMax;
    if (x>xMax+d) {
      x = xMax + d;
    } else if (x<xMax-d) {
      x = xMax - d;
    }
    xMax = x;

    y = (1.0+a*(yyPress()-contPress))*yMax;
    if (y>yMax+d) {
      y = yMax + d;
    } else if (y<yMax-d) {
      y = yMax - d;
    }
    yMax = y;

  }

  //

  private void scaleLength(double xlen) {
    int i;
    double rr;

    rr = xlen/xMax;
    xMax = xlen;
    yMax = rr*yMax;

    for(i=0;i<Nmt;i++) {
      xx[i] = rr*xx[i];
      yy[i] = rr*yy[i];
    }
  }

  //

  private int numberOfDefect() {
    int i,n;

    n = 0;
    for(i=0;i<Nmt;i++) {
      if (meanNear[i]<5.5) n += 1;
    }
    return(n);
  }

  //

  private int regMax() {
    int i,m;

    m = 0;
    for(i=0;i<Nmt;i++) {
      if (reg[i][0]>m) m = reg[i][0];
    }
    return(m);
  }


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

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

  public static double getxMax() { return xMax; }
  public static double getyMax() { return yMax; }

  public void setTempMode(int tmode) { tempMode = tmode; }
  public int getTempMode() { return tempMode; }
  public void setContTemp(double tmp) { contTemp = tmp; vAjustment(); }
  public double getContTemp() { return contTemp; }
  public void setContPress(double p) { contPress = p; }
  public double getContPress() { return contPress; }
  public void setxLength(double len) { xLength = len; scaleSW = 1; }
  public double getxLength() { return xLength; }

  public void setMAT1(int m) { MAT1 = m; resetSW = 1; }
  public int getMAT1() { return MAT1; }
  public void setMAT2(int m) { MAT2 = m; resetSW = 1; }
  public int getMAT2() { return MAT2; }
  public void setMat1content(int m2c) { mat2content = m2c; resetSW = 1; }
  public int getMat2content() { return mat2content; }

  // get data
  public double getTime() { return t; }
  public double getdt() { return dt; }
  public int getNmt() { return Nmt; }
  public float getColorOf(int knd) { return (float)particle.getColorOf(knd); }

  public double getDia(int i) { return particle.getDiaOf(kind[i]); }
  public double getMass(int i) { return particle.getMassOf(kind[i]); }
  public float getColor(int i) { return (float) particle.getColorOf(kind[i]); }

  public int getKind(int i) { return kind[i]; }
  public double getxx(int i) { return xx[i]; }
  public double getyy(int i) { return yy[i]; }
  public double getvx(int i) { return vx[i]; }
  public double getvy(int i) { return vy[i]; }
  public double getTemp(int i) { return temp(i); }

  public double getTemperature() { return temperature(); }
  public double getPressure() { return pressure(); }
  public double getxxPress() { return xxPress(); }
  public double getyyPress() { return yyPress(); }
  public double getxyPress() { return xyPress(); }

  public int getnnReg(int i) { return reg[i][0]; }
  public int getReg(int i, int k) { return reg[i][k]; }
  public int getRegMax() { return regMax(); }

  public double getMeanNear(int i) { return meanNear[i]; }
  public int getNumberOfDefect() { return numberOfDefect(); }
  public double getMeanGrainSize() {
    return xMax/(numberOfDefect()/Math.sqrt(Nmt*1.0)/4.0); 
  }
}