The information presented here is placed in the public domain, and was written by Doug Burke. The notebook used to create this page is available, and questions can be asked using the GitHub issues page or via Twitter: https://twitter.com/doug_burke.
It will probably make a lot-more sense if you have already read the previous entries in this "series".
In today's installment I want to look at the K correction term I mentioned in my first post when calculating the luminosity of a source given a flux and a redshift. However, there's a few things I have to get through first - such as how to represent a spectral model, integration of these models, converting between units - before I can get to the actual K correction. As that leaves me with about one third of a notebook to fill up, I repeat the process, but this time tracking quantities using the dimensional-tf package.
In [1]:
import Data.Time
getCurrentTime
Before I get to calculating the K correction for a model, I play around with creating a model that represents the spectrum of a source, as it turns out I'm going to want this later. Since I work for the Chandra X-ray Center, I am going to working with models of X-ray emitting sources, concentrating on the 0.1 to 10 keV pass band (I'll explain in a second why I am talking about wavelengths in terms of energy, but the short answer is because Astronomers).
In [2]:
-- Load up modules and set up the code to automatically display plots
-- from Charts.
--
import Numeric.Integration.TanhSinh
import IHaskell.Display
import Graphics.Rendering.Chart.Backend.Diagrams
import Graphics.Rendering.Chart.Easy
import qualified Numeric.Units.Dimensional.TF.Prelude as P
import Numeric.Units.Dimensional.TF.NonSI (electronVolt)
import qualified Numeric.NumType.TF as N
instance IHaskellDisplay (Renderable a) where
-- Pretend you didn't notice that in previous notebooks this
-- was written slightly differently - e.g.
--
-- do
-- (svg, _) <- renderableToSVG renderable 450 300
-- display svg
--
-- The following code is the same as above, but just in a more
-- compact (or, equivalently, obtuse) form.
--
display renderable = renderableToSVG renderable 450 300 >>= display . fst
Theoretically, you might expect for spectral models of X-ray sources - that is, the amount of "flux" from a source - to be
in SI units
a function of wavelength
However, via a combination of history, the fact that in X-ray astronomy we generally deal with the particle nature of light rather than its wave-like nature$^\dagger$, and the need to compare the model to data, X-ray Astrophysical models generally calculate the flux of a source in photon/cm$^2$/s as a function of energy. That is, if $S(E)$ is the source model, then the photon flux, $f$, is defined as
$$ f = \int_{\rm{elo}}^{\rm{ehi}} S(E) \rm{d}E $$where elo
and ehi
are the energy band we are interested in
(that is, the wavelength range we are using and $E = hc/\lambda$
to convert between wavelength and energy).
For sources observable by the
Chandra X-ray Observatory,
this energy range is close to 0.1 to 10 kilo electronVolts (referred
to as keV from now on), which corresponds to 12 to 0.12 nanometers.
This means that the units of $S(E)$ are going to be
photons/cm$^2$/s/keV, so it is a "flux density". Although not needed
here, I will later on need a conversion factor from keV to Joules, so let's
write it now.
$^\dagger$ The detectors on Chandra count photons, although some data is taken with a transmission grating in the light path, for which wavelength is the "natural unit", rather than energy.
In [3]:
-- Convert an energy in keV into ergs using dimensional-tf, but
-- note that the input and output values are plain Doubles.
--
enConv :: Double -> Double
enConv x = x P.*~ P.kilo electronVolt P./~ erg
-- Set up a unit to represent an erg.
erg :: P.Unit P.DEnergy Double
erg = P.prefix 1.0e-7 P.joule
As the comments indicate, this routine uses the
dimensional-tf
package,
but its use is "hidden", in that both the input and output values of enConv
are Double
values. As a reminder, the *~
operator from the package converts
a value into a type which "knows" about its dimensions, and /~
is its
inverse, "extracting" the value whilst checking that the dimensions match, and
applying any unit conversion (and, due to the precedence of the operators,
x *~ a /~ b
is equivalent to (x *~ a) /~ b
).
A quick check that this is working: I know that 1 keV should be about $1.602 \times 10^{-16}\ \rm{J}$, or $1.602 \times 10^{-9}\ \rm{erg}$:
In [4]:
enConv 1.0
Whilst I'm playing around with the dimensional-tf
package (more can be found at the end of the notebook),
I can also check that I got the conversion between energy and wavelength correct,
taking the product of the Plank constant and the speed of light, c
, from the
Wikipedia page:
In [5]:
hc = 1.98644568e-25 P.*~ (P.joule P.* P.meter)
l1 = hc P./ (1.0 P.*~ P.kilo electronVolt)
putStrLn ("The wavelength corresponding to 1 keV is " ++ show l1)
putStrLn (" which is " ++ show (l1 P./~ P.nano P.meter) ++ " nm")
So, if 1 keV is 1.2 nm then 0.1 keV is 12 nm and 10 keV is 0.12 nm, which is what I claimed. Success!
In the example I used
putStrLn
to display the string, using
show
to convert a value into a string and
++
to append strings.
If I had just left it to the notebook to display a string
then it would have also displayed the quotes - e.g.
"The wavelength corresponding to ..."
In general I let the notebook display values, but occasionally I'll use putStrLn
.
After that little diversion, I can get back to thinking about how I want to represent the spectral models by Haskell code. The obvious interface - at least to me - is a function that accepts a single parameter (an energy, in keV) and returns the flux density for that energy (in units of photon/cm$^2$/s/keV). That is, I can define the following type synonym to represent such a function:
In [6]:
type Model = Double -> Double
Now, on its own this turns out to be rather useless, since there is no way to select the model or change any parameters of a model! To get around this, I'll expect models to be partially applied, so that a model function will accept the model parameters and then end with the energy at which to evaluate.
Hopefully an example will help. A common model - both because of its simplicity, and because there are many physical processes which behave like this - is a power law, written as (at least, by X-ray astronomers)
$$S(E) = S_0 E^{-\gamma}$$where $\gamma$ is the slope of the power law$^\dagger$ and a normalisation factor $S_0$, which represents the flux density of the source at 1 keV (when $E$ is in units of $keV$). The power law gets steeper (so more signal at lower energies) as $\gamma$ increases.
$^\dagger$ it is often referred to as the photon index by X-ray astronomers
In [7]:
-- Create a power-law model of two parameters - the
-- normalisation at 1 keV and the slope (gamma) - which
-- is evaluated at an energy e (in keV).
--
plMdl :: Double -> Double -> Model
plMdl s0 gamma e = s0 * e**(-gamma)
Using this I can create a function of type Model
by specifying the two parameters:
for example
In [8]:
myModel = plMdl 1 2
:type myModel
which can then be evaluated at one or more energies:
In [9]:
myModel 5
map myModel [0.1, 1, 10]
Now, I don't have to use this approach of creating a separate function for
each model, since I can also use plMdl
directly; that is
In [10]:
plMdl 1 2 5
map (plMdl 1 2) [0.1, 1, 10]
and get the same result, but there are times when it's useful.
If I had not given a type signature, then the compiler would have inferred the type
plMdl :: Double -> Double -> Double -> Double
(well, actually it would have come up with something a bit-more general, involving the
Floating
type class,
but let's ignore that for now). Since Model
just means Double -> Double
, this is
the same as the signature I gave - namely
Double -> Double -> Model
- although I do admit that it can be surprising when
looking at a type signature that looks like it takes two arguments when it actually
takes three.
Since Model
is a type synonym - so that its just a textual substitution,
just like a #define
statement with the C pre-processor - I haven't gained
any type safety by making this substitution. Hopefully it has provided
some clarity by better documenting the intent of my code.
Now I want to visualize the model, so I create a function that plots a model over an energy range (with a hard-coded step size of 0.01 keV) and sets up the axis labels (I don't plan to explain how this code works, since it is quite complex; the Chart wiki contains examples and documentation, but I think it's perfectly fine at the moment to just think of this as "and now, the magic of plotting happens"):
In [11]:
-- plotModel :: Double -> Double -> Model -> Renderable ()
plotModel elo ehi mdl = toRenderable (do
layout_title .= "Model"
layout_x_axis . laxis_title .= "E (keV)"
layout_y_axis . laxis_title .= "photon / cm^2 / s / keV"
plot (line "" [cs])
)
where
de = 0.01
es = [elo, elo + de .. ehi]
cs = map (\e -> (e, mdl e)) es
Without an explicit type signature - such as the one given above as a comment - the inferred type is rather long:
In [12]:
:t plotModel
This says that the function takes two values of the same type (x
) and a function which converts from that type
to another (y
) and creates a Renderable
(I'm going to skip over the ()
part). The constraints on the x
type are
that it is
enumerable
(because the values are used in the statement [elo, elo+de .. ehi]
),
fractional
(since de
is not an integer),
and belongs to the mysterious
PlotValue
type class, which y
also has to. This typeclass is used
by Chart to indicate that a type can be included in a plot and, fortunately
for me, Double
values are an instance of it. I should also add that the choice of x
and y
as names of types here has nothing to do with X and Y axes of the graph, in case you were wondering.
In [13]:
plotModel 0.1 10 myModel
Hmmm, that's not particularly spectacular now, is it?
What I really want is to display the data on a log-log
axis (that is, use log scales for both axes). One way to do this is to "wrap" the values to be plotted
in the
LogValue
constructor,
which tells Chart to use a log scale. That is, in plotModel
, instead of
cs = map (\e -> (e, mdl e)) es
use
cs = map (\e -> (LogValue e, LogValue (mdl e))) es
To make this configurable, I create a helper function and then two
"plot" routines: plotModel
for linear axes and plotModelLL
for
log axes. I could also create linear-log and log-linear variants,
but I'll leave that as an exercise for the interested reader!
In [14]:
-- The xconv and yconv values allow the user to specify how
-- to convert the values: the obvious choices for these parameters,
-- at least if you have used Chart for a while, are id (linear)
-- or LogValue (log). I have given a signature just to be
-- explicit; I could have let the compiler infer this.
--
plotModelHelper :: (PlotValue x, PlotValue y) => (Double -> x) -> (Double -> y) -> Double -> Double -> Model -> Renderable ()
plotModelHelper xconv yconv elo ehi mdl = toRenderable (do
layout_title .= "Model"
layout_x_axis . laxis_title .= "E (keV)"
layout_y_axis . laxis_title .= "photon / cm^2 / s / keV"
plot (line "" [cs])
)
where
de = 0.01
es = [elo, elo + de .. ehi]
cs = map (\e -> (xconv e, yconv (mdl e))) es
-- plotModel and plotModelLL have the following type:
-- Double -> Double -> Model -> Renderable ()
--
plotModel = plotModelHelper id id
plotModelLL = plotModelHelper LogValue LogValue
The plotModel
variant uses the id
function, which has an interesting type:
:type id
a -> a
Due to parametricity, this function can only return the value it is sent (the compiler does not know anything about the type, so it can't create or change a value, so the only thing it can do is return the input$^\dagger$), which makes it seem rather pointless. However, it is surprisingly useful; in this particular case it is used where a conversion function is used but no conversion is actually needed.
$^\dagger$ well, it could also return an
error
or
undefined
value, but as both of these
would end up in a run-time error I'll just look slightly embarassed, mumble a bit, and forget I even
mentioned it.
The two versions of the plots for this model look like
In [15]:
plotModel 0.1 10 myModel
plotModelLL 0.1 10 myModel
How about a model that looks a bit-more interesting, such as a gaussian, defined as
$$S(E) = S_0 e^{-0.5 (E-E_0)^2 / \sigma^2}$$
In [16]:
-- The values E, E0, and sigma are in keV
--
gMdl :: Double -> Double -> Double -> Model
gMdl s0 e0 sigma e = let x = (e-e0)/sigma in s0 * exp (-0.5 * x * x)
In [17]:
plotModel 0.1 10 (gMdl 10 6.7 1.8)
The plot appears cut off at the high-energy end because I chose an upper limit of 10 keV. Increasing this upper limit leads to more of the Gaussian being drawn:
In [18]:
plotModel 0.1 15 (gMdl 10 6.7 1.8)
Compound models can also be created; cMdl
also has a type of
Model
and can be used in the same way as the individual components:
In [19]:
cMdl e = gMdl 1.89 4.3 1.2 e + plMdl 0.023 2 e
plotModel 0.1 15 cMdl
plotModelLL 0.1 15 cMdl
In order to calculate the K-correction for a model I need the flux
for a model over a set energy range; that is, I want
$\int S(E) dE$.
The models I have used so far have
analytic solutions for this integral,
but instead of actually doing any thinking$^\dagger$,
I'm just going to use the same
integration technique as in previous notebooks and use
an approximation (this time with simpson
rather than
trap
in case the models are not well behaved).
$^\dagger$ also known knowadays as 'using Wolfram Alpha'
In [20]:
-- Integrate up a model, in units of photons/cm^2/s/keV, over
-- the given energy range, in keV, with a fixed absolute
-- tolerance. The flux is in units of photons/cm^2/s.
--
iModel :: Double -> Double -> Model -> Double
iModel elo ehi mdl = result (absolute 1.0e-6 (simpson mdl elo ehi))
Using the power law model, this calculates the 0.1 to 10 keV photon flux to be:
In [21]:
iModel 0.1 10 myModel
As a quick check, in part to verify that the integration package is working but mainly to show off some simple LaTeX, the analytic solution for the power law model is
$$ \begin{align} \int_{e_1}^{e_2} S_0 E^{-\gamma} {\rm d}E && = && \frac{S_0}{1-\gamma} (e_2^{1-\gamma} - e_1^{1-\gamma}) && \gamma \ne 1 \\ && = && S_0 \rm{log}_e(e_2 / e_1) && \gamma = 1 \end{align} $$This can be written as (neglecting the fact that I should not be using an equality check for floating-point numbers, so this code is not robust):
In [22]:
exact s0 gamma e1 e2 | gamma == 1 = s0 * log (e2 / e1)
| otherwise = let p = 1-gamma in s0 * (e2**p - e1**p) / p
exact 1 2 0.1 10
So the approximate value calculated by iModel
is essentially the same as the exact solution
(for the accuracy I need). The exact solution could be calculated for the gaussian model, but I'll leave that
for the interested reader (the
erf
or
hmatrix-special
packages may be of interest, although I'm sure there's other packages that provide
an erf
function).
Actually, what I really need is not the flux in photon/cm$^2$/s but in erg/cm$^2$/s, that is
$\int E S(E) \rm {d}E$. First up, I want a simple way to convert an
energy from keV to erg, which is just what enConv
, created way back at the
start of the notebook, does. Using this, I can modify iModel
so that
it integrates up energy * photon flux density (i.e. energy flux density)
instead:
In [23]:
iFModel :: Double -> Double -> Model -> Double
iFModel elo ehi mdl =
-- I create a simple "wrapper" function to integrate that
-- evaluates the supplied model and then multiplies it by
-- the energy.
let wmdl x = enConv x * mdl x
in result (absolute 1.0e-6 (simpson wmdl elo ehi))
Note that there's nothing in the types that tell me what the units are; that is,
the compiler can't tell me if I've made a mistake and used iFModel
when
I meant iModel
, or used energies in eV
rather than keV
. It's left to
the user - along with the documentation - to ensure things are correct.
Once I've got to calculating the K correction, I show how the
dimensional-tf
package can be used to rewrite this code and ensure
that the dimensions match. There are also alternative approaches,
such as creating types to represent the values of interest, which are
less general but may be easier to use. I may play around with this idea
in a future notebook.
With all this, I can compare the two fluxes (photon and energy):
In [24]:
putStrLn ("flux = " ++ show (iModel 0.1 10 myModel) ++ " photon/cm^2/s")
putStrLn (" = " ++ show (iFModel 0.1 10 myModel) ++ " erg/cm^2/s")
I can compare these results to those from the software that X-ray astronomers use. In this case I use Sherpa, which is a Python-based analysis system. The following commands set up a power-law model with $\gamma=2$, normalisation of 1 and then calculate the photon and energy fluxes over the range 0.1 to 10 keV:
sherpa> dataspace1d(start=0.01, stop=11, step=0.01, id=1)
sherpa> set_source(powlaw1d.plmdl)
sherpa> plmdl.gamma = 2
sherpa> plmdl.ampl = 1
sherpa> calc_photon_flux(lo=0.1, hi=10, id=1)
9.9000999000999279
sherpa> calc_energy_flux(lo=0.1, hi=10, id=1)
7.3812306496266669e-09
sherpa> calc_photon_flux(lo=0.1, hi=10, id=1) / 9.899999999978695
1.0000100909213367
sherpa> calc_energy_flux(lo=0.1, hi=10, id=1) / 7.3833833576195705e-9
0.99970843881610427
The numbers are very close$^\dagger$, but it isn't quite a fair comparison, since the Sherpa models are designed to be compared to binned data, which means that the energy flux is approximated by calculating the photon flux within a bin (i.e. range of energies) and then multiplied by the mid-point energy of the bin. That is, the energy flux, $f_{\rm E}$, is calculated using
$$ f_{\rm E} \simeq \sum_{i} \frac{e_{\rm{lo};i} + e_{\rm{hi};i}}{2} s(e_{\rm{lo};i}, e_{\rm{hi};i}) $$where $s(a,b)$ represents the model integrated over the energy range $a$ to $b$,
and the bins cover the desired energy range (in the example above
they are defined by the dataspace1d
call). In early drafts of this
notebook I went off and showed that I get similar results using
this approximation scheme, but I think I'll leave that for another time,
and get back to the K correction.
$^\dagger$ This level of agreement is significantly better than the accuracy we can get for most Astronomical measurements!
The K correction is used in the equation to convert a flux ($f_X$) into a luminosity ($L_X$) for a source at a redshift $z$:
$$ L_X = 4 \pi d_L(z)^2 k(z) f_X $$where $d_L$ is the luminosity distance and $k(z)$ is the k correction. The reason it is needed is that if you have a source far enough away that you need to worry about including General Relativistic effects when calculating $d_L$, then you also need to account for the fact that - due to the expansion of space - the light we observe has been "redshifted". That is, if a source has a feature at a energy $E_0$ then we would observe that feature at an energy of $E_0/(1+z)$ when the source is at a redshift of $z$. Hopefully an image, created with 1000 characters or so, will be better than 1000 words:
In [25]:
-- Plot up a model spectrum, showing the "observed" and "rest frame"
-- bands if the source is at z=2 for measurements made over
-- the 1 to 5 keV pass band.
--
import Data.Colour.Names (orange, red)
toRenderable (do
let (elo, ehi, z) = (1, 5, 2)
es = [0.5, 0.51 .. 30]
f e = (LogValue e, LogValue (cMdl e))
g e = (LogValue e, (LogValue 1e-5, LogValue (cMdl e)))
cs = map f es
zp1 = 1 + z
er = [elo, ehi]
getE e1 e2 = takeWhile (<= e2) . dropWhile (< e1)
rest = map g (getE (elo*zp1) (ehi*zp1) es)
obs = map g (getE elo ehi es)
t = "Observation at z=" ++ show z ++ " for " ++
show elo ++ " to " ++ show ehi ++ " keV"
fillBetween color title vs = liftEC $ do
plot_fillbetween_title .= title
plot_fillbetween_style .= solidFillStyle (color `withOpacity` 0.5)
plot_fillbetween_values .= vs
layout_title .= t
layout_x_axis . laxis_title .= "E (keV)"
layout_y_axis . laxis_title .= "photon / cm^2 / s / keV"
plot (fillBetween orange "observed" rest)
plot (fillBetween red "wanted" obs)
plot (line "model" [cs])
)
So, although we are observing in the 1 to 5 keV range (in this example), the emission we measure is actually from the 3 to 15 keV range (i.e. $1+z=3$), which will not (for virtually all objects we are interested in) have the same flux. In this case this is indicated by the two different pass bands (i.e. 1-5 and 3-15 keV) having different shapes (and, if you summed up the model) fluxes. If this shift is not corrected for then you can not compare luminosities for objects at different redshifts. What we need is a correction for this, which is where the K correction comes in.
Following Appendix B of Jones et al. 1998, ApJ, 495, 100-114 (or the calculating K corrections document I wrote for people analysing Chandra data), the K correction is calculated as the ratio of the flux of the source in the rest and observed frames, namely
$$ k(z) = \frac{\int_{e_1}^{e_2} E S(E) {\rm d}E}{\int_{e_1 (1+z)}^{e_2 (1+z)} E S(E) {\rm d}E} $$One way of writing this is to send in the observed frame energy
band (elo
to ehi
), a redshift (z
), and a model:
In [26]:
-- The arguments are energy band (low, high) in keV,
-- the redshift, and the model.
--
kcorr :: Double -> Double -> Double -> Model -> Double
kcorr elo ehi z mdl =
let obs = iFModel (elo * zp1) (ehi * zp1) mdl
rest = iFModel elo ehi mdl
zp1 = 1 + z
in rest / obs
As a quick test:
In [27]:
kcorr 0.1 10 0.5 myModel
kcorr 0.1 10 0.5 (plMdl 1 3)
Interestingly enough, the "default" model I have been using - a power-law with $\gamma=2$ - actually has a K correction of unity here. Is this a fluke, or did I pick a special value? Well, for this model the K correction can be calculated analytically by expanding out $S(E) = S_0 E^{-\gamma}$:
$$ \begin{align} k(z) && = && \frac{\int_{e_1}^{e_2} E^{1-\gamma} {\rm d}E}{\int_{e_1 (1+z)}^{e_2 (1+z)} E^{1-\gamma} {\rm d}E} \\ && = && (1+z)^{\gamma-2} \end{align} $$So, for a power law, the k correction:
myModel
What does it look like for other models?
In [28]:
-- Plot the k correction as a function of redshift over the
-- range 0 to 2.
--
plotKZ :: Double -> Double -> Model -> Renderable ()
plotKZ elo ehi mdl = toRenderable (do
layout_title .= ("K correction: " ++ show elo ++ ":" ++ show ehi ++ " keV")
layout_x_axis . laxis_title .= "z"
layout_y_axis . laxis_title .= "k(z)"
plot (line "" [cs])
)
where
zs = [0, 0.02 .. 2]
cs = map (\z -> (z, kcorr elo ehi z mdl)) zs
The value of the K correction depends, in general, on the redshift of the object, the energy band being used, and the model parameters. It need not be a nice monotonic function, as I show here, when looking at a 3 keV gaussian line with a sigma of 1.2 keV for two different energy bands.
In [29]:
plotKZ 0.1 10 (gMdl 1 3 1.2)
plotKZ 1 5 (gMdl 1 3 1.2)
So, with this, I can finally calculate luminosities of sources, assuming I have their redshift and flux.
As I mentioned earlier, using a routine like iFModel
, which has a signature
of
iFModel :: Double -> Double -> (Double -> Double) -> Double
(where I have expanded out the Model
synonym)
requires the user to make sure that they understand what units
the parameters are. There is no way for the compiler to spot
when a user accidentally gives an energy in Joules rather than
keV, or a model that calculates an energy, rather than flux, density.
I started
these notebooks
by seeing how easy it was to track quantities and units in Physical
calculations, so let's try that approach here
by repeating the above steps but using types
which track dimensions.
I originally tried to do this with the
units package,
but ended up getting a bit stuck, so I ended up with
dimensional-tf
instead. Hopefully I can use the following to
guide me when I find time to try again with units
.
First I start with defining a bunch of types and functions;
they aren't always required$^\dagger$, but provide documentation and
simplify some code, and I based the naming scheme on the
existing dimensional-tf
code.
$^\dagger$ I actually developed the code without explicit types, checked what the compiler had inferred to make sure it made sense, and then created these types.
In [30]:
-- Flux is measured in photon/cm^2/s and the flux density in photon/cm^2/s/keV.
--
-- The Dim constructor arguments are for mass, length, time, and then others that
-- are not relevant here, so I can set the remaineders to 0 (for technical
-- reasons we have to use symbolic values such as Pos1, Neg1, and Zero,
-- rather than integers).
--
type DPhotonFluxDensity = P.Dim N.Neg4 N.Neg1 N.Pos1 N.Zero N.Zero N.Zero N.Zero
type PhotonFluxDensity = P.Quantity DPhotonFluxDensity Double
type DPhotonFlux = P.Dim N.Neg2 N.Zero N.Neg1 N.Zero N.Zero N.Zero N.Zero
type PhotonFlux = P.Quantity DPhotonFlux Double
-- I want to be able to easily refer to fluxes, so I create
-- units that can be used explicitly or in routines like
-- "fdens", defined below.
--
photonFluxDensity :: P.Unit DPhotonFluxDensity Double
photonFluxDensity = P.one P./ (P.centi P.meter P.* P.centi P.meter) P./ P.second P./ P.kilo electronVolt
photonFlux :: P.Unit DPhotonFlux Double
photonFlux = P.one P./ (P.centi P.meter P.* P.centi P.metre) P./ P.second
-- Helper routines to get values "into" and "out of" a particular unit.
-- A function x takes a value (Double, normal Haskell type) and converts it into that type
-- and unX takes the type and converts it back into a Double.
--
keV = (P.*~ P.kilo electronVolt)
unKeV = (P./~ P.kilo electronVolt)
fdens = (P.*~ photonFluxDensity)
unFdens = (P./~ photonFluxDensity)
pflux = (P.*~ photonFlux)
unPFlux = (P./~ photonFlux)
I am going to use the convention that names that end in P
are versions
of the "unit-less" code that have been converted to deal with
photon fluxes or photon flux densisites. So ModelP
is the
version of
type Model = Double -> Double
In [31]:
type Energy = P.Energy Double
type ModelP = Energy -> PhotonFluxDensity
The conversion of the power-law model is relatively easy, but there is extra effort required to deal with the units, in particular with operator precedence. As an example, defining
alpha = -gamma P.*~ P.one
results in a type error because -
has a lower precedence
than P.*~
; that is, the above is taken to mean
alpha = - (gamma P.*~ P.one)
and the error looked like
No instance for (Num (P.Quantity P.DOne Double)) arising from a use of syntactic negation
In the expression: - gamma P.*~ P.one
In an equation for ‘alpha’: alpha = - gamma P.*~ P.one
One thing I do like about the following is that it makes it explicit that the normalization refers to the value at 1 keV:
In [32]:
plMdlP :: PhotonFluxDensity -> Double -> ModelP
plMdlP s0 gamma e =
let alpha = (-gamma) P.*~ P.one
eterm = (e P./ keV 1) P.** alpha
in s0 P.* eterm
Evaluating the model returns both the numeric value and units:
In [33]:
plMdlP (fdens 1) 2 (keV 1.2)
Since I'm going to use this set of parameters, I can create an
instance of the model in the same way I created myModel
above:
In [34]:
myModelP = plMdlP (fdens 1) 2
:t myModelP
Values can be "extracted" from the result either explicitly:
In [35]:
myModelP (keV 1.2) P./~ photonFluxDensity
or by using routines such as unFdens
(this also shows
an example of using different units, in this case
eV rather than keV):
In [36]:
unFdens (myModelP (1200 P.*~ electronVolt))
As examples of the "type safety" that this approach brings, trying to use a wavelength as an argument causes an error:
In [37]:
myModelP (12 P.*~ P.nano P.meter)
I have to explicitly convert the length to an energy:
In [38]:
myModelP (hc P./ (12 P.*~ P.nano P.meter))
It's also an error to treat the output value incorrectly. As an example, if I try to calculate an energy flux (by trying to treat the answer as if it has units of W/m$^2$/s) then I also get a compile-time error:
In [39]:
myModelP (keV 1.2) P./~ (P.watt P./ P.meter P./ P.meter P./ P.second)
A version of the gaussian model can also be created:
In [40]:
gMdlP :: PhotonFluxDensity -> Energy -> Energy -> ModelP
gMdlP s0 e0 sigma e =
let x = (e P.- e0) P./ sigma
in s0 P.* P.exp ((-0.5) P.*~ P.one P.* x P.* x)
and used:
In [41]:
gMdlP (fdens 1) (keV 4.3) (keV 1.2) (keV 5.8)
gMdlP (fdens 1) (keV 4.3) (keV 1.2) (keV 5.8) P./~ photonFluxDensity
These models cal also be combined, although you have to remember to use the
"correct" mathematical operators (in this case, P.+
rather than just +
),
otherwise you get a type error:
In [42]:
cMdlP e = gMdlP (fdens 1.89) (keV 4.3) (keV 1.2) e P.+ plMdlP (fdens 0.023) 2 e
Unfortunately, because of the types, I can't use my plotModel
routine to display
this model. It's a relatively simple conversion, but this notebook is long enough
already, so I'll leave it for another time.
Given that I can create a model with units of
photon/cm$^2$/s/keV (photonFluxDensity
),
how do I integrate this to get photon/cm$^2$/s
(photonFlux
)? That is, I want a version of
iModel :: Double -> Double -> Model -> Double
iModel elo ehi mdl = result (absolute 1.0e-6 (simpson mdl elo ehi))
Since the integration code only works with Double
values,
i.e. it would be a type error to send simpson
a model with
type ModelP
,
the code will have to manually convert between types with dimensions
and Double
values:
In [43]:
iModelP :: Energy -> Energy -> ModelP -> PhotonFlux
iModelP elo ehi mdlP =
let x1 = unKeV elo
x2 = unKeV ehi
wmdl x = mdlP (keV x) P./~ photonFluxDensity
in iModel x1 x2 wmdl P.*~ photonFlux
So, does this work? The photon flux for myModel
over the range 0.1 to 1 keV was
9.9, so what do I get now?
In [44]:
iModelP (keV 0.1) (keV 10) myModelP
Yay - that's only a factor of $10^4$ out, which is quite close for Astronomy!
The reason for the discrepancy is that
the "native" units used by iModelP
are SI, so the output is
in photon/m$^2$/s rather than photon/cm$^2$/s,
which is what iModel
returns. This
explains the difference of $(100)^2$. Converting the
output to my units of choice, using photonFlux
,
returns the expected value of 9.9 (give or take a
little bit of floating-point accuracy):
In [45]:
unPFlux (iModelP (keV 0.1) (keV 10) myModelP)
Using iModelP I can calculate the photon flux, that is, the flux in photon/cm$^2$/s, but what I really need is the energy flux, i.e. erg/cm$^2$/s. First stop: more types and helper functions!
In [46]:
type DFluxDensity = P.Dim N.Neg2 N.Zero N.Neg1 N.Zero N.Zero N.Zero N.Zero
type FluxDensity = P.Quantity DFluxDensity Double
type DFlux = P.Dim N.Zero N.Pos1 N.Neg3 N.Zero N.Zero N.Zero N.Zero
type Flux = P.Quantity DFlux Double
energyFluxDensity :: P.Unit DFluxDensity Double
energyFluxDensity = erg P./ (P.centi P.meter P.* P.centi P.meter) P./ P.second P./ P.kilo electronVolt
energyFlux :: P.Unit DFlux Double
energyFlux = erg P./ (P.centi P.meter P.* P.centi P.metre) P./ P.second
Calculating the energy flux should look familiar:
In [47]:
iModelE :: Energy -> Energy -> ModelP -> Flux
iModelE elo ehi mdlP =
let x1 = unKeV elo
x2 = unKeV ehi
wmdl x = let e = keV x in e P.* mdlP e P./~ energyFluxDensity
in iModel x1 x2 wmdl P.*~ energyFlux
As a reminder, the photon flux was calculated with
iModelP :: Energy -> Energy -> ModelP -> PhotonFlux
iModelP elo ehi mdlP =
let x1 = unKeV elo
x2 = unKeV ehi
wmdl x = mdlP (keV x) P./~ photonFluxDensity
in iModel x1 x2 wmdl P.*~ photonFlux
Finally I can calculate the "energy flux" and compare with value from Sherpa:
In [48]:
iModelE (keV 0.1) (keV 10) myModelP
iModelE (keV 0.1) (keV 10) myModelP P./~ energyFlux
putStrLn "Sherpa:\n7.3796292733846489e-9"
So I agree with Sherpa, since
In [49]:
7.3796292733846489e-9 / 7.3833833576195705e-9
The K correction for a model can then be calculated using iModelE
:
In [50]:
-- The return value could be P.Dimensionless Double rather than
-- just Double, but that would be annoying to deal with in this
-- notebook, so return a Double.
--
kcorrE :: Energy -> Energy -> Double -> ModelP -> Double
kcorrE elo ehi z mdlP =
let obs = iModelE (elo P.* zp1) (ehi P.* zp1) mdlP
rest = iModelE elo ehi mdlP
zp1 = (1 + z) P.*~ P.one
in (rest P./ obs) P./~ P.one -- convert a Dimensionless Double to a Double
This is very similar to the original code:
kcorr :: Double -> Double -> Double -> Model -> Double
kcorr elo ehi z mdl =
let obs = iFModel (elo * zp1) (ehi * zp1) mdl
rest = iFModel elo ehi mdl
zp1 = 1 + z
in rest / obs
Does this version give the same answer as the original version?
In [51]:
cMdl e = gMdl 1 4.3 1.2 e + plMdl 0.01 2 e
cMdlP e = gMdlP (fdens 1) (keV 4.3) (keV 1.2) e P.+ plMdlP (fdens 0.01) 2 e
:t cMdl
:t cMdlP
In [52]:
kcorr 0.1 10 0.5 cMdl
kcorrE (keV 0.1) (keV 10) 0.5 cMdlP
So, yes it does (at least to the level of accuracy I care about). So, why bother with all this effort? The main one is because I can, but it can also reduce mistakes such as trying to integrate an invalid model: in the following example I try to calculate the K correction from a model defined as
$$S_1(E) + E S_2(E)$$(in this case $S_1$ is a power law and $S_2$ is a gaussian)
which the "unitless" version accepts, but the dimensional-tf version correctly
identifies as having incompatible units (that is, contrivedModelP
can not
even be created as the types do not match up):
In [53]:
contrivedModel e = m1 + m3
where
m1 = plMdl 1 2 e
m2 = gMdl 1 4.3 1.2 e
m3 = e * m2
:type contrivedModel
kcorr 0.1 10 0.5 contrivedModel
In [54]:
contrivedModelP e = m1 P.+ m3
where
m1 = plMdlP (fdens 1) 2 e
m2 = gMdlP (fdens 1) (keV 4.3) (keV 1.2) e
m3 = e P.* m2
Note that the equations for contrivedModel
and contrivedModelP
were broken
up into multiple lines to make the error message from contrivedModelP
somewhat
meaningful (and to make them a bit-more readable). If you are interested in
inscrutable Haskell errors, then try:
In [55]:
contrivedModelP e = plMdlP (fdens 1) 2 e P.+ e P.* gMdlP (fdens 1) (keV 4.3) (keV 1.2) e
which isn't a claim that the types do not match, but that the compiler has got itself stuck in some sort of loop trying to match up terms, and that it might work if given enough memory (or, at least, that's my understanding of the situation). By breaking up the definition over multiple lines the compiler is able to infer enough about the individual components that it can avoid this trap and error out.
To be "fair" to ghc, this is by no means the worst error message it can create, I just don't have a handy link to an example showing how deep the hole can go!
There you go; I hope you enjoyed it. If you have any questions, then please use the GitHub issues page or contact me on Twitter at https://twitter.com/doug_burke.