integrating a light pulse that has gaussian temporal and spatial characteristic

I am trying to solve this integration of a gaussian pulse, it is a integration over all wave numbers i.e. from -inf to inf. I have tried it with the infinite bounds and also with definate bounds, and consistently get an infinite result. The goal is to have a pulse wich is E_ref and compare it to E_g, but I can't get a correct result for E_g any help would be greatly appreciated. (Also this is a replication of work from a paper GROUP VELOCITY OF GAUSSIAN BEAMS M.V. VASNETSOV, V.A. PAS’KO, where you can see the equations I have used). The following is my code so far.

% estimation of the group velocity as compared to a reference pulse

%for a gaussian pulse

clear, clc

%parameters of the wave

T=10e-15; %duration of pulse

K=1/(10e-6); %wave# of carrier frequency

c=3e8; %speed of light

w=100e-6; %initial beam waist

zr=50e-3; %Rayleigh length

%z=4.99e-3:1e-8:5.01e-3;

z=linspace(zr-1e-5,zr+1e-5,2000)

n=1;

tr=n*zr/c;

E_0=1;

%k=K-5*(2/(c*T)):K+5*(2/(c*T));

%k=linspace(K-5*(2/(c*T)),K+5*(2/(c*T)),length(z));

%Reference pulse

E_ref=exp(-1.*(z-c*tr).^2/(c^2*T^2)).*exp(i*K.*(z-c*tr));

for i=1:length(z)

F = @(k) [E_0./sqrt(1+(2*z(i)./(k*w^2)).^2)].*exp(-c^2*T^2*[(k-K).^2/4]...

-i*atan(2*z(i)./(k*w^2))+i*k*(z(i)-c*tr));

E_g(i)=integral(F,K-(2/(c*T)),K+(2/(c*T)));

end

E_1=abs(E_ref).^2/max(E_ref);

E_2=abs(E_g).^2/max(E_g);

figure(1)

plot(z,E_1)

figure(2)

plot(z,E_2)

Answer 1

[Leepakshi]

Hi Noah,

The issue with your code is that the integration bounds for the Gaussian in k-space are too narrow, given the very short pulse duration (T), which makes the Gaussian extremely broad. As a result, most of the pulse's spectral energy lies outside your chosen integration range, causing the integral to return an infinite or unreasonably large value. To fix this, widen your integration bounds to at least ±5 standard deviations (sigma) from the central wavenumber (K), where sigma is 2/(c*T), so you capture the entire relevant spectrum and obtain a correct, finite result.

Try replacing

E_g(i)=integral(F,K-(2/(c*T)),K+(2/(c*T))); 

with

E_g(i)=integral(F, K-5*(2/(c*T)), K+5*(2/(c*T)), 'ArrayValued', true, 'AbsTol', 1e-12, 'RelTol', 1e-10); 

Hope this helps!

Thanks