vid=VideoWriter('outmov');
open(vid)
%
%　変数の配列領域確保
%
p1=zeros(50,50);p2=zeros(50,50);p3=zeros(50,50);
u1=zeros(50,50);u2=zeros(50,50);w1=zeros(50,50);w2=zeros(50,50); 
vin=zeros(200);
%     
%　定数～pai2=２π，vpl=vl/vp，それ以外は吸収境界に必要なもの
%
       pai2=3.14*2.0; 
       vpl=0.5; 
       dd=199./200.; 
       a=2.*(vpl-1.)/(vpl+1.); 
       b=((vpl-1.)^2)/((vpl+1.)^2); 
       c=2.*(vpl-1.)/(vpl+1.); 
       d=dd*2.; 
       e=dd^2; 
%
%　入力値（振幅１の音圧１周期分）
% 
       for it=1:15 
            vin(it)=sin(it/15.*pai2);
       end
%
%　計算開始
%
 for it=1:150
%
%　圧力 p の計算
%
  for ix=2:41 
   for iz=2:41 
    p1(ix,iz)=p2(ix,iz)-vpl*(u2(ix,iz)-u2(ix-1,iz)+w2(ix,iz)-w2(ix,iz-1)); 
   end
  end 
%
%　描画とビデオフレームへの取り込み
%
 figure;
 surf(p1);
 axis([2 40 2 40 -0.2 0.2]);
 view(30,60);
 frame=getframe(gcf);
 writeVideo(vid,frame);
 close;
%
%　吸収境界における圧力の計算
%
 ix=41; 
 for iz=2:41 
  p1(ix,iz)=a*(p1(ix-1,iz)-p2(ix,iz))-b*(p1(ix-2,iz)-2.*p2(ix-1,iz)...
     +p3(ix,iz)); 
  p1(ix,iz)=p1(ix,iz)-c*(p2(ix-2,iz)-p3(ix-1,iz))+d*p2(ix-1,iz)...
    -e*p3(ix-2,iz); 
 end
% 
 iz=41; 
 for ix=2:41 
  p1(ix,iz)=a*(p1(ix,iz-1)-p2(ix,iz))-b*(p1(ix,iz-2)-2.*p2(ix,iz-1)+...
     p3(ix,iz)); 
  p1(ix,iz)=p1(ix,iz)-c*(p2(ix,iz-2)-p3(ix,iz-1))+d*p2(ix,iz-1)-e*...
    p3(ix,iz-2); 
 end
%
%　中心に圧力を与える（音源）
% 
 p1(2,2)=vin(it); 
%
%　新旧圧力値の入れ替え
% 
 for ix=2:41 
  for iz=2:41 
   p3(ix,iz)=p2(ix,iz); 
   p2(ix,iz)=p1(ix,iz); 
  end 
 end 
%
%　x 方向粒子速度 u の計算
% 
 for ix=2:40 
  for iz=2:41 
   u1(ix,iz)=u2(ix,iz)-vpl*(p2(ix+1,iz)-p2(ix,iz)); 
  end 
 end 
%
%　z 方向粒子速度 w の計算
%
 for ix=2:41 
  for iz=2:40 
   w1(ix,iz)=w2(ix,iz)-vpl*(p2(ix,iz+1)-p2(ix,iz)); 
  end 
 end 
%
%　対称条件の設定
% 
 for iz=2:41 
  u1(1,iz)=-u1(2,iz); 
 end 
% 
 for ix=2:41

  w1(ix,1)=-w1(ix,2); 
 end 
% 
%　新旧粒子速度値の入れ替え
%
 for ix=2:41 
  for iz=2:41 
   u2(ix,iz)=u1(ix,iz); 
   w2(ix,iz)=w1(ix,iz);
  end
 end 
% 
 end  
 close(vid)