I am currently writing a small Python script to estimate the pulse duration from the optical spectrum.
At the end, the idea is to observe the effects of the spectral phase on the pulse duration and shape. But for now, I cannot retrieve a well-known result without taking into account spectral phase.
Indeed, a time-bandwidth limited Gaussian pulse with an optical spectrum of 30 nm width centered at 800 nm has a duration of around 31 fs. This is a basic result for a rule of thumb. You can check the result here and here, for instance.
Sadly, my Python script, using
numpy/matplotlib, gives me a pulse duration which is twice too large:
from numpy import power, linspace, sqrt, log, exp, abs, max, conj from numpy.fft import fft, fftshift, fftfreq from scipy.constants import c from matplotlib import pyplot as plt # Number of points for the computations n_points = 2**9 # Spectral width and center in wavelength (nm) wl_center = 800e-9 wl_width = 30e-9 # Conversion of spectral parameters in frequency (s-1) freq_center = c/wl_center freq_width = wl_width*c/power(wl_center, 2) # Building of the x axes frequencies = linspace(freq_center - freq_width*8, freq_center + freq_width*8, n_points) wavelengths = c/frequencies # Building the Gaussian spectrum according to the parameters diff_freq = frequencies - freq_center sig_freq = freq_width/(sqrt(8*log(2))) spectrum_freq_abs = exp(-(power(diff_freq, 2))/(2*power(sig_freq, 2))) # FFT step = frequencies - frequencies time = fftshift(fftfreq(n_points, d=step)) field = fftshift(fft(spectrum_freq_abs)) # Plotting # The commented line plot the spectrum vs the wavelength. The second line plot the pulse. # The second line plot the optical pulse envelope vs the time in fs. It is twice too large (60 fs fwhm instead of 30 fs.) # plt.plot(wavelengths*1e9, spectrum_freq_abs) plt.plot(time*1e15, abs(field*conj(field))) plt.show()
I checked the time array created thanks to the
fftfreq function of
numpy. It seems to be good according to the frequency array size and step. So I don't think it's an x axis problem.
I cannot find where I am wrong in this small and easy Python script, which is very frustrating. Maybe the community could give me an hint somewhere?
PS: This is a repost from physics community, where I've been kindly redirected to scicomp.