3次元ガウス分布の自由空間伝搬

はじめに
以下の式で初期条件が与えられる3次元ガウシアンパルスの自由空間伝搬について考えます。
http://www.geocities.jp/mhashimo0724/blog/equation/gauss3d1.gif
波源中心から距離rの点で観測する音圧時間波形の厳密解は、
http://www.geocities.jp/mhashimo0724/blog/equation/gauss3d2.gif
となります。
厳密解は積分となっており、実際に解くには数値積分を用いなければなりません。オーソドックスな台形積分でも解くことはできるけれど、ちょっと高級な数値積分であるRomberg積分を使ってみました。

解説
初めの方で定義した定数の説明です。円周率PIやガウス分布の広がりSIGMAです。KMINとKMAXは積分範囲です。KMINは0というのは式の通りですが、KMAXは本来、無限であるところを700に設定しました。というのもKが大きくなると被積分関数のexpが非常に小さくなるからです。
次に関数を設定します。一つ目の関数は被積分項の関数で、二つ目の関数はRomberg積分用の関数です。
main以降はコメントに書いておきます。

Cプログラム

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
//===define===//
#define PI 3.141592653589793238
#define SIGMA 20e-3/sqrt(2)
#define KMIN 0
#define KMAX 700
#define RR 36e-3*sqrt(3)
//===function===//
double f(double k,double c,double t){
  return pow(k,2)*cos(k*c*t)*sin(k*RR)*exp(-pow(SIGMA*k,2)/2)/(k*RR);
}
double h(int m1){
  return (KMAX-KMIN)/pow(2,m1-1);
}
//===main===//
int main(){
  //===Parameter===//
  int    NT=54;                  //ステップ数
  int    n=15;                    //Romberg積分の精度を決める。
  int    m0,m1,m2,T;           //ループ用変数
  double rho=1.21,K=1.42e5; //室温20度の空気密度と体積弾性率
  double c=sqrt(K/rho);      //音速
  double CFL=1/sqrt(3);     //安定条件
  double ds=4e-3;             //空間格子間隔4mm
  double dt=CFL*ds/c;       //時間刻み
  double tt,tmp;                //時間ttとtmp
  double AA=SIGMA*SIGMA*SIGMA*sqrt(2/PI);    //積分の外の係数
  double pinc[NT];            //時間波形保存用メモリ
  double Rom[n+1][n+1];    //Romberg積分用のメモリ
  //===Romberg Integration===//
  for (T=0;T<NT;T++){
    tt=T*dt;
    Rom[1][1]=0;
    for(m1=2;m1<=n;m1++){
      //==Rom[1][1]--->Rom[n][1]
      tmp=0;
      for(m2=1;m2<=pow(2,m1-2);m2++){
        tmp=tmp+f(KMIN+(2*m2-1)*h(m1),c,tt);
      }
      Rom[m1][1]=(Rom[m1-1][1]+h(m1-1)*tmp)/2;
      //==Rom[n][1]--->Rom[n][n]
      for(m2=2;m2<=m1;m2++){
        Rom[m1][m2]=Rom[m1][m2-1]+(Rom[m1][m2-1]-Rom[m1-1][m2-1])/(pow(4,m2-1)-1);
      }
    }
    pinc[T]=AA*Rom[n][n];
  }
  //===Output Dat Files===//
  for (T=0;T<NT;T++){
    printf("%4.16e\n",pinc[T]);
  }
  return(0);
}

参考
Romberg積分