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 .