# -*- coding: utf-8 -*-
"""
Module for computing star formation rates with time delays.
This module provides functions for computing the star formation rate at a given
redshift, accounting for time delays between formation and observation. The time
delay distribution follows a 1/t power-law form, and the formation redshift is
computed using the cosmological age-redshift relation.
Key Features: \n
- Time-delayed star formation rate computation \n
- Integration with Madau & Fragos (2017) SFR model \n
- Monte Carlo integration for time delay averaging \n
- Cosmological calculations using Astropy \n
Copyright (C) 2026 Hemanta Ph. Distributed under MIT License.
"""
import numpy as np
from numba import set_num_threads
from scipy.integrate import quad
from scipy.optimize import fsolve
from astropy.cosmology import LambdaCDM
from .prior_functions import sfr_madau_fragos2017
[docs]
def sfr_with_time_delay_function(input_args):
"""
Compute star formation rate at observed redshift with time delay.
The star formation rate is time-delayed relative to the observed redshift,
with a time delay uniformly distributed between td_min and td_max. The
formation redshift is computed using the cosmological age-redshift relation.
Parameters
----------
input_args : ``list``
List containing the following elements in order: \n
- z (``float``): Observed redshift \n
- idx (``int``): Index identifier for the computation \n
- td_min (``float``): Minimum time delay (Gyr) \n
- td_max (``float``): Maximum time delay (Gyr) \n
- H0 (``float``): Hubble constant (km/s/Mpc) \n
- Omega_M (``float``): Matter density parameter \n
- Omega_Lambda (``float``): Dark energy density parameter \n
- a (``float``): Madau-Fragos SFR normalization parameter \n
- b (``float``): Madau-Fragos low-z power-law slope \n
- c (``float``): Madau-Fragos turnover parameter \n
- d (``float``): Madau-Fragos high-z power-law slope \n
Returns
-------
idx : ``int``
Index identifier (same as input).
result : ``float``
Time-averaged star formation rate at observed redshift z.
Examples
--------
>>> from ler.gw_source_population.sfr_with_time_delay import sfr_with_time_delay
>>> args = [0.5, 0, 0.02, 13.0, 70.0, 0.3, 0.7, 0.01, 2.6, 3.2, 6.2]
>>> idx, sfr = sfr_with_time_delay(args)
"""
set_num_threads(1)
z = input_args[0]
idx = input_args[1]
td_min = input_args[2]
td_max = input_args[3]
H0 = input_args[4]
Omega_M = input_args[5]
Omega_Lambda = input_args[6]
a = input_args[7]
b = input_args[8]
c = input_args[9]
d = input_args[10]
def _E(z_prime):
"""Compute dimensionless Hubble parameter E(z)."""
return np.sqrt(Omega_M * (1 + z_prime)**3 + Omega_Lambda)
def _integrand(z_prime):
"""Compute integrand for lookback time calculation."""
return 1 / (H0 * (1 + z_prime) * _E(z_prime)) * 977.813
def _time_delay(zform, z):
"""Compute time delay between formation and observation redshifts."""
integral, _ = quad(_integrand, z, zform)
return integral
def _equation_to_solve(zform, z, td):
"""Residual equation for finding formation redshift."""
return td - _time_delay(zform, z)
def _find_zform(z, td):
"""Find formation redshift for given observation redshift and time delay."""
zform_solution = fsolve(_equation_to_solve, z, args=(z, td))
return zform_solution
def _integrand_rates(z, size=100000, zform_max=1000.):
"""
Compute time-averaged SFR using Monte Carlo integration.
Parameters
----------
z : ``float``
Observed redshift.
size : ``int``
Number of Monte Carlo samples. \n
default: 100000
zform_max : ``float``
Maximum formation redshift cutoff. \n
default: 1000.0
Returns
-------
integral : ``float``
Time-averaged star formation rate.
"""
td = np.random.uniform(td_min, td_max, size)
td_max_allowed = _time_delay(zform_max, z)
idx_valid = np.where(td < td_max_allowed)[0]
P_td = np.zeros_like(td)
P_td[idx_valid] = 1 / (np.log(td_max / td_min) * td[idx_valid])
zform = np.zeros_like(td)
for idx_ in idx_valid:
zform[idx_] = _find_zform(z, td[idx_])
psi = sfr_madau_fragos2017(zform, a, b, c, d)
integral = 1 / (td_max - td_min) * np.sum(P_td * psi)
return integral
result = _integrand_rates(z)
return int(idx), result