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".

K corrections

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.

Last time the notebook was run


In [1]:
import Data.Time
getCurrentTime


2015-03-04 01:01:57.459731 UTC

Spectral models

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

  1. in SI units

  2. 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


1.60217733e-9

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")


The wavelength corresponding to 1 keV is 1.2398413351660643e-9 m
  which is 1.2398413351660642 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


myModel :: Model

which can then be evaluated at one or more energies:


In [9]:
myModel 5
map myModel [0.1, 1, 10]


4.0e-2
[99.99999999999999,1.0,1.0e-2]

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]


4.0e-2
[99.99999999999999,1.0,1.0e-2]

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


plotModel :: forall x y. (PlotValue y, PlotValue x, Fractional x, Enum x) => x -> x -> (x -> y) -> Renderable ()

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


Calculating fluxes

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


9.899999999978695

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


9.9

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).

Calculating useful fluxes

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")


flux = 9.899999999978695     photon/cm^2/s
     = 7.3833833576195705e-9 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!

K corrections

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)


1.0000000000000018
1.5000000000000122

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:

  • depends on the slope and redshift, but not the energy range
  • the K correction will be 1 when $\gamma-2=0$, i.e. $\gamma=2$ as used in 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.

What would it mean to include dimensions and units?

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)


4.334379418815298e19 m^-4 kg^-1 s

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


myModelP :: ModelP

Values can be "extracted" from the result either explicitly:


In [35]:
myModelP (keV 1.2) P./~ photonFluxDensity


0.6944444444444445

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))


0.6944444444444445

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)


Couldn't match type ‘N.Z’ with ‘N.S N.Z’
Expected type: P.Unit P.DEnergy Double
Actual type: P.Unit P.DLength Double
In the first argument of ‘P.nano’, namely ‘P.meter’
In the second argument of ‘(P.*~)’, namely ‘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))


5.846819506755189e21 m^-4 kg^-1 s

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)


Couldn't match type ‘N.Z’ with ‘N.N N.Pos4’
Expected type: P.Unit DPhotonFluxDensity Double
Actual type: Numeric.Units.Dimensional.TF.Dimensional P.DUnit (P.Div (P.Dim N.Z (N.S N.Z) (N.N (N.S (N.S (N.S N.Z)))) N.Z N.Z N.Z N.Z) P.DTime) Double
In the second argument of ‘(P./~)’, namely ‘(P.watt P./ P.meter P./ P.meter P./ P.second)’
In the expression: myModelP (keV 1.2) P./~ (P.watt P./ P.meter P./ P.meter P./ P.second)
In an equation for ‘it’: it = 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


2.857569840734262e19 m^-4 kg^-1 s
0.4578333617716145

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


98999.99999978695 m^-2 s^-1

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)


9.899999999978695

Fluxes

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"


7.383383357619569e-12 kg s^-3
7.3833833576195705e-9
Sherpa:
7.3796292733846489e-9

So I agree with Sherpa, since


In [49]:
7.3796292733846489e-9 / 7.3833833576195705e-9


0.999491549598187

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


cMdl :: Double -> Double
cMdlP :: Energy -> Quantity DPhotonFluxDensity Double

In [52]:
kcorr 0.1 10 0.5 cMdl
kcorrE (keV 0.1) (keV 10) 0.5 cMdlP


0.7051545397316024
0.7051545397316025

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


contrivedModel :: Double -> Double
0.6949509062138369

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


Couldn't match type ‘N.Z’ with ‘N.S N.Pos1’
Expected type: P.Quantity DPhotonFluxDensity Double
Actual type: Numeric.Units.Dimensional.TF.Dimensional P.DQuantity (P.Dim (N.N (N.S (N.S N.Z))) N.Z (N.N (N.S N.Z)) N.Z N.Z N.Z N.Z) Double
In the second argument of ‘(P.+)’, namely ‘m3’
In the expression: m1 P.+ m3

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


Type function application stack overflow; size = 201
Use -ftype-function-depth=N to increase stack size to N
N.Add (N.Succ (N.N (N.S N.Z))) (N.Pred t0) ~ N.Pos1
In the second argument of ‘(P.+)’, namely ‘e P.* gMdlP (fdens 1) (keV 4.3) (keV 1.2) e’
In the expression: plMdlP (fdens 1) 2 e P.+ e P.* gMdlP (fdens 1) (keV 4.3) (keV 1.2) e
In an equation for ‘interactive:IHaskell1131.contrivedModelP’: interactive:IHaskell1131.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!

The end

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.