389 implement split atmosphere interpolation for 4 dimensional atmospheres

[Basdorsman]

I hope doing this in a pull request makes sense, I thought of writing a status update in the issue itself but there you cannot easily see my code updates as neatly.

What is the goal?

The goal is to implement the split interpolation methodology for 4D (NSX) atmospheres, similar to how it is already implemented for 5D atmospheres.

What is the status?

The implementation is now finished, but there appears to be a small offset between the interpolation methods in split vs combined 4D that I don't see in split vs combined 5D. The small offset is -68144 vs -68147 in ln-space when I run TestRun_Num_test.py using split vs combined methods respectively. The difference in pulse profile is not visible by eye. Selecting one energy-phase bin the difference in intensity is 0.007%. Small, but I'm still suspicious because I don't see any difference in 5D. Regardless of the small difference in the result, it runs faster. in the same testrun, I see that the photosphere.integrate() is almost twice as fast. The potential gain in likelihood speedup is not as impressive as there are other parts in the code that contribute heavily too, such as the ray construction(?).

I have also moved all unit conversions to hot_wrapper.pyx, taking place around the relevant eval_hot_I functions. My idea here is that I think it is cleaner if all interpolations are cleared of unit conversions, which makes the interpolators unit agnostic and easily testable with a generic case.

Possible reasons for difference in pulse

I am a bit at a loss about what is causing the difference, but here is one idea.

• Changing the inpolation from one step to two exposes the process to different rounding errors (or otherwise numerical artefacts). I wonder if perhaps due to unit conversion (there is a unit conversion in the NSX atmosphere that involves a log) could enhance such numerical artefacts and introduce a slight error.

How could that difference be addressed

I would test if the split interpolator is healthy by comparing results from the combined interplator preferably in the test run that I have used so far. A broader parameter sweep could later identify if it is somehow parameter dependent.

To dos related to clean-up

• I think hot_num_4D.pyx can already be safely removed, to be superseded by hot_num_4D_split.pyx. Maybe the name can be shortened.
• Same for hot_num_2D.pyx, but in this case some edits are left in the hot_wrapper.pyx
• Commented lines that serve no function can be removed.
• Docstrings.

[thjsal]

Nice work! Need to look more into that. But I realized that the CI tests probably fail because you named the new example so that it ends "_test.py" and therefore GitHub Actions try to execute it (and it fails since it does not find the required input files).

[thjsal]

@Basdorsman: I did a bit testing and realized why the 2+2 split integration does not exactly work for the NSX model. The reason is that one of the parameters is surface gravity (log g), which actually is not constant within the spot. So when you apply the usual 4D interpolation, you have different surface gravity values for each co-latitude (since because of the oblate shape of the star, each co-latitude corresponds to different distance from the star center). But when you apply the split interpolation technique, you assume that log g is constant within the whole spot. This approximation is worse the bigger the spot is or the more oblate the star is.

So if you don't want to loose any accuracy in the current 4D calculation, you could only do 1+3 split interpolation. But I don't know if that would have a significant speed benefit, and is it worth to do.

Or if loosing this bit in the accuracy is ok (to get the higher speed benefit), then maybe even this 2+2 method could be useful (I used this same simplification in Salmi+2019 when simulating RMP data, discussed in Section 2 there).