我想画(R,x2,z2)的图,但是我现在只会画(R,x2)的图,我想把x2,z2都看成变量,R=integral(h1(x1,xs,z2)*h2(x2,xs)*f(xs)),其中x1=0
111.PNG (32.43 KB, 下载次数: 0)
2020-10-26 09:50 上传
微信图片_20201026091523.jpg (30.59 KB, 下载次数: 0)
2020-10-26 09:51 上传
lamda=800*10^(-6);
k=2.*pi./lamda;
b=180*10^(-3);
d=400*10^(-3);
zi=400;
zs1=400;
Zb=zs2*(zs1+zi)/(zi+zs1+zs2);
w=1.8;
N=100;
n=-100:1:100;
deltxs=0.01;
deltx0=0.002;
x0=-1:deltx0:1;
xs=-1:deltxs:1;
x2=-2:0.02:2;
T=@(x0,n) sum(0*(abs(x0-n.*d)/b>0.5)+1*(abs(x0-n.*d)/b<=0.5))
h1=@(x1,xs) k/(2*pi*1i*sqrt(zs1*zs2))*exp(1i*k*(zs1+zs2))*sum(arrayfun(@(x0)T(x0,n),x0).*exp(1i.*k./(2.*zs1).*(x0-xs).^2+1i.*k./(2.*zs2).*(x1-x0).^2)*deltx0)
x1=0
h1_mat=arrayfun(@(xs) h1(x1,xs),xs);
[xsm,x2m]=meshgrid(xs,x2);
h2=@(x2,xs) sqrt(k/(2*pi*zi))*exp(-1i*pi/4+1i*k*zi+1i*k/(2*zi)*(x2-xs).^2);
f=@(xs) exp(-2.*xs.^2./(w.^2));
f_mat=arrayfun(@(xs) f(xs),xsm);
h2_mat=arrayfun(@(x2,xs) h2(x2,xs),x2m,xsm)
G=sum(h2_mat.*f_mat.*repmat(h1_mat,size(h2_mat,1),1)*deltxs,2);
G=abs(G').^2;
G_max=max(G);
plot(x2,G1)