Troubles with drPopulate()

[Alexander]

Hello!

I'm looking into using ChiantiPy version 0.15.0 in my research.

I noticed that the drPopulate() method from ChiantiPy.core.ion doesn't work. I used FeXIV, FeXXIII and CrXXIII as ions.

I tried to fix the errors in the code in order to restore the functionality of the method. Below are a few pieces of code by which you can identify areas where changes have been made. Also at the end there will be a link to the file where the entire drPopulate() method is written with my edits.

After I made changes, the method seems to work for FeXIV and FeXXIII. However, for CrXXIII I get the following error:

Traceback (most recent call last):
  File "C:\Python_programms\CHIANTI_PROJECT\venv\lib\site-packages\IPython\core\interactiveshell.py", line 3460, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-4-2c91415f15e9>", line 1, in <module>
    cr23.drPopulate()
  File "C:\Python_programms\CHIANTI_PROJECT\venv\lib\site-packages\ChiantiPy\core\Ion.py", line 1689, in drPopulate
    l1 = rrlvl['lvl1'][itrans] - 1
IndexError: index 123 is out of bounds for axis 0 with size 123

As I was able to understand, the problem occurs when executing this piece of drPopulate():

                if self.Nrrlvl:
                    for itrans in range(self.Nrrlvl):
                        l1 = rrlvl['lvl1'][itrans] - 1
                        l2 = rrlvl['lvl2'][itrans] - 1
                        popmat[l2+ci, -1] += self.EDensity[itemp]*hPop[itemp, l1]*self.RrlvlRate['rate'][itrans, itemp]
                        popmat[-1, -1] -= self.EDensity[itemp]*hPop[itemp, l1]*self.RrlvlRate['rate'][itrans,itemp]
                        rrTot[itemp] += hPop[itemp, l1]*self.RrlvlRate['rate'][itrans, itemp]
#                    else:
#                        rrTot[itemp] = 0.

And, as I understand it, the root cause for it is the code section from ChiatiPy.core.Ion.setup(), where self.Nrrlvl is defined, which for CrXXIII turns out to be 127.
Relevant section of code:

        if os.path.isfile(rrlvlfile):
            self.Rrlvl = io.cireclvlRead(self.IonStr, filetype='rrlvl')
            self.Nrrlvl = len(self.Rrlvl['lvl1'])
            nRrlvl = max(self.Rrlvl['lvl2'])
            self.Nrrlvl = nRrlvl
            nlvlList.append(nRrlvl)
        else:
            self.Nrrlvl = 0
        #  psplups file may not exist

After that, I decided to contact you. Please tell me what you can say about the current situation? Perhaps there is some mistake on my part?

Code used to test drPopulate():

import ChiantiPy.core as ch
import numpy as np
import matplotlib.pyplot as plt

temperature = np.logspace(5.8,6.8,21)
fe14 = ch.ion('fe_14',temperature=temperature,eDensity=1.e+9,em=1.e+27)
fe14.drPopulate()

fe23 = ch.ion('fe_23',temperature=temperature,eDensity=1.e+9,em=1.e+27)
fe23.drPopulate()

cr23 = ch.ion('cr_23',temperature=temperature,eDensity=1.e+9,em=1.e+27)
cr23.drPopulate()

Link to the file with my edits: https://drive.google.com/file/d/11G5w_LuPJLY2I7NmiGurutIedhCcPa1u/view?usp=share_link

Excerpts with edits from the file I modified:

FIX1

        #  Changed because there was an error:
        # Traceback (most recent call last):
        #   File "C:\Python_programms\CHIANTI_PROJECT\venv\lib\site-packages\IPython\core\interactiveshell.py", line 3460, in run_code
        #     exec(code_obj, self.user_global_ns, self.user_ns)
        #   File "<ipython-input-6-1a981ab8ece7>", line 1, in <module>
        #     fe14.drPopulate()
        #   File "C:\Python_programms\CHIANTI_PROJECT\venv\lib\site-packages\ChiantiPy\core\Ion.py", line 1601, in drPopulate
        #     temp = temperature
        # NameError: name 'temperature' is not defined
        #
        # Before:
        # temp = temperature
        # After:
        temp = self.Temperature
        ntemp = temp.size
        dens = self.EDensity
        ndens = dens.size
#        cc = const.collision*self.EDensity
        #
        # (4 pi a0^2)^(3/2) = 6.6011e-24 (Badnell et al, 2003, A&A 406, 1151
        coef1 = 6.6011e-24*(const.hartree/(2.*const.boltzmann*self.Temperature))**1.5
        coef2 = (const.planck)**3/(2.*const.pi*const.emass*const.boltzmann*self.Temperature)**1.5

FIX2

            if self.Nrrlvl:
                rrlvl = self.Rrlvl
                if hasattr(self, 'RrlvlRate'):
                    rrlvlRate = self.RrlvlRate
                else:
                    self.rrlvlDescale('rrlvl')
                    rrlvlRate = self.RrlvlRate
            # Changed because there was an error:
            # Traceback (most recent call last):
            #   File "C:\Python_programms\CHIANTI_PROJECT\venv\lib\site-packages\IPython\core\interactiveshell.py", line 3460, in run_code
            #     exec(code_obj, self.user_global_ns, self.user_ns)
            #   File "<ipython-input-7-1a981ab8ece7>", line 1, in <module>
            #     fe14.drPopulate()
            #   File "C:\Python_programms\CHIANTI_PROJECT\venv\lib\site-packages\ChiantiPy\core\Ion.py", line 1714, in drPopulate
            #     if higher.RecombRate['rate'][itemp] > (rrTot[itemp] + drTot[itemp]):
            # UnboundLocalError: local variable 'higher' referenced before assignment
            #
            # Before:
            #if rec:
            #After:
            if 1:

                # get ionization rate of this current ion
                self.ionizRate()
                #  get the higher ionization stage and its recombination rates
                highers = util.zion2name(self.Z, self.Ion+1)
                higher = ion(highers, temperature=self.Temperature, eDensity=self.EDensity)

FIX3

# Changed because there were errors where temperature and eDensity variables were not known.
        # Before:
        # self.Population = {"temperature":temperature,"eDensity":eDensity,"population":pop, "protonDensity":protonDensity, "ci":ci, "rec":rec, 'popmat':popmat, 'method':'drPopulate'}
        # After:
        self.Population = {"temperature":self.Temperature,"eDensity": self.EDensity,"population":pop, "protonDensity":protonDensity, "ci":ci, "rec":rec, 'popmat':popmat, 'method':'drPopulate'}
        if len(errorMessage) > 0:
            self.Population['errorMessage'] = errorMessage

[Ken Dere]

do not use drPopulate - just use ion.populate()

drPopulate is no longer in the current version for good reasons

[Alexander]

Thank you very much for your answer!

I would like to clarify. Do I understand correctly that now the ChiantiPy.core.Ion.populate() method includes the dielectronic recombination from all levels specified by the .auto file?
Because I wanted to use the drPopulate() method precisely because of this possibility, which is indicated here: ChiantiPy.core package — ChiantiPy 0.15.0 documentation

четверг, 27 апреля 2023 г. в 00:18:38 UTC+3, Ken Dere:

[Ken Dere]

ion.populate includes everything, included dielectronic and radiative recombinations where the necessary files are available

if you have further questions, please ask

regards,

Ken

[Alexander]

Accepted. Thanks a lot!

четверг, 27 апреля 2023 г. в 19:55:02 UTC+3, Ken Dere: