/* template - periodic smoothed particle hydrodynamics 3D */
/* coded by Ikeuchi Mitsuru   */
/* ver 0.0.1   2005.05.21     */
/* ver 0.0.2   2005.07.10 bug of dwwvdr(r,h) fixed */

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

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

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

  Button bt_reset, bt_start, bt_stop;
  Choice ch_view;
  Scrollbar sc_rate,sc_alpha;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 30;

  int dispMode = 0;
  int incRate = 10;
  int resetSW = 0;
  int started = 1;
  double alpha = 1.0;

  double xMax = 10.0;
  double yMax = 10.0;
  double zMax = 10.0;
  
  int dgX, dgY, dgXb,dgYb;
  double cx = 0.5*xMax;
  double cy = 0.5*yMax;
  double cz = 0.5*zMax;
  double theta = -20.0*3.14/180.0;
  double fai = 10.0*3.14/180.0;
  double dtheta = 0.0*3.14/180.0;
  double pai = 3.1415926536;

  double t = 0.0;
  double dt = 0.01;
  double hh0 = 0.25;
  double rBond = hh0*1.2;

  int NNp = 1000+100;
  int kind[] = new int [NNp];
  int bond[][] = new int[NNp][20];
  double xx[] = new double [NNp];
  double yy[] = new double [NNp];
  double zz[] = new double [NNp];
  double vx[] = new double [NNp];
  double vy[] = new double [NNp];
  double vz[] = new double [NNp];
  double ax[] = new double [NNp];
  double ay[] = new double [NNp];
  double az[] = new double [NNp];

  double hh[] = new double [NNp];
  double mass[] = new double [NNp];
  double density[] = new double [NNp];
  double pressure[] = new double [NNp];
  double energy[] = new double [NNp];
  double power[] = new double [NNp];

  int Nsx = 20;
  int Nsy = 20;
  int Nsz = 20;
  int section[][][][] = new int[Nsx][Nsy][Nsz][200];

  int srtz[][] = new int[2][NNp];

  double px[] = new double[NNp]; 
  double py[] = new double[NNp]; 
  double pz[] = new double[NNp]; 

  double wkx[] = new double[8]; 
  double wky[] = new double[8]; 
  double wkz[] = new double[8]; 

  double pwkx[] = new double[8]; 
  double pwky[] = new double[8]; 
  double pwkz[] = new double[8]; 

  int boxp[][] = {
   {0,1},{0,2},{0,4},{1,3},{1,5},{2,3},{2,6},{4,5},{4,6},{3,7},{5,7},{6,7}
  };

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

  public void init() {

    resize(630,320);
    setBackground(Color.white);

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

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

    bt_start= new Button("start");
    bt_start.addActionListener(this); 

    bt_stop = new Button("stop");
    bt_stop.addActionListener(this);

    ch_view = new Choice();
    ch_view.addItem("particle"); ch_view.addItem("buble"); 
    ch_view.addItem("x-colored"); ch_view.addItem("y-colored");
    ch_view.addItem("z-colored"); ch_view.addItem("density");
    ch_view.addItem("pressure"); ch_view.addItem("velocity");
    ch_view.addItemListener(this);
    ch_view.select("particle");

    sc_rate= new Scrollbar(Scrollbar.HORIZONTAL,10,10,1,110);
    sc_rate.addAdjustmentListener(this);

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

    addMouseListener(this);  
    addMouseMotionListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
    pnl.add(bt_reset);
    pnl.add(bt_start);
    pnl.add(bt_stop);
    pnl.add(new Label("  "));
    pnl.add(sc_alpha);
    pnl.add(ch_view);
    add(pnl,"North");

    setInitialCondition();
    setWaku();
  }

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

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

  public void itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_view) {
      dispMode = ch_view.getSelectedIndex();
    }
  }

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

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if(ev.getSource() == sc_rate) { 
      incRate = sc_rate.getValue();
    } else if(ev.getSource() == sc_alpha) { 
      alpha = 0.1*sc_alpha.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;
    }
    theta += 0.5*3.14/180.0*(dgX-dgXb);
    fai += 0.5*3.14/180.0*(dgY-dgYb);
  }


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

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

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

