Some brainstorming about ideas for a package for computing differential emission measures using a unified API and data container

dem_package_ideas.md

sunkit-dem

A package for computing differential emission measures using multiple methods with a common API

DEM inversion methods

The table below is a (probably not comprehensive) list of the different methods for computing the differential emission measure.

Name	Reference	Notes	Source code
MCMC	Vinay & Kashyap (1998)	Easily leverage many Python MCMC implementations	PINTofALE
demreg	Hannah & Kontar (2012)		https://github.com/ianan/demreg
sparse AIA	Cheung et al. (2015)		http://www.lmsal.com/~cheung/AIA/tutorial_dem/
sparse Bayesian	Warren et al. (2017)
xrt_dem_iterative2.pro		XRT DEM model. Does not seem to have a reference in literature?	SSW
Fast DEM Inversion	Plowman et al. (2012)		http://solar.physics.montana.edu/plowman/firdems.tgz
Pottasch method	Pottasch (1963)
EM loci method	see Veck et al. (1984)	Find loci of curves defined by ratio of intensity to contribution function; these are upper limits for DEM
A. Leonard Python method	Leonard and Morgan (2014)	In Python	https://github.com/SolarDrew/CoronaTemps

Also see Aschwanden et al. (2015) for a comparison of many of these methods as well as section 7 of Del Zanna and Mason (2018) for a brief overview of some methods.

Implementation

The idea would be that the package provides the infrastructure and the different methods are implemented as plugins. This could be, e.g. as subclasses of a base method class. Some would be provided by the package, but any arbitrary method could be used by a user.

Initially, we would need to develop the infrastructure to allow for this plugin architecture. Ideally, both the input and output would be an NDCube object.

A rough outline of the steps one would go through to calculate the DEM is as follows:

1. Input a multiwavelength observation. This could be a single pixel, an average over many pixels, or a a grid of pixels (0, 1, or 2 dimensions). The 3rd (or 2nd or 1st) dimension would be wavelength. This could also be a channel in the case of a broadband imager, though I don't think that distinction need be made.
2. Input the response functions or contributions for each channel or wavelength. Again this distinction is not that important. The only requirement is that the number of response (contribution) functions is equal to the number of wavelengths (channels).
3. Run the inversion method. This is the plugin part. To select the method, this could either be a string when creating the model instance (for a builtin method) or by providing a class or a function. How does astropy handle plugins? This could be a good guide...
4. Return an NDCube object, this time with the 3rd axis as temperature.