/** applet No. 1071
 *
 * sound propagation in water - sound wave 3D
 * - multi thread  SW3D :: display - asynchronous
 * - grid 192 x 140 x 140, source point/line/plane
 * - float precision - material table
 *
 * Created by Ikeuchi Mitsuru on November 21 2006.
 * Copyright (c) 2006-2007 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2006.12.08  created
 * ver 0.0.2   2007.06.08  improved code
 *
 */

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

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

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

  Thread th = null; // for run()-paint() thread
  SoundWave3D sw = null;
  DispSW3D dg = new DispSW3D();

  // for event
  Choice ch_source, ch_mat2, ch_view;
  Button bt_reset, bt_startStop, bt_step;
  Scrollbar sc_freq, sc_mag;

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

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

  int sleepTime = 100;
  int dispMode = 0;
  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();

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

    addMouseListener(this);  
    addMouseMotionListener(this);

    sc_freq= new Scrollbar(Scrollbar.HORIZONTAL,100,10,20,110);
    sc_freq.addAdjustmentListener(this);

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

    ch_source = new Choice();
    ch_source.add("point");
    ch_source.add("line");
    ch_source.add("plane");
    ch_source.addItemListener(this);
    ch_source.select("plane");

    ch_mat2 = new Choice();
    ch_mat2.add("water");
    ch_mat2.add("polystyrene");
    ch_mat2.add("glass");
    ch_mat2.add("iron");
    ch_mat2.addItemListener(this);
    ch_mat2.select("glass");

    ch_view = new Choice();
    ch_view.add("p density");
    ch_view.add("v vector");
    ch_view.add("p (at z=0)");
    ch_view.addItemListener(this);
    ch_view.select("p density");

    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(new Label("sound in water"));
    pnl.add(new Label("  "));
    pnl.add(new Label("  "));

    pnl.add(ch_source);
    pnl.add(sc_freq);
    pnl.add(new Label("  "));//pnl.add(ch_mat2);
    pnl.add(new Label("  "));
    pnl.add(sc_mag);
    pnl.add(ch_view);
    add(pnl,"North");
  }

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

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


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

  public void itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_source){
      sw.setSourceMode( ch_source.getSelectedIndex() );
    } else if (ev.getSource() == ch_mat2){
      sw.setMAT2( 1 + ch_mat2.getSelectedIndex() );
    } else if (ev.getSource() == ch_view){
      dispMode = ch_view.getSelectedIndex();
    }
  }

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

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_freq) { 
      sw.setFrequency( 100.0*(double)sc_freq.getValue() );
    } else if(ev.getSource() == sc_mag) { 
      dg.setMag( 0.01*(double)(sc_mag.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;
    }
    dg.setViewDelta( (dgX-dgXb), (dgY-dgYb) );
  }



  // ========================= 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);
  }

  public void offPaint() {
    if (sw != null) {
      dg.plotGraph(gOff, sw, dispMode, thCount);
    }
  }

}



// ================================ display SW3D class ===========

class DispSW3D {

  private Graphics gOff;
  private SoundWave3D sim;

  private int NNx;
  private int NNy;
  private int NNz;
  private double dx;
  private double dy;
  private double dz;
  private double xMax;
  private double yMax;
  private double zMax;
  private double cxx;
  private double cyy;
  private double czz;

  private double dispScale;
  private double viewScale = 1.0;
  private double scale;
  private double magn = 1.0;

  private int xImageLoc = 40;
  private int yImageLoc = 120;

  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 double wkx[] = new double[8]; 
  private double wky[] = new double[8]; 
  private double wkz[] = new double[8]; 

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

