sirocco-rt / sirocco

This is the repository for Sirocco, the radiative transfer code used to winds in AGN and other syatems
GNU General Public License v3.0
27 stars 24 forks source link

Various forms of sigma_phot #86

Closed kslong closed 9 years ago

kslong commented 10 years ago

We have three routines, sigma_phot, sigma_phot_topbase, and sigma_phot_verner in radiation.c which calculate x-sections given certain information. sigma_phot and sigma_phot_verner refer to two different verner datasets. Apparently the difference is that sigma_phot uses a version of the verner data that includes inner shell ionization, but sigma_phot_verner does not.

I would have thought that sigma_phot should get the photoionization x-section regardless of the type of data that is supplied, but that is not the case.

As a documentation issue, at least, we need to understand better the reason for the various versions of sigma_phot, and why they are used in different situations. It is possible that we rename sigma_phot to something else, and that we should create a new sigma_phot which can calculate the x-section in all situations. (note that there are separate routines to calcuate partial x-sections).

Comments & explanations are invited.

kslong commented 10 years ago

All of the places that various versions of sigma_phot are used

atomic.c: double sigma_phot(x_ptr,freq) calculates the photionization crossection due to the transition atomic.c: sigma_phot uses Verner et al.'s interpolation formulae for the photoionization crossection atomic.c:sigma_phot (x_ptr, freq) balance_abso.c: double sigma_phot (); balance_abso.c: x = (f - ft) * f * f * exp (fb_alpha * (f - ft)) * sigma_phot (fb_xptr, f); estimators.c: double sigma_phot_topbase (); estimators.c: x = sigma_phot_topbase (&phot_top[n], freq_av); //this is the cross section estimators.c: x = sigma_phot_verner (&augerion[n], freq_av); //this is the cross section estimators.c: x = sigma_phot_topbase (cont_ext_ptr2, freq); //this is the cross-section estimators.c: x = sigma_phot_topbase (cont_ext_ptr2, freq); //this is the cross-section estimators.c: x = sigma_phot_topbase (cont_ext_ptr2, freq); //this is the cross-section estimators.c: x = sigma_phot_topbase (cont_ext_ptr2, freq); //this is the cross-section get_atomicdata.c:The routine sigma_phot(xptr, freq) calculates the crossection based on this. get_atomicdata.c: double sigma_phot(); get_atomicdata.c: xphot_tab[j].x[n]=sigma_phot (xver, freq); matom.c: x = sigma_phot_topbase (cont_ext_ptr, freq); //this is the cross-section

Radiation is where the various versions of sigma_phot are defined radiation.c: sigma_phot_topbase (x_top_ptr, radiation.c: sigma_phot (x_ptr, freq_xs) * density * frac_path; radiation.c: x = sigma_phot_verner (&augerion[n], freq); //this is the cross section radiation.c: double sigma_phot(x_ptr,freq) calculates the photionization crossection due to the transition radiation.c: sigma_phot uses Verner et al.'s interpolation formulae for the photoionization crossection radiation.c:sigma_phot (x_ptr, freq) radiation.c: double sigma_phot_topbase(x_ptr,freq) calculates the radiation.c: sigma_phot uses the Topbase x-sections to calculate the bound free radiation.c:sigma_phot_topbase (x_ptr, freq) radiation.c: double sigma_phot_verner(x_ptr,freq) calculates the photionization crossection due to the transition radiation.c: Same as sigma_phot but using the older compitation from Verner that includes inner shells radiation.c: 08dec SS Coded (actually mostly copied from sigma_phot) radiation.c:sigma_phot_verner (x_ptr, freq)

recomb.c: double sigma_phot (); recomb.c: x = sigma_phot (fb_xver, freq); recomb.c: double sigma_phot_topbase (); recomb.c: x = sigma_phot_topbase (fb_xtop, freq); recomb.c: double sigma_phot (), sigma_phot_topbase (); resonate.c: kap_bf[nn] = x = sigma_phot_topbase (&phot_top[n], freq) * density; //stimulated recombination? (SS) variable_temperature.c: answer = sigma_phot_topbase (xtop, freq); variable_temperature.c: answer = sigma_phot (xver, freq); variable_temperature.c:// answer = sigma_phot_topbase (xtop, freq); // and finally multiply by the cross section. variable_temperature.c: answer = sigma_phot_topbase (xtop, freq); // and finally multiply by the cross section. variable_temperature.c: answer = sigma_phot (xver, freq); // and finally multiply by the cross section. variable_temperature.c: answer = sigma_phot (xver, freq); // and finally multiply by the cross section. variable_temperature.c: answer *= sigma_phot_topbase (xtop, freq); // and finally multiply by the cross section.

Higginbottom commented 10 years ago

A few months ago, mainly in order to improve performance, I wrote some code that creates a tabulated version of the verner cross sections. This is used in the integrations in the variable temperature subroutines, and that now only uses sigma_phot_topbase (the calls listed in the previous comment to sigma_phot_ver are actually commented out in the code).

This seems to suggest a way forward. If we are happy that we can always tabulate the verner data to capture the cross section properly (and I think we should be) then we could simply use the tabulated version in all cases.

So, during read_atomic data, we decide which is the best photoionization cross section data for that ion (presumably topbase first, then verner) tabulate it if necessary, and store it. We use the same array whatever the origin of the data. This might need a little care....

Then, wherenver we want to get a cross section, we would use the same subroutine - essentially sigma_phot_topbase, but probably renamed.

If we get some new cross section data in the future, all we need to do is write code in get_atomic_data to read it in, and tabulate it.

kslong commented 10 years ago

Is the suggestion just to tabulate all of these x-sections, whatever the source. If so, seems reasonable. I assume if we do this we would maintain helper routines to do the tabulations. Is there any issue regarding level descriptions?

jhmatthews commented 9 years ago

Potential solution here is probably to have one structure for photoionization cross sections which contains all the variables one might need for either topbase or VFKY cross sections, instead of a separate one for both. This is then passed to the governing sigma_phot routine.

This may take up a bit more space, but we can use dynamic allocation to address this. We would then need a 'type' variable in the structure in order to distinguish between the type of cross section read in.

If we tabulate the data, this is easier on space, because we simply have all data in the same format as our topbase files.

Higginbottom commented 9 years ago

Sounds good. The tabulate_verner routine at the end of get_atomic_data currently puts verner data into a separate structure called xphot_tab. If we only put topbase or verner into a general phot_tab structure, I think this would save space since currently we store both sets and then only use the VFKY tabulation in the absence of topbase. the get_pi_rates routine would require minimal modification to make use of a new structure, the main effort would be in the recomb routines. Should be conceptually pretty simple tho, one would just need to be careful....

kslong commented 9 years ago

No objections, though this has grown from my original issue, which was simply that we should have a single call to sigma-phot, so that the only place there were if statements were in the sigma-phot routine itself. Note that get-atomic-data takes care of the order in which different cross sections are chosen

jhmatthews commented 9 years ago

This has been fixed by #143 so can be closed. See also discussions in #141