Add calculation for two-photon continuum about fiasco HOT 15 CLOSED

The Chianti routine has lots of options in terms of inputs, outputs, and keywords. What input/output (variables, units, etc.) would be preferred here?

from fiasco.

Good question. Looking at the list of input parameters, I would say most are not needed due to the design of fiasco, e.g. sngl_ion or element do not need to be specified because this will be a method on an Ion class.

In an effort to keep the API as similar to that of the free-free and free-bound methods, I would say we should aim to just have wavelength and density be the inputs.

We could consider implementing a minimum abundance optional input too but that would be applied at the level of the method on IonCollection and if we go that route, I'd rather implement that on all of the continuum methods at once.

from fiasco.

Semi unrelated but the only reference I've found w.r.t. how this is implemented in CHIANTI is Eq 4 of this paper: http://adsabs.harvard.edu/abs/1999A%26AS..135..339L.

from fiasco.

Actually, equation 11 of this paper and the accompanying text are much more helpful: http://adsabs.harvard.edu/abs/2003ApJS..144..135Y

from fiasco.

The hydrogenic calculation uses both A_sum and A in the IDL routine, but this isn't mentioned in the Equations (or the Chianti user's guide):

y=wvl0/wvl[w_ind]

               y2=spl_init(y0,psi0[*,iz-1])
               distr=y*spl_interp(y0,psi0[*,iz-1],y2,y)/asum[iz-1]/wvl[w_ind]

               FOR i=0,nt-1 DO BEGIN
                  IF lookup_swtch EQ 0 THEN BEGIN
                     pop_solver, input, temps[i],edens[t_ind[i]],pop
                  ENDIF ELSE BEGIN
                     pop=reform(popx[i,i,*])   ; popx from lookup table.
                  ENDELSE 

                  IF keyword_set(photons) THEN BEGIN
                     distr1=rescale/4d0/!pi*avalue[iz-1]*this_abund* $
                            distr * $
                            (ioneq1[i]*pop[pop_idx]/edens[t_ind[i]]) * dem_arr[t_ind[i]]
                  ENDIF ELSE BEGIN
                     distr1=rescale*factor*1d8*avalue[iz-1]*this_abund* $
                            (distr/wvl[w_ind]) * $
                            (ioneq1[i]*pop[pop_idx]/edens[t_ind[i]]) * dem_arr[t_ind[i]]
                  ENDELSE
                  h_rad[w_ind,t_ind[i]]=h_rad[w_ind,t_ind[i]]+distr1

I'm assuming this is folded into the definition of the spectral distribution function, but I'm not sure I know what that function is. Do you have any idea?

from fiasco.

I think all of this information is available in the hseq_2photon.dat and heseq_2photon.dat files. Conveniently, I already wrote a parser for these files and so the information is available on each ion already. Have a look at the _hseq and heseq attributes on the ion object. Just taking a cursory look at the IDL code you pasted above I think those should have everything you need.

from fiasco.

Oh, for sure I can translate the code. I'm just wondering what it actually is, and what these equations actually are. There's no mention of A_sum in any of the Chianti papers I've seen, for example.

from fiasco.

It's certainly quite mysterious. The only thing I can find is this note in the IDL docstring:

; PROGRAMMING NOTES
;
;    For He 2-photon decays, the distribution function is from Table II 
;    of Drake et al. (1969), except that the values have been divided by 
;    the A-value from Drake (1986)

This maybe explains why the spectral distribution function for the H-like ions is normalized by A_sum but the He-like ions are not. Why it's done at all, I do not know.

@pryoung do you know why these A_sum terms appear in the H spectral distribution function?

from fiasco.

I can't access the Drake papers, but this one has a simple definition of the spectral distribution function:

No mention of A_sum though, or at least from what I'm understanding . . . but perhaps what you say is correct about its origin.

One more inquiry: should A_sum be unitless? It's currently given by the parser as 1/s, but I'm not sure that's correct from simple unit analysis. The hydrogenic calculation in Chianti has a factor of 1/A_sum which helium-like doesn't, and the units of everything else are otherwise identical.

H: distr=y*spl_interp(y0,psi0[*,iz-1],y2,y)/asum[iz-1]/wvl[w_ind]
He: distr=y*spl_interp(y0,psi0[*,iz-1],y2,y)/wvl[w_ind]

from fiasco.

Yeah, something about the units is not right. I think when I wrote the parser I assumed the units of A_sum were 1/s because those are the units of A but if that's the case, then psi can't be unitless for the H files.

We should just elevate this to the CHIANTI mailing list as I'm sure either Ken or Peter knows. Maybe some normalization was applied when the files were prepared.

from fiasco.

I'll send an e-mail. Not having much luck figuring it out otherwise.

from fiasco.

Cross referencing chianti-atomic/ChiantiPy#458 for completeness

from fiasco.

Think I sort of figured out $A_{\text{sum}}$: it's a normalization factor of the integral of $\psi(y)$. That is,

$\frac{1}{A_{\text{sum}}} \int_{0}^{1} \psi(y) dy = 2$

The IDL code has the following comments:

;                      psi for H- and He-like are normalized differently
;   psi is properly normalized to give an integral of 2.0

which seems to check out for any element. This would indeed be unitless, then.

In any case, I think we should come up with a better name. It has nothing to do with the $A_{ji}$ Einstein coefficients.

from fiasco.

Oh great! I don't know how I missed that comment. I suppose me naively assuming A_ij and A_sum are related and thus must have the same units is what led to all of this confusion!

I completely agree re: a new name. Maybe something like psi_norm? And then in the description field we can give more details. Let's also try to preserve some of this discussion in a comment somewhere in the code since it doesn't seem to appear anywhere else in the CHIANTI papers, software, or docs!

from fiasco.

from fiasco.