diffraction_circle.m

用matlab模拟圆孔衍射图像

%圆孔衍射的数值模拟
function diffraction_circle()
distance=input('请输入光波传输距离（毫米）:');%=100000;%观察屏到衍射屏的距离，单位mm
length=20; %衍射屏宽度，单位mm
lambda=0.0006328;%波长，单位mm

%input('衍射屏的长度(mm)','length')
%input('光波传播的距离(mm)','distance')
%input('波长(mm)','lambda')
%规定
%0代表物平面的参数
%1代表菲涅耳解析DFFT
%2代表菲涅耳数值DFFT
%3代表瑞利-索末非衍射DFFT
%4代表角谱理论DFFT
%5代表菲涅耳衍射强度SFFT

N0=N_object();%求解物函数的采样点的函数
N1=length^2/(lambda*distance);%N<=
N2=length^2/(lambda*distance);%N>=
N3=length^2/(lambda*sqrt(distance^2+length^2/2));%N>=

N4=length^2/(lambda*sqrt(distance^2+length^2/2));%N<=
%N3=N4
N5=length^2/(lambda*distance);%N>=
%N1=N2=N5
choice=0;
if  N0<=N4
    choice=4;%一般情况下，满足N0<N4,都有N0<N1;只有当distance短时才会有差别
    N=N0;
else
   choice=3;
   N=N0;
end
if distance>150000
    choice=5;
    N=N0;
end
%判断N属于那个2的几次幂段   
i=0;
while N>=2
    N=N/2;
    i=i+1;
end
N=2^(i+1);   

switch choice
    case 1
        U=Frensel_DFFT(length,distance,lambda,N);
    case 2
        U=Frensel_DFFT_S(length,distance,lambda,N);
    case 3
        U=Rayleigh_DFFT(length,distance,lambda,N);
    case 4
        U=Angular_DFFT(length,distance,lambda,N);
    case 5
        U=Frensel_SFFT_I(length,distance,lambda,N);%不使用Frensel_SFFT
    otherwise
        U=Frensel_DFFT_S(length,distance,lambda,N);
end
U=abs(U);
U=U.^2;
%归一化
U=U/max(max(U));
imshow(U);
%物函数的抽样点数
function N=N_object()
for i=5:10
    N=2^i;
    U=Object(N);
    U=fftshift(fft2(U));
    U=abs(U);
    eps=U(N,N)/max(max(U));
    if eps<0.0002
        break;
    end
end
%物函数
function U0=Object(N)
U0=zeros(N,N);
for x=1:N
    for y=1:N
        if((x-N/2)^2+(y-N/2)^2)<=(N/4)^2
            U0(x,y)=1;
        end
    end
end
%数值传递函数的菲涅耳衍射积分公式的数值模拟
function U=Frensel_DFFT_S(length,distance,lambda,N)
U=zeros(N,N);
dl=1/length;

k=2*pi/lambda;
temp=exp(i*k*distance)/(i*lambda*distance);
for m=1:N
    for n=1:N
        U(m,n)=temp*exp(i*k/2/distance*dl^2*((m-N/2)^2+(n-N/2)^2));
    end
end
U=fftshift(fft2(U));
U=U.*fftshift(fft2(Object(N)));
U=fftshift(ifft2(U));
%解析传递函数的菲涅耳衍射积分公式的数值模拟
function U=Frensel_DFFT(length,distance,lambda,N)
U=Object(N);
U=fftshift(fft2(U));
deltafx=1/length;
deltafy=1/length;
k=2*pi/lambda;
for m=1:N
    for n=1:N
        U(m,n)=U(m,n)*exp(i*k*distance*(1-lambda^2/2*((deltafx^2*(m-N/2)^2+deltafy^2*(n-n/2)^2))));
    end
end
U=ifft2(U);
%数值传递函数形式的瑞利—索末非衍射积分公式的数值模拟，它的公式与基尔霍夫衍射公式很接近，省略了基尔霍夫衍射公式的模拟
function U=Rayleigh_DFFT(length,distance,lambda,N)
U=Object(N);
U=fftshift(fft2(U));
h=zeros(N,N);
k=2*pi/lambda;
dx=length/N;
dy=length/N;
for m=1:N
    for n=1:N
        temp=distance^2+(dx*(m-N/2))^2+(dy*(n-N/2))^2;
        h(m,n)=distance*exp(i*k*(temp)^0.5)/(i*lambda*temp);
    end
end
h=fftshift(fft2(h));
U=U.*h;
U=fftshift(ifft2(U));
%角谱理论的频域衍射积分公式的数值模拟
function U=Angular_DFFT(length,distance,lambda,N)
U=Object(N);
U=fftshift(fft2(U));
k=2*pi/lambda;
deltafx=1/length;
deltafy=1/length;
for m=1:N
    for n=1:N
        U(m,n)=U(m,n)*exp(i*k*distance*(1-(lambda*deltafx*(m-N/2))^2-(lambda*deltafy*(n-N/2))^2)^0.5);
    end
end
U=ifft2(U);
%一次傅里叶变换的菲涅耳衍射积分的数值模拟
function U=Frensel_SFFT(length,distance,lambda,N)
U=Object(N);
k=2*pi/lambda;
dx0=length/N;
dy0=length/N;
for m=1:N
    for n=1:N
        U(m,n)=U(m,n)*exp(i*k/(2*distance)*((dx0*(m-N/2))^2+(dy0*(n-N/2))^2));
    end
end
U=fftshift(fft2(U));
dx=lambda*distance/length;
dy=lambda*distance/length;
for m=1:N
    for n=1:N
        U(m,n)=U(m,n)*exp(i*k/(2*distance)*((dx*(m-N/2))^2+(dy*(n-N/2))^2));
    end
end
%一次傅里叶变换菲涅耳衍射积分的数值模拟，忽略了积分式的二次位相因子
function U=Frensel_SFFT_I(length,distance,lambda,N)
U=Object(N);
k=2*pi/lambda;
dl=length/N;
temp=i*k/(2*distance)*dl^2;
for m=1:N
    for n=1:N
        U(m,n)=U(m,n)*exp(temp*((m-N/2)^2+(n-N/2)^2));
    end
end
U=fftshift(fft2(U));