/* ball flow in SI unit - periodic smoothed particle hydrodynamics 2D */
/* coded by Ikeuchi Mitsuru   */
/* ver 0.0.1   2005.08.13     */

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

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

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

  Button bt_reset, bt_startStop;
  Choice ch_view;
  Scrollbar sc_vx0, sc_turns, sc_alpha, sc_dispScale;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 30;

  int dispMode = 1;
  int incRate = 10;
  int resetSW = 0;
  int started = 1;
  double dispScale = 1.0;
  double scale = 20.0;
  double alpha = 0.1;
  double thermalDiffusivity = 3.4e-5; /* m2/s */
  double vx0 = 50.0; /* m/s */
  double turns = 10.0; 

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

  double xMax = 20.0; /* m */
  double yMax = 10.0; /* m */

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

  int NNp = 5000;
  int kind[] = new int [NNp];
  int bond[][] = new int[NNp][8];
  double xx[] = new double [NNp];
  double yy[] = new double [NNp];
  double vx[] = new double [NNp];
  double vy[] = new double [NNp];
  double ax[] = new double [NNp];
  double ay[] = new double [NNp];

  double hh[] = new double [NNp];
  double mass[] = new double [NNp];
  double cp[] = 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 = 40;
  int Nsy = 20;
  int section[][][] = new int[Nsx][Nsy][50];

