Question

大家谁有pudn下载账号啊，帮帮忙啊！

http://www.pudn.com/downloads579/sourcecode/math/detail2374231.html
http://www.pudn.com/downloads579/sourcecode/math/detail2374230.html

Answer 1

% Script file amplitude_ps.m
%
% for students
% all kinds of distrabutiona of photon sieve figure
% 孔密度:选用不同窗函数，半窗
% d/w=1.53
clear
close all                                               %close all open figure windows

%lambda=632.8e-9;                                          %wave length=6328nm
lambda=400e-9;  
p=inf;                                                   %distance from source to ps plane
q=1;
f=1/(1/p+1/q);                                          %focal length
An=1;                                                   %real amplitude
k=2*pi/lambda;                                          %wave number
d_min=100e-6;                                           %smallest dtructure size
an_min=d_min./2;                                        %radii of smallest hole
ratio=1.53;                                              %prearrange d/w
 %ratio=3.51;
%ratio=5.52;
%ratio=7.53;
variable=lambda*f;
w_min=d_min/ratio;
% rn_max=9.6e-04;
% N=floor(rn_max^2/variable);
N=floor((w_min.^2/variable-2+variable/w_min.^2)/4);
if rem(N,2)==1
    N=N-1;
end
num1=2:2:N;                                             %numbers of zones
rho1=sqrt(num1*variable);                               %radii of transparence zones
rho2=sqrt((num1-1)*variable);                           %radii of opacifying zones
w=rho1-rho2;                                            %width of zones
rn=rho2+w/2;                                            %center of pinholes
length_rn=length(rn);
an=ratio.*w/2;                                          %radii of hole
U=0;                                                    %inatial complexx amplitude
gridnumber=200;                                         %focal place gridnumber*guidnumber
XYlimit=0.2e-3;                                         %observation slope in screen
nX=linspace(-XYlimit,XYlimit,gridnumber);               %coodinate in the XY plane
nY=nX;                                                  %coodinate in the XY plane
%======figure of window function===========================================
double_length_rn=2*length_rn;
%win=rectwin(double_length_rn);                        %Tectangular Window
%win=Triang(double_length_rn);                         %Triangle Window
%win=Gausswin(double_length_rn,3);                                %Gauss Window
 %win=hamming(double_length_rn);                                %Hamming Window
   %win=hann(double_length_rn+2);win=win(2:end-1);                           %Hanning Window
   win=chebwin(double_length_rn+2,60);win=win(2:end-1);                              %Chebv Window
    %win=blackman(double_length_rn+2);win=win(2:end-1);                                %Blackman Window
% win=tukeywin(double_length_rn,0.5);                                %Window function
win=win(length_rn+1:end);
figure(1)
plot(win)
title('Window Function')
%======dentity of holes====================================================
angle_min=2.*asin(an./rn);                              %minmum angle between two holes
Nf=pi.*an.^2./(lambda*q);                               %Fresnel number
M=floor(2*pi./angle_min.*win');
% density_zone=M./(2*pi);
holenumber=sum(M);                                      %total number of pinholes

for m=1:length_rn                                       %mth layer zone
    alfa=2*pi/M(m);
    theta=zeros(1,M(m));
    theta(1)=unifrnd(0,2*pi-angle_min(1));              %interval of per two centers of a hole
    for miu=1:M(m)
        if miu>=2
            theta(miu)=theta(miu-1)+delta;
            if  miu<M(m)
                alfa=(2*pi-(theta(miu)-theta(1))-2*angle_min(m))/(M(m)-miu);
            else
                alfa=2*pi-(theta(miu)-theta(1))-4*angle_min(m);
            end
        end
        if angle_min(m)<alfa
            delta=unifrnd(angle_min(m),alfa);
        else
            delta=angle_min(m);
        end
        xn=rn(m).*cos(theta(miu));
        yn=rn(m).*sin(theta(miu));
        [xout,yout]=scircle1(xn,yn,an(m));
       % figure(3)
       % fill(xout,yout,'w')
        %hold on
        %         ========calculation of diffraction field Un(X,Y)==================
        Ln=(xn.^2+yn.^2)./(2*p);                        %eikonal
        gn=xn/p;
        hn=yn/p;
        X1=nX-xn;                                      %1_D X'
        Y1=nY-yn;                                       %3-D Y'
        R1=sqrt(X1.^2+Y1.^2);                           %R'
        X2=X1-gn*q;                                     %X''
        Y2=Y1-hn*q;                                     %Y''
        rou=sqrt(X2.^2+Y2.^2);
        t=k*an(m).*rou/q;
        Jinc=(besselj(0,t)+besselj(2,t))/2;
        F0=Nf(m).*Jinc;
        F1=Nf(m).^2*(3*besselj(0,t)+2*besselj(2,t)-besselj(4,t))/12;
        F=F0+j*F1;
        UnXY=2*An.*exp(j*k*(Ln+R1.^2/(2*q))).*F;
        if isfinite(UnXY)==0
            disp('Warning: the value is NaN or inf or -inf')
            break
        end
        U=UnXY+U;
    end
    disp('N=')
    disp(m)
end
axis off
axis square
set(gcf,'Units','centimeters','position',[9 5 10 10],'color','k')
title('Photon sieve','color','w')
%=======distribution of intensity==========================================
I=abs(U).^2;
Imax=max(I);
I_guiyi=I./Imax;
Airy_center_ps=nX(find(I_guiyi==max(I_guiyi)));                  %position of maxmum, peak
figure(4)
plot(nX.*1e3,I_guiyi)
xlabel('x/mm'),ylabel('Intensity')
figure(5)
semilogy(nX.*1e3,I_guiyi)             %halp_log axis
xlabel('x/mm'),ylabel('Intensity')
grid on

直接给你贴出来了

追问

谢谢！一共两个程序的，可以帮我贴下另外一个吗？

追答

太长贴不下，我再试试