Erica's Task List

TIP Last Update: JeffMangum - 2021-05-04

NOTE: Completed tasks shown as strike through

2020-10-15: Analysis path worked example

(HCN/HNC reduction and analysis).

2020-11-01: Heating model review

(see Meijerink, Bayet, Kazandjian, and de Mijolla articles above). When complete should have a good general understanding of how physical+chemical models of starburst galaxies are designed.

2021-05-01: Heating source investigation

This will involve looking for correlations between several potential heating sources (radio continuum point sources, vibrational spectral line point sources) and the NGC253 GMCs.
  • Radio Continuum: Review Ulvestad and Antonucci (1997) VLA imaging of NGC253. Use spectral index to identify supernovae/remnants and HII regions.
  • Vibrationally-Excited Spectral Line: Start with vib-HC3N positions and intensities from Rico-Villas et al.(2020). Theory here is that the high kinetic temperatures within the NGC 253 GMCs could be driven by the same physics that produces vib-HC3N rather than by mechanical heating.

2021-05-01: Complete General Integrated Intensity Image Ratio Analysis Routine

Recall that we have largely been studying integrated spectral line intensities and their ratios toward just the GMC positions in NGC253 for the HCN/HNC and H3Op/SO projects. In order to derive a more general analysis of the integrated intensity ratios within the ALCHEMI archive it would be useful to have a tool (python scripts) that can take two integrated intensity images, calculate a ratio image using proper signal-to-noise statistics, and extract integrated intensity line ratio values from specific positions or areas in the integrated intensity ratio image. For example, one might like to know what the integrated intensity ratio is along a NE-SW line along the CMZ where many of the radio continuum sources are found. Or perhaps one would also want to know what the integrated intensity image ratio is toward just the radio continuum sources that are identified as possible supernova remnants. It would be most useful if this tool (likely one or more python scripts) was written in a general way so as to allow anyone in the ALCHEMI collaboration to use it. This will likely require extensive use of astropy and associated python libraries. You can also leverage existing tools (CubeLineMoment: See the "Spectral Line Extraction" section of the ALCHEMIBestPractices page) and a FITS image ratio script that I have. The tricky part will be the generalizing feeding of positions or areas to extract ratio information and how to collect that information for later processing.

One of the first applications of this integrated intensity ratio investigative tool will be to extend our HCN/HNC and H3Op/SO to go beyond just the NGC253 GMCs.

2021-07-01: Generalized pyradex Modeling Package

There are two possible approaches to this task:

  • This first approach will involve starting with an adaptation of Adam Ginsburg's pyradex modeling construct applied to HCN and H13CN. The process that I went through to develop this modeling process is described in my Evernote for this task, which contains quite a bit of useful reference information on pyradex and model fitting. Leveraging the first-step generalization that I developed, Adam and I have spelled-out a path forward for developing a general RADEX model fitting routine. Note that subsequent discussions honed this "path forward" a bit, and clarified how the linewidth parameter is used in RADEX. I will provide the entire discussion thread once you start working on this task:
Hi Adam,

After working with your "adaptable" pyradex_{mol} (to create RADEX model cubes) and constrain_parameters_{mol}.py (to constrain the resultant LVG model (tr, tau, tex) cubes) to model H2CO, HCN, and H13CN, I think that I can see a path forward to a rather generalized set of model scripts that Joe Astronomer would be able to pick-up and use with little to no modification.  It assumes that no matter if Joe Astronomer uses taline or ratio to constrain the LVG models, the actual chi2 minimization is done on the ratio (i.e. if taline is entered, I calculate ratios).  The way I have done it at the moment, though, begs a few questions on exactly what that path is:

    The integrated intensities I used (extracted from CubeLineMoment mom0 images) are in units of Jy/beam*km/s, with subsequent conversion to TR*km/s before using them in constrain_parameters_{mol}.py.  Since I fitting to line ratios I think that these are equivalent to line temperatures, where of course I am making the simplification that all line widths are the same.  Since RADEX takes an input linewidth parameter (I have set it to 10 km/s), is it correct that this linewidth parameter is used to turn the modeled column densities per unit line width to modeled column densities for the assumed line width?
    Option 1 above may be too simplistic in that all linewidths are not really equal for all lines modeled.  Maybe I should try using mom1 from CubeLineMoment extracted as part of  [NOTE: While looking at one of the mom1 images, for NGC253_H13CN_21, there seemed to be some negative values near Region 3…  I thought that mom1 could not be negative?].  Assuming that we believe that the mom1 images are reliable (I have not done much with them), it would be easy to modify my script to extract the mom1 values associated with each input region.

To further explain how I have done the fitting with ratios, consider a linear molecule.  If you enter taline values (and their uncertainties) for the 1-0, 2-1, and 3-2, I calculate 2-1/1-0 and 3-2/2-1 and send those to the chi2 fitter to derive best-fits to (n,N,T).  One thing I should mention.  Since one is fitting to ratios, is the number of degrees of freedom N-3 or N-3-1, where N is the number of taline values entered?  Using the example I just gave, I gave it 3 taline values, but only 2 ratios were used to do the chi2 constraining.

