# coding: utf-8
#
# Copyright (C) 2007, 2016, 2018, 2019, 2021 - 2024
# 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.
#
import numpy
from sherpa.models.basic import clean_kwargs1d, clean_kwargs2d
from sherpa.models.parameter import Parameter, tinyval
from sherpa.models.model import ArithmeticModel, RegriddableModel2D, RegriddableModel1D, modelCacher1d
from sherpa.astro.utils import apply_pileup
from sherpa.utils import bool_cast, lgam
from sherpa.utils.err import ModelErr
from sherpa.utils.guess import _guess_ampl_scale, get_fwhm, \
get_peak, get_position, guess_amplitude, guess_amplitude2d, \
guess_amplitude_at_ref, guess_fwhm, guess_position, \
guess_radius, guess_reference, param_apply_limits
import sherpa.models._modelfcts
from . import _modelfcts
__all__ = ('Atten', 'BBody', 'BBodyFreq', 'Beta1D', 'BPL1D', 'Dered', 'Edge',
'LineBroad', 'Lorentz1D', 'Voigt1D', 'PseudoVoigt1D', 'NormBeta1D', 'Schechter',
'Beta2D', 'DeVaucouleurs2D', 'HubbleReynolds', 'Lorentz2D',
'JDPileup', 'Sersic2D', 'Disk2D',
'Shell2D')
[docs]
class Atten(RegriddableModel1D):
"""Model the attenuation by the Inter-Stellar Medium (ISM).
This model calculates the transmission of the interstellar medium
using the description of the ISM absorption of [1]_. It includes
neutral He autoionization features. Between 1.2398 and 43.655
Angstroms (i.e. in the 0.28-10 keV range) the model also accounts
for metals as described in [2]_. It should only be used when the
independent axis has units of Angstroms.
Attributes
----------
hcol
The column density of HI in atoms cm^-2.
heiRatio
The ratio of the HeI to HI column densities.
heiiRatio
The ratio of the HeII to HI column densities.
Notes
-----
The code uses the best available photoionization cross-sections to
date from the atomic data literature and combines them in an
arbitrary mixture of the three ionic species: HI, HeI, and HeII.
This model provided courtesy of Pat Jelinsky.
The grid version is evaluated by numerically intgerating the
function over each bin using a non-adaptive Gauss-Kronrod scheme
suited for smooth functions [3]_, falling over to a simple
trapezoid scheme if this fails.
References
----------
.. [1] Rumph, Bowyer, & Vennes 1994, AJ 107, 2108
http://adsabs.harvard.edu/abs/1994AJ....107.2108R
.. [2] Morrison & McCammon 1983, ApJ 270, 119
http://adsabs.harvard.edu/abs/1983ApJ...270..119M
.. [3] https://www.gnu.org/software/gsl/manual/html_node/QNG-non_002dadaptive-Gauss_002dKronrod-integration.html
"""
def __init__(self, name='atten'):
self.hcol = Parameter(name, 'hcol', 1e+20, 1e+17, 1e+24, tinyval)
self.heiRatio = Parameter(name, 'heiRatio', 0.1, 0, 1)
self.heiiRatio = Parameter(name, 'heiiRatio', 0.01, 0, 1)
ArithmeticModel.__init__(self, name,
(self.hcol, self.heiRatio, self.heiiRatio))
[docs]
@modelCacher1d
def calc(self, p, *args, **kwargs):
kwargs = clean_kwargs1d(self, kwargs)
# atten should act like xsphabs, never integrate.
kwargs['integrate'] = False
return _modelfcts.atten(p, *args, **kwargs)
[docs]
class BBody(RegriddableModel1D):
"""A one-dimensional Blackbody model.
A model representing the ideal blackbody function. It can be
used when the independent axis is in energy or wavelength space.
Attributes
----------
space
Switch to select whether the independent axis is energy or
wavelength. This parameter is not fit (``alwaysfrozen`` is
set), and should be set to either 0, when the independent axis
is energy with units of keV, or 1 when the axis is wavelength
with units of Angstrom.
kT
The temperature if the blackbody, in keV.
ampl
The amplitude of the blackbody component.
See Also
--------
BBodyFreq
Notes
-----
The blackbody emission is calculated as a function of energy using
the expression::
f(E) = ampl * E^2 / (exp(E / kT) - 1)
where E is the photon energy, and kT is the blackbody temperature
(both in keV). The amplitude, ampl, is related to the ratio of
source radius to distance by::
ampl = (2 * pi / (c^2 * h^3)) (R / d)^2
= 9.884 x 10^31 (R / d)^2
with Planck's constant (h) specified in keV-s and the speed of
light (c) specified in cm/s, and with R and d representing the
radius of, and distance to, the source respectively.
There are two conditions when the above equation is not used:
- if E/kt < 10^-4 then f(E) = ampl * E * kT
- if E/kT > 60, f(E) = 0.
"""
def __init__(self, name='bbody'):
self.space = Parameter(name, 'space', 0, 0, 1, 0, 1,
'0 - energy | 1 - wave', alwaysfrozen=True)
self.kT = Parameter(name, 'kT', 1, 0.1, 1000, 0, 1e+10, 'keV')
self.ampl = Parameter(name, 'ampl', 1, 1e-20, 1e+20, 1e-20, 1e+20)
ArithmeticModel.__init__(self, name, (self.space, self.kT, self.ampl))
[docs]
def guess(self, dep, *args, **kwargs):
if args[0][0] > args[0][-1]:
self.space.val = 1
Emax = get_peak(dep, *args)
tMax = Emax / 1.594
kt = {'val': tMax,
'min': tMax / _guess_ampl_scale,
'max': tMax * _guess_ampl_scale}
param_apply_limits(kt, self.kt, **kwargs)
norm = guess_amplitude(dep, *args)
modampl = norm['val'] * (numpy.exp(Emax / tMax) - 1) / \
numpy.square(Emax)
mod = {'val': modampl,
'min': modampl / _guess_ampl_scale,
'max': modampl * _guess_ampl_scale}
param_apply_limits(mod, self.ampl, **kwargs)
[docs]
@modelCacher1d
def calc(self, p, *args, **kwargs):
kwargs = clean_kwargs1d(self, kwargs)
return _modelfcts.bbody(p, *args, **kwargs)
[docs]
class BBodyFreq(RegriddableModel1D):
"""A one-dimensional Blackbody model (frequency).
This model can be used when the independent axis is in frequency
space.
Attributes
----------
T
The temperature if the blackbody, in Kelvin.
ampl
The amplitude of the blackbody component.
See Also
--------
BBody
Notes
-----
The blackbody emission is calculated as a function of frequency
(v) using Wien's law (hv >> kT):
f(v) = ampl * 2 * h * v^3 * exp(-h * v / (k * T)) / c^2
where T is the blackbody temperature in Kelvin, h is Planck's
constant, k is Boltzmann's constant, and c the speed of light.
"""
def __init__(self, name='bbodyfreq'):
self.T = Parameter(name, 'T', 1e+06, 1000, 1e+10, 1000, 1e+10)
self.ampl = Parameter(name, 'ampl', 1, 0, hard_min=0)
ArithmeticModel.__init__(self, name, (self.T, self.ampl))
[docs]
def guess(self, dep, *args, **kwargs):
vmax = get_peak(dep, *args)
tMax = vmax / 5.88e+10
t = {'val': tMax,
'min': tMax / _guess_ampl_scale,
'max': tMax * _guess_ampl_scale}
norm = guess_amplitude(dep, *args)
c_cm = 2.99792458e+10
h_erg = 6.6260693e-27
factor = numpy.exp(2.82) * numpy.square(c_cm) / h_erg / 2.
modampl = norm['val'] * factor / numpy.power(vmax, 3.)
mod = {'val': modampl,
'min': modampl / _guess_ampl_scale,
'max': modampl * _guess_ampl_scale}
param_apply_limits(mod, self.ampl, **kwargs)
param_apply_limits(t, self.t, **kwargs)
[docs]
@modelCacher1d
def calc(self, p, *args, **kwargs):
kwargs = clean_kwargs1d(self, kwargs)
return _modelfcts.bbodyfreq(p, *args, **kwargs)
[docs]
class Beta1D(RegriddableModel1D):
"""One-dimensional beta model function.
The beta model is a Lorentz model with a varying power law.
Attributes
----------
r0
The core radius.
beta
This parameter controls the slope of the profile at large
radii.
xpos
The reference point of the profile. This is frozen by default.
ampl
The amplitude refers to the maximum value of the model, at
x = xpos.
See Also
--------
Beta2D, Lorentz1D, NormBeta1D
Notes
-----
The functional form of the model for points is::
f(x) = ampl * (1 + ((x - xpos) / r0)^2)^(0.5 - 3 * beta)
The grid version is evaluated by numerically intgerating the
function over each bin using a non-adaptive Gauss-Kronrod scheme
suited for smooth functions [1]_, falling over to a simple
trapezoid scheme if this fails.
References
----------
.. [1] https://www.gnu.org/software/gsl/manual/html_node/QNG-non_002dadaptive-Gauss_002dKronrod-integration.html
"""
def __init__(self, name='beta1d'):
self.r0 = Parameter(name, 'r0', 1, tinyval, hard_min=tinyval)
self.beta = Parameter(name, 'beta', 1, 1e-05, 10, 1e-05, 10)
self.xpos = Parameter(name, 'xpos', 0, 0, frozen=True)
self.ampl = Parameter(name, 'ampl', 1, 0)
ArithmeticModel.__init__(self, name,
(self.r0, self.beta, self.xpos, self.ampl))
[docs]
def get_center(self):
return (self.xpos.val,)
[docs]
def set_center(self, xpos, *args, **kwargs):
self.xpos.set(xpos)
[docs]
def guess(self, dep, *args, **kwargs):
pos = get_position(dep, *args)
param_apply_limits(pos, self.xpos, **kwargs)
ref = guess_reference(self.r0.min, self.r0.max, *args)
param_apply_limits(ref, self.r0, **kwargs)
norm = guess_amplitude_at_ref(self.r0.val, dep, *args)
param_apply_limits(norm, self.ampl, **kwargs)
[docs]
@modelCacher1d
def calc(self, p, *args, **kwargs):
kwargs = clean_kwargs1d(self, kwargs)
return _modelfcts.beta1d(p, *args, **kwargs)
[docs]
class BPL1D(RegriddableModel1D):
"""One-dimensional broken power-law function.
Attributes
----------
gamma1
The slope of the power law below the break point.
gamma2
The slope of the power law above the break point.
eb
The position of the break point.
ref
The reference position for the amplitude.
ampl
The amplitude, defined with respect to the reference position.
See Also
--------
Beta1D, Lorentz1D, NormBeta1D
Notes
-----
The functional form of the model for points is::
f(x) = ampl * (x / ref)^(-gamma1) x <= eb
= ampl' * (x / ref)^(-gamma2) otherwise
ampl' = ampl (eb / ref)^(gamma2 - gamma1)
and for an integrated grid it is the integral of this over
the bin.
"""
def __init__(self, name='bpl1d'):
self.gamma1 = Parameter(name, 'gamma1', 0, -10, 10)
self.gamma2 = Parameter(name, 'gamma2', 0, -10, 10)
self.eb = Parameter(name, 'eb', 100, 0, 1000, 0)
self.ref = Parameter(name, 'ref', 1, frozen=True)
self.ampl = Parameter(name, 'ampl', 1, 1e-20)
ArithmeticModel.__init__(self, name,
(self.gamma1, self.gamma2,
self.eb, self.ref, self.ampl))
[docs]
def guess(self, dep, *args, **kwargs):
pos = get_position(dep, *args)
param_apply_limits(pos, self.eb, **kwargs)
ref = guess_reference(self.ref.min, self.ref.max, *args)
param_apply_limits(ref, self.ref, **kwargs)
norm = guess_amplitude_at_ref(self.ref.val, dep, *args)
param_apply_limits(norm, self.ampl, **kwargs)
[docs]
@modelCacher1d
def calc(self, p, *args, **kwargs):
kwargs = clean_kwargs1d(self, kwargs)
return _modelfcts.bpl1d(p, *args, **kwargs)
# TODO: what are the units of the independent axis: Angstrom?
[docs]
class Dered(RegriddableModel1D):
"""A de-reddening model.
De-reddening model applied multiplicatively to a spectrum.
The integrate flag of this model should be set to False when
used with an integrated grid.
Attributes
----------
rv
The ratio of total to selective extinction.
nhgal
The absorbing column density of H_gal in units of 10^20 cm^-2.
Notes
-----
This dereddening model uses the analytic formula for the mean
extension law described in [1]_::
A(lambda) = E(B-V) * (a * rv + b)
= 1.086 tau(lambda)
where tau(lambda) is the wavelength-dependent optical depth::
I(lambda) = I(0) * exp(-tau(lambda))
and a and b are computed using wavelength-dependent formulae,
which are not reproduced here, for the wavelength range 1000
Angstroms to 3.3 microns. The relationship between the color
excess and the column density (nhgal) is [2]_::
E(B-V) = nhgal / 58.0
where the units of nhgal is 10^20 cm^-2. The value of the ratio
of total to selective extinction, rv, is initially set to 3.1, the
standard value for the diffuse ISM. The final model form is::
I(lambda) = I(0) exp(-nhgal * (a * ev + b) / (58.0 * 1.086)
This model provided courtesy of Karl Forster.
References
----------
.. [1] Cardelli, Clayton, & Mathis 1989, ApJ 345, 245
http://adsabs.harvard.edu/abs/1989ApJ...345..245C
.. [2] Bohlin, Savage, & Drake 1978, ApJ 224, 132
http://adsabs.harvard.edu/abs/1978ApJ...224..132B
"""
def __init__(self, name='dered'):
self.rv = Parameter(name, 'rv', 10, 1e-10, 1000, tinyval)
self.nhgal = Parameter(name, 'nhgal', 1e-07, 1e-07, 100000)
ArithmeticModel.__init__(self, name, (self.rv, self.nhgal))
[docs]
@modelCacher1d
def calc(self, p, *args, **kwargs):
kwargs = clean_kwargs1d(self, kwargs)
return _modelfcts.dered(p, *args, **kwargs)
[docs]
class Edge(RegriddableModel1D):
"""Photoabsorption edge model.
This model can be used when the independent axis is in energy
or wavelength space.
Attributes
----------
space
Switch to select whether the independent axis is energy or
wavelength. This parameter is not fit (``alwaysfrozen`` is
set), and should be set to either 0, when the independent axis
is energy with units of keV, or 1 when the axis is wavelength
with units of Angstrom.
thresh
The edge position (in energy or wavelength units matching
the data grid).
abs
The absorption coefficient.
Notes
-----
A phenomenological photoabsorption edge model as a function of
energy::
f(x) = exp(-abs * (x / thresh)^-3) if x >= thresh
= 1.0 otherwise
or, as a function of wavelength::
f(x) = exp(-abs * (x / thresh)^3) if x <= thresh
= 1.0 otherwise
"""
def __init__(self, name='edge'):
self.space = Parameter(name, 'space', 0, 0, 1, 0, 1,
'0 - energy | 1 - wave')
self.thresh = Parameter(name, 'thresh', 1, 0, hard_min=0)
self.abs = Parameter(name, 'abs', 1, 0)
ArithmeticModel.__init__(self, name,
(self.space, self.thresh, self.abs))
[docs]
@modelCacher1d
def calc(self, p, *args, **kwargs):
kwargs = clean_kwargs1d(self, kwargs)
return _modelfcts.edge(p, *args, **kwargs)
# DOC-NOTE:
# The equation is different to the ahelp file, but matches the code
#
# DOC-TODO:
# add note about integrated version
# convince myself that it is doing something useful.
#
[docs]
class LineBroad(RegriddableModel1D):
"""A one-dimensional line-broadening profile.
Attributes
----------
ampl
The amplitude of the line.
rest
The rest wavelength.
vsini
The rotation velocity (v sin(i)), in km/s.
Notes
-----
The model is::
f(lambda) = 2 * ampl * c * sqrt(x) / (pi * rest * vsini)
x = 1 - ((lambda - rest) * c / (rest * vsini))^2
"""
def __init__(self, name='linebroad'):
self.ampl = Parameter(name, 'ampl', 1, 0, hard_min=0)
self.rest = Parameter(name, 'rest', 1000, tinyval, hard_min=tinyval)
self.vsini = Parameter(name, 'vsini', tinyval, tinyval,
hard_min=tinyval)
ArithmeticModel.__init__(self, name,
(self.ampl, self.rest, self.vsini))
[docs]
def guess(self, dep, *args, **kwargs):
ref = guess_reference(self.rest.min, self.rest.max, *args)
param_apply_limits(ref, self.rest, **kwargs)
norm = guess_amplitude_at_ref(self.rest.val, dep, *args)
fwhm = get_fwhm(dep, *args)
c_km = 2.99792458e+5
if self.rest.val != 0:
vsini = 2. * fwhm * c_km / numpy.sqrt(3.) / self.rest.val
vs = {'val': vsini,
'min': vsini / _guess_ampl_scale,
'max': vsini * _guess_ampl_scale}
param_apply_limits(vs, self.vsini, **kwargs)
modampl = norm['val'] * numpy.pi * self.vsini.val * \
self.rest.val / 2. / c_km
mod = {'val': modampl,
'min': modampl / _guess_ampl_scale,
'max': modampl * _guess_ampl_scale}
param_apply_limits(mod, self.ampl, **kwargs)
[docs]
@modelCacher1d
def calc(self, p, *args, **kwargs):
kwargs = clean_kwargs1d(self, kwargs)
return _modelfcts.linebroad(p, *args, **kwargs)
# DOC-NOTE: for some reason the division in the equation in the notes
# section confuses sphinx (it thinks it is a section title).
#
[docs]
class Lorentz1D(RegriddableModel1D):
"""One-dimensional normalized Lorentz model function.
Attributes
----------
fwhm
The full-width half maximum of the line.
pos
The center of the line.
ampl
The amplitude refers to the integral of the model.
See Also
--------
Beta1D, NormBeta1D, NormGauss1D, Voigt1D
Notes
-----
The functional form of the model for points is::
f(x) = ampl * fwhm
--------------------------------------
2 * pi * (0.25 * fwhm^2 + (x - pos)^2)
and for an integrated grid it is the integral of this over
the bin.
The area under the function as defined above is 1 if ampl is 1.
"""
def __init__(self, name='lorentz1d'):
self.fwhm = Parameter(name, 'fwhm', 10, 0, hard_min=0)
self.pos = Parameter(name, 'pos', 1)
self.ampl = Parameter(name, 'ampl', 1)
ArithmeticModel.__init__(self, name,
(self.fwhm, self.pos, self.ampl))
[docs]
def get_center(self):
return (self.pos.val,)
[docs]
def set_center(self, pos, *args, **kwargs):
self.pos.set(pos)
[docs]
def guess(self, dep, *args, **kwargs):
pos = get_position(dep, *args)
fwhm = guess_fwhm(dep, *args)
param_apply_limits(pos, self.pos, **kwargs)
param_apply_limits(fwhm, self.fwhm, **kwargs)
norm = guess_amplitude(dep, *args)
if fwhm != 10:
aprime = norm['val'] * self.fwhm.val * numpy.pi / 2.
ampl = {'val': aprime,
'min': aprime / _guess_ampl_scale,
'max': aprime * _guess_ampl_scale}
param_apply_limits(ampl, self.ampl, **kwargs)
else:
param_apply_limits(norm, self.ampl, **kwargs)
[docs]
@modelCacher1d
def calc(self, p, *args, **kwargs):
kwargs = clean_kwargs1d(self, kwargs)
return _modelfcts.lorentz1d(p, *args, **kwargs)
[docs]
class Voigt1D(RegriddableModel1D):
"""One dimensional Voigt profile.
The Voigt profile is a convolution between a Gaussian distribution
a Cauchy-Lorentz distribution [Voigt_1912]_, [wiki_voigt]_. It is often used in
analyzing spectroscopy data.
.. versionadded:: 4.12.2
Attributes
----------
fwhm_g
The full-width half-maximum (FWHM) of the Gaussian distribution.
fwhm_l
The full-width half-maximum of the Lorentzian distribution.
pos
The center of the profile.
ampl
The amplitude of the profile.
See Also
--------
NormGauss1D, Lorentz1D, PseudoVoigt1D
Notes
-----
Following [wiki_voigt]_, the Voigt profile can be written as::
f(x) = ampl * Re[w(z)] / (sqrt(2 * PI) * sigma)
where Re[w(z)] is the real part of the Faddeeva function [wiki_faddeeva]_
and sigma and gamma are parameters of the Gaussian and
Lorentzian model respectively::
z = (x - pos + i * gamma) / (sqrt(2) * sigma)
sigma = fhwm_g / sqrt(8 * log(2))
gamma = fwhm_l / 2
One common simplification is to tie the sigma and gamma
parameters together, which can be achieved by linking the
fwhm_l parameter to fwhm_g with the following equation::
fwhm_l = fwhm_g / sqrt(2 * log(2))
An approximation for the FWHM of the profile, taken from [wiki_voigt]_,
is
0.5346 fwhm_l + sqrt(0.2166 fwhm_l^2 + fwhm_g^2)
References
----------
.. [Voigt_1912] https://publikationen.badw.de/de/003395768
.. [wiki_voigt] https://en.wikipedia.org/wiki/Voigt_profile
.. [wiki_faddeeva] https://en.wikipedia.org/wiki/Faddeeva_function
Examples
--------
Force the widths of the Gaussian and Lorentzian components
to be the same:
>>> mdl = Voigt1D()
>>> mdl.fwhm_l = mdl.fwhm_g / np.sqrt(2 * np.log(2))
"""
def __init__(self, name='voigt1d'):
self.fwhm_g = Parameter(name, 'fwhm_g', 10, tinyval, hard_min=tinyval)
self.fwhm_l = Parameter(name, 'fwhm_l', 10, 0, hard_min=0)
self.pos = Parameter(name, 'pos', 0.0)
self.ampl = Parameter(name, 'ampl', 1.0)
ArithmeticModel.__init__(self, name,
(self.fwhm_g, self.fwhm_l, self.pos, self.ampl))
return
[docs]
def get_center(self):
return (self.pos.val,)
[docs]
def set_center(self, pos, *args, **kwargs):
self.pos.set(pos)
[docs]
def guess(self, dep, *args, **kwargs):
"""Guess the parameter values.
The fwhm_g and fwhm_l parameters are set to the
same value.
"""
pos = get_position(dep, *args)
fwhm = guess_fwhm(dep, *args)
param_apply_limits(pos, self.pos, **kwargs)
param_apply_limits(fwhm, self.fwhm_g, **kwargs)
param_apply_limits(fwhm, self.fwhm_l, **kwargs)
# I am using an approximate conversion to get the
# amplitude from guess_amplitude. The conversion
# between height (as returned - roughly - by guess_amplitude)
# and amplitude is, from lmfit
# amplitude = height * (max(2.220446049250313e-16, sigma*sqrt(2*pi))))*wofz((1j*gamma)/(max(2.220446049250313e-16, sigma*sqrt(2))).real'
# but I'm just going to use the lorentz1d amplitude guess
#
norm = guess_amplitude(dep, *args)
aprime = norm['val'] * fwhm['val'] * numpy.pi / 2.
ampl = {'val': aprime,
'min': aprime / _guess_ampl_scale,
'max': aprime * _guess_ampl_scale}
param_apply_limits(ampl, self.ampl, **kwargs)
[docs]
@modelCacher1d
def calc(self, p, *args, **kwargs):
kwargs = clean_kwargs1d(self, kwargs)
return _modelfcts.wofz(p, *args, **kwargs)
[docs]
class PseudoVoigt1D(RegriddableModel1D):
"""A weighted sum of a Gaussian and Lorentzian distribution.
Unlike the Voigt1D model, which is a convolution between a
Gaussian and Lorentz distribution, this approximates the
Voigt profile with a linear combination of the two profiles [1]_.
It is often used in spectroscopy.
.. versionadded:: 4.12.2
Attributes
----------
frac
The fraction of the model composed of the Gaussian profile
(0 to 1).
fwhm
The full-width half-maximum (FWHM) of each component.
pos
The center of the profile.
ampl
The amplitude of the profile.
See Also
--------
NormGauss1D, Lorentz1D, Voigt1D
Notes
-----
The model can be written as::
f(x) = frac * g(x) + (1 - frac) * l(x)
where g(x) and l(x) are NormGauss1D and Lorentz1D models with
the fwhm, pos, and ampl values taken from this model.
References
----------
.. [1] https://en.wikipedia.org/wiki/Voigt_profile#Pseudo-Voigt_approximation
"""
def __init__(self, name='pseudovoigt1d'):
self.frac = Parameter(name, 'frac', 0.5, 0, 1, hard_min=0, hard_max=1)
self.fwhm = Parameter(name, 'fwhm', 10, tinyval, hard_min=tinyval)
self.pos = Parameter(name, 'pos', 0.0)
self.ampl = Parameter(name, 'ampl', 1.0)
ArithmeticModel.__init__(self, name,
(self.frac, self.fwhm, self.pos, self.ampl))
return
[docs]
def get_center(self):
return (self.pos.val,)
[docs]
def set_center(self, pos, *args, **kwargs):
self.pos.set(pos)
[docs]
def guess(self, dep, *args, **kwargs):
"""Guess the parameter values.
The frac parameter is set to 0.5.
"""
self.frac.set(0.5, min=0, max=1)
pos = get_position(dep, *args)
fwhm = guess_fwhm(dep, *args)
param_apply_limits(pos, self.pos, **kwargs)
param_apply_limits(fwhm, self.fwhm, **kwargs)
# See Voigt1D for the amplitude guess
# amplitude from guess_amplitude. although I've
# added another factor of 1/2
#
norm = guess_amplitude(dep, *args)
aprime = norm['val'] * fwhm['val'] * numpy.pi / 4.
ampl = {'val': aprime,
'min': aprime / _guess_ampl_scale,
'max': aprime * _guess_ampl_scale}
param_apply_limits(ampl, self.ampl, **kwargs)
[docs]
@modelCacher1d
def calc(self, p, *args, **kwargs):
kwargs = clean_kwargs1d(self, kwargs)
frac = p[0]
gmdl = sherpa.models._modelfcts.ngauss1d(p[1:], *args, **kwargs)
lmdl = _modelfcts.lorentz1d(p[1:], *args, **kwargs)
return frac * gmdl + (1 - frac) * lmdl
[docs]
class NormBeta1D(RegriddableModel1D):
"""One-dimensional normalized beta model function.
This is the same model as the ``Beta1D`` model but with a
different slope parameter and normalisation.
Attributes
----------
pos
The center of the line.
w
The line width.
alpha
The slope of the profile at large radii.
ampl
The amplitude refers to the integral of the model.
See Also
--------
Beta1D, Lorentz1D
Notes
-----
The functional form of the model for points is::
f(x) = A * (1 + ((x - pos) / w)^2)^(-alpha)
A = ampl / integral f(x) dx
The grid version is evaluated by numerically intgerating the
function over each bin using a non-adaptive Gauss-Kronrod scheme
suited for smooth functions [1]_, falling over to a simple
trapezoid scheme if this fails.
References
----------
.. [1] https://www.gnu.org/software/gsl/manual/html_node/QNG-non_002dadaptive-Gauss_002dKronrod-integration.html
"""
def __init__(self, name='normbeta1d'):
self.pos = Parameter(name, 'pos', 0)
self.width = Parameter(name, 'width', 1, tinyval, hard_min=tinyval)
self.index = Parameter(name, 'index', 2.5, 0.5, 1000, 0.5)
self.ampl = Parameter(name, 'ampl', 1, 0)
ArithmeticModel.__init__(self, name,
(self.pos, self.width, self.index, self.ampl))
[docs]
def get_center(self):
return (self.pos.val,)
[docs]
def set_center(self, pos, *args, **kwargs):
self.pos.set(pos)
[docs]
def guess(self, dep, *args, **kwargs):
ampl = guess_amplitude(dep, *args)
pos = get_position(dep, *args)
fwhm = guess_fwhm(dep, *args)
param_apply_limits(pos, self.pos, **kwargs)
norm = fwhm['val'] * numpy.sqrt(numpy.pi) * \
numpy.exp(lgam(self.index.val - 0.5) - lgam(self.index.val))
for key in ampl.keys():
ampl[key] *= norm
param_apply_limits(ampl, self.ampl, **kwargs)
[docs]
@modelCacher1d
def calc(self, p, *args, **kwargs):
kwargs = clean_kwargs1d(self, kwargs)
return _modelfcts.nbeta1d(p, *args, **kwargs)
[docs]
class Schechter(RegriddableModel1D):
"""One-dimensional Schechter model function.
This model is for integrated data grids only.
Attributes
----------
alpha
The slope of the power-law component.
ref
The reference position, which controls the switch from the
power-law to the exponential.
norm
The normalisation of the model.
Notes
-----
The functional form of the model for grids is::
f(xlo,xhi) = norm * (xlo / ref)^alpha
* exp(-xlo / ref)
* (xhi - xlo) / ref
and for points the model is::
f(x) = 0
"""
def __init__(self, name='schechter'):
self.alpha = Parameter(name, 'alpha', 10)
self.ref = Parameter(name, 'ref', 1)
self.norm = Parameter(name, 'norm', 1)
ArithmeticModel.__init__(self, name, (self.alpha, self.ref, self.norm))
[docs]
def guess(self, dep, *args, **kwargs):
norm = guess_amplitude(dep, *args)
param_apply_limits(norm, self.norm, **kwargs)
[docs]
@modelCacher1d
def calc(self, p, *args, **kwargs):
if not self.integrate:
raise ModelErr('alwaysint', self.name)
kwargs = clean_kwargs1d(self, kwargs)
return _modelfcts.schechter(p, *args, **kwargs)
[docs]
class Beta2D(RegriddableModel2D):
"""Two-dimensional beta model function.
The beta model is a Lorentz model with a varying power law.
Attributes
----------
r0
The core radius.
xpos
X0 axis coordinate of the model center (position of the peak).
ypos
X1 axis coordinate of the model center (position of the peak).
ellip
The ellipticity of the model.
theta
The angle of the major axis. It is in radians, measured
counter-clockwise from the X0 axis (i.e. the line X1=0).
ampl
The model value at the peak position (xpos, ypos).
alpha
The power-law slope of the profile at large radii.
See Also
--------
Beta1D, DeVaucouleurs2D, HubbleReynolds, Lorentz2D, Sersic2D
Notes
-----
The functional form of the model for points is::
f(x0,x1) = ampl * (1 + r(x0,x1)^2)^(-alpha)
r(x0,x1)^2 = xoff(x0,x1)^2 * (1-ellip)^2 + yoff(x0,x1)^2
-------------------------------------------
r0^2 * (1-ellip)^2
xoff(x0,x1) = (x0 - xpos) * cos(theta) + (x1 - ypos) * sin(theta)
yoff(x0,x1) = (x1 - ypos) * cos(theta) - (x0 - xpos) * sin(theta)
The grid version is evaluated by adaptive multidimensional
integration scheme on hypercubes using cubature rules, based
on code from HIntLib ([1]_) and GSL ([2]_).
References
----------
.. [1] HIntLib - High-dimensional Integration Library
http://mint.sbg.ac.at/HIntLib/
.. [2] GSL - GNU Scientific Library
http://www.gnu.org/software/gsl/
"""
def __init__(self, name='beta2d'):
self.r0 = Parameter(name, 'r0', 10, tinyval, hard_min=tinyval)
self.xpos = Parameter(name, 'xpos', 0)
self.ypos = Parameter(name, 'ypos', 0)
self.ellip = Parameter(name, 'ellip', 0, 0, 0.999, 0, 0.9999,
frozen=True)
self.theta = Parameter(name, 'theta', 0, -2*numpy.pi, 2*numpy.pi, -2*numpy.pi,
4*numpy.pi, 'radians', True)
self.ampl = Parameter(name, 'ampl', 1)
self.alpha = Parameter(name, 'alpha', 1, -10, 10)
ArithmeticModel.__init__(self, name,
(self.r0, self.xpos, self.ypos, self.ellip,
self.theta, self.ampl, self.alpha))
self.cache = 0
[docs]
def get_center(self):
return (self.xpos.val, self.ypos.val)
[docs]
def set_center(self, xpos, ypos, *args, **kwargs):
self.xpos.set(xpos)
self.ypos.set(ypos)
[docs]
def guess(self, dep, *args, **kwargs):
xpos, ypos = guess_position(dep, *args)
norm = guess_amplitude2d(dep, *args)
rad = guess_radius(*args)
param_apply_limits(xpos, self.xpos, **kwargs)
param_apply_limits(ypos, self.ypos, **kwargs)
param_apply_limits(norm, self.ampl, **kwargs)
param_apply_limits(rad, self.r0, **kwargs)
[docs]
def calc(self, p, *args, **kwargs):
kwargs = clean_kwargs2d(self, kwargs)
return _modelfcts.beta2d(p, *args, **kwargs)
[docs]
class DeVaucouleurs2D(RegriddableModel2D):
"""Two-dimensional de Vaucouleurs model.
A spatial de Vaucouleurs profile which is a formulation of the
R^(1/4) law introduced by [1]_. It is a special case of the
``Sersic2D`` model with ``n=4``, as described in [2]_, [3]_,
and [4]_.
Attributes
----------
r0
The core radius.
xpos
The center of the model on the x0 axis.
ypos
The center of the model on the x1 axis.
ellip
The ellipticity of the model.
theta
The angle of the major axis. It is in radians, measured
counter-clockwise from the X0 axis (i.e. the line X1=0).
ampl
The amplitude refers to the maximum peak of the model.
See Also
--------
Beta2D, HubbleReynolds, Lorentz2D, Sersic2D
Notes
-----
The model used is the same as the ``Sersic2D`` model with ``n=4``.
References
----------
.. [1] de Vaucouleurs G., 1948, Ann. d'Astroph. 11, 247
http://adsabs.harvard.edu/abs/1948AnAp...11..247D
.. [2] http://ned.ipac.caltech.edu/level5/March05/Graham/Graham2.html
.. [3] Graham, A. & Driver, S., 2005, PASA, 22, 118
http://adsabs.harvard.edu/abs/2005PASA...22..118G
.. [4] Ciotti, L. & Bertin, G., A&A, 1999, 352, 447-451
http://adsabs.harvard.edu/abs/1999A%26A...352..447C
"""
def __init__(self, name='devaucouleurs2d'):
self.r0 = Parameter(name, 'r0', 10, 0, hard_min=0)
self.xpos = Parameter(name, 'xpos', 0)
self.ypos = Parameter(name, 'ypos', 0)
self.ellip = Parameter(name, 'ellip', 0, 0, 0.999, 0, 0.9999)
self.theta = Parameter(name, 'theta', 0, -2*numpy.pi, 2*numpy.pi, -2*numpy.pi,
4*numpy.pi, 'radians')
self.ampl = Parameter(name, 'ampl', 1)
ArithmeticModel.__init__(self, name,
(self.r0, self.xpos, self.ypos, self.ellip,
self.theta, self.ampl))
self.cache = 0
[docs]
def get_center(self):
return (self.xpos.val, self.ypos.val)
[docs]
def set_center(self, xpos, ypos, *args, **kwargs):
self.xpos.set(xpos)
self.ypos.set(ypos)
[docs]
def guess(self, dep, *args, **kwargs):
xpos, ypos = guess_position(dep, *args)
norm = guess_amplitude2d(dep, *args)
rad = guess_radius(*args)
param_apply_limits(xpos, self.xpos, **kwargs)
param_apply_limits(ypos, self.ypos, **kwargs)
param_apply_limits(norm, self.ampl, **kwargs)
param_apply_limits(rad, self.r0, **kwargs)
[docs]
def calc(self, p, *args, **kwargs):
kwargs = clean_kwargs2d(self, kwargs)
return _modelfcts.devau(p, *args, **kwargs)
[docs]
class HubbleReynolds(RegriddableModel2D):
"""Two-dimensional Hubble-Reynolds model.
Attributes
----------
r0
The core radius.
xpos
The center of the model on the x0 axis.
ypos
The center of the model on the x1 axis.
ellip
The ellipticity of the model.
theta
The angle of the major axis. It is in radians, measured
counter-clockwise from the X0 axis (i.e. the line X1=0).
ampl
The amplitude refers to the maximum peak of the model.
See Also
--------
Beta2D, DeVaucouleurs2D, Lorentz2D, Sersic2D
Notes
-----
The functional form of the model for points is::
f(x0,x1) = ampl / (1 + r(x0,x1))^2
r(x0,x1)^2 = xoff(x0,x1)^2 * (1-ellip)^2 + yoff(x0,x1)^2
-------------------------------------------
r0^2 * (1-ellip)^2
xoff(x0,x1) = (x0 - xpos) * cos(theta) + (x1 - ypos) * sin(theta)
yoff(x0,x1) = (x1 - ypos) * cos(theta) - (x0 - xpos) * sin(theta)
The grid version is evaluated by adaptive multidimensional
integration scheme on hypercubes using cubature rules, based
on code from HIntLib ([1]_) and GSL ([2]_).
References
----------
.. [1] HIntLib - High-dimensional Integration Library
http://mint.sbg.ac.at/HIntLib/
.. [2] GSL - GNU Scientific Library
http://www.gnu.org/software/gsl/
"""
def __init__(self, name='hubblereynolds'):
self.r0 = Parameter(name, 'r0', 10, 0, hard_min=0)
self.xpos = Parameter(name, 'xpos', 0)
self.ypos = Parameter(name, 'ypos', 0)
self.ellip = Parameter(name, 'ellip', 0, 0, 0.999, 0, 0.9999)
self.theta = Parameter(name, 'theta', 0, -2*numpy.pi, 2*numpy.pi, -2*numpy.pi,
4*numpy.pi, 'radians')
self.ampl = Parameter(name, 'ampl', 1)
ArithmeticModel.__init__(self, name,
(self.r0, self.xpos, self.ypos, self.ellip,
self.theta, self.ampl))
self.cache = 0
[docs]
def get_center(self):
return (self.xpos.val, self.ypos.val)
[docs]
def set_center(self, xpos, ypos, *args, **kwargs):
self.xpos.set(xpos)
self.ypos.set(ypos)
[docs]
def guess(self, dep, *args, **kwargs):
xpos, ypos = guess_position(dep, *args)
norm = guess_amplitude2d(dep, *args)
rad = guess_radius(*args)
param_apply_limits(xpos, self.xpos, **kwargs)
param_apply_limits(ypos, self.ypos, **kwargs)
param_apply_limits(norm, self.ampl, **kwargs)
param_apply_limits(rad, self.r0, **kwargs)
[docs]
def calc(self, p, *args, **kwargs):
kwargs = clean_kwargs2d(self, kwargs)
return _modelfcts.hr(p, *args, **kwargs)
[docs]
class Lorentz2D(RegriddableModel2D):
"""Two-dimensional un-normalised Lorentz function.
Attributes
----------
fwhm
The full-width half maximum.
xpos
The center of the model on the x0 axis.
ypos
The center of the model on the x1 axis.
ellip
The ellipticity of the model.
theta
The angle of the major axis. It is in radians, measured
counter-clockwise from the X0 axis (i.e. the line X1=0).
ampl
The amplitude refers to the maximum peak of the model.
See Also
--------
Beta1D, DeVaucouleurs2D, HubbleReynolds, Lorentz1D, Sersic2D
Notes
-----
The functional form of the model for points is::
f(x0,x1) = ampl / (1 + 4 * r(x0,x1)^2)
r(x0,x1)^2 = xoff(x0,x1)^2 * (1-ellip)^2 + yoff(x0,x1)^2
-------------------------------------------
fwhm^2 * (1-ellip)^2
xoff(x0,x1) = (x0 - xpos) * cos(theta) + (x1 - ypos) * sin(theta)
yoff(x0,x1) = (x1 - ypos) * cos(theta) - (x0 - xpos) * sin(theta)
and for an integrated grid it is the integral of this over
the bin.
"""
def __init__(self, name='lorentz2d'):
self.fwhm = Parameter(name, 'fwhm', 10, tinyval, hard_min=tinyval)
self.xpos = Parameter(name, 'xpos', 0)
self.ypos = Parameter(name, 'ypos', 0)
self.ellip = Parameter(name, 'ellip', 0, 0, 0.999, 0, 0.9999,
frozen=True)
self.theta = Parameter(name, 'theta', 0, -2*numpy.pi, 2*numpy.pi,
-2*numpy.pi, 4*numpy.pi, 'radians', frozen=True)
self.ampl = Parameter(name, 'ampl', 1)
ArithmeticModel.__init__(self, name,
(self.fwhm, self.xpos, self.ypos, self.ellip,
self.theta, self.ampl))
self.cache = 0
[docs]
def get_center(self):
return (self.xpos.val, self.ypos.val)
[docs]
def set_center(self, xpos, ypos, *args, **kwargs):
self.xpos.set(xpos)
self.ypos.set(ypos)
[docs]
def guess(self, dep, *args, **kwargs):
xpos, ypos = guess_position(dep, *args)
norm = guess_amplitude2d(dep, *args)
param_apply_limits(xpos, self.xpos, **kwargs)
param_apply_limits(ypos, self.ypos, **kwargs)
param_apply_limits(norm, self.ampl, **kwargs)
[docs]
def calc(self, p, *args, **kwargs):
kwargs = clean_kwargs2d(self, kwargs)
return _modelfcts.lorentz2d(p, *args, **kwargs)
[docs]
class JDPileup(RegriddableModel1D):
"""A CCD pileup model for the ACIS detectors on Chandra.
This model is based on the work by John Davis ([1]_). It is
intended only for modeling the pileup in one-dimensional spectra
obtained in imaging mode (i.e. no grating data), but can be used
with the zeroth-order spectrum of a grating data set.
Attributes
----------
alpha
The alpha parameter parameterizes "grade migration" in the
detector and represents the fraction of piled-up events that
result in a good grade.
g0
The probability of assigning a grade of zero. This should
remain frozen.
f
The fraction of flux falling into the pileup region. This
should remain frozen.
n
The number of detection cells. This parameter can not be fit.
ftime
The frame time in seconds (as given by the ``EXPTIME``
keyword of the event file). This parameter can not be fit.
fracexp
The fractional exposure that the source experienced while
dithering on the chip (as given by the ``FRACEXPO`` keyword
in the ARF file). This parameter can not be fit.
nterms
The maximum number of photons considered for pileup in a
single readout frame. This should not be changed from its
default value of 30. This parameter can not be fit.
Notes
-----
As the pileup model is inherently non-linear, it is strongly
advised that multiple optimization methods are used to thoroughly
investigate the search space for the model.
The alpha parameter should vary with photon energy and detector
position, but for simplicity it is treated as independent of
energy and location.
An example of using this model to fit a Chandra spectrum is
provided in [2]_.
References
----------
.. [1] Davis, J, 2001, ApJ, 562, 575-582.
http://adsabs.harvard.edu/abs/2001ApJ...562..575D
.. [2] "Using A Pileup Model",
http://cxc.harvard.edu/sherpa/threads/pileup/
"""
def __init__(self, name='jdpileup'):
self.alpha = Parameter(name, 'alpha', 0.5, 0, 1, 0, 1)
self.g0 = Parameter(name, 'g0', 1, tinyval, 1, tinyval, 1, frozen=True)
self.f = Parameter(name, 'f', 0.95, 0.9, 1, 0, 1)
self.n = Parameter(name, 'n', 1, tinyval, 100, tinyval, 2048,
alwaysfrozen=True)
self.ftime = Parameter(name, 'ftime', 3.241, tinyval, 5, tinyval, 100,
'sec', alwaysfrozen=True)
self.fracexp = Parameter(name, 'fracexp', 0.987, 0, 1, 0, 1,
alwaysfrozen=True)
self.nterms = Parameter(name, 'nterms', 30, 1, 100, 1, 100,
alwaysfrozen=True)
self._results = None
ArithmeticModel.__init__(self, name,
(self.alpha, self.g0, self.f, self.n,
self.ftime, self.fracexp, self.nterms))
def __str__(self):
s = ArithmeticModel.__str__(self)
if self._results is not None:
pileup_fractions, integral_ae, num_terms = self._results
sum = 0.0
if num_terms > 0:
sum = pileup_fractions[1]
sum_piled = pileup_fractions[2:num_terms + 1].sum()
sum += sum_piled
pn = numpy.exp(-integral_ae)
s += '\n\n'
for i in range(1, num_terms + 1):
pn *= integral_ae / float(i)
s += ' %d: %g %g\n' % (i, pn,
pileup_fractions[i] / sum)
s += ' *** pileup fraction: %g' % (sum_piled / sum)
return s
[docs]
def calc(self, p, arf_source, exposure_time, min_energy, max_energy,
specresp, model):
(alpha, g0, psf_frac, num_regions, frame_time, fracexpo,
max_num_terms) = p
out = apply_pileup(arf_source, exposure_time, int(max_num_terms),
min_energy, max_energy, specresp, fracexpo,
frame_time, alpha, g0, num_regions, psf_frac, model)
self._results = out[1:]
return out[0]
[docs]
class Sersic2D(RegriddableModel2D):
"""Two-dimensional Sersic model.
This is a generalization of the ``DeVaucouleurs2D`` model,
in which the exponent ``n`` can vary ([1]_, [2]_, and [3]_).
Attributes
----------
r0
The core radius.
xpos
The center of the model on the x0 axis.
ypos
The center of the model on the x1 axis.
ellip
The ellipticity of the model.
theta
The angle of the major axis. It is in radians, measured
counter-clockwise from the X0 axis (i.e. the line X1=0).
ampl
The amplitude refers to the maximum peak of the model.
n
The Sersic index (n=4 replicates the ``DeVaucouleurs2D``
model).
See Also
--------
Beta2D, DeVaucouleurs2D, HubbleReynolds, Lorentz2D
Notes
-----
The functional form of the model for points is can be
expressed as the following::
f(x0,x1) = ampl * exp(-b(n) * (r(x0,x1)^(1/n) - 1))
b(n) = 2 * n - 1 / 3 + 4 / (405 * n) + 46 / (25515 * n^2)
r(x0,x1)^2 = xoff(x0,x1)^2 * (1-ellip)^2 + yoff(x0,x1)^2
-------------------------------------------
r0^2 * (1-ellip)^2
xoff(x0,x1) = (x0 - xpos) * cos(theta) + (x1 - ypos) * sin(theta)
yoff(x0,x1) = (x1 - ypos) * cos(theta) - (x0 - xpos) * sin(theta)
The grid version is evaluated by adaptive multidimensional
integration scheme on hypercubes using cubature rules, based
on code from HIntLib ([4]_) and GSL ([5]_).
References
----------
.. [1] http://ned.ipac.caltech.edu/level5/March05/Graham/Graham2.html
.. [2] Graham, A. & Driver, S., 2005, PASA, 22, 118
http://adsabs.harvard.edu/abs/2005PASA...22..118G
.. [3] Ciotti, L. & Bertin, G., A&A, 1999, 352, 447-451
http://adsabs.harvard.edu/abs/1999A%26A...352..447C
.. [4] HIntLib - High-dimensional Integration Library
http://mint.sbg.ac.at/HIntLib/
.. [5] GSL - GNU Scientific Library
http://www.gnu.org/software/gsl/
"""
def __init__(self, name='sersic2d'):
self.r0 = Parameter(name, 'r0', 10, 0, hard_min=0)
self.xpos = Parameter(name, 'xpos', 0)
self.ypos = Parameter(name, 'ypos', 0)
self.ellip = Parameter(name, 'ellip', 0, 0, 0.999, 0, 0.9999)
self.theta = Parameter(name, 'theta', 0, -2*numpy.pi, 2*numpy.pi, -2*numpy.pi,
4*numpy.pi, 'radians')
self.ampl = Parameter(name, 'ampl', 1)
self.n = Parameter(name, 'n', 1, .1, 10, 0.01, 100, frozen=True)
ArithmeticModel.__init__(self, name,
(self.r0, self.xpos, self.ypos, self.ellip,
self.theta, self.ampl, self.n))
self.cache = 0
[docs]
def get_center(self):
return (self.xpos.val, self.ypos.val)
[docs]
def set_center(self, xpos, ypos, *args, **kwargs):
self.xpos.set(xpos)
self.ypos.set(ypos)
[docs]
def guess(self, dep, *args, **kwargs):
xpos, ypos = guess_position(dep, *args)
norm = guess_amplitude2d(dep, *args)
rad = guess_radius(*args)
param_apply_limits(xpos, self.xpos, **kwargs)
param_apply_limits(ypos, self.ypos, **kwargs)
param_apply_limits(norm, self.ampl, **kwargs)
param_apply_limits(rad, self.r0, **kwargs)
[docs]
def calc(self, p, *args, **kwargs):
kwargs = clean_kwargs2d(self, kwargs)
return _modelfcts.sersic(p, *args, **kwargs)
# ## disk2d and shell2d models
# ## Contributed by Christoph Deil of the HESS project
# ## Added to CIAO 4.6 for Dec. 2013 release SMD
[docs]
class Disk2D(RegriddableModel2D):
"""Two-dimensional uniform disk model.
Two-dimensional step function model consisting of a uniform
intensity disk.
Attributes
----------
xpos
The center of the disk on the x0 axis.
ypos
The center of the disk on the x1 axis.
ampl
The amplitude of the signal within the disk.
r0
The radius of the disk.
See Also
--------
Shell2D
Notes
-----
The functional form of the model for points is::
f(x0,x1) = ampl if (x0 - xpos)^2 + (x1 - ypos)^2 <= r0^2
= 0 otherwise
For grids, the x0lo,x1lo values are used in the above equation.
This model was provided by Christoph Deil.
"""
def __init__(self, name='disk2d'):
self.xpos = Parameter(name, 'xpos', 0) # p[0]
self.ypos = Parameter(name, 'ypos', 0) # p[1]
self.ampl = Parameter(name, 'ampl', 1) # p[2]
self.r0 = Parameter(name, 'r0', 1, 0) # p[3]
ArithmeticModel.__init__(self, name, (self.xpos, self.ypos, self.ampl, self.r0))
[docs]
def calc(self, p, x, y, *args, **kwargs):
# Compute radii
r2 = (x - p[0]) ** 2 + (y - p[1]) ** 2
# Return ampl when r2 <= r0 else return 0
return numpy.select([r2 <= p[3] ** 2], [p[2]])
# DOC-NOTE: TODO finish the functional form description
[docs]
class Shell2D(RegriddableModel2D):
"""A homogeneous spherical 3D shell projected onto 2D.
Attributes
----------
xpos
The center of the shell on the x0 axis.
ypos
The center of the shell on the x1 axis.
ampl
The amplitude.
r0
The line-of-sight distance.
width
The width of the shell.
See Also
--------
Disk2D
Notes
-----
The functional form of the model for points is::
f(x0,x1) = ampl * r(x0,x1)
= 0 otherwise
r(x0,x1) = ?
For grids, the x0lo,x1lo values are used in the above equation.
It is not correct for very large shells on the sky.
This model was provided by Christoph Deil.
"""
def __init__(self, name='shell2d'):
self.xpos = Parameter(name, 'xpos', 0) # p[0]
self.ypos = Parameter(name, 'ypos', 0) # p[1]
self.ampl = Parameter(name, 'ampl', 1) # p[2]
self.r0 = Parameter(name, 'r0', 1, 0) # p[3]
self.width = Parameter(name, 'width', 0.1, 0)
ArithmeticModel.__init__(self, name, (self.xpos, self.ypos, self.ampl, self.r0, self.width))
[docs]
def calc(self, p, x, y, *args, **kwargs):
"""Homogeneously emitting spherical shell,
projected along the z-direction
(this is not 100% correct for very large shells on the sky).
"""
(xpos, ypos, ampl, r_0, width) = p
r2 = (x - xpos) * (x - xpos) + (y - ypos) * (y - ypos)
r_out = r_0 + width
r_in2, r_out2 = r_0 * r_0, r_out * r_out
# r_in3, r_out3 = r_in * r_in2, r_out * r_out2
# We only call abs() in sqrt() to avoid warning messages.
sphere_out = numpy.sqrt(numpy.abs(r_out2 - r2))
sphere_in = numpy.sqrt(numpy.abs(r_in2 - r2))
# Note: for r > r_out 'numpy.select' fills automatically zeros!
non_normalized = numpy.select([r2 <= r_in2, r2 <= r_out2],
[sphere_out - sphere_in, sphere_out])
# Computed with Mathematica:
integral = 2 * numpy.pi / 3 * (r_out ** 3 - r_0 ** 3)
# integral = 1
return ampl * integral * non_normalized