Wakeflow User API

Module contents

Wakeflow allows users to generate and manipulate semi-analytic models of planet wakes.

Wakeflow module

Contains the WakeflowModel class, intended for use by users to generate, configure and run models of planet wakes.

class wakeflow.wakeflow.WakeflowModel

Bases: object

Model object allowing you to configure, generate and save planet wake model results.

model_params

Dictionary containing the parameters to be used in the model as specified by the user.

Type

dict

default_params

Dictionary containing the default WakeflowModel parameters.

Type

dict

_run_wakeflow(params: _Parameters) → None

Internal use method for generating the planet wake by calling other parts of Wakeflow.

Parameters

params (Parameters) – Paramaters object specfiying options for the model to be generated.

configure(name: str = 'default_results_dir', system: str = 'default_parent_dir', m_star: float = 1.0, m_planet: float = 0.1, r_outer: float = 500, r_inner: float = 100, r_planet: float = 250, phi_planet: float = 0, r_ref: Optional[float] = None, r_c: float = 0, z_max: float = 3, q: float = 0.25, p: float = 1.0, hr: float = 0.1, dens_ref: float = 1.0, cw_rotation: bool = False, grid_type: float = 'cartesian', n_x: int = 400, n_y: int = 400, n_r: int = 200, n_phi: int = 160, n_z: int = 50, r_log: bool = False, make_midplane_plots: bool = True, show_midplane_plots: bool = True, dimensionless: bool = False, include_planet: bool = True, include_linear: bool = True, save_perturbations: bool = True, save_total: bool = True, write_FITS: bool = False, run_mcfost: bool = False, inclination: float = -225, PA: float = 45, PAp: float = 45, temp_star: float = 9250, distance: float = 101.5, v_max: float = 3.2, n_v: float = 40) → None

Configure the model according the specifications by the user.

Parameters

• name (str) – results are saved a directory called name.

• system (str) – name of parent directory of results.

• m_star (float) – central star mass in solar masses.

• m_planet (float) – planet mass in Jupiter masses. Can also give list(float) to run multiple planet masses.

• r_outer (float) – outer disk radius in au.

• r_inner (float) – inner disk radius in au.

• r_planet (float) – orbital radius of planet in au.

• phi_planet (float) – azimuthal position of the planet in radians.

• r_ref (float) – reference radius r_ref in au.

• r_c (float) – critical radius r_c in au, used for exponentially tapered density profile. ignored if set to 0.

• z_max (float) – height of the disk in units of pressure scale height at r_outer.

• q (float) – q index for sound speed profile, defined as c_s propto r^{-q}.

• p (float) – p index for density profile (NOT surface density), defined as rho propto r^{-p} or propto r^{-p}*exp{-r/r_c}^{2-p}.

• hr (float) – disk aspect ratio at r_ref.

• dens_ref (float) – density at r_ref in g/cm^3.

• cw_rotation (bool) – clockwise disk rotation (True) or anticlockwise disk rotation (False).

• grid_type (float) – grid type, either “cylindrical” or “cartesian” or “mcfost”. Must choose “mcfost” to write a .FITS file.

• n_x (int) – number of grid points in x.

• n_y (int) – number of grid points in y.

• n_r (int) – number of grid points in radius.

• n_phi (int) – number of grid points in azimuth.

• n_z (int) – number of grid points in z.

• r_log (bool) – True for logarithmically spaced radii, False for linear. Ignored for “mcfost” grid type.

• make_midplane_plots (bool) – Create midplane plots of density and velocity perturbations?

• show_midplane_plots (bool) – Display midplane plots of density and velocity perturbations to user?

• dimensionless (bool) – True for dimensionless results otherwise False.

• include_planet (bool) – Include the planet-induced perturbations? Set to False to generate unperturbed disk.

• include_linear (bool) – Include the results from the linear regime?

• save_perturbations (bool) – Save the perturbations?

• save_total (bool) – Save the totals (perturbations + background disk)?

• write_FITS (bool) – Generate a .FITS file to run in MCFOST? Requires “mcfost” grid type.

• run_mcfost (bool) – Call MCFOST to generate line emission cubes on results, requires write_FITS=True.

