notebooks/k corrections with a dimensional twist.ipynb

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