Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

adding phase dependence to prop. effects #359

Merged
merged 8 commits into from
Apr 25, 2023
Merged
Show file tree
Hide file tree
Changes from 5 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
42 changes: 31 additions & 11 deletions sncosmo/models.py
Original file line number Diff line number Diff line change
Expand Up @@ -1546,26 +1546,31 @@ def _flux(self, time, wave):
"""Array flux function."""

a = 1. / (1. + self._parameters[0])
phase = (time - self._parameters[1]) * a
obsphase = (time - self._parameters[1])
restphase = obsphase * a
restwave = wave * a

# Note that below we multiply by the scale factor to conserve
# bolometric luminosity.
f = a * self._source._flux(phase, restwave)
f = a * self._source._flux(restphase, restwave)

# Pass the flux through the PropagationEffects.
for effect, frame, zindex in zip(self._effects, self._effect_frames,
self._effect_zindicies):
if frame == 'obs':
effect_wave = wave
effect_phase = obsphase
elif frame == 'rest':
effect_wave = restwave
effect_phase = restphase
else: # frame == 'free'
effect_a = 1. / (1. + self._parameters[zindex])
effect_wave = wave * effect_a

f = effect.propagate(effect_wave, f)

effect_phase = obsphase * effect_a
try:
f = effect.propagate(effect_wave, f, phase=effect_phase)
except:
f = effect.propagate(effect_wave, f)
return f

def flux(self, time, wave):
Expand Down Expand Up @@ -1948,15 +1953,24 @@ def minwave(self):
def maxwave(self):
return self._maxwave

def minphase(self):
return self._minphase

def maxphase(self):
return self._maxphase
benjaminrose marked this conversation as resolved.
Show resolved Hide resolved

@abc.abstractmethod
def propagate(self, wave, flux):
def propagate(self, wave, flux, phase=None):
pass

def _headsummary(self):
summary = """\
class : {0}
wavelength range: [{1:.6g}, {2:.6g}] Angstroms"""\
.format(self.__class__.__name__, self._minwave, self._maxwave)
wavelength range: [{1:.6g}, {2:.6g}] Angstroms
phase range : [{3:.2g}, {4:.2g}]"""\
.format(self.__class__.__name__,
self._minwave, self._maxwave,
self._minphase, self._maxphase)
return dedent(summary)


Expand All @@ -1966,11 +1980,13 @@ class CCM89Dust(PropagationEffect):
param_names_latex = ['E(B-V)', 'R_V']
_minwave = 1000.
_maxwave = 33333.33
_minphase = np.nan
_maxphase = np.nan

def __init__(self):
self._parameters = np.array([0., 3.1])

def propagate(self, wave, flux):
def propagate(self, wave, flux, phase=None):
"""Propagate the flux."""
ebv, r_v = self._parameters
return extinction.apply(extinction.ccm89(wave, ebv * r_v, r_v), flux)
Expand All @@ -1982,11 +1998,13 @@ class OD94Dust(PropagationEffect):
param_names_latex = ['E(B-V)', 'R_V']
_minwave = 909.09
_maxwave = 33333.33
_minphase = np.nan
_maxphase = np.nan

def __init__(self):
self._parameters = np.array([0., 3.1])

def propagate(self, wave, flux):
def propagate(self, wave, flux, phase=None):
"""Propagate the flux."""
ebv, r_v = self._parameters
return extinction.apply(extinction.odonnell94(wave, ebv * r_v, r_v),
Expand All @@ -1997,6 +2015,8 @@ class F99Dust(PropagationEffect):
"""Fitzpatrick (1999) extinction model dust with fixed R_V."""
_minwave = 909.09
_maxwave = 60000.
_minphase = np.nan
_maxphase = np.nan

def __init__(self, r_v=3.1):
self._param_names = ['ebv']
Expand All @@ -2005,7 +2025,7 @@ def __init__(self, r_v=3.1):
self._r_v = r_v
self._f = extinction.Fitzpatrick99(r_v=r_v)

def propagate(self, wave, flux):
def propagate(self, wave, flux, phase=None):
"""Propagate the flux."""
ebv = self._parameters[0]
return extinction.apply(self._f(wave, ebv * self._r_v), flux)
59 changes: 57 additions & 2 deletions sncosmo/tests/test_models.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,8 @@

from io import StringIO

import scipy

import numpy as np
from numpy.testing import assert_allclose, assert_approx_equal

Expand All @@ -18,9 +20,48 @@ def flatsource():
return sncosmo.TimeSeriesSource(phase, wave, flux)


class AchromaticMicrolensing(sncosmo.PropagationEffect):
"""
An achromatic microlensing object. Tests phase dependence.
"""
_param_names = []
param_names_latex = []
_minwave = np.nan
_maxwave = np.nan
_minphase = np.nan
_maxphase = np.nan

def __init__(self, time, magnification):
"""
Parameters
----------
time: :class:`~list` or :class:`~numpy.array`
A time array for your microlensing
dmag: :class:`~list` or :class:`~numpy.array`
microlensing magnification
"""
self._parameters = np.array([])
self.mu = scipy.interpolate.interp1d(time,
magnification,
bounds_error=False,
fill_value=1.)
self._minphase = np.min(time)
self._maxphase = np.max(time)

def propagate(self, wave, flux, phase):
"""
Propagate the magnification onto the model's flux output.
"""

mu = np.expand_dims(np.atleast_1d(self.mu(phase)), 1)
return flux * mu


class StepEffect(sncosmo.PropagationEffect):
"""Effect with transmission 0 below cutoff wavelength, 1 above.
Useful for testing behavior with redshift."""
"""
Effect with transmission 0 below cutoff wavelength, 1 above.
Useful for testing behavior with redshift.
"""

_param_names = ['stepwave']
param_names_latex = [r'\lambda_{s}']
Expand Down Expand Up @@ -312,3 +353,17 @@ def test_effect_frame_free():

wave = np.array([12000., 13000., 14000., 14999., 15000., 16000.])
assert_allclose(model.flux(0., wave), [0., 0., 0., 0., 0.5, 0.5])


def test_effect_phase_dependent():
"""Test Model with PropagationEffect with phase dependence"""

model = sncosmo.Model(source=flatsource())
flux = model.flux(50., 5000.).flatten()

model_micro = sncosmo.Model(source=flatsource(),
effects=[AchromaticMicrolensing([40., 60.],
[10., 10.])],
effect_frames=['rest'],
effect_names=['micro'])
assert_approx_equal(model_micro.flux(50., 5000.).flatten(), flux*10.)