Revamp some of the warm, thick target code

[settwi]

I think "legacy" is the only place where the warm, thick target code appears (described in Kontar+2015).

The legacy code is is basically just a direct port from IDL and isn't very descriptive:

sunkit-spex/sunkit_spex/legacy/fitting/photon_models_for_fitting.py

Lines 100 to 161 in 7f350be

	def thick_warm(total_eflux, ind, e_c, plasma_d, loop_temp, length, energies=None):
	"""Calculates the warm thick-target bremsstrahlung radiation as seen from Earth.
	[1] Kontar et al, ApJ 2015 (http://adsabs.harvard.edu/abs/2015arXiv150503733K)
	[2] https://hesperia.gsfc.nasa.gov/ssw/packages/xray/idl/f_thick_warm.pro
	Parameters
	----------
	energies : 2d array
	Array of energy bins for the model to be calculated over.
	E.g., [[1,1.5],[1.5,2],[2,2.5],...].
	total_eflux : int or float
	Total integrated electron flux in units of 10^35 e^- s^-1.
	ind : int or float
	Power-law index of the electron distribution.
	e_c : int or float
	Low-energy cut-off of the electron distribution in units of keV.
	plasma_d : int or float
	Plasma number density in units of 1e10 cm^-3.
	loop_temp : int or float
	Plasma temperature in megakelvin.
	length : int or float
	Plasma column length (usually half the full loop length?) in Mm.
	Returns
	-------
	A 1d array of the warm thick-target bremsstrahlung radiation in units
	of ph s^-1 cm^-2 keV^-1.
	"""
	ME_KEV = 511 # [keV]
	CC = 2.99e10 # speed of light [cm/s]
	KK = 2.6e-18
	n_p = plasma_d * 1e10 # was in 1e10 cm^-3, now in cm^-3
	tloop = loop_temp * 8.6173e-2 # was in MK, now in keV
	l = length * 1e8 # was in Mm, now in cm
	ll = tloop**2 / (2 * KK * n_p) # collisional stopping distance for electrons of Tloop energy
	emin = tloop * 3 * (5 * ll / l) ** 4
	if emin > 0.1:
	print(
	f"The loop_temp ({loop_temp}), plasma density ({plasma_d}), and loop length ({length}) make emin ({emin}) >0.1. Fixing emin to 0.1."
	)
	emin = 0.1
	lmin = e_c**2 / (2 * KK * n_p) / 3
	if lmin > l:
	print("Minimum length>length")
	em_add = 3 * np.pi / 2 / KK / CC * np.sqrt(ME_KEV / 8.0) * tloop**2 / np.sqrt(emin) * total_eflux * 1e35
	em46 = em_add * 1e-46 # get EM in units of 10^46 cm^(-3)
	return thick_fn(total_eflux, ind, e_c, energies=energies) + f_vth(loop_temp, em46, energies=energies)

I would suggest moving forward that we factor out some computations into some separate functions. The bulk of the computation relates to computing the E_min and effective emission measure parameters. The computations could use more comments to help describe what's actually going on, and use better variable names to help with code readability.

Here's a crack at that portion. These returned parameters can be used to compute the emission from the thermalized warm target electrons. Strictly speaking, the e_min isn't necessary outside of computing the emission measure, so we could drop that return value.

https://gist.github.com/settwi/b5a398ec3cb2b2a73303344af151870f