reason of getting zero answer

Question

Habib 2014-11-11

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/162236-reason-of-getting-zero-answer

评论： Habib 2014-11-12

I write below code for solve diffraction integral by trapz order but in final I have got zero answer for u2. what is problem? please help to find defect.

%gaussian propagation%
clear
clc
r_1=linspace(-1,1,50);
r_2=linspace(-1,1,50);
[r1,r2]=meshgrid(linspace(-1,1,100));
z2=4.2;                % is the axial distance from the beam's narrowest point
z1=1;
L=z2-z1;
wvl=500;                   % is wavelength
k=2*pi./wvl;
w0=sqrt(wvl.*z1./pi);                           % is the waist size
w1=w0.*sqrt(1+wvl.*z1./pi.*w0.^2);         % is the radius at which the field amplitude and intensity drop to 1/e and 1/e2 of...
                                                         %...their axial values, respectively,
w2=w0.*sqrt(1+wvl.*z2./pi.*w0.^2);                                                      
R1=z1.*(1+(pi.*w0.^2./wvl.*z1).^2);            % is the radius of curvature of the beam's wavefronts
R2=z2.*(1+(pi.*w0.^2./wvl.*z2).^2); 
u1=w0./w1.*exp((-1./w1.^2+1i.*pi./wvl.*R1).*(r1.^2+r2.^2)-1i.*k.*z1); % mathematical form of Gaussian beam
figure(1)
mesh(r1,r2,real(u1))
K=zeros(1,length(r1));
u1=zeros(1,length(r1));
u2=zeros(1,length(r1));
for nn=1:50
    for mm=1:50
        r1=linspace(0,1,length(r1));
         r2=linspace(0,1,length(r1));
    K=exp(-1i.*pi.*(r1(nn).^2+z1.^2+r2(mm).^2+z1.^2-2.*(r1(nn).*r2(mm)+z1.*z2))./(wvl.*L)-1i.*k.*L);
    u1=w0./w1.*exp((-1./w1.^2+1i.*pi./wvl.*R1).*(r1(nn).^2+r2(mm).^2)-1i.*k.*z1);
    u2(nn,mm)=trapz(K.*u1);
    end
end

2 个评论
显示无隐藏无

Zoltán Csáti 2014-11-11

Please format your code so that we can read it easier.

the cyclist 2014-11-11

Did it for him.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

the cyclist 2014-11-11

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/162236-reason-of-getting-zero-answer#answer_158506

在 MATLAB Online 中打开

I have not tried to fully understand your code, but the specific reason that u2 is all zeros is that when you calculate

u2(nn,mm)=trapz(K.*u1)

both K and u1 are scalars, and the numerically evaluated integral of a scalar is zero.

5 个评论
显示 3更早的评论隐藏 3更早的评论

the cyclist 2014-11-12

在 MATLAB Online 中打开

The following will run to completion. It seems that it might do what you intend, but I am not certain of that.

clear
clc
r_1=linspace(-1,1,50);
r_2=linspace(-1,1,50);
[r1,r2]=meshgrid(linspace(-1,1,100));
z2=4.2;                % is the axial distance from the beam's narrowest point
z1=1;
L=z2-z1;
wvl=500;                   % is wavelength
k=2*pi./wvl;
w0=sqrt(wvl.*z1./pi);                           % is the waist size
w1=w0.*sqrt(1+wvl.*z1./pi.*w0.^2);         % is the radius at which the field amplitude and intensity drop to 1/e and 1/e2 of...
                                                         %...their axial values, respectively,
w2=w0.*sqrt(1+wvl.*z2./pi.*w0.^2);                                                      
R1=z1.*(1+(pi.*w0.^2./wvl.*z1).^2);            % is the radius of curvature of the beam's wavefronts
R2=z2.*(1+(pi.*w0.^2./wvl.*z2).^2); 
u1=w0./w1.*exp((-1./w1.^2+1i.*pi./wvl.*R1).*(r1.^2+r2.^2)-1i.*k.*z1); % mathematical form of Gaussian beam
figure(1)
mesh(r1,r2,real(u1))
K=zeros(1,length(r1));
u1=zeros(1,length(r1));
u2=zeros(1,length(r1));
r1=linspace(0,1,length(r1));
r2=linspace(0,1,length(r1));
for nn=1:50
    for mm=1:50
    K(nn,mm)=exp(-1i.*pi.*(r1(nn).^2+z1.^2+r2(mm).^2+z1.^2-2.*(r1(nn).*r2(mm)+z1.*z2))./(wvl.*L)-1i.*k.*L);
    u1(nn,mm)=w0./w1.*exp((-1./w1.^2+1i.*pi./wvl.*R1).*(r1(nn).^2+r2(mm).^2)-1i.*k.*z1);
    end
    u2(nn)=trapz(r1,K(nn,:).*u1(nn,:));
end