astropy / astroquery
1 1
import numpy as np
2 1
import re
3 1
from astropy import constants
4 1
from astropy import units as u
5

6

7 1
def ncrit(lamda_tables, transition_upper, transition_lower, temperature, OPR=3,
8
          partners=['H2', 'OH2', 'PH2']):
9
    """
10
    Compute the critical density for a transition given its temperature.
11

12
    The critical density is defined as the Einstein A value divided by the sum
13
    of the collision rates into the state minus the collision rates out of that
14
    state.  See Shirley et al 2015, eqn 4
15
    (http://esoads.eso.org/cgi-bin/bib_query?arXiv:1501.01629)
16

17
    Parameters
18
    ----------
19
    lamda_tables : list
20
        The list of LAMDA tables returned from a Lamda.query operation.
21
        Should be [ collision_rates_dict, Avals/Freqs, Energy Levels ]
22
    transition_upper : int
23
        The upper transition number as indexed in the lamda catalog
24
    transition_lower: int
25
        The lower transition number as indexed in the lamda catalog
26
    temperature : float
27
        Kinetic temperature in Kelvin.  Will be interpolated as appropriate.
28
        Extrapolation uses nearest value
29
    OPR : float
30
        ortho/para ratio of h2 if para/ortho h2 are included as colliders
31
    partners : list
32
        A list of valid partners.  It probably does not make sense to include
33
        both electrons and H2 because they'll have different densities.
34

35
    Returns
36
    -------
37
    ncrit : astropy.units.Quantity
38
        A quantity with units cm^-3
39
    """
40

41 0
    fortho = (OPR) / (OPR - 1)
42

43
    # exclude partners that are explicitly excluded
44 0
    crates = {coll: val for coll, val in lamda_tables[0].items()
45
              if coll in partners}
46 0
    avals = lamda_tables[1]
47 0
    enlevs = lamda_tables[2]
48

49 0
    aval = avals[(avals['Upper'] == transition_upper) &
50
                 (avals['Lower'] == transition_lower)]['EinsteinA'][0]
51

52 0
    temperature_re = re.compile(r"C_ij\(T=([0-9]*)\)")
53 0
    crate_temperatures = np.array(
54
        [int(temperature_re.search(cn).groups()[0])
55
         for cn in crates[list(crates.keys())[0]].keys()
56
         if temperature_re.search(cn)])
57

58 0
    if temperature < crate_temperatures.min():
59 0
        crates_ji_all = {
60
            coll: cr['C_ij(T={0})'.format(crate_temperatures.min())]
61
            for coll, cr in crates.items()}
62 0
    elif temperature > crate_temperatures.max():
63 0
        crates_ji_all = {
64
            coll: cr['C_ij(T={0})'.format(crate_temperatures.max())]
65
            for coll, cr in crates.items()}
66 0
    elif temperature in crate_temperatures:
67 0
        crates_ji_all = {coll: cr['C_ij(T={0})'.format(temperature)]
68
                         for coll, cr in crates.items()}
69
    else:  # interpolate
70 0
        nearest = np.argmin(np.abs(temperature - crate_temperatures))
71 0
        if crate_temperatures[nearest] < temperature:
72 0
            low, high = (crate_temperatures[nearest],
73
                         crate_temperatures[nearest + 1])
74
        else:
75 0
            low, high = (crate_temperatures[nearest - 1],
76
                         crate_temperatures[nearest])
77 0
        crates_ji_all = {coll:
78
                         (cr['C_ij(T={0})'.format(high)] -
79
                          cr['C_ij(T={0})'.format(low)]) *
80
                         (temperature - low) / (high - low) +
81
                         cr['C_ij(T={0})'.format(low)]
82
                         for coll, cr in crates.items()}
83

84 0
    transition_indices_ji = {
85
        coll: np.nonzero(cr['Upper'] == transition_upper)[0]
86
        for coll, cr in crates.items()}
87

88 0
    crates_ji = {coll: crates_ji_all[coll][transition_indices_ji[coll]]
89
                 for coll in crates}
90

91
    # i > j: collisions from higher levels
92 0
    transition_indices_ij = {
93
        coll: np.nonzero(cr['Lower'] == transition_upper)[0]
94
        for coll, cr in crates.items()}
95 0
    crates_ij = {}
96 0
    for coll in crates.keys():
97 0
        crates_ind = crates[coll][transition_indices_ij[coll]]
98 0
        degeneracies_i = enlevs['Weight'][crates_ind['Upper'] - 1]
99 0
        degeneracies_j = enlevs['Weight'][crates_ind['Lower'] - 1]
100 0
        energy_i = enlevs['Energy'][crates_ind['Upper'] - 1] * u.cm ** -1
101 0
        energy_j = enlevs['Energy'][crates_ind['Lower'] - 1] * u.cm ** -1
102
        # Shirley 2015 eqn 4:
103 0
        crates_ij[coll] = (
104
            crates_ji_all[coll][transition_indices_ij[coll]] *
105
            degeneracies_i / degeneracies_j.astype('float') *
106
            np.exp((-energy_i - energy_j).to(u.erg, u.spectral()) /
107
                   (constants.k_B * temperature * u.K)))
108

109 0
    crates_tot_percollider = {
110
        coll: (np.sum(crates_ij[coll]) + np.sum(crates_ji[coll])) *
111
        u.cm ** 3 / u.s for coll in crates}
112

113 0
    if 'OH2' in crates:
114 0
        crates_tot = (fortho * crates_tot_percollider['OH2'] +
115
                      (1 - fortho) * crates_tot_percollider['PH2'])
116 0
    elif 'PH2' in crates:
117 0
        crates_tot = crates_tot_percollider['PH2']
118 0
    elif 'H2' in crates:
119 0
        crates_tot = crates_tot_percollider['H2']
120

121 0
    return ((aval * u.s ** -1) / crates_tot).to(u.cm ** -3)

Read our documentation on viewing source code .

Loading