  private 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}
  };

  private int xpts[] = new int[200];
  private int ypts[] = new int[200];
  private double zval[] = new double[200];

  private int drawCount = 0;


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

  DispSW3D() {
  }


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

  public void plotGraph(Graphics gg, SoundWave3D sw, int dispMode, int thCount) {
    int simCount,gb;

    gOff = gg;
    sim = sw;
    if (sw==null) return;
    simCount = sim.getCount();
    if (simCount==drawCount) return;
    drawCount = simCount;
    setGlobal();

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

    if (dispMode==0) {
      densityPlot(0,0.1*magn,2); // P density

    } else if (dispMode==1) {
      densityPlot(1,500000.0*magn,4);  // v vector

    } else if (dispMode==2) {
      gridPlot(0, NNz/2, 50.0*magn); // P (@z=0) wave;

    } else if (dispMode==3) {
      ;

    } else if (dispMode==4) {
      ;

    }

    gOff.setColor(Color.black);
    gOff.drawString("source",630/6*0+20,70);
    gOff.drawString("freq="+(int)(sim.getFrequency()+0.5)+" ",630/6*1+10,70);
    //gOff.drawString("object",630/6*2+20,70);
    gOff.drawString("mag="+(int)(magn*100.0+0.5)+"%",630/6*4+10,70);
    gOff.drawString("view",630/6*5+10,70);

    gOff.drawString("grid = "+(int)(xMax*10000+0.5)/100.0+" x "+(int)(yMax*10000+0.5)/100.0+" x "+(int)(zMax*10000+0.5)/100.0+" (cm)",400,100);
    gOff.drawString("t="+(int)(sim.getTime()*1.0e7+0.5)/10.0+" us",400,120);

    //gOff.drawString("scattering object = "+sim.getMAT2String()+"",400,140);
    //gOff.drawString("density = "+sim.getMAT2Density()+"kg/m3",400,160);
    //gOff.drawString("sound velocity = "+sim.getMAT2Velocity()+"m/s",400,180);

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

  }

  private void setGlobal() {

    NNx = sim.getNNx();
    NNy = sim.getNNy();
    NNz = sim.getNNz();
    dx = sim.getdx();
    dy = dx;
    dz = dx;
    xMax = NNx*dx;
    yMax = NNy*dy;
    zMax = NNz*dz;
    cxx = xMax/2.0;
    cyy = yMax/2.0;
    czz = zMax/2.0;
    dispScale = 320.0/xMax;
    scale = dispScale*viewScale;
  }


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

  private void densityPlot(int mode, double mag,int idiv) {
    int i,j,k,ix,iy,ir,fp;
    int iBegin, iEnd, iInc, jBegin, jEnd, jInc, kBegin, kEnd, kInc;
    double x,y,z,dens,col;

    fp = drawBox(0);

    iBegin = 0; iEnd = NNx; iInc = idiv;
    if (wkx[fp]>10.0) {
      iBegin = NNx-1; iEnd = 0; iInc = -idiv;
    }
    jBegin = 0; jEnd = NNy; jInc = idiv;
    if (wky[fp]>10.0) {
      jBegin = NNy-1; jEnd = 0; jInc = -idiv;
    }
    kBegin = 0; kEnd = NNz; kInc = idiv;
    if (wkz[fp]>10.0) {
      kBegin = NNz-1; kEnd = 0; kInc = -idiv;
    }

    for (i=iBegin; ((i<iEnd && iInc==idiv) || (i>=iEnd && iInc==-idiv)) ; i += iInc) {
      x = dx*i;
      for (j=jBegin; ((j<jEnd && jInc==idiv) || (j>=jEnd && jInc==-idiv)) ; j += jInc) {
        y = dy*j;
        for (k=kBegin; ((k<kEnd && kInc==idiv) || (k>=kEnd && kInc==-idiv)) ; k += kInc) {
          z = dz*k;

          if (mode==0) {
            dens = mag*sim.getSoundPressure(i,j,k);
            col = 0.33f;
            if (dens<0) {
              col = 0.01f;
              dens = -dens;
            }
            if (dens>0.999) dens = 0.999;
            ir = (int)(10.0*dens+0.9);
            if (ir>0) plotDot(x,y,z,ir, (float)col);
	
          } else if (mode==1) {
            if (i<=0 || i>=NNx-1 || j<=0 || j>=NNy-1 || k<=0 || k>=NNy-1 ) continue;
            plotVector(x,y,z, 
              mag*sim.getMediaVelocityX(i,j,k),
              mag*sim.getMediaVelocityY(i,j,k),
              mag*sim.getMediaVelocityZ(i,j,k) );

          } else if (mode==2) {
            ;
          }
        }
      }
    }
    drawBox(1);
  }

  private void plotDot(double x, double y, double z, int ir, float col) {
    int ix,iy;
    double ppx,ppy,ppz;

    ppy = cosFi*(y-cyy)+sinFi*(z-czz) + cyy; 
    ppz = -sinFi*(y-cyy)+cosFi*(z-czz) + czz;
    ppx = cosTh*(x-cxx)+sinTh*(ppz-czz) + cxx; 
    //ppz = -sinTh*(x-cxx)+cosTh*(ppz-czz) + czz;

    ix = (int)(scale*ppx)+xImageLoc;
    iy = (int)(scale*ppy)+yImageLoc;

    gOff.setColor(Color.getHSBColor(col,0.9f, 0.9f ));
    gOff.fillOval(ix-ir/2,iy-ir/2,ir,ir);
  }

  private void plotVector(double x, double y, double z, double u, double v, double w) {
    double x2,y2,z2;
    float col;

    x2 = x + u;
    y2 = y + v;
    z2 = z + w;
    if (u>=0.0) {
      col = 0.65f;
    } else {
      col = 0.05f;
    }
    plotLine(x,y,z, x2,y2,z2, col);
  }

  private void plotLine(double x, double y, double z, double x2, double y2, double z2, float col) {
    int ix,iy, ix2,iy2;
    double ppx,ppy,ppz, ppx2,ppy2,ppz2;

    ppy = cosFi*(y-cyy)+sinFi*(z-czz) + cyy; 
    ppz = -sinFi*(y-cyy)+cosFi*(z-czz) + czz;
    ppx = cosTh*(x-cxx)+sinTh*(ppz-czz) + cxx; 
    //ppz = -sinTh*(x-cxx)+cosTh*(ppz-czz) + czz;

    ix = (int)(scale*ppx)+xImageLoc;
    iy = (int)(scale*ppy)+yImageLoc;

    ppy2 = cosFi*(y2-cyy)+sinFi*(z2-czz) + cyy; 
    ppz2 = -sinFi*(y2-cyy)+cosFi*(z2-czz) + czz;
    ppx2 = cosTh*(x2-cxx)+sinTh*(ppz2-czz) + cxx; 
    //ppz2 = -sinTh*(x2-cxx)+cosTh*(ppz2-czz) + czz;

    ix2 = (int)(scale*ppx2)+xImageLoc;
    iy2 = (int)(scale*ppy2)+yImageLoc;

    gOff.setColor(Color.getHSBColor(col,0.9f, 0.9f ));
    gOff.drawLine(ix,iy,ix2,iy2);
  }

  // draw box

  private int drawBox(int imode) {
    int i,gx,gy,g2x,g2y,farPoint;
    double sc;

    if (imode==0) {
      setBox();
      rotateBox();
    }
    sc = dispScale*viewScale;
    farPoint = findFarPoint();

    for (i=0; i<12; i++) {
      gx = (int)(sc*pwkx[boxp[i][0]])+xImageLoc;
      gy = (int)(sc*pwky[boxp[i][0]])+yImageLoc;
      g2x = (int)(sc*pwkx[boxp[i][1]])+xImageLoc;
      g2y = (int)(sc*pwky[boxp[i][1]])+yImageLoc;
      if (imode==0) {
        if (boxp[i][0]==farPoint || boxp[i][1]==farPoint) { 
          gOff.setColor(Color.lightGray);
          gOff.drawLine(gx,gy, g2x, g2y);
        }
      } else if (imode==1) {
        if (boxp[i][0]!=farPoint && boxp[i][1]!=farPoint) { 
          gOff.setColor(Color.darkGray);
          gOff.drawLine(gx,gy, g2x, g2y);
        }
      }
    }
    return farPoint;
  }

  private int findFarPoint() {
    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);
  }

  private void setBox() {
    int i,ix,iy,iz;
    double ballSize;

    i=0;
    ballSize = 0.0*xMax;

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

  private void rotateBox() {
    int i;
 
    for (i=0; i<8; i++) {
      pwky[i] = cosFi*(wky[i]-cyy)+sinFi*(wkz[i]-czz) + cyy; 
      pwkz[i] = -sinFi*(wky[i]-cyy)+cosFi*(wkz[i]-czz) + czz;
      pwkx[i] = cosTh*(wkx[i]-cxx)+sinTh*(pwkz[i]-czz) + cxx; 
      pwkz[i] = -sinTh*(wkx[i]-cxx)+cosTh*(pwkz[i]-czz) + czz;
    }  
  }

  //

  void gridPlot(int mode, int k, double mag) {
    int i,j,gbx,gby;
    int cx,cy,cz;
    double ff,z,px,py,pz,sc;

    cx = NNx/2;
    cy = NNy/2;
    cz = 0;
    sc = 2.5;

    gbx = 40+(int)(cx/sc);
    gby = 60+(int)(cy/sc);

    for(j=0; j<NNy; j += 2) {
      for (i=0; i<NNx; i++) {
        ff = sim.getSoundPressure(i,j,k);
        z = mag*ff + 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;
        zval[i] = ff;
      }
      for (i=0; i<NNx-1; i++) {
        gOff.setColor(Color.getHSBColor(lineColor(i),0.6f,0.9f));
        gOff.drawLine(xpts[i],ypts[i],xpts[i+1],ypts[i+1]);
      }
    }

    for(i=0; i<NNx; i += 2) {
      for (j=0; j<NNy; j++) {
        ff = sim.getSoundPressure(i,j,k);
        z = mag*ff + 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;
        zval[j] = ff;
      }
      for (j=0; j<NNy-1; j++) { 
        gOff.setColor(Color.getHSBColor(lineColor(j),0.6f,0.9f));
        gOff.drawLine(xpts[j],ypts[j],xpts[j+1],ypts[j+1]);
      }
    }
  }

  private float lineColor(int i) {
	
    if (zval[i]*zval[i+1]<0.0) return 0.66f;
    if (zval[i]+zval[i+1]>0.02) return 0.33f;
    if (zval[i]+zval[i+1]<-0.02) return 0.05f;
    return 0.66f;
  }


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

  public void setViewDelta(double xDisplacement, double yDisplacement) {
    viewTheta += 0.5*3.14159/180.0*xDisplacement;
    viewFai += 0.5*3.14159/180.0*yDisplacement;
    cosTh = Math.cos(viewTheta); 
    sinTh = Math.sin(viewTheta);
    cosFi = Math.cos(viewFai);
    sinFi = Math.sin(viewFai);
  }
  public void setMag(double mg) { magn = mg; }
  public double getMag() { return magn; }
}



