AmmoniaLevelPopulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notes on the Ammonia model ##"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Exact equation from the [source code](https://github.com/pyspeckit/pyspeckit/blob/master/pyspeckit/spectrum/models/ammonia.py#L210):\n",
    "\n",
    "            population_upperstate = lin_ntot * orthoparafrac * partition/(Z.sum())\n",
    "\n",
    "            tau_dict[linename] = (population_upperstate /\n",
    "                                  (1. + np.exp(-h*frq/(kb*tkin) ))*ccms**2 /\n",
    "                                  (8*np.pi*frq**2) * aval *\n",
    "                                  (1-np.exp(-h*frq/(kb*tex))) /\n",
    "                                  (width/ckms*frq*np.sqrt(2*np.pi)) )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation}\n",
    "    \\tau = N_{tot} g_{opr} Z_{upper} \\frac{A_{ij} c^2}{8\\pi\\nu^2} \n",
    "    \\left(1-\\exp{ \\frac{-h \\nu}{k_B T_{ex}} } \\right)\n",
    "    \\left(1+\\exp{\\frac{-h \\nu}{k_B T_K}}\\right)\n",
    "    \\left((2\\pi)^{1/2} \\nu \\sigma_\\nu / c\\right)^{-1}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Equation 16 from Rosolowsky et al 2008:\n",
    "\n",
    "$$N(1,1) = \\frac{8 \\pi k \\nu_0^2}{h c^3} \\frac{1}{A_{1,1}} \\sqrt{2\\pi}\\sigma_\\nu (T_{ex}-T_{bg})\\tau$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rearranges to:\n",
    "\n",
    "$$\\tau = N(1,1) \\frac{h c^3}{8\\pi k \\nu_0^2} A_{1,1} \\frac{1}{\\sqrt{2 \\pi} \\sigma_\\nu} \\left(T_{ex}-T_{bg})\\right)^{-1}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Equation A4 of Friesen et al 2009:\n",
    "\n",
    "$$N(1,1) = \\frac{8\\pi\\nu^2}{c^2} \\frac{g_1}{g_2} \\frac{1}{A_{1,1}} \\frac{1+\\exp\\left(-h\\nu_0/k_B T_{ex}\\right)}{1-\\exp\\left(-h \\nu_0/k_B T_{ex}\\right)} \\int \\tau(\\nu) d\\nu$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Equation 98 of Mangum & Shirley 2015:\n",
    "\n",
    "\n",
    "$$N_{tot} = \\frac{3 h}{8 \\pi \\mu^2 R_i} \\frac{J_u(J_u+1)}{K^2} \n",
    "\\frac{Q_{rot}}{g_J g_K g_I} \\frac{\\exp{E_u/k_B T_{ex}}}{\\exp{h \\nu/k_B T_{ex}} - 1} \n",
    "\\left[\\frac{\\int T_R dv}{f\\left(J_\\nu(T_{ex})-J_\\nu{T_B}\\right) }\\right]$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## From Scratch ##"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\tau_\\nu = \\int \\alpha_\\nu ds$$\n",
    "$$\\alpha_\\nu = \\frac{c^2}{8\\pi\\nu_0^2} \\frac{g_u}{g_l} n_l A_{ul} \\left(1-\\frac{g_l n_u}{g_u n_l}\\right) \\phi_\\nu$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Excitation temperature:\n",
    "$$T_{ex} \\equiv \\frac{h\\nu_0/k_b}{\\ln \\frac{n_l g_u}{n_u g_l} } $$\n",
    "\n",
    "$\\nu_0$ = rest frequency of the line\n",
    "\n",
    "Rearranges to:\n",
    "$$ \\frac{n_l g_u}{n_u g_l} = \\exp\\left(\\frac{h \\nu_0}{k_B T_{ex}}\\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Boltzman  distribution:\n",
    "$$ \\frac{n_u}{n_l} = \\frac{g_u}{g_l} \\exp\\left(\\frac{-h \\nu_0}{k_B T}\\right)$$\n",
    "where T is a thermal equilibrium temperature"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rearranges to:\n",
    "$$ 1-\\frac{n_u g_l}{n_l g_u} = 1-\\exp\\left(\\frac{-h \\nu_0}{k_B T}\\right)$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Column Density \n",
    "$$N_u \\equiv \\int n_u ds$$\n",
    "$$N_l \\equiv \\int n_l ds$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Starting to substitute previous equations into each other:\n",
    "\n",
    "$$\\tau_\\nu d\\nu= \\alpha_\\nu d\\nu = \\frac{c^2}{8\\pi\\nu_0^2} \\frac{g_u}{g_l} n_l A_{ul} \\left(1-\\frac{g_l n_u}{g_u n_l}\\right) \\phi_\\nu d\\nu$$\n",
    "\n",
    "$$\\frac{g_u}{g_l}N_l = N_u\\exp\\left(\\frac{h \\nu_0}{k_B T_{ex}}\\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First substitution is the Boltzmann distribution, with $T_{ex}$ for T\n",
    "$$\\int \\tau_\\nu d\\nu = \\int \\frac{c^2}{8\\pi\\nu_0^2} \\frac{g_u}{g_l} n_l A_{ul} \\left[  1-\\exp\\left(\\frac{-h \\nu_0}{k_B T_{ex}}\\right) \\right] \\phi_\\nu d\\nu $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Second is the $N_l$ - $N_u$ relation:\n",
    "$$\\int \\tau_\\nu d\\nu = \\frac{c^2}{8\\pi\\nu_0^2}  A_{ul} N_u\\left[\\exp\\left(\\frac{h \\nu_0}{k_B T_{ex}}\\right)\\right] \\left[  1-\\exp\\left(\\frac{-h \\nu_0}{k_B T}\\right) \\right] \\int \\phi_\\nu d\\nu $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then some simplification:\n",
    "$$\\int \\tau_\\nu d\\nu = \\frac{c^2}{8\\pi\\nu_0^2}  A_{ul} N_u \\left[  \\exp\\left(\\frac{h \\nu_0}{k_B T}\\right) - 1 \\right] \\int \\phi_\\nu d\\nu $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$A_{ul} = \\frac{64\\pi^4\\nu_0^3}{3 h c^3} \\left|\\mu_{lu}\\right|^2$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Becomes, via some manipulation, equation 29 of Mangum & Shirley 2015:\n",
    "$$N_u = \\frac{3 h c}{8\\pi^3 \\nu \\left|\\mu_{lu}\\right|^2} \\left[\\exp\\left(\\frac{h\\nu}{k_B T_{ex}}\\right) -1\\right]^{-1} \\int \\tau_\\nu d\\nu$$\n",
    "where I have used $T_{ex}$ instead of $T$ here because that is one of the substitutions invoked (quietly) in their derivation.  There is some sleight-of-hand regarding assuming $N_l = n_l$ that essentially assumes $T_{ex}$ is constant along the line of sight, but that is fine.\n",
    "\n",
    "(Equation 30 is the same as this one, but with $dv$ instead of $d\\nu$ units)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve for tau again (because that's what's implemented in the code):\n",
    "\n",
    "$$\\mathrm{\"tau\"} = \\int \\tau_\\nu d\\nu = N_u \\frac{\n",
    "c^2 A_{ul}}{8\\pi\\nu_0^2} \\left[\\exp\\left(\\frac{h\\nu}{k_B T_{ex}}\\right) -1\\right] $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The key difference from Erik's derivation is that this is $N_u$, but he has defined $N_{(1,1)}= N_u + N_l$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, we get $N_l$ the same way as above:\n",
    "$$N_l = \\frac{8\\pi\\nu_0^2}{c^2} \\frac{g_l}{g_u} A_{ul}^{-1} \\left[  1-\\exp\\left(\\frac{-h \\nu_0}{k_B T_{ex}}\\right) \\right]^{-1} \\int \\tau d\\nu$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$N_l = \\frac{3 h c}{8 \\pi^3 \\nu \\left|\\mu_{lu}\\right|^2} \\frac{g_l}{g_u}  \\left[  1-\\exp\\left(\\frac{-h \\nu_0}{k_B T_{ex}}\\right) \\right]^{-1} \\int \\tau d\\nu$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Added together:\n",
    "$$N_u + N_l = \\frac{3 h c}{8 \\pi^3 \\nu \\left|\\mu_{lu}\\right|^2} \\frac{\\frac{g_l}{g_u}  +\\exp\\left(\\frac{-h \\nu_0}{k_B T_{ex}}\\right)}{1-\\exp\\left(\\frac{-h \\nu_0}{k_B T_{ex}}\\right)} \\int \\tau d\\nu$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can solve that back for tau, which is what Erik has done:\n",
    "\n",
    "$$\\int \\tau d\\nu = (N_u + N_l) \\frac{8 \\pi^3 \\nu \\left|\\mu_{lu}\\right|^2}{3 h c} \\frac{1-\\exp\\left(\\frac{-h \\nu_0}{k_B T_{ex}}\\right)} {\\frac{g_l}{g_u}  +\\exp\\left(\\frac{-h \\nu_0}{k_B T_{ex}}\\right)} $$\n",
    "\n",
    "$$=(N_u + N_l) \\frac{g_u}{g_l}\\frac{8 \\pi^3 \\nu \\left|\\mu_{lu}\\right|^2}{3 h c} \\frac{1-\\exp\\left(\\frac{-h \\nu_0}{k_B T_{ex}}\\right)} {1  +\\frac{g_u}{g_l}\\exp\\left(\\frac{-h \\nu_0}{k_B T_{ex}}\\right)} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$=(N_u + N_l) \\frac{g_u}{g_l}\\frac{A_{ul}c^2}{8\\pi\\nu_0^2} \\frac{1-\\exp\\left(\\frac{-h \\nu_0}{k_B T_{ex}}\\right)} {1  +\\frac{g_u}{g_l}\\exp\\left(\\frac{-h \\nu_0}{k_B T_{ex}}\\right)} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "now identical to Erik's equation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is actually a problem, because $N_u$ is related to $N_{tot}$ via the partition function, but there is some double-counting going on if we try to relate $N_{(1,1)}$ to $N_{tot}$ with the same equation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, to reformulate the equations in pyspeckit using the appropriate values, we want to use both the partition function (calculated using $T_{kin}$) and $N_u$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Eqn 31:\n",
    "    \n",
    "\n",
    "$$N_u = N_{tot} \\frac{g_u}{Q_{rot}} \\exp\\left(\\frac{-E_u}{k_B T_{kin}}\\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "is implemented correctly in pyspeckit:\n",
    "\n",
    "    population_upperstate = lin_ntot * orthoparafrac * partition/(Z.sum())\n",
    "\n",
    "where ``partition`` is \n",
    "\n",
    "$$Z_i(\\mathrm{para}) = (2J + 1) \\exp\\left[ \\frac{ -h (B_0 J (J+1) + (C_0-B_0)J^2)}{k_B T_{kin}}\\right]$$\n",
    "\n",
    "$$Z_i(\\mathrm{ortho}) = 2(2J + 1) \\exp\\left[ \\frac{ -h (B_0 J (J+1) + (C_0-B_0)J^2)}{k_B T_{kin}}\\right]$$\n",
    "\n",
    "\n",
    "...so I'm assuming (haven't checked) that $E_u = h (B_0 J (J+1) + (C_0-B_0)J^2)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the leading \"2\" above cancels out in the Z/sum(Z), so it doesn't matter if it's right or not.  I suspect, though, that the 2 belongs in front of both the para and ortho states, but it should be excluded for the J=0 case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Erik's diff here)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "EWR: The above equation is problematic because it relates the total column density to the $(J,J)$ state which is the equivalent of the $N_{(1,1)}$ term.  In the notation above $N_{(1,1)} = N_u + N_l$, so to get this right, you need to consider the inversion transition splitting on top of the total energy of the state so that \n",
    "$$ E_u = h (B_0 J (J+1) + (C_0-B_0)J^2) + \\Delta E_{\\mathrm{inv}}, g_u = 1 $$ \n",
    "and \n",
    "$$ E_l = h (B_0 J (J+1) + (C_0-B_0)J^2) - \\Delta E_{\\mathrm{inv}}, g_l = 1 $$ \n",
    "or, since the splitting is small compared to the rotational energy (1 K compared to > 20 K), then\n",
    "\n",
    "$$Z_J \\approx 2 (2J + 1) \\exp\\left[ \\frac{ -h (B_0 J (J+1) + (C_0-B_0)J^2)}{k_B T_{\\mathrm{rot}}}\\right]$$\n",
    "\n",
    "where the leading 2 accounts for the internal inversion states.  Since this 2 appears in all the terms, it cancels out in the sum.  Note that I have also changed the $T_{\\mathrm{kin}}$ to $T_{\\mathrm{rot}}$ since these two aren't the same and it is the latter which establishes the level populations.\n",
    "\n",
    "Returning to the above, I would then suggest \n",
    "\n",
    "$$N_{(J,J)} = N_{tot} \\frac{Z_J}{\\sum_j Z_j} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some numerical checks: How bad was the use of Tkin instead of Tex in the $\\tau$ equation?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$(N_u + N_l) \\frac{g_u}{g_l}\\frac{A_{ul}c^2}{8\\pi\\nu_0^2} \\frac{1-\\exp\\left(\\frac{-h \\nu_0}{k_B T_{ex}}\\right)} {1  +\\frac{g_u}{g_l}\\exp\\left(\\frac{-h \\nu_0}{k_B T_{ex}}\\right)} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from astropy import units as u\n",
    "from astropy import constants\n",
    "freq = 23*u.GHz\n",
    "def tau_wrong(tkin, tex):\n",
    "    return (1-np.exp(-constants.h * freq/(constants.k_B*tkin)))/(1+np.exp(-constants.h * freq/(constants.k_B*tex)))\n",
    "def tau_right(tex):\n",
    "    return (1-np.exp(-constants.h * freq/(constants.k_B*tex)))/(1+np.exp(-constants.h * freq/(constants.k_B*tex)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tkin = np.linspace(5,40,101)*u.K\n",
    "tex = np.linspace(5,40,100)*u.K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "grid = np.array([[tau_wrong(tk,tx)/tau_right(tx) for tx in tex] for tk in tkin])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import pylab as pl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.imshow(grid, cmap='hot', extent=[5,40,5,40])\n",
    "pl.xlabel(\"Tex\")\n",
    "pl.ylabel(\"Tkin\")\n",
    "pl.colorbar()\n",
    "pl.contour(tex, tkin, grid, levels=[0.75,1,1/0.75], colors=['w','w','k'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the error could be 50%-700% over a somewhat reasonable range.  That's bad, and it affects the temperature estimates.  However, the effect on temperature estimates should be pretty small, since each line will be affected in the same way.  The biggest effect will be on the column density."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But, is this error at all balanced by the double-counting problem?\n",
    "\n",
    "Because we were using the partition function directly, it's not obvious.  I was assuming that we were using the equation with $N_u$ as the leader, but we were using $N_u+N_l$.  i.e., I was using this equation:\n",
    "\n",
    "$$\\int \\tau d\\nu =(N_u + N_l) \\frac{g_u}{g_l}\\frac{A_{ul}c^2}{8\\pi\\nu_0^2} \\frac{1-\\exp\\left(\\frac{-h \\nu_0}{k_B T_{ex}}\\right)} {1  +\\frac{g_u}{g_l}\\exp\\left(\\frac{-h \\nu_0}{k_B T_{ex}}\\right)} $$\n",
    "\n",
    "but with $N_u$ in place of $N_u + N_l$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The magnitude of the error can therefore be estimated by computing $(N_u+N_l)/N_u = 1 + \\frac{N_l}{N_u}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the Boltzmann distribution to compute this error, then:\n",
    "$$ \\frac{n_u}{n_l} = \\frac{g_u}{g_l}\\exp\\left(\\frac{-h \\nu_0}{k_B T}\\right)$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def nunlnu_error(Tkin):\n",
    "    return 1+np.exp(-constants.h * freq / (constants.k_B * Tkin))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.plot(tkin.value, nunlnu_error(tkin))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we were always off by a factor very close to 2.  The *relative* values of $\\tau$ should never have been affected by this issue."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It will be more work to determine exactly how much the T_K and column estimates were affected."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# New work in May 2016: T_{rot} #"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Comparing Trot and Tkin.  If we start with the equation that governs level populations,\n",
    "\n",
    "$$N_u = N_{tot} \\frac{g_u}{Q_{rot}} \\exp\\left(\\frac{-E_u}{k_B T_{kin}}\\right)$$\n",
    "\n",
    "we get \n",
    "\n",
    "$$N_u / N_l = \\frac{g_u}{g_l} \\exp\\left(\\frac{-E_u}{k_B T_{kin}} + \\frac{E_l}{k_B T_{kin}}\\right)$$\n",
    "\n",
    "where we really mean $T_{rot}$ instead of $T_{kin}$ here as long as we're talking about just two levels.  This gives us a definition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$T_{rot} = \\left(\\frac{E_l-E_u}{k_B}\\right)\\left[\\ln\\left(\\frac{N_u g_l}{N_l g_u}\\right)\\right]^{-1}$$\n",
    "\n",
    "which is the rotational temperature for a two-level system... which is just a $T_{ex}$, but governing non-radiatively-coupled levels."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, for example, if we want to know $T_{rot}$ for the 2-2 and 1-1 lines at $n=10^4$ and $T_{kin}=20$ K:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/adam/repos/pyradex/pyradex/core.py:852: RuntimeWarning: invalid value encountered in true_divide\n",
      "  frac_level_diff = level_diff/self.level_population\n"
     ]
    }
   ],
   "source": [
    "from pyradex import Radex\n",
    "from astropy import constants, units as u\n",
    "R = Radex(species='p-nh3', column=1e13, collider_densities={'pH2':1e4}, temperature=20)\n",
    "tbl = R(collider_densities={'ph2': 1e4}, temperature=20, column=1e13)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<Table length=2>\n",
       "<table id=\"table6188359864\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>Tex</th><th>tau</th><th>frequency</th><th>upperstateenergy</th><th>upperlevel</th><th>lowerlevel</th><th>upperlevelpop</th><th>lowerlevelpop</th><th>brightness</th><th>T_B</th></tr></thead>\n",
       "<thead><tr><th>K</th><th></th><th>GHz</th><th>K</th><th></th><th></th><th></th><th></th><th>erg / (cm2 Hz s sr)</th><th>K</th></tr></thead>\n",
       "<thead><tr><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>bytes6</th><th>bytes6</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th></tr></thead>\n",
       "<tr><td>6.78936052483</td><td>0.0929414309528</td><td>23.6944955</td><td>1.14</td><td>01_01_</td><td>01_01_</td><td>0.393279209876</td><td>0.464987835394</td><td>6.17629302965e-17</td><td>0.358064405386</td></tr>\n",
       "<tr><td>6.53340348878</td><td>0.0212672724409</td><td>23.7226333</td><td>42.32</td><td>02_02_</td><td>02_02_</td><td>0.064641582585</td><td>0.0769470693495</td><td>1.37491511187e-17</td><td>0.0795203502905</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=2>\n",
       "     Tex            tau       frequency  ...      brightness           T_B      \n",
       "      K                          GHz     ... erg / (cm2 Hz s sr)        K       \n",
       "   float64        float64      float64   ...       float64           float64    \n",
       "------------- --------------- ---------- ... ------------------- ---------------\n",
       "6.78936052483 0.0929414309528 23.6944955 ...   6.17629302965e-17  0.358064405386\n",
       "6.53340348878 0.0212672724409 23.7226333 ...   1.37491511187e-17 0.0795203502905"
      ]
     },
     "execution_count": 97,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tbl[8:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# we're comparing the upper states since these are the ones that are emitting photons\n",
    "trot = (u.Quantity(tbl['upperstateenergy'][8]-tbl['upperstateenergy'][9], u.K) *\n",
    "        np.log((tbl['upperlevelpop'][9] * R.statistical_weight[8]) /\n",
    "               (tbl['upperlevelpop'][8] * R.statistical_weight[9]))**-1\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$22.806041 \\; \\mathrm{K}$"
      ],
      "text/plain": [
       "<Quantity 22.806041117736964 K>"
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.6613820068044447"
      ]
     },
     "execution_count": 104,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tbl['Tex'][8:10].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pause here ##\n",
    "\n",
    "$T_{rot} = 60$ K for $T_{kin}=25$ K?  That doesn't seem right.  Is it possible RADEX is doing something funny with level populations?\n",
    "\n",
    "# ERIK I SOLVED IT #\n",
    "I had left out the $^{-1}$ in the code.  Oops!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "delta-T for 1-1_upper - 1-1_lower: -1.137156499241892 K\n",
      "Excitation temperature computed is 6.789341498213386 K and should be 6.789360524825584\n"
     ]
    }
   ],
   "source": [
    "dT_oneone = -(constants.h * u.Quantity(tbl['frequency'][8], u.GHz)/constants.k_B).to(u.K)\n",
    "print(\"delta-T for 1-1_upper - 1-1_lower: {0}\".format(dT_oneone))\n",
    "tex = (dT_oneone *\n",
    "        np.log((tbl['upperlevelpop'][8] * R.statistical_weight[8]) /\n",
    "               (tbl['lowerlevelpop'][8] * R.statistical_weight[8]))**-1\n",
    "        )\n",
    "print(\"Excitation temperature computed is {0} and should be {1}\".format(tex.to(u.K), tbl['Tex'][8]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Moving on: comparison to Swift et al 2005 ##"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Swift et al 2005 eqn A6\n",
    "\n",
    "$$T_R = T_K \\left[ 1 + \\frac{T_K}{T_0} \\ln \\left[1+0.6\\exp\\left( -15.7/T_K \\right)\\right] \\right]^{-1}$$\n",
    "\n",
    "where $T_0=-41.18$ K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-41.18"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T0=tbl['upperstateenergy'][8]-tbl['upperstateenergy'][9]\n",
    "T0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def tr_swift(tk):\n",
    "    return tk*(1+tk/T0 * np.log(1+0.6*np.exp(-15.7/tk)))**-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "22.662533151044173"
      ]
     },
     "execution_count": 113,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tr_swift(20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def trot_radex(column=1e13, density=1e4, tkin=20):\n",
    "    tbl = R(collider_densities={'ph2': density}, temperature=tkin, column=column)\n",
    "    trot = (u.Quantity(tbl['upperstateenergy'][8]-tbl['upperstateenergy'][9], u.K) *\n",
    "        np.log((tbl['upperlevelpop'][9] * R.statistical_weight[8]) /\n",
    "               (tbl['upperlevelpop'][8] * R.statistical_weight[9]))**-1\n",
    "        )\n",
    "    return trot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import pylab as pl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x1794621d0>"
      ]
     },
     "execution_count": 116,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAETCAYAAADtZdsKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYFNXZ9/HvzaYRgrgDQlDcoyyKYBQ084oLKopBwYgb\n7gmXYPSJ0YS8AWPyuCQuaKJhcQEDbhhBNBhFHJAoiig4Iu6KEgUFBVERkbmfP05NbMeerZeq7unf\n57rqmu6qOtV3U87cnqXOMXdHREQkDk2SDkBEREqHko6IiMRGSUdERGKjpCMiIrFR0hERkdgo6YiI\nSGxiTzpm1sHMZpvZEjOrMLPh0f5rzGypmS0ys/vNrHUN5S8ys5fM7EUzm2xmLaL9o8xsuZk9H239\n4vxeIiJSN4v7OR0zawu0dfdFZtYKWAgMADoAs9290syuAtzdf12tbHtgHrCnu39lZvcAD7v7JDMb\nBaxz9+ti/UIiIlJvsdd03H2Fuy+KXn8GLAV2dPdZ7l4ZnTafkITSaQq0NLNmwBbA+ynHLE9hi4hI\nDiTap2NmOwHdgWeqHToLmFn9fHd/H7gWeBf4D7DG3WelnHJB1Dw3wcy2zEvQIiKSsdib1/77waFp\nrRy4wt2np+wfCezn7iekKdMGuB8YBKwFpgL3ufsUM9sOWOXubmZ/ANq5+9lprqF5f0REMuDuWbcm\nJVLTiZrGpgJ3Vks4Q4GjgSE1FD0MeMvdP3b3TcA/gIMA3P0j/yaDjgd61vT57p74NmrUqMSv1ZBy\n9Tm3rnNqOt6Q/bn8dyuEe1cs96+hxwr13hXj/SuU371cSap57TbgZXcfU7UjGm12CXCcu2+oody7\nwI/MbHMzM6AvoU+oaoBClYHAS3mJPEfKysoSv1ZDytXn3LrOqel4Q/cnLddxFcP9a+ixQr13UHz3\nr7H97iUxeq03MBeoADzaRgI3Ai2A1dGp8919mJm1A8a7e/+o/Cjgp8BG4AXgHHffaGaTCP1DlcA7\nwPnuvjLN53vc31lyZ/To0YwePTrpMCQDunfFzczwHDSvJdankxQlneJWXl5e0P8XLTXTvStuSjoZ\nUtIREWm4XCUdTYMjIiKxUdIREZHYKOmIiEhslHRERCQ2SjoiIhIbJR0REYmNko6IiMRGSUdERGKj\npCMiIrFR0hERkdgo6YiISGyUdEREJDZKOiIiEhslHRERiY2SjoiIxKZZ0gEk4Xe/a9j5+Vh+Z7vt\n4Oc/h+bNc39tEZFCVZJJp1kG39qyXrro22bMgAcegPvug223ze21RUQKlVYOTcimTfDb38Ldd8O0\nadCtW9IRiYjUTMtVZ6hQkk6Vu+6CESPg5pth0KCkoxERSU9JJ0OFlnQAnn8efvITOP10uPxyaKLh\nHSJSYJR0MlSISQfgww/hxBOhTRv4+9+hdeukIxIR+YaSToYKNekAfPUVXHhhSDo77wydOn1323nn\nMPIt1wMbRERqo6SToUJOOlU++QSWLfv29s474efbb4fktOuusMsu3/65++7Qvr0SkojknpJOhooh\n6dTlk0/gzTfD9sYb32yvvgobNsCee4Ztr72+2XbZBZo2TTpyESlWSjoZagxJpzarV8Mrr8DSpWF7\n5RV4+eXQZ/TDH4ah2VVb166hD0lEpC5KOhlq7EmnJp9+ChUVsHgxvPhi+FlRER5M3X9/6NkTevWC\nHj00iEFEvktJJ0OlmnTSqawMzXLPPQfPPgsLFsCiRWHAQq9eYevTB/bZR8O4RUqdkk6GlHRqt3Ej\nLFkSktD8+fDkk6HJrndvOPjgsPXoAS1aJB2piMRJSSdDSjoN98EHMG9eSEBPPgmvvw4HHABHHBG2\nbt1UExJp7JR0MqSkk701a2DuXHjsMXj00fD+8MNDAjr8cGjXLukIRSTXlHQypKSTe++8800Cevzx\n0Cd0/PFh69pVzw2JNAZKOhlS0smvTZvg6afDzNkPPBAGK1QloN69M1tWQkSSp6STISWd+LiHQQnT\npoVt2bKQfE49NQxIqKsfaPVqmDo11KJ694aTTgozLohI/JR0MqSkk5x334V77oHJk+Hjj+Hkk0MC\n6tLlm3M++wymT4cpU+Df/4Z+/eDII8MAhmnToHv3UO6EE2DrrZP7LiKlRkknQ0o6heGll0LymTw5\nzIpw4olh5oRHHgnPBp18MgwYAK1afVPmyy9h5syQkB59FA45BE45JdSeNt88ue8iUgqKNumYWQdg\nErADUAmMc/ebzOwa4FhgA/AmcKa7f5qm/EXA2VHZiui8r8xsK+AeoBPwDjDY3demKa+kU0AqK8Nw\n7H/8I8wRd+KJsM02dZf79NNQI5o0KTzQesopcO65sPfe+Y9ZpBQVc9JpC7R190Vm1gpYCAwAOgCz\n3b3SzK4C3N1/Xa1se2AesGeUaO4BHnb3SWZ2NbDa3a8xs0uBrdz9sjSfr6TTyLz1Ftx6K9x+O+y0\nE5x3HgweDFtskXRkIo1HrpJO7I/0ufsKd18Uvf4MWArs6O6z3L0yOm0+IQml0xRoaWbNgC2A/0T7\nBwATo9cTgePzEb8Uns6d4Y9/DH1Gl14aBh906AAXXBAeZBWRwpHoc+RmthPQHXim2qGzgJnVz3f3\n94FrgXcJyWaNuz8eHd7e3VdG560Ats9P1FKomjUL/UAPPRQmNN1ySzjooNDnM3duGE0nIslK7KmJ\nqGltKnBhVOOp2j8S2OjuU9KUaUOo0XQC1gJTzWxIunOBGv/EjB49+r+vy8rKKCsry/BbSKHq2DHU\nfkaOhIkT4ZxzQhK6+OLQb9S8edIRihS28vJyysvLc37dREavRU1jDwEz3X1Myv6hwLnAoe6+IU25\nE4Ej3f3c6P1pwAHufoGZLQXK3H1l1G/0hLvvleYa6tMpQZWV8PDDcN11YfG7iy6C889Xv49IfRVt\nn07kNuDlagmnH3AJcFy6hBN5F/iRmW1uZgb0JfQJATwIDI1enwFMz0fgUpyaNIFjj4UnnggzJcyb\nF5b4vuEGWL8+6ehESkcSo9d6A3MJw5092kYCNwItgNXRqfPdfZiZtQPGu3v/qPwo4KfARuAF4Bx3\n32hmWwP3Ah2BZYQh02vSfL5qOgKEoda//31YwuHSS8Oot+99L+moRApT0Q6ZTpqSjlT3wgtw+eVh\nEbuq5KOHTUW+rdib10QKxr77hil2ZsyAWbNgzz3h7rs12k0kH1TTEammvDyMcttsM7j+evjRj5KO\nSCR5qumI5ElZGTz3XBjddsIJMGRImCFbRLKnpCOSRpMmMHQovPYa7L477Lcf/OY3YRZsEcmcko5I\nLVq2hNGj4cUXwzQ7e+8dnvcRkcyoT0ekAWbNgp/9LNR8xoyBdu2SjkgkHurTEUnAYYdBRUV4sLRr\nVxg7Nsx2ICL1o5qOSIYqKsIzPU2bhuSjtXykMVNNRyRhXbqE6XSGDAkj3q6+GjZtSjoqkcKmmo5I\nDixbBqefHl5PmgSdOiUbj0iuqaYjUkA6dYLZs+Hoo2H//WHyZM1oIJKOajoiOfbCC3DKKdCtG9x8\nM2y1VdIRiWRPNR2RArXvvrBwIWy3XUg8TzyRdEQihUM1HZE8+te/4MwzYdiwMKNBE/1vnhQpLW2Q\nISUdidv778OgQbD11mGQgZrbpBipeU2kSLRvH5rYOncOgwwWL046IpHkKOmIxKBFizBtzh/+EGY1\nmDQp6YhEkqHmNZGYLVkCAwfCoYfCDTeEdXtECp2a10SK1N57h6WxP/wwzGTw4YdJRyQSHyUdkQS0\nbg1Tp8IRR4SVSV9+OemIROKh5jWRhE2aBL/8Jdx1F/Ttm3Q0IumpeU2kkTj9dLj33jBx6G23JR2N\nSH6ppiNSIF59FY45Bk46Ca64Qg+SSmHRw6EZUtKRQvbRR3D88dCxI9xxB2y+edIRiQRqXhNphLbb\nDh5/PLw+4ghYuzbZeERyTUlHpMBsvjlMmRIWievbF1atSjoikdypV/OamTUDBgEHRrtaApuAL4AX\ngSnu/mW+gswlNa9JsXCHkSNh+nR49FHYccekI5JSFlufjpn1BA4GHnP3ijTHdwGOARa7+5xsA8o3\nJR0pNldfDWPHwqxZYf42kSTEmXS6pEs2ac7rDCx396+yDSqflHSkGN1yC/zxj2GphL33TjoaKUUF\nMXrNzLZw9y+yDSJOSjpSrCZPhv/5H3jooTBbtUicCmX02nnZBiAi9XPKKaGZ7eijYd68pKMRyUx9\nmteuAw4BPq1+CNjT3dvlKba8UE1Hit1jj4UE9OCDYd42kTjE2adjwC/c/fo0x37h7jdkG0SclHSk\nMfjnP2Ho0PBTTW0Sh1j7dMxsK3f/JM3+lu7+ebZBxElJRxqL6dPhvPPC4ILu3ZOORhq7XCWdZvU5\nKV3CifYXVcIRaUwGDICNG+Goo0KT2z77JB2RSN1in5HAzDqY2WwzW2JmFWY2PNp/jZktNbNFZna/\nmbVOU3Z3M3vBzJ6Pfq41sxHRsVFmtjw69ryZ9Yv7u4nE7cQT4dprw5Q5r7ySdDQidctoyLSZneDu\n92f0gWZtgbbuvsjMWgELgQFAB2C2u1ea2VWAu/uva7lOE2A50Mvdl5vZKGCdu19Xx+ereU0anYkT\nw+wFTzwBu+2WdDTSGMXavJZGq0w/0N1XACui15+Z2VJgR3eflXLafOCEOi51GPCmuy9P2Zf1P4hI\nMTrjDPjqqzBX25NPQqdOSUckkl6mzWs5qSqY2U5Ad+CZaofOAmbWUfwk4K5q+y6ImucmmNmWuYhR\npFice254ePTIIzVJqBSuxGaZjprWpgIXuvtnKftHAhvdfUotZZsDxwH3pey+Gejs7t0JNalam9lE\nGqMLLwzr8fTvD59rmI8UoEyb17ISzVo9FbjT3aen7B8KHA0cWscljgIWuvtHVTtSXwPjgRk1FR49\nevR/X5eVlVFWVlb/4EUK3JVXhmd4Bg+GadOgefOkI5JiVF5eTnl5ec6vm+lAgtPdfVLGH2o2CVjl\n7hen7OsHXAsc4u6r6yh/F/CIu09M2dc26i/CzC4Cerr7kDRlNZBAGr2NG8OQ6rZt4dZbwdTbKVlK\ndMLPbCb6NLPewFyggtA35MBI4EagBVCVcOa7+zAzaweMd/f+VZ8NLCM0pa1Lue4kQv9QJfAOcL67\nr0zz+Uo6UhI+/xwOPRQOOyzMUC2SjYKYZboYKelIKVm1Cnr3hgsugOHDk45GilnSQ6ZFpAhsuy08\n8ggcfDDssEPo5xFJUlZJpxjX0xEpNTvvDA8/DIcfDu3bQ58+SUckpUzr6YiUgG7dYNIkGDQI3n47\n6WiklGk9HZESctNNYSG4p56C1t+Z3VCkZlpPJ0NKOlLK3MOggrffDovANVOvrtST1tPJkJKOlLqv\nvw5LXv/wh3BDUf0voyQpV0mn1j4dM9vMzLapz3o6ZtYx22BEJP+aNYN77w2Lv40dm3Q0UmpqrVy7\n+wYzO9zMvg9Mc/f11c8xszbAYOBl4L38hCkiudSmDcyYEUay7bprmJ1aJA71bV5rS5j5eXtgc0Ky\n2gR8QVjTZoK7r81jnDmj5jWRb8yZE57dmTsX9tgj6WikkGlGggwp6Yh82223wVVXwbPPhhqQSDqx\nJB0z2wPY5O5vZPtBhUJJR+S7RowII9qmT4cmiS14IoUsrqTTDCgD9iA0py1w94XZfmiSlHREvmvj\nxtCvc9hh8LvfJR2NFKJEmtfMrBfQg/Bg6KtAubtvyjaIOCnpiKS3YgXsv38Y0XbMMUlHI4Um8T6d\nqOmtDGgO/Af4VzHMw6akI1Kzp58OK4/++99hVJtIlcSTTrVg2gEHu/u9WV8sz5R0RGp3yy1he/pp\naNky6WikUBRU0ikmSjoitXOHs86CL7+EKVO06qgEcU+D0wwYBBwY7WrJN8/pvAhMcfcvsw0mDko6\nInVbvz48OHrqqXDRRUlHI4Ugzgk/exJmmX7M3V9Mc3wX4BhgsbvPyTagfFPSEamfZcvggAPg7ruh\nrCzpaCRpcSadLu5eUY+AOgPL3f2rbIPKJyUdkfqbNQtOPx2efx7atk06GklS3M1rvYAOwDPu/p9o\nX19ghbsvyTaIOCnpiDTM5ZeH6XIeewyaNk06GklKnDWdK4A9gbeArsBsd/9T1M+zwt23zTaIOCnp\niDTMpk1wxBFw8MEwenTS0UhScpV06rOE0xp3H5TywX3M7DLgGkB/vUUauaZNYfJk6NEjJB7NSC3Z\nqM8sS1+a2VZm9jMz+567zwP+BpxH/ZKWiBS5tm3hzjvhtNPCzAUimapP0hkH9APaApUA7r7G3f8G\n/CqPsYlIATn0UDjvPBgyJDS5iWRCD4eKSL1t2gSHHw6HHKL+nVKjGQkypKQjkp0VK2C//UJzm/p3\nSkeukk5GK2eY2QnZfrCIFCf170g2Ml2uqVVOoxCRotK3b+jfOfVUqKxMOhopJpkmHbVPiZS43/4W\nNmyAP/856UikmGhhWhHJSLNm8Pe/h6Tz3HNJRyPFQklHRDLWqRP85S9w8smwbl3S0UgxUNIRkawM\nHhyGUI8YkXQkUgwyTTpTcxqFiBS1MWPgqafCMggitdFzOiKSEwsXQr9+sGAB7LRT0tFIriX6nI6I\nSHU9esCll4Zpcr7+OulopFDVK+mYWS8zG2hmO6bs62tme+cvNBEpNhdfDK1awRVXJB2JFKrY19Mx\nsw7AJGAHwgSi49z9JjO7BjgW2AC8CZzp7p9WK7s7cA/hOSEDOgP/391vNLOtomOdgHeAwe6+Ns3n\nq3lNJI8++AD23Rfuuy8shSCNQ5yLuP2Pu1+b8r4P0Iewns5Kd9+uQR9o1hZo6+6LzKwVsBAYQFiZ\ndLa7V5rZVYC7+69ruU4TYDnQy92Xm9nVwGp3v8bMLgW2cvfL0pRT0hHJsxkzwmi2RYtgyy2TjkZy\nIc4+nZyup+PuK9x9UfT6M2ApsKO7z3L3qgk15hOSUG0OA9509+XR+wHAxOj1ROD4hsYmIrlx7LFh\ntVENo5bqEl1Px8x2AroDz1Q7dBYws47iJwF3pbzf3t1XRvGtALbPJjYRyc6114Zh1FP1gIWkqLWm\nYmabAa3c/a50x919fMq5Hd39vfp+cNS0NhW4MKrxVO0fCWx09ym1lG0OHAd8p/ksNbyaDoxOWQik\nrKyMsrKy+oYtIvXUqlWYJue44+Cgg6B9+6QjkoYoLy+nvLw859etT59Of+D7wDR3X5/meBtgMPBy\n1PRW94eGQQgPATPdfUzK/qHAucCh7r6hlvLHAcPcvV/KvqVAmbuvjPqNnnD3vdKUVZ+OSIxGj4an\nn4aZM6GJHtIoWrEu4hb9ET+L0GS1OaGGtAn4gtCZPyHdSLFarjcJWOXuF6fs6wdcCxzi7qvrKH8X\n8Ii7T0zZdzXwsbtfrYEEIoVj40bo0ycsgzB8eNLRSKaKduVQM+sNzAUqCE1gDowEbgRaAFUJZ767\nDzOzdsB4d+8fld8CWAZ0dvd1KdfdGrgX6BgdH+zua9J8vpKOSMxefx0OPBDmzoUf/jDpaCQTiScd\nM2sNfA40cfeN2QYSFyUdkWSMGwe33ALPPAMtWiQdjTRUYknHzHoBRxNqKHcAHdz939kGEhclHZFk\nuIdBBfvsA1demXQ00lCFUNM5CtgW2MzdJ2QbSFyUdESS8+GH0K0b3HuvZisoNoWQdLYhdNa/kW0Q\ncVLSEUnWgw/CL34BixfD97+fdDRSX4WQdH4ONAU+cffJ2QYSFyUdkeSdfXYYPj1+fN3nSmEohKUN\nPnL3vxAm5xQRqbfrr4dZs8IcbVJaskk6z5rZGMIM1CIi9da6NUycCOefDx99lHQ0Eqf6zEjQJGUi\nztT9v3D3G/IWWZ6oeU2kcFxyCbz1VpifzbJuuJF8irN57VkzG2RmXart72hmZ6Qu7CYi0hBXXAGv\nvRbmaJPSUJ+azvnuPjbN/kMIi6X1rmlC0EKkmo5IYVm0KCyD8Nxz8IMfJB2N1CTOmk7PaMXOb3H3\nucD6Yko4IlJ4uneHiy6CoUOh8jsN+dLY1CfpdAJuNLO3zOxRM/uDmR1rZjsAJ+Y5PhEpAb/6FWzY\nADfdlHQkkm/1aV47z93HRa93BXpFW0+gq7sX1eNdal4TKUxvvBEmBZ0zR5OCFqLYHg41s5nAcekm\n9TSzX7r7n7MNIk5KOiKFa+zY8MDo009D8+ZJRyOp4kw6uwH7Aq+5+6Jqx7q5++Jsg4iTko5I4XKH\n/v2hRw/4/e+TjkZSJT4NTrFS0hEpbB98EAYXPPggHHBA0tFIlUKYBkdEJOfatYO//AVOOw2++CLp\naCTXVNMRkYJ06qnQpk1IQJI8Na9lSElHpDisWQNdu8KECeHhUUmWmtdEpFFr0wZuvz0sg/Dxx0lH\nI7mimo6IFLQLLwwrjt6luU8SpZqOiJSEK6+EF16Au+9OOhLJBdV0RKTgLVgAxxwDzz8PHTokHU1p\nUk1HREpGz54wYoQmBW0MlHREpChcdhmsXw9jxiQdiWRDzWsiUjTeeivMUjB7NnSpvqyk5JWa10Sk\n5HTuDH/6EwwZAl9+mXQ0kgnVdESkqLjDoEFhldHrrks6mtKhGQkypKQjUvxWr4Zu3eCOO+Cww5KO\npjSoeU1EStY228Btt8GZZ2q2gmJTkjWdpGMQkVy5HtgRGJx0ICUhFzWdZrkIpNiUWqIVaazWrw/P\n8FxyiXPGGUlH07iZZZ1vwnVK7Q+w+nREGpfFi0O/zlNPwW67JR1N46U+HRERwoCCUaPg5JPhq6+S\njkbqopqOiBQ9d/jJT2CXXeDaa5OOpnHSkOkMKemINE6rV8O++8Lf/gZHH510NI2Pkk6GlHREGq+5\nc+Gkk2DhQmjfPuloGpei7dMxsw5mNtvMlphZhZkNj/ZfY2ZLzWyRmd1vZq1rKL+lmd0XnbvEzA6I\n9o8ys+Vm9ny09Yvze4lI8g45BH72MzjtNNi0KeloJJ3Yazpm1hZo6+6LzKwVsBAYAHQAZrt7pZld\nBbi7/zpN+TuAOe5+u5k1A7Zw90/NbBSwzt1rnRhDNR2Rxm3TJjj0UDjySPjNb5KOpvEo2pqOu69w\n90XR68+ApcCO7j7L3atWyphPSELfEtV+Dnb326PyX7v7p6mn5Dd6ESl0TZvC5MlhCYSnnko6Gqku\n0SHTZrYT0B14ptqhs4CZaYrsDKwys9ujJrRxZva9lOMXRM1zE8xsy7wELSIFr0MHGDcuzEa9Zk3S\n0UiqxAYSRE1r5cAV7j49Zf9IYD93PyFNmR6EWtCB7v6cmd0ArHX3UWa2HbDK3d3M/gC0c/ez01zD\nR40a9d/3ZWVllJWV5fjbiUghGDEC3nsP/vEPyNED9SWjvLyc8vLy/76//PLLi3f0WtQX8xAw093H\npOwfCpwLHOruG9KU2wF42t07R+/7AJe6+7HVzusEzHD3rmmuoT4dkRKxYQP8+McwcCD86ldJR1Pc\nirZPJ3Ib8HK1hNMPuAQ4Ll3CAXD3lcB7ZrZ7tKsv8HJUvm3KqQOBl/IRuIgUj802g/vug+uvD6uN\nSvKSGL3WG5gLVAAebSOBG4EWwOro1PnuPszM2gHj3b1/VL4bMAFoDrwFnOnua81sEqF/qBJ4Bzg/\nSlLVP181HZES8/jjcOqpsGBB6O+RhtPDoRlS0hEpTVdeCQ8+CHPmQIsWSUdTfJR0MqSkI1KaKitD\n307HjnDTTUlHU3yKvU9HRCRWTZrAxInwyCPhOR5Jhmo6IlJSXnwR+vYNAwu6dEk6muKhmo6ISAa6\ndg2j2QYOhLVrk46m9KimIyIlafhweOMNmDEDmjVLOprCp5qOiEgWrr8+/BwxIiwCJ/FQ0hGRktSs\nGdxzDzz5JNxwQ9LRlA5VKkWkZLVuDQ89BAcdBJ07w4ABSUfU+KlPR0RK3oIFYYnrRx6BHj2SjqYw\nqU9HRCRHevaEsWNDTee995KOpnFT85qICGEI9ZtvQv/+MG8efP/7SUdUONaty921VNMREYn88pdw\nwAHw05/C118nHU1h+PprOOmk3F1PSUdEJGIGf/0rbNoEQ4eGn6XMPQwpz+W/g5KOiEiK5s3DSqMf\nfABnnVXaieeGG0JT43335e6aGr0mIpLGF1/AMcfAzjvDhAlhwtBSMm0aXHABPPUU/OAHWtogY0o6\nIlJfn38ORx0Fe+wRRreVSuJ57rnwvVOHkGvItIhInrVsCf/8JyxdCsOGhTV5Grtly8LQ8QkT8vPM\nkpKOiEgtWrWCmTPDkgjDhzfuedrWrg1Dxi+5JH+zM6h5TUSkHj79FI44Anr1Ch3sja2pbePGkHB2\n2y2srGrVGtLUp5MhJR0RydSaNeEPc9u2YRXSli2Tjig3vvgCTj45NB8+8ED6pR7UpyMiErM2beDx\nx0OyOfjgxjFlzscfw+GHhxkY7r8//2sLKemIiDTAZpvBHXfAkCHwox/B/PlJR5S5996DPn3gwANh\n0iRo0SL/n6mkIyLSQGZhypyxY+G44+DOO5OOqOGWLIHeveHss+HPf46vj0p9OiIiWViyBI49FgYP\nhv/93+IYYDBvHpxwAlx3HZxySv3KaCBBhpR0RCTXVq2CE08Mw6vHjYP27ZOOqGbTpsF558HkyaEv\np740kEBEpEBsuy08+ijsuy906wY331x4D5J+/nloEhw2LDx31JCEk0tKOiIiOdCiBVxxBZSXw5Qp\nob+koiLpqIJHHoF99oGVK2Hx4mRXR1XzmohIjlVWwq23wsiRYabq3/0Ottgi/jhWroSLLgoj7P72\nt/Bwa6bUvCYiUqCaNIFzzw1T5yxbBl26hL6UuJZJcIfbbguf27EjvPRSdgknl1TTERHJs5kzYdSo\nUPM499xQ+8nHYIMvv4QZM8I0NuvXw/jx0L17bq6tmo6ISJE46ih49tlQ21m+PPSvDBwI//pX9gMO\nKithzhw455yQyMaODYlt/vzcJZxcUk1HRCRm69aFwQZjx4b53I4/Hvba65ttm21qLltZCatXw7vv\nhhVOJ08OU9icdlqYJaFDh/zErOd0MqSkIyKFwj0smDZ7dlizp2rbfPOQfPbcM8x+8MEH32wrV4Yk\n07596Kc57bQwTLv6rNC5pqSTISUdESlk7vD++/DKKyEBNWkC7dp9s7VtG+Z/i5uSToaUdEREGq5o\nBxKYWQfiHDHRAAAGLElEQVQzm21mS8yswsyGR/uvMbOlZrbIzO43s9Y1lN/SzO6Lzl1iZgdE+7cy\ns0fN7FUz+5eZbRnn95J4lJeXJx2CZEj3TiCZ0WtfAxe7+97AgcAFZrYn8Ciwt7t3B14Hfl1D+THA\nP919L6AbsDTafxkwy933AGbXUl6KmP5wFS/dO4EEko67r3D3RdHrzwhJY0d3n+XuVYMH5wPfGYMR\n1X4Odvfbo/Jfu/un0eEBwMTo9UTg+Dx+jazl8hcw02s1pFx9zq3rnJqON3R/0nIdVzHcv4YeK9R7\nB8V3/xrb716iz+mY2U5Ad+CZaofOAmamKbIzsMrMbjez581snJl9Lzq2vbuvhJDYgO3zE3VuKOlk\nvj9pxfZHq77nKunEe71STTqJDSQws1ZAOXCFu09P2T8S2M/dT0hTpgehFnSguz9nZjcAa919lJl9\n4u5bpZy72t2/M9rdzDSKQEQkA7kYSJDn1bDTM7NmwFTgzmoJZyhwNHBoDUWXA++5+3PR+6nApdHr\nFWa2g7uvNLO2wIfpLpCLfzQREclMUs1rtwEvu/uYqh1m1g+4BDjO3TekKxQ1n71nZrtHu/oCL0ev\nHwSGRq/PAKYjIiIFJfbmNTPrDcwFKgCPtpHAjUALYHV06nx3H2Zm7YDx7t4/Kt8NmAA0B94CznT3\ntWa2NXAv0BFYBgx29zXxfTMREalLyT0cKiIiydEs0yIiEhslHRERiU1JJx0z29nMJpjZvdH7AdGz\nP3eZ2eFJxyc1S3Pv9jSzW8zsXjP7WdLxSe2q379o3xZmtsDMjk4yNqldmt+9H5vZ3Oj375C6ypd0\n0nH3t939nJT30939PODnwODkIpO6pLl3r7j7z4GTgIOSi0zqo/r9i1wK3JNEPFJ/ae6dA+uAzQiP\ntdSqUSUdM7vVzFaa2YvV9vczs1fM7DUzu7Sm8il+C/w1P1FKOrm4d2Z2LPAQ8M98xirfle39M7PD\nCI8/fAToWboYZXvv3H2uux9DmP/y93V9XqNKOsDtwJGpO8ysCfCXaP/ewMnRBKPfOi3l/KsIE4ou\nynOs8m1Z3zt3nxH9x39qnmOV78r2/pUBBwBDgOo1IMmvrH/3ImsIj73UKpEZCfLF3eeZWadqu3sB\nr7v7MgAzu5swOegr0bM9fwS6R5n8C8IDp63NbFd3Hxdj+CUtB/duPjCQUMV/OL7IBbK/f+7+2+ic\n04FVMYZe8nLwu/caITltSUhUtWpUSacGOwLvpbxfTvgHxd0/JvTfpLopprikbg29d3Niikvqp6H3\nD3efFE9oUoeG3rsH6nvhxta8JiIiBawUks5/gB+kvO8Q7ZPCp3tX3HT/ilfe7l1jTDrGtzu4FgC7\nmlknM2sB/JQwOagUHt274qb7V7xiu3eNKumY2RTgKWB3M3vXzM50903AcMJy2EuAu919aW3Xkfjp\n3hU33b/iFfe904SfIiISm0ZV0xERkcKmpCMiIrFR0hERkdgo6YiISGyUdEREJDZKOiIiEhslHRER\niY2SjkgOmNkO0Yqzr0erXz5kZrvWcG4nM6uIO0aRQlAKs0yLxOEB4HZ3PxnAzLoAOwBv1HC+nsqW\nkqSajkiWzOz/AV+5+/iqfe5e4e7/NrM/mVmFmS02s+8sgW5mZ5jZTSnvZ1StM29m68zsGjN7ycwe\nNbOeZvaEmb1hZv1Tyt9vZjPN7FUzuzqGryySMSUdkeztAyysvtPMBgJd3b0LcDjwJzPbIU35mmo9\nLYFZ7r4P8BlwBWGRwYHR6yrdgEFAV+AkM9sx0y8ikm9KOiL50we4C8DdPwTKgZ4NKL/B3R+NXlcA\nc9y9MnqdutLj4+7+mbtvAF6udkykoCjpiGRvCbB/Pc6rvqY8wNd8+/dw85TXG1NeVwIbADzM0pva\nH7sh5fUm1FcrBUxJRyRL7j4baGFm51TtiwYSrCE0dzUxs+2Ag4Fnq06Jfr5DWGvezKwj0ZLA1c5J\np7ZjIgVL/0ckkhs/AcaY2WXAekIy+QWhX2YxoaZyibt/aGadiPpxosEG7xBqS0v5dt9QbSPcajqm\nUXFS0LSejoiIxEbNayIiEhslHRERiY2SjoiIxEZJR0REYqOkIyIisVHSERGR2CjpiIhIbP4PuA+G\nH/FVPbgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x179462240>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cols = np.logspace(12,15)\n",
    "trots = [trot_radex(column=c).to(u.K).value for c in cols]\n",
    "pl.semilogx(cols, trots)\n",
    "pl.hlines(tr_swift(20), cols.min(), cols.max(), color='k')\n",
    "pl.xlabel(\"Column\")\n",
    "pl.ylabel(\"$T_{rot} (2-2)/(1-1)$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x1794987b8>"
      ]
     },
     "execution_count": 121,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAETCAYAAAD+spv+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYHFW9//H3JwuQhU1AggnrZZPtFxEiGITRoEZkEVQW\nuaJyAeHKZkBBQBMWHyIIyCqXACJIQEHCEgEJwhDhagghgQATlssaSIhsISErM9/fH1VjOpPuWau7\npns+r+fpp6ura/meTKa+c86pOkcRgZmZWRZ65R2AmZnVDicVMzPLjJOKmZllxknFzMwy46RiZmaZ\ncVIxM7PMVDypSBoi6SFJz0qaKenEgu9OkNSQrh9bYv9XJT0labqkxysXuZmZtaVPDuf8GBgVETMk\nDQSmSXoAGATsB+wYER9LWr/E/k1AXUS8X6F4zcysnSqeVCJiLjA3XV4oqQEYDBwDjI2Ij9Pv3ilx\nCOFmOzOzbinXi7OkzYChwBRga2BPSf+U9LCkXUrsFsAkSVMlHV2ZSM3MrD3yaP4CIG36uh04Ka2x\n9AHWjYjdJO0K/AnYosiuwyNijqQNSJJLQ0Q8WuT4Hn/GzKwTIkKd3TeXmkqaQG4HboqIu9LVbwB3\nAETEVKBJ0not942IOen7v4AJwLBS54mILr1Gjx7dpW1KfddyfWuf27NcrrLlVb6ulq2c5WvPulou\nX3vKWk3ly+N3r73HyOva0lV5NX9dDzwXEZcWrLsT+BKApK2BvhHxbuFOkvqnNRwkDQC+AjxTriDr\n6uq6tE2p71qub+1zqeWuau+xXL6Or6vl8rW3rF1VqfLl8bNr7/Gq9Xevyxm3oy9gONAIzACmA08C\nI4G+wE3ATOAJYK90+42Aieny5gX7zQROb+U8UctGjx6ddwhlU8tli3D5ql2tly+9dnb6Gp/H3V+P\nAb1LfP3dItvPAfZNl18h6djv8bL+y6k7qeWygctX7Wq9fF2lyKANrTuSFLVaNjOzcpFEVFtHvZmZ\n1SYnFTMzy4yTipmZZcZJxczMMuOkYmZmmXFSMTOzzDipmJlZZpxUzMwsM04qZmaWGScVMzPLjJOK\nmZllxknFzMwA+Oijrh/DScXMzGhogO226/pxnFTMzHq4p5+GESPg3HO7fiwnFTOzHuzJJ+ErX4FL\nLoEjjuj68So+SZeZmXUPU6bA/vvD1VfDgQdmc0wnFTOzHujRR+Ggg+B3v4Ovfz274zqpmJn1MA8/\nDAcfDDffnDR9ZanifSqShkh6SNKzkmZKOrHguxMkNaTrx5bYf6SkWZJekHRa5SI3M6t+kybBIYfA\nbbdln1AghznqJQ0CBkXEDEkDgWnAAcAg4Axgn4j4WNL6EfFOi317AS8AI4C3gKnAoRExq8h5PEe9\nmVmBp56CvfeGCRNgjz2Kb1N1c9RHxNyImJEuLwQagMHAccDYiPg4/e6dIrsPA16MiNciYjlwK0lC\nMjOzVsyZk3TKX3ll6YSShVxvKZa0GTAUmAJsDewp6Z+SHpa0S5FdBgNvFHyena4zM7MSFi+Gb3wD\njjoq6Uspp9w66tOmr9uBkyJioaQ+wLoRsZukXYE/AVt05Rxjxoz593JdXR11dXVdOZyZWdVpaoLv\nfx+23BLOOmvV7+vr66mvr8/sfBXvUwFIE8hE4L6IuDRddy/wq4h4JP38EvC5iHi3YL/dgDERMTL9\nfDoQEfGrIudwn4qZ9XijRyed8w89BGus0fb2VdenkroeeK45oaTuBL4EIGlroG9hQklNBbaUtKmk\n1YBDgbsrEbCZWbUZPx5+//ukY749CSULedxSPBw4HPiSpOmSnpQ0EvgdsIWkmcB44Ih0+40kTQSI\niEbgeOAB4Fng1ohoqHQZzMy6u3/8A04+Ge65BzbcsHLnzaX5qxLc/GVmPdVrr8Huu8O4cR1/Wr6r\nzV9OKmZmNWT58uSW4W9/G049teP7V2ufipmZlcEvfwnrrAOjRuVzfo/9ZWZWI/7xj2TE4SefhF45\nVRlcUzEzqwELFsB3vwtXXQWf+lR+cbhPxcysBvzXfyXv113XteN0tU/FzV9mZlXujjugvh5mzMg7\nEtdUzMyq2ltvwWc+A3femdxG3FW++8vMrIdqaoIf/ACOOy6bhJIFJxUzsyp1xRUwfz6ceWbekazg\n5i8zsyr03HOw557JbcRbbZXdcf1EfQlOKmZWqyKgri55av7447M9tvtUzMx6mJtvhoULk76U7sY1\nFTOzKvLBB7Dddslw9p/7XPbHd/NXCU4qZlaLTjwRliyBa64pz/H98KOZWQ8xfTr88Y/w7LN5R1Ka\n+1TMzKpAUxP86Edw3nmw/vp5R1Oak4qZWRW44YYksTSP8dVduU/FzKybe++9pHP+3nth553Ley53\n1JfgpGJmteLYY6FPn+QJ+nKruo56SUOAG4ENgSbgmoi4XNJo4GhgXrrpGRFxf5H9XwXmp/suj4hh\nFQnczCwHU6fCXXclT9BXg4rXVCQNAgZFxAxJA4FpwAHAIcCCiLi4jf1fBj4bEe+3sZ1rKmZW1Rob\nYbfdkqfmv/e9ypyz6moqETEXmJsuL5TUAAxOv25PQYRvMDCzHuAPf4C+feGII/KOpP1yvThL2gwY\nCkxJVx0vaYakayWtXWK3ACZJmirp6AqEaWZWcUuWwC9+ARdeCOp0vaHycnv4MW36uh04Ka2xXAWc\nExEh6TzgYqDYzXPDI2KOpA1IkktDRDxa7Bxjxoz593JdXR11dXVZF8PMrCyuugqGDoXhw8t7nvr6\neurr6zM7Xi53f0nqA0wE7ouIS4t8vylwT0Ts1MZxRlOiH8Z9KmZWrebPh623hocegu23r+y5q3WU\n4uuB5woTStqB3+wg4JmWO0nqn9ZwkDQA+Eqx7czMqtmFF8I++1Q+oWQhj7u/hgOTgZkk/SMBnAF8\nh6R/pQl4FfhhRLwtaSNgXETsK2lzYEK6Tx/g5ogYW+I8rqmYWdWZMwd22CEZ52uTTSp/fj/8WIKT\niplVo+OOg/794aKL8jm/k0oJTipmVm1efBF23x2efx7WWy+fGKq1T8XMzFo46ywYNSq/hJIF11TM\nzLqBJ56A/fdPaisDBuQXh2sqZmY14PTTk4cd80woWXBSMTPL2aRJ8Prr3X+ulPZwUjEzy1FEUkv5\n5S+Tcb6qnZOKmVmO7r47mdHxW9/KO5JsOKmYmeUkAs45J+lLqaZBI1vjpGJmlpN774Xly+GAA/KO\nJDtOKmZmOWiupfz859Crhq7ENVQUM7PqMWkSLFwI3/xm3pFky0nFzKzCIuDss5Mn6GuplgJOKmZm\nFffww/DOO3DwwXlHkj0nFTOzCjvnnKSW0rt33pFkz0nFzKyCHnkEZs+Gww7LO5LyaNcc9en0v98G\ndk9XDQAagUXA08D4iFhSlgjNzGrIuefCmWdCn3ZdfatPm6MUS9oV+AIwKSJmFvn+P4CvA09FxCNl\nibITPEqxmXU3//u/cPjh8MIL3XdIlrJP0iVpx2LJpMh2WwCzI2JZZ4PJkpOKmXU3X/saHHggHHNM\n3pGUluvMj5L6R8SiTh+gjJxUzKw7efzxZHyvF1+E1VfPO5rS8p5PpcP5VtIQSQ9JelbSTEknpOtH\nS5ot6cn0NbLE/iMlzZL0gqTTuhi/mVlFnHtuMhpxd04oWWhP89fFwJ7Ahy2/AraNiI06dEJpEDAo\nImZIGghMAw4ADgEWRMTFrezbC3gBGAG8BUwFDo2IWUW2dU3FzLqFp5+GkSPh5ZdhjTXyjqZ1Xa2p\ntOf+g1OAkyPikiInP7mjJ4yIucDcdHmhpAZgcPMh29h9GPBiRLyWnv9WkoS0SlIxM+suLrgATj65\n+yeULLTZ/JX+uX9Dia/HdeXkkjYDhgJT0lXHS5oh6VpJaxfZZTDwRsHn2axISGZm3c4rr8D998Ox\nx+YdSWW0607piHi/xPqPOnvitOnrduCktMZyFXBORISk84CLgS5NrjlmzJh/L9fV1VFXV9eVw5mZ\nddhFFyV3e621Vt6RFFdfX099fX1mx+vS3V+dPmnyMOVE4L6IuLTI95sC90TETi3W7waMiYiR6efT\nSSpTvypyDPepmFmu5s2DbbeFhgbYcMO8o2mfXO7+ktTVwZqvB54rTChpB36zg4Bniuw3FdhS0qaS\nVgMOBe7uYixmZmVx2WVw6KHVk1Cy0NmBAgZ29oSShgOHAzMlTQcCOAP4jqShQBPwKvDDdPuNgHER\nsW9ENEo6HniAJCFeFxENnY3FzKxcPvwQrr46eT6lJ+lU85ekIyLixjLEkxk3f5lZnn79a5g2DW65\nJe9IOqYStxSbmVkHLF0Kl1wCf/lL3pFUnoe+NzPL2E03wU47wdCheUdSea6pmJllqLExedhxXJee\n4qtena2p3J5pFGZmNWLCBFhvPdhzz7wjyUenkkp3HZnYzCxPETB2bDJwpDrd1V3d3KdiZpaRv/0N\nFi2C/fbLO5L8dCmpSOqfVSBmZtVu7Fj46U+hVw/+c73i86mYmdWiadOSaYK/8528I8lXxedTqRQ/\n/GhmlXTwwbD77vDjH+cdSddUYo560cp8KhHxm86evJycVMysUl56CXbbLRnmfs01846mayoyR72k\ndYsNfy9pQFeGvy8nJxUzq5TjjktuIz7vvLwj6bqyJhVJqwMDI+LddgSycUS80dZ2leKkYmaV8Pbb\n8OlPw6xZ8MlP5h1N15V16PuIWArsLukwSf1KBLCOpGOATTsbhJlZtWoe3r4WEkoW2tv8NQg4Evgk\nsAbJ8C6NwCKSKX2vjYj5ZYyzw1xTMbNyW7AANt88Gd5+iy3yjiYbFelTqUZOKmZWbhddBFOnwq23\n5h1Jdsrdp7IN0BgRL3X2BHlxUjGzclq2LKmd3H037Lxz3tFkp9zzqfwfUCfpqyTNXVMjYlpnT2Zm\nVivGj0866GspoWShQ81fkoYBnyV58PF5oD4iGssUW5e4pmJm5dLUBDvskHTS77133tFkq6IzP0bE\n48Dj6Ym3AY6S1Bd4E/irRy82s55g4kTo1w9GjMg7ku4nk456SRsBX4iIP7Vj2yHAjcCGQBMwLiIu\nK/j+FOBCYP2IeK/I/q8C89N9l0fEsBLncU3FzMpi+HA46aRkaJZa0y3mqI+IOUCbCSX1MTAqImZI\nGghMk/RARMxKE86Xgdda2b8JqCv2hL+ZWbk9+ijMnQsHHZR3JN1Tu5KKpD7At4Hd01UDWPGcytPA\n+IhY0p5jRcRcYG66vFBSAzAYmAVcAvwEuLu1cPA8MGaWk7Fj4dRToY8nYy+qzX8WSbuSjFI8KSJu\nKfL9fwDHSHoqIh7pyMklbQYMBaZI2h94IyJmqvUp0wKYJKkRuCYieuhM0GZWaTNmwPTpcLsnVC+p\nPbl2SURcVOrLiPg/4DJJW0haLSKWtefEadPX7cBJJLWeM0iavv69SYldh0fEHEkbkCSXhoh4tNiG\nY8aM+fdyXV0ddXV17QnNzKyo88+HUaNgjTXyjiQ79fX11NfXZ3a89g7TMgwYAkyJiDfTdSOAuRHx\nbIdPmjSnTQTui4hLJe0APEjSnKb0XG8CwyJiXivHGQ0siIiLi3znjnozy8wLL8Aee8DLL8PAgXlH\nUz6VmE/lXGBb4GVgJ+ChiLgwTQxzI2L9Dp9UuhF4JyJGlfj+FWDnlp3x6fTFvdK+mAHAA8DZEfFA\nkWM4qZhZZo48EjbbDH7xi7wjKa9K3P31QUR8u+CEe0g6HbiApH+jQyQNBw4HZkqanh7jjIi4v2Cz\nIG3+Sm9XHhcR+5LchjxBUqSx31wsoZiZZen11+HOO5PJuKx17amp/AgYDxwC/D4iFktaBzgUOD8i\n1i1/mB3nmoqZZeXEE5N+lAsuyDuS8qtETeUa4FvAIJJnRIiID4Cr0zuwzMxq1ttvwx/+AM89l3ck\n1cFD35uZteJnP4MPP4Qrr8w7ksrwfColOKmYWVe9/z5suSU8+SRs2kPmti3rdMKtnPSbnT2hmVm1\nuPJK2G+/npNQstDZgQZq+C5tMzP46KNkaPvJk/OOpLp0dgwttyuZWU275hrYay/Ydtu8I6kuHhLN\nzKyFpUuT+efvuSfvSKqPR/s1M2vhd7+DnXaCz3wm70iqj2sqZmYFliyBX/4S/vznvCOpTp2tqXjg\nZzOrSePGJTWUYUXnlLW2+DkVM7PUokXJcyl/+UvPbfrK5TkVM7Na9Nvfwu6799yEkoVc5lOpBNdU\nzKwjFi5MaikPPgg77JB3NPmpyvlUKsFJxcw64vzz4emn4ZZVJk3vWapuPhUzs+5m/ny45BL4+9/z\njqT6tadPZYmkdSUdK6lfOh/81cAx+JZkM6sBv/kNfO1rsM02eUdS/Tyfipn1aO+9B5dfDlOm5B1J\nbWi1T0XS6sDAiHi3zQNJG0fEG1kG1xXuUzGz9jjrLJg7F669Nu9IuodKdNTvC6wJ3BkRi4t8vw5w\nMPBc2jTWLTipmFlb3nknafLqSfOltKUik3RJGgQcCXwSWIOk2awRWATMBq6NiPntDHgIcCOwIUlz\n2riIuKzg+1OAC4H1I+K9IvuPBH5D0h90XUT8qsR5nFTMrFU//WlyK/FVV+UdSfdRdTM/pglqUETM\nkDQQmAYcEBGz0oRzLbAN8NmWSUVSL+AFYATwFjAVODQiZhU5j5OKmZU0Z07yPMrTT8PgwXlH033k\n9kS9pLUk9ZbUtyP7RcTciJiRLi8EGoDmH+klwE9a2X0Y8GJEvBYRy4FbgQM6Hr2Z9XQ//zkceaQT\nStY6fEtw+nT9PiTPqNxA8qT9Y505uaTNgKHAFEn7A29ExEypZJIcDBTeDDCbJNGYmbXbjBkwcSI8\n/3zekdSeDieViHgceFzS14A9gdXpRFJJm75uB04i6Z85A/hy4SYdPWZLY8aM+fdyXV0ddXV1XT2k\nmVW5CBg1CkaPhrXXzjua/NXX11NfX5/Z8TrdpyJpPWDdiHipE/v2ASYC90XEpZJ2AB4k6fgXSe3n\nTWBYRMwr2G83YExEjEw/nw5Esc5696mYWTF33QVnnpnUVvr48e1V5NZRL+k4oDfwfkTc3MF9bwTe\niYhRJb5/Bdg5It5vsb438DxJR/0c4HHgsIhoKHIMJxUzW8myZbD99nDFFfDVr+YdTfeU59D3/4qI\nK4D/68hOkoYDhwNfkjRd0pPpbcKFgrT5S9JGkiYCREQjcDzwAPAscGuxhGJmVsyVV8JWWzmhlFNX\naiqbAKcA0yPihiyDyoJrKmZW6N13YdttYfJk+PSn846m+6rEE/W9IqKpyPqTI+I3nT1xuTmpmFmh\nE05IOumvuCLvSLq3Sgx9/7ikXwGzImJmwfqNJX0PeLB54i4zs+6ooQFuvTV5t/JqT5/KuIi4rUVC\nAbgLeJjktmIzs27r1FPhZz+D9bvllIK1pT1JZVdJW7dcGRGTgcUR0cPnSTOz7uyBB5KHHH/0o7wj\n6Rna06cyieThxK2Bl0hu452Svh8UEb8td5Cd4T4VM1u+HHbeGc45Bw48MO9oqkMlOuqPiYhr0uUt\nSYZFGQbsCuwUEWt29uTl5KRiZuefD/X1cP/9UHr0JytUiaRyH7B/OoBjy+9OjYhfd/bk5eSkYtaz\nPf88DB8OTzwBm22WdzTVoxJJZSvgM8ALzaMLF3z3/yLiqc6evJycVMx6rqYm2GsvOPjg5FZia7+q\nm0+lUpxUzHquq66CP/wB/v536N0772iqi5NKCU4qZj3T668nnfOTJ8N22+UdTfXJc+wvM7NuJQKO\nPRZOPtkJJS9OKmZWM8aPhzffTOaet3y4+cvMasK8ebDjjvCXv8Auu+QdTfVyn0oJTipmPcthh8GQ\nIXDhhXlHUt0qMaCkmVm3ds89MHUqXHdd3pGYk4qZVbV58+C44+Cmm6B//7yjMXfUm1nVamyE73wH\nvvc9+OIX847GwEnFzKrYOeckieXss/OOxJq5+cvMqtIDD8C118K0adDHV7Juo+I/CklDgBuBDYEm\n4JqIuFzSOcAB6bq3ge9HxNwi+78KzE+3Wx4RwyoVu5l1D7NnJ01et9wCgwblHY0VqvgtxZIGAYMi\nYoakgcA0kmQyOyIWptucAGwXEccV2f9l4LMR8X4b5/EtxWY1aPlyqKuDffdNZnO0bFXdMC0RMbd5\ntOM0iTQAg5sTSmoASU2kGOG+ILMe62c/g3XWgdNOyzsSKybXlkhJmwFDSWaSRNJ5wBHAB0CpezkC\nmCSpkaTpbFz5IzWz7uCuu+D225N+lF7+07Jbyi2ppE1ftwMnNddSIuIs4CxJpwEnAGOK7Do8IuZI\n2oAkuTRExKPFzjFmzIrd6+rqqKury7QMZlY5L78MRx+dPOi43np5R1M76uvrqa+vz+x4uQzTIqkP\nMBG4LyIuLfL9xsC9EbFjG8cZDSyIiIuLfOc+FbMa8f77yaRbRx/tSbfKrer6VFLXA88VJhRJWxZ8\n/w2SvpaVSOqf1nCQNAD4CvBMmWM1sxwtWgT77Qd77w3HH593NNaWPO7+Gg5MBmaS9I8EcAZwFLAN\n0Ai8BhybNnNtBIyLiH0lbQ5MSPfpA9wcEWNLnMc1FbMqt3w5HHggfOITcMMN7kepBI9SXIKTill1\na2qCH/wA3n0XJkyAvn3zjqhn8CjFZlZzIuAnP4GXXoJJk5xQqomTipl1OxdckAzDMnmyRx6uNk4q\nZtatXHcdXH01PPYYrLtu3tFYRzmpmFm3cdtt8POfwyOPwKc+lXc01hm+l8LMuoUrroCTT4Z774Wt\ntso7Guss11TMLFdNTcl4XnfeCY8+CptvnndE1hVOKmaWm2XL4MgjkyFYHnsM1l8/74isq9z8ZWa5\n+PBD2Gcf+OgjePBBJ5Ra4aRiZhX31luw556wzTbJqMO+bbh2OKmYWUU99RR8/vNwyCFJ53zv3nlH\nZFlyUjGzimhqgksuSQaGPP/8pHNenR4MxLord9SbWdm99RZ8//uwcCFMmQJbbJF3RFYurqmYWVnd\neSfsvDMMH54Mu+KEUttcUzGzsvjoIxg1KhkQ8o47kn4Uq32uqZhZ5v72N/jsZ2HxYpgxwwmlJ3FN\nxcwy8/zzyZD1zzwDv/41HHRQ3hFZpbmmYmZd9s47ydzxe+yRzCXf0OCE0lM5qZhZpy1dChddBJ/+\ndPK5oQFOOQVWXz3fuCw/FU8qkoZIekjSs5JmSjohXX+OpKckTZd0v6RBJfYfKWmWpBcknVbZ6M0M\nkluDr7wySSYPP5zc1XX55R5qxXKYoz5NFoMiYoakgcA04ABgdkQsTLc5AdguIo5rsW8v4AVgBPAW\nMBU4NCJmFTmP56g3y9jrrydPwV9/fdLM9eMfJ01eVju6Okd9xWsqETE3ImakywuBBmBwc0JJDQCa\niuw+DHgxIl6LiOXArSQJyczK6J//TIZVGToUli+HqVPhz392QrFV5Xr3l6TNgKHAlPTzecARwAfA\nF4vsMhh4o+DzbJJEY2YZe/PNZCbG8eOTjvgTT4Rx42CttfKOzLqz3Drq06av24GTmmspEXFWRGwC\n3AyckFdsZj3VvHlw1VVJ09aOOybPmJx9Nrz4YjIroxOKtSWXmoqkPiQJ5aaIuKvIJuOBe4ExLda/\nCWxS8HlIuq6oMWNW7F5XV0ddXV2n4jWrVRHw0kvJfCZ33JE0a33968kdXF/9qu/i6gnq6+upr6/P\n7HgV76gHkHQj8E5EjCpYt2VEvJQunwB8ISIObrFfb+B5ko76OcDjwGER0VDkHO6oNyti7tzkiffm\nV2MjjBgB+++fTJrVr1/eEVqeutpRn8fdX8OBycBMINLXGcBRwDZAI/AacGxEzJG0ETAuIvZN9x8J\nXErSdHddRIwtcR4nFevxli+HZ5+FadPgiSfg739P+kq++MUkkYwYkUyU5SHorVnVJZVKcVKxnubD\nD5NhUp55Jkkg06bBzJmw6aawyy7JWFyf/3wyYrAnxrJSnFRKcFKxWrRoEbzxBrzySpJAZs1a8T5/\nPmy9NWy/fZJAdtkluQV4zTXzjtqqiZNKCU4qVk0aG5Pbdt9+e8Vr7twkgbz++orXggWw8cZJ7WPb\nbZOmq+b3IUOglwdesi5yUinBScUqKQKWLEku+h9+uOr7Bx/A++/De++t+j5vXvK+zjqw4YYrvzbZ\nZOXXBhs4cVh5OamU4KRSmyKSuc4//jjphG5+L3y1XLdsWfHPy5YlAyIuXbpiufl9yZLir8WLk8mn\nFi1KXoXLvXsnz3Gsueaq7+usA5/4BKy77qrvn/xkkiz6eCIK6wacVEpwUumcpqaVL5yLFyev5gtq\n4fuSJStfiAsv0IWv5gt44efWkkFb7716JRfg5lffvqu+mtevttqKdc3Lze+rr54sr7baiuXm9379\nYI01ir8GDID+/Ve8N7/69s37p2fWdU4qJdRiUolILujLlydt8C1fy5YlTS3z5694NX9ubopZuHDl\n9wULkgTSnESWLFn5wtm/f/K5+SLb8mLbfCEufBVenJsv4sUu7qUSQbH3wmU3/5iVj5NKCZLihz8M\nBg5klVe/fiu/+vdfsVzsYtee2y8jVv3Le/ny5C/3wr/4Fy9e8fmjj5KLe+GFvvBi35wImpcXLlxx\nQe7de9VX376w9trJa621ViyvvXbSBFP4GjhwxfvAgSuSSL9+vmib9WRdTSo13Yo7dOiKC/XcuSuW\nW17kCy/0xRIDJH8hl3pArLmmUKwpZvXVV1ysC5NXv34rJ7p11006Ypsv8MXa5gcOdLu7mXVvNV1T\nyapsjY1Je34pzW38firZzKqdayoV0Ny8ZGZmrXPruZmZZcZJxczMMuOkYmZmmXFSMTOzzNR0R718\nO5aZWUXVdFKp1dulzczKpat/jLv5y8zMMuOkYmZmmXFSMTOzzFQ8qUgaIukhSc9KminphHT9BZIa\nJM2Q9GdJa5XY/1VJT0maLunxykbffdTX1+cdQtnUctnA5at2tV6+rsqjpvIxMCoitgd2B46XtC3w\nALB9RAwFXgR+VmL/JqAuIj4TEcMqEnE3VMv/sWu5bODyVbtaL19XVTypRMTciJiRLi8EGoDBEfFg\nRDSlm/0TGFLiEKJCcbfnP09r25T6ruX61j6XWu6q9h7L5ev4ulouX3vL2lWVKl8eP7v2Hq9af/dy\n7VORtBnPX+IFAAAHP0lEQVQwFJjS4qsjgftK7BbAJElTJR1dvuhq+wfvpNL2dk4qq65zUslGLV9b\nchv6XtJAoB44NyLuKlh/JrBzRHyzxH4bRcQcSRsAk4DjI+LRItv5IRUzs06ouqHvJfUBbgduapFQ\nvg/sA3yp1L4RMSd9/5ekCcAwYJWk0pV/FDMz65y8mr+uB56LiEubV0gaCfwE2D8ilhbbSVL/tIaD\npAHAV4BnKhCvmZm1Q8WbvyQNByYDM0n6RwI4E7gMWA14N930nxHx35I2AsZFxL6SNgcmpPv0AW6O\niLEVLYCZmZVUs9MJm5lZ5fmJejMzy4yTipmZZabHJBVJ20r6raQ/STo273jKIb2RYaqkffKOJWuS\n9pI0Of0Z7pl3PFlT4jxJl0n6bt7xZEnSHunPbZykVe7UrHaSNpY0QdK1kk7LO56sSfq0pD9KulJS\n0Uc9CvWYpBIRsyLiOOAQ4PN5x1MmpwF/zDuIMglgAbA6MDvnWMrhAJJRJJZRY+WLiEfT372JwO/z\njqcMdgRui4ijSB7mrjVfAy6LiB8BR7S1cdUmFUnXSXpb0tMt1o+UNEvSCy3/apC0H8l/7HsrGWtn\ndLR8kvYGngP+RTKUTbfW0fJFxOSI+DpwOnBOpePtqE78/9wGeCwiTgX+u6LBdlBnfvdS3wHGVybK\nzutE+f4JHCXpQeD+igbbCZ0o303AoZIuAD7R5gkioipfwB4kfxU8XbCuF/ASsCnQF5gBbFtk34l5\nx591+YDzgIuBvwIT8o6/XD8/ktvO/5R3/GX4+R0OfCtdvjXv+LP+2QEbA/+Td+xl+tmdAuyRLt+W\nd/zl+PkVbNPmtaVqpxOOiEclbdpi9TDgxYh4DUDSrSTNCrMk7QUcRNJ88peKBtsJHS1fRJyVrjsC\neKeiwXZCJ35+BwJfBdYGrqhosJ3Q0fIBdwCXS/oC8EhFg+2gTpQN4L+A31Uuys7rRPnuB8ZIOhx4\npaLBdkInfvc2Bc4A+gMXtnX8qk0qJQwG3ij4PJvkH4uIeIRu/svaDiXL1ywibqxoRNlq7ec3geTB\n12rWWvkWA0flEVRGWv2/GRFjKh1Qxlr72T0LfDuPoDLUWvleA37Y3gNVbZ+KmZl1P7WWVN4ENin4\nPCRdVytcvupWy+Wr5bKBy9du1Z5UxMp3Ok0FtpS0qaTVgEOBu3OJLBsun8vXXdVy2cDl63z58r4T\noQt3MIwH3gKWAq8DP0jXfw14nmRK4tPzjtPlc/lqrXy1XDaXr+vl84CSZmaWmWpv/jIzs27EScXM\nzDLjpGJmZplxUjEzs8w4qZiZWWacVMzMLDNOKmZmlhknFatJkh6S9OUW606SdGUr+2wqaWb5o1vp\nnI2SnpT0jKTpkkZJynQ+HEk/lPSf6fL3JA3K8vhmhWptlGKzZuOBw4BJBesOBU5tY79KPw38UUTs\nDCBpfeAWYC1gTFYniIj/Kfj4feAZYG5Wxzcr5JqK1ao/A/tI6gNJLQTYKCIeSz9fKGmmpKckHdxy\n5/Qv+ssLPt8jac90eYGkC9LaxQOSdpX0sKSXJO2bbtMr3WaKpBmSjm4r4Ih4BzgGOL61Y0jaKz3f\nbZIaJN1UEOfYNK4Z6Ux9SBot6RQl84vvAvwhrR3tI2lCwb57S7qjo//QZoVcU7GaFBHvS3qcZDyj\ne0hqKX8CkHQQsFNE7Cjpk8BUScXm2ilVaxkAPBgRP00vwucCI4AdSOZgn0gyKdUHEfG5dIC+xyQ9\nEOkkSK3E/UqaTDYAvlHsGOmmQ4HtSGocj0n6PMmEUd+IiG3Tcq618qHjz5KOB0ZFxPR0m19LWi8i\n3gV+AFzXWnxmbXFNxWrZrSTJhPS9eX70PUiamYiIeUA9sGsHjrs0Ipov7jOBRyKiKV1unlHvK8AR\nkqYDU0jm9t6qg/G3dozHI2JOJIP3zQA2A+YDiyVdm86UubjEcQv7bG4C/lPS2sBuwH0djNFsJa6p\nWC27C7hY0meAfhExo8R2xTrGP2blP7rWKFheXrDcRDLaKxERzc1t6TFPiIjCPp02SdoCaIyIf6Ud\n9qscI50ae2nBqkagT0Q0ShpGUmv6Nkkz2og2TnkDSU1uKcn86k0didesJddUrGZFxEcktZDrSWsm\nqb8DhxQ0M30BeDz9rjnBvAoMVWJjVp62ubW7s5q/+yvw3wV9OltJ6tfK9qSx/BZo7sspdoz+JU+c\nfLdORNwPjAJ2KrLZApIbAQCIiDkkw6CfSZXMIW/dm2sqVutuAe4ADmleERETJO0GPEVS0/hJRMxL\nO/Mj3eYxSa8CzwINwLSCY7Z2h1jzd9eSNEk9mdY45pH0kbS0hqQngdVIakA3RsQlHTxG8znXAu6S\n1Fyr+nGRbW8Arpa0CNg9IpYCNwPrR8TzrZTLrF08n4pZD5fe5fZkRLimYl3mpGLWg0l6AlgIfDki\nlre1vVlbnFTMzCwz7qg3M7PMOKmYmVlmnFTMzCwzTipmZpYZJxUzM8vM/wfi8gbjItZRTwAAAABJ\nRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x179b85c50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "densities = np.logspace(3,9)\n",
    "trots = [trot_radex(density=n).to(u.K).value for n in densities]\n",
    "pl.semilogx(densities, trots)\n",
    "pl.hlines(tr_swift(20), densities.min(), densities.max(), color='k')\n",
    "pl.xlabel(\"Volume Density\")\n",
    "pl.ylabel(\"$T_{rot} (2-2)/(1-1)$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x17a4398d0>"
      ]
     },
     "execution_count": 118,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEPCAYAAABY9lNGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcjvX+x/HXRymVJSlxTqWkkmQrJMmcFO2kRsuxVDqo\nLIcS1a/DyTmVDmXN2HeiZDlOkmKIsiYmImnTNlqQMGOZz++P61bSDDNjZq77nnk/H48eZq77vq/r\n87hi3nN9V3N3REREDlco7AJERCQ6KSBERCRdCggREUmXAkJERNKlgBARkXQpIEREJF2hBISZlTCz\nV8zsIzNbZ2a1zaykmb1pZhvNbK6ZlQijNhERCYT1BNEfeN3dLwaqAhuA7sBb7n4RMB94PKTaREQE\nsLyeKGdmxYHV7n7+Ycc3APXdPdnMygCJ7l4xT4sTEZFfhfEEcR7wg5mNNrP3zWyYmZ0MnOnuyQDu\n/h1QOoTaREQkIoyAOB6oAQx29xrALoLmpcMfZbQGiIhIiI4P4ZpfAVvcfWXk+2kEAZFsZmce0sS0\nNb0Pm5mCQ0QkG9zdsvL+PH+CiDQjbTGzCyOHGgDrgFnAvZFjrYCZRzhHzP7Xo0eP0GtQ/eHXURDr\nj+XaY7X+1NRU4uPjufbaa7P18zqMJwiAjsBEMysMfArcBxwHTDWz+4EvgGYh1SYiEvP27NlDfHw8\nxx13HP/973856aSTsnyOUALC3dcANdN5KXsxJyIiv/rll1+49dZbKVOmDGPHjqVw4cLZOo9mUuex\nuLi4sEs4Jqo/XLFcfyzXDrFT//bt22nYsCHnn38+48ePz3Y4QAjzII6VmXms1Swikhe+//57GjZs\nSP369XnxxRcx+61P2szwaO+kFhGRnPf1119Tv359br755j+EQ3YpIEREYtwnn3xCvXr1aNWqFb16\n9cqRcAAFhIhITFuzZg3169fn8ccfp1u3bjl67rCGuYqIyDFasmQJTZs2ZeDAgTRrlvMzAxQQIiIx\n6I033qBFixZMmDCBRo0a5co11MQkIhJjpkyZQqtWrZg5c2auhQMoIEREYsrQoUPp0qUL8+bN48or\nr8zVa6mJSUQkBrg7zzzzDCNHjmThwoVUqFAh16+pgBARiXJpaWl06tSJRYsWsWTJEsqWLZsn11VA\niIhEsdTUVFq2bElycjKLFi2iRIkSeXZt9UGIiESpnTt3ctNNN7F//37eeOONPA0HUECIiESlrVu3\nEhcXR4UKFZg6dSpFihTJ8xoUECIiUebTTz+lbt263HLLLQwZMoTjjjsulDoUECIiUWT16tXUq1eP\nLl260LNnzxxbVyk71EktIhIl5s6dS4sWLUhISKBp06Zhl6MnCBGRaDBmzBhatWrF9OnToyIcQE8Q\nIiKhOjgBbvjw4SxYsICLL7447JJ+pYAQEQnJ/v37ad++PcuWLeO9997LswlwmaWAEBEJwa5du7j7\n7rtJSUlh4cKFFC9ePOyS/kB9ECIieWzr1q00aNCAkiVLMnv27KgMB1BAiIjkqQ0bNnDFFVfQsGFD\nxowZwwknnBB2SRlSE5OISB5ZuHAhzZo147nnnuO+++4Lu5yjUkCIiOSBiRMn0rlzZyZNmsS1114b\ndjmZooAQEclF7s6///1vhg8fzvz586lcuXLYJWWaAkJEJJfs27ePtm3bsmbNGpYuXRp1w1iPRgEh\nIpILtm3bRnx8PEWKFGHhwoUULVo07JKyLJRRTGb2uZmtMbPVZrY8cqykmb1pZhvNbK6Z5e3C5yIi\nOWTz5s3UqVOHypUrM3PmzJgMBwhvmGsaEOfu1d29VuRYd+Atd78ImA88HlJtIiLZ9s4771C3bl06\ndepEv379QluqOyeEFRCWzrUbA2MjX48FmuRpRSIix2j8+PHcfvvtjB07lgcffDDsco5ZWH0QDswz\nswPAUHcfAZzp7skA7v6dmZUOqTYRkSxJS0vjH//4B5MmTWLBggVccsklYZeUI8IKiLru/q2ZnQG8\naWYbCULjUId/LyISdfbs2cO9997LV199xbJlyzjjjDPCLinHhBIQ7v5t5M/vzWwGUAtINrMz3T3Z\nzMoAWzP6fM+ePX/9Oi4ujri4uNwtWEQkHd988w1NmjThggsu4O233w5l3+iMJCYmkpiYeEznMPe8\n/UXdzE4GCrn7L2Z2CvAm8E+gAfCTu/c2s25ASXfvns7nPa9rFhE53MqVK7ntttt46KGH6N69e6hb\ng2aGmeHuWSoyjIA4D5hO0IR0PDDR3Z8zs9OAqcDZwBdAM3ffns7nFRAiEqopU6bQoUMHhg0bRpMm\nsTGeJiYC4lgpIEQkLGlpafTs2ZNx48Yxc+ZMqlatGnZJmZadgNBMahGRTNi1axctW7YkOTmZ5cuX\nU7p0/h9oqf0gRESO4ssvv+Sqq66iRIkSvP322wUiHEABISJyRIsWLaJ27do0b96ckSNHcuKJJ4Zd\nUp5RE5OISAaGDBlCz549GT9+PA0bNgy7nDyngBAROczevXvp0KEDixcvZsmSJVSoUCHskkKhgBAR\nOURycjK33347p59+OkuXLqVYsWJhlxQa9UGIiESsXLmSmjVr0qBBA1577bUCHQ6gJwgREQDGjRvH\nI488wtChQ2natGnY5UQFBYSIFGj79u2jS5cuvPHGGyxYsCCm9ozObQoIESmwvvvuO+Lj4zn11FNZ\nsWIFp556atglRRX1QYhIgfTee+/92t8wc+ZMhUM69AQhIgWKuzNs2DCeeuopRo4cyS233BJ2SVFL\nASEiBUZKSgrt27fnvffeY/HixVx44YVhlxTV1MQkIgXC559/zlVXXcWOHTtYunSpwiETFBAiku/N\nmTOH2rVr89e//pWpU6cWuPkNe/Zk73NqYhKRfCstLY2nn36a4cOH8+qrr1KvXr2wS8pT69fD0KEw\ncWL2Pq8nCBHJl3788UduuukmFixYwKpVqwpMOKSkwIQJUK8eXHstFC8Oq1Zl71wKCBHJd1auXMll\nl11G5cqVeeuttyhTpkzYJeW6jRvhkUfg7LNh/Hjo3Bm++AJ69YJy5bJ3TjUxiUi+4e4MGTKEHj16\nkJCQwO233x52Sblq716YPj1oRlq3Du67D5Ytg/Llc+b8CggRyRd27txJmzZtWL9+PUuWLMnXo5Q+\n+wyGDYPRo6FSJWjXDpo0gRNOyNnrqIlJRGJeUlISNWvWpGjRovl2COv+/TBzJtxwA9SqBampkJgI\n8+dDs2Y5Hw6gJwgRiXFjxoyha9eu9OnTh1atWoVdTo77+msYORKGDw/6F9q1g9deg5NOyv1rKyBE\nJCbt3r2bDh068O677+a7VVjT0uDttyEhIXhCuOsumD0bqlbN2zoUECISczZs2ECzZs2oXLkyK1as\noGjRomGXlCN+/BHGjAmC4eST4cEHg+/DmtenPggRiSnjxo2jXr16dOjQgYkTJ8Z8OLjD0qXQqhVU\nqABr1sC4cfDBB0FzUpiTvvUEISIxYdeuXbRv356lS5cyf/58Lr300rBLOia//AKTJsGQIbBzZxAG\nffvC6aeHXdlv9AQhIlEvKSmJyy+/HHdnxYoVMR0O69dDhw7B5LU5c6B3b/j4Y3j00egKB1BAiEgU\nc3dGjBjBNddcQ/fu3RkzZkxMNint3QtTp0JcXLD8xWmnBU1J06dDw4ZQKEp/EquJSUSi0o4dO2jX\nrh0ffvghixYt4uKLLw67pCzbsiWY0DZiBFx8MTz8cDChrXDhsCvLnNByy8wKmdn7ZjYr8n1JM3vT\nzDaa2VwzKxFWbSISrmXLllG9enVKlizJ8uXLYyoc0tJg3jy47TaoVg22bw+GrM6fD/HxsRMOAObu\n4VzYrDNwGVDc3W81s97Aj+7+vJl1A0q6e/d0Pudh1SwiuSstLY3nn3+eF154gYSEBJo2bRp2SZm2\nbRuMHRt0OhcpEjwt3HMPREuLmJnh7paVz4TSxGRmZwE3Av8GukQONwbqR74eCyQCfwgIEcmfvv32\nW1q2bElKSgorV67knHPOCbukTFm9GgYPhmnT4MYbYdQouPJKsCz9KI5OYTUxvQh0BQ59FDjT3ZMB\n3P07oHQYhYlI3pszZw41atSgbt26LFiwIOrDITU12HPhyiuhceNg9dSNG4ONeerWzR/hACE8QZjZ\nTUCyu39gZnFHeGuG7Ug9e/b89eu4uDji4o50GhGJVikpKTz++ONMmzaNl19+mfr16x/9QyH64otg\nae2RI4NlL7p1g5tuguOjcLhPYmIiiYmJx3SOTPVBmNnxQDxQJ3LoFOAAsBtYC0xy95RMXdDsGaA5\nsB84CSgGTAcuB+LcPdnMygAL3P0PPVPqgxDJH9atW8c999zDBRdcwLBhwzjttNPCLildaWnw1ltB\nM9LixdCiRbAExkUXhV1Z1mSnD+KoAWFmNYF6wDx3T0rn9fOBm4A17r4wSxc3qw88Eumkfp6gk7q3\nOqlF8i93Z/DgwfTs2ZPevXtz//33Y1HYJrN9e9Dp/NJLv3U6//WvcMopYVeWPbnVSZ3i7i9k9KK7\nbwYGmFl5MzvB3fdmpYBDPAdMNbP7gS+AZtk8j4hEqa1bt3LfffexdetW3n333ajct2Ht2uBpYepU\nuP76oDkpP/UrZMUxDXM1s5PdfXcO1pOZa+oJQiQGzZkzh9atW9OqVSv++c9/ckJu7HCTTfv2BbOa\nBw2CzZuhbVv429+gbNmwK8s5YQxzbQP0O8ZziEg+tnv3brp168bMmTOZNGlSVA0q+e67YKbz0KFw\nwQXBGkmxNNM5tx01IMzsBeBq4OfDXwIqooAQkQysWrWK5s2bU716ddasWUPJkiXDLgl3eO+94Glh\nzhy480544w2I4fX/ck1mniAeAf7u7i8e/oKZ/T3nSxKRWHfgwAGee+45+vfvT//+/bn77rvDLok9\ne2Dy5CAYfvkFHnoo6IA+9dSwK4temR3mWtLdt6Vz/BR335UrlWVci/ogRKLYp59+SosWLShSpAhj\nxozh7LPPDrWezz8PgmD0aKhdG9q3j+4VVHNLdvogMnWL0guHyPE8DQcRiV7uzujRo6lduzZ33HEH\n8+bNCy0c3IO5C40bw+WXw4EDwa5ts2cHI5MKWjhkVxTO/xORWJOcnEybNm344osvWLBgAZUrVw6l\njp07g+06Bw0KZjd36BDs2harcxfClq0cNbPbc7oQEYlN06ZNo2rVqlSuXJnly5eHEg4ffwwdOwa7\ntC1YAAkJwXyGNm0UDsciu08QUbKArYiEZfv27XTo0IFly5YxY8YMrrjiijy9flpaMApp4MBgRdUH\nHgh2aQu5yyNfyW5AqJdYpAB78803eeCBB2jcuDGrV6/mlDz8NX3HjqDDedCgYARShw4wY0awHIbk\nLPVBiEim/fLLLzz22GPMnj2bUaNGce211+bZtdevD0Lh5ZeDjubx4+GKKwrmEhh5RX35IpIpCxcu\npGrVquzevZu1a9fmSTgcOACzZsF110GDBlC6NKxbF3Q816mjcMhteoIQkSPatWvXr3s2DB06lJtv\nvjnXr7l9e7Az26BBcMYZQQf0HXfAiSfm+qXlENl9gng1R6sQkaj0zjvvULVqVbZt20ZSUlKuh8P6\n9cFeC+edB++/H8x8XrYsWGZb4ZD3svUEkdcruIpI3tq9ezdPPvkkU6ZMYciQITRu3DjXrnXgAPzv\nf8FopKQkaNcuCIr8tJJqrFITk4j8zjvvvEPr1q257LLLSEpKolSpUrlynUObkUqVgk6dID5eTwrR\n5JgCIoz9IEQkd/zyyy88/vjjvPbaawwePJgmTZrkynU2bAieFiZNghtuCP6sXVsdztHoWEcxtcmR\nKkQkVPPnz6dKlSrs3LmTpKSkHA+HtDR4/fVgeGr9+sETw8HRSBqqGr0ysyf1EfeDcPc8bSnUaq4i\nOefnn3+ma9euzJkzh4SEBG688cYcPf/OnTBmTPDEULRo0Ix0552a1BaG3NpRTvtBiORDc+bMoV27\ndjRq1IikpCRKlCiRY+fevDnoWxg3Dq65JuhrKKj7OseyowaEu7uZjcng5eE5W46I5LYffviBzp07\ns2TJEkaOHJljE97cYf586N8/2LGtdetgjaRzzsmR00sIjtgHYWYnmlmpzOwHYWZaIkskirk7kyZN\nonLlypxxxhkkJSXlSDjs3g3Dh0OVKkET0s03wxdfwHPPKRxi3RGfINw91cyuM7NiwAx333P4e8zs\nVKAZsB7Ykjtlisix2LJlC+3atePLL79k1qxZ1KpVKwfOGezUNmJEsOxFv35Bc5KakfKPzDQxzTaz\nMkBnMysNFIl87gCwG/gKGOHuO3K1UhHJsrS0NBISEvjHP/5Bp06dmD59OieccEK2z+ceNB/17x/s\n2NaiRfB9hQo5WLREjUztSR1NNIpJJHPWrVtHmzZtcHdGjBhBpUqVsn2uvXvhlVeCp4Rt24K1ke69\nF4oXz7l6JXfl+J7UZnaRmel3A5EYkpKSwlNPPUVcXBzNmzdn8eLF2Q6HrVuhVy8499xgD4YePX7b\nvU3hkP8drYlpMxBnZo0ImpRWuPuq3C9LRLIjMTGRtm3bUrlyZT744AP+/Oc/Z+s8a9cGTwvTpwer\nqM6dC5demsPFStTLUhOTmdUCLiOYJLcRSHT3A7lUW0Y1qIlJ5DA//fQTjz32GHPnzmXgwIHZmgl9\ncNG8fv1g40Z46CFo2xZOPz0XCpY8l1sT5X7l7suB5ZGLXQQ8YGaFga+BuVqXSSRvHRy6+uijj3LH\nHXewbt06imex7efnn4PmowEDgiUw/v734KnhGPqyJZ/I9mJ97r6R4CkCMysL3AxMPdrnzOxEYBFw\nQuT6r7r7P82sJDAFKAd8DjTTyCiRjG3atImHHnqI77//nhkzZlC7du0sff7TT4MlMMaODXZsmzBB\n6yLJ7+XIlqPu/q27HzUcIu9NBf7i7tWBasANkaar7sBb7n4RMB94PCdqE8lvUlNTefrpp6lTpw7X\nX389K1euzHQ4uMPChXDbbVCrFhQuDB98AFOmaAtP+aNMPUGY2fFAPFAncugUfpsHsRaY5O4pmb3o\nIU1RJ0ZqcKAxUD9yfCyQSBAaIhKxYMECHnzwQS666CLef/99zsnkVOXU1CAE+vULZj536hQ8MZxy\nSi4XLDEtM6u51iRYzXWeu69N5/XzgZuANe6+MFMXNSsErALOBwa7++Nmts3dSx7ynp/c/bR0PqtO\nailwtm7dSteuXVmwYAEDBgzIdCf01q0wdCgMGRKMQvr736FRIyiUI20HEktyq5M6xd37ZvSiu28G\nBphZeTM7wd33Hu2E7p4GVDez4sB0M7uE4Cnid2/L6PM9e/b89eu4uDji4uKOdkmRmHTgwAGGDRtG\njx49aNGiBevWraNYsWJH/VxSUvC08NprwS5t8+bBJZfkQcESNRITE0lMTDymc2RqmGukj+AsYJm7\nfx051gD4zt3XHVMBZk8RNFU9AMS5e3JkaY8F7n5xOu/XE4QUCCtXruTBBx/kxBNPZMiQIVx6lIkI\naWkwZw68+GKwp/PDD2uYqvwmO08QmWli6gVUBD4FqgDz3f0/kX6J79w9S3/9zOx0YJ+77zCzk4C5\nwHME/Q8/uXtvM+sGlHT3P/RBKCAkv9u2bRv/93//x7Rp03juuedo2bIlhY7QJrRrVzASqX//YFOe\nzp2hWTMNU5Xfy60mpu3uHn/IRa4ys+7A8xyhGegIygJjI/0QhYAp7v66mS0FpprZ/cAXBCvEihQY\n7s748ePp1q0bTZo0Yf369Zx22h+64X61ZUuwKc/IkXD11cGqqlddpZFIknMy1QcRmaNwJzDW3Reb\n2YcE+1FneR6FuycBNdI5/hOQMzuXiMSYNWvW8PDDD5OSksLMmTOPuBz3smVB/8Kbb0LLlrB8OZQv\nn4fFSoGRmbEMw4DrgTJAGoC7b3f3BOCxXKxNJN/bvn07HTt2pGHDhrRo0YJly5alGw779gWrqV55\nJdx1VzCH4bPPgv4GhYPklszsB7EPmJzBa9pyVCQb0tLSGDt2LE888QS33nor69evp1SpUn94X3Jy\nsFtbQkIQBI8+Co0bw3HHhVC0FDjZXmpDRLLn/fff5+GHHyYtLY1Zs2ZRs2bN373uHjQbDRwYLJ4X\nHx/8WbVqSAVLgZWt6TJmdntOFyKS3/3www+0bduWG264gdatW/Pee+/9LhxSUoLRSLVqwd13Q/Xq\nwXpJw4YpHCQc2Z1PWTRHqxDJx/bv38+gQYOoVKkSRYoUYcOGDTzwwAO/Dl396it48kkoVw4mTQo2\n5dm0CR55BEqWPMrJRXJRdpuYNBFBJBMWLFhAx44dKV26NPPnz6dy5cpA0Iy0eHHQjDRvHjRvDosW\nwUUXhVywyCHUByGSC7788kseffRRli9fTt++fWnatClmRkoKTJ4c7L2waxe0bx/MX9D2nRKNtGSX\nSA7atWsXPXr0oEaNGlSqVIn169dz++238/XXxpNPwjnnBMNVn3kGNmzQ3s4S3RQQIjnA3Zk4cSIV\nK1Zk06ZNrF69mh49evLBBydz111QpUqwc9vixfD663DDDVpRVaJfdpuYXs3RKkRi2IoVK+jUqRN7\n9+7l5Zdf5vLL6zJ1arA20vbt0KFDsOR2iRJhVyqSNdkKCO09LQLffPMNTzzxBG+++SbPPPMMjRq1\nZNiwQtxxB1SuDD17wo036klBYpf+6opk0e7du3n66aepUqUKZcuWZerUjSQm3kulSoX45ht4661g\nZNLNNyscJLZl6q+vmdUys6Zm9udDjjWIbPQjUiCkpaUxYcIEKlasyLp16/n3v1fy7rvPcvfdxahY\nET75JGhK0sY8kl/k+X4Qx0r7QUgYFi9eTJcuXUhLM6688kVmz76SMmWCvZ2bNoXChcOuUOTIYmU/\nCJGY8emnn9K9e3cWL15KlSrPsWzZXVx4YSEmT4batcOuTiR3ZaaJKcXMSppZOzM7yd0XAwlkcz8I\nkVjw008/0aXLI1SvXoukpCrs3buBGjXuISmpEJMmKRykYMjMD/hhwB0cth8EkGBmB3KxNpE8l5qa\nyoABL9Gr17OccEJTSpdeR6dOZ9KyJZx8ctjVieStI/ZBmNmJQFF3//GoJzI729235GRxGVxHfRCS\n49ydUaNepWvX7uzefTE1ajzPU09VolEjjUSS/CHH+yDcPdXMrjOzYsAMd9+TzkVPJdg/ej2Q6wEh\nktMmTFjCI4905ccf99Cw4TCef74BkTX1RAq0o45iAjCzMsD9QGmgCEGwHAB2A18BI9x9Ry7WeWgt\neoKQY+YOo0Zt4IknuvPDD6u59dZeDBnSnDJl9Lgg+VN2niAyFRDRRAEhxyIlBV566Vv+9a9/8vPP\n02ja9DGGDevAqacWCbs0kVyVW8NcM7pYcWAXUCiyb7VI1Pr+e+jXbyf9+/dh795BNG58PwkJGylV\n6rSwSxOJWlkOCDOrBdxIMAdiDHAWsCRnyxLJGR99BH367GXy5KHAM1x77XUMHPg+5cqVC7s0kaiX\n7SYmM7sBOB040d1H5GhVR76umpjkiNxhwQLo0yeNJUsmU6jQU1x2WUX69n2WqtrcWQqoPO2DMLNS\nQEl3/yRbJ8gmBYRkZO9emDIF+vZ1fvzxDeBxypYtwn/+05v69euHXZ5IqPK0D4JgaOtxZrbN3Sce\nw3lEjsmOHTBsWLD/QpkyS0lL607Rosk888wzNGnSBLMs/ZsQkYhjGdP3vbsPAjbnVDEiWfHFF9Cl\nC5x3Hixc+CEXXNCE5OR4OnZsTlJSErfddpvCQeQYHEtALDez/gQrvYrkmVWr4J57oEYN2LnzU+Li\nWrBiRQNuvbU+mzZt4oEHHuD447VMmMixOmpAmFlG72nq7p3cfUxWLmhmZ5nZfDNbZ2ZJZtYxcryk\nmb1pZhvNbK6ZaYNG+VVaGsyeDX/5C9x2G1xwwbc0bfow06fXomrVCmzatInOnTtTpIjmM4jklMw8\nQSw3s3gzu/Sw42ebWatDNxHKpP1AF3e/BKgDPGxmFYHuwFvufhEwH3g8i+eVfCglBYYPDzbh6dED\n7rrrR+68szuDBlWmePEibNiwgR49elC8ePGwSxXJdzITEMPd/RV3Tzrs+ExgAXB1Vi7o7t+5+weR\nr38BPiKYS9EYGBt521igSVbOK/nLDz/A00/DuefCzJnw/PM7uPnmnjz55EXs3LmdNWvW0LdvX04/\nPU/3qxIpUDITEDXN7MLDD7r7ImCPu0/O7sXN7FygGrAUONPdkyPn/o5g3ScpYDZtgoceggsugC+/\nhP/9bxf16vWmdesL+Pzzz1i+fDkJCQmcddZZYZcqku9lpievHDAgEhKfAMuBZZE/mwJDsnNhMysK\nvAp0cvdfzOzwyQ0ZTnbo2bPnr1/HxcURFxeXnRIkirz7LvTpA++8A+3awerVKcyaNYybb36Wq6++\nmoULF3LxxReHXaZIzEhMTCQxMfGYzpGZPanbuPuwyNcVgFqR/2oCVdy9WJYvGuxnPRuY4+79I8c+\nAuLcPTmyeuwCd//DTwRNlMs/DhwImo/69IHk5GDI6t13pzJlyiieeeYZqlWrRq9evahWrVrYpYrE\nvNyaKHebmY12932RWdOfAJMiF3w0G3UCjALWHwyHiFnAvUBvoBVBH4fkQ7t3w9ix8MILUKoUdO0K\nN920jwkTxlKjxr+oVKkS06ZNo1atWmGXKlKgZSYgOhKExMcHO5cPMS+rFzSzusBfgSQzW03QlPQE\nQTBMNbP7gS8IZmpLPvL99zB4MLz0EtSpA6NHQ+3a+5k4cQKVKj1NhQoVmDx5MnXq1Am7VBEhEwHh\n7puATRm8tiarF3T3JcBxGbx8bVbPJ9Fv06bgaWHKFIiPD/oZzj9/P5MnT+b++3vxpz/9iTFjxnD1\n1VkaECciuUzTTSXXvPce/Oc/v3U8f/QRlCoVBMMtt/SibNmyJCQkcM0114RdqoikQwEhOSotDWbN\nCoLh22+Djufx4+HEE/czadIkevUKnhiGDRum0WciUU4BITliz54gCPr2heLFg47npk0Bfh8Mw4cP\nVzCIxAgFhByTH38MOp0HD4bLLw+W3b76ati/fx/jxo3nmWee4ayzzlIwiMQgBYRky2efBR3PEyYE\ni+e9/XbCxGGXAAARMElEQVSwXlJqairDho3h2Wef5YILLmDUqFHqfBaJUQoIyZJVq4L+hbfeggce\ngHXr4E9/gpSUFAYPHknv3r2pXLmyhquK5AMKCDkqd5g3D3r3ho8/hs6dg6ak4sVh165dvPDCUPr2\n7ctll13GtGnTqFmzZtgli0gOUEBIhvbvh1dfheefD/Z7fuwxuOsuOOEE2LFjB//+9yD69+9P/fr1\nmT17NtWrVw+7ZBHJQQoI+YM9e4JZzn37Bs1HTz8NN94IhQrBDz/8QL9+/UhISODGG2/UInoi+ZgC\nQn61bVswImngQKhVC8aNg7p1g9e+/vprXnjhBUaPHk18fDzLly+nfPny4RYsIrnqWPaklnxiy5Zg\nQtv558MnnwQjkmbNCsLhk08+oU2bNlx66aWkpaWxdu1ahg4dqnAQKQAUEAXYunVw771QtSqYwdq1\nQdPSJZfA2rVrueeee7jiiisoU6YMH3/8MS+++KI26hEpQBQQBYw7LF4Mt9wCDRrAhRfC5s1Bf8NZ\nZ8G7777LLbfcQqNGjahWrRqffvopTz/9tLb2FCmA1AdRQKSlwX//GwxV3boVHn0Upk6Fk04Cd+f1\n1+fw3HPPsWXLFrp27crUqVM56aSTwi5bREKkgMjnUlOD2c7/+Q8UKwbdugUzn487Dvbv38/EiVPo\n3bs3Zkb37t2Jj4/n+OP110JEFBD51o4dMHQo9O8Pl14KQ4ZAXFzQ17B7925Gjx5Nnz59KFeuHL17\n9+b666/HLEu7EYpIPqeAyGe++Qb69YORI+H66+H114NOaAjmMAwePJiXXnqJOnXqMGnSJC2HISIZ\nUid1PrF+Pdx/P1SuHMx6fv99mDgxCIfPPvuMDh06cOGFF/LVV1+xcOFCZsyYoXAQkSNSQMSwgyOS\nbr0V/vIXOO+8YHvPfv2gXDlYtWoVd911FzVr1qRYsWKsW7eO4cOHU7FixbBLF5EYoCamGJSWBjNn\nBh3PycnBiKQpU34bkTRnzhv06dOHjz/+mM6dOzN8+HCKFSsWdtkiEmMUEDEkJSXYta1PHyhRIlg8\n7+CIpNTUVEaNmkjfvn0pXLgwXbt2pVmzZhQuXDjsskUkRikgYsC2bcEopIEDoUaNYHRS/frBiKSf\nfvqJhIQEBg0aRNWqVRkwYADXXHONRiSJyDFTH0QU+/LLYO+F88+HjRvhzTfhf/8Lhqtu3vwJHTp0\noEKFCmzatIm5c+cyZ84cGjRooHAQkRyhgIhCa9ZA8+ZQvXrQfLRmDYwdC5UrO++88w633XYbderU\noUSJEqxbt47Ro0dz6aWXhl22iOQzamKKEu7BNp7/+U+wiF6nTjB4cNDXsG/fPl5+eRp9+/Zlx44d\ndO7cmQkTJnDKKaeEXbaI5GPm7mHXkCVm5rFW85Hs2xesidSnTzB/4dFH4Z574MQTYdu2bQwfPpxB\ngwZRvnx5unTpws0330yhQnrwE5GsMTPcPUvtz3qCCMnOnTBiRDBn4bzz4F//ghtuCHZt27hxI/37\n92fy5MnccsstzJgxgxo1aoRdsogUMAqIPPbNNzBgQBAODRoEez7XrBnMX3jrrbfo168fK1asoG3b\ntqxfv56yZcuGXbKIFFChBISZjQRuBpLdvUrkWElgClAO+Bxo5u47wqgvN6xbF+y5MGNG0AG9fDmU\nLx8snDds2AQGDBgAQOfOnXn11Ve11LaIhC6sxuzRQKPDjnUH3nL3i4D5wON5XlUOc4f58+HGG4On\nhfLlg6UwBgyA44//km7dulGuXDlmz55N//79SUpKonXr1goHEYkKoQSEuy8Gth12uDEwNvL1WKBJ\nnhaVg/bvh8mT4fLL4aGHoGlT+PxzePJJ56OPFhMfH0+1atVITU3lvffeY9asWZq/ICJRJ5r6IEq7\nezKAu39nZqXDLiirDnY89+8fLJbXsyfcdBPs3ZvCyy+/zMCBA/n555/p2LEjI0eOpHjx4mGXLCKS\noWgKiMPFzFjWr7/+fcfzlClQuzZs2bKF//u/IYwYMYLLLruMf/3rXzRq1EjDVEUkJkRTQCSb2Znu\nnmxmZYCtGb2xZ8+ev34dFxdHXFxc7leXjjVrgo7n2bOhRQtYuRLOPTeY7RwfP5D58+fTvHlzFi9e\nzIUXXhhKjSJSMCUmJpKYmHhM5whtopyZnQv8190vjXzfG/jJ3XubWTegpLt3T+dzoU6Ucw/WROrb\nFz78EDp2hLZt4YQTdjFx4kQGDRrE3r176dChAy1bttQy2yISFbIzUS6UgDCzSUAcUApIBnoAM4BX\ngLOBLwiGuW5P57OhBERqKkyaBC+8EHzfpUsw4/nLLzfx0ksvMX78eOrVq8fDDz+sDmcRiToxM5Pa\n3e/J4KVr87SQTPjpJ0hIgEGDgu08+/SBBg0O8MYbc2jSZDCrVq2idevWrFq1inLlyoVdrohIjomm\nPoio8sknwTIYEydC48bwxhtQtuz3jBo1inbtEjjjjDNo374906dPp0iRImGXKyKS4zSc5hAH93hu\n2hTq1IFixeDDD5127ZbSp09LLrzwQjZu3MjUqVNZvnw5LVu2VDiISL6lgCCY2DZ1KlxxBdx7bzBU\ndd26XZx//ghuueUymjdvTtWqVdm8eTOjRo2iZs2aYZcsIpLrCnQT088/w8iRwcS2s8+Gxx+H889f\nz/DhCVx88UTq1q3Ls88+y3XXXae5CyJS4BTIgPj882B/5zFj4LrrYOLEvWzZ8hr9+iWwceNGHnjg\nAVavXs0555wTdqkiIqEpUAGxdGkwTPXtt+H++2HmzM94/fXh3H77KCpVqsTDDz9MkyZNKFy4cNil\nioiELt+3m+zfD6+8AldeGcxbqFNnP4MHz+TDD2+gSZOa7Nmzh8TERObPn098fLzCQUQkIt9uObp9\ne9C/MHAgnHUWtGz5FV9/PZKRI4dTrlw52rZtS3x8vJbWFpECITsT5fLdE8TmzcHyF+XLw4oVB/j7\n31/n9NOb0L17FX74YSuvv/46S5YsoWXLlgoHEZEjyBdPEO6waBG8+GIwj+HOO7/hpJNG8sorIzjz\nzDNp06YNd911F0WLFg2pahGRcMXMWkzH4tCASE2Fl18OZjzv3n2Aa699ky+/HMqSJYu48847adOm\nDdWrVw+5YhGR8BWYgPjuOychAYYMgYoVv+Lss0exaNFISpcuTdu2bfW0ICJymJhZrO9YXXTRfq64\nYg4VKw5j7dolXHzxnUyfPp0aNWqEXZqISL4RkwFx8snnsmPH2fztb3/jzjtf5pRTTgm7JBGRfCcm\nm5jWrFlDlSpVwi5FRCRmFJg+iFirWUQkbJoHISIiOUYBISIi6VJAiIhIuhQQIiKSLgWEiIikSwEh\nIiLpUkCIiEi6FBAiIpIuBYSIiKRLASEiIulSQIiISLoUECIikq6oCggzu97MNpjZx2bWLex6REQK\nsqgJCDMrBAwCGgGXAHebWcVwq8p5iYmJYZdwTFR/uGK5/liuHWK//uyImoAAagGb3P0Ld98HvAw0\nDrmmHBfrf8lUf7hiuf5Yrh1iv/7siKaA+DOw5ZDvv4ocExGREERTQIiISBSJmh3lzOwKoKe7Xx/5\nvjvg7t77sPdFR8EiIjEmZrccNbPjgI1AA+BbYDlwt7t/FGphIiIF1PFhF3CQux8ws/bAmwRNXyMV\nDiIi4YmaJwgREYkuMdVJbWafm9kaM1ttZsvDrudozGykmSWb2dpDjpU0szfNbKOZzTWzEmHWeCQZ\n1N/DzL4ys/cj/10fZo0ZMbOzzGy+ma0zsyQz6xg5HhP3P536O0SOx8r9P9HMlkX+rSaZWY/I8Vi5\n/xnVHxP3H4K5ZZEaZ0W+z/K9j6knCDP7FLjM3beFXUtmmNlVwC/AOHevEjnWG/jR3Z+PzBYv6e7d\nw6wzIxnU3wPY6e4vhFrcUZhZGaCMu39gZkWBVQTzau4jBu7/Eeq/kxi4/wBmdrK77470Ly4BOgK3\nEwP3HzKs/wZi5/53Bi4Dirv7rdn52RNTTxCAEUM1u/ti4PAwawyMjXw9FmiSp0VlQQb1Q/D/Iaq5\n+3fu/kHk61+Aj4CziJH7n0H9B+cFRf39B3D33ZEvTyTo73Ri5P5DhvVDDNx/MzsLuBEYccjhLN/7\nmPlhG+HAPDNbYWZ/C7uYbCrt7skQ/BAASodcT3a0N7MPzGxEtDYRHMrMzgWqAUuBM2Pt/h9S/7LI\noZi4/5EmjtXAd8A8d19BDN3/DOqH2Lj/LwJd+S3UIBv3PtYCoq671yBIxocjTSCxLnba+AIvAeXd\nvRrBP5yoftSONM+8CnSK/CZ++P2O6vufTv0xc//dPc3dqxM8udUys0uIofufTv2ViIH7b2Y3AcmR\nJ9AjPe0c9d7HVEC4+7eRP78HphOs3xRrks3sTPi1nXlryPVkibt/7791XA0HaoZZz5GY2fEEP1zH\nu/vMyOGYuf/p1R9L9/8gd/8ZSASuJ4bu/0GH1h8j978ucGukz3YycI2ZjQe+y+q9j5mAMLOTI79N\nYWanAA2BD8OtKlOM36f4LODeyNetgJmHfyDK/K7+yF+sg5oS3f8PRgHr3b3/Icdi6f7/of5Yuf9m\ndvrB5hczOwm4jqAfJSbufwb1b4iF++/uT7j7Oe5eHrgLmO/uLYD/ksV7HzOjmMzsPIKnBifoMJro\n7s+FW9WRmdkkIA4oBSQDPYAZwCvA2cAXQDN33x5WjUeSQf1/IWgPTwM+B9oebNeMJmZWF1gEJBH8\nnXHgCYIZ+lOJ8vt/hPrvITbu/6UEHaGFIv9Ncfd/m9lpxMb9z6j+ccTA/T/IzOoDj0RGMWX53sdM\nQIiISN6KmSYmERHJWwoIERFJlwJCRETSpYAQEZF0KSBERCRdCggREUlX1GwYJJIdkbHdbxPMEygL\nHAC+j3xfy933h1heuszsPuB/7h71s4ilYNM8CMk3zOwfwC/RsBSzmRVy97QMXnsHaO/ua7JwvuPc\n/UCOFSiSCWpikvzkdwuTmVnLyKYv75vZoMix48xsm5n1NbMPzWyOmdUys0Qz++TgBjBm1trMXosc\n32hmT2byvC+a2QdATTPraWbLzWytmb0UeV8zgpm4L0c+X9jMtphZ8cjrtc1sXuTrXmY21swWA6Mj\n1+hrZksjq4neH3nfn8zsncj51prZFbl9o6VgUEBIvhRZOfQ2oE5kBeDCZnZX5OUSBE08lYF9BEuI\nXAM0A3odcpqawK1AdeAeM6uSifMmuns1d18G9HP3WpHNlk41s0buPhX4gGCZgxruvo8jr3B6EfAX\nd28JtCFYpfMKgoUq20fW/W8OzIrUUxVYi0gOUB+E5FfXApcDK83MgCIE688A7Hb3+ZGvk4Dt7p5m\nZklAuUPOMTeykidmNh24Cih8hPOmHrJqLMB1ZvZo5D2lgJXA3Mhrhz7tHGlJ5pmREIFggcqKZnZ3\n5PviwAXACmComRWJvF8BITlCASH5lQGj3L3H7w4G20fuPeRQGpB6yNeH/ps49Dd5O+T7jM6755Dv\nTwIGAtXc/Tsz60UQFOnZz29P84e/Z9dhNTzk7gsOP0FkUbabgHFm1tvdJ2dwLZFMUxOT5FdvAc3M\nrBQEo50izTFw5N/YD32toZkVN7OTCbZrXEIwYioz5z2JYETVj2ZWjGAv5oN2Evz2f9BnBHsHc9j7\nDjeXYKOs4yLXvtDMipjZOQRNTyOA0QRNYiLHTE8Qki+5+4dm9k/gLTMrRPDU0A74liPvpHXoaysI\n9i8oC4w52HSTmfO6+09mNpZgD4RvCLY7PWg0MMLMdhP0JfwTGG5m2wiW+M7IUOAc4IOgdYtkguBq\nAHQxs33Az0CLI5xDJNM0zFUkHWbWGrjE3buEXYtIWNTEJCIi6dIThIiIpEtPECIiki4FhIiIpEsB\nISIi6VJAiIhIuhQQIiKSLgWEiIik6/8BxBJBlAAbcGIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x17a4a7128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "temperatures = np.linspace(5,40)\n",
    "trots = [trot_radex(tkin=t).to(u.K).value for t in temperatures]\n",
    "pl.plot(temperatures, trots)\n",
    "pl.plot(temperatures, tr_swift(temperatures), color='k')\n",
    "pl.xlabel(\"Temperatures\")\n",
    "pl.ylabel(\"$T_{rot} (2-2)/(1-1)$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tests of cold_ammonia reproducing pyspeckit ammonia spectra ##"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: No header given.  Creating an empty one.\n",
      "WARNING: No header given.  Creating an empty one.\n",
      "WARNING: No header given.  Creating an empty one.\n",
      "INFO:astropy:Creating spectra\n",
      "INFO:astropy:Concatenating data\n",
      "INFO:astropy:Left region selection unchanged.  xminpix, xmaxpix: 0,753\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: Creating spectra [pyspeckit.spectrum.classes]\n",
      "INFO: Concatenating data [pyspeckit.spectrum.classes]\n",
      "INFO: Left region selection unchanged.  xminpix, xmaxpix: 0,753 [pyspeckit.spectrum.interactive]\n",
      "For Tkin=20 -> Trot=22.80, pyspeckit's cold_ammonia fitter got:\n",
      "Param #0        tkin0 =      27.1194 +/-         140.857   Range:[2.7315,inf)\n",
      "Param #1         tex0 =      6.65999 +/-         2721.09   Range:[2.7315,inf)\n",
      "Param #2        ntot0 =           13 +/-         128.091   Range:    [5,25]\n",
      "Param #3       width0 =          0.5 +/-         2.66699   Range:   [0,inf)\n",
      "Param #4      xoff_v0 = -2.94186e-05 +/-               0 \n",
      "Param #5      fortho0 =            0 (fixed)  Range:     [0,1]\n"
     ]
    }
   ],
   "source": [
    "from pyspeckit.spectrum.models.tests import test_ammonia\n",
    "from pyspeckit.spectrum.models import ammonia\n",
    "spc = test_ammonia.make_synthspec(lte=False, tkin=None, tex=6.66, trot=22.80)\n",
    "spc.specfit.Registry.add_fitter('cold_ammonia',ammonia.cold_ammonia_model(),6)\n",
    "spc.specfit(fittype='cold_ammonia', guesses=[23, 5, 13.1, 1, 0.5, 0],\n",
    "            fixed=[False,False,False,False,False,True])\n",
    "\n",
    "print(\"For Tkin=20 -> Trot=22.80, pyspeckit's cold_ammonia fitter got:\\n{0}\".format(spc.specfit.parinfo))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:astropy:Left region selection unchanged.  xminpix, xmaxpix: 0,753\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: Left region selection unchanged.  xminpix, xmaxpix: 0,753 [pyspeckit.spectrum.interactive]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Param #0        tkin0 =      27.1194 +/-      7.2437e-05   Range:[2.7315,inf)"
      ]
     },
     "execution_count": 120,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spc.specfit(fittype='cold_ammonia', guesses=[22.80, 6.6, 13.1, 1, 0.5, 0],\n",
    "            fixed=[False,False,False,False,False,True])\n",
    "spc.specfit.parinfo[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}