clear
%条件設定
lx=0.1;%x方向の膜の大きさ(m)
ly=0.1;%y方向の膜の大きさ(m)
lt=0.01;%計算時間(s)
dx=1e-3;%x方向の距離刻み(m)
dy=1e-3;%y方向の距離刻み(m)
dt=1e-5;%時間刻み(s)
c=70;%音速(m/s)

%条件から定数等を決定
nx=round(lx/dx)+1;%x方向の配列の数
ny=round(ly/dy)+1;%y方向の配列の数
nt=round(lt/dt)+1;%t方向の配列の数
Kx=(c*dt/dx)^2;%差分計算の時に用いる係数
Ky=(c*dt/dy)^2;%差分計算の時に用いる係数
p=zeros(nx,ny,nt);%音圧値p(x,y,t)
p_p=zeros(nx,ny);%過去の値
n_p=zeros(nx,ny);%現在の値
f_p=zeros(nx,ny);%未来の値
x_1=1:nx-2; %x-1のベクトルを生成
x0=2:nx-1; %xの（以下同文）
x1=3:nx; %x+1
y_1=1:ny-2;%y-1
y0=2:ny-1;%y
y1=3:ny;%y+1

%%%%%%%%%%% Initial condition %%%%%%%%%%%%%%     
[xn yn]=meshgrid(-(nx-1)/2:(nx-1)/2,-(ny-1)/2:(ny-1)/2);
w=nx/20;%ガウス分布のパルス幅（分散）                           
shift=0;%ガウス分布の中心(平均値)
p(:,:,1)=exp(-((xn.^2+yn.^2).^0.5 - shift).^2/(2*w^2));%ガウス分布型の初期値を与える
p(:,:,2)=p(:,:,1);%%初期条件から静かに離す
p_p=p(:,:,1);
n_p=p(:,:,2);

%%%%%%%%%%%%%% FDTD %%%%%%%%%%%%%%%
for t=2:nt-1          
    f_p(x0,y0)=Kx*(n_p(x1,y0)+n_p(x_1,y0)-2*n_p(x0,y0))+Ky*(n_p(x0,y1)+n_p(x0,y_1)-2*n_p(x0,y0))-p_p(x0,y0)+2*n_p(x0,y0);
    p_p=n_p;
    n_p=f_p;
    p(:,:,t+1)=f_p;
end

%%%%%%%% Animation %%%%%%%%%%%%%%%%%
for t=1:nt
    mesh(p(:,:,t));
    axis([1 nx 1 ny -1 1]);
    M(t)=getframe;
end
