# Copyright (C) 2025 LMU Munich
# Author: Luca Tortorelli
# created: Apr 2025
import os
import numpy as np
import scipy.interpolate
from scipy.integrate import trapezoid
from scipy.ndimage import gaussian_filter1d
import galsbi
from galsbi.ucat.spectrum_util import extinction_coefficient
[docs]
def apply_extinction(spectra, observed_wavelengths, excess_b_v, extinction_spline):
"""
This function applies Milky Way extinction to galaxy spectra.
:param spectra: (array_like[n_gal, n_lambda]) input spectra.
:param observed_wavelengths: (array_like[n_gal, n_lambda]) observed- frame
wavelengths array.
:param excess_b_v: (array_like[n_gal,]) E(B-V) Milky Way extinction.
:param extinction_spline: (obj) scipy InterpolatedUnivariateSpline object.
:return ext_spectra: (array_like[n_gal, n_lambda]) extincted spectra.
"""
observed_wavelengths_micrometer = observed_wavelengths / 10000
ext_spectra = np.array(
[
spectra[i]
* 10
** (
-2
/ 5
* extinction_coefficient(
observed_wavelengths_micrometer[i], excess_b_v[i], extinction_spline
)
)
for i in range(len(spectra))
]
)
return ext_spectra
[docs]
def apply_velocity_broadening(
spectra, lam, sigma_smooth, min_wave=3000, max_wave=20000
):
"""
This function applies velocity broadening to the spectra based on the
velocity dispersion sampled from GalSBI.
:param spectra: (array_like[n_gal, n_lambda]) input spectra.
:param lam: (array_like[n_lambda,]) rest-frame wavelength array.
:param sigma_smooth: (array_like[n_gal,]) Velocity dispersions in km/s.
:param min_wave: (float) Minimum wavelength.
:param max_wave: (float) Maximum wavelength.
:return vel_disp_smooth_spectra: (array_like[n_gal, n_lambda]) Velocity
dispersion smoothed spectra.
"""
def process_spectrum(spec, lambda_arr, minl, maxl, sigma, ckms=299792.458, m=4):
"""Process a single galaxy spectrum with smoothing and interpolation.
:param spec: (array_like[n_lambda,]) input spectrum.
:param lambda_arr: (array_like[n_lambda,]) rest-frame wavelength array.
:param minl: (float) Minimum wavelength.
:param maxl: (float) Maximum wavelength.
:param sigma: (float) Velocity dispersion in km/s.
:param ckms: (float) Speed of light in km/s.
:param m: (int) Multiplier for kernel range (default: 4).
:return spec: (array_like[n_lambda,]) Velocity dispersion smoothed
spectrum.
"""
nspec = len(lambda_arr)
dlstep = (np.log(maxl) - np.log(minl)) / nspec
lnlam = np.arange(1, nspec + 1) * dlstep + np.log(minl)
# Interpolate to logarithmic wavelength space
interp_spec = scipy.interpolate.interp1d(
np.log(lambda_arr), spec, kind="linear", bounds_error=False, fill_value=0
)
tspec = interp_spec(lnlam)
# Gaussian smoothing
fwhm = sigma * 2.35482 / ckms / dlstep
psig = fwhm / 2.0 / np.sqrt(-2.0 * np.log(0.5))
grange = int(np.floor(m * psig))
# Smooth spectrum using Gaussian kernel
tnspec = gaussian_filter1d(tspec, psig, mode="constant", truncate=m)
# Preserve unsmoothed ends
tnspec[:grange] = tspec[:grange]
tnspec[-grange:] = tspec[-grange:]
# Interpolate back to original wavelength space
il = np.searchsorted(lambda_arr, minl)
ih = np.searchsorted(lambda_arr, maxl)
interp_tnspec = scipy.interpolate.interp1d(
np.exp(lnlam), tnspec, kind="linear", bounds_error=False, fill_value=0
)
spec[il:ih] = interp_tnspec(lambda_arr[il:ih])
return spec
def process_all_galaxies(
specs, lambda_arr, minl, maxl, sigma, ckms=299792.458, m=4
):
"""Process spectra for all galaxies.
:param specs: (array_like[n_gal, n_lambda]) input spectra.
:param lambda_arr: (array_like[n_lambda,]) rest-frame wavelength array.
:param minl: (float) Minimum wavelength.
:param maxl: (float) Maximum wavelength.
:param sigma: (array_like[n_gal,]) Velocity dispersions in km/s.
:param ckms: (float) Speed of light in km/s.
:param m: (int) Multiplier for kernel range (default: 4).
:return results: (array_like[n_gal, n_lambda]) Velocity dispersion
smoothed spectra.
"""
num_galaxies = specs.shape[0]
results = np.empty_like(specs)
for i in range(num_galaxies):
results[i, :] = process_spectrum(
specs[i, :], lambda_arr, minl, maxl, sigma[i], ckms, m
)
return results
vel_disp_smooth_spectra = process_all_galaxies(
spectra, lam, min_wave, max_wave, sigma_smooth
)
return vel_disp_smooth_spectra
[docs]
def compute_Temple21_AGN_contributions(
par, rest_frame_fluxes_erg_s_A, rest_frame_wave_A, fagn, redshift
):
"""
This function computes the quasar spectra from the Temple+21 AGN templates.
As specified in the paper and through private communication, we compute the AGN
bolometric luminosity as fraction fagn of the galaxy bolometric luminosity and
apply simple prescription to compute the luminosity at 3000 Å and the absolute
magnitude in the i-band. These quantities are needed for qsogen package to work.
:param par: (obj) par objects containing the Ucat parameters.
:param rest_frame_fluxes_erg_s_A: (array_like[n_gal, n_lambda]) galaxy
rest-frame flux in erg/s/Å from ProSpect.
:param rest_frame_wave_A: (array_like[n_lambda,]) rest-frame wavelength array
in Ångstrom.
:param fagn: (array_like[n_gal,]) fraction of AGN luminosity sampled from
GalSBI-SPS.
:param redshift: (array_like[n_gal,]) galaxy redshift.
:return quasar_fluxes_erg_s_cm2_A: (array_like[n_gal, n_lambda]) quasar
observed-frame flux in erg/s/cm2/Å.
"""
import sys
if os.path.isabs(par.abs_path_temple2021):
sys.path.append(par.abs_path_temple2021)
else:
os.makedirs(os.path.join(galsbi.__path__[0], "resources/qsogen"), exist_ok=True)
os.system(
f'git clone https://github.com/MJTemple/qsogen.git '
f'{os.path.join(galsbi.__path__[0], "resources/qsogen")}'
)
sys.path.append(os.path.join(galsbi.__path__[0], "resources/qsogen"))
from qsosed import Quasar_sed
quasar_fluxes_erg_s_cm2_A = np.zeros_like(rest_frame_fluxes_erg_s_A)
bol_wave_mask = np.where(
(rest_frame_wave_A >= 50000) & (rest_frame_wave_A <= 200000)
)
for i in range(len(rest_frame_fluxes_erg_s_A)):
rest_frame_Lbol_erg_s = trapezoid(
rest_frame_fluxes_erg_s_A[i][bol_wave_mask],
x=rest_frame_wave_A[bol_wave_mask],
)
LogLbol = np.log10(fagn[i] * rest_frame_Lbol_erg_s)
LogL3000 = LogLbol - 0.7
M_i = (LogLbol - 36) / (-0.4)
# controlling the balance between strong, peaky, systemic emission and
# weak, highly skewed emission
# emline_type = np.random.uniform(low=-2, high=3, size=1)[0]
# L3000 re-scale the output model flux appropriately, M_i control the
# emission line properties
# Baldwin effect, EW of strong emission lines anticorrelated with
# intrinsic luminosity
quasar_sed = Quasar_sed(
z=redshift[i],
LogL3000=LogL3000,
wavlen=rest_frame_wave_A,
ebv=0.0,
M_i=M_i,
gflag=False,
) # , emline_type=emline_type)
quasar_fluxes_erg_s_cm2_A[i] = quasar_sed.flux
return quasar_fluxes_erg_s_cm2_A