/* ----------------------------- 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_startStop= new Button("start/stop");
    bt_startStop.addActionListener(this); 

    ch_view = new Choice();
    ch_view.addItem("particle");
    ch_view.addItem("velocity");
    ch_view.addItem("density");
    ch_view.addItem("temperature");
    ch_view.addItem("pressure");
    ch_view.addItemListener(this);
    ch_view.select("velocity");

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

    sc_turns= new Scrollbar(Scrollbar.HORIZONTAL,100,10,0,210);
    sc_turns.addAdjustmentListener(this);

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

    sc_dispScale= new Scrollbar(Scrollbar.HORIZONTAL,100,10,50,210);
    sc_dispScale.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(sc_vx0);
    pnl.add(sc_turns);
    pnl.add(sc_alpha);
    pnl.add(ch_view);
    add(pnl,"North");

    setInitialCondition();
  }

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

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

  public void 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_startStop){
	  if (started==0) {
        started = 1;
	  } else {
        started = 0;
	  }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_vx0) { 
      vx0 = 1.0*sc_vx0.getValue();
    } else if (ev.getSource() == sc_turns) { 
      turns = 0.1*sc_turns.getValue();
    } else if (ev.getSource() == sc_alpha) { 
      alpha = 0.01*sc_alpha.getValue();
    } else if (ev.getSource() == sc_dispScale) { 
      dispScale = 0.01*sc_dispScale.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;
    }
    xShift += (1.0/scale/dispScale)*(dgX-dgXb);
    yShift += (1.0/scale/dispScale)*(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 gb,gbx,gby;
    double sc;

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

    sc = scale*dispScale;
    gbx = 30+(int)(sc*xShift);
    gby = 80+(int)(sc*yShift);

    gOff.setColor(Color.gray);
    gOff.drawRect(gbx,gby,(int)(sc*xMax+0.5),(int)(sc*yMax+0.5));

    if (dispMode==0) {

      ovalPlot(gbx,gby,sc);

	} else if (dispMode==1) {

      linePlot(gbx,gby,sc,0.2);

	} else if (dispMode==2) {

      fieldPlot(gbx,gby,sc,density,0.1);

    } else if (dispMode==3) {

      tempPlot(gbx,gby,sc,273.0,1.0/200.0);

    } else if (dispMode==4) {

      fieldPlot(gbx,gby,sc,pressure,0.00005);

    }

    gOff.setColor(Color.black);
    gOff.drawString("t="+(int)(t*1000*100.0+0.5)/100.0+" ms",630/6*0+10,40);
	if (started==1) {
      gOff.drawString("started",630/6*1+10,40);
	} else {
	  gOff.drawString("stopped",630/6*1+10,40);
	}
    gOff.drawString("vx0="+(int)(vx0*10.0)/10.0+" m/s",630/6*2+10,40);
    gOff.drawString("turns="+(int)(turns*10.0)/10.0+" /s",630/6*3+10,40);
    gOff.drawString("alpha="+(int)(alpha*100.0)/100.0+" ",630/6*4+10,40);
    gOff.drawString("view",630/6*5+10,40);

    gOff.drawString("maxSec="+maxSection()+" ",10,60);
    gOff.drawString("Lx=20m, Ly=10m, ball flow",200,60);
  }

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

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

    ir = (int)(2.0*hh0*sc+0.5);
    for (i=0; i<NNp; i++) {
	  if (kind[i]>0) {
	    gOff.setColor(Color.cyan);
	  } else if (kind[i]==0) {
	    gOff.setColor(Color.green);
	  } else {
		gOff.setColor(Color.orange);
	  }
      gOff.drawOval(gx+(int)(xx[i]*sc+0.5)-ir/2,gy+(int)(yy[i]*sc+0.5)-ir/2, ir, ir);
    }
	
	gOff.setColor(Color.red);
    for (i=0; i<NNp; i++) {
	  if (kind[i]<0) {
	    for (j=0; j<NNp; j++) {
		  if (bondQ(i,j)==1) {
		    if ((xx[i]-xx[j])*(xx[i]-xx[j])+(yy[i]-yy[j])*(yy[i]-yy[j])<2.0*rBond*rBond) {
		      gOff.setColor(Color.red);
		      gOff.drawLine(gx+(int)(xx[i]*sc+0.5),gy+(int)(yy[i]*sc+0.5),gx+(int)(xx[j]*sc+0.5),gy+(int)(yy[j]*sc+0.5));
		    }
		  }
		}
	  }
    }
  }

  void linePlot(int gx, int gy, double sc, double mag) {
    int i,ix,iy,ir;

    ir = (int)(2.0*hh0*sc+0.5);
    for (i=0; i<NNp; i++) {
      ix = gx+(int)(xx[i]*sc+0.5);
      iy = gy+(int)(yy[i]*sc+0.5);
      if (vx[i]>0.0) {
        gOff.setColor(Color.getHSBColor(0.66f,0.8f,0.8f));
      } else {
        gOff.setColor(Color.getHSBColor(0.01f,0.8f,0.8f));
      }
      gOff.drawLine(ix,iy,ix+(int)(vx[i]*dispScale*mag+0.5),iy+(int)(vy[i]*dispScale*mag+0.5));
      if (kind[i]<0) {
	    gOff.setColor(Color.orange);
		gOff.drawOval(ix-ir/2,iy-ir/2, ir, ir);
      } else if (kind[i]==0) {
	    gOff.setColor(Color.green);
		gOff.drawOval(ix-ir/2,iy-ir/2, ir, ir);
      }
    }
  }

  void fieldPlot(int gx, int gy, double sc, double field[], double mag) {
    int i,j,ir;
	double c;

    ir = (int)(2.0*hh0*sc+0.5);
    for (i=0; i<NNp; i++) {
	  c = mag*field[i]; if (c<0.0) c = 0.0;
	  if (c>1.0) c = 1.0;
	  gOff.setColor(Color.getHSBColor((float)(0.66-0.66*c),0.8f,0.8f));
      gOff.drawOval(gx+(int)(xx[i]*sc+0.5)-ir/2,gy+(int)(yy[i]*sc+0.5)-ir/2, ir, ir);
    }
  }

  void tempPlot(int gx, int gy, double sc, double offset, double mag) {
    int i,j,ir;
	double c;

    ir = (int)(2.0*hh0*sc+0.5);
    for (i=0; i<NNp; i++) {
	  c = mag*(energy[i]/(mass[i]*cp[i])-offset); if (c<0.0) c = 0.0;
	  if (c>1.0) c = 1.0;
	  gOff.setColor(Color.getHSBColor((float)(0.66-0.66*c),0.8f,0.8f));
      gOff.drawOval(gx+(int)(xx[i]*sc+0.5)-ir/2,gy+(int)(yy[i]*sc+0.5)-ir/2, ir, ir);
    }
  }

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

  void setInitialCondition() {
    int i,j,ix,iy,ith,nx,ny,nth;
    double x,y,th,x0,y0,r0,r2;

    t = 0.0;
    i = 0;

    x0 = 0.25*xMax;
	y0 = 0.5*yMax;
	r0 = 1.0;
	r2 = (r0+0.3)*(r0+0.3);

    nx = (int)(xMax/0.3);
    ny = (int)(yMax/0.3/0.866);

    for (ix=0; ix<nx; ix++) {
      for (iy=0; iy<ny; iy++) {
	    x = 0.15+0.3*(ix)+0.15*(iy%2);
	    y = 0.15+0.3*(0.866*iy);
		if ((x-x0)*(x-x0)+(y-y0)*(y-y0)>r2) {
          kind[i] = 1;
          mass[i] = (0.029/0.0224)/4.0;
          cp[i] = 1000.0;
	      xx[i] = x;
	      yy[i] = y;
	      vx[i] = 0.5*(Math.random()-0.5)+vx0;
	      vy[i] = 0.5*(Math.random()-0.5);
	      ax[i] = 0.0;
	      ay[i] = 0.0;
	      hh[i] = hh0;
          i += 1;
		}
	  }
    }

    nth = (int)(2.0*3.1416*r0/0.3);
    for (ith=0; ith<nth; ith++) {
	  th = (2.0*3.1416/nth)*ith;
      kind[i] = 0;
      mass[i] = (0.029/0.0224)/4.0;
      cp[i] = 1000.0;
      xx[i] = r0*Math.cos(th)+x0;
      yy[i] = r0*Math.sin(th)+y0;
      vx[i] = 0.0;
      vy[i] = 0.0;
      ax[i] = 0.0;
      ay[i] = 0.0;
      hh[i] = hh0;
      i += 1;
    }

    NNp = i;
    bonding();

    for (i=0; i<NNp; i++) {
      energy[i] = mass[i]*cp[i]*300;
      power[i] = 0.0;
    }

    dividingSection();
    setDensity();

  }

  void rotate(double dth) {
    int i;
	double x,y,x0,y0;

    x0 = 0.25*xMax;
	y0 = 0.5*yMax;

    for (i=0; i<NNp; i++) {
      if (kind[i]==0) {
        x = xx[i]-x0;
		y = yy[i]-y0;
		xx[i] = x*Math.cos(dth)-y*Math.sin(dth)+x0;
		yy[i] = x*Math.sin(dth)+y*Math.cos(dth)+y0;
      }
    }
  }

/* --------------------------- 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<8; j++) {
        bond[i][j]=-1;
      }
    }
  }

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

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

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

    for (j=0; j<8; 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<8; 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;

    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; 
    return (xij*xij+yij*yij);
  }

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

  void timeEvolution() {
    int i;

    if (resetSW==1) {
      setInitialCondition();
      resetSW = 0;
    }
    if (started==1) {
      for (i=0; i<incRate; i++) {
        timeStep();
		rotate(2.0*3.1416*turns*dt);
      }
    }
  }

/*-------------------------------  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++) {
      if (kind[i]>0) {
        pressure[i] = 8.31*density[i]*(energy[i]/(mass[i]*cp[i]));

      } else {
        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];
        xx[i] += vx[i]*dt;
        yy[i] += vy[i]*dt;
	  }
    }
    changeCalc();
    for (i=0; i<NNp; i++) {
      if (kind[i]!=0) {
        vx[i] += dtv2*ax[i];
        vy[i] += dtv2*ay[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 (xx[i] < 0.3) {
        vx[i] = vx0;
      }
    }

  }

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

    cv = 400.0;
    rcc = hh0*2.0;
    rcc2 = rcc*rcc;

    for (i=0; i<NNp; i++) {
      ax[i] = 0.0; ay[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);
      for(ii=i0;ii<=i1;ii++) {
        iii = (ii+Nsx)%Nsx;
        for(jj=j0;jj<=j1;jj++) {
		  jjj = (jj+Nsy)%Nsy;
          for(iq=1;iq<=section[iii][jjj][0];iq++) {
            j = section[iii][jjj][iq];

            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<rcc2) {
              rij = Math.sqrt(r2);
              if (rij>0) {
			    if (kind[j]>0 && kind[i]>0) {
                aj = pressure[j]/(density[j]*density[j]);
                  vr = xij*(vx[i]-vx[j])+yij*(vy[i]-vy[j]);
                  pij = 0;
                  if (vr<0.0) {
                    mu = 0.5*vr/(r2+0.001*0.5*0.5);
                    pij = (-alpha*cv*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;
                  ax[i] += -b*gradWx;
				  ay[i] += -b*gradWy;
                  lapW = d2wwvdr2(rij, (hh[i]+hh[j])/2.0);
                  power[i] += 0.5*b*((vx[i]-vx[j])*gradWx+(vy[i]-vy[j])*gradWy) + thermalDiffusivity*mass[j]/density[j]*(energy[j]-energy[i])*lapW;
                } else {
                  a = force(rij,i,j)/mass[i];
                  ax[i] += a*xij/rij;
                  ay[i] += a*yij/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 = 50000000.0*ra*ra;
		return ( ret );
	  }

    }

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

      kc = 50000000.0;
      ra = rBond;
      ra1 = 1.2*ra;
      ra2 = 1.5*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);
        return ( 0.0 );
      } 

      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 boundaryForce(int i, double rc, double k) {
    double r;

    if (xx[i] < rc) {
	  r = rc-xx[i];
      ax[i] += k*r*r; 
    }
    if (xx[i] > xMax-rc) {
	  r = xx[i]-(xMax-rc);
      ax[i] += -k*r*r; 
    }
    if (yy[i] < rc) {
	  r = rc-yy[i];
      ay[i] += k*r*r; 
    }
    if (yy[i] > yMax-rc) {
	  r = yy[i]-(yMax-rc);
      ay[i] += -k*r*r; 
    }
  }
*/

  void setDensity() {
    int i,j,i0,i1,j0,j1,ii,jj,iq,iii,jjj,ikind;
    double s,rij,r2,xij,yij;
    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);
      for(ii=i0;ii<=i1;ii++) {
        iii = (ii+Nsx)%Nsx;
		ikind = kind[iii];
        for(jj=j0;jj<=j1;jj++) {
          jjj = (jj+Nsy)%Nsy;
          for(iq=1;iq<=section[iii][jjj][0];iq++) {
            j = section[iii][jjj][iq];
			if (kind[j]==ikind) {
              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<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,ip,iq;

    for(i=0;i<Nsx;i++) {
      for(j=0;j<Nsy;j++) {
        section[i][j][0] = 0;
      }
    }

    for(ip=0;ip<NNp;ip++) {
      i = (int)(Nsx*xx[ip]/xMax); if (i>=Nsx) i = Nsx-1;
	  if (i<0) i = 0;
      j = (int)(Nsy*yy[ip]/yMax); if (j>=Nsy) j = Nsy-1;
	  if (j<0) j = 0;
      iq = section[i][j][0]+1;
      section[i][j][0] = iq;
      section[i][j][iq] = ip;
    }
  }

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

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

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

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

    ret = 0.0;
	q = r/h;
	a = 10.0/(7.0*3.141592)/(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 = 10.0/(7.0*3.141592)/(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 );
  }

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

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

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

}