衝撃波を描く

衝撃波の動画を描いて見ようと思います.下記に示すような画像を何枚も描くことになります.

下記に,そのプログラムを示します。

　まず,パルス波を1秒間に1個の割合でだすことにします.8秒間の動画を作る予定ですので,
波の数は8個です.　 8個の波を波紋の中心の位置を変えながら描いていきます. 標準では1秒間に
15コマの動画となりますので, 8x15枚の画像を描けば,8秒間の動画ができあがります。

keika=8;(* 経過時間 *)
bunkatu=keika*15;(* 描く画像の総枚数 *)
kankaku=keika/bunkatu;(* 一枚一枚の画像の時間間隔 *)
naminokazu=8;(* 波の数の総数 *)

　まず,それぞれの波をどの時刻から描き始めるかを決定します. kaishilist という
リストにそれぞれの波を描き出す時刻を入れておきます. 全体の時間(keika)を波の個数で
割って、均等に（下記では1秒間隔で）描き始める時刻を設定します。

kaishilist=Table[keika/naminokazu*(i-1),{i,naminokazu}];

    　音源(船でしょうか)の速度は波の速度の1.414倍で,x軸方向に移動します. (cx, cy)は
    波紋の中心の位置です。それぞれの波の開始時刻に合わせて中心をずらしています。
    
    cx=x0+1.414*velo*kaishilist[[j]];
  cy=y0;

    
  　波の速度は前と同じく, 4 m/sec です.周期も前と同様に1Hz です。
  
    　nagasa は計算する領域を表します. 下記では50 m x 50 m の領域を計算します.
    kizami はメッシュの切り方を指定しています. 50mを100に分割します。
    これを増加させると計算時間が長くなります。
    
  nagasa=50;
  kizami=100;
  dl=nagasa/kizami;

    snaplist は計算したその地点での変位を保存しておくリストです.
    このリストの内容はプログラムの最後の方で,test.txt に書き出しています。
    
     snaplist >>"test.txt";
   
    gsnaplist は描いた画像を保存しておくリストです. これも最後にtest.avi として
書き出しています。

  Export["test.avi",gsnaplist,"AVI"];

  計算するときに一枚一枚の画像が表示されるのですが,計算途中の画像表示を
DisplayFunction → Identity を用いて表示を抑制すればメモリ-を節約できます。

      後は,それぞれの時刻での各波の波形を計算して,重ね合わせるだけです.
    Doル-プを用いて計算しています.それぞれの波は開始の時刻と中心の座標が異なるので,
そこに注意します。2重のループになっていますが、外側のループは各時刻のループです。
つまり、時間の経過を示します。中のループの役目は、各波の描きはじめの予定時刻と
現在の時刻を比較して、書き始めの時刻が過去であるかどうか判断し、そうであれば描く
というものです。

Remove[" Global `*"];
Timing[
  frq=1;(* 周波数 *)
  shuki=1/frq;
  r[x_,y_,x0_,y0_]=Sqrt[(x-x0)^2 + (y-y0)^2];
  warur[x_,y_,x0_,y0_]=
  　　　　If[Sqrt[(x-x0)^2 + (y-y0)^2]==0,1,
  　　　　(1-Exp[-Sqrt[(x-x0)^2 + (y-y0)^2]])/Sqrt[(x-x0)^2 + (y-y0)^2]];
  
  velo=4;
  hachou=velo/frq;(* 波長 *)
  
  keika=8;(* 経過時間 *)
  bunkatu=keika*15;(* 描く画像の枚数 *)
  kankaku=keika/bunkatu;(* 一枚一枚の画像の時間間隔 *)
  naminokazu=8;(* 波の数の総数 *)
  
  t0list=Table[kankaku*(i-1),{i,bunkatu+1}];
  kaishilist=Table[keika/naminokazu*(i-1),{i,naminokazu}];
  (* それぞれのパルス波を描き始める時刻 *)
  
  nagasa=50;(* 計算する範囲, 50m x 50 m の範囲を計算する *)
  kizami=100;(* 上記のnagasaを何個に分割して計算するか, その分割数. *)
  dl=nagasa/kizami;
  
  snaplist={Null};
  gsnaplist={Null};
  
  nami[t_,x_,y_,x0_,y0_]=
    If[(velo*(t-shuki)<r[x,y,x0,y0])&&(r[x,y,x0,y0]<velo*t),
    Sin[-2*Pi*(r[x,y,x0,y0]-velo*(t-shuki))/hachou]*warur[x,y,x0,y0],0];
  
  x0=-10;(* 最初に描く波の中心の座標(x0,y0) *)
  y0=0;
  
  
  
  Do[
    ima=t0list[[k]];(* 現在の時刻, イマ *)
    amplist=Table[0,{i,kizami+1},{j,kizami+1}];
    
    Do[
      If[kaishilist[[j]]<ima,(* 現在の時刻とそれぞれの波の描きはじめの時刻とを比較する. *)
          cx=x0+1.414*velo*kaishilist[[j]];(* (cx,cy)はそれぞれの波を描き始める中心の座標 *)
          cy=y0;
          t=ima-kaishilist[[j]];(* それぞれの波の開始時刻から何秒経過しているかを求める *)
          tmplist=Table[
          　N[nami[t,dl*(i-(kizami/2)),dl*(j-(kizami/2)),cx,cy]],
          　　　　　　　　　　　　　　　　　　　{i,0,kizami},{j,0,kizami}];
          
          amplist=amplist+tmplist,Continue[]];
      ,{j,naminokazu}];
    amplist=Transpose[amplist];
    snaplist=Append[snaplist,amplist];
    
    Print[" k = ";k];
    gsnap=ListPlot3D[amplist,PlotRange→{-1,1},Mesh→False,Boxed→False,Axes→False];
    gsnaplist=Append[gsnaplist,gsnap];
    ,{k,1,bunkatu+1,1}];
  
  snaplist=Drop[snaplist,1];
  gsnaplist=Drop[gsnaplist,1];
  
  SetDirectory["d:\\today"];
  snaplist >>"test.txt";
  Export["test.avi",gsnaplist,"AVI"];
  ]

この様にして,この動画が出来上がりました.

動画を見ると,波面の共通接線が強めあっている様子は分かりますね.

計算に用いたプログラムはこれです.

パラメ-タ-を変更して、波の数を増やしてみると、衝撃波の内側の円形の波が打ち消されて、
波の数の増加と共に消えていきます。下記の図は波の数を1秒間に4個としたものです。

音源の速度を波の速度よりも小さくするとドップラー効果となります。ただし、この場合には波の数を
減らす必要があります。1秒間に1個の波だと波どおしが重なってしまうので、2秒間に1個ぐらいに
した方が良いです。

Created by Mathematica  (May 31, 2007)