paranoya / population-synthesis-toolkit

Spetral synthesis, photometry, SED fitting, etc.
BSD 3-Clause "New" or "Revised" License
3 stars 1 forks source link

SSP Kinematics #20

Open PabloCorcho opened 2 months ago

PabloCorcho commented 2 months ago

SSPs should include a method for convolving the spectra with an input kernel to include kinematic effects.

This is an example that could be adapted to this package

def convolve_ssp(config, los_sigma, los_vel, los_h3=0., los_h4=0.):
    velscale = config["velscale"]
    oversampling = config["oversampling"]
    extra_pixels = config["extra_pixels"]
    ssp_sed = config["ssp_sed"]
    flux = config["flux"]
    # Kinematics
    sigma_pixel = los_sigma / (velscale / oversampling)
    veloffset_pixel = los_vel / (velscale / oversampling)
    x = np.arange(-5 * sigma_pixel, 5 * sigma_pixel) - veloffset_pixel
    losvd_kernel = specBasics.losvd(x, sigma_pixel=sigma_pixel,
                                    h3=los_h3, h4=los_h4)
    sed = fftconvolve(ssp_sed, np.atleast_2d(losvd_kernel), mode="same", axes=1)
    # Rebin model spectra to observed grid
    sed = (
        sed[:, extra_pixels * oversampling : -(extra_pixels * oversampling + 1)]
        .reshape((sed.shape[0], flux.size, oversampling))
        .mean(axis=2)
    )
    ### Mask pixels at the edges with artifacts produced by the convolution
    mask = np.ones_like(flux, dtype=bool)
    mask[: int(5 * sigma_pixel)] = False
    mask[-int(5 * sigma_pixel) :] = False
    return sed, mask

def losvd(vel_pixel, sigma_pixel, h3=0, h4=0):
     y = vel_pixel / sigma_pixel                                                                                                                                                      
     g = (
         np.exp(-(y**2) / 2)
         / sigma_pixel
         / np.sqrt(2 * np.pi)
         * (
             1
             + h3 * (y * (2 * y**2 - 3) / np.sqrt(3))  # H3
             + h4 * ((4 * (y**2 - 3) * y**2 + 3) / np.sqrt(24))  # H4
         )
     )
    return g