/* ----------------------------- offPaint -------------------- */

  void offPaint() {
    int i;
    int gbx,gby;
    double sc;

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

    rotate(theta,fai);
    zSort();

    gbx = 70;
    gby = 80;
    sc = 20.0;

    plotBoxInit(gbx, gby, sc);

    if (dispMode==0) {
	  plotParticle(gbx, gby, sc, 1);

	} else if (dispMode==1) {
      plotParticle(gbx, gby, sc, 0);

    } else if (dispMode==2) {
      plotColorParticle(gbx, gby, sc, 0);

	} else if (dispMode==3) {
      plotColorParticle(gbx, gby, sc, 1);

	} else if (dispMode==4) {
      plotColorParticle(gbx, gby, sc, 2);

	} else if (dispMode==5) {
      plotColorOf(density, 0.01, gbx, gby, sc);

	} else if (dispMode==6) {
      plotColorOf(pressure, 0.02, gbx, gby, sc);

	} else if (dispMode==7) {
      plotVelocity(100.0, gbx, gby, sc);

	} else if (dispMode==8) {

	}

    plotBoxFin(gbx, gby, sc);


    gOff.setColor(Color.black);
    gOff.drawString("t="+(int)(t*10.0+0.5)/10.0+" ",630/6*0+10,40);
	gOff.drawString("viscosity="+(int)(alpha*10.0)/10.0+" ",630/6*4+10,40);
    gOff.drawString("view",630/6*5+10,40);
/*
    gOff.drawString("maxSec="+maxSection()+" ",150,40);
*/
  }

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

  void plotParticle(int gbx, int gby, double scale, int filled) {
    int i,j,gx,gy,gz;
	double c;

    for (i=0; i<NNp; i++) {
        j = srtz[0][i];
        gx = (int)(scale*px[j])+gbx;
        gy = (int)(scale*py[j])+gby;
        gz = (int)(10.0*(0.4*(pz[j]/zMax)+0.6));
		c = 0.66;
		if (kind[j]!=0) c= 0.15;
        gOff.setColor(Color.getHSBColor((float)c,0.6f,(float)(0.4*pz[j]/zMax+0.4)));
		if (filled==1) {
          gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
		} else {
          gOff.drawOval(gx-gz/2,gy-gz/2, gz, gz);
		}
    }
  }

  void plotColorParticle(int gbx, int gby, double scale, int direction) {
    int i,j,gx,gy,gz;

    for (i=0; i<NNp; i++) {
      j = srtz[0][i];
      gx = (int)(scale*px[j])+gbx;
      gy = (int)(scale*py[j])+gby;
      gz = (int)(10.0*(0.4*(pz[j]/zMax)+0.6));
      if (direction==0) {
        gOff.setColor(Color.getHSBColor((float)(0.66*(1.0-xx[j]/xMax)),0.8f,(float)(0.4*pz[j]/zMax+0.4)));
      } else if (direction==1) {
        gOff.setColor(Color.getHSBColor((float)(0.66*yy[j]/yMax),0.8f,(float)(0.4*pz[j]/zMax+0.4)));
      } else if (direction==2) {
        gOff.setColor(Color.getHSBColor((float)(0.66*(1.0-zz[j]/zMax)),0.8f,(float)(0.4*pz[j]/zMax+0.4)));
      }
      gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
    }
  }

  void plotColorOf(double a[], double mag, int gbx, int gby, double scale) {
    int i,j,gx,gy,gz;
	double ai;

    for (i=0; i<NNp; i++) {
      j = srtz[0][i];
      gx = (int)(scale*px[j])+gbx;
      gy = (int)(scale*py[j])+gby;
      gz = (int)(10.0*(0.4*(pz[j]/zMax)+0.6));
	  ai = mag*a[j]; if (ai<=0.0) ai=0.0;
	  if (ai>1.0) ai=1.0;
      gOff.setColor(Color.getHSBColor((float)(0.66*(1.0-ai)),0.8f,(float)(0.4*pz[j]/zMax+0.4)));
      gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
    }
  }

  void plotVelocity(double mag, int gbx, int gby, double scale) {
    int i,j,gx,gy,gz,g2x,g2y,g2z;
	double xx2j,yy2j,zz2j,px2,py2,pz2;
    double cosTh,sinTh,cosFi,sinFi;

    cosTh = Math.cos(theta); sinTh = Math.sin(theta);
    cosFi = Math.cos(fai); sinFi = Math.sin(fai);

    for (i=0; i<NNp; i++) {
      j = srtz[0][i];
      gx = (int)(scale*px[j])+gbx;
      gy = (int)(scale*py[j])+gby;
      gz = (int)(10.0*(0.4*(pz[j]/zMax)+0.6));
	  
	  xx2j = xx[j]+vx[j]*dt*mag;
	  yy2j = yy[j]+vy[j]*dt*mag;
	  zz2j = zz[j]+vz[j]*dt*mag;
      py2 = cosFi*(yy2j-cy)+sinFi*(zz2j-cz) + cy; 
      pz2 = -sinFi*(yy2j-cy)+cosFi*(zz2j-cz) + cz;
      px2 = cosTh*(xx2j-cx)+sinTh*(pz2-cz) + cx; 
      pz2 = -sinTh*(xx2j-cx)+cosTh*(pz2-cz) + cz;
      g2x = (int)(scale*px2)+gbx;
      g2y = (int)(scale*py2)+gby;
      g2z = (int)(10.0*(0.4*(pz[j]/zMax)+0.6));

	  if (vy[j]<0.0) {
	    gOff.setColor(Color.getHSBColor(0.01f,0.8f,(float)(0.4*pz[j]/zMax+0.4)));
	  } else {
        gOff.setColor(Color.getHSBColor(0.66f,0.8f,(float)(0.4*pz[j]/zMax+0.4)));
	  }
      gOff.drawLine(gx,gy, g2x, g2y);
    }
  }


  void plotBoxInit(int gbx, int gby, double scale) {
    int i,gx,gy,g2x,g2y,farPoint;

    farPoint = findMin();
    for (i=0; i<12; i++) {
      gx = (int)(scale*pwkx[boxp[i][0]])+gbx;
      gy = (int)(scale*pwky[boxp[i][0]])+gby;
      g2x = (int)(scale*pwkx[boxp[i][1]])+gbx;
      g2y = (int)(scale*pwky[boxp[i][1]])+gby;
      if (boxp[i][0]==farPoint || boxp[i][1]==farPoint) { 
        gOff.setColor(Color.lightGray);
        gOff.drawLine(gx,gy, g2x, g2y);
      } 
    }
  }

  void plotBoxFin(int gbx, int gby, double scale) {
    int i,gx,gy,g2x,g2y,farPoint;

    farPoint = findMin();
    for (i=0; i<12; i++) {
      gx = (int)(scale*pwkx[boxp[i][0]])+gbx;
      gy = (int)(scale*pwky[boxp[i][0]])+gby;
      g2x = (int)(scale*pwkx[boxp[i][1]])+gbx;
      g2y = (int)(scale*pwky[boxp[i][1]])+gby;
      if (boxp[i][0]!=farPoint && boxp[i][1]!=farPoint) { 
        gOff.setColor(Color.gray);
        gOff.drawLine(gx,gy, g2x, g2y);
      }
    }
  }

  int findMin() {
    int i,im;
    double m;

    im = 0; m = pwkz[im];
    for (i=0; i<8; i++) {
      if (pwkz[i]<m) {
        m = pwkz[i];
        im = i;
      }
    }
    return(im);
  }

