reason of getting zero answer

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

[Zoltán Csáti]

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

[the cyclist]

Did it for him.

Answer 1

[the cyclist]

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.

[Habib]

thank you for your attention. why both K and u1 become scalar? so, how I can do it that u1 and k don't be scalar and i get right answer?

[the cyclist]

For example here:

u1=w0./w1.*exp((-1./w1.^2+1i.*pi./wvl.*R1).*(r1(nn).^2+r2(mm).^2)-1i.*k.*z1);

you are select one value of r1 and one value of r2, and defining the scalar u1 from them.

Did you perhaps mean to write

u1(nn,mm)=w0./w1.*exp((-1./w1.^2+1i.*pi./wvl.*R1).*(r1(nn).^2+r2(mm).^2)-1i.*k.*z1);

so that u1 is also a matrix of values?

I am not certain what you are trying to integrate, if you have two-dimensional matrices.

[Habib]

thank you again,

I tried to integrate respect to r1(nn) or r1. but when program run, the below error appears:

??? Error using ==> trapz at 59 LENGTH(X) must equal the length of the first non-singleton dimension of Y.

Error in ==> gaussianpropagation at 33 u2(nn,mm)=trapz(r1(nn),K.*u1);

whats means error?

how i can eliminate this error?

[the cyclist]

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]

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)