Updating lines structure in the grid to give correct continua

[WillJRoper]

This PR restructures the line calculation on a grid.

• Lines are now stored under the same keys as spectra with continua derived from the spectra.
• Spectra are always read now.
• Line luminosities, wavelengths and continua are stored in separate attributes on a Grid for better "ergonomics" internally.
• A plot has been introduced to the galaxy_lines docs to demonstrate the different line luminiosities/continua for different emission models.

I'm not 100% sure this works yet and the luminosities appear to be very wrong in the their magnitude but I'm raising for comment since I'll be off the rest of the day.

Issue Type

• Bug
• Enhancement

Checklist

• [x] I have read the [CONTRIBUTING.md]() -->
• [x] I have added docstrings to all methods
• [x] I have added sufficient comments to all lines
• [x] I have made corresponding changes to the documentation
• [x] My changes generate no pep8 errors
• [x] I have added tests that prove my fix is effective or that my feature works
• [x] New and existing unit tests pass locally with my changes

[stephenmwilkins]

Something still isn't working right.

• The escaped line luminosity (which should be zero) is the same as the nebular for some reason.
• It's also confusing that while the escaped line luminosity is non-zero, the EW is zero. I thought the EW was calculated when the Line object is initialised?
• The nebular and reprocessed continuum are identical but reprocessed should be the sum of the nebular and transmitted?!
• The intrinsic EW should be non-zero, at the moment it looks like it's just a copy of escaped.

This is:

from synthesizer.emission_models import IntrinsicEmission

emission_model = IntrinsicEmission(
    grid=grid,
    fesc=0.2,
    )

line_ids = [
    ','.join([line_ratios.Hb, line_ratios.O3b, line_ratios.O3r]),
    line_ratios.Ha]

galaxy.stars.get_lines(line_ids, emission_model)

[stephenmwilkins]

I think this is the part I don't understand.

        # Now we have the spectra we can get the emission at each line
        # and include it
        lines[this_model.label] = {}
        for line_id in line_ids:
            # Get the emission at this lines wavelength
            lam = lines[this_model.lum_intrinsic_model.label][
                line_id
            ].wavelength

            # Get the luminosity at this wavelength
            luminosity = spectra.get_lnu_at_lam(lam)

            # Create the line (luminoisty = continuum)
            lines[this_model.label][line_id] = Line(
                line_id=line_id,
                wavelength=lam,
                luminosity=luminosity * Hz,
                continuum=luminosity,
            )

Perhaps I'm missing something but isn't this just setting the line luminosity (luminosity) to the same as the continuum, multiplied by Hz?

However, a quick test reveals this is not the cause of the issue I'm seeing.

[stephenmwilkins]

Something still isn't working right.

• The escaped line luminosity (which should be zero) is the same as the nebular for some reason.
• It's also confusing that while the escaped line luminosity is non-zero, the EW is zero. I thought the EW was calculated when the Line object is initialised?
• The nebular and reprocessed continuum are identical but reprocessed should be the sum of the nebular and transmitted?!
• The intrinsic EW should be non-zero, at the moment it looks like it's just a copy of escaped.

This is:

from synthesizer.emission_models import IntrinsicEmission

emission_model = IntrinsicEmission(
    grid=grid,
    fesc=0.2,
    )

line_ids = [
    ','.join([line_ratios.Hb, line_ratios.O3b, line_ratios.O3r]),
    line_ratios.Ha]

galaxy.stars.get_lines(line_ids, emission_model)

Ah... but when I just use EscapedEmission it works as expected.

Is the combination doing something weird like updating escaped with intrinsic in someway?

[stephenmwilkins]

Something still isn't working right.

• The escaped line luminosity (which should be zero) is the same as the nebular for some reason.
• It's also confusing that while the escaped line luminosity is non-zero, the EW is zero. I thought the EW was calculated when the Line object is initialised?
• The nebular and reprocessed continuum are identical but reprocessed should be the sum of the nebular and transmitted?!
• The intrinsic EW should be non-zero, at the moment it looks like it's just a copy of escaped.

