wtbarnes / fiasco

Python interface to the CHIANTI atomic database
http://fiasco.readthedocs.io/
BSD 3-Clause "New" or "Revised" License
21 stars 17 forks source link

Add calculation for two-photon continuum #43

Closed wtbarnes closed 6 months ago

jwreep commented 9 months ago

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

wtbarnes commented 9 months ago

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.

wtbarnes commented 9 months ago

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.

wtbarnes commented 9 months ago

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

jwreep commented 9 months ago

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?

wtbarnes commented 9 months ago

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.

jwreep commented 9 months ago

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.

wtbarnes commented 9 months ago

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?

jwreep commented 9 months ago

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

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]

wtbarnes commented 9 months ago

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.

jwreep commented 9 months ago

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

wtbarnes commented 8 months ago

Cross referencing chianti-atomic/ChiantiPy#458 for completeness

jwreep commented 8 months ago

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.

wtbarnes commented 8 months ago

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!

pryoung commented 8 months ago

Thanks for sorting this out – I made a start on a technical report for the two-photon continuum, so I’ll give this a mention. Peter.

From: Will Barnes @.> Reply-To: wtbarnes/fiasco @.> Date: Thursday, February 15, 2024 at 12:11 AM To: wtbarnes/fiasco @.> Cc: Peter Young @.>, Mention @.***> Subject: [EXTERNAL] [BULK] Re: [wtbarnes/fiasco] Add calculation for two-photon continuum (#43)

CAUTION: This email originated from outside of NASA. Please take care when clicking links or opening attachments. Use the "Report Message" button to report suspicious messages to the NASA SOC.

Oh great! I don't know how I missed that comment.

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!

— Reply to this email directly, view it on GitHubhttps://github.com/wtbarnes/fiasco/issues/43#issuecomment-1945382321, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEL4UCGJ2WEXSGBIJGCS3H3YTWKIBAVCNFSM4GQQRLN2U5DIOJSWCZC7NNSXTN2JONZXKZKDN5WW2ZLOOQ5TCOJUGUZTQMRTGIYQ. You are receiving this because you were mentioned.Message ID: @.***>