• inclination (float) – Inclination in degrees for MCFOST.

• PA (float) – Position angle in degrees for MCFOST.

• PAp (float) – Planet position in degrees angle for MCFOST.

• temp_star (float) – Star temperature in Kelvin for MCFOST.

• distance (float) – System distance in parsecs for MCFOST.

• v_max (float) – maxmimum magnitude velocity for cube in MCFOST.

• n_v (float) – number of velocity channels for MCFOST.

configure_from_file(param_file: str) → None

Configure the model by reading in .yaml file provided by user.

Parameters

param_file (str) – path to .yaml file provided by user specifying run parameters.

run(overwrite: bool = False) → None

Generate results for model, requires user to have called either configure or configure_from_file first.

Parameters

overwrite (bool) – Overwrite previous results with identical name?

Spiralmoments User API

Sub-module contents

Spiralmoments allows users to map and project spiral-shapes and velocity fields in 3 dimensions to predict peak-velocity maps and line-of-sight spiral shapes.

Spiralmoments sub-module

Contains the Spiral class that allows uses to manipulate, project and rotate different spiral shapes easily.

class spiralmoments.spirals.Spiral(radii: ndarray, phi_function: Callable)

Bases: object

Class representing spiral shape, to be instantiated by user. Allows for easy 3D projection and rotation to find the line-of-sight shape for a spiral on the surface of a disk.

_phi = None

_rad = None

_x = None

_y = None

_z = None

_z_disk = None

get_height(height_func: Callable) → None

Calculate the height of the spiral on the disk-surface using function that gives the disk height as a function of radius.

Parameters:

height_func : Callable function that takes disk radius as an argument and returns the respective disk height

returns : None

property height: ndarray

Get height of the spiral on the disk surface. Requires that get_height has been called previously.

Parameters:

returns : numpy.ndarray containing the z-coordinates of the spiral in the disk frame.

perform_rotation(position_angle: float, inclination: float, azimuth_offset: float) → None

Rotate the spiral shape to match the position on the sky, given a disk position angle and inclination. Also rotates the disk around its own z-axis by azimuth_offset (to match a planet position or similar).

Parameters:

position_angle : disk position angle in radians, defined counter-clockwise from the positive x-axis direction inclincation : disk inclination in radians azimuth_offset : angle through which to rotate the spiral around the disk z-axis.

returns : None

property rad_phi: Tuple[ndarray, ndarray]

Get radial and azimuthal disk coordinates of the spiral.

Parameters:

returns : Tuple containing numpy.ndarrays, where the first is the radial disk coordinates of the spiral, and the second is the azimuthal.

property sky_coords: Tuple[ndarray, ndarray, ndarray]

Get the sky coordinates of the spiral. Requires that perform_rotation has been called.

Parameters:

returns : Sky coordinates of spiral.

property xy: Tuple[ndarray, ndarray]

Get the cartesian disk coordinates of the spiral.

Parameters:

returns: Tuple containing numpy.ndarrays, where the first is the disk x coordinates of the spiral, and the second is the y coordinates.

Contains the HeightFunctions class, which allows users to generate functions that return the disk height as a function of radius, using the parameterisation of their choice.

class spiralmoments.height.HeightFunctions

Bases: object

Class containing methods to generate functions that allow users to get disk height as a function of radius (most-likely used for the height of the emitting layer for a certain line).

static interpolate_tabulated() → Callable

Height function defined by interpolating over user-specified values. NOT IMPLEMENTED.

static powerlaw(z_ref: float, r_ref: float, power_index: float) → Callable

Power-law height profile height = z_ref * (r / r_ref)**power_index.

Parameters

• z_ref (disk height at r_ref) –

• r_ref (reference radius) –

• power_index (power-law index for the height profile) –

• returns (Callable function that will return the height value at a provided radius, for the power-law profile.) –

static powerlaw_tapered(z_ref: float, r_ref: float, power_index: float, r_taper: float, power_index_taper: float) → Callable

Power-law height profile height = z_red * (r / r_ref)**power_intex * exp(- (r / r_taper)**power_index_taper).

Parameters

• z_ref (disk height at r_ref) –

• r_ref (reference radius) –

