大家谁有pudn下载账号啊,帮帮忙啊!
发布网友
发布时间:2022-07-26 22:50
我来回答
共1个回答
热心网友
时间:2023-10-23 04:30
% 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
直接给你贴出来了
追问谢谢!一共两个程序的,可以帮我贴下另外一个吗?
追答太长贴不下,我再试试