# This is mostly from mosfit/transforms but rewritten to be modular
import numpy as np
from scipy.interpolate import interp1d
from redback.constants import *
[docs]
class Diffusion(object):
[docs]
def __init__(self, time, dense_times, luminosity, kappa, kappa_gamma, mej, vej, **kwargs):
"""
:param time: source frame time in days
:param dense_times: dense time array in days
:param dense_luminosity: luminosity
:param kappa: opacity
:param kappa_gamma: gamma-ray opacity
:param mej: ejecta mass
:param vej: ejecta velocity
Adds new attributes for tau_diffusion and new luminosity accounting for the
interaction process at the time values
"""
self.kappa = kappa
self.kappa_gamma = kappa_gamma
self.luminosity = luminosity
self.time = time
self.dense_times = dense_times
self.m_ejecta = mej
self.v_ejecta = vej
self.reference = 'https://ui.adsabs.harvard.edu/abs/1982ApJ...253..785A/abstract'
self.tau_d = []
self.new_luminosity = []
self.tau_d, self.new_luminosity = self.convert_input_luminosity()
def convert_input_luminosity(self):
timesteps = 100
minimum_log_spacing = -3
diffusion_constant = 2.0 * solar_mass / (13.7 * speed_of_light * km_cgs)
trapping_constant = 3.0 * solar_mass / (4*np.pi * km_cgs ** 2)
tau_diff = np.sqrt(diffusion_constant * self.kappa * self.m_ejecta / self.v_ejecta) / day_to_s
trap_coeff = (trapping_constant * self.kappa_gamma * self.m_ejecta / (self.v_ejecta ** 2)) / day_to_s ** 2
min_te = np.min(self.dense_times)
tb = max(0.0, min_te)
luminosity_interpolator = interp1d(self.dense_times, self.luminosity, copy=False, assume_sorted=True)
uniq_times = np.unique(self.time[(self.time >= tb) & (self.time <= self.dense_times[-1])])
lu = len(uniq_times)
num = int(round(timesteps / 2.0))
lsp = np.logspace(np.log10(tau_diff /self.dense_times[-1]) + minimum_log_spacing, 0, num)
xm = np.unique(np.concatenate((lsp, 1 - lsp)))
int_times = np.clip(tb + (uniq_times.reshape(lu, 1) - tb) * xm, tb, self.dense_times[-1])
int_te2s = int_times[:, -1] ** 2
int_lums = luminosity_interpolator(int_times)
int_args = int_lums * int_times * np.exp((int_times ** 2 - int_te2s.reshape(lu, 1)) / tau_diff**2)
int_args[np.isnan(int_args)] = 0.
uniq_lums = np.trapezoid(int_args, int_times, axis=1)
uniq_lums *= -2.0 * np.expm1(-trap_coeff / int_te2s) / tau_diff**2
new_lums = uniq_lums[np.searchsorted(uniq_times, self.time)]
return tau_diff, new_lums
[docs]
class AsphericalDiffusion(object):
[docs]
def __init__(self, time, dense_times, luminosity, kappa, kappa_gamma, mej, vej, area_projection, area_reference, **kwargs):
"""
:param time: source frame time in days
:param dense_times: dense time array in days
:param luminosity: luminosity
:param kappa: opacity
:param kappa_gamma: gamma-ray opacity
:param mej: ejecta mass
:param vej: ejecta velocity
:param area_projection: projected area of cocoon/polar ejecta
:param area_reference: remaining reference area i.e., the equitorial ejecta
Adds new attributes for tau_diffusion and new luminosity accounting for the interaction process
"""
self.kappa = kappa
self.kappa_gamma = kappa_gamma
self.luminosity = luminosity
self.time = time
self.dense_times = dense_times
self.m_ejecta = mej
self.v_ejecta = vej
self.area_projection = area_projection
self.area_reference = area_reference
self.reference = 'https://ui.adsabs.harvard.edu/abs/2020ApJ...897..150D/abstract'
self.tau_d = []
self.new_luminosity = []
self.tau_d, self.new_luminosity = self.convert_input_luminosity()
def convert_input_luminosity(self):
timesteps = 100
minimum_log_spacing = -3
diffusion_constant = 2.0 * solar_mass / (13.7 * speed_of_light * km_cgs)
trapping_constant = 3.0 * solar_mass / (4*np.pi * km_cgs ** 2)
tau_diff = np.sqrt(diffusion_constant * self.kappa * self.m_ejecta / self.v_ejecta) / day_to_s
trap_coeff = (trapping_constant * self.kappa_gamma * self.m_ejecta / (self.v_ejecta ** 2)) / day_to_s ** 2
min_te = min(self.dense_times)
tb = max(0.0, min_te)
luminosity_interpolator = interp1d(self.dense_times, self.luminosity, copy=False,assume_sorted=True)
uniq_times = np.unique(self.time[(self.time >= tb) & (self.time <= self.dense_times[-1])])
lu = len(uniq_times)
num = int(round(timesteps / 2.0))
lsp = np.logspace(np.log10(tau_diff /self.dense_times[-1]) + minimum_log_spacing, 0, num)
xm = np.unique(np.concatenate((lsp, 1 - lsp)))
int_times = np.clip(tb + (uniq_times.reshape(lu, 1) - tb) * xm, tb, self.dense_times[-1])
int_te2s = int_times[:, -1] ** 2
int_lums = luminosity_interpolator(int_times)
int_args = int_lums * int_times * np.exp((int_times ** 2 - int_te2s.reshape(lu, 1)) / tau_diff**2)
int_args[np.isnan(int_args)] = 0.0
uniq_lums = np.trapezoid(int_args, int_times, axis=1)
uniq_lums *= -2.0 * np.expm1(-trap_coeff / int_te2s) / tau_diff**2
uniq_lums *= (1 + 1.4 * (2 + uniq_times/tau_diff/0.59) / (1 + np.exp(uniq_times/tau_diff/0.59)) *
(self.area_projection/self.area_reference - 1))
new_lums = uniq_lums[np.searchsorted(uniq_times, self.time)]
return tau_diff, new_lums
[docs]
class CSMDiffusion(object):
[docs]
def __init__(self, time, dense_times, luminosity, kappa, r_photosphere, mass_csm_threshold, csm_mass, **kwargs):
"""
:param time: source frame time in days
:param dense_times: dense time array in days
:param luminosity: luminosity in erg/s
:param kappa: opacity in cm^2/g
:param r_photosphere: photosphere radius in cm (CGS)
:param mass_csm_threshold: mass of the optically-thick CSM (tau > 2/3) in grams (CGS),
as returned directly by _csm_engine. Used in CGS form in the diffusion timescale formula.
:param csm_mass: total csm mass in solar masses
:param csm_diffusion_method: Method used for the CSM diffusion integral.
Default is "cumulative". Use "quadrature" for the legacy per-time integral.
:param csm_diffusion_timesteps: Number of quadrature points used for each
requested time when csm_diffusion_method="quadrature". Default is 3000.
Adds new attribute for luminosity accounting for the interaction process
"""
self.time = time
self.dense_times = dense_times
self.luminosity = luminosity
self.kappa = kappa
self.r_photosphere = r_photosphere
self.mass_csm_threshold = mass_csm_threshold
self.csm_mass = csm_mass * solar_mass
self.timesteps = kwargs.get("csm_diffusion_timesteps", 3000)
self.diffusion_method = kwargs.get("csm_diffusion_method", "cumulative")
self.reference = 'https://ui.adsabs.harvard.edu/abs/2013ApJ...773...76C/abstract'
self.new_luminosity = []
self.new_luminosity = self.convert_input_luminosity()
def convert_input_luminosity(self):
if self.diffusion_method == "cumulative":
return self._convert_input_luminosity_cumulative()
if self.diffusion_method in ["quadrature", "legacy"]:
return self._convert_input_luminosity_quadrature()
raise ValueError("csm_diffusion_method must be 'cumulative' or 'quadrature'")
def _diffusion_timescale(self):
beta = 4. * np.pi ** 3. / 9.
return self.kappa * self.mass_csm_threshold / (beta * speed_of_light * self.r_photosphere) / day_to_s
def _convert_input_luminosity_cumulative(self):
t0 = self._diffusion_timescale()
dense_times = np.asarray(self.dense_times)
time = np.asarray(self.time)
valid_times = time[(time >= dense_times[0]) & (time <= dense_times[-1])]
integration_times = np.unique(np.concatenate((dense_times, valid_times)))
luminosity = np.interp(integration_times, dense_times, self.luminosity)
new_lums = np.zeros_like(integration_times, dtype=float)
dt = np.diff(integration_times)
luminosity_slopes = np.diff(luminosity) / dt
scaled_dt = dt / t0
exp_decay = np.exp(-scaled_dt)
one_minus_exp_decay = -np.expm1(-scaled_dt)
for ii in range(len(dt)):
new_lums[ii + 1] = (
new_lums[ii] * exp_decay[ii]
+ luminosity[ii] * one_minus_exp_decay[ii]
+ luminosity_slopes[ii] * (dt[ii] - t0 * one_minus_exp_decay[ii])
)
return np.interp(time, integration_times, new_lums)
def _convert_input_luminosity_quadrature(self):
timesteps = self.timesteps
minimum_log_spacing = -3
t0 = self._diffusion_timescale()
min_te = min(self.dense_times)
tb = max(0.0, min_te)
uniq_times = np.unique(self.time[(self.time >= tb) & (self.time <= self.dense_times[-1])])
lu = len(uniq_times)
num = int(round(timesteps / 2.0))
log_start = min(np.log10(t0 / self.dense_times[-1]) + minimum_log_spacing, 0.0)
lsp = np.logspace(log_start, 0, num)
xm = np.unique(np.concatenate((lsp, 1 - lsp)))
xm = xm[(xm >= 0.0) & (xm <= 1.0)]
xm = np.unique(np.concatenate(([0.0], xm, [1.0])))
int_times = tb + (uniq_times.reshape(lu, 1) - tb) * xm
int_tes = uniq_times.reshape(lu, 1)
int_lums = np.interp(int_times.ravel(), self.dense_times, self.luminosity).reshape(int_times.shape)
int_args = int_lums * np.exp((int_times - int_tes) / t0)
int_args[np.isnan(int_args)] = 0.0
uniq_lums = np.trapezoid(int_args, int_times, axis=1)
uniq_lums /= t0
new_lums = uniq_lums[np.searchsorted(uniq_times, self.time)]
return new_lums
[docs]
class Viscous(object):
[docs]
def __init__(self, time, dense_times, luminosity, t_viscous, **kwargs):
"""
:param time: source frame time in days
:param dense_times: dense time array in days
:param luminosity: luminosity
:param t_viscous: viscous timescale
Adds new attribute for luminosity accounting for the interaction process
"""
self.luminosity = luminosity
self.time = time
self.dense_times = dense_times
self.tvisc = t_viscous
self.reference = ''
self.new_luminosity = []
self.new_luminosity = self.convert_input_luminosity()
def convert_input_luminosity(self):
timesteps = 1000
minimum_log_spacing = -3
min_te = min(self.dense_times)
tb = max(0.0, min_te)
luminosity_interpolator = interp1d(self.dense_times, self.luminosity, copy=False,assume_sorted=True)
uniq_times = np.unique(self.time[(self.time >= tb) & (self.time <= self.dense_times[-1])])
lu = len(uniq_times)
num = int(round(timesteps / 2.0))
lsp = np.logspace(np.log10(self.tvisc /self.dense_times[-1]) +minimum_log_spacing, 0, num)
xm = np.unique(np.concatenate((lsp, 1 - lsp)))
int_times = np.clip(tb + (uniq_times.reshape(lu, 1) - tb) * xm, tb, self.dense_times[-1])
int_tes = int_times[:, -1]
int_lums = luminosity_interpolator(int_times)
int_args = int_lums * np.exp((int_times - int_tes.reshape(lu, 1)) / self.tvisc)
int_args[np.isnan(int_args)] = 0.0
uniq_lums = np.trapezoid(int_args, int_times, axis=1)/self.tvisc
new_lums = uniq_lums[np.searchsorted(uniq_times, self.time)]
return new_lums