/** applet No. 1067
 *
 * interference - sound in air - sound wave 2D
 * - multi thread  SW2D :: display - asynchronous
 * - grid 640 x 640, source point/plane
 *
 * Created by Ikeuchi Mitsuru on December 02 2006.
 * Copyright (c) 2006-2007 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2006.12.02  created
 * ver 0.0.2   2006.12.05  bug fixed
 * ver 0.0.3   2007.06.03  improved code
 *
 */

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

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

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

  Thread th = null; // for run()-paint() thread
  SoundWave2D sw = null;
  DrawGraph2D dg2d = new DrawGraph2D();

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

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

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

  int sleepTime = 50;
  int dispMode = 1;
  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);

    bt_mclear = new Button("MEM clear");
    bt_mclear.addActionListener(this);

    addMouseListener(this);  
    addMouseMotionListener(this);

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

    sc_mypos= new Scrollbar(Scrollbar.HORIZONTAL,180,10,40,620);
    sc_mypos.addAdjustmentListener(this);

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

    ch_source = new Choice();
    ch_source.add("point 1");
    ch_source.add("point 2");
    ch_source.add("point 1+2");
    ch_source.addItemListener(this);
    ch_source.select("point 1+2");

    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 grid");
    ch_view.add("P density");
    ch_view.add("vx-vy");
    ch_view.add("vx-vy & P");
    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 air"));
    pnl.add(new Label("  "));
    pnl.add(bt_mclear);

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

  public void start() {
    if (sw == null) {
      sw = new SoundWave2D();
      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();
      }
    } else if (ev.getSource() == bt_mclear){
      dg2d.oscMemoryClear();
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_freq) { 
      sw.setFrequency( 10.0*(double)sc_freq.getValue() );
    } else if(ev.getSource() == sc_mypos) { 
      dg2d.setMikeYPos( sc_mypos.getValue() );
    } else if(ev.getSource() == sc_mag) { 
      dg2d.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;
    }
    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();
      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) {
      dg2d.plotGraph2D(gOff, sw, dispMode, thCount);
    }
  }

}



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

class DrawGraph2D {

  private Graphics gOff;
  private SoundWave2D sim;

  private int NNx;
  private int NNy;
  private double dx;
  private double magn = 1.0;

  double viewTheta = -15.0*3.14/180.0;
  double viewFai = -60.0*3.14/180.0; //-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 xpts[] = new int[2000];
  private int ypts[] = new int[2000];

  private int micNN = 250;
  private double mic[][] = new double[micNN][4];
  private int micp = 0;
  private int micx[] = { 0, 480, 320, 600 };
  private int micy[] = { 0, 180, 320, 320 };
  private float micColor[] = { 0, 0.01f, 0.30f, 0.60f };  
  private int mikeSW = 1;

  private int drawnCount = 0;


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

  DrawGraph2D() {
  }


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

  public void plotGraph2D(Graphics gg, SoundWave2D sw, int dispMode, int thCount) {
    int simCount,gb;

    gOff = gg;
    sim = sw;

    if (sim==null) return;
    simCount = sim.getCount();
    if (simCount==drawnCount) return;
    drawnCount = simCount;

    NNx = sim.getNNx(); NNy = sim.getNNy();
    dx = sim.getdx();

    microphone();

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

    if (dispMode==0) { // sound Pressure
      gridPlot((NNx/10.0)*magn,NNx/40);

    } else if (dispMode==1) { // sound Pressure density
      densityPlot((NNx/20.0)*magn, 1, NNx/160);

    } else if (dispMode==2) { // vx-vy field
      densityPlot((NNx/20.0)*magn, 2, NNx/160);

    } else if (dispMode==3) { // sound Pressure + vx-vy
      densityPlot((NNx/20.0)*magn, 3, NNx/160);

    } else if (dispMode==4) {
      ;

    }

    gOff.setColor(Color.black);
    gOff.drawString("source",630/6*0+20,70);
    gOff.drawString("freq="+(int)(sim.getFrequency()+0.5)+" Hz",630/6*1+10,70);
    //gOff.drawString("object",630/6*2+20,70);
    gOff.drawString("micY="+micy[1]+"",630/6*3+10,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)(NNx*dx*100+0.5)/100.0+" x "+(int)(NNy*dx*100+0.5)/100.0+" (m)",400,100);
    gOff.drawString("t="+(int)(sim.getTime()*100000.0+0.5)/100.0+" ms",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 SW2D = "+sim.getCount()+" ",20,440);
    gOff.drawString("thread disp = "+thCount+" ",220,440);
  }


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

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

    cx = NNx/2;
    cy = NNy/2;
    cz = 0;
    sc = 1.5*(320.0/NNx);
    ic = 2;

    gbx = 60;
    gby = 20;

