reason of getting zero answer

1 次查看（过去 30 天）
Habib
Habib 2014-11-11
评论： 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
Zoltán Csáti 2014-11-11
Please format your code so that we can read it easier.

请先登录，再进行评论。

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

回答（1 个）

the cyclist
the cyclist 2014-11-11
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更早的评论
the cyclist
the cyclist 2014-11-12
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
Habib
Habib 2014-11-12
thank you for your reply and your kindness
your codes worked, but if we can plot u2, it is understandable that the Trapz worked correctly or not!!!
if it is not bothered you,I have another question!!!
my question is about plot u2, how we can plot u2 by mesh, surf or another plot commands?
(I want plot it as a three-dimensional surface)

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by