%圓孔衍射的數(shù)值模擬
function diffraction_circle()
distance=input('請輸入光波傳輸距離（毫米）:');%=100000;%觀察屏到衍射屏的距離，單位mm
length=20; %衍射屏寬度，單位mm
lambda=0.0006328;%波長，單位mm

%input('衍射屏的長度(mm)','length')
%input('光波傳播的距離(mm)','distance')
%input('波長(mm)','lambda')
%規(guī)定
%0代表物平面的參數(shù)
%1代表菲涅耳解析DFFT
%2代表菲涅耳數(shù)值DFFT
%3代表瑞利-索末非衍射DFFT
%4代表角譜理論DFFT
%5代表菲涅耳衍射強(qiáng)度SFFT

N0=N_object();%求解物函數(shù)的采樣點(diǎn)的函數(shù)
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;只有當(dāng)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);
%物函數(shù)的抽樣點(diǎn)數(shù)
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
%物函數(shù)
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
%數(shù)值傳遞函數(shù)的菲涅耳衍射積分公式的數(shù)值模擬
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));
%解析傳遞函數(shù)的菲涅耳衍射積分公式的數(shù)值模擬
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);
%數(shù)值傳遞函數(shù)形式的瑞利—索末非衍射積分公式的數(shù)值模擬，它的公式與基爾霍夫衍射公式很接近，省略了基爾霍夫衍射公式的模擬
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));
%角譜理論的頻域衍射積分公式的數(shù)值模擬
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);
%一次傅里葉變換的菲涅耳衍射積分的數(shù)值模擬
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
%一次傅里葉變換菲涅耳衍射積分的數(shù)值模擬，忽略了積分式的二次位相因子
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));