matom line + bf luminosity bug

[jhmatthews]

Current Status: Unsolved.
Version discovered in: Python 76, although seemingly as far back as 68.
Brief statement of problem: As described on the macro atoms page, our line luminosities and bound free luminosities
How to reproduce: Macro atom yso.pf model
Discussion:

When we compare python 76 and python 58 of observations, we make a number of observations

• After one cycle, the temperatures are the same
• After one cycle, the line luminosities are different in py_wind, and printed on screen.
• This is partly due to concerns about how we report and save cooling/luminosity #36 - we were confused as to how an identical temperature could lead to a different bf luminosity ... this issue confused things somewhat- the heating and cooling that are actually used to minimise the temperature are the same
• After one cycle, the level populations are different, which is what causes the luminosities to be different after one cycle
• After more cycles, and once the model has converged, the temperature and ionization structure are different, and so is the spectrum accordingly.

This possibly suggests that it is something in the macro_pops routine.

The macro_pops routine should compute the fractional level populations for macro atoms and store them in "levden" array. The ion fractions are also computed and stored in w[n].density[nion]
This routine uses a matrix inversion method to get the level populations, the matrices being inverted are the rate matrices.

Initial debugging of this routine reveals that the rate_matrices are different, but I haven't nailed down where the actual problem occurs. Will resume on my return.

[jhmatthews]

A little progress here.

• We know from the above that something is causing macro_pops() to work out the level populations differently, which would lead to line and bf luminosities coming out differently as we observe
• First, that means understanding how the macro_pops routine works.
  • Broadly speaking, it calculates rates between different levels of a macro-atom and constructs a matrix which stores these rates
  • This matrix is then inverted to obtain the level populations via LU decomposition (a method which is not that important to understand for the purposes of this bug, because...
• Inspection of 58 and 76 reveals that the routines are close to identical, and it is the calculation of the rates which is different, rather than anything in the matrix method
• to calculate the rate of bound-bound transitions we use the normal formula

rate = (a21 (line_ptr) * p_escape (line_ptr, xplasma));

• if any of the rates are wrong in this routine that would cause different level populations, and even though the problem may only be in a bb calculation it could propagate through the bb and bf emission as everything is linked to the levels.
• in 58 and 76, the escape probabilities are different
• why? because the density of neutral hydrogen is different when macro pops is called- so the bb opacities are different in this calculation!
• when I plotted up, and even compared numbers, ne was very very similar (to one part in 10^7) in all cells. However, because this plasma is almost all ionised then this means that the n_h values (around 10^5) can be out by a factor of 2 while hardly affecting ne...
• And from what I can tell, the reason the densities are different is because of partition functions.
• the saha equation is called just before you call macro_pops, which alters the ion densities which are fed to macro pops
• this returns different values, and the reason it returns different values is because the partition function calculations are different.

So, I have a number of concluding points here.

• First, how to fix, and how to tell which is right. There are quite a lot of differences in code between the way partition functions are done in 58 and 76 (for example, partition.c doesn't exist in older versions). Inspection of the confluence site reveals a few notes on partition functions around python 62:

Python_62 - A ksl update to the code
0808 - 62
- Rewrote the partition function sections of the code to be more like calls to other similar routines in this section of the program
- Removed the forttran routine in the code, which contained Ivan's calculations of partition functions
- Need to fix up the partition function problem
- Fixed the partition problem; it was due to the difference between the indexing into the levden array and the nlte levels.
- Modified get_atomicdata so that it will only read in one type of configuration information for each ion ... i.e. macro-records, topbase records, etc.

Python_62 08Aug
- Fixes problems with the calculation of partiotion functions that arose because the the levden array in the plasma structure is not exactly parallel to the nlte levels
- Changes the behavior of get_atomicdata so that one has to use a single type of configuration for a single ion. One can no longer, for example, mix kurucz and topbase levels. Now if you allow say 5 nlte levels (by which one means that the level density is tracked in the plasma structure), the first five levels from the configuration file are nlte levels and the remainder for that ion are treated as simple ions.
080810 - Still need to verify that everything is handle properly in den_config

Haven't really been able to find any info on what the 'partition function problem' was, maybe someone can shed some light on that?

• Second, to understand what is going on here a little better.
  • Is calling saha before macro_pops everytime a good idea? It means that this escape probability calculation will always be done with saha ion densities..
  • Is the macro pops routine as well as we can do? It uses an ion density as an input and also alters the ion density itself. Further more, the ion densities it uses as inputs are not those calculated in the previous cycle...sort of similar to the concerns about how we report and save cooling/luminosity #36 bug in that there could a kind of iteration here as it's not entirely consistent, but perhaps it's good enough.
• Third, is it the only problem...experience tells me probably not...

Still, getting there, I think.

[jhmatthews]

This bug is now mainly discussed in #40

[jhmatthews]

I'm closing this issue with the caveat that issue #40 is still open and we don't fully understand the convergence on macro populations yet, but our line strengths now agree pretty well. Bug 40 will be the place for future discussion of this issue.