Source code for slugpy.cloudy.read_cloudy_hcon

"""
Routine to ready cloudy hydrogen conditions and hydrogen ionization
files, produced by 'save last hydrogen conditions' and 'save last
hydrogen ionization'.
"""

import numpy as np
try:
    # Various constants in cgs
    c = physcons.c * 1e2
    a0 = physcons.physical_constants['Bohr radius'][0]*1e2
    alpha = physcons.alpha
    sigma_PI = 2.**9*np.pi**2/(3.0*np.exp(1.0)**4)*alpha*a0**2
except:
    # This exception is here to deal with readthedocs not having scipy
    c = 3.0e10
    sigma_PI = 6.3e-18
fi = 1.1

[docs]def read_cloudy_hcon(hcon_file, r0 = 0.0): """ Reads cloudy outputs produce by the 'save last hydrogen conditions' and 'save last hydrogen ionization' file, and uses these to return various HII region diagnostic parameters. Parameters hcon_file : str Name of hydrogen conditions file to be read r0 : float Inner radius for the calculation Returns r1 : float outer radius, in cm nII : float average number density of H nuclei Omega : float wind parameter, defined as r0^3 / (r1^3 - r0^3) Notes All averages are averages over the ionized volume, i.e., the average is taken with a weighting factor 4 pi r^2 (n_H+/n_H) dV. """ # Load the hydrogen conditions file data = np.atleast_2d(np.loadtxt(hcon_file)) # Form the radii and the weighting factor for volume integrations; # be sure to handle pathological case of only 1 or 2 zones of # cloudy output r = data[:,0] + r0 if data.shape[0] >= 3: dr = 0.5*(data[2:,0] - data[:-2,0]) dr = np.append(dr, data[-1,0]-data[-2,0]) dr = np.insert(dr, 0, 0.5*(r[1]+r[0])-r0) elif data.shape[0] == 2: dr = [ 0.5*(r[1]+r[0])-r0, data[1,0] - data[0,0] ] else: dr = 1.0 # Doesn't matter in this case, since there's # nothing to weight w = r**2*dr*data[:,4] wsum = np.sum(w) # Get r1 and Omega r1 = r[-1] Omega = r0**3 / (r1**3 - r0**3) # Get average number density nII = np.sum(data[:,2]*w)/wsum # Return return r1, nII, Omega