• power_index (power-law index for the height profile) –

• r_taper (taper radius) –

• power_index_taper (power-law index for the exponential taper) –

• returns (Callable function that will return the height value at a provided radius, for the tapered power-law profile.) –

Contains the PhiFunctions class, which allows users to generate functions that return the azimuthal position of a spiral as a function of radius, using the parameterisation of their choice. Includes options to use physical spirals from the Lin-Shu dispersion relation, including the planet wake shape (Ogilvie and Lubow 2002), as well as parameterised spirals such as logarithmic or power-law spirals.

class spiralmoments.phi.PhiFunctions

Bases: object

Class containing methods to generate functions that allow users to get spiral shape functions phi(r) using either simple or physically-motivated parameterisations. In the latter case, phi(r) can be generated from arbitrary rotation and sound speed curves.

static log(b: float, c: float) → Callable

Logarithmic spiral, defined by phi = b * log(r) + c. It has constant pitch angle arctan(1 / b).

Parameters

• b (spiral scaling factor (pitch angle = arctan(1 / b))) –

• c (constant offset in phi) –

• returns (Callable function that will return the phi value at a provided radius, given the defined spiral shape.) –

static m_mode_spiral_numerical(phi_0: float, r_0: float, m: int, Omega_func: float, cs_func: float, rotation_sign: int) → Callable

Individual spiral mode with azimuthal wave number m, defined in terms of the rotation curve Omega(r) and the sounds speed cs(r). Keep in mind that there are evanescent regions interior to the locations of the Lindblad radii. TODO: add a useful reference here.

Parameters

• phi_0 (azimuthal launching location for the spiral wave) –

• r_0 (radial launching location for the spiral wave) –

• m (azimuthal wave number of the spiral wave) –

• Omega_func (function that takes only r as an argument and returns the rotation Omega at r) –

• cs_func (function that takes only r as an argument and returns the sound speed cs at r) –

• rotation_sign (sign of the rotation. +1 for clockwise and -1 for anticlockwise.) –

• returns (Callable function that will numerically evaluate and return the phi value at a provided radius, given the origin and mode of the wave.) –

static planet_wake_numerical(phi_planet: float, r_planet: float, Omega_func: float, cs_func: float, rotation_sign: int) → Callable

Planet wake form given in Rafikov (2002), defined in terms of the rotation curve Omega(r) and the sounds speed cs(r), calculated numerically.

Parameters

• phi_planet (azimuthal location of the planet in radians) –

• r_planet (orbital radius of the planet) –

• Omega_func (function that takes only r as an argument and returns the rotation Omega at r) –

• cs_func (function that takes only r as an argument and returns the sound speed cs at r) –

• rotation_sign (sign of the rotation. +1 for clockwise and -1 for anticlockwise.) –

• returns (Callable function that will numerically evaluate and return the phi value at a provided radius, given the defined planet wake shape.) –

static planet_wake_powerlaw_disk(phi_planet: float, r_planet: float, hr_planet: float, index_cs: float, rotation_sign: int) → Callable

Ogilvie and Lubow (2002) analytical planet wake shape for a power-law disk with Keplerian rotation.

Parameters

• phi_planet (azimuthal location of the planet in radians) –

• r_planet (orbital radius of the planet) –

• hr_planet (disk aspect ratio H/r at the planet orbital radius r=r_planet. H is the pressure scale height H = cs / Omega.) –

• index_cs (power law index for the radial power-law parameterisation of the sound speed in the disk cs = cs_planet * (r / r_planet)**-index_cs) –

• rotation_sign (sign of the rotation. +1 for clockwise and -1 for anticlockwise.) –

• returns (Callable function that will return the phi value at a provided radius, given the defined planet wake shape.) –

static powerlaw(a: float, c: float, power_index: float) → Callable

Power-law spiral, defined by phi = a * r**power_index + c. The pitch-angle evolves as a power-law. TODO: add form for pitch-angle

Parameters

• a (spiral scaling factor) –

• c (constant offset in phi) –

• returns (Callable function that will return the phi value at a provided radius, given the defined spiral shape.) –

static powerlaw_plus_log(a: float, b: float, c: float, power_index: float) → Callable

