Welcome to concord’s documentation!
concord
implements some commonly-used approaches for determining system
parameters for neutron stars from thermonuclear bursts,
and fully accounting for astrophysical uncertainties.
The code includes the capability for quantitative comparison of observed burst
properties with the predictions of numerical models,
and the tools presented here are intended
to make comprehensive model-observation comparisons straightforward.
This code is under active development, but the v1.0.0 release is associated with a companion paper accepted by Astrophysical Journal Supplements (see Galloway et al. 2022, also available at arXiv:2210.03598). A preprint of the paper is available in the doc subdirectory of the repository
Getting started
Clone the repository and install using pip:
python3 -m pip install
You can then import the repository and use the functions. Here’s a very simple example, to find the peak luminosity of a burst from 4U 0513+40 measured by RXTE, as part of the MINBAR sample. The first part calculates the isotropic luminosity, neglecting the uncertainty in both the peak flux and the distance:
import concord as cd
import astropy.units as u
F_pk, e_F_pk = 21.72, 0.6 # 1E-9 erg/cm^2/s bolometric; MINBAR #3443
d = (10.32, 0.24, 0.20) # asymmetric errors from Watkins et al. 2015
l_iso = cd.luminosity( F_pk, dist=d[0], isotropic=True )
print (l_iso)
# 2.767771097997098e+38 erg / s
The second part takes into account both the uncertainties in the peak flux and distance (including the asymmetric errors), and also includes the model-predicted effect of the high system inclination (>80 degrees):
l_asym = cd.luminosity( (F_pk, e_F_pk), dist=d, burst=True, imin=80, imax=90, fulldist=True)
lc = l_asym['lum'].pdf_percentiles([50, 50 - cd.CONF / 2, 50 + cd.CONF / 2])
l_unit = 1e38*u.erg/u.s
print ('''\nIsotropic luminosity is {:.2f}e38 erg/s
Taking into account anisotropy, ({:.2f}-{:.2f}+{:.2f})e38 erg/s'''.format(l_iso/l_unit, lc[0]/l_unit, (lc[0]-lc[1])/l_unit, (lc[2]-lc[0])/l_unit))
# Isotropic luminosity is 2.77e38 erg/s
# Taking into account anisotropy, (4.85-0.43+0.51)e38 erg/s
With the fulldist=True option, the function returns a dictionary with the Monte-Carlo generated distribution of the result (key lum) and all the intermediate quantities.
Check the Inferring burster properties
notebook, which includes the
examples used in the paper, for additional demonstrations of usage.
Obtaining burst data
You’ll also need some burst measurements or data; at their simplest, most of the functions operate on measurable quantities (with uncertainties), like the peak burst flux, burst fluence, persistent flux etc.
You can access data from over 7000 bursts from 85 sources via MINBAR, by also installing the Python repository from https://github.com/outs1der/minbar.
You can also use your own time-resolved spectroscopic data to generate an
concord.burstclass.ObservedBurst
object, and then run the
functions on that. You might need to define routines to read in your own
data, depending on the format.
You can download some reference bursts from
http://burst.sci.monash.edu/reference, and put the data (including the
.dat
files and the table2.tex
file) in the concord/data
directory.
You can use the example concord.burstclass.KeplerBurst
class
to read in some model bursts from
http://burst.sci.monash.edu/kepler.
Alternatively you can use that class as an example to write your own appropriate for your model results.
How the functions work
There are three principal components; the functions in utils.py
, the
anisotropy treatment in diskmodel.py
and the observed and model
burst classes in burstclass.py
.
The functions in utils.py
provide the basic functionality as
described in Functions. In order to treat the uncertainties in
the observed quantities, we adopt a Monte-Carlo (MC) approach via the
astropy
Distribution package.
Input values (including measured
quantities) can be provided in four different ways:
scalar
value with symmetric error, represented as a tuple
(value, error)
value with asymmetric error, as a 3-element tuple
(value, lo_error, hi_error)
an arbitrary array.
Quantities with error estimates will be converted to a Distribution
object of user-defined size, and with a symmetric or asymmetric normal
distribution. In this way we can represent a wide range of probability
distribution functions (PDFs) for the input parameters. These
Distribution
objects can then be used in calculations as for scalars,
hence providing uncertainty propagation at the cost of additional
computation.
Input values can be provided with units, or if units are absent, will be assumed to have the standard units for MINBAR quantities; for flux, \(10^{-9}\ \rm{erg\,cm^{-2}\,s^{-1}}\), burst fluence \(10^{-6}\ \rm{erg\,cm^{-2}}\), recurrence time hr, and so on.
All functions can operate assuming isotropic distributions for the burst
or persistent emission (isotropic=True
), but by default will include
the possible effects of anisotropic emission via the diskmodel.py
routine. This routine incorporates the modelling of He & Keek (2016, ApJ
819, #47), via ASCII
tables provided with the code.
There are three options for the treatment of emission anisotropy:
isotropic=True
- assumes both the burst and persistent emission is isotropic (\(\xi_b=\xi_p=1\))isotropic=False
- incorporates anisotropy, by drawing a set of inclination values and calculating the H-fraction for each value. Inclinations are uniform on the sphere (up to a maximum value of \(i=72^\circ\), based on the absence of dips) and the value returned is the mean and 1-sigma error, or the complete distributions (with the fulldist option)inclination=<i>
- calculate for a specific value of inclination
Where observed or model burst lightcurves are available, the classes in
burstclass.py
provide a way to represent those observations and
perform various standard analyses, including observation-model
comparisons. The ObservedBurst
class can be instantiated from an
ASCII file giving the burst flux as a function of time, or in a number of
other ways.
An example KeplerBurst
model burst class is provided, which offers a
number of ways of reading in Kepler burst runs, and which can be
adapted to outputs of different codes.
The reference bursts provided by Galloway et al. (2017, PASA 34, e109) can be read
in directly, provided the data is included in the data
directory.
Functions
- concord.utils.L_Edd(F_pk, dist=<Quantity 8. kpc>, nsamp=1000, isotropic=False, inclination=None, imin=0.0, imax=75.0, dip=False, conf=68.0, fulldist=False)
This routine estimates the Eddington luminosity of the source, based on the measured peak flux of a radius-expansion burst and the distance (or some distribution thereof. The default isotropic=False option also takes into account the likely effect of anisotropic burst emission, based on the models provided by concord
This routine has been supplanted by
concord.utils.luminosity()
, which it now calls; the only thing not incorporated into the new routine is the “simple” estimate of the isotropic luminosity error,L_Edd_iso_err = L_Edd_iso * np.sqrt((F_pk_err / F_pk) ** 2 + (2. * dist_err / dist) ** 2)
Burst emission is assumed (
burst=True
); apart from the F_pk value, parameters are as forconcord.utils.luminosity()
- Parameters:
F_pk – peak burst flux, MINBAR units assumed if not present
- Returns:
dictionary including calculation results and assumed distributions (as for luminosity)
- concord.utils.Q_nuc(Xbar, quadratic=True, old_relation=False, coeff=False)
This function implements the approximation to the nuclear energy generation rate Q_nuc given the mean hydrogen fraction in the fuel layer, Xbar, determined by Goodwin et al. 2019 (ApJ 870, #64). The value returned is in units of MeV/nucleon:
>>> q = Q_nuc(0.73) 5.35517578
There are three versions of the approximation, which can be selected by using the quadratic and old_relation flags; by default the most precise (and recent) version is used:
>>> q = Q_nuc(0.73, quadratic=False) 5.758715 >>> q = Q_nuc(0.73, old_relation=True) 4.52
If you want the coefficients, rather than the actual Q_nuc value, you can use the coeff flag (in which case the first argument is ignored):
>>> q_0, q_1 = Q_nuc(0.0,quadratic=False,coeff=True) [1.3455, 6.0455]
- Parameters:
Xbar – average hydrogen mass fraction of the burst column
quadratic – set to
True
to use the (more accurate) quadratic approximationold_relation – set to
True
to use the earlier approximationcoeff – set to
True
to return the coefficients used, in which caseXbar
is ignored
- Returns:
Q_nuc value, in MeV/nucleon, or the coefficients if
coeff=True
- concord.utils.X_0(Xbar, zcno, tdel, opz=1.259, debug=False, old_relation=False)
Routine (extracted from
concord.utils.hfrac()
) to determine the accreted H fraction X_0 given the average H-fraction at ignition, Xbar, the CNO metallicity zcno, and the burst recurrence time, tdel.The time to exhaust the accreted hydrogen is calculated according to the formula derived by Lampe et al. (2016, ApJ 819, #46); if the
old_relation=True
flag is set, the alternative prefactor quoted by Galloway et al. (2004 ApJ 601, 466) is used instead- Parameters:
Xbar – average hydrogen mass fraction of the burst column
zcno – CNO mass fraction in the burst fuel
tdel – burst recurrence time
opz – NS surface redshift, 1+z
old_relation – use the older expression for the time to exhaust the accreted hydrogen
debug – display debugging messages
- Returns:
inferred hydrogen mass fraction X_0 of the accreted fuel
- concord.utils.alpha(_tdel, _fluen, _fper, c_bol=None, nsamp=1000, conf=68.0, fulldist=False)
Routine to calculate alpha, the observed ratio of persistent to burst energy (integrated over the recurrence time, tdel) from the input measurables
Usage:
>>> alpha = cd.alpha(2.681, 0.381, 3.72, 1.45) >>> print (alpha) 136.64233701
>>> alpha = cd.alpha((2.681, 0.007), (0.381, 0.003), (3.72, 0.18), (1.45, 0.09)) >>> print (alpha) [136.23913536 9.88340505 11.21502717]
- Parameters:
tdel – burst recurrence time
fluen – burst fluence
fper – persistent flux
c_bol – bolometric correction on persistent flux
- Returns:
alpha-values, either a scalar (if all the inputs are also scalars); a central value and confidence range; or (if
fulldist=True
) a dictionary with the distributions of the alpha value and the input or adopted distributions of the intermediate values
- concord.utils.asym_norm(m, sigm=None, sigp=None, nsamp=1000, statistics='max', positive=False, model=1)
Draw samples from an asymmetric error distribution characterised by a mean and upper and lower 68% confidence intervals (sigp and sigm, respectively). Used primarily to generate distributions from quantities with asymmetric errors, by
concord.utils.value_to_dist()
.Follows the treatment of Barlow (2003)
With the positive flag set, it will continue drawing samples until all are > 0. Note that the resulting distribution may not quite have the required shape
- Parameters:
m – mean (central) value for the distribution
sigm – lower 68th-percentile uncertainty
sigp – upper 68th-percentile uncertainty
nsamp – number of samples required
statistics – convention for input stats for asymmetric distributions
positive – ensure all the samples are positive (may affect the distribution)
model – implementation of the distribution function; only one option is currently implemented
- Returns:
distribution with
nsamp
values having the desired shape
- concord.utils.calc_mr(g, opz)
‘ Calculates neutron star mass and radius given a surface gravity and redshift. This routine is required when (for example) comparing a model computed at a fixed gravity, with observations requiring a particular value of redshift (to stretch the model lightcurve to match the observed one). The combination of surface gravity and redshift implies unique mass and radius, and the constraints on 1+z can be translated to constraints on M, R
- Parameters:
g – surface gravity, in units equivalent to cm/s**2
opz – surface redshift 1+z
- Returns:
neutron star mass in g, radius in cm
- concord.utils.check_M_R_opz(M, R, opz)
Utility routine to check the consistency of the mass, radius and redshift passed to various routines (e.g. mdot)
- Parameters:
M – neutron star mass, or None
R – neutron star radius, or None
opz – 1+z where z is the surface gravitational redshift, or None
- Returns:
boolean,
True
if the values are consistent,False
otherwise
- concord.utils.create_logger()
Create a logger instance where messages are sent. See https://docs.python.org/3/library/logging.html
- Returns:
logger instance
- concord.utils.decode_LaTeX(string, delim='pm')
This function converts a LaTeX numerical value (with error) to one or more floating point values. It is expected that the number is formatted as
$3.1\pm1.2$
, and the default delimiter assumes this scenarioWill not work on asymmetric errors, expressed (for example) as
$1.4_{-0.1}^{+0.2}$
Example usage:
>>> cd.decode_LaTeX('$1.4\pm0.3') (1.4, 0.3)
- Parameters:
string – string with LaTeX numerical value and error
delim – separator of number and uncertainty; defaults to
\pm
- Returns:
tuple giving the value and error(s)
- concord.utils.dist(_F_pk, nsamp=None, dip=False, empirical=False, X=0.0, M=<Quantity 1.4 solMass>, opz=1.259, T_e=0.0, isotropic=False, inclination=None, imin=0.0, imax=75.0, model='he16_a', conf=68.0, fulldist=False, plot=False)
This routine estimates the distance to the source, based on the measured peak flux of a radius-expansion burst. Two main options are available;
empirical=False
(default), in which case the Eddington luminosity is calculated as equation 7 from Galloway et al. (2008, ApJS 179, 360); orempirical=True
, in which case the estimate of Kuulkers et al. 2003 (A&A 399, 663) is used. The defaultisotropic=False
option also takes into account the likely effect of anisotropic burst emission, based on the models provided by concordExample usage:
>>> cd.dist((30., 3.), isotropic=True, empirical=True) <Quantity [10.27267949, 0.50629732, 0.59722829] kpc> >>> cd.dist((30., 2., 5.), isotropic=False, empirical=True) <Quantity [10.75035916, 1.39721634, 1.461003 ] kpc> >>> cd.dist((30., 3.), isotropic=False, empirical=True, fulldist=True, nsamp=100) {'dist': <QuantityDistribution [11.12852694, 11.4435149 , 9.79809775, 9.82146105, 9.34573615, 10.4330567 , 10.99421133, 10.66816313, 8.83667458, 12.83249538, 14.15062276, 13.26883293, 9.80702496, 8.72915705, 12.04546987, . . .
- Parameters:
_F_pk – peak burst flux
nsamp – number of samples required for the distributions
dip – set to True if the source is a “dipper”
empirical – flag to use the empirical Eddington luminosity
X – hydrogen mass fraction in the atmosphere
M – neutron-star mass
opz – neutron-star surface redshift
T_e – temperature factor used in the theoretical L_Edd expression
isotropic – set to True to assume isotropic emission
inclination – inclination value (or distribution)
imin – minimum of allowed inclination range
imax – maximum of allowed inclination range
conf – confidence interval for output limits
fulldist – set to True to output full distributions for all parameters
plot – set to True to generate a simple plot
- Returns:
inferred distance values, either a scalar (if all the inputs are also scalars); a central value and confidence range; or (if
fulldist=True
) a dictionary with the distributions of the distance and the input or adopted distributions of the other parameters
- concord.utils.g(M, R, Newt=False, units='cm/s^2')
This function calculates the surface gravity given a mass M and radius R, using the Newtonian expression if the flag Newt is set to True
The result is returned in units of cm/s^2 by default
- Parameters:
M – neutron-star mass
R – neutron-star radius
Newt – default
False
returns the GR value, set to True to calculate the Newtonian valueunits – the required units
- Returns:
the surface gravity, in the required units
- concord.utils.hfrac(_tdel, _alpha=None, fper=None, fluen=None, c_bol=None, zcno=0.02, opz=1.259, old_relation=False, isotropic=False, inclination=None, imin=0.0, imax=75.0, model='he16_a', conf=68.0, fulldist=False, nsamp=None, debug=False)
Estimates the H-fraction at burst ignition, based on the burst properties In the absence of the alpha-value(s), you need to supply the persistent flux and burst fluence, along with the recurrence time, so that alpha can be calculated
There is also a mode for calculating a single value of the parameters, based on a single value of alpha, tdel and the inclination; this mode is used by the function itself, within a loop over the inclination values.
Example usage:
>>> import astropy.units as u >>> cd.hfrac(2.5*u.hr, 140., inclination=30.) ** WARNING ** assuming inclination in degrees (<Quantity 0.19564785>, <Quantity 0.26656582>, 30.0)
>>> cd.hfrac((2.681, 0.007), fluen=(0.381, 0.003), fper=(3.72, 0.18), c_bol=(1.45, 0.09),nsamp=100,fulldist=True) {'xbar': NdarrayDistribution([ 0.15246606, 0.20964612, 0.14169592, 0.11998812, 0.14293726, 0.13757633, -0.0850957 , 0.27077167, 0.01620781, 0.13673547, -0.00821256, 0.10193499, -0.08151998, -0.0370918 , 0.26400253, 0.1813988 , -0.04371571, 0.14432454, 0.28422351, -0.04962202, . . . 'X_0': NdarrayDistribution([ 0.22852126, 0.28567535, 0.21799341, 0.19609461, 0.21895133, . . . 'i': <QuantityDistribution [45.58518673, 15.61321579, 36.30632308, 48.18062101, 45.46515182, . . .
- Parameters:
tdel – burst recurrence time
alpha – burst alpha, which (if not supplied) can be calculated instead from the fluence, fper, and c_bol
fper – persistent flux
fluen – burst fluence
c_bol – bolometric correction on persistent flux
opz – neutron-star surface redshift
zcno – CNO mass fraction in the burst fuel
old_relation – flag to use the old (incorrect) relation for Q_nuc
isotropic – set to True to assume isotropic emission
inclination – inclination value (or distribution)
imin – minimum of allowed inclination range
imax – maximum of allowed inclination range
conf – confidence interval for output limits
model – model string of He & Keek (2016), defaults to “A”
conf – confidence % level for limits, ignored if fulldist=True
fulldist – set to True to return the distributions of each parameter
debug – display debugging messages
- Returns:
a tuple of values of Xbar, X_0 and inclination; or a dictionary including distributions of each of the values, as well as the distributions adopted for the observables
- concord.utils.homogenize_params(theta, nsamp=None)
This method is used by the MC routines, to ensure consistency of the parameter set. We want to make sure that each parameter (if a distribution) has the same length as any other parameter (unless it’s a scalar). Each parameter is labeled and provided with a value (and possibly also error), and units. Inclination is also provided with the isotropy flag, and the angle limits.
If none of the input parameters are already distributions, the provided number of samples nsamp will be used as the dimensions of the output distributions (this may be redundant, since we now also do the inclination distributions here)
Usage:
(par1, par2, … nsamp ) = homogenize_params( dictionary_of_input_params_and_units, nsamp )
Example usage:
>>> F_pk, _nsamp = homogenize_params( {'fpeak': ((33., 2.), cd.MINBAR_FLUX_UNIT)}, 100)
>>> F_pk, _inclination, _nsamp = cd.homogenize_params( {'fpeak': ((33., 2.), cd.MINBAR_FLUX_UNIT), 'incl': (None, u.deg, isotropic, 0.0, 75.)}, 100)
- Parameters:
theta – dictionary with parameters and units
nsamp – desired size of distribution objects, if not already set by one or more of the parameters
- Returns:
tuple with all the parameters in the same order they’re passed, plus the actual number of samples per parameter
- concord.utils.intvl_to_errors(perc)
This function converts a confidence interval defined by a 3-element array to the more useful version with errors .
- Parameters:
perc – the array of percentiles (central, lower limit, upper limit)
- Returns:
3-element array with (central value, lower uncertainty, upper uncertainty)
- concord.utils.iso_dist(nsamp=1000, imin=0.0, imax=75.0)
Routine to generate an isotropic distribution of inclinations (angle from the system rotation axis to the line of sight) from imin up to some maximum value, defaulting to 75 degrees. This value corresponds to the likely maximum possible for a non-dipping source.
- Parameters:
nsamp – number of samples required
imin – lower limit of distribution in degrees (default is zero)
imax – upper limit of distribution in degrees (default is IMAX_NDIP)
- Returns:
distribution of inclinations between
imin
,imax
- concord.utils.len_dist(d)
Utility routine to replace the len function for distributions, and also return sensible values for scalars rather than throwing an error
- Parameters:
d – object (scalar, array or distribution)
- Returns:
length of the object (0 if empty/
None
)
- concord.utils.lum_to_flux(_lum, dist=None, c_bol=None, nsamp=None, burst=True, dip=False, isotropic=False, inclination=None, imin=0.0, imax=75.0, model='he16_a', conf=68.0, fulldist=False, plot=False)
This routine converts luminosity to flux, based on the provided distance Basically the inverse of the
concord.utils.luminosity()
function- Parameters:
lum – luminosity to convert to flux
dist – distance to the source
c_bol – bolometric correction to apply
nsamp – number of samples
burst – set to True for burst emission, to select the model anisotropy for bursts
dip – set to True for a dipping source (deprecated)
isotropic – set to True if isotropic value required
inclination – system inclination or distribution thereof
imin – minimum inclination for generated distribution
imax – maximum inclination for generated distribution
model – model string of He & Keek (2016), defaults to “A”
conf – confidence % level for limits, ignored if fulldist=True
fulldist – set to True to return the distributions of each parameter
plot – plots the resulting distributions
- Returns:
inferred flux, either a scalar (if all the inputs are also scalars); a central value and confidence range; or (if
fulldist=True
) a dictionary with the distributions of the flux and the input or adopted distributions of the intermediate values
- concord.utils.luminosity(_F_X, dist=None, c_bol=None, nsamp=None, burst=True, dip=False, isotropic=False, inclination=None, imin=0.0, imax=75.0, model='he16_a', conf=68.0, fulldist=False, plot=False)
Calculate the inferred luminosity given a measured X-ray flux and distance The flux, error and distance can be single values, or arrays (in which case it’s assumed that they’re distributions)
This is a more sophisticated version of the L_Edd routine (forgot I had that)
Example usage:
Calculate the isotropic luminosity corresponding to a flux of 3e-9 erg/cm^2/s at 7.3 kpc
>>> import concord as cd cd.luminosity(3e-9,7.3,isotropic=True)
Calculate the range of luminosities corresponding to a persistent flux of 3e-9 erg/cm^2/s at 7.3 kpc, and assuming isotropic inclination distribution (i < 75 deg)
>>> cd.luminosity(3e-9,7.3,burst=False)
Calculate the range of luminosities corresponding to a burst flux of 3e-8, with uncertainty 1e-9, and an inclination of 45-60 degrees, plot (and return) the resulting distribution
>>> cd.luminosity((3e-8,1e-9),7.3,imin=45,imax=60,plot=True,fulldist=True)
- Parameters:
F_X – X-ray flux, MINBAR units assumed if not present
dist – distance to the source
c_bol – bolometric correction to apply
nsamp – number of samples
burst – set to True for burst emission, to select the model anisotropy for bursts
dip – set to True for a dipping source (deprecated)
isotropic – set to True if isotropic value required
inclination – system inclination or distribution thereof
imin – minimum inclination for generated distribution
imax – maximum inclination for generated distribution
model – model string of He & Keek (2016), defaults to “A”
conf – confidence % level for limits, ignored if fulldist=True
fulldist – set to True to return the distributions of each parameter
plot – plots the resulting distributions
- Returns:
luminosity, either a scalar (if all the inputs are also scalars); a central value and confidence range; or (if
fulldist=True
) a dictionary with the distributions of the luminosity and the input or adopted distributions of the intermediate values
- concord.utils.mdot(_F_per, _dist, c_bol=None, M=None, R=None, opz=None, isotropic=False, inclination=None, imin=0.0, imax=75.0, dip=False, model='he16_a', nsamp=None, conf=68.0, fulldist=False)
Routine to estimate the mdot given a (bolometric) persistent flux, distance, and inclination. This calculation was adapted initially from equation 2 of Galloway et al. (2008, ApJS 179, 360), and uses an approximation to \(Q_\rm{grav} = c^2 z/(1+z) \approx GM_{\rm NS}/R_{\rm NS}\) which is good to about 10% for a typical neutron star
Usage:
>>> import concord as cd # calculate the mdot corresponding to a flux of 1e-9 at 10kpc, for a 10km # 1.4 M_sun neutron star (this is the prefactor for equation 2 in Galloway # et al. 2008) >>> import astropy.constants as c >>> cd.mdot(1., 10., M=1.4*c.M_sun, R=10.*u.km, isotropic=True) WARNING:homogenize_params:no bolometric correction applied <Quantity 6691.30392224 g / (cm2 s)>
- Parameters:
F_per – persistent flux, MINBAR units assumed if not present
dist – distance to the source
c_bol – bolometric correction to apply
M – neutron-star mass
R – neutron-star radius
opz – 1+z where z is the surface gravitational redshift, or None
isotropic – set to True if isotropic value required
inclination – system inclination or distribution thereof
imin – minimum inclination for generated distribution
imax – maximum inclination for generated distribution
dip – set to True for a dipping source (deprecated)
model – model string of He & Keek (2016), defaults to “A”
nsamp – number of samples to generate
conf – confidence % level for limits, ignored if fulldist=True
fulldist – set to True to return the distributions of each parameter
- Returns:
mdot-values, either a scalar (if all the inputs are also scalars); a central value and confidence range; or (if
fulldist=True
) a dictionary with the distributions of the mdot value and the input or adopted distributions of the intermediate values
- concord.utils.redshift(M, R)
This function calculates the gravitational redshift 1+z
- Parameters:
M – neutron-star mass
R – neutron-star radius
- Returns:
1+z
- concord.utils.solve_radius(M, R_Newt, eta=1e-06)
This routine determines the GR radius given the NS mass and Newtonian radius, assuming the GR and Newtonian masses are identical
Solving is tricky so we just use an iterative approach
- Parameters:
M – neutron-star mass
R_Newt – Newtonian neutron-star radius
eta – tolerance for convergence
- Returns:
neutron-star radius R in GR
- concord.utils.tdel_dist(nburst, exp, nsamp=1000)
Function to generate a synthetic distribution of recurrence times, given nburst observed events over a total exposure of exp. The resulting distribution may approximate the PDF of the underlying burst rate, but does assume that the bursts are independent, which is generally not the case
- Parameters:
nburst – number of events detected
exp – total exposure time
nsamp – number of samples to generate
- Returns:
distribution giving the PDF of the average recurrence time
- concord.utils.value_to_dist(_num, nsamp=1000, unit=None, statistics='max', verbose=False)
This method converts a measurement to a distribution, to allow flexibility in how values are implemented in the various routines. Primarily used by
concord.utils.homogenize_distributions()
.Intended to be very flexible in terms of the way the units are provided, but unfortunately np.shape for a tuple with elements having units, does not work! So have to pass the units on the entire object
Also have an issue with combining distributions with and without units; “…the NdarrayDistribution will not combine with Quantity objects containing units” see https://docs.astropy.org/en/stable/uncertainty for more details
- Parameters:
num – scalar/array to convert to distribution
statistics – convention for input stats for asymmetric distributions
- Returns:
astropy distribution object
Example usage:
>>> y = cd.value_to_dist(3.) # scalars are kept as single values >>> z = cd.value_to_dist((3.,0.1),nsamp=10) # generate 10 samples from a normal distribution around 3.0 with st. dev 0.1 >>> a = cd.value_to_dist((3.,0.5,0.1)) # generate default number of samples from an asymmetric Gaussian >>> b = cd.value_to_dist(a.distribution) # convert an array to a distribution
with units:
>>> z = cd.value_to_dist((3.,0.1)*u.hr,nsamp=10) # generate 10 samples from a normal distribution around 3.0 with st. dev 0.1
- concord.utils.yign(_E_b, dist=None, nsamp=None, R=<Quantity 11.2 km>, opz=1.259, Xbar=0.7, quadratic=False, old_relation=False, isotropic=False, inclination=None, imin=0.0, imax=75.0, dip=False, model='he16_a', conf=68.0, fulldist=False)
Calculate the burst column from the burst fluence, adapted initially from equation 4 of Galloway et al. (2008, ApJS 179, 360)
- Parameters:
_E_b – burst fluence
_dist – distance ot bursting source
R – neutron star radius
opz – 1+z, surface redshift
Xbar – mean H-fraction at ignition
quadratic – flag to use the quadratic expression for Q_nuc
old_relation – flag to use the old (incorrect) relation for Q_nuc
isotropic – assume isotropic flux (or not)
inclination – inclination value (or array; in degrees)
imin – minimum inclination (degrees)
imax – maximum inclination (degrees)
dip – whether or not the source is a dipper
model – model for anisotropy calculation
nsamp – number of samples to generate
conf – confidence interval for uncertainties; or
fulldist – output the full distribution of calculated values
- Returns:
inferred ignition column values, either a scalar (if all the inputs are also scalars); a central value and confidence range; or (if
fulldist=True
) a dictionary with the distributions of the ignition column and the input or adopted distributions of the other parameters
Classes
- class concord.burstclass.KeplerBurst(filename=None, run_id=None, path=None, source='gs1826', grid_table=None, batch=None, run=None, **kwargs)
Example simulated burst class. Apart from the lightcurve (which is defined with time, lumin, and lumin_err columns), the additional (minimal) attributes requred are:
- Filename:
source file name
- Tdel, tdel_err:
recurrence time and error (hr)
- Lacc:
accretion luminosity (in units of Mdot_Edd)
- R_NS:
model-assumed radius of the neutron star (GR)
- G:
surface gravity assumed for the run
An example call is as follows:
>>> c_loZ=KeplerBurst(filename='mean1.data',path='kepler', lAcc=0.1164,Z=0.005,H=0.7, tdel=4.06/opz,tdel_err=0.17/opz, g = 1.858e+14*u.cm/u.s**2, R_NS=11.2*u.km)
If you provide a batch and run value, like so:
>>> c = KeplerBurst(batch=2,run=9)
or (for the older results) a run_id value, like so:
>>> c = KeplerBurst(run_id='a028',path='../../burst/reference')
the appropriate model lightcurve will be read in, and the model parameters will be populated from the table
- info()
Display information about the KeplerBurst
- class concord.burstclass.Lightcurve(*args, **kwargs)
The fundamental burst object, which has attributes like time (array), flux/luminosity etc. and methods including plot. Minimal attributes are:
- Time, dt:
array of times and time bin duration
- Timepixr:
set to 0.0 (default) if time is the start of the bin, 0.5 for midpoint
- Flux, flux_err:
flux and error (no units assumed)
- Lumin, lumin_err:
luminosity and error
Normally only one of flux or luminosity would be supplied
Example: simulated burst (giving time and luminosity)
>>> Lightcurve(self, time=d['time']*u.s, lumin=d['luminosity']*u.erg/u.s, lumin_err=d['u_luminosity']*u.erg/u.s)
Example: observed burst giving flux
>>> Lightcurve(self, time=d['col1']*u.s, dt=d['col2']*u.s, flux=d['col3']*1e-9*u.erg/u.cm**2/u.s, flux_err=d['col4']*1e-9*u.erg/u.cm**2/u.s)
- fluence(plot=False, warnings=True)
Calculate the fluence for a lightcurve, following the approach from
get_burst_data
, as used for MINBARThis code relies on the presence of the
dt_nogap
andgood
attributes, which should be calculated when the Lightcurve object is defined- Parameters:
plot – set to
True
to display a diagnostic plotwarnings – set to
False
to suppress display of warnings
- Returns:
the burst fluence and uncertainty
- observe(param=[<Quantity 6.1 kpc>, <Quantity 60. deg>, 1.26, <Quantity -10. s>], obs=None, disc_model='he16_a', c_bol=1.0)
Convert a luminosity profile to a simulated observation, with plausible errors. See Keek & Heger (2011, ApJ 743, #189) for the background calculations
This method might make more sense as part of the :py:class:`KeplerBurst`_ class (or a parent
SimulatedBurst
class), but for now we include this method as part of the lightcurve class, so that future classes of model bursts can access it- Parameters:
param – tuple with the distance, inclination, redshift and time offset
obs – template observed burst to use for the time bins. If not supplied, uniform 0.25 s bins will be used, covering the burst
disc_model – choice of disc model for the anisotropy
c_bol – bolometric correction to adjust the burst fluxes
- Returns:
- plot(yerror=True, obs_color='b', model_color='g', **kwargs)
Plot the lightcurve, accommodating both flux and luminosities
- Parameters:
yerror – display flux/luminosity errors (default
True
)obs_color – set the colour for an observed burst
model_color – set the colour for a simulated burst
- Returns:
the plot
- print()
Print the basic parameters of the lightcurve. Can be called in turn by the info commands for
concord.burstclass.KeplerBurst
,concord.burstclass.ObservedBurst
to display the properties of the parent class… once those are written- Returns:
- write(filename='lightcurve.csv', addhdr=None)
Write the lightcurve to a file
- Parameters:
filename – filename to save the data to
addhdr – additional text to include as a header (default
None
)
- class concord.burstclass.ObservedBurst(time, dt, flux, flux_err, **kwargs)
Observed burst class. Apart from the lightcurve (which is defined with
time
,flux
, andflux_err
columns), the additional attributes set are:- Filename:
source file name
- Tdel, tdel_err:
recurrence time and error (hr)
- Comments:
ASCII file header text
- Table_file:
LaTeX file of table 2 from Galloway et al. (2017)
- Table:
contents of table 2
- Row:
entry in the table corresponding to this burst
- Fper, fper_err:
persistent flux level (1e-9 erg/cm^2/s)
- C_bol:
bolometric correction
- Mdot:
accretion rate (total, not per unit area)
- Fluen_table:
burst fluence
- F_pk:
burst peak flux
- Alpha:
alpha value
The units for selected parameters are all stored as astropy.units objects
- compare(mburst, param=[<Quantity 6.1 kpc>, <Quantity 60. deg>, 1.26, <Quantity -10. s>], breakdown=False, plot=False, subplot=True, weights={'fluxwt': 1.0, 'tdelwt': 2500.0}, disc_model='he16_a', debug=False)
This is the key method for burst observation-model lightcurve comparisons, and generates a likelihood that can be used for MCMC (for example). It can be used to plot the observations with the models rescaled by the appropriate parameters.
weights
give the relative weight to the recurrence time and persistent flux for the likelihood. Since you have many more points in the lightcurve, you may want to weight these greater than one so that the MCMC code will try to match those preferentially- Parameters:
mburst – model burst for comparison, for example an
concord.burstclass.KeplerBurst
objectparam – tuple with the distance, inclination, redshift and time offset
breakdown – display the breakdown of the likelihood between the recurrence time, persistent flux and lightcurve
plot – set to
True
to show a plot of the comparisonsubplot – include a subplot with the recurrence time comparison (default
True
)weights – dict giving the relative weight to the recurrence time (key
tdelwt
) and persistent flux (fluxwt
) for the likelihood.disc_model – choice of disc model for the anisotropy
debug – display debugging information
- classmethod from_file(filename, path='.', **kwargs)
Method to read in a file and populate the relevant arrays to create an
concord.burstclass.ObservedBurst
object This routine will read in any ascii lightcurve file which matches the format of the “reference” bursts, i.e. with columns:'Time [s]' 'dt [s]' 'flux [10^-9 erg/cm^2/s]' 'flux error [10^-9 erg/cm^2/s]' 'blackbody temperature kT [keV]' 'kT error [keV]' 'blackbody normalisation K_bb [(km/d_10kpc)^2]' 'K_bb error [(km/d_10kpc)^2]' chi-sq -1.750 0.500 1.63 0.054 1.865 0.108 13.891 4.793 0.917 -1.250 0.500 2.88 1.005 1.862 0.246 22.220 4.443 1.034 -0.750 0.500 4.38 1.107 1.902 0.093 30.247 2.943 1.089 -0.250 0.500 6.57 0.463 1.909 0.080 46.936 6.969 0.849 . . .
Calling approach: >>> obs = ObservedBurst.from_file(filename) >>> obs = ObservedBurst.from_file(‘gs1826-24_3.530h.dat’,’example_data’)
- info()
Display information about the ObservedBurst
- Returns:
- classmethod minbar(id, remote=False)
Method to generate an observed burst from MINBAR, and populate the relevant arrays. Requires the minbar repo to be present (obviously)
Calling approach:
>>> obs = ObservedBurst.minbar(2203)
- Parameters:
id – the burst ID in MINBAR DR1
remote – if true, force download from the remote repo
- Returns:
- classmethod ref(source, dt)
Method to read in a reference burst and populate the relevant arrays to create an
concord.burstclass.ObservedBurst
Calling approach:
>>> obs = ObservedBurst.ref('GS 1826-24', 3.5)
- concord.burstclass.fper(mburst, param, c_bol=1.0)
Calculates the persistent flux, based on the supplied mdot, redshift etc. Earlier we passed the individual parameters, but easiest to just provide the model burst
- concord.burstclass.lhoodClass(params, obs, model, **kwargs)
Calculate the likelihood related to one or more model-observation comparisons The corresponding call to emcee will (necessarily) look something like this:
>>> sampler = emcee.EnsembleSampler(nwalkers, ndim, lhoodClass, args=[obs, models, weights])
Use the kwargs construction to pass the additional parameters weights, disc_model to the compare method if required
- Parameters:
params – tuple with the distance, inclination, redshift and as many time offsets as there are bursts to compare
obs – one or more :py:class:`concord.burstclass.ObservedBurst`s
model – one or more model bursts, one for each of the observed bursts; for example a
concord.burstclass.KeplerBurst
- Returns:
likelihood value
- concord.burstclass.modelFunc(p, obs, model, disc_model)
This function performs the stretching and rescaling of the (model predicted) burst lightcurve based on the input parameter set The calculations are based on those in Appendix B of Keek & Heger (2011, ApJ 743, #189)
- Parameters:
p – parameter array, (for now) a tuple with the appropriate units; ( distance, inclination, redshift, time offset )
obs – observed burst to use for the time bins
model – model burst to be stretched and rescaled
disc_model – choice of disc model for the anisotropy
- concord.burstclass.plot_comparison(obs, models, param=None, sampler=None, ibest=None)
This routine will plot each observation vs. each model, for a given set of parameters. If you supply an emcee sampler object, it will find the highest-likelihood sample and use that
- concord.burstclass.plot_contours(sampler, parameters=['$d$', '$i$', '$1+z$'], ignore=10, plot_size=6)
Simple routine to plot contours of the walkers, ignoring some initial fraction of the steps (the “burn-in” phase)
Documentation is here https://samreay.github.io/ChainConsumer/index.html