problem with free-bound calculation about chiantipy HOT 13 CLOSED

yes, that is what I meant, I was speaking in svn

from chiantipy.

What file and line number are being referenced specifically? is this in the Continuum module or spectrum module?

from chiantipy.

I just ran the code again, just as the user supplied, and got no errors. I may have done a pull from my origin master in between the two runs.

The problem had been in the Spectrum module after line 170

from chiantipy.

I just checked and I was using the v0.7.0 release so th bug may still be there in Spectrum after line 170

from chiantipy.

went back to the trunk version (~0.7.1) and the problem still seemed to be gone

from chiantipy.

When you say v0.7.1, do you mean the master branch at the most current commit?

from chiantipy.

I ran the code snippet from the mailing list on the most current commit of master and it didn't error, but I didn't think that was the issue. I thought the problem was that the free-bound contribution to spectrum gives zero over that particular wavelength range and that this result is not correct.

from chiantipy.

yes, your are right. I just forgot

from chiantipy.

Just to consolidate the issue to this thread, the code to calculate the spectrum in ChiantiPy is,

temp=[1.e+7,2.e+7]
dens=1.e+9
wvl = 1. + 0.002*np.arange(4501)
emeas=[1.e+27,1.e+27]
s3 = ch.spectrum(temperature=temp,eDensity=dens,wavelength=wvl,filter = (chfilters.gaussian,.015),doContinuum=1, em=emeas,minAbund=1.e-5,verbose=0)
plt.plot(wvl,s3.FreeFree['intensity'][0])
plt.plot(wvl,s3.FreeBound['intensity'][0])

while the corresponding IDL output gives,

(solid is free-free, dashed is free-bound, dot-dashed is two-photon)

The next step I guess would just be to make sure that the Continuum module alone is giving the expected output, e.g. something like

mport numpy as np
import ChiantiPy.core as ch
import ChiantiPy.tools.data as chdata

temperature = np.array([1e7,2e7])
wavelength = np.linspace(1,10,4501)
total_fb_emission = np.zeros((temperature.shape[0],wavelength.shape[0]))

for ion in chdata.MasterList:
    tmp = ch.Continuum(ion,temperature)
    tmp.calculate_free_bound_emission(wavelength)
    total_fb_emission += tmp.free_bound_emission

If it does, then we know the problem is confined to spectrum

from chiantipy.

the problem is in Continuum, at line 449, in inserted

        if self.emission_measure is not None:
            for i, em in enumerate(self.emission_measure):
                fb_emiss[i] *= em

for this problem, fb has a shape of [2,4501] and emeas has a shape of [2] and they can not be directly multiplied together. I have put a correction and have pushed to to my fork and then will do it to the upstream master

also, line 421: requires an array of temperatures, not a float - it tries to get len(temperature)
I will let you work on this because the code with regard to this has changed from what I had
in v

from chiantipy.

Ok yes that's a problem. There are a couple of other corner cases I've hit where things cannot be multiplied together because of weird dimension mismatches. Note that this can be solved with numpy.newaxis. See my comment in #139. As they come up, I'll try to fix them. I think most of the time it is just a matter of using the numpy.newaxis fix. This is probably a case where more tests with a larger range of inputs would be useful...

Regarding line 421, I don't think this is an issue because of this line in the constructor which makes sure temperature is a Numpy array so calling len() on it shouldn't be a problem.

from chiantipy.

there is another problem with the fix and that is when emission_measure is not interable.
never used numpy.newaxis
I think I should leave this to you

from chiantipy.

the free-bound calculation now works with a single temperature. However, the syntax has mostly reverted to the original style

`import numpy as np
import ChiantiPy.core as ch
import ChiantiPy.tools.data as chdata

temperature = np.array([1e7,2e7])
wavelength = np.linspace(1,10,4501)
total_fb_emission = np.zeros((temperature.shape[0],wavelength.shape[0]))

for ion in chdata.MasterList:
tmp = ch.continuum(ion,temperature)
tmp.freeBound(wavelength)
total_fb_emission += tmp.FreeBound['intensity']`

previously, freeBound would through an error if the data file necessary to do the calculation was missing. Now it just returns a correctly shaped array of zeros.

from chiantipy.