Composite spiral with both a power-law and a logarithmic component. Defined by phi = a * r**power_index + b * log(r) + c. Corresponds to a spiral that evolves with a pitch-angle with both a constant and a power-law component. TODO: add form for pitch-angle

Parameters

• a (power-law spiral scaling factor) –

• b (log spiral scaling factor (pitch angle approaches arctan(1 / b) as r approaches zero)) –

• c (constant offset in phi) –

• returns (Callable function that will return the phi value at a provided radius, given the defined spiral shape.) –

Wakeflow Dev API

wakeflow.burgers

Contains the solve_burgers function used in the wake non-linear propagation calculations.

wakeflow.burgers._solve_burgers(eta, profile, gamma, beta_p, C, CFL, eta_tilde, t0, linear_solution, linear_t, show_teta, tf_fac)

Propagate the wake in (t,eta,chi) space by solving Eq. 10 from Bollati et al. 2021 using Godunov scheme.

wakeflow.grid

Contains the Grid class on which all Wakeflow models are run/stored.

class wakeflow.grid._Grid(parameters: _Parameters)

Bases: object

Grid object to generate and store Wakeflow results on.

Note that all 3D arrays in the Grid have dimensions (x,z,y) or (phi,z,r).

p

Parameters object (model_setup.py) containing the options used for the model.

Type

Parameters

info

Contains the type, size and contents of the grid.

Type

dict

height

Scale height of the disk.

Type

ndarray

v_r

Contains radial velocities.

Type

ndarray

v_phi

Contains azimuthal velocities, positive is counterclockwise.

Type

ndarray

rho

Contains densities.

Type

ndarray

x

1D. Contains x coordinates of Cartesian grid.

Type

ndarray

y

1D. Contains y coordinates of Cartesian grid.

Type

ndarray

z_xy

1D. Contains z coordinates of Cartesian grid.

Type

ndarray

X

3D. Contains x coordinates of Cartesian grid points.

Type

ndarray

Y

3D. Contains y coordinates of Cartesian grid points.

Type

ndarray

Z_xy

3D. Contains z coordinates of Cartesian grid points.

Type

ndarray

R_xy

3D. Contains radial coordinates of Cartesian grid points.

Type

ndarray

PHI_xy

3D. Contains azimuthal coordinates of Cartesian grid points.

Type

ndarray

r

1D. Contains radial coordinates of Cylindrical or Mcfost grid.

Type

ndarray

phi

1D. Contains azimuthal coordinates of Cylindrical or Mcfost grid.

Type

ndarray

z

1D. Contains z coordinates of Cylindrical grid.

Type

ndarray

R

3D. Contains radial coordinates of Cylindrical or Mcfost grid points.

Type

ndarray

PHI

3D. Contains azimuthal coordinates of Cylindrical or Mcfost grid points.

Type

ndarray

Z

3D. Contains z coordinates of Cylindrical or Mcfost grid points.

Type

ndarray

_add_linear_perturbations(LinearPerts: _LinearPerts, rho_background: _Grid.rho) → None

Add results from the linear regime, stored in LinearPerts object, to the grid

Parameters

• LinearPerts (LinearPerts) – LinearPerts object containing results from the linear regime.

• rho_background (Grid.rho) – unperturbed density from Grid object where make_keplerian_disk has been used.

_add_non_linear_perturbations(NonLinearPerts: _NonLinearPerts, rho_background: _Grid.rho) → None

Add results from the non-linear regime, stored in NonLinearPerts object, to the grid

Parameters

• NonLinearPerts (NonLinearPerts) – NonLinearPerts object containing results from the non-linear regime.

• rho_background (Grid.rho) – unperturbed density from Grid object where make_keplerian_disk has been used.

_flip_results() → None

Go from counterclockwise to clockwise rotation; should be used last

_get_r_phi_coords() → None

Find the (r,phi) coordinates corresponding to each (x,y) point

_make_cartesian_grid() → None

Generates a Cartesian grid using parameters in self.p

_make_cylindrical_grid() → None

Generates a Cylindrical grid using parameters in self.p

_make_empty_disk() → None

fill v_r, v_phi and rho with zeros

