Strange/Incorrect values returned by get_phi2()

[TimoEichner]

According to its documentation, the utils.laser_utils.get_phi2() function should return the group delay dispersion (GDD) of the pulse. However, this does not seem to be the case. See for example the following code:

from lasy.laser import Laser
from lasy.profiles import GaussianProfile
from lasy.optical_elements import PolynomialSpectralPhase
from lasy.utils.laser_utils import get_phi2

w0 = 1e-3
tau = 30e-15/np.sqrt(2*np.log(2))
energy = 1
wavelength = 0.8e-6
pol = (1,0)
t_peak = 0

profile = GaussianProfile(wavelength, pol, laser_energy=energy, w0=w0, tau=tau, t_peak=t_peak)

dim='xyt'
lo = (-5e-3, -5e-3, -300e-15)
hi = (5e-3, 5e-3, 300e-15)
npoints = (100, 100, 200)

laser = Laser(dim=dim, lo=lo, hi=hi, npoints=npoints, profile=profile)

gdd_value = 1000e-30

gdd = PolynomialSpectralPhase(laser.profile.omega0, gdd=gdd_value)

laser.apply_optics(gdd)

phi2_output = get_phi2(grid=laser.grid, dim=laser.dim)

print('Calculated:', phi2_output, ', Expected:', gdd_value)

which gives the output: Calculated: (-9.046750219782097e+26, (-5.526846523370059e-28+3.3601099512172264e-27j)) , Expected: 1e-27.

Looking into the function, it looks like the first of the returned values is the intensity-weighted mean of the second derivative of the temporal phase, rather than the second derivative of the spectral phase (which would be the GDD).
What the second value is, I haven't yet understood, but it doesn't seem right to me that it is a complex value while the documentation says it should be the GDD in units of second^2.

I would suggest to either adapt the documentation to more accurately describe what the function does, or adapt the function to return the correct group delay dispersion.

[huixingjian]

Thanks for reporting this bug. Previously get_phi2 function was tested and designed based on a longitudinally Gaussian-shaped laser envelope. PolynomialSpectralPhase seems to distort it drastically. This is fixed and more detailed explain in PR #381 . You can try to use it as following before the PR get merged.

PHI2, BUG_phi2_to_be_fixed = get_phi2(grid=laser.grid, dim=laser.dim)
GOOD_phi2 = np.max(np.roots([4 * np.abs(PHI2), -4, tau**4 *np.abs(PHI2)]))
#tau is the original duration before PolynomialSpectralPhase