"Marcin bug" bugfix

[aprsa]

Marcin discovered that extrapolation/blending can cause haywire intensities for high Teffs and low mus. This was the consequence of limb darkening extrapolation defaulting to ATM instead of LD extrapolation method, causing the intensities to go crazy. The fix is to use LD extrapolation instead of ATM extrapolation.

[Raindogjones]

I'm not sure I understand the fix. The options for both methods are the same, so is it a problem of mixing none, nearest and linear? Is it not just the case that for some regions the linear extrapolation will drive things into non-physical regimes?

[mwrona77]

It fixes the problem in the test case I ran. However, when I disable blending entirely using b.set_value_all('blending_method', 'none') to force the extrapolation method, the issue persists and I am still seeing extremely high intensities.

I also want to make sure I fully understand the underlying mechanics of this fix. Why is the atm_extrapolation_method parameter suddenly being assigned the value of ld_extrapolation_method? Why does atm_extrapolation_method not hold a valid value on its own here?

Additionally, why do we apply this change specifically to result = self.interpolate_imus and inorms_atm = self.interpolate_inorms in this block, but leave it unchanged in the other inorms_bb = self.interpolate_inorms call just above (around line 1653 in phoebe/atmospheres/passbands.py)?

[aprsa]

Extrapolation (and subsequently blending) rely on extrapolated atmospheric limb darkening coefficients to apply LD correction to blackbody intensities. The reason for this is that blackbody intensities have no inherent limb darkening so we need to appropriate model atmosphere limb darkening. While on-grid, this is straight-forward because the values are all defined. But when off-grid, we need to extrapolate and the question, then, is how. The default was linear extrapolation in LD space, which causes blow-up. LDs should always use 'nearest' extrapolation except when just marginally off-grid; the fix thus ought to work whenever ld_extrapolation_mode is set to 'nearest', irrespective of blending. This is especially pronounced when mu is small (close to the limb) and Teff is above grid values.

[mwrona77]

So, after the fix, when I tried to use extrapolation by setting:
b.set_value_all('blending_method', value="none")
the intensities still blow up.

However, when I went back to the pre-fix version and run the code with:
b.set_value_all('atm_extrapolation_method', value="nearest")
it works for both blending and extrapolation (no super-high intensities)

Also, in the fixed version, when someone decides to use:
b.set_value_all('ld_extrapolation_method', value="linear")
the problem with intensities remains.

Shouldn't we then just hard-code:
atm_extrapolation_method='nearest' instead of atm_extrapolation_method=atm_extrapolation_method in the places in the passband file where you changed it? Or we can just check if the out-of-grid triangles are on the limb, and if so, use the "nearest" method for both atm_extrapolation_method and ld_extrapolation_method regardless of what was set by the user.

[aprsa]

This might be a bigger discussion to have than it seems at first glance.

What would you expect to happen if atm parameters are off-grid, LD method is set to 'interp', and blending method is set to 'blackbody' or to 'none'? How exactly should limb darkening be handled?

The current implementation uses ld_extrapolation_method to determine LD coefficients, i.e. when the LD method is non-interp. This fix uses ld_extrapolation_method for the case when LD method is set to 'interp'. But that implies that we must always extrapolate intensities (specific and normal), the ratio of which is then the limb darkening correction. Perhaps we shouldn't allow 'interp' for off-grid extrapolation, but how do we prescribe that from the user's point of view? We would need yet another set of parameters to give us the fallback option for LDs, which sounds messy.

Forcing Inorm and Imu to 'nearest' is indeed another possibility, but that should happen only if extrapolation isn't 'none'. So perhaps something like:

ld_extrapolation_method = 'none' if atm_extrapolation_method == 'none' else 'nearest'

But then we're essentially ignoring what the user passes as ld_extrapolation_method.

None of these seem ideal to me.

[Raindogjones]

Maybe not ideal but I think the second option while throwing a warning makes most sense.

[mwrona77]

Sorry you have to explain this to me in detail, I cannot say I fully understand this method.

I see your point, but let me ask this question first: do we know what exactly causes the intensity to blow up? Is it division by something close to zero? In this case, it has to be inorms_atm in imus_bb = imus_atm / inorms_atm * inorms_bb. You said earlier that:

The default was linear extrapolation in LD space, which causes blow-up.

and I am just trying to understand why this is so drastic.

Now for the fix. Your patch only works when blending is on and when ld_extrapolation_method is "nearest". As soon as someone turns off blending and switches to extrapolation, the problem reappears. To make it work when the blending is off, I have to use:
b.set_value_all('atm_extrapolation_method', value="nearest").

So the question is, are we missing a change from atm_extrapolation_method to ld_extrapolation_method somewhere? You changed it in two places, maybe there is another one that has to be changed.

I like Dave's idea with a warning. Maybe when a user decides to use "linear" instead of "nearest", we can show a warning that it may cause unwanted extrapolation if the off-grid triangles are on the limb.

EDIT

elif atm.name in self.ndp.keys():  # atm in one of the model atmospheres
    result = InterpResult.from_ndpolator(
        self.ndp[atm.name].ndpolate(
            f'imu@{intens_weighting}',
            query_pts=query.subset(atm.basic_axis_names + ['mus']).pts,
            extrapolation_method=atm_extrapolation_method
        )
    )

I just changed extrapolation_method=atm_extrapolation_method to extrapolation_method=ld_extrapolation_method, and now it also works when the blending is off. Of course, only because ld_extrapolation_method is nearest.

[kecnry]

pretty sure this is in the repo (although I don't think we use it anymore and probably should clean it up).