integrating a light pulse that has gaussian temporal and spatial characteristic

9 次查看（过去 30 天）
Noah Cowper
Noah Cowper 2021-3-10
回答： Leepakshi 2025-9-2，19:09
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)

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

回答（1 个）

Leepakshi
Leepakshi 2025-9-2，19:09
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

类别

在 Help Center 和 File Exchange 中查找有关 MATLAB 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by