// ===================================== sound wave 3D ===========

class SoundWave3D extends Thread {

  private float t = 0.0f;
  private float dt = 1.0e-6f;
  private float dx = 1.0e-2f;

  private int NNx = 160+32;
  private int NNy = 108+32;
  private int NNz = 108+32;
  private byte kind[][][] = new byte[NNx][NNy][NNz];
  private float fai[][][][] = new float[3][NNx][NNy][NNz]; // v-potential

  private int pbf=0;
  private int ppr=1;
  private int pnx=2;

  private float frequency = 10000.0f; //Hz
  private float omega = (float)(2.0*Math.PI*frequency);
  private int sourceMode = 2; // 0-point 1-plane source

  private float material[][] = {
  // c(m/s), rho(kg/m3), alpha(atenuation)
    {  340.0f,    1.2f, 0.0001f }, // 0 air
    { 1500.0f, 1000.0f, 0.01f   }, // 1 water
    { 2350.0f, 1056.0f, 0.01f   }, // 2 polystyrene
    { 5100.0f, 2400.0f, 0.01f   }, // 3 glass
    { 5950.0f, 7860.0f, 0.01f   }, // 4 iron

    { 1500.0f, 1000.0f, 10000.0f }, // 5 absorber1
    { 1500.0f, 1000.0f, 40000.0f }, // 6 absorber2
    { 1500.0f, 1000.0f, 200000.0f } // 7 absorber3
  };
  private String materialString[] = {
    "air", "water", "polystyrene", "glass", "iron" };

