Add albedo correction

[samaloney]

PR Description

A a function to read in precomputed greens table and return spectrum with albedo correction, would close #128

• Split function at the moment every time the function is called we would make web request and load a file from disk both very slow not ideal for fitting. The idea would be split up into at least two parts depending on the use case/s. As far as I know often $\theta$ is fix so the loading and interpolation between $\theta$ can be cached. So one function to load the appropriate greens matrix file and interpolate for given angle. A second function would call the previous function and do the energy intepotalation (as this could change in a fit e.g. gain correction) and fold the given spectrum through obtained green matrix.
• Turn in to astropy model e.g. spec | albedo (still waiting on proper units support in astroy for this work Should the model composition operator, |, support units changing astropy/astropy#17302)
• Docs and tests

Recreated Fig 2 upper from the paper

[edkontar]

I have found IDL code used in 2005 to generate Green Matrices for OSPEX and put it on github:

https://github.com/edkontar/albedo

[natsat0919]

I have added the changes we made to the albedo model during the workshop. For now we made the model "unitless" until we resolve the astropy model unit problems with SRM. All the units-related lines are commented so it can be easily updated.

[samaloney]

One of the main reason to separate things is to make them testable with these change it will be harder to test so I wonder are they necessary?

[samaloney]

Why move onto the object?

[samaloney]

Don't think this is needed as already cached by the decorators

[jajmitchell]

Perhaps we can remove the decorator then? referencing in the init step seems to be more inline with the Astropy API. From looking over the Astropy API all model classes have their reliant functions handled internally, the previous method was encountering issues with fitting theta as a free parameter, I think at least in part because of the functions being defined outside of the class, although I could be wrong about that. Either way, this version solves that issue and fully functions with the Astropy fitters, and seemingly produces correct results, also it is more in line with the Astropy model class API. I guess it comes back to the main question of to what degree we want to align ourselves with the API

[edkontar]

In OSPEX/drm_correct_albedo.pro, the anisotropy parameter is “a coefficient showing the ratio of the flux in observer direction to the flux downwards." A smaller value means more of an effect from albedo. A value of 1000 means very little effect. if anisotropy=1 (default) the source is isotropic. However, in the paper, the coefficient α that accounts for anisotropy of the source in the Eddington approximation – that is, "distinct but constant specific intensities in the downward and in the upward hemispheres, in the ratio α” i.e. the anisotropy described as the reciprocal of the anisotropy parameter in the OSPEX.

Should the description state the definition of "α" was assumed as in OSPEX?

class Albedo(FittableModel):
r"""
Aldedo model which adds albdeo correction to input spectrum.
Following [Kontar2006]_ using precomputed green matrices distributed as part of [SSW]_.
.. [Kontar2006] https://doi.org/10.1051/0004-6361:20053672
.. [SSW] https://www.lmsal.com/solarsoft/
Parameters
==========
energy_edges :
Energy edges associated with input spectrum
theta :
Angle between Sun-observer line and X-ray source
anisotropy :
Ratio of the flux in observer direction to the flux downwards, 1 for an isotropic source

[jajmitchell]

I have reverted the code back to its state before the previous commit with some small edits, the functions are now defined outside of the class and the class is used a wrapper, therefore the tests should now work, is it ok for me to push these changes? The albedo model is now functional with astropy fitters.

[samaloney]

Go for it, sounds good

[samaloney]

Suggested change

	if hasattr(theta, "unit"):
	theta = Quantity(theta.value,theta.unit)
	else:
	theta = theta*u.deg
	if not isinstance(theta, Quantity)
	theta = theta*u.deg

[settwi]

hey it looks like y'all are making good progress on this, i am gonna put a couple style suggestions along these lines if that's ok. feel free to reject them if you want but yea