/*-------------------------------------  sort z  ----------*/

  void zSort() {
    int i,j,w;

    for (i=0; i<NNp; i++) {
      j = srtz[0][i];
      srtz[1][i] = (int)(32767.0*pz[j]/zMax);
    } 
    for (i=0; i<100; i++) {
      for (j=0; j<NNp-1; j++) {
        if (srtz[1][j]>srtz[1][j+1]) {
          w = srtz[0][j]; srtz[0][j]=srtz[0][j+1]; srtz[0][j+1]=w;
          w = srtz[1][j]; srtz[1][j]=srtz[1][j+1]; srtz[1][j+1]=w;
        }
      }
    }    
  }

  void zSortInit() {
    int i,j,w;

    for (i=0; i<NNp; i++) {
      srtz[0][i] = i;
      srtz[1][i] = (int)(32767.0*pz[i]/zMax);
    } 
    for (i=0; i<NNp; i++) {
      for (j=0; j<NNp-1; j++) {
        if (srtz[1][j]>srtz[1][j+1]) {
          w = srtz[0][j]; srtz[0][j]=srtz[0][j+1]; srtz[0][j+1]=w;
          w = srtz[1][j]; srtz[1][j]=srtz[1][j+1]; srtz[1][j+1]=w;
        }
      }
    }    
  }

