Source code for sherpa.astro.optical

#
#  Copyright (C) 2011, 2016, 2017, 2018, 2019, 2020, 2021, 2023
#  Smithsonian Astrophysical Observatory
#
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 3 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  You should have received a copy of the GNU General Public License along
#  with this program; if not, write to the Free Software Foundation, Inc.,
#  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#

"""
Optical models intended for SED Analysis

The models match those used by the SpecView application [1]_,
and are intended for un-binned one-dimensional data sets defined
on a wavelength grid, with units of Angstroms. When used with
a binned data set the lower-edge of each bin is used to evaluate
the model. This module does not contain all the spectral
components from SpecView ([2]_).

References
----------

.. [1] http://www.stsci.edu/institute/software_hardware/specview/

.. [2] http://specview.stsci.edu/javahelp/Components.html

"""

import numpy
from sherpa.models.parameter import Parameter, tinyval
from sherpa.models.model import ArithmeticModel, RegriddableModel1D
from sherpa.utils import sao_fcmp
from sherpa.utils.err import ModelErr
from sherpa.utils.numeric_types import SherpaFloat

_tol = numpy.finfo(numpy.float32).eps

# Optical Models for SED Analysis
#
# This Sherpa Python module contains optical models for fitting to SEDs.
# These models are Python versions of models found in the Specview
# application for analyzing spectra and SEDs.  These models are meant
# to be used in conjunction with Specview to serve the VAO SED project.
#
# These models work in wavelength space (Angstroms).
#

__all__ = ('AbsorptionEdge', 'AccretionDisk', 'AbsorptionGaussian',
           'AbsorptionLorentz', 'EmissionLorentz', 'OpticalGaussian',
           'EmissionGaussian', 'BlackBody',
           'AbsorptionVoigt', 'EmissionVoigt',  # Use Voigt1D instead
           'Bremsstrahlung', 'BrokenPowerlaw', 'CCM', 'LogAbsorption',
           'LogEmission', 'Polynomial', 'Powerlaw', 'Recombination',
           'XGal', 'FM', 'LMC', 'SM', 'SMC', 'Seaton')

# The speed of light in km/s
c_km = 2.99792458e+5


# Helper function, to do a particular sort of crude interpolation
# from tables supplied for certain extinction curves.  No general
# applicability outside this function, so we do not make it available
# to code outside this module.  Ported from the Specview file
# AbstractExtinction.java for Iris use.  5/26/11 SMD
def _extinct_interp(xtable, etable, x):
    out = numpy.zeros_like(x)
    last = len(xtable) - 1

    for i in range(len(x)):
        xval = x[i]
        if (xval <= xtable[0]):
            out[i] = etable[0]
        elif (xval >= xtable[last]):
            out[i] = etable[last]
        else:
            index = 0
            for j in range(last + 1):
                if xval < xtable[j]:
                    index = j
                    break
            x1 = xtable[index - 1]
            x2 = xtable[index]
            e1 = etable[index - 1]
            e2 = etable[index]

            out[i] = ((e2 - e1) / (x2 - x1)) * (xval - x1) + e1

    return out


