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



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; this argument is used for older
          verisons of cloudy (pre-C16), which output the emission
          as a flux per unit area at the cloud inner radius, so the
          inner radius was required to convert this to a luminosity;
          more recent versions of cloudy output the luminosity
          directly, so this is not needed. If left as None, the
          outputs are assumed to be from C16 or later, and so are
          interpreted as luminosities

    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 frequency to wavelength
    nu = en * physcons.Rydberg*physcons.c  # Rydbergs to Hz
    wl = 1e8 * c/nu[::-1]                  # Wavelength in Angstrom
    emission = emission[:,::-1]            # Reverse order to match wl

    # Convert from nu L_nu to L_lambda or 
    L_lambda = emission / wl
    if r0 is not None:
        L_lambda *= 4.0 * np.pi * r0**2

    # Save and return output
    out_type = namedtuple('cloudy_continuum',
                              ['wl', 'L_lambda'])
    out = out_type(wl, L_lambda)        
    return out