据所学内容编写一个程序，演示多个振动源(>2 个)的水波干涉现象

%水波干涉的动画演示

a=4;w=4;A0=0.3;v=1;%水波源的位置及方程各参量

[x,y]=meshgrid([-4*pi:pi/20:4*pi],[-4*pi:pi/20:4*pi]);%产生网络坐标

r1=sqrt((x-a).^2+y.^2)+5;%空间点（x，y）到波源1（a，0）的距离

c1=1./r1;%波源1传播到各点的振幅

r2=sqrt((x+a).^2+y.^2)+5;%空间点（x，y）到波源2（-a，0）的距离

c2=1./r2;%波源2传播到各点的振幅

t=0:0.1:10;%设定运行时间

for j=1:41%此循环产生动画

    z=A0*c1.*cos(w*(t(j)-r1/v))+0.3*c2.*cos(w*(t(j)-r2/v));%波动叠加

surf(x,y,z)%绘制波动叠加曲面

colormap(gray);%创建线性灰阶色图给图形表面着色

axis equal;%使每个坐标轴的刻度增量相同

shading interp;%用细致添加描影法渲染图形

view([180 90])%在方位角180°，仰角90°观察图形

axis([-15 15 -15 15 -0.1 0.1])%设定x、y、z、轴的范围

axis off%去掉坐标轴

m(j)=getframe;%保存绘制好的图形到变量数组m（j）

end

movie(m,10)%放映储存在m矩阵中电影动画10遍

%水波干涉的动画演示

a=4;w=4;A0=0.3;v=1;%水波源的位置及方程各参量

[x,y]=meshgrid([-4*pi:pi/20:4*pi],[-4*pi:pi/20:4*pi]);%产生网络坐标

r1=sqrt((x-a).^2+(y-a).^2)+5;%空间点（x，y）到波源1（a，a）的距离

c1=1./r1;%波源1传播到各点的振幅

r2=sqrt((x-a).^2+(y+a).^2)+5;%空间点（x，y）到波源2（a，-a）的距离

c2=1./r2;%波源2传播到各点的振幅

r3=sqrt((x+a).^2+(y+a).^2)+5;%空间点（x，y）到波源3（-a，-a）的距离

c3=1./r3;%波源3传播到各点的振幅

r4=sqrt((x+a).^2+(y-a).^2)+5;%空间点（x，y）到波源3（-a，a）的距离

c4=1./r4;%波源4传播到各点的振幅

t=0:0.1:10;%设定运行时间

for j=1:41%此循环产生动画

    z=A0*c1.*cos(w*(t(j)-r1/v))+0.3*c2.*cos(w*(t(j)-r2/v))+0.3*c3.*cos(w*(t(j)-r3/v))+0.3*c4.*cos(w*(t(j)-r4/v));%波动叠加

surf(x,y,z)%绘制波动叠加曲面

colormap(gray);%创建线性灰阶色图给图形表面着色

axis equal;%使每个坐标轴的刻度增量相同

shading interp;%用细致添加描影法渲染图形

view([180 90])%在方位角180°，仰角90°观察图形

axis([-15 15 -15 15 -15 15])%设定x、y、z、轴的范围

axis off%去掉坐标轴

m(j)=getframe;%保存绘制好的图形到变量数组m（j）

end

movie(m,10)%放映储存在m矩阵中电影动画10遍