[docs] class AbsorptionVoigt(ArithmeticModel): """This model has been replaced by Voigt1D .. versionchanged:: 4.12.2 The AbsorptionVoigt model has been removed. The Voigt1D model should be used instead. """ def __init__(self, *args): raise ModelErr("The AbsorptionVoigt model has been replaced by Voigt1D")
[docs] class EmissionVoigt(ArithmeticModel): """This model has been replaced by Voigt1D .. versionchanged:: 4.12.2 The EmissionVoigt model has been removed. The Voigt1D model should be used instead. """ def __init__(self, *args): raise ModelErr("The EmissionVoigt model has been replaced by Voigt1D")
# This model sets in edge (in Angstroms) beyond which absorption # is a significant feature to the spectrum or SED. #
[docs] class AbsorptionEdge(RegriddableModel1D): """Optical model of an absorption edge. This model is intended to be used to modify another model (e.g. by multiplying the two together). It is for use when the independent axis is in wavelength units (e.g. Angstrom). Attributes ---------- egdew The location of the edge. Above this value the model is set to 1. tau The optical depth of the edge. index The exponent used for the relative distance from the edge. It is a hidden parameter, with a value fixed at 3. See Also -------- AbsorptionGaussian, AbsorptionLorentz, OpticalGaussian Notes ----- The functional form of the model for points is:: f(x) = exp(-tau * (x / edgew)^index) for x <= edgew = 1 otherwise and for integrated data sets the low-edge of the grid is used. """ def __init__(self, name='absorptionedge'): self.edgew = Parameter(name, 'edgew', 5000., tinyval, frozen=True, units='angstroms') self.tau = Parameter(name, 'tau', 0.5) self.index = Parameter(name, 'index', 3.0, alwaysfrozen=True, hidden=True) ArithmeticModel.__init__(self, name, (self.edgew, self.tau, self.index)) # We can turn on model caching with this commented-out feature, # if we find we need it. # @modelCacher1d
[docs] def calc(self, p, x, xhi=None, **kwargs): x = numpy.asarray(x, dtype=SherpaFloat) y = numpy.ones_like(x) if 0.0 == p[0]: raise ValueError('model evaluation failed, ' + '%s edgew cannot be zero' % self.name) idx = (x <= p[0]) y[idx] = numpy.exp(-(p[1] * numpy.power(x[idx] / p[0], p[2]))) return y
# This model is an accretion disk continuum function.
[docs] class AccretionDisk(RegriddableModel1D): """A model of emission due to an accretion disk. It is for use when the independent axis is in Angstroms. Attributes ---------- ref The reference wavelength, in Angstroms. beta The index of the power-law component. ampl The amplitude of the disk. norm The normalization value for the position. It is a hidden parameter, with a value fixed at 20000. Notes ----- The functional form of the model for points is:: f(x) = ampl * (x / norm)^(-beta) * exp(-ref / x) and for integrated data sets the low-edge of the grid is used. """ def __init__(self, name='accretiondisk'): self.ref = Parameter(name, 'ref', 5000., frozen=True, units='angstroms') self.beta = Parameter(name, 'beta', 0.5, -10, 10) self.ampl = Parameter(name, 'ampl', 1.) self.norm = Parameter(name, 'norm', 20000.0, tinyval, alwaysfrozen=True, hidden=True) ArithmeticModel.__init__(self, name, (self.ref, self.beta, self.ampl, self.norm)) # @modelCacher1d
[docs] def calc(self, p, x, xhi=None, **kwargs): x = numpy.asarray(x, dtype=SherpaFloat) if 0.0 in x: raise ValueError('model evaluation failed, ' + 'x cannot be zero') if p[3] == 0.0: raise ValueError('model evaluation failed, ' + 'norm cannot be zero') return p[2] * numpy.power(x / p[3], -p[1]) * numpy.exp(-p[0] / x)
# This model calculates a Gaussian function expressed in # equivalent width, and models absorption due to this Gaussian.
[docs] class AbsorptionGaussian(RegriddableModel1D): """Gaussian function for modeling absorption (equivalent width). This model is intended to be used to modify another model (e.g. by multiplying the two together). It is for use when the independent axis is in Angstroms. Attributes ---------- fwhm The full-width half-maximum of the model in km/s. pos The center of the gaussian, in Angstroms. ewidth The equivalent width of the model, in Angstroms. limit The model is only evaluated for points that lie within limit sigma of pos. It is a hidden parameter, with a value fixed at 4. See Also -------- AbsorptionEdge, AbsorptionLorentz, EmissionGaussian, OpticalGaussian Notes ----- The functional form of the model for points is:: f(x) = 1 - ampl * exp(-0.5 * ((x - pos)/sigma)^2) ampl = ewidth / sigma / 2.50662828 sigma = pos * fwhm / (2.354820044 * c) and for integrated data sets the low-edge of the grid is used. The calculation is only done for those points that are in the range:: |x - pos| < limit * sigma Outside this range the model is set to 1. """ def __init__(self, name='absorptiongaussian'): self.fwhm = Parameter(name, 'fwhm', 100., tinyval, hard_min=tinyval, units="km/s") self.pos = Parameter(name, 'pos', 5000., tinyval, frozen=True, units='angstroms') self.ewidth = Parameter(name, 'ewidth', 1.) self.limit = Parameter(name, 'limit', 4., alwaysfrozen=True, hidden=True) ArithmeticModel.__init__(self, name, (self.fwhm, self.pos, self.ewidth, self.limit)) # @modelCacher1d
[docs] def calc(self, p, x, xhi=None, **kwargs): x = numpy.asarray(x, dtype=SherpaFloat) if 0.0 == p[0]: raise ValueError('model evaluation failed, ' + '%s fwhm cannot be zero' % self.name) if 0.0 == p[1]: raise ValueError('model evaluation failed, ' + '%s pos cannot be zero' % self.name) y = numpy.ones_like(x) sigma = p[1] * p[0] / 705951.5 # = 2.9979e5 / 2.354820044 ? delta = numpy.abs((x - p[1]) / sigma) ampl = p[2] / sigma / 2.50662828 # document this constant idx = (delta < self.limit.val) y[idx] = 1.0 - ampl * numpy.exp(- delta[idx] * delta[idx] / 2.0) return y
# This model calculates a Lorentzian function expressed in # equivalent width, and models absorption due to this Lorentzian.
[docs] class AbsorptionLorentz(RegriddableModel1D): """Lorentz function for modeling absorption (equivalent width). This model is intended to be used to modify another model (e.g. by multiplying the two together). It is for use when the independent axis is in Angstroms. Attributes ---------- fwhm The full-width half-maximum of the model in km/s. pos The center of the lorentzian, in Angstroms. ewidth The equivalent width of the model, in Angstroms. See Also -------- AbsorptionEdge, AbsorptionGaussian, EmissionLorentz, OpticalGaussian Notes ----- The functional form of the model for points is:: f(x) = 1 - ewidth * c / (1.571 * fwhm * pos * l(x)) l(x) = 1 + 4 * ((1 / x - 1 / pos) * pos * c / fwhm)^2 c = speed of light in km/s and for integrated data sets the low-edge of the grid is used. The speed of light can be found by inspecting the module variable ``sherpa.astro.optical.c_km``. """ def __init__(self, name='absorptionlorentz'): self.fwhm = Parameter(name, 'fwhm', 100., tinyval, hard_min=tinyval, units="km/s") self.pos = Parameter(name, 'pos', 5000., tinyval, frozen=True, units='angstroms') self.ewidth = Parameter(name, 'ewidth', 1.) ArithmeticModel.__init__(self, name, (self.fwhm, self.pos, self.ewidth)) # @modelCacher1d
[docs] def calc(self, p, x, xhi=None, **kwargs): x = numpy.asarray(x, dtype=SherpaFloat) if 0.0 == p[0]: raise ValueError('model evaluation failed, ' + '%s fwhm cannot be zero' % self.name) if 0.0 == p[1]: raise ValueError('model evaluation failed, ' + '%s pos cannot be zero' % self.name) y = (1.0 / x - 1.0 / p[1]) * p[1] * c_km / p[0] y = 1.0 + 4.0 * y * y y *= 1.571 * p[0] * p[1] / c_km y = 1.0 - p[2] / y return y
# This model computes a Lorentzian profile for emission features.
[docs] class EmissionLorentz(RegriddableModel1D): """Lorentz function for modeling emission. It is for use when the independent axis is in Angstroms. Attributes ---------- fwhm The full-width half-maximum of the model in km/s. pos The center of the lorentzian, in Angstroms. flux The normalisation of the lorentzian. kurt The kurtosis of the lorentzian. See Also -------- AbsorptionLorentz, EmissionGaussian Notes ----- The functional form of the model for points is:: f(x) = flux * 2 * pi * s / l(x) l(x) = abs(x - pos)^kurt + (0.5 * s)^2 s = pos * fwhm / c c = speed of light in km/s and for integrated data sets the low-edge of the grid is used. The speed of light can be found by inspecting the module variable ``sherpa.astro.optical.c_km``. """ def __init__(self, name='emissionlorentz'): self.fwhm = Parameter(name, 'fwhm', 100., tinyval, hard_min=tinyval, units="km/s") self.pos = Parameter(name, 'pos', 5000., tinyval, frozen=True, units='angstroms') self.flux = Parameter(name, 'flux', 1.) self.kurt = Parameter(name, 'kurt', 2., frozen=True) ArithmeticModel.__init__(self, name, (self.fwhm, self.pos, self.flux, self.kurt)) # @modelCacher1d
[docs] def calc(self, p, x, xhi=None, **kwargs): x = numpy.asarray(x, dtype=SherpaFloat) if 0.0 == p[0]: raise ValueError('model evaluation failed, ' + '%s fwhm cannot be zero' % self.name) if 0.0 == p[1]: raise ValueError('model evaluation failed, ' + '%s pos cannot be zero' % self.name) sigma = p[0] * p[1] / c_km arg = numpy.power(numpy.abs(x - p[1]), p[3]) + \ sigma / 2.0 * sigma / 2.0 arg[arg < 1.0e-15] = 1.0e-15 return p[2] * sigma / arg / (numpy.pi * 2)
# This model computes an absorption Gaussian feature expressed in # optical depth.
[docs] class OpticalGaussian(RegriddableModel1D): """Gaussian function for modeling absorption (optical depth). This model is intended to be used to modify another model (e.g. by multiplying the two together). It is for use when the independent axis is in Angstroms. Attributes ---------- fwhm The full-width half-maximum of the model in km/s. pos The center of the gaussian, in Angstroms. tau The optical depth of the model. limit The model is only evaluated for points that lie within limit sigma of pos. It is a hidden parameter, with a value fixed at 4. See Also -------- AbsorptionEdge, AbsorptionGaussian, AbsorptionLorentz, EmissionGaussian Notes ----- The functional form of the model for points is:: f(x) = exp(-tau * g(x)) g(x) = exp(-0.5 * ((x - pos) / sigma)^2) sigma = pos * fwhm / (2.9979e5 * 2.354820044) and for integrated data sets the low-edge of the grid is used. The calculation is only done for those points that are in the range:: |x - pos| < limit * sigma Outside this range the model is set to 1. """ def __init__(self, name='opticalgaussian'): self.fwhm = Parameter(name, 'fwhm', 100., tinyval, hard_min=tinyval, units="km/s") self.pos = Parameter(name, 'pos', 5000., tinyval, frozen=True, units='angstroms') self.tau = Parameter(name, 'tau', 0.5) self.limit = Parameter(name, 'limit', 4., alwaysfrozen=True, hidden=True) ArithmeticModel.__init__(self, name, (self.fwhm, self.pos, self.tau, self.limit)) # @modelCacher1d
[docs] def calc(self, p, x, xhi=None, **kwargs): x = numpy.asarray(x, dtype=SherpaFloat) if 0.0 == p[0]: raise ValueError('model evaluation failed, ' + '%s fwhm cannot be zero' % self.name) if 0.0 == p[1]: raise ValueError('model evaluation failed, ' + '%s pos cannot be zero' % self.name) y = numpy.ones_like(x) sigma = p[1] * p[0] / 705951.5 # = 2.9979e5 / 2.354820044 ? delta = numpy.abs((x - p[1]) / sigma) idx = (delta < self.limit.val) y[idx] = numpy.exp(-p[2] * numpy.exp(- delta[idx] * delta[idx] / 2.0)) return y
# This model computes a Gaussian profile for emission features.
[docs] class EmissionGaussian(RegriddableModel1D): """Gaussian function for modeling emission. It is for use when the independent axis is in Angstroms. Attributes ---------- fwhm The full-width half-maximum of the model in km/s. pos The center of the gaussian, in Angstroms. flux The normalisation of the gaussian. skew The skew of the gaussian. limit The model is only evaluated for points that lie within limit sigma of pos. It is a hidden parameter, with a default value of 4. See Also -------- AbsorptionGaussian, EmissionLorentz, LogEmission Notes ----- The functional form of the model for points is:: f(x) = flux * exp(-0.5 * d(x)^2) / s2 if skew = 1 = 2 * flux * exp(-0.5 * d2(x)^2) / (s2 * (1 + skew)) otherwise d(x) = (x - pos) / s d2(x) = d(x) if x <= pos = d(x) / skew otherwise s2 = 2.50662828 * s s = pos * fwhm / (2.9979e5 * 2.354820044) and for integrated data sets the low-edge of the grid is used. The calculation is only done for those points that are in the range:: |x - pos| < limit * sigma Outside this range the model is set to 0. """ def __init__(self, name='emissiongaussian'): self.fwhm = Parameter(name, 'fwhm', 100., tinyval, hard_min=tinyval, units="km/s") self.pos = Parameter(name, 'pos', 5000., tinyval, frozen=True, units='angstroms') self.flux = Parameter(name, 'flux', 1.) self.skew = Parameter(name, 'skew', 1., tinyval, frozen=True) self.limit = Parameter(name, 'limit', 4., alwaysfrozen=True, hidden=True) ArithmeticModel.__init__(self, name, (self.fwhm, self.pos, self.flux, self.skew, self.limit)) # @modelCacher1d
[docs] def calc(self, p, x, xhi=None, **kwargs): x = numpy.asarray(x, dtype=SherpaFloat) fwhm = p[0] pos = p[1] flux = p[2] skew = p[3] if 0.0 == fwhm: raise ValueError('model evaluation failed, ' + '%s fwhm cannot be zero' % self.name) if 0.0 == pos: raise ValueError('model evaluation failed, ' + '%s pos cannot be zero' % self.name) if 0.0 == skew: raise ValueError('model evaluation failed, ' + '%s skew cannot be zero' % self.name) y = numpy.zeros_like(x) sigma = pos * fwhm / 705951.5 # = 2.9979e5 / 2.354820044 delta = numpy.abs((x - pos) / sigma) idx = (delta < self.limit.val) arg = - delta * delta / 2.0 s2 = 2.50662828 * sigma if sao_fcmp(skew, 1.0, _tol) == 0: y[idx] = flux * numpy.exp(arg[idx]) / s2 else: left = (x <= pos) arg[left] = numpy.exp(arg[left]) right = ~left arg[right] = numpy.exp(arg[right] / skew / skew) y[idx] = 2.0 * flux * arg[idx] / s2 / (1.0 + skew) return y
# This model computes continuum emission as a blackbody function.
[docs] class BlackBody(RegriddableModel1D): """Emission from a black body as a function of wavelength. It is for use when the independent axis is in Angstroms. Attributes ---------- refer The reference point, in Angstroms. ampl The amplitude of the emission; it is defined at the reference point but its numerical value there also depends on the temperature. temperature The temperature in Kelvin. See Also -------- Bremsstrahlung Notes ----- The functional form of the model for points is:: f(x) = ampl * g(refer) / g(x) g(x) = x^5 * (exp(1.438786E8 / temperature / x) - 1) and for integrated data sets the low-edge of the grid is used. """ def __init__(self, name='blackbody'): self.refer = Parameter(name, 'refer', 5000., tinyval, hard_min=tinyval, frozen=True, units="angstroms") self.ampl = Parameter(name, 'ampl', 1., tinyval, hard_min=tinyval, units="angstroms") self.temperature = Parameter(name, 'temperature', 3000., tinyval, hard_min=tinyval, units="Kelvin") self._argmin = 1.0e-3 self._argmax = 1000.0 ArithmeticModel.__init__(self, name, (self.refer, self.ampl, self.temperature)) # @modelCacher1d
[docs] def calc(self, p, x, xhi=None, **kwargs): x = numpy.asarray(x, dtype=SherpaFloat) c1 = 1.438786e8 efactor = c1 / p[2] if ((efactor / p[0]) > self._argmax): # raise error exp too big raise ValueError('model evaluation failed, either ' + 'temperature or reference wavelength too small') numer = p[1] * numpy.power(p[0], 5.0) * \ (numpy.exp(efactor / p[0]) - 1.0) y = numpy.zeros_like(x) x0 = numpy.where(x > 0.0)[0] if (len(x0) > 0): arg = numpy.zeros_like(x) arg[x0] = efactor / x[x0] denon = numpy.zeros_like(x) denon[x0] = numpy.power(x[x0], 5) argmin_slice = numpy.where(arg < self._argmin)[0] if (len(argmin_slice) > 0): denon[argmin_slice] *= arg[argmin_slice] * \ (1.0 + 0.5 * arg[argmin_slice]) arg = numpy.where(arg > self._argmax, self._argmax, arg) arg_slice = numpy.where(arg >= self._argmin)[0] if (len(arg_slice) > 0): denon[arg_slice] *= numpy.exp(arg[arg_slice]) - 1.0 y[x0] = numer / denon[x0] return y
# This model computes continuum emission with the bremsstrahlung function.
[docs] class Bremsstrahlung(RegriddableModel1D): """Bremsstrahlung emission. It is for use when the independent axis is in Angstroms. Attributes ---------- refer The reference point, in Angstroms. ampl The amplitude of the emission. Note that the model does not equal ``ampl`` at the reference point, as shown in the functional form below. temperature The temperature in Kelvin. See Also -------- BlackBody Notes ----- The functional form of the model for points is:: f(x) = ampl * (refer/x)^2 * exp(-1.438779E8 / temperature / x) and for integrated data sets the low-edge of the grid is used. """ def __init__(self, name='bremsstrahlung'): self.refer = Parameter(name, 'refer', 5000., tinyval, hard_min=tinyval, frozen=True, units="angstroms") self.ampl = Parameter(name, 'ampl', 1., tinyval, hard_min=tinyval, units="angstroms") self.temperature = Parameter(name, 'temperature', 3000., tinyval, hard_min=tinyval, units="Kelvin") ArithmeticModel.__init__(self, name, (self.refer, self.ampl, self.temperature)) # @modelCacher1d
[docs] def calc(self, p, x, xhi=None, **kwargs): x = numpy.asarray(x, dtype=SherpaFloat) if 0.0 == p[0]: raise ValueError('model evaluation failed, ' + '%s refer cannot be zero' % self.name) return p[1] * numpy.power((p[0] / x), 2) * \ numpy.exp(-1.438779e8 / x / p[2])
# This model computes continuum emission with a broken power-law; # that is, the power-law index changes after a break at a particular # wavelength.
[docs] class BrokenPowerlaw(RegriddableModel1D): """Broken power-law model. It is for use when the independent axis is in Angstroms. Attributes ---------- refer The reference point at which the amplitude is defined, with units of Angstroms. ampl The amplitude at the reference point. index1 The index for the power law below the reference point. index2 The index for the power law above the reference point. See Also -------- Polynomial, Powerlaw Notes ----- The functional form of the model for points is:: f(x) = ampl * (x / refer)^index1 x < refer = ampl * (x / refer)^index2 x >= refer and for integrated data sets the low-edge of the grid is used. """ def __init__(self, name='brokenpowerlaw'): self.refer = Parameter(name, 'refer', 5000., tinyval, hard_min=tinyval, frozen=True, units="angstroms") self.ampl = Parameter(name, 'ampl', 1., tinyval, hard_min=tinyval, units="angstroms") self.index1 = Parameter(name, 'index1', 0.1, -10.0, 10.0) self.index2 = Parameter(name, 'index2', -0.1, -10.0, 10.0) ArithmeticModel.__init__(self, name, (self.refer, self.ampl, self.index1, self.index2)) # @modelCacher1d
[docs] def calc(self, p, x, xhi=None, **kwargs): if 0.0 == p[0]: raise ValueError('model evaluation failed, ' + '%s refer cannot be zero' % self.name) x = numpy.asarray(x, dtype=SherpaFloat) arg = x / p[0] arg = p[1] * (numpy.where(arg > 1.0, numpy.power(arg, p[3]), numpy.power(arg, p[2]))) return arg
# This model computes extinction using the function published by # Cardelli, Clayton, and Mathis # (ApJ, 1989, vol 345, pp 245)
[docs] class CCM(RegriddableModel1D): """Galactic extinction: the Cardelli, Clayton, and Mathis model. The interstellar extinction is calculated using the formula from [1]_. This model is intended to be used to modify another model (e.g. by multiplying the two together). It is for use when the independent axis is in Angstrom. Attributes ---------- ebv E(B-V) r R_v See Also -------- FM, LMC, Seaton, SM, SMC, XGAL Notes ----- When evaluated on a binned grid, the lower-edges of the bins are used for the calculation. References ---------- .. [1] Cardelli, Clayton, & Mathis, 1989, ApJ, 345, 245-256. http://adsabs.harvard.edu/abs/1989ApJ...345..245C """ def __init__(self, name='ccm'): self.ebv = Parameter(name, 'ebv', 0.5) self.r = Parameter(name, 'r', 3.2) ArithmeticModel.__init__(self, name, (self.ebv, self.r)) # @modelCacher1d
[docs] def calc(self, p, x, xhi=None, **kwargs): x = numpy.asarray(x, dtype=SherpaFloat) y = numpy.zeros_like(x) y2 = numpy.zeros_like(x) y3 = numpy.zeros_like(x) a = numpy.zeros_like(x) b = numpy.zeros_like(x) x = 1000.0 / (x / 10.0) # Infrared wavelengths xp = numpy.zeros_like(x) ir_slice = numpy.where((x >= 0.3) & (x <= 1.1))[0] if (len(ir_slice) > 0): xp[ir_slice] = numpy.power(x[ir_slice], 1.61) a[ir_slice] = 0.574 * xp[ir_slice] b[ir_slice] = -0.527 * xp[ir_slice] # Optical opt_slice = numpy.where((x > 1.1) & (x <= 3.3))[0] if (len(opt_slice) > 0): y[opt_slice] = x[opt_slice] - 1.82 a[opt_slice] = 1.0 + 0.17699 * y[opt_slice] \ - 0.50477 * y[opt_slice] * y[opt_slice] \ - 0.02427 * numpy.power(y[opt_slice], 3) \ + 0.72085 * numpy.power(y[opt_slice], 4) \ + 0.01979 * numpy.power(y[opt_slice], 5) \ - 0.77530 * numpy.power(y[opt_slice], 6) \ + 0.32999 * numpy.power(y[opt_slice], 7) b[opt_slice] = 0.0 + 1.41338 * y[opt_slice] \ + 2.28305 * y[opt_slice] * y[opt_slice] \ + 1.07233 * numpy.power(y[opt_slice], 3) \ - 5.38434 * numpy.power(y[opt_slice], 4) \ - 0.62551 * numpy.power(y[opt_slice], 5) \ + 5.30260 * numpy.power(y[opt_slice], 6) \ - 2.09002 * numpy.power(y[opt_slice], 7) # Near-UV nuv_slice = numpy.where((x > 3.3) & (x <= 8.0))[0] if (len(nuv_slice) > 0): a[nuv_slice] = 0.0 b[nuv_slice] = 0.0 nuv_slice2 = numpy.where((x >= 5.9) & (x <= 8.0))[0] if (len(nuv_slice2) > 0): y[nuv_slice2] = x[nuv_slice2] - 5.9 y2[nuv_slice2] = y[nuv_slice2] * y[nuv_slice2] y3[nuv_slice2] = y2[nuv_slice2] * y[nuv_slice2] a[nuv_slice2] = -0.04473 * y2[nuv_slice2] \ - 0.009779 * y3[nuv_slice2] b[nuv_slice2] = 0.21300 * y2[nuv_slice2] \ + .120700 * y3[nuv_slice2] a[nuv_slice] = a[nuv_slice] + 1.752 \ - 0.316 * x[nuv_slice] \ - 0.104 / (0.341 + numpy.power((x[nuv_slice] - 4.67), 2)) b[nuv_slice] = b[nuv_slice] - 3.090 \ + 1.825 * x[nuv_slice] \ + 1.206 / (0.263 + numpy.power((x[nuv_slice] - 4.62), 2)) # Far-UV fuv_slice = numpy.where((x > 8.0) & (x <= 20.0))[0] if (len(fuv_slice) > 0): y[fuv_slice] = x[fuv_slice] - 8.0 y2[fuv_slice] = y[fuv_slice] * y[fuv_slice] y3[fuv_slice] = y2[fuv_slice] * y[fuv_slice] a[fuv_slice] = -1.073 - 0.628 * y[fuv_slice] \ + 0.137 * y2[fuv_slice] - 0.070 * y3[fuv_slice] b[fuv_slice] = 13.670 + 4.257 * y[fuv_slice] \ - 0.420 * y2[fuv_slice] + 0.374 * y3[fuv_slice] # Final extinction curve aext = p[1] * a + b return numpy.power(10.0, (-0.4 * p[0] * aext))
# This model computes absorption using a Gaussian function expressed # in optical depth, and using the log of the FWHM.
[docs] class LogAbsorption(RegriddableModel1D): """Gaussian function for modeling absorption (log of fwhm). This model is intended to be used to modify another model (e.g. by multiplying the two together). It is for use when the independent axis is in Angstroms. Attributes ---------- fwhm The full-width half-maximum of the feature in km/s. pos The center of the feature, in Angstroms. tau The optical depth of the feature. See Also -------- AbsorptionEdge, AbsorptionGaussian, AbsorptionLorentz, EmissionGaussian, LogEmission, OpticalGaussian Notes ----- The functional form of the model for points is:: f(x) = exp(-tau * (x / pos)^(term * alpha)) term = -1 if x >= pos = 1 otherwise alpha = log(2) / log(1 + 0.5 * fwhm / c) and for integrated data sets the low-edge of the grid is used. """ def __init__(self, name='logabsorption'): self.fwhm = Parameter(name, 'fwhm', 100., tinyval, hard_min=tinyval, units="km/s") self.pos = Parameter(name, 'pos', 5000., tinyval, frozen=True, units='angstroms') self.tau = Parameter(name, 'tau', 0.5) ArithmeticModel.__init__(self, name, (self.fwhm, self.pos, self.tau)) # @modelCacher1d
[docs] def calc(self, p, x, xhi=None, **kwargs): x = numpy.asarray(x, dtype=SherpaFloat) if 0.0 == p[0]: raise ValueError('model evaluation failed, ' + '%s fwhm cannot be zero' % self.name) if 0.0 == p[1]: raise ValueError('model evaluation failed, ' + '%s pos cannot be zero' % self.name) y = numpy.ones_like(x) alpha = 0.69314718 / numpy.log(1.0 + p[0] / 2.9979e5 / 2.0) if (alpha <= 1.0): alpha = 1.0001 y = numpy.where(x >= p[1], p[2] * numpy.power((x / p[1]), -alpha), p[2] * numpy.power((x / p[1]), alpha)) return numpy.exp(-y)
# This model computes emission using a Gaussian function expressed # in optical depth, and using the log of the FWHM. # # DOC NOTE: the specview docs and ahelp file claim that fmax # requires c but the code uses the pos parameter. # WHAT IS CORRECT? See # https://github.com/sherpa/sherpa/issues/220 #
[docs] class LogEmission(RegriddableModel1D): """Gaussian function for modeling emission (log of fwhm). It is for use when the independent axis is in Angstroms. Attributes ---------- fwhm The full-width half-maximum of the model in km/s. pos The center of the gaussian, in Angstroms. flux The normalisation of the gaussian. skew The skew of the gaussian. See Also -------- EmissionGaussian, EmissionLorentz, LogAbsorption Notes ----- The functional form of the model for points is:: f(x) = fmax * (x / pos)^arg if x <= pos fmax * (x / pos)^(-alpha) otherwise arg = log(2) / log(1 + 0.5 * fwhm / c) arg1 = log(2) / log(1 + 0.5 * skew * fwhm / c) alpha = arg if skew == 1 = arg1 otherwise fmax = (arg - 1) * flux / (2 * pos) if skew == 1 = (arg - 1) * flux / (pos * (1 + (arg - 1) / (arg1 - 1))) otherwise c = 2.9979e5 and for integrated data sets the low-edge of the grid is used. """ def __init__(self, name='logemission'): self.fwhm = Parameter(name, 'fwhm', 100., tinyval, hard_min=tinyval, units="km/s") self.pos = Parameter(name, 'pos', 5000., tinyval, frozen=True, units='angstroms') self.flux = Parameter(name, 'flux', 1.) self.skew = Parameter(name, 'skew', 1., tinyval, frozen=True) ArithmeticModel.__init__(self, name, (self.fwhm, self.pos, self.flux, self.skew)) # @modelCacher1d
[docs] def calc(self, p, x, xhi=None, **kwargs): x = numpy.asarray(x, dtype=SherpaFloat) if 0.0 == p[0]: raise ValueError('model evaluation failed, ' + '%s fwhm cannot be zero' % self.name) if 0.0 == p[1]: raise ValueError('model evaluation failed, ' + '%s pos cannot be zero' % self.name) if 0.0 == p[3]: raise ValueError('model evaluation failed, ' + '%s skew cannot be zero' % self.name) arg = 0.69314718 / numpy.log(1.0 + p[0] / 2.9979e5 / 2.0) if (arg <= 1.0): arg = 1.0001 fmax = (arg - 1.0) * p[2] / p[1] / 2.0 if (p[3] == 1.0): return numpy.where(x >= p[1], fmax * numpy.power((x / p[1]), -arg), fmax * numpy.power((x / p[1]), arg)) arg1 = 0.69314718 / numpy.log(1.0 + p[3] * p[0] / 2.9979e5 / 2.0) fmax = (arg - 1.0) * p[2] / p[1] / (1.0 + (arg - 1.0) / (arg1 - 1.0)) return numpy.where(x <= p[1], fmax * numpy.power((x / p[1]), arg), fmax * numpy.power((x / p[1]), -arg1))
# This model computes continuum emission as a polynomial, # y = c0 + # c1 * (x - offset) + # c2 * (x - offset)^2 + # c3 * (x - offset)^3 + # c4 * (x - offset)^4 + # c5 * (x - offset)^5 #
[docs] class Polynomial(RegriddableModel1D): """Polynomial model of order 5. This model can be used with any one-dimensional data set since there are no units on the parameters. Attributes ---------- c0 The constant term. c1 The amplitude of the (x-offset) term. c2 The amplitude of the (x-offset)^2 term. c3 The amplitude of the (x-offset)^3 term. c4 The amplitude of the (x-offset)^4 term. c5 The amplitude of the (x-offset)^5 term. offset There is a degeneracy between ``c0`` and ``offset``, so it is recommended that at least one of these remains frozen. See Also -------- Powerlaw Notes ----- The functional form of the model for points is:: f(x) = sum_(i=0)^(i=8) c_i * (x - offset)^i and for integrated data sets the low-edge of the grid is used. """ def __init__(self, name='polynomial'): pars = [] for i in range(6): pars.append(Parameter(name, 'c%d' % i, 0, frozen=True)) pars[0].val = 1 pars[0].frozen = False for p in pars: setattr(self, p.name, p) self.offset = Parameter(name, 'offset', 0, frozen=True) pars.append(self.offset) ArithmeticModel.__init__(self, name, pars) # @modelCacher1d
[docs] def calc(self, p, x, xhi=None, **kwargs): x = numpy.asarray(x, dtype=SherpaFloat) y = numpy.zeros_like(x) xtemp = x - p[6] y += p[5] for i in [4, 3, 2, 1, 0]: y = y * xtemp + p[i] return y
# This model computes continuum emission using a power-law.
[docs] class Powerlaw(RegriddableModel1D): """Power-law model. It is for use when the independent axis is in Angstroms. Attributes ---------- refer The reference point at which the amplitude is defined, with units of Angstroms. ampl The amplitude at the reference point. index The index for the power law. See Also -------- BrokenPowerlaw, Polynomial Notes ----- The functional form of the model for points is:: f(x) = ampl * (x / refer)^index and for integrated data sets the low-edge of the grid is used. """ def __init__(self, name='powerlaw'): self.refer = Parameter(name, 'refer', 5000., tinyval, hard_min=tinyval, frozen=True, units="angstroms") self.ampl = Parameter(name, 'ampl', 1., tinyval, hard_min=tinyval, units="angstroms") self.index = Parameter(name, 'index', -0.5, -10.0, 10.0) ArithmeticModel.__init__(self, name, (self.refer, self.ampl, self.index)) # @modelCacher1d
[docs] def calc(self, p, x, xhi=None, **kwargs): if 0.0 == p[0]: raise ValueError('model evaluation failed, ' + '%s refer cannot be zero' % self.name) x = numpy.asarray(x, dtype=SherpaFloat) arg = x / p[0] arg = p[1] * numpy.power(arg, p[2]) return arg
# This model computes the continuum with an optically thin # recombination function.
[docs] class Recombination(RegriddableModel1D): """Optically-thin recombination continuum model. It is for use when the independent axis is in Angstroms. Attributes ---------- refer The reference point, in Angstroms. ampl The amplitude of the emission; it is defined at the reference point but its numerical value there also depends on the temperature. temperature The temperature in Kelvin. fwhm The full-width half-maximum of the model in km/s. Notes ----- The functional form of the model for points is:: f(x) = ampl * (refer / x)^2 * exp(-1.440E8 * (1 / x - 1 / refer) / temperature) if x < refer = ampl * exp(-0.5 * (x - refer)^2 / g(fwhm, refer)^2) otherwise g(fwhm, refer) = refer * fwhm / (2.354820044 * c) where c is the speed of light in km/s. For integrated data sets the low-edge of the grid is used. """ def __init__(self, name='recombination'): self.refer = Parameter(name, 'refer', 5000., tinyval, hard_min=tinyval, frozen=True, units="angstroms") self.ampl = Parameter(name, 'ampl', 1., tinyval, hard_min=tinyval, units="angstroms") self.temperature = Parameter(name, 'temperature', 3000., tinyval, hard_min=tinyval, units="Kelvin") self.fwhm = Parameter(name, 'fwhm', 100., tinyval, hard_min=tinyval, units="km/s") ArithmeticModel.__init__(self, name, (self.refer, self.ampl, self.temperature, self.fwhm)) # @modelCacher1d
[docs] def calc(self, p, x, xhi=None, **kwargs): if 0.0 == p[0]: raise ValueError('model evaluation failed, ' + '%s refer cannot be zero' % self.name) x = numpy.asarray(x, dtype=SherpaFloat) sigma = p[0] * p[3] / 705951.5 # = 2.9979e5 / 2.354820044 delta = 1.440e8 * (1.0 / x - 1.0 / p[0]) / p[2] return numpy.where(delta < 0.0, p[1] * numpy.exp(-numpy.power((x - p[0]), 2.0) / numpy.power(sigma, 2.0) / 2.0), p[1] * numpy.power((p[0] / x), 2.0) * numpy.exp(-delta) )
# This model computes the extragalactic extinction function of # Calzetti, Kinney and Storchi-Bergmann, 1994, ApJ, 429, 582
[docs] class XGal(RegriddableModel1D): """Extragalactic extinction: Calzetti, Kinney and Storchi-Bergmann The extragalactic extinction is calculated using the formula from [1]_. This model is intended to be used to modify another model (e.g. by multiplying the two together). It is for use when the independent axis is in Angstrom. Attributes ---------- ebv E(B-V) See Also -------- CCM, FM, LMC, Seaton, SM, SMC Notes ----- When evaluated on a binned grid, the lower-edges of the bins are used for the calculation. References ---------- .. [1] Calzetti, Kinney, Storchi-Bergmann, 1994, ApJ, 429, 582 http://adsabs.harvard.edu/abs/1994ApJ...429..582C """ def __init__(self, name='xgal'): self.ebv = Parameter(name, 'ebv', 0.5) ArithmeticModel.__init__(self, name, (self.ebv,)) # @modelCacher1d
[docs] def calc(self, p, x, xhi=None, **kwargs): x = numpy.asarray(x, dtype=SherpaFloat) if 0.0 in x: raise ValueError('model evaluation failed, ' + 'x cannot be zero') x = 1000.0 / x # Formula from paper with zero point moved to (x = 0) ext = ((0.011 * x - 0.198) * x + 1.509) * x # Normalize the result according to Kailash Sahu's calculations ext *= 2.43 return numpy.power(10.0, (-0.4 * p[0] * ext))
# This model computes the extinction curve for wavelengths for UV # spectra below 3200 A. Fitzpatrick and Massa 1988 extinction curve # with Drude UV bump. # See Fitzpatrick and Massa (ApJ, 1988, vol. 328, p. 734)
[docs] class FM(RegriddableModel1D): """UV extinction curve: Fitzpatrick and Massa 1988. The UV extinction is calculated using [1]_. This model is intended to be used to modify another model (e.g. by multiplying the two together). It is for use when the independent axis is in Angstrom. Attributes ---------- ebv E(B-V) x0 Position of the Drude bump. width Width of the Drude bump. c1 The intercept of the linear term. c2 The slope of the linear term. c3 Normalization of the Drude bump. c4 Normalization of the FUV curvature. See Also -------- CCM, LMC, Seaton, SM, SMC, XGAL Notes ----- When evaluated on a binned grid, the lower-edges of the bins are used for the calculation. References ---------- .. [1] Fitzpatrick and Massa 1988 http://adsabs.harvard.edu/abs/1988ApJ...328..734F """ def __init__(self, name='fm'): self.ebv = Parameter(name, 'ebv', 0.5) # E(B-V) self.x0 = Parameter(name, 'x0', 4.6) # Position of Drude bump self.width = Parameter(name, 'width', 0.06) # Width of Drude bump self.c1 = Parameter(name, 'c1', 0.2) self.c2 = Parameter(name, 'c2', 0.1) self.c3 = Parameter(name, 'c3', 0.02) self.c4 = Parameter(name, 'c4', 0.1) ArithmeticModel.__init__(self, name, (self.ebv, self.x0, self.width, self.c1, self.c2, self.c3, self.c4)) # @modelCacher1d
[docs] def calc(self, p, x, xhi=None, **kwargs): x = numpy.asarray(x, dtype=SherpaFloat) x = 10000. / x av = 3.14 dru = x * x / (p[2] * p[2] * x * x + numpy.power((x - p[1]), 2.0)) dx = x - 5.9 fuv = 0.5392 * dx * dx + 0.0564 * numpy.power(dx, 3.0) ext = av + p[3] + p[4] * x + p[5] * dru ext = numpy.where(x > 5.9, p[6] * fuv + ext, ext) return numpy.power(10.0, (-0.4 * ext * p[0]))
# This model computes the extinction curve using the # LMC extinction curve from Howarth 1983 MNRAS, 203, 301
[docs] class LMC(RegriddableModel1D): """LMC extinction: the Howarth model. The interstellar extinction is calculated using the formula from [1]_. This model is intended to be used to modify another model (e.g. by multiplying the two together). It is for use when the independent axis is in Angstrom. Attributes ---------- ebv E(B-V) See Also -------- CCM, FM, Seaton, SM, SMC, XGAL Notes ----- When evaluated on a binned grid, the lower-edges of the bins are used for the calculation. References ---------- .. [1] Howarth 1983 MNRAS, 203, 301 http://adsabs.harvard.edu/abs/1983MNRAS.203..301H """ def __init__(self, name='lmc'): self.ebv = Parameter(name, 'ebv', 0.5) self._R = 3.1 ArithmeticModel.__init__(self, name, (self.ebv,)) # @modelCacher1d
[docs] def calc(self, p, x, xhi=None, **kwargs): x = numpy.asarray(x, dtype=SherpaFloat) # convert from wavelength in Angstroms to 1/microns x = 10000.0 / x extmag = numpy.zeros_like(x) # Infrared - extend optical results linearly to 0 at 1/lambda = 0 slice1 = numpy.where(x <= 1.83)[0] slice2 = numpy.where((x > 1.83) & (x <= 2.75))[0] slice3 = numpy.where(x > 2.75)[0] x = numpy.where(x > 10.96, 10.96, x) if (len(slice1) > 0): extmag[slice1] = ((1.86 - 0.48 * x[slice1]) * x[slice1] - 0.1) * \ x[slice1] if (len(slice2) > 0): extmag[slice2] = self._R + 2.04 * (x[slice2] - 1.83) + \ 0.094 * (x[slice2] - 1.83) * (x[slice2] - 1.83) # continue out to lambda = 912 A if (len(slice3) > 0): extmag[slice3] = self._R - 0.236 + 0.462 * x[slice3] + \ 0.105 * x[slice3] * x[slice3] + \ 0.454 / ((x[slice3] - 4.557) * (x[slice3] - 4.557) + 0.293) return numpy.power(10.0, (-0.4 * extmag * p[0]))
# This model computes the galactic interstellar extinction function # from Savage & Mathis (1979, ARA&A, 17, 73-111)
[docs] class SM(RegriddableModel1D): """Galactic extinction: the Savage & Mathis model. The interstellar extinction is calculated using the formula from [1]_. This model is intended to be used to modify another model (e.g. by multiplying the two together). It is for use when the independent axis is in Angstrom. Attributes ---------- ebv E(B-V) See Also -------- CCM, FM, LMC, Seaton, SMC, XGAL Notes ----- When evaluated on a binned grid, the lower-edges of the bins are used for the calculation. References ---------- .. [1] Savage & Mathis, 1979, ARA&A, 17, 73-111 http://adsabs.harvard.edu/abs/1979ARA%26A..17...73S """ def __init__(self, name='sm'): self.ebv = Parameter(name, 'ebv', 0.5) self._xtab = numpy.array([0.00, 0.29, 0.45, 0.80, 1.11, 1.43, 1.82, 2.27, 2.50, 2.91, 3.65, 4.00, 4.17, 4.35, 4.57, 4.76, 5.00, 5.26, 5.56, 5.88, 6.25, 6.71, 7.18, 8.00, 8.50, 9.00, 9.50, 10.00]) self._extab = numpy.array([0.00, 0.16, 0.38, 0.87, 1.50, 2.32, 3.10, 4.10, 4.40, 4.90, 6.20, 7.29, 8.00, 8.87, 9.67, 9.33, 8.62, 8.00, 7.75, 7.87, 8.12, 8.15, 8.49, 9.65, 10.55, 11.55, 12.90, 14.40]) ArithmeticModel.__init__(self, name, (self.ebv,)) # @modelCacher1d
[docs] def calc(self, p, x, xhi=None, **kwargs): x = numpy.asarray(x, dtype=SherpaFloat) # convert from wavelength in Angstroms to 1/microns x = 10000.0 / x extmag = numpy.zeros_like(x) extmag = _extinct_interp(self._xtab, self._extab, x) return numpy.power(10.0, (-0.4 * extmag * p[0]))
# This model computes the SMC extinction function of # Prevot et al. 1984, A&A, 132, 389-392
[docs] class SMC(RegriddableModel1D): """SMC extinction: the Prevot et al. 1984 model. The interstellar extinction is calculated using the formula from [1]_. This model is intended to be used to modify another model (e.g. by multiplying the two together). It is for use when the independent axis is in Angstrom. Attributes ---------- ebv E(B-V) See Also -------- CCM, FM, LMC, Seaton, SM, XGAL Notes ----- When evaluated on a binned grid, the lower-edges of the bins are used for the calculation. References ---------- .. [1] Prevot et al. 1984, A&A, 132, 389-392 http://adsabs.harvard.edu/abs/1984A%26A...132..389P """ def __init__(self, name='smc'): self.ebv = Parameter(name, 'ebv', 0.5) self._xtab = numpy.array([0.00, 0.29, 0.45, 0.80, 1.11, 1.43, 1.82, 2.35, 2.70, 3.22, 3.34, 3.46, 3.60, 3.75, 3.92, 4.09, 4.28, 4.50, 4.73, 5.00, 5.24, 5.38, 5.52, 5.70, 5.88, 6.07, 6.27, 6.48, 6.72, 6.98, 7.23, 7.52, 7.84]) self._extab = numpy.array([-3.10, -2.94, -2.72, -2.23, -1.60, -0.78, 0.00, 1.00, 1.67, 2.29, 2.65, 3.00, 3.15, 3.49, 3.91, 4.24, 4.53, 5.30, 5.85, 6.38, 6.76, 6.90, 7.17, 7.71, 8.01, 8.49, 9.06, 9.28, 9.84, 10.80, 11.51, 12.52, 13.54]) ArithmeticModel.__init__(self, name, (self.ebv,)) # @modelCacher1d
[docs] def calc(self, p, x, xhi=None, **kwargs): x = numpy.asarray(x, dtype=SherpaFloat) # convert from wavelength in Angstroms to 1/microns x = 10000.0 / x extmag = numpy.zeros_like(x) extmag = _extinct_interp(self._xtab, self._extab, x) return numpy.power(10.0, (-0.4 * extmag * p[0]))
# This model computes Seaton's interstellar extinction function. # The formulae are based on an adopted value of R = 3.20. # # This function implements Seaton's function as originally # implemented in STScI's Synphot program. # # For wavelengths > 3704 Angstrom, the function interpolates # linearly in 1/lambda in Seaton's table 3. For wavelengths # < 3704 Angstrom, the class uses the formulae from Seaton's # table 2. The formulae match at the endpoints of their respective # intervals. There is a mismatch of 0.009 mag/ebmv at nu=2.7 # (lambda=3704 Angstrom). Seaton's tabulated value of 1.44 mag at # 1/lambda = 1.1 may be in error; 1.64 seems more consistent with # his other values. # # Wavelength range allowed is 0.1 to 1.0 microns; outside this # range, the class extrapolates the function. # # References: # # lambda < 1000 same as lambda = 1000. # 1000 < lambda < 3704 Seaton (1979) MNRAS 187,73p. # 3704 < lambda < 10,000 Nandy (1975) A+A 44, 195. (corrected to R=3.2) # 10000 < lambda extrapolate linearly in 1/lambda
[docs] class Seaton(RegriddableModel1D): """Galactic extinction: the Seaton model from Synphot. The interstellar extinction is calculated using the formula from [1]_ as implemented in STSCI's Synphot program [2]_. The supported wavelength range is 1000 to 10000 Angstroms, and the Notes section describes the changes from [1]_. This model is intended to be used to modify another model (e.g. by multiplying the two together). It is for use when the independent axis is in Angstrom. Attributes ---------- ebv E(B-V) See Also -------- CCM, FM, LMC, SM, SMC, XGAL Notes ----- The formulae are based on an adopted value of R=3.20. For wavelengths above 3704 Angstrom, the function interpolates linearly in 1/lambda in table 3 [1]_. For wavelengths below this, the formulae from table 2 [1]_ are used (see also [3]_, corrected to R=3.2). The formulae match at the endpoints of their respective intervals. There is a mismatch of 0.009 mag/ebmv at nu=2.7 (lambda=3704 Angstrom). Seaton's tabulated value of 1.44 mag at 1/lambda = 1.1 may be in error; 1.64 seems more consistent with his other values. For wavelengths below 1000 Angstrom, a constant value equal to the value at 1000 Angstrom is used. For wavelengths above 10000 Angstroms a linear extrapolation (in 1/lambda) is used. When evaluated on a binned grid, the lower-edges of the bins are used for the calculation. References ---------- .. [1] Seaton, M. J. 1979, MNRAS, 187, 73 http://adsabs.harvard.edu/abs/1979MNRAS.187P..73S .. [2] http://www.stsci.edu/institute/software_hardware/stsdas/synphot .. [3] Nandy et al., 1975, A&A, 44, 195-203. http://adsabs.harvard.edu/abs/1975A%26A....44..195N """ def __init__(self, name='seaton'): self.ebv = Parameter(name, 'ebv', 0.5) self._xtab = numpy.array([0.0, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7]) self._extab = numpy.array([0.0, 1.36, 1.64, 1.84, 2.04, 2.24, 2.44, 2.66, 2.88, 3.14, 3.36, 3.56, 3.77, 3.96, 4.15, 4.26, 4.40, 4.52, 4.64]) ArithmeticModel.__init__(self, name, (self.ebv,)) # @modelCacher1d
[docs] def calc(self, p, x, xhi=None, **kwargs): x = numpy.asarray(x, dtype=SherpaFloat) # convert from wavelength in Angstroms to 1/microns x = 10000.0 / x extmag = numpy.zeros_like(x) ir_slice = numpy.where(x <= 1.0)[0] opt_slice = numpy.where((x > 1.0) & (x < 2.7))[0] uv1_slice = numpy.where((x >= 2.7) & (x < 3.65))[0] uv2_slice = numpy.where((x >= 3.65) & (x < 7.14))[0] uv3_slice = numpy.where((x >= 7.14) & (x <= 10.0))[0] uv_extra_slice = numpy.where(x > 10.0)[0] # Infrared - extend optical results linearly to 0 at 1/lambda = 0 if (len(ir_slice) > 0): extmag[ir_slice] = self._extab[1] * x[ir_slice] * x[ir_slice] # Optical - interpolate in Seaton's table 3 if (len(opt_slice) > 0): extmag[opt_slice] = _extinct_interp(self._xtab, self._extab, x[opt_slice]) # UV - use analytic formulae from Seaton's table 2 if (len(uv1_slice) > 0): extmag[uv1_slice] = 1.56 + 1.048 * x[uv1_slice] + \ 1.01 / ((x[uv1_slice] - 4.6) * (x[uv1_slice] - 4.6) + 0.280) # UV again if (len(uv2_slice) > 0): extmag[uv2_slice] = 2.29 + 0.848 * x[uv2_slice] + \ 1.01 / ((x[uv2_slice] - 4.6) * (x[uv2_slice] - 4.6) + 0.280) # and more UV still if (len(uv3_slice) > 0): extmag[uv3_slice] = 16.17 + \ x[uv3_slice] * (-3.20 + 0.2975 * x[uv3_slice]) # Extrapolate beyond 1/lambda = 10.0 if (len(uv_extra_slice) > 0): x[uv_extra_slice] = numpy.where(x[uv_extra_slice] < 50.0, x[uv_extra_slice], 50.0) extmag[uv_extra_slice] = 16.17 + \ x[uv_extra_slice] * (-3.20 + 0.2975 * x[uv_extra_slice]) return numpy.power(10.0, (-0.4 * extmag * p[0]))