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.
- 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
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
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 dict
s. For each of those 6 keys, their logical purpose and description of their associated values' dict
's keys are given in separate tables below.
'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 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.
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 |
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 |
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 |
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']) |
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 |
Parameter/key | Description | Type | Example |
---|---|---|---|
'model_dcy' |
Root directory for all products and directories associated with Pipeline execution |
str |
/my/rajepy/exec/directory |
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.
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.
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
│ └──...
└──...
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 isRaJePy/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 isTrue
.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
.
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 uponparams
: Full path to the pipeline parameters file, an example for which isRaJePy/files/example-pipeline-params.py
(see a full description below)
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 |
Data class containing all observational properties derived from pipeline param file, for a continuum run.
Data class containing all observational properties derived from pipeline param file, for an RRL run. Child class of ContinuumRun
.
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
.
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)
Python >= 3.6 (f-string dependence)
- collections.abc
- errno
- os
- pickle (developed with 4.0)
- shutil
- sys
- time
- warnings
- 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)
- Working casa installation (developed with 6.1.2.7) with casa's executable located in
$PATH
Incorporate inclination into jet modelIncorporate position angle into jet modelImplement more than one channel across bandwidth for more accurate multi-frequency synthesisIncorporate 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