Skip to content

Radiative transfer and synthetic radio imaging of ionised jets from young forming stars

Notifications You must be signed in to change notification settings

SimonP2207/RaJePy

Repository files navigation

RaJePy

Overview

Radio Jets in Python (RaJePy) is a Python library which conducts radiative transfer (RT) calculations towards a power-law-based, physical model of an ionised jet. Data products from those calculations are subsequently used as the models to conduct synthetic, interferometric radio imaging upon. Both continuum and radio recombination line (RRL) observations are supported, though be aware that RRL RT calculations currently assume local thermodynamic equilibrium (LTE). Incorporating non-LTE effects into the RRL RT calculations is the next priority in RaJePy's development so that physically accurate synthetic RRL images can be derived.

Purpose

  • Inform radio astronomers on the significance of the detrimental effects of the interferometric imaging of ionised jets
  • Allow observers to determine the best telescope configurations and observing frequencies for their science target
  • Determine the spectral and morphological evolution of jets with variable ejection events/mass loss rates
  • Determine RRL-observability with current and future observers

Instructions for use

RaJePy's operation is based upon the interplay between the JetModel and Pipeline classes defined in classes.py and initialised in the namespace of the RaJePy module. Basic use can be achieved via the command line with a command of the following syntax:

python /home/user/RaJePy/main.py -rt -so -v example-model-params.py example-pipeline-params.py

Help for using the command line can be found by using the command python /home/user/RaJePy/main.py --help, which gives the following output:

usage: main.py [-h] [-v] [-rt] [-so] [-r] [-c] model_param_file pipeline_param_file

positional arguments:
  model_param_file      Full path to model parameter file
  pipeline_param_file   Full path to pipeline parameter file

optional arguments:
  -h, --help            show this help message and exit
  -v, --verbose         Increase output verbosity
  -rt, --radiative-transfer
                        Compute radiative transfer solutions
  -so, --simobserve     Conduct synthetic observations using CASA
  -r, --resume          Resume previous pipeline run if present
  -c, --clobber         Overwrite any data products/files present

Model parameter file

For the model parameter file, a full-working example is given in RaJePy/files/example-model-params.py which contains the following:

import numpy as np

params = {
    "target": {"name": "test2",  # Jet/YSO/Model name
               "ra": "04:31:34.07736",  # R.A. [HH:MM:SS.SS...]
               "dec": "+18:08:04.9020",  # Declination [DD:MM:SS.SS...]
               "epoch": "J2000",  # Coordinate epoch, J2000 or B1950
               "dist": 120.,  # Distance from observer [pc]
               "v_lsr": 6.2,  # System's local standard of rest velocity [km/s]
               "M_star": 0.55,  # Central object mass [M_sol]
               "R_1": .25,  # inner disc radii sourcing the jet [au]
               "R_2": 2.5,  # outer disc radii sourcing the jet [au]
               },
    "grid": {"n_x": 50,  # No. of cells in x
             "n_y": 400,  # No. of cells in y
             "n_z": 50,  # No. of cells in z
             "l_z": 2.,  # Length of z-axis. Overrides n_x/n_y/n_z
             "c_size": 0.5,  # Cell size [au]
             },
    "geometry": {"epsilon": 7. / 9.,  # Jet width index
                 "opang": 25.,  # Jet opening angle [deg]
                 "w_0": 1.,  # Half-width of jet base [au]
                 "r_0": 1.,  # Launching radius [au]
                 "inc": 90.,  # Inclination angle where 0 <= i <= 90 [deg]
                 "pa": 0.,  # BLUE (approaching) jet position angle [deg]
                 "rotation": "CCW",  # Rotation sense, either "CCW" or "CW"
                 },
    "power_laws": {"q_v": 0.,  # Velocity index
                   "q_T": 0.,  # Temperature index
                   "q_x": 0.,  # Ionisation fraction index
                   "q^d_n": 0.,  # Cross-sectional density index
                   "q^d_T": 0.,  # Cross-sectional temperature index
                   "q^d_v": 0.,  # Cross-sectional jet-velocity index
                   "q^d_x": 0.  # Cross-sectional ionisation fraction index
                   },
    "properties": {"v_0": 150.,  # Ejection velocity [km/s]
                   "x_0": 0.1,  # Initial HII fraction
                   "T_0": 1E4,  # Temperature [K]
                   "mu": 1.3,  # Mean atomic weight [u]
                   "mlr_bj": 1e-7,  # BLUE jet steady-state MLR [Msol/yr]
                   "mlr_rj": 5e-8,  # RED jet steady-state MLR [Msol/yr]
                   },
    "ejection": {"t_0": np.array([0.5, 0.75, 1., 2.]),  # Peak times of bursts [yr]
                 "hl": np.array([0.15, 0.15, 0.45, 0.5]),  # Half-lives of bursts [yr]
                 "chi": np.array([5., 5., 2.5, 10.]),  # Burst factors
                 "which": np.array(["R", "B", "B", "RB"])  # Which jet the burst relates to
                 }
    }

