johnson-c
Added capability to model energy evolution of a species in a plasma
This adds the capability to track the time evolution of the temperature of a species during the course of an ionization balance. It runs the ionization balance as normal, and in addition tracks the change in total energy of each charge state over time by evaluating the collisional heating between each charge state and background ions.
The rate equations for energy transfer are non-linear, so the matrix exponentiation method cannot be used for solving the system. Instead the scipy.integrate.solve_ivp method is used, which numerically integrates the system of ODE's.
The temperature of each charge state is calculated using the equation E = 3/2 * N * T, where E is the energy in each charge state, N is the population of each charge state, and T is the temperature of each charge state. During times where the population is small, the division of two very small numbers can cause T to become numerically unstable.
This modelling was motivated by trying to model the temperature of carbon charge states sputtered off a material into a background plasma (with the goal of understanding how the temperature of each charge state measured by a spectrometer relates to the background ion temperature). Therefore, it is not tested for atoms other than carbon.
@johnson-c I mostly fixed the numerical instability issue with the solver. For the vast majority of cases, the energy balance runs in less than a second, but there are still an extremely small number of cases (below 0.01% of the total over the density/temperature grid I'm using) where the solver becomes numerically unstable and takes hours to run (although it still seems to get the right answer in the end). I'm not sure why this happens but I think in the end it is fine since it happens so rarely.
Looking at the results, such as in the plot below, I notice inflection points in the curves at 50 eV and 70 eV. These seem to correlate with the Te grid points in the acd/scd adf11 files. Is it possible something weird is going on with the interpolation of the rate coefficients from the acd/scd files? This is not really a problem, but if there is an easy fix (maybe doing log interpolation instead of linear interpolation?) then I would be interested.
From my side, I think this PR is ready to merge.
@johnson-c I added the example file. And I also had to fix a small bug in read_adf11 that was causing it to fail for windows-style paths of adf11 files if the path contained the character 'r' (which made the file get interpreted as metastable resolved even if it wasn't).