Note that this will not meet every LVG modeling need, but I think that the ratio-chi2-fitting approach is about as general as one can get.

-- Jeff
  • In this second approach, Jon Holdship has developed a very similar package called SpectralRadex. Here is Jon's description of what he has created and why he created it:
Hi Jeff,

I'm glad you emailed! I've spent far too much time on this so I'm happy someone else might be interest in a similar tool. I've tried a litany of python-radex implementations over the past 6 years and ndradex was just the latest version I had found. Nothing has ever satisfied me completely.

    pyradex cannot be run in parallel (you can't have more than one copy of that class according to the developer)
    ndradex's parallelization is internal so you can run a massive grid but you can't set independent samplers off with parallel prcoessing. It's demand for required lines is also weird.
    subprocessing is messy and requires some file i/o stuff that I don't like

So as a result, I have actually created SpectralRadex ( for the C2H project and have been using it for H3O+/SO. I'll be "releasing" it as part of the C2H paper with a pypi distribution and documentation site. However, you can install it yourself by cloning the repo if you want to use it now.

SpectralRadex has the following features:

    RADEX rewritten in modern fortran to remove pyradex's common block problem
    An F2PY wrapper for that code so you can run radex through a simple function in python with dictionary inputs
    Non-LTE spectral modelling using the RADEX outputs which I benchmarked against MADCUBA and CASSIS

Parallelization is left to the user so you can run a grid in parallel and set samplers off with different MPI jobs as you like.

If this fits the bill for what you wanted (ability to run radex through python and spectral modelling in python) then I'd be keen to get you up and running.

So, what Jon has created is a spectrum generator, while what we need is a spectrum fitter. Here is what I described to Jon as what I wanted for a generalized pyradex modeler:
Thanks Jon.  Really glad I asked you this question as I had wondered if someone had developed a non-LTE spectrum generator based on RADEX.  What I have been working on works similarly but produces visualization of the modeled fit ranges rather than spectra:

    You have measurements of integrated intensities for several transitions from a molecule that you want to use to constrain (n,N,T).
    Run RADEX models over a pre-defined range of (n,N,T) to determine the best-fit and range (in a chi-squared sense) for (n,N,T).  A plot of an example where I fit 5 transitions from HCN (four ratios) to constrain (n,T) is attached.

I often have the need to just get a quick estimate of the (n,N,T) range implied by a set of spectral line measurements, if only to justify going down a certain interpretation path.  pyradex has built-in the modeling plus fitting capabilities, but it looks like SpectralRadex has the underpinnings to do this better (with a less obtuse implementation of RADEX, for example).  Do I have that right?

This is a task that I was going to hand-off to Erica, which she will likely get to in the spring or early summer.

-- Jeff
To which Jon replied:
Hi Jeff,

That's something we can do easily. I just need to add a function to do the actual fitting and that's definitely something I'd like to include. I can get this to you some time in early January if that's useful?

On 2021-01-21 Jon let me know that he has made the changes noted above...
In fact, after a successful day of navigating the world of python development, I can update my instructions.

You should be able to "pip install spectralradex" to get the python module and our docs are live at

On Jan 22 2021, at 11:31 am, Jon Holdship <> wrote:

    Hi Jeff,

    A little before Christmas, you asked about running RADEX grids in python. A PhD student in our group and I have added a few features to our python module to do just that!

    I'm struggling to get the documentation to compile but you may eventually find it at . For now, you can find the attached page as it compiled locally.

    We'll also eventually get to the point where you can "pip install spectralradex". However, currently you need to follow the instructions on our github: . It's super simple just install it to a location following the command in the README and then make sure you add it to the PYTHONPATH environmental variable.

    If you have any questions, let me know!

Looks like this is ready to experiment with.

2021-07-01: VLA Continuum Source Extraction and Spectral Index Measurement

Reduced images located at /lustre/cv/users/jmangum/NGC253/VLACont. The text file ReductionFileLocationDirectory.txt summarizes locations and names for final images. There is also an Overleaf file with data reduction notes that can be used as a starter for a potential future article.

2021-07-01: Get Familiar with pyBDSF for Generalized Source Extraction

pyBDSF information, installation, and examples. Will need this for VLA continuum data reduction and, in the future, for a more generalized molecular cloud identification in NGC253 and elsewhere.

2021-10-01: Generalized Rotation Diagram Package

Starting with a generalized LTE model of a molecule's emission that Adam Ginsburg wrote a while back develop a tool which generates a rotation diagram, which is used to get best-fit estimate of the LTE molecular column density and excitation temperature for a molecule. Since the LTE column density for a molecule is given by N = constant*IntIntensity*exp(Eu/Tex), you can take the logarithm of N/constant, which results in a linear function involving constant*IntIntensity and 1/Tex, so that the y-intercept is the column density and the slope is the excitation temperature. You can flesh this out with the appropriate column density equations in Mangum & Shirley 2015.

-- JeffMangum - 2020-12-25

This topic: Main > TWikiUsers > JeffMangum > EricaResearchProjects > EricaTaskList
Topic revision: 2021-05-04, JeffMangum
This site is powered by FoswikiCopyright © by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding NRAO Public Wiki? Send feedback