求一个用MATLAB做圆孔或圆屏衍射的实验，要求程序正确，能出图像，非常感谢

%%傅里叶变换法

clear all

z1=20

r=4

N=512

d=180

lambda=0.0039

k=2*pi/lambda

f1=16*pi/粗芹lambda

R=z1+f1^2/z1

k0=-N/2:N/2-1

a=sqrt(lambda*d/岩拍毕N)

[x0,y0]=meshgrid(k0*a)

w=4*sqrt(1+(z1/f1)^2)

Ein=100/w*exp(-(x0.^2+y0.^2)/w^2).*exp(-i*(k*(z1+(x0.^2+y0.^2)/2/R)-atan(z1/f1)))

for m=1:N%%%%%%%%%%%%%%%%%%%圆孔

for n=1:N

if (x0(m,n)^2+y0(m,n)^2>(r)^2)

Ein(m,n) = 0

end

end

end

E=fft2(Ein.*exp(i*k/2/d*(x0.^2+y0.^2)))%%%%傅里叶变换

Eout=fftshift(abs(E))

x=x0

y=y0

U=exp(i*k*d)/i/lambda/d.*exp(i*k/贺枣2/d*(x.^2+y.^2)).*Eout

I=U.*conj(U)

figure,mesh(x,y,I)

figure,plot(x,I)

先写这么多，待会儿再补充，吃饭去了，回来接着写。。

function [uf,vf] = sfn(u,v)%此函数为自定义的sfn函数，u,v为输入变量

uf、vf为输出变量

global e b k%定义全局变量 e b k；

ny = size(u,1)%求出矩阵u的行逗铅瞎的数目，赋给ny

nx = size(u,2)%求出矩阵u的列的数目，赋给nx；

uer = [u(:,1),u,u(:,nx)]

%构造一个矩阵uer，u(:,1）为矩阵u的第一列； u(:,nx) 为矩阵第最后一列,u就不解释了。

uec = [u(1,:)uu(ny,:)]

%构造一个矩阵uec，u(1,:)为矩阵u的第一行；u(ny,:)为u最后一行

ul = uec(3:ny+2,:)+uec(1:ny,:)+uer(:,1:nx)+uer(:,3:nx+2)-4*u

%用矩阵uec的第三行到第（ny+2）行+uec的第一行到第ny行+uer的第一列到第nx列+uer的第三行到（nx+2）列-4*u

u3 = u.*u.*u%把u的三次方赋给u3

uf = (u-u3/3-v)/e + k^2*ul%这行不用激罩解释吧，简单的矩阵运算。

vf = e*(u+b-0.5*v)%同上。。山空

程序2：

clc %清屏

global k b e %定义全局变量k b e

nsteps = 2000%将2000赋给nsteps

nxy = 100%100赋给nxy

u = 5*ones(nxy,1)*([1:nxy]-0.4*nxy)/nxy

%ones(nxy,1)是创建一个nxy*1的矩阵，[1:nxy]是行向量[1 2 3 4 .. nxy]

v = 5*(([1:nxy]-0.4*nxy)')*ones(1,nxy)/nxy%同上

h = 0.04

e = 0.1

b = 0.67

colormap hsv% 创建一个hsv颜色映射表

set(0,'DefaultSurfaceEdgecolor','none')%设置为缺省边缘颜色为none。

k=1

for i=1:nsteps

[uf,vf] = sfn(u,v)%调用上面写的sfn函数

u = u+h*uf%矩阵运算

v = v+h*vf

if mod(i,3) == 1 %如果i对3求余数等于1

pcolor(u)drawnow%按u中元素着色，画图。。

end

end

---