    for(j=0; j<NNy; j += inc) {
      for (i=0; i<NNx; i+=ic) { 
        z = mag*sim.getSoundPressure(i,j)+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;
      }
      for (i=0; i<NNx-ic; i+=ic) {
        if (sim.getAttenuation(i,j)>1.0) {
          gOff.setColor(Color.getHSBColor(colorOf(i,j),0.9f,0.5f));
        } else {
          gOff.setColor(Color.getHSBColor(colorOf(i,j),0.9f,0.9f));
        }
        gOff.drawLine(xpts[i],ypts[i],xpts[i+ic],ypts[i+ic]);
      }
    }

    for(i=0; i<NNx; i += inc) {
      for (j=0; j<NNy; j+=ic) { 
        z = mag*sim.getSoundPressure(i,j)+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;
      }
      for (j=0; j<NNy-ic; j+=ic) { 
        if (sim.getAttenuation(i,j)>1.0) {
          gOff.setColor(Color.getHSBColor(colorOf(i,j),0.9f,0.5f));
        } else {
          gOff.setColor(Color.getHSBColor(colorOf(i,j),0.9f,0.9f));
        }
        gOff.drawLine(xpts[j],ypts[j],xpts[j+ic],ypts[j+ic]);
      }
    }
  }

  private void densityPlot(double mag, int mode,int inc) {
    int i,j,gx,gy,ix,iy,ir,ch;
    double cs,cb,sc,vx,vy;

    gx = 30; gy = 100;
    sc = 1.0*(320.0/NNx);
    ir = (int)(sc*inc+0.99999);

    for(i=0+inc/2; i<NNx-inc; i += inc) {
      for (j=0+inc/2; j<NNy-inc; j += inc) { 
        cb = 0.5+0.25*sim.getSoundPressure(i,j); if (cb>0.999) cb = 0.999;
        if (cb<0.0) cb = 0.0;
        if (mode==1 || mode==3) {
          gOff.setColor(Color.getHSBColor( colorOf(i,j), 0.8f, (float)cb ));
        } else if (mode==2) {
          gOff.setColor(Color.getHSBColor( colorOf(i,j), 0.3f, 0.9f ));
        }
        ix = (int)(gx+sc*i);
        iy = (int)(gy+sc*j);
        gOff.fillRect(ix,iy,ir,ir);
      }
    }

    if (mikeSW==1) {
      for (ch=1; ch<4; ch++) {
        gOff.setColor(Color.getHSBColor( micColor[ch], 0.8f, 0.8f ));
        ix = (int)(gx+sc*micx[ch]);
        iy = (int)(gy+sc*micy[ch]);
        gOff.fillRect(ix-2,iy-2,4,4);
      }
    }

    if (mode==2 || mode==3) {
      for(i=0+inc; i<NNx-inc*2; i += inc*2) {
        for (j=0+inc; j<NNy-inc*2; j += inc*2) {
          vx = sim.getMediaVelocityX(i,j);
          vy = sim.getMediaVelocityY(i,j);
          if (vx*vx+vy*vy<1.0e-16) continue;
          if (vx>0) {
            gOff.setColor(Color.red);
          } else {
            gOff.setColor(Color.blue);
          }
          ix = (int)(gx+sc*(i+inc));
          iy = (int)(gy+sc*(j+inc));
          gOff.drawLine(ix,iy,ix+(int)(mag*500.0*vx),iy+(int)(mag*500.0*vy) );
        }
      }
    }

    if (mikeSW==1) {
      plotMicrophone(200.0*magn);
    }
  }

  private float colorOf(int i, int j) {
    float c;

    c = 0.33f;
    if (sim.getSoundVelocity(i,j)>400.0) c = 0.10f;
    if (sim.getAttenuation(i,j)>1.0) c = c + 0.1f;
    return c;
  }

  //

  private void microphone() {
    if (sim.getStartSW()==1) {
      micp = (micp+1)%micNN;
      mic[micp][0] = sim.getTime();
      mic[micp][1] = sim.getSoundPressure(micx[1],micy[1]);
      mic[micp][2] = sim.getSoundPressure(micx[2],micy[2]);
      mic[micp][3] = sim.getSoundPressure(micx[3],micy[3]);
    }
  }

  private void plotMicrophone(double mag) {
    int n,ch,i,ip1,ip2,gx1,gy1,gx2,gy2;

    gOff.setColor(Color.black);
    gOff.drawString("microphone (2.5 ms/div)",420,230);
    gOff.setColor(Color.lightGray);
    for (i=0; i<=200; i+=50) {
      gOff.drawLine(400+i,240,400+i,440);
    }
    for (i=0; i<=200; i+=50) {
      gOff.drawLine(400,240+i,600,240+i);
    }

    for (ch=1; ch<4; ch++) {
      gOff.setColor(Color.getHSBColor( micColor[ch], 0.8f, 0.6f ));
      for (i=1; i<micNN; i++) {
        ip1 = (micp+i)%micNN;
        ip2 = (micp+i+1)%micNN;
        gx1 = 400+(int)(20000.0*(mic[micp][0]-mic[ip1][0]));
        gy1 = 340-(int)(mag*mic[ip1][ch]);
        gx2 = 400+(int)(20000.0*(mic[micp][0]-mic[ip2][0]));
        gy2 = 340-(int)(mag*mic[ip2][ch]);
        if (gx1<400 || gx1>600 || gx2<400 || gx2>600) continue;
        gOff.drawLine(gx1,gy1,gx2,gy2);
      }
    }
  }

  private void memoryClear() {
    int i;

    micp = 0;
    for (i=0; i<micNN; i++) {
      mic[i][0] = -1.0;
      mic[i][1] = 0.0;
      mic[i][2] = 0.0;
      mic[i][3] = 0.0;
    }
  }


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

  public void setMikeYPos(int yp) { micy[1] = yp; }

  public void oscMemoryClear() { memoryClear(); }
}



