求一个MATLAB做杨氏干涉实验的程序，只要程序没错 能出图形都可以！

程序：

clear all

lam=500e-9%输入波长

a=2e-3D=1

ym=5*lam*D/axs=ym%设定光屏的范围

n=101ys=linspace(-ym,ym,n)%把光屏的y方向分成101点

for i=1:n

r1=sqrt((ys(i)-a/2).^2+D^2)

r2=sqrt((ys(i)+a/2).^2+D^2)

phi=2*pi*(r2-r1)/lam

B(i,:)=4*cos(phi/2).^2

end

N=225%确定用灰度等级为255级

Br=(B/4.0)*N%使最大光强对应于最大灰度级（白色）

subplot(1,2,1)

image(xs,ys,Br)%面干涉条纹

colormap(gray(N))

subplot(1,2,2)

plot(B,ys)%画出光强变化曲线

你是两个波函数的叠加。

phi_1=A1 cos(w1*t-(x-x1)/lambda1)和phi_2=A2 cos(w2*t-(x-x2)/lambda2).

当 w1=w2, lambda1=lambda2时发生相干，产生干涉图样。

格点的干涉图样为 abs(phi1+phi2)^2. 我想这个程序matlab很好编，你可以试试看。

P.S. 我这里把两束光简化为两个点光源。如果是两束平行线光源，你需要引入一个y分量。

function photo_diffraction

%光学衍射仿真 矩形

lmda=632.8e-9% 波长

xmax=0.05% 观察屏所取范围

a=1e-3f=1

ymax=xmax

def=0.0001

x=-xmax:def:xmax

y=-ymax:def:ymax

lenm=length(x)

lenn=length(y)

for m=1:lenm

for n=1:lenn

alpha=pi*x(m)*a/(lmda*f)%0.5*k*l*a

beta=pi*y(n)*a/(lmda*f)%bb=0.5*k*w*b

I(m,n)=((sin(alpha))/(alpha))^2*((sin(beta))/(beta))^2

end

end

I=I/(max(max(I)))

[X,Y]=meshgrid(x,y)

figure

mesh(X,Y,I)

xlabel('x')

ylabel('y')

zlabel(' 光强')

rotate3D

hold on

figure

imshow(255*I)

xlabel('x')

ylabel('y')

%---------光学衍射仿真 单缝------------

clear

lam=500e-9

a= 1e-3f=1

xm= 3*lam*f/a

nx= 51

xs=linspace(-xm,xm,nx)

np=51

xp=linspace(0,a,np)

for i=1:nx

sinphi= xs(i)/f

alpha=2*pi*xp*sinphi/lam

sumcos=sum(cos(alpha))

sumsin=sum(sin(alpha))

B(i,:)=(sumcos^2+sumsin^2)/np^2

end

N=255

Br=(B/max(B))*N

figure

subplot(1,2,1)

image(xm,xs,Br)

colormap(gray(N))

subplot(1,2,2)

plot(B,xs)

%--------光学衍射仿真 多缝--------------

clear

lam=500e-9N=2

a= 2e-4z=5d=5*a

xm=2*lam*z/ay0=xm

n=1001

x0=linspace(-xm,xm,n)

for i= 1: n

sinphi=x0(i)/z

alpha=pi*a*sinphi/lam

beta=pi*d*sinphi/lam

B(i,:)=(sin(alpha)./alpha).^2.*(sin(N*beta)./sin(beta)).^2

B1=B/max(B)

end

NC=255

Br=(B/max(B))*NC

figure

subplot(1,2,1)

image(y0,x0,Br)

colormap( gray(NC) )

subplot(1,2,2)

plot(B1,x0)

%--------光学衍射仿真 圆孔--------------

clear

lam=500% '请输入光的波长:')

lam=lam*1e-9

a=2e-3

f=1

m=200

ym=2000*lam*f

ys=linspace(-ym,ym,m)

xs=ys

n=255

for i=1:m

r=xs(i)^2+ys.^2

sinth=sqrt(r./(r+f^2))

x=2*pi*a*sinth./lam

hh=(2*BESSELJ(1,x)).^2./x.^2

b(:,i)=(hh)'.*5000

end

figure

subplot(1,2,1)

image(xs,ys,b)

colormap(gray(n))

subplot(1,2,2)

b(:,m/2)

plot(ys,b(:,m/2))

---