This is:

from synthesizer.emission_models import IntrinsicEmission

emission_model = IntrinsicEmission(
    grid=grid,
    fesc=0.2,
    )

line_ids = [
    ','.join([line_ratios.Hb, line_ratios.O3b, line_ratios.O3r]),
    line_ratios.Ha]

galaxy.stars.get_lines(line_ids, emission_model)

Ah... but when I just use EscapedEmission it works as expected.

Is the combination doing something weird like updating escaped with intrinsic in someway?

Replacing:

# Loop over lines copying over the first set of lines
        for line_id in line_ids:
            lines[this_model.label][line_id] = copy.copy(
                lines[this_model.combine[0].label][line_id]
            )

with

        # Loop over lines copying over the first set of lines
        for line_id in line_ids:

            lum = 0.
            cont = 0.

            # Combine the lines
            for combine_model in this_model.combine[0:]:
                lum += lines[combine_model.label][line_id]._luminosity
                cont += lines[combine_model.label][line_id]._continuum

            lines[this_model.label][line_id] = Line(
                line_id=line_id,
                wavelength=lines[this_model.combine[0].label][line_id].wavelength,
                luminosity=lum*erg/s,
                continuum=cont*erg/s/Hz)

makes things make more sense now.

type line_luminosity continuum_luminosity equivalent_width
escaped -inf 28.35 0.00 Å
transmitted -inf 28.35 0.00 Å
nebular 42.89 28.51 3443.68 Å
reprocessed 42.89 28.74 2034.58 Å
intrinsic 42.89 28.89 1443.80 Å

However, there is a remaining issues that for since this model assumes fesc=0.2 escaped and transmitted should not be identical.

[WillJRoper]

I think this is the part I don't understand.

        # Now we have the spectra we can get the emission at each line
        # and include it
        lines[this_model.label] = {}
        for line_id in line_ids:
            # Get the emission at this lines wavelength
            lam = lines[this_model.lum_intrinsic_model.label][
                line_id
            ].wavelength

            # Get the luminosity at this wavelength
            luminosity = spectra.get_lnu_at_lam(lam)

            # Create the line (luminoisty = continuum)
            lines[this_model.label][line_id] = Line(
                line_id=line_id,
                wavelength=lam,
                luminosity=luminosity * Hz,
                continuum=luminosity,
            )

Perhaps I'm missing something but isn't this just setting the line luminosity (luminosity) to the same as the continuum, multiplied by Hz?

However, a quick test reveals this is not the cause of the issue I'm seeing.

This is indeed wrong. I didn't really think about it when I was writing it. Any dust emission only contributes to the continuum so the luminosity should be 0.

[WillJRoper]

Your issue that escaped and transmitted are the same is because of how the lines are defined in the grid generation this is guaranteed. There are no "transmitted" lines and since escaped is derived from transmitted these are guaranteed to be the same for lines.

This needs. ironing out at the grid level.

[stephenmwilkins]

Your issue that escaped and transmitted are the same is because of how the lines are defined in the grid generation this is guaranteed. There are no "transmitted" lines and since escaped is derived from transmitted these are guaranteed to be the same for lines.

This needs. ironing out at the grid level.

I don't think that's the case. When I use just the EscapedEmission and ReprocessedEmission independently I think it made sense. I assumed the issue was something to do with the combining linking the two somehow.

[WillJRoper]

Your issue that escaped and transmitted are the same is because of how the lines are defined in the grid generation this is guaranteed. There are no "transmitted" lines and since escaped is derived from transmitted these are guaranteed to be the same for lines. This needs. ironing out at the grid level.

I don't think that's the case. When I use just the EscapedEmission and ReprocessedEmission independently I think it made sense. I assumed the issue was something to do with the combining linking the two somehow.