  private int MAT1 = 1; // water
  private int MAT2 = 3; // scattering object - glass

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


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

  SoundWave3D() {
    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,j,k;

    t = 0.0f;
    count = 0;

    int pbf=0;
    int ppr=1;
    int pnx=2;
   
    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        for (k=0; k<NNz; k++) {
          fai[0][i][j][k] = 0.0f;
          fai[1][i][j][k] = 0.0f;
          fai[2][i][j][k] = 0.0f;
          kind[i][j][k] = (byte)MAT1;

          if ( i<16 || i>NNx-16 || j<16 || j>NNy-16 || k<16 || k>NNz-16) {
            kind[i][j][k] = 5;
            if ( i<8 || i>NNx-8 || j<8 || j>NNy-8 || k<8 || k>NNz-8) {
              kind[i][j][k] = 6;
              if ( i<4 || i>NNx-4 || j<4 || j>NNy-4 || k<4 || k>NNz-4) {
                kind[i][j][k] = 7;
              }
            }
          }
          /*
          if ((i-NNx/2)*(i-NNx/2)+(j-NNy/2)*(j-NNy/2)<(NNy*NNy/16)) {
            velocity[i][j][k] = material[MAT2][0];
            rho[i][j][k] = material[MAT2][1];
            alpha[i][j][k] = material[MAT2][2];
          }
          */
        }
      }
    }
  }


  // ---------------------------------- 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,j,k;
    float a,b,v,al;

    t = t + dt;

    pbf = (pbf+1)%3;
    ppr = (ppr+1)%3;
    pnx = (pnx+1)%3;

    a = dt*dt/(dx*dx);
    for (i=1; i<NNx-1; i++) {
      for (j=1; j<NNy-1; j++) {
        for (k=1; k<NNz-1; k++) {
          v = material[kind[i][j][k]][0];
          al = material[kind[i][j][k]][2];
          b = v*v*a;
          fai[pnx][i][j][k] = 2.0f*fai[ppr][i][j][k] - fai[pbf][i][j][k]
            -al*dt*(fai[ppr][i][j][k] - fai[pbf][i][j][k])
            +b*(fai[ppr][i+1][j][k]+fai[ppr][i][j+1][k]+fai[ppr][i][j][k+1]
            +fai[ppr][i-1][j][k]+fai[ppr][i][j-1][k]+fai[ppr][i][j][k-1]
            -6.0f*fai[ppr][i][j][k]);
        }
      }
    }

    soundSource(sourceMode,t,omega);
  }

  private void soundSource(int sMode,double tt,double w) {
    int ix,iy,iz;
    double th;

    if (t>2.0f/frequency) return;

    if (sMode==0) { // point source
      ix = NNx/5; iy = NNy/2; iz = NNz/2;
      fai[pnx][ix][iy][iz] += (float)(0.025*Math.sin(w*tt)*dt);

    } else if (sMode==1) { // line wave
      ix = NNx/5; iz = NNz/2;
      for (iy=16; iy<NNy-16; iy++) {
        fai[pnx][ix][iy][iz] += (float)(0.001*Math.sin(w*tt)*dt);
      }

    } else if (sMode==2) { // plane wave
      ix = NNx/5;
      for (iy=16; iy<NNy-16; iy++) {
        for (iz=16; iz<NNz-16; iz++) {
          fai[pnx][ix][iy][iz] += (float)(0.0001*Math.sin(w*tt)*dt);
        }
      }

    }
  }


  // ------------------------------------- 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 int getCount() { return count; }

  public void setSourceMode(int sm) { sourceMode = sm; resetSW = 1; }
  public int getSourceMode() { return sourceMode; }
  public void setMAT2(int m2) { MAT2 = m2; resetSW = 1; }
  public double getMAT2Velocity() { return material[MAT2][0]; }
  public double getMAT2Density() { return material[MAT2][1]; }
  public String getMAT2String() { return materialString[MAT2]; }
  public void setFrequency(double f) {
    frequency = (float)f;
    omega = (float)(2.0*Math.PI*frequency);
  }
  public double getFrequency() { return frequency; }

  // get data
  public double getTime() { return t; }
  public int getNNx() { return NNx; }
  public int getNNy() { return NNy; }
  public int getNNz() { return NNz; }
  public double getdx() { return dx; }

  public double getSoundVelocity(int i, int j, int k) { return material[kind[i][j][k]][0]; }
  public double getAttenuation(int i, int j, int k) { return material[kind[i][j][k]][2]; }
  public double getVPotential(int i, int j, int k) { return fai[ppr][i][j][k]; } // m^2/s
  public double getSoundPressure(int i, int j, int k) { // Pa
    return material[kind[i][j][k]][1]*(fai[ppr][i][j][k]-fai[pbf][i][j][k])/dt;
  }
  public double getMediaVelocityX(int i, int j, int k) { // m/s
    return(-0.5*(fai[ppr][i+1][j][k]-fai[ppr][i-1][j][k])/dx);
  }
  public double getMediaVelocityY(int i, int j, int k) { // m/s
    return(-0.5*(fai[ppr][i][j+1][k]-fai[ppr][i][j-1][k])/dx);
  }
  public double getMediaVelocityZ(int i, int j, int k) { // m/s
    return(-0.5*(fai[ppr][i][j][k+1]-fai[ppr][i][j][k-1])/dx);
  }
  public double getSoundEnergyDensity(int i, int j, int k) { // J/m^3
    double p;

    p = material[kind[i][j][k]][1]*(fai[ppr][i][j][k]-fai[pbf][i][j][k])/dt;
    return( p*p/(material[kind[i][j][k]][1]*material[kind[i][j][k]][0]) );
  }

}


/* sound wave 3D
 *
 * velocity potential (fai)
 * particle velocity v = - grad(fai)
 * sound pressure p = rho d(fai)/dt
 * media density rho (kg/m3)
 * sound speed c (m/s)
 *
 * wave propagation equation
 *  d^2(fai)/dr^2 - (1/c^2) d^2(fai)/dt^2 = 0
 *
 * finite difference approximation
 *  d^2(fai)/dt^2 -> { f(k+1) + f(k-1) -2f(k) }/dt^2
 *  d^2(fai)/dr^2 ->
 *    { f(i+1) + f(i-1) -2f(i) }/dx^2 + { f(j+1) + f(j-1) -2f(j) }/dx^2
 *      + {f(k+1) + f(k-1) -2f(k) }/dx^2
 */