This shows a python dict containing 6 keys associated with more nested dicts. For each of those 6 keys, their logical purpose and description of their associated values' dict's keys are given in separate tables below.

Model-parameter section 'target'

'Science' target information.

Parameter/key Description Type Example
"name" Target object name str "S255 IRS3"
"ra" Target object Right Ascension (HH:MM:SS.S) str "01:23:45.67"
"dec" Target object Declination (DD:MM:SS.S) str "-87:65:43.2"
"epoch" Epoch of RA/Dec str "J2000"
"dist" Distance (pc) float 1780.
"v_lsr" Systemic velocity with respect to the Local Standard of Rest (km/s) float -4.2
"M_star" Central protostellar mass (solar masses). float 10.0
"R_1" Inner disc radius from which jet material is sources (au) float 1.0
"R_2" Outer disc radius from which jet material is sources (au) float 10.0

Model-parameter section 'grid'

Model grid dimensions

Parameter/key Description Type Example
"n_x" Number of cells in x-axis int 100
"n_y" Number of cells in y-axis int 100
"n_z" Number of cells in z-axis int 400
"l_z" Plane-of-the-sky length of bi-polar jet (arcsec). Overrides "n_x"/"n_y"/"n_z" parameters float, None 2.0
"c_size" Grid cell size (au) float 2.0

NB - If not None, "l_z" calculates (using supplied "dist" and "c_size") and updates "n_x"/"n_y"/"n_z" parameters to fully encompass a jet, plane-of-the-sky length of "l_z" arcseconds.

Model-parameter section 'geometry'

Jet geometry parameters