Nah, you can see in the source:

class EscapedEmission(StellarEmissionModel):
    """
    An emission model that extracts the escaped radiation field.

    This defines an extraction of key "escaped" from a SPS grid.

    This is a child of the EmissionModel class for a full description of the
    parameters see the EmissionModel class.

    Note: The escaped model is the "mirror" to the transmitted model. What
    is transmitted is not escaped and vice versa. Therefore
    EscapedEmission.fesc = 1 - TransmittedEmission.fesc. This will be
    checked to be true at the time of spectra generation.

    Attributes:
        grid (synthesizer.grid.Grid): The grid object to extract from.
        label (str): The label for this emission model.
        extract (str): The key to extract from the grid.
        fesc (float): The escape fraction of the emission. Note that this will
                        be set to 1 - fesc during generation, i.e. it is the
                        fraction that does not escape that should be passed.
    """

    def __init__(
        self,
        grid,
        label="escaped",
        fesc=0.0,
        **kwargs,
    ):
        """
        Initialise the EscapedEmission object.

        Args:
            grid (synthesizer.grid.Grid): The grid object to extract from.
            label (str): The label for this emission model.
            fesc (float): The escape fraction of the emission. (Note that,
                          1-fesc will be used during generation).
        """
        StellarEmissionModel.__init__(
            self,
            grid=grid,
            label=label,
            extract="transmitted",
            fesc=(1 - fesc),
            **kwargs,
        )

i.e. it's the inversion of "transmitted" because that is how escaped is defined (since Malta). Since we attach lines to the nebular model (which makes sense IMO) you will never get transmitted being non-zero by construction and thus transmitted=escaped.

I'm not arguing for whether this is correct, this is just "how it is" based on the definitions we use. It would be better if we did have lines for each extraction component (incident, nebular, transmitted) but I'm guessing cloudy doesn't work like that.

[WillJRoper]

Agreed that the combination wasn't working properly. That has now been changed.

[stephenmwilkins]

There are two things here... leaving aside the lines.

• First, I'm fairly sure the extract for escaped should be incident and not transmitted. However, for most lines that will make very little difference.

• Second, if I choose fesc = 0.2 and calculate IntrinsicEmission the continuum of escaped and transmitted are identical. However, transmitted should be about 4x escaped.

[stephenmwilkins]

Unfortunately, I need to leave soon. However, it looks, in my example code that actually fesc is doing nothing to the continuum only the lines. That would cause the issue.

[WillJRoper]

Ok. Well your first point could be true (although I'm sure there's a reason it was decided it should be transmitted). This all boils down to exactly what escape fractions are being defined. Is it the escape of the incident or the escape of the incident reprocessed by the gas immediately around the star?

Second, if I choose fesc = 0.2 and calculate IntrinsicEmission the continuum of escaped and transmitted are identical. However, transmitted should be about 4x escaped.

After fixing the combination that is true:

escaped/H 1 4861.32A: lum=0.0 erg/s cont=2.5761699526441095e+27 erg/(Hz*s)
transmitted/H 1 4861.32A: lum=0.0 erg/s cont=1.0304679810576438e+28 erg/(Hz*s)

[WillJRoper]

Unfortunately, I need to leave soon. However, it looks, in my example code that actually fesc is doing nothing to the continuum only the lines. That would cause the issue.

Fair, but see the message above. fesc is now effecting the continuum as expected and all I've done is fix the combination operation.

[WillJRoper]

Turns out that the parametric generate_line function wasn't applying fesc to the continuum. This is also now fixed and I think everything is working. These are the plots the docs generates:

[stephenmwilkins]

Code wise it all looks fine and it passes all my current sanity checks so I think it can safely be merged. However, we might discover something else when writing the lines the paper but that can be addressed later.

[WillJRoper]

Agreed, I'm 90% sure there will be a gremlin but let us face it when someone feeds it after midnight.