// ===================================== sound wave 2D ===========

class SoundWave2D extends Thread {

  private double t = 0.0;
  private double dt = 1.0e-5;
  private double dx = 1.0e-2;

  private int NNx = 640;
  private int NNy = 640;
  private double fai[][][] = new double[3][NNx][NNy]; // v-potential
  private double velocity[][] = new double[NNx][NNy]; // sound speed
  private double alpha[][] = new double[NNx][NNy]; // attenuation
  private double rho[][] = new double[NNx][NNy]; // density  

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

  private double frequency = 500.0; //Hz
  private double omega = 2.0*Math.PI*frequency;
  private int sourceMode = 2; // 0-point1,1-point2,2-point1+2

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

  private int MAT1 = 0; // air
  private int MAT2 = 3; // scattering object

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

  SoundWave2D() {
    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;

    t = 0.0;
    count = 0;

    int pbf=0;
    int ppr=1;
    int pnx=2;
   
    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        fai[0][i][j] = 0.0;
        fai[1][i][j] = 0.0;
        fai[2][i][j] = 0.0;
        velocity[i][j] = material[MAT1][0];
        rho[i][j] = material[MAT1][1];
        alpha[i][j] = material[MAT1][2]; // atenuation

        if ( i<32 || i>NNx-32 || j<32 || j>NNy-32) {
          alpha[i][j] = 1000.0;
          if ( i<16 || i>NNx-16 || j<16 || j>NNy-16) {
            alpha[i][j] = 3000.0;
            if ( i<8 || i>NNx-8 || j<8 || j>NNy-8) {
              alpha[i][j] = 10000.0;
            }
          }
        }
        /*
        if ((i-NNx/2)*(i-NNx/2)+(j-NNy/2)*(j-NNy/2)<NNy*NNy/64) {
          velocity[i][j] = material[MAT2][0];
          rho[i][j] = material[MAT2][1];
          alpha[i][j] = 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;
    double a,b;

    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++) {
        b = velocity[i][j]*velocity[i][j]*a;
        fai[pnx][i][j] = 2.0*fai[ppr][i][j] - fai[pbf][i][j]
          -alpha[i][j]*dt*(fai[ppr][i][j] - fai[pbf][i][j])
          +b*(fai[ppr][i+1][j]+fai[ppr][i][j+1]+fai[ppr][i-1][j]+fai[ppr][i][j-1]-4.0*fai[ppr][i][j]);
      }
    }

    soundSource(sourceMode,t,omega);
  }

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

    //if (t>2.0/frequency) return;

    if (sMode==0) { // point 1
      a = 10.0*Math.sin(w*tt)*dt;

      ix = NNx/2; iy = 3*NNy/8;
      fai[pnx][ix][iy] += a;

    } else if (sMode==1) { // point 2
      a = 10.0*Math.sin(w*tt)*dt;

      ix = NNx/2; iy = 5*NNy/8;
      fai[pnx][ix][iy] += a;

    } else if (sMode==2) { // point 1+2
      a = 10.0*Math.sin(w*tt)*dt;

      ix = NNx/2; iy = 3*NNy/8;
      fai[pnx][ix][iy] += a;

      ix = NNx/2; iy = 5*NNy/8;
      fai[pnx][ix][iy] += a;
    }
  }


  // ------------------------------------- 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 = f;
    omega = 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 double getdx() { return dx; }

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

    p = rho[i][j]*(fai[ppr][i][j]-fai[pbf][i][j])/dt;
    return( p*p/(rho[i][j]*velocity[i][j]) );
  }

}

/* sound wave 2D
 *
 * 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
 *
 */