Parameter/key Description Type Example
"epsilon" Power-law coefficient for jet width float +1.0
"opang" Opening angle (deg) float 20.0
"w_0" Jet half-width at jet base (au) float 2.0
"r_0" Launching radius (au) float 4.0
"inc" Jet inclination (deg) float 90.
"pa" Jet position angle (deg) float 0.
"rotation" Rotation sense of the jet and counter-jet, around the propagation (r) axis. Must be either "CCW" (counter-clockwise) or "CW" (clockwise) str `CCW

Model-parameter section 'power_laws'

Power-law coefficients defining the physical model in velocity, temperature, ionisiation fraction and number density

Parameter/key Description Type Example
"q_v" Power-law coefficient for jet velocity as function of r float -0.5
"q_T" Power-law coefficient for jet temperature as function of r float -0.5
"q_x" Power-law coefficient for jet ionisation fraction as function of r float 0.0
"q^d_n" Cross-sectional power-law coefficient for jet number density as function of w float -2.0
"q^d_T" Cross-sectional power-law coefficient for jet electron temperature as function of w float -2.0
"q^d_v" Cross-sectional power-law coefficient for jet velocity as function of w float -0.5
"q^d_x" Cross-sectional power-law coefficient for jet ionisation fraction as function of w float -2.0

Model-parameter section 'properties'

Jet physical parameter values

Parameter/key Description Type Example
"v_0" Jet initial velocity (km/s) float 500.
"x_0" Initial jet ionisation fraction (0 < x_0 <= 1) float 0.1
"T_0" Initial jet temperature (K) float 1e4
"mu" Mean atomic weight of jet (hydrogen atom mass) float 1.3
"mlr_bj" Blue-jet mass loss rate (solar masses per yr) float, None 1e-5
"mlr_rj" Red-jet mass loss rate (solar masses per yr) float, None 1e-5

Model-parameter section 'ejection'

Jet mass loss variability parameters. The number of bursts defined here does not impact code-performance and so complex mlr(t) profiles can be defined dependent upon requirements.

Parameter/key Description Type Example
"t_0" Burst(s) peak times (yr) np.array (dtype=float) np.array([0., 1., 2.])
"hl" Burst(s) 'half-lives', i.e. FWHM in time (yr) np.array (dtype=float) np.array([0.2, 0.1, 0.8])
"chi" Burst(s) factors (multiple of jet's steady state mass loss rate np.array (dtype=float) np.array([10., 5., 2.])
"which" Whether the burst is in the blue ('B'), red ('R'), or both ('RB') jets np.array (dtype=str) np.array(['R', 'B', 'RB'])

Pipeline parameter file

For the pipeline parameter file, a full-working example is given in RaJePy/files/example-pipeline-params.py which contains the following:

import os
import numpy as np

params = {'min_el':    20.,    # Minimum elevation for synthetic observations, deg
          'dcys':      {"model_dcy": os.sep.join([os.getcwd(), 'test_output_dcy'])},  # Output root directory
          # Continuum observations
          'continuum': {'times':  np.array([0.]),  # Model times, yr
                        'freqs':  np.array([1.5, 3.0, 6., 10., 22., 33., 43., 5.05]) * 1e9,  # Frequencies of observations, Hz
                        't_obs':  np.array([1200, 1200, 3600, 1200, 1200, 1800, 2400, 59400]),  # Total on-source times, s
                        'tscps':  np.array([('VLA', 'A'), ('VLA', 'A'),
                                            ('VLA', 'A'), ('VLA', 'A'),
                                            ('VLA', 'A'), ('VLA', 'A'),
                                            ('VLA', 'A'), ('EMERLIN', '0')]),  # List of 2-tuples of (telescope, configuration)
                        't_ints': np.array([5, 5, 5, 5, 5, 5, 5, 5]),  # Visibility integration times, s
                        'bws':    np.array([1e9, 2e9, 2e9, 4e9, 4e9, 4e9, 8e9, .5e9]),  # Observational bandwidth, Hz
                        'chanws': np.array([1e8, 1e8, 1e8, 1e8, 1e8, 1e8, 1e8, 2.e8])},  # Channel widths, Hz
          # Radio recombination line observations
          'rrls':      {'times':  np.array([0.]),  # Model times, yr
                        'lines':  np.array(['H58a']),  # RRL lines to observe (Element+n+dn)
                        't_obs':  np.array([30000]),  # Total on-source times, s
                        'tscps':  np.array([('VLA', 'A')]),  # List of 2-tuples of (telescope, configuration)
                        't_ints': np.array([60]),  # Visibility integration times, s
                        'bws':    np.array([1e8]),  # Observational bandwidth, Hz
                        'chanws': np.array([1e5])},  # Channel widths, Hz
          }

This shows a python dict with 4 associated keys which are:

Parameter/key Description Type Example
"min_el" Minimum elevation to conduct synthetic observations (deg) float 20.
"dcys " Contains relevant directories needed float 20.
"continuum" Contains all continuum run specifications dict See below
"rrls" Contains all radio recombination line run specifications dict See below

Pipeline-parameter section 'dcys'

Parameter/key Description Type Example
'model_dcy' Root directory for all products and directories associated with Pipeline execution str /my/rajepy/exec/directory

Pipeline-parameter section 'continuum'

Imaging/pipeline parameters defining continuum radiative-transfer/synthetic observations are defined here. When executed, radiative transfer and/or synthetic observations are conducted at each value of "times", for each value of "freqs". For example, if "times" specified are 0, 1 and 10 years, whilst the observing frequencies are specified at 2, 10, 20 and 50 GHz, a total of 12 radiative transfer/synthetic observations will be conducted.

Parameter/key Description Type Example
"times" Observational times (yr) np.array(dtype=float) np.linspace(0., 5., 4)
"freqs" Observational intermediate frequencies (Hz) np.array(dtype=float) np.array([1e9, 5e9, 2e10, 5e10])
"t_obs" Total scan times (s) np.array(dtype=int) np.array([3600, 1800, 900, 450])
"tscps" Telescope names and configurations np.array(dtype=tuple(str, str)) np.array([('VLA', 'A'), ('EMERLIN', '0'), ('VLA', 'B'), ('VLA', 'C')])
"t_ints" Visibility integration times (s) np.array(dtype=int) np.array([5, 3, 3, 2])
"bws" Bandwidths (Hz) np.array(dtype=float) np.array([0.5e9, 2e9, 2e9, 4e9])
"chanws" Channel widths within bandwidth (Hz) np.array(dtype=int) np.array([1, 1, 1, 1])

NB - Specified arrays for "freqs", "t_obs", "tscps", "t_int", "bws" and "chanws" should all have the same length. Array values at the same indices in each of those arrays are used in combination to define the observational parameters for one run.

Pipeline-parameter section 'rrls'

Imaging/pipeline parameters defining radio recombination line radiative-transfer/synthetic observations are defined here. When executed, radiative transfer and/or synthetic observations are conducted at each value of "times", for each value of "lines". For example, if "times" specified are 0, 1 and 10 years, whilst the observed RRL lines are specified as H48a and H58a, a total of 6 radiative transfer/synthetic observations will be conducted.

Parameter/key Description Type Example
"times" Observational times (yr) np.array with dtype=float np.linspace(0., 5., 4)
"lines" Radio recombination lines for observing np.array with dtype=str np.array(['H56a', 'H42a', 'H76a'])
"t_obs" Total scan times (s) np.array with dtype=int np.array([3600, 1800, 900, 450])
"tscps" Telescope names and configurations np.array with dtype=tuple(str, str) np.array([('VLA', 'A'), ('EMERLIN', '0'), ('VLA', 'B'), ('VLA', 'C')])
"t_ints" Visibility integration times (s) np.array with dtype=int np.array([5, 3, 3, 2])
"bws" Bandwidths (Hz) np.array with dtype=float np.array([0.5e9, 2e9, 2e9, 4e9])
"chanws" Channel widths within bandwidth (Hz) np.array with dtype=int np.array([1, 1, 1, 1])

RaJePy automatically calculates the central observing frequency based upon the RRL(s) defined above. See the splatalogue for useful RRL information.

NB - Specified arrays for "lines", "t_obs", "tscps", "t_int", "bws" and "chanws" should all have the same length. Array values at the same indices in each of those arrays are used in combination to define the observational parameters for one run.

Data products

If /example/dcy was given for params['dcys']['model_dcy'] in the pipeline parameter file, the output from the above code would be:

/example/dcy/  # Set by user in pipeline parameter file's params['dcys']['model_dcy']
├──ModelRun_YYYY-MM-DD-HH:MM:SS.log  # Generated log file by Pipeline
├──jetmodel.save  # Saved JetModel instance (large file size!)
├──modelrun.save  # Saved Pipeline instance
├──pointings.ptg  # Synthetic observation pointing file
├──GridPlot.pdf  # Plot of defined jet-grid profile from each aspect (x, y and z)
├──JMLPlot.pdf  # Plot of jet mass loss rate as a function of time
├──Day0/  # First epoch directory for modelling/observations
│  ├──ModelPlot.pdf  # Plot of resulting physical model
│  ├──RadioPlot.pdf  # Plot of model fluxes, optical depths and emission measures
│  ├──2GHz/  # Directory containing data products for first radio frequency
│  │  ├──DDMMYYYY_HHMMSS.py  # Casa script to be executed, generated by Pipeline
│  │  ├──DDMMYYYY_HHMMSS.log  # Casa logfile for execution of DDMMYYYY_HHMMSS.py
│  │  ├──EM_Day0_1000MHz.fits  # Model emission measure image
│  │  ├──Flux_Day0_1000MHz.fits  # Model radio flux image
│  │  ├──Tau_Day0_1000MHz.fits  # Model optical depth image
│  │  ├──SynObs.vla.a.noisy.imaging.fits  # Final clean image
│  │  ├──SynObs.vla.a.noisy.imaging.estimates  # Input to casa.imfit (continuum only)
│  │  ├──SynObs.vla.a.noisy.imaging.imfit  # Results of Gaussian fits to jet (continuum only)
│  │  └──SynObs/  # All files generated by CASA
│  │     ├──SynObs.vla.a.ms/  # No-noise measurement set (casa.simobserve)
│  │     ├──SynObs.vla.a.noisy.ms/  # Noisy measurement set (casa.simobserve)
│  │     ├──SynObs.vla.a.noisy.imaging.image/  # Clean image (casa.tclean)
│  │     ├──SynObs.vla.a.noisy.imaging.mask/  # Clean mask (casa.tclean)
│  │     ├──SynObs.vla.a.noisy.imaging.model/  # Deconvolved model image (casa.tclean)
│  │     ├──SynObs.vla.a.noisy.imaging.pb/  # Primary beam image (casa.tclean)
│  │     ├──SynObs.vla.a.noisy.imaging.psf/  # Point spread function image (casa.tclean)
│  │     ├──SynObs.vla.a.noisy.imaging.residual/  # Residual map (casa.tclean)
│  │     ├──SynObs.vla.a.noisy.imaging.sumwt/  # Sum-of-weights image (casa.tclean)
│  │     ├──SynObs.vla.a.skymodel/  # Input model image (casa.simobserve)
│  │     └──SynObs.vla.a.skymodel.flat/  # Flattened input model image (casa.simobserve)
│  ├──5GHz/  # Directory containing data products for second radio frequency|  └──...
│  ├──20GHz/  # Directory containing data products for third radio frequency|  └──...
│  ├──50GHz/  # Directory containing data products for fourth radio frequency|  └──...
│  └──H58a/  # Directory containing data products for RRL line H58a
│     └──...
├──Day100/  # Second epoch directory for modelling/observations
│  └──...
└──...

classes.py

classes.JetModel

The purpose of the JetModel class is to initialise/calculate the physical model grid and conduct radiative transfer calulations.

JetModel's __init__ method takes one compulsory argument, params, and two optional keyword-arguments, verbose and log:

  • params : Full path to the model parameters file, an example for which is RaJePy/files/example-model-params.py (see a full description below)
  • verbose: Boolean determining if on an execution of JetModel's class/instance methods, verbose output to the terminal is wanted. Default is True.
  • log: Full path to an optional log file containing a detailed log of all of an instance's method outputs. None by default.

For the initialisation of the jet model's physical parameters and their values in each grid cell, the analytical, power law-based approach of Reynolds (1986) is used. Those parameters described in that work form the basis of the input parameters required to define a jet model in RaJePy.

classes.Pipeline class

The purpose of the Pipeline class is to handle directory/file manipulation and perform synthetic observations via casa, their subsequent measurements and other analyses conducted on a JetModel instance.

Pipeline's __init__ method takes two compulsory arguments, jetmodel and params:

  • jetmodel: JetModel instance to conduct all analyses and observations upon
  • params : Full path to the pipeline parameters file, an example for which is RaJePy/files/example-pipeline-params.py (see a full description below)

Pipeline.execute method

The execute method runs the complete pipeline, producing relevant plots, .fits model files, measurement sets and final clean image .fits files. It takes 5, optional, keyword-arguments:

kwarg Description Type
simobserve Whether to conduct synthetic observations on the model .fits file bool
verbose Verbose output to the terminal? bool
dryrun Whether to execute a dry run, without actually running any calculations (for testing) bool
resume Whether to resume a previously saved run (if saved model file and saved pipeline file exist) bool
clobber Whether to redo and overwrite previously written files [soon to be deprecated] bool

classes.ContinuumRun class

Data class containing all observational properties derived from pipeline param file, for a continuum run.

classes.RRLRun class

Data class containing all observational properties derived from pipeline param file, for an RRL run. Child class of ContinuumRun.

casa.Script class

Class containing all casa tasks to be executed with given parameters, in a defined order. Methods include add_task and execute. The latter method writes the .py file, based upon instance's parameter values, which is executed via the command casa --nogui --nologger --agg --logfile XXX.log -c XXX.py.

Use RaJePy as its own module to execute your own pipeline

After creating instances of the JetModel and Pipeline classes, execution of the desired synthetic observations etc. takes place via the Pipeline class' execute method. A complete script for execution of a synthetic observing and model calculation run would be:

import os
import RaJePy as rjp

model_params = 'model-params.py'
pline_params = 'pipeline-params.py'
log_file = 'TestJet.log'

jm = rjp.JetModel(model_params, log=log_file)
pline = rjp.Pipeline(jm, pline_params)

pline.execute(simobserve=True, dryrun=False)

Requirements:

Python >= 3.6 (f-string dependence)

Python standard library packages:

  • collections.abc
  • errno
  • os
  • pickle (developed with 4.0)
  • shutil
  • sys
  • time
  • warnings

Other Python packages:

  • astropy (developed with 4.0)
  • imageio (developed with 2.8.0)
  • matplotlib (developed with 3.2.1)
  • mpmath (developed with 1.1.0)
  • numpy (developed with 1.18.1)
  • pandas (developed with 1.0.5)
  • scipy (developed with 1.3.1)
  • tabulate (developed with 0.8.9)

System-based installations

  • Working casa installation (developed with 6.1.2.7) with casa's executable located in $PATH

Future work and direction

  • Incorporate inclination into jet model
  • Incorporate position angle into jet model
  • Implement more than one channel across bandwidth for more accurate multi-frequency synthesis
  • Incorporate asymmetric mass-loss for bursts and steady-state mass loss rates
  • Incorporate non-LTE effects into RRL radiative transfer calculations
  • Parallelise code, especially different synthetic observations and model calculations

About

Radiative transfer and synthetic radio imaging of ionised jets from young forming stars

Topics

Resources

Stars

Watchers

Forks

Releases

No releases published

Packages

No packages published

Languages