/*-------------------------------------  rotate  ----------*/

  void rotate(double th,double fi) {
    int i;
    double cosTh,sinTh,cosFi,sinFi;

    cosTh = Math.cos(th); sinTh = Math.sin(th);
    cosFi = Math.cos(fi); sinFi = Math.sin(fi);
    for (i=0; i<NNp; i++) {
      py[i] = cosFi*(yy[i]-cy)+sinFi*(zz[i]-cz) + cy; 
      pz[i] = -sinFi*(yy[i]-cy)+cosFi*(zz[i]-cz) + cz;
    } 
    for (i=0; i<NNp; i++) {
      px[i] = cosTh*(xx[i]-cx)+sinTh*(pz[i]-cz) + cx; 
      pz[i] = -sinTh*(xx[i]-cx)+cosTh*(pz[i]-cz) + cz;
    } 

    for (i=0; i<8; i++) {
      pwky[i] = cosFi*(wky[i]-cy)+sinFi*(wkz[i]-cz) + cy; 
      pwkz[i] = -sinFi*(wky[i]-cy)+cosFi*(wkz[i]-cz) + cz;
    } 
    for (i=0; i<8; i++) {
      pwkx[i] = cosTh*(wkx[i]-cx)+sinTh*(pwkz[i]-cz) + cx; 
      pwkz[i] = -sinTh*(wkx[i]-cx)+cosTh*(pwkz[i]-cz) + cz;
    }  
  }

/*-------------------------------------  setWaku  ---------*/

  void setWaku() {
    int i,ix,iy,iz;
    double sgm;
	
    sgm = 0.25;
    i=0;

    for (ix=0; ix<2; ix++) {
      for (iy=0; iy<2; iy++) {
        for (iz=0; iz<2; iz++) {
          wkx[i] = -0.5*sgm+(xMax+sgm)*(double)(ix);
          wky[i] = -0.5*sgm+(yMax+sgm)*(double)(iy);
          wkz[i] = -0.5*sgm+(zMax+sgm)*(double)(iz);
          i += 1;
        }
      }
    }
  }


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

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

    t = 0.0;

    i = 0;
    for (ix=0; ix<10; ix++) {
      for (iy=0; iy<10; iy++) {
        for (iz=0; iz<10; iz++) {
          kind[i] = 0;
          mass[i] = 1.0;
          xx[i] = 0.2*(ix)+2.0;
          yy[i] = 0.2*(iy)+2.0;
          zz[i] = 0.2*(iz)+2.0;
          vx[i] = 0.5*(Math.random()-0.5);
          vy[i] = 0.5*(Math.random()-0.5);
          vz[i] = 0.5*(Math.random()-0.5);
          ax[i] = 0.0;
          ay[i] = 0.0;
          az[i] = 0.0;
          hh[i] = hh0;
          i += 1;
        }
      }
    }

    for (ix=0; ix<10; ix++) {
      for (iy=0; iy<10; iy++) {
        for (iz=0; iz<1; iz++) {
          kind[i] = 1;
          mass[i] = 2.0;
          xx[i] = 2.0+0.866*rBond*ix;
          yy[i] = 2.0+rBond*iy+0.5*rBond*(ix%2);
          zz[i] = 1.0;
          vx[i] = 0.0;
          vy[i] = 0.0;
          vz[i] = 0.0;
          ax[i] = 0.0;
          ay[i] = 0.0;
          az[i] = 0.0;
          hh[i] = hh0;
          i += 1;
        }
      }
    }

	NNp = i;
    bonding();

    for (i=0; i<NNp; i++) {
      energy[i] = 1.0;
	  power[i] = 0.0;
    }

    dividingSection();
    setDensity();

    rotate(theta,fai);
    zSortInit();
  }