_make_grid() → None

Generates grid according to parameters in self.p and define disk scale height

_make_keplerian_disk() → None

Generate unperturbed accretion disk in vertical hydrostatic equilibrium, according to parameters in self.p

_make_mcfost_grid() → None

Generates an Mcfost grid using parameters in self.p

_merge_grids(grid_to_merge: _Grid) → None

Merge v_r, v_phi and rho from two compatible Grid objects (ie. with the same _Parameters)

_remove_dimensions(scale_dens: bool = False) → None

Convert results to be dimensionless by appropriate scaling; should be used last

Parameters

scale_dens (bool, optional) – If true, the density is scaled by the rho_ref value in the parameters. Should be left as default for scaling perturbations (delta rho / rho) and set to True otherwise.

_save_results(label: str, printed: str) → None

Save results and grid to .npy files in results directory

Parameters

• label (str) – Prepended to saved file names (used in wakeflow.py to give “delta” for perturbations and “total” for perturbations applied on top of an unperturbed disk).

• printed (str) – Message to be printed to the user in format f”{printed} saved to {savedir}”.

_show_disk2D(z_slice: int, show: bool = False, save: bool = False, dimless: bool = False, vphi_lim: float = 1.0, cw: bool = False) → None

Generate plots of the disk at constant z value with index z_slice

For example, mid-plane plots are generated with z_slice=0.

Parameters

• z_slice (int) – Index in Grid of z_slice to plot.

• show (bool, optional) – True to show plots to user using matplotlib.pyplot.show(). Default is False.

• save (bool, optional) – True to save plots in directory with results. Default is False.

• dimless (bool, optional) – Set to True if using dimensionless results to get accurately labelled plots.

• vphi_lim (float, optional) – Modifier for the limits on the v_phi plot.

• cw (bool, optional) – True for clockwise rotation and False for anticlockwise. Since anticlockwise is default, False is default.

_smooth_box(big_box_grid: _Grid) → None

Under development. Smooths the solution between the linear and non-linear regimes. Currently only smooths in v_r (it would be easy to extend to the other components if you need it).

_write_fits_file() → None

Write a .FITS file from the results compatible with being run in MCFOST

wakeflow.linear_perts

Contains the LinearPerts class responsible for handling the linear regime of the models.

class wakeflow.linear_perts._LinearPerts(parameters: _Parameters)

Bases: object

Class to store the results from the linear regime nearby the planet.

_cut_box_annulus_segment() → None

Extract the part of the linear solution needed for the model, and interpolate onto appropriate grid

_cut_box_square() → None

Deprecated. This is only used in the case where you want to plot the linear solution in t, eta space

wakeflow.mcfost_interface

Contains functions for generating MCFOST grid information and writing MCFOST parameter files. Requires pymcfost which may be found here: https://github.com/cpinte/pymcfost

wakeflow.mcfost_interface._make_mcfost_grid_data(parameters: _Parameters) → None

Take grid parameters specified by Wakeflow user and call MCFOST to generate grid data.

wakeflow.mcfost_interface._make_mcfost_parameter_file(parameters: _Parameters) → None

Take Wakeflow parameters and write an MCFOST .para file by calling pymcfost.

wakeflow.mcfost_interface._read_mcfost_grid_data(parameters: _Parameters) → Tuple[ndarray, ndarray]

Read in grid data generated by MCFOST and return radial and z coordinates.

wakeflow.model_setup

Contains the Parameters class responsible for storing the parameters used in a Wakeflow model. Additionally, contains functions used read in said parameters, check their validity, and create/replace the results directory.

class wakeflow.model_setup._Constants

Bases: object

Constants class for storing physical constants needed by Wakeflow. Inherited by _Parameters.

class wakeflow.model_setup._Parameters(config: dict)

Bases: _Constants

Class for storing Wakeflow model parameters to be easily handed around by various parts of Wakeflow.

_do_sanity_checks() → bool

Called to check that the parameters specified by the user will not break the results/run. Some parameter combinations are incompatible.

wakeflow.model_setup._load_config_file(config_file: str, default_config_dict: Optional[dict] = None) → dict

Reads .yaml parameter file into a dictionary so that Wakeflow may parse it.

