py: photo ionisation and heating rates for given UV spectral slope

calculate_rates.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""Calculate photo ionisation and photo heating rates for J21 = 1.0 and a given
UV spectral slope value, α.

Sorry for the current lack of comments!
"""

__author__ = "Simon Mutch"
__date__ = "2019-06-07"

import numpy as np
from scipy.integrate import quad
from astropy import units as U, constants as C
from box import Box
import click

# Using parameters from Table 1 of Verner+ (1996)

params_HI = Box(
    dict(Z=1, N=1, Eth=1.360e1, Emax=5.000e4, E0=4.29e-1, sigma0=5.475e4, ya=3.288e1, P=2.963, yw=0.0, y0=0.0, y1=0.0),
    frozen=True,
)

params_HeI = Box(
    dict(
        Z=2,
        N=2,
        Eth=2.459e1,
        Emax=5.000e4,
        E0=1.361e1,
        sigma0=9.492e2,
        ya=1.469,
        P=3.188,
        yw=2.039,
        y0=4.434e-1,
        y1=2.136,
    ),
    frozen=True,
)

params_HeII = Box(
    dict(Z=2, N=1, Eth=5.442e1, Emax=5.000e4, E0=1.720, sigma0=1.369e4, ya=3.288e1, P=2.693, yw=0.0, y0=0.0, y1=0.0),
    frozen=True,
)

params_HeIII = Box(
    dict(Z=2, N=1, Eth=5.442e1, Emax=5.000e4, E0=1.720, sigma0=1.369e4, ya=3.288e1, P=2.693, yw=0.0, y0=0.0, y1=0.0),
    frozen=True,
)


def _x(E: U.Quantity, params: dict):
    return E / (params.E0 * U.eV) - params.y0


def _F(x: U.Quantity, y: U.Quantity, params: dict):
    return ((x - 1.0) ** 2 + params.yw ** 2) * y ** (0.5 * params.P - 5.5) * (1 + np.sqrt(y / params.ya)) ** (-params.P)


def _y(x: U.Quantity, params: dict):
    return np.sqrt(x * x + params.y1 ** 2)


def sigma(ν: U.Quantity, params: dict):
    E = (C.h * ν).to("eV")
    x = _x(E, params)
    y = _y(x, params)
    F = _F(x, y, params)

    return params.sigma0 * U.Unit("1e-18 cm2") * F


def ion_rate(alpha: float, params: dict):
    ν_HI = (params_HI.Eth * U.eV / C.h).to("1e6 GHz").value
    ν_0 = (params.Eth * U.eV / C.h).to("1e6 GHz").value

    def E_i(ν: float):
        J = (ν / ν_HI) ** -alpha
        return J / ν * sigma(ν * U.GHz * 1e6, params).value

    E, err = quad(E_i, ν_0, np.inf)
    E *= U.Unit("1e-39 erg")

    return (4.0 * np.pi / C.h * E).decompose()


def photoheat_rate(alpha: float, params: dict):
    ν_HI = (params_HI.Eth * U.eV / C.h).to("1e6 GHz").value
    ν_0 = (params.Eth * U.eV / C.h).to("1e6 GHz").value

    def Edot_i(ν: float):
        J = (ν / ν_HI) ** -alpha
        return J / ν * (ν - ν_0) * sigma(ν * U.GHz * 1e6, params).value

    Edot, err = quad(Edot_i, ν_0, np.inf)
    Edot *= U.Unit("1e-33 erg GHz")

    return (4.0 * np.pi * Edot).to("eV s-1")


@click.command()
@click.argument("α", type=click.FLOAT)
def cli(α):
    for name, params in (("HI", params_HI), ("HeI", params_HeI), ("HeII", params_HeII)):
        print(f"Γ_{name}(α={α}) = {ion_rate(α, params):.3e}")
    for name, params in (("HI", params_HI), ("HeI", params_HeI), ("HeII", params_HeII)):
        print(f"g_{name}(α={α}) = {photoheat_rate(α, params):.3e}")


if __name__ == "__main__":
    cli()