/* --------------------------- bonding manager ---------- */

  void bonding() {
    int i,j;
    double r2,rb,rb2;

    rb = rBond*1.2;
	rb2=rb*rb;

    bondInit();

    for (i=0; i<NNp; i++) {
	  if (kind[i]!=0) {
        for (j=i+1; j<NNp; j++) {
		  if (kind[j]!=0) {
            r2 = rSq(i,j);
            if (r2<=rb2) {
              bondReg(i, j); bondReg(j, i);
			}
		  }
        }
      }
    }
  }

  void bondInit() {
    int i,j;

    for (i=0; i<NNp; i++) {
      for (j=0; j<20; j++) {
        bond[i][j]=-1;
      }
    }
  }

  void bondReg(int ii, int jj) {
    int j;

    for (j=0; j<20; j++) {
      if (bond[ii][j]==-1) {
	    bond[ii][j] = jj;
		break;
	  }
	}
  }

  void bondUnreg(int ii, int jj) {
    int j;

    for (j=0; j<20; j++) {
      if (bond[ii][j]==jj) {
	    bond[ii][j] = -1;
		break;
	  }
	}
  }

  int bondQ(int ii, int jj) {
    int j,q;
    double r2,rb,rb2;

    rb = rBond*1.2;
    rb2 = rb*rb;

    q = 0;
    for (j=0; j<20; j++) {
      if (bond[ii][j]==jj) {
        q = 1;
        r2 = rSq(ii,jj);
        if (r2>rb2) {
          bondUnreg(ii,jj); bondUnreg(jj,ii);
          q = 0;
        }
        break;
      }
	}
	return (q);
  }

  double rSq(int i, int j) {
    double xij,yij,zij;

    xij = xx[i] - xx[j];
    if (xij>0.5*xMax) xij -= xMax;
    if (xij<-0.5*xMax) xij += xMax; 
    yij = yy[i] - yy[j];
    if (yij>0.5*yMax) yij -= yMax;
    if (yij<-0.5*yMax) yij += yMax;
    zij = zz[i] - zz[j];
    if (zij>0.5*zMax) zij -= zMax;
    if (zij<-0.5*zMax) zij += zMax; 

    return (xij*xij+yij*yij+zij*zij);
  }

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

  void timeEvolution() {
    int i;

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

/*-------------------------------  particle motion   ---------*/

  void timeStep() {
    int i,j;
	double dtv2,s,rr;

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

    dividingSection();

    setDensity();

    for (i=0; i<NNp; i++) {
	  pressure[i] = 0.5*density[i];
    }

    for (i=0; i<NNp; i++) {
	  if (kind[i]>=0) {
        vx[i] += dtv2*ax[i];
        vy[i] += dtv2*ay[i];
        vz[i] += dtv2*az[i];
        xx[i] += vx[i]*dt;
        yy[i] += vy[i]*dt;
        zz[i] += vz[i]*dt;
      }
    }
    changeCalc();
    for (i=0; i<NNp; i++) {
      if (kind[i]>=0) {
        vx[i] += dtv2*ax[i];
        vy[i] += dtv2*ay[i];
        vz[i] += dtv2*az[i];
        energy[i] += power[i]*dt;
      }
    }

    rr = 1.0;
    for (i=0; i<NNp; 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;
      }
      if (zz[i] < 0.0) { 
        zz[i] += zMax; 
      }
      if (zz[i] > zMax) {
        zz[i] -= zMax; 
      }
    }
  }

  void changeCalc() {
    int i,j,k,i0,i1,j0,j1,k0,k1,ii,jj,kk,iq,iii,jjj,kkk;
    double ai,aj,b,r2,rij,gradW,gradWx,gradWy,gradWz,a,xij,yij,zij;
    double rcc,rcc2,vr,mu,pij;

    rcc = hh0*2.0;
    rcc2 = rcc*rcc;

    for (i=0; i<NNp; i++) {
	  ax[i] = 0.0; ay[i] = 0.0; az[i] = 0.0;
	  power[i] = 0.0;
	  ai = pressure[i]/(density[i]*density[i]);

      i0 = (int)(Nsx*(xx[i]-rcc)/xMax+(double)Nsx)-Nsx;
      i1 = (int)(Nsx*(xx[i]+rcc)/xMax);
      j0 = (int)(Nsy*(yy[i]-rcc)/yMax+(double)Nsy)-Nsy;
      j1 = (int)(Nsy*(yy[i]+rcc)/yMax);
      k0 = (int)(Nsz*(zz[i]-rcc)/zMax+(double)Nsz)-Nsz;
      k1 = (int)(Nsz*(zz[i]+rcc)/zMax);
      for(ii=i0;ii<=i1;ii++) {
        iii = (ii+Nsx)%Nsx;
        for(jj=j0;jj<=j1;jj++) {
          jjj = (jj+Nsy)%Nsy;
          for(kk=k0;kk<=k1;kk++) {
            kkk = (kk+Nsz)%Nsz;
            for(iq=1;iq<=section[iii][jjj][kkk][0];iq++) {
              j = section[iii][jjj][kkk][iq];
              xij = xx[i] - xx[j]; yij = yy[i] - yy[j]; zij = zz[i] - zz[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;
              if (zij>0.5*zMax) zij -= zMax;
              if (zij<-0.5*zMax) zij += zMax; 
              r2=xij*xij+yij*yij+zij*zij;
              if (r2<rcc2) {
                rij = Math.sqrt(r2);
                if (rij>0) {
                  aj = pressure[j]/(density[j]*density[j]);
                  vr = xij*(vx[i]-vx[j])+yij*(vy[i]-vy[j])+zij*(vz[i]-vz[j]);
                  pij = 0;
                  if (vr<0.0) {
                    mu = 0.5*vr/(r2+0.001*0.5*0.5);
                    pij = (-alpha*0.05*mu+0.0*mu*mu)/((density[i]+density[j])/2.0);
                  }
                  b = mass[j]*(ai+aj+pij);
                  gradW = dwwvdr(rij, (hh[i]+hh[j])/2.0);
                  gradWx = gradW*xij/rij;
                  gradWy = gradW*yij/rij;
                  gradWz = gradW*zij/rij;
                  ax[i] += -b*gradWx;
                  ay[i] += -b*gradWy;
                  az[i] += -b*gradWz;
                  power[i] += 0.5*b*((vx[i]-vx[j])*gradWx+(vy[i]-vy[j])*gradWy+(vz[i]-vz[j])*gradWz);
                  if (kind[j]!=0 || kind[i]!=0) {
                    a = force(rij,i,j)/mass[i];
                    ax[i] += a*xij/rij;
                    ay[i] += a*yij/rij;
                    az[i] += a*zij/rij;
                  }
                }
              }
            }
          }
        }
      }
    }
  }

  double force(double r, int ii, int jj) {
    int i,j;
    double ret,rc,ra,ra1,ra2,ram,fa,f,r1,r3,kc;


    if (kind[ii]==0 || kind[jj]==0) {

      rc = hh0*1.2;
      ret = 0.0;
      if (r<rc) {
	    ra = (rc-r)/rc;
	    ret = 1000.0*ra*ra;
		return ( ret );
	  }

	}

    if( kind[ii]>0 && kind[jj]>0 ) {

      kc = 2000.0;
      ra = rBond;
      ra1 = 1.2*ra;
      ra2 = 1.4*ra;
      ram = (ra1+ra2)/2.0;
      fa = 2.0/((ra2-ra1)*(ra2-ra1));

      if (ii<jj) {
        i= ii; j=jj;
      } else {
        i= jj; j=ii;
      }
      if (r>ra2) {
        bondUnreg(i,j); bondUnreg(j,i);
      } 

      f = 0.0;
      if (bondQ(i,j)==0) {
        if (r<ra) {
          f = kc*(r-ra);
        }
      } else {
        if (r<=ra1) {
          f = kc*(r-ra);
        } else if (r>=ra2) {
          f = 0.0;
          bondUnreg(i,j); bondUnreg(j,i);
        } else if (ra1<=r && r<ram) {
          f = kc*(r-ra)*(1.0-fa*(r-ra1)*(r-ra1));
        } else if (ram<=r && r<ra2) {
          f = kc*(r-ra)*fa*(r-ra2)*(r-ra2);
        }
      }
      return ( -f );

	}

    return (0.0);
  }


  void setDensity() {
    int i,j,k,i0,i1,j0,j1,k0,k1,ii,jj,kk,iq,iii,jjj,kkk;
    double s,r2,rij,xij,yij,zij;
    double rcc,rcc2;

    rcc = hh0*2.0;
    rcc2 = rcc*rcc;

    for (i=0; i<NNp; i++) {
      s = 0.0;
      i0 = (int)(Nsx*(xx[i]-rcc)/xMax+(double)Nsx)-Nsx;
      i1 = (int)(Nsx*(xx[i]+rcc)/xMax);
      j0 = (int)(Nsy*(yy[i]-rcc)/yMax+(double)Nsy)-Nsy;
      j1 = (int)(Nsy*(yy[i]+rcc)/yMax);
      k0 = (int)(Nsz*(zz[i]-rcc)/zMax+(double)Nsz)-Nsz;
      k1 = (int)(Nsz*(zz[i]+rcc)/zMax);
      for(ii=i0;ii<=i1;ii++) {
        iii = (ii+Nsx)%Nsx;
        for(jj=j0;jj<=j1;jj++) {
          jjj = (jj+Nsy)%Nsy;
          for(kk=k0;kk<=k1;kk++) {
            kkk = (kk+Nsz)%Nsz;
            for(iq=1;iq<=section[iii][jjj][kkk][0];iq++) {
              j = section[iii][jjj][kkk][iq];
              xij = xx[i] - xx[j]; yij = yy[i] - yy[j]; zij = zz[i] - zz[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;
              if (zij>0.5*zMax) zij -= zMax;
              if (zij<-0.5*zMax) zij += zMax; 
              r2=xij*xij+yij*yij+zij*zij;
              if (r2<rcc2) {
                rij = Math.sqrt(r2);
                s += mass[j]*kernel(rij, (hh[i]+hh[j])/2.0);
              }
            }
          }
        }
      }
      density[i] = s;
    }
  }

/*-------------------------------  registration   -------------*/

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

    for(i=0;i<Nsx;i++) {
      for(j=0;j<Nsy;j++) {
        for(k=0;k<Nsz;k++) {
          section[i][j][k][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;
      k = (int)(Nsz*zz[ip]/zMax); if (k>=Nsz) k = Nsz-1;
      iq = section[i][j][k][0]+1;
      section[i][j][k][0] = iq;
      section[i][j][k][iq] = ip;
    }
  }

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

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

/*----------------------------  smoothed particle ---------*/

  double kernel(double r, double h) {
    double a,q,ret;

    ret = 0.0;
	q = r/h;
	a = 1.0/3.141592/(h*h*h);
	if (q<1.0) {
      ret = a*(1.0-1.5*q*q+0.75*q*q*q);
    } else if (q<2.0) {
      ret = a*0.25*(2.0-q)*(2.0-q)*(2.0-q);
	}
    return ( ret );
  }

  double dwwvdr(double r, double h) {
    double a,q,ret;

    ret = 0.0;
	q = r/h;
	a = 1.0/3.141592/(h*h*h*h);
	if (q<1.0) {
      ret = a*(-3.0*q+2.25*q*q);
    } else if (q<=2.0) {
      ret = -a*0.75*(2.0-q)*(2.0-q);
	}
    return ( ret );
  }

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

}