wakeflow.model_setup._run_setup(param_dict: dict, overwrite: bool = False) → _Parameters

Perform setup by generating parameters object, checking the parameters are okay, creating results directory, writing .yaml file with parameters used to results dir, and calling MCFOST to generate grid data if needed.

wakeflow.model_setup._write_config_file(config_dict: dict, directory: str, filename: str) → None

Writes parameters dictionary to .yaml file.

wakeflow.non_linear_perts

Contains the NonLinearPerts class responsible for handling the non-linear regime of the models.

class wakeflow.non_linear_perts._NonLinearPerts(parameters: _Parameters, Grid: _Grid)

Bases: object

Class to store the results from the non-linear wake propagation away from the planet.

_extract_ICs(LinearPerts: _LinearPerts) → None

Extract the initial conditions for the non-linear evolution from the linear regime.

_extract_ICs_ann(LinearPerts: _LinearPerts) → None

Alternate initial condition extraction where the IC is read from the edges of the box as included in the final solution. Usually, far more y-extent of the linear regime is used than this. Using this method will invalidate the solution and is meant for developer use only.

_get_non_linear_perts() → None

Calculate the perturbations in the non-linear regime by propagating the wake away from the planet by solving Burger’s eqn and then transforming the results to physical coordinates.

wakeflow.transformations

Contains functions for mapping the results between (t,eta,chi) coordinates and physical coordinates, see Bollati et al. 2021 for details.

wakeflow.transformations._Eta(r, phi, Rp, hr, q, cw)

Eq. (14) Bollati et al. 2021

wakeflow.transformations._Eta_vector(r, phi, Rp, hr, q, cw)

Eq. (14) Bollati et al. 2021

Vectorised version of _Eta. Simply replaces the modular arithmetic which involved if statements with modulus operators and constant offsets.

wakeflow.transformations._Lambda_fu(r, Rp, csp, hr, gamma, q, p)

Eq. (28) Bollati et al. 2021

wakeflow.transformations._Lambda_fv(r, Rp, csp, hr, gamma, q, p)

Eq. (29) Bollati et al. 2021

wakeflow.transformations._g(r, Rp, hr, q, p)

Equation (12) Bollati et al. 2021

wakeflow.transformations._get_chi(pphi, rr, time_outer, time_inner, eta_outer, eta_inner, eta_tilde_outer, eta_tilde_inner, C_outer, C_inner, solution_outer, solution_inner, t0_outer, t0_inner, tf_outer, tf_inner, Rp, x_match, l, cw, hr, q, p, t1)

wakeflow.transformations._get_chi_vector(pphi, rr, tt, time_outer, time_inner, eta_outer, eta_inner, eta_tilde_outer, eta_tilde_inner, C_outer, C_inner, solution_outer, solution_inner, t0_outer, t0_inner, tf_outer, tf_inner, Rp, x_match_l, x_match_r, l, cw, hr, q, p)

This is a vectorised version of the previous _get_chi function. Uses masks rather than logical statements to replace if statements. Uses one RectBivariateSpline interpolation function for all points needing interpolation, rather than recreating a function for every point.

wakeflow.transformations._get_dens_vel(rr, Chi, gamma, Rp, cw, csp, hr, q, p, use_old_vel)

wakeflow.transformations._mod2pi(phi)

wakeflow.transformations._phi_wake(r, Rp, hr, q, cw)

Eq. (4) Bollati et al. 2021

wakeflow.transformations._plot_r_t(params)

wakeflow.transformations._t(r, Rp, hr, q, p)

Equation (43) Rafikov 2002 (Eq. 13 Bollati et al. 2021)

wakeflow.transformations._t_integral(up, q, p)

wakeflow.transformations._t_integrand(x, q, p)

wakeflow.transformations._t_vector(rr, Rp, hr, q, p)

Equation (43) Rafikov 2002 (Eq. 13 Bollati et al. 2021) This is a vectorised version of _t. _t computes an integral where the integrand is independent of the radius r. The radius r only changes the end point of the integral. Instead of using quad for each end point (which is not as easily vectorisable), we can utilise an ODE solver to obtain the result at every end point without doing redundant work.