jzuhone / pyxsim

Simulating X-ray observations from astrophysical sources.
http://hea-www.cfa.harvard.edu/~jzuhone/pyxsim
Other
20 stars 8 forks source link

The emission_measure derived from yt #50

Closed Arppe closed 1 year ago

Arppe commented 1 year ago

Hi John,

I notice that the H_nuclei_density used to calculate the emission_measure field is fixed to primordial_H_mass_fraction = 0.76. However, some simulations (e.g. IllustrisTNG) indeed offer the correct values of the metals. Therefore, maybe we can get a better estimate of the emission_measure.

Thank you.

jzuhone commented 1 year ago

Hi @Arppe,

This is not the case in codes which track elemental abundances, such as TNG which you mentioned. This value of the H fraction is only used if there are no species fields defined in the dataset at all. For TNG (and similar situations), it should find the value of the H fraction and compute the emission measure field accordingly.

You may be looking at an older version of yt, pyXSIM, or both--if so, I highly recommend upgrading.

Arppe commented 1 year ago

I'm afraid that's the case.

https://github.com/yt-project/yt/blob/main/yt/frontends/arepo/fields.py

        if (ptype, "NeutralHydrogenAbundance") not in self.field_list:
            amu_cgs = self.ds.units.physical_constants.amu_cgs
            muinv = _primordial_mass_fraction["H"] / ChemicalFormula("H").weight

            def _h_number_density(field, data):
                return data["gas", "density"] * muinv / amu_cgs

            self.add_field(
                (ptype, "H_number_density"),
                sampling_type="particle",
                function=_h_number_density,
                units=self.ds.unit_system["number_density"],
            )
            self.alias(("gas", "H_number_density"), (ptype, "H_number_density"))
            self.alias(("gas", "H_nuclei_density"), ("gas", "H_number_density"))

https://github.com/yt-project/yt/blob/main/yt/fields/astro_fields.py

def _emission_measure(field, data):
    dV = data[ftype, "mass"] / data[ftype, "density"]
    nenhdV = data[ftype, "H_nuclei_density"] * dV
    nenhdV *= data[ftype, "El_number_density"]
    return nenhdV

registry.add_field(
    (ftype, "emission_measure"),
    sampling_type="local",
    function=_emission_measure,
    units=unit_system["number_density"],
)

https://github.com/yt-project/yt/blob/main/yt/utilities/physical_ratios.py

primordial_H_mass_fraction = 0.76

_primordial_mass_fraction = {
    "H": primordial_H_mass_fraction,
    "He": (1 - primordial_H_mass_fraction),
}
jzuhone commented 1 year ago

@Arppe Does your dataset not have NeutralHydrogenAbundance? Most TNG datasets should. This only happens if that field is not included. Immediately above the lines you quoted in yt/frontends/arepo/fields.py, it says this:

# If we have ElectronAbundance but not NeutralHydrogenAbundance, assume the
# cosmic value for hydrogen to generate the H_number_density

Otherwise, you get the TNG value for the H fraction.

jzuhone commented 1 year ago

There is a potential corner case here where you have the metal fractions defined in GFM_Metals but not the NeutralHydrogenAbundance, but this case is rare.

Arppe commented 1 year ago

Maybe I didn't notice the other defination of ("gas", "H_nuclei_density") anywhere else. I probably do not need to worry about it. Thanks.

By the way, is the temperature field derived from GFM_Metals or primordial_H_mass_fraction = 0.76?

jzuhone commented 1 year ago

The temperature is derived from GFM_Metals

jzuhone commented 1 year ago

The corner case I referred to above was resolved in https://github.com/yt-project/yt/pull/4419.