PAHFIT / pahfit

Model Decomposition for Near- to Mid-Infrared Spectroscopy of Astronomical Sources
https://pahfit.readthedocs.io/
18 stars 26 forks source link

Blackbody continuum sum incorrectly implemented #171

Closed jdtsmith closed 2 years ago

jdtsmith commented 2 years ago

Looking closely at the model implementation, it seems the sum of blackbody continua is incorrectly implemented. Here is the PAHFIT model in simple form:

httptir astro utoledo edujdsmithresearchdownloadpah_template_final_apj pdf

Notice the B_nu/lam^2, which is a sum of emissivity=2 MBB's. I believe the current model continuum sources are unmodified (which makes them "peakier" and thus a bit harder to define smooth continuum with).

This brings up the need to validate against IDL PAHFIT not just in terms of the fit quality, but also the fit parameters including the tau's in the above.

I will take a crack at fixing this, probably by adding a "Modified" flag to the Blackbody1D model component and an associated column to the scipacks. Or perhaps just make a ModifiedBlackBody1D derived model.

Can someone take a look and confirm my suspicion here?

jdtsmith commented 2 years ago

Update: on even closer look, it seems the BlackBody1D model is not in fact a Blackbody, but a MBB. This means that the starlight which is not supposed to be modified (stars don't have strong frequency dependence of their emissivity) is the one that is incorrect. Which is less severe, but still incorrect. I'll put together a PR.

karllark commented 2 years ago

Great that you have caught this bug. Yes, I remember just implementing on BlackBody. Definitely sounds like we need two. Or one BB and one MBB (maybe with the latter using the first?).

jdtsmith commented 2 years ago

I am evaluating this, which seems to work fine:

class ModifiedBlackBody1D(BlackBody1D):
    """
    Modified blackbody with an emissivity propoportional to nu^2
    """

    @staticmethod
    def evaluate(x, amplitude, temperature):
        return BlackBody1D.evaluate(x, amplitude, temperature) * ((9.7 / x) ** 2)
karllark commented 2 years ago

Can you comment on your implementation? In a previous comment you say that the current "BlackBody1D" is actually a modified blackbody. But then in your next comment you are modifying the currenct BlackBody1D. Feeling a bit confused.

jdtsmith commented 2 years ago

Yep I altered BlackBody1D to be an actual BB, then added an MBB which calls it.