Improve reporting of macro bound-free heating

[jhmatthews]

Currently, we include collisional bound-free processes in macro-atoms as ncontribution to the bound-free heating and cooling (i.e. colisional recombination is calculated as "bound-free heating" and reported as "photo" heating). This change preserves this functionality but separates out the calculation into separate routines so that improved diagnostics are accessible.

This allows the user to see when macro-atom collisional processes are dominating the bound-free heating and cooling. The printout goes from looking like this:

!!wind_update: Absorbed flux    0.00e+00  (photo 0.00e+00 ff 2.75e+41 compton 7.70e+41 auger 0.00e+00 induced_compton 0.00e+00 lines 2.94e+41)
!!wind_update: Wind heating     4.09e+42  (photo 2.75e+42 ff 2.75e+41 compton 7.70e+41 auger 0.00e+00 induced_compton 0.00e+00 lines 2.94e+41 adiabatic 0.00e+00)
!!wind_update: Wind luminosity  5.42e+46 (recomb 6.24e+44 ff 2.73e+44 lines 5.33e+46) after update
init_freebound: Creating recombination emissivities between 1.000000e+14 and 1.000000e+20 in structure element 1
!!wind_update: Rad luminosity  5.42e+46 (recomb 6.24e+44 ff 2.56e+44 lines 5.33e+46) after update
!!wind_update: Wind cooling     5.41e+46 (recomb 4.20e+44 ff 2.56e+44 compton 2.42e+41 DR 0.00e+00 DI 5.73e+43 lines 5.33e+46 adiabatic 3.07e+40) after update

to this:

!!wind_update: Absorbed flux    0.00e+00  (photo 0.00e+00 ff 2.75e+41 compton 7.70e+41 auger 0.00e+00 induced_compton 0.00e+00 lines 2.94e+41)
!!wind_update: Wind heating     4.09e+42  (photo 2.75e+42 ff 2.75e+41 compton 7.70e+41 auger 0.00e+00 induced_compton 0.00e+00 lines 2.94e+41 adiabatic 0.00e+00)
!!wind_update: Wind luminosity  5.42e+46 (recomb 6.24e+44 ff 2.73e+44 lines 5.33e+46) after update
init_freebound: Creating recombination emissivities between 1.000000e+14 and 1.000000e+20 in structure element 1
!!wind_update: Rad luminosity  5.42e+46 (recomb 6.24e+44 ff 2.56e+44 lines 5.33e+46) after update
!!wind_update: Wind cooling     5.41e+46 (recomb 4.20e+44 ff 2.56e+44 compton 2.42e+41 DR 0.00e+00 DI 5.73e+43 lines 5.33e+46 adiabatic 3.07e+40) after update
!!wind_update: macro-atom heating: photoionization 4.53e+40 three body recomb 5.17e+41 lines 2.94e+41
!!wind_update: macro-atom cooling: photoionization 3.42e+44 collisional ionization 4.77e+43 lines 5.19e+44
!!wind_update: macro-atom energy flow: in 1.18e+41 out 1.78e+41

Corresponding columns have been added to the heat.txt files created by windsave2table which allows inspection of the quantities on a cell-by-cell basis. This is especially useful for high density issues, see #1152

[jhmatthews]

Here's the current regression test result from this. My confusion is that I would expect no changes to the spectra, but there are small "noise" level changes I am confused by.

summary.pdf

[kslong]

The noise differences in the regression tests imply that the random number secquence has changed and don't seem excessive to me. Claude identified sthe following:

Looking at this carefully, the split of the heating/cooling terms looks mathematically equivalent so it shouldn't change physics. But there
  are two things that stand out:

  1. trans_phot.c — this is the most suspicious for changing the random number sequence:
  Previously, BACK_RAD_SPECULAR photons hitting the star fell through the if block entirely — their direction was not changed. Now they enter
  the block and flip lmn[2]. If any regression test model uses BACK_RAD_SPECULAR, those photons now travel a completely different path,
  consuming random numbers differently from that point on. That would produce exactly the kind of noise-like differences you're seeing.

  2. Floating point ordering in normalise_macro_estimators:
  The old code incremented heat_tot twice in sequence:
  heat_tot += macro_bb_heating(...)
  heat_tot += macro_bf_heating(...)
  The new code does it in one combined increment:
  heat_tot += heat_lines_macro + heat_photo_macro + heat_qrecomb_macro
  The sum is algebraically identical but floating point addition is not associative — tiny rounding differences propagate through the
  temperature iteration, changing the converged temperature slightly, which then changes photon paths and the random number sequence in
  subsequent cycles.
  
  The trans_phot.c change is the more likely culprit if any regression test uses specular reflection. Do you know what BACK_RAD_SPECULAR
  corresponds to in the .pf files, and whether any regression test uses it?

Number 2 seems the likely culprit to me

[jhmatthews]

Ok, I'll quickly test that.

[jhmatthews]

Ok, thanks Knox, this was the cause. I checked that taking Sirocco 1.3 and just changing the order, which nicely illustrates how a very minimal change can result in noise-level differences in spectra. Merging this PR now.