Source code for slugpy.cloudy.read_cloudy_continuum

"""
Routine to read a cloudy continuum output, produced by save last continuum
"""

import numpy as np
from collections import namedtuple
import scipy.constants as physcons
try:
    c = physcons.c * 1e2   # c in CGS units
except:
    # This exception is here to deal with readthedocs not having scipy
    c = 3.0e10

[docs]def read_cloudy_continuum(filename, r0=None): """ Reads a cloudy continuum output, produced by save last continuum Parameters filename : string name of the file to be read r0 : float inner radius, in cm; if included, the quantities returned will be total energies instead of energy emission rates instead of rates per unit area Returns A namedtuple containing the following fields: wl : array wavelengths in Angstrom incident : array incident radiation field intensity """ # Open file fp = open(filename, 'r') # Burn the header line fp.readline() # Prepare storage en = [] emission = [] line_label = [] cont_label = [] line_density = [] # Read the data for line in fp: # Skip comment lines that start with # if line.strip()[0] == '#': continue linesplit = line.split() en.append(float(linesplit[0])) incident = float(linesplit[1]) transmitted = float(linesplit[2]) emitted = float(linesplit[3]) transmitted_plus_emitted = float(linesplit[4]) emission.append([incident, transmitted, emitted, transmitted_plus_emitted]) fp.close() # Convert to arrays en = np.array(en) emission = np.transpose(np.array(emission)) # Convert units to match convension in SLUG -- emission is # L_lambda, units of erg / s / Ang, flux is F_lambda, units of erg # / cm^2 / s / Ang; cloudy gives us nu L_nu / 4 pi r_0^2, units of # erg / s / cm^2 nu = en * physcons.Rydberg*physcons.c # Rydbergs to Hz F_nu = emission / nu # Convert nu L_nu to L_nu F_lambda = F_nu * nu**2/c / 1e8 # Convert L_nu to L_lambda, # and units to erg / s / Ang wl = 1e8 * c/nu[::-1] # Wavelength in Angstrom F_lambda = np.copy(F_lambda[:,::-1]) # Reverse order if r0 is not None: L_lambda = F_lambda * 4.0*np.pi*r0**2 out_type = namedtuple('cloudy_continuum', ['wl', 'L_lambda']) out = out_type(wl, L_lambda) else: out_type = namedtuple('cloudy_continuum', ['wl', 'F_lambda']) out = out_type(wl, F_lambda) # Return output return out