diskpy.ICgen package

Submodules

diskpy.ICgen.AddBinary module

@dflemin3 David Fleming 2015

This module is useful when considering a two body system in a Keplerian potential. Uses include binary stars, a satellite orbit a massive primary (m1 >> m2), a satellite orbiting far from a massive central binary, etc. Functions included in this module can be used to compute orbital elements, convert between Cartersian coordinates and orbital elements and other uses.

Note on units: Functions are optimized to take pynbody.array.SimArray as input types so that they may handle all the different, weird units and array types/shapes. If you don’t use SimArrays, they should work if your inputs are in cgs, but use at your own peril.

diskpy.ICgen.AddBinary.aToP(a=1, M=1)[source]

Given a semimajor axis (au), convert into period (in days) given Kepler law

Parameters:

a : float

semimajor axis [au]

M : float

mass of system [Msol]
Returns:

Period : float

[days]
diskpy.ICgen.AddBinary.accretionEDot(Binary, M_dot, dt)[source]

Given a Binary object and an accretion rate (assumed to be in M_sol/yr), compute the theoretical rate of change of the binary’s eccentricity in 1/second. Assumptions: -radius, velocity of binary nearly constant over accretion (found to more or less apply via empirical measurements)

Parameters:

Binary: binary object class

Mdot: float

accretion rate in M_sol/yr

dt: float

time interval
Returns:

de/dt : float

change in eccentricity in 1/second
diskpy.ICgen.AddBinary.binaryPrecession(s, r_in, r_out)[source]

Computes the period of the binary precession due to an axisymmetric disk given by Rafikov 2013.

Parameters:

s : Tipsy snapshot

r_in, r_out : float

inner and outer radii of the circumbinary disk [AU]
Returns:

T : SimArray

Period of binary argument of periastron precession in yr
diskpy.ICgen.AddBinary.calcArgPeri(x1, x2, v1, v2, m1, m2, flag=True)[source]

Given pynbody arrays for positions and velocity of primary and secondary bodies and masses, calculates the eccentricity vector in the reduced two body system. Usage note: Intended for binary system, but pass x2 = v2 = 0 to use with any location in the disk.

Parameters:

Assumed as pynbody SimArrays [preferred units]

x1,x2 : SimArrays

Primary and secondary position arrays [AU]

v1,v2 : SimArrays

Primary and secondary velocity arrays [km/s]

m1, m2 : SimArrays

Primary and secondary masses [Msol]

Flag : bool

Whether or not to internally convert to cgs units If flag == false, input assumed to be in cgs
Returns:

w : float

Argument of pericenter in degrees
diskpy.ICgen.AddBinary.calcCOM(m1=0.5, m2=0.5, x1=1, x2=1)[source]

Given pynbody arrays for the mass and position of the two binary stars, function calculates the CoM of the two particles in order to check that it’s still roughly 0.

Deprecated

Parameters:

(as pynbody SimArrays)

m1, m2 : SimArrays

Primary and secondary mass arrays (in Msol)

x1, x2 : SimArrays

Primary and secondary position arrays (in AU)
Returns:

center of mass: array

Center of mass position vector (numpy array in AU)
diskpy.ICgen.AddBinary.calcCircularFrequency(x1, x2, v1, v2, m1, m2, flag=True)[source]

Given pynbody arrays for positions and velocity of primary and secondary bodies and masses, calculates the circular frequency in the reduced two body system. Usage note: Intended for binary system, but pass x2 = v2 = 0 to use with any location in the disk.

omega = (L_z)/(R^2) which assumes spherical symmetry (ok assumption here) L = sqrt(G*M*a*(1-e^2) for ~Keplerian

Deprecated (just use sqrt(G*M/a^3))…why did I make this?

Parameters:

Assumed as pynbody SimArrays in simulation units (AU, scaled velocity, etc)

x1,x2 : arrays

Primary and secondary position arrays (in AU)

v1, v2 : arrays

Primary and secondary velocity arrays (km/s)

m1, m2 : floats

Primary and secondary masses (Msol)

Flag : Whether or not to internally convert to cgs units
Returns:

omega: float

Circular frequency in 1/days
diskpy.ICgen.AddBinary.calcCriticalRadius(a=1, e=0, m1=0.5, m2=0.5)[source]

Calculates the approximate bounds for where we would expect a planet to form given the conditions of a circumbinary disk around a short-period binary system. Calculates based on best fit from Holman&Wiegert+1999 (outer/P-type region) Assumes m1 ~ m2 and NOT m1 >> m2 or m2 >> m1

Parameters:

a : float

Semimajor axis a of the binary system (au)

e : float

Eccentricity e of binary system

m1, m2 : floats

Masses of binary components m1, m2 (Msol)
Returns:

ac, pmac : floats

Lower bounds for circumbinary planet’s distace from binary COM and error terms of the following form: ac, pmac (error bounds are symmetric) both in au
diskpy.ICgen.AddBinary.calcEcc(x1, x2, v1, v2, m1, m2, flag=True)[source]

Given as pynbody arrays the masses of the binary system, arrays of the positions and velocities, compute its orbital eccentricity.

Calculates e using following: e = sqrt(1 + (2*e*h^2)/(mu^2) for h = r x v, mu = G*(m1+m2), e = (v^2)/2 - (mu/|r|)

Parameters:

All inputs expected to be pynbody simArrays

x1, x2 : SimArrays

Position arrays of primary and secondary x1, x2 (in AU)

v1, v2 : SimArrays

Velocity arrays of primary and secondary v1, v2 (in km/s)

m1, m2 : SimArrays

Masses of primary and secondary m1, m2 (in Msol)

Flag: bool

Whether or not to internally convert to cgs units
Returns:

e : float

Scalar eccentricity of binary system.
diskpy.ICgen.AddBinary.calcEccVector(x1, x2, v1, v2, m1, m2, flag=True)[source]

Given pynbody arrays for positions and velocity of primary and secondary bodies and masses, calculates the eccentricity vector in the reduced two body system. Usage note: Intended for binary system, but pass x2 = v2 = 0 to use with any location in the disk.

Parameters:

x1,x2 : SimArrays

Primary and secondary position arrays [length]

v1, v2 : SimArrays

Primary and secondary velocity arrays [velocity]

m1,m2 : SimArrays

Primary and secondary masses [mass]

Flag : bool

Whether or not to internally convert to cgs units
Returns:

Ecc : SimArray

Eccentricity vector in cgs
diskpy.ICgen.AddBinary.calcEccentricAnomaly(x1, x2, v1, v2, m1, m2, flag=True)[source]

Given pynbody arrays for positions and velocity of primary and secondary bodies and masses, calculates the eccentric anomaly in the reduced two body system. Usage note: Intended for binary system, but pass x2 = v2 = 0 to use with any location in the disk.

Parameters:

Assumed as pynbody SimArrays [preferred units]

x1,x2 : SimArrays

Primary and secondary position arrays [AU]

v1,v2 : SimArrays

Primary and secondary velocity arrays [km/s]

m1, m2 : SimArrays

Primary and secondary masses [Msol]

Flag : bool

Whether or not to internally convert to cgs units if flag is false, input assumed to be in cgs

Returns

——-

E : float

Eccentric anomaly in degrees
diskpy.ICgen.AddBinary.calcInc(x1=1, x2=0, v1=1, v2=0, flag=True)[source]

Given pynbody arrays for positions and velocities of primary and secondaries bodies and masses in a binary orbit, calculate’s the orbit’s inclination. (Given: Orbit starts in xy plane)

i = arccos(h_z/|h|)

Parameters:

as pynbody SimArrays [preferred units]

x1,x2 : SimArrays

Primary and secondary position arrays [AU]

v1, v2 : SimArrays

Primary and secondary velocity arrays [km/s]

Flag : bool

Whether or not to internally convert to cgs units
Returns:

i : float

Inclination in degrees
diskpy.ICgen.AddBinary.calcLongOfAscNode(x1, x2, v1, v2, flag=True)[source]

Given pynbody arrays for positions and velocity of primary and secondary bodies and masses, calculates the binary’s longitude of the ascending node Omega. Usage note: Intended for binary system, but pass x2 = v2 = 0 to use with any location in the disk.

Calculates Omega using the following: Omega = arccos(n_x/|n|) n_y > 0 Omega = 2*pi - arccos(n_x/|n|) n_y < 0 where n = (0,0,1) x h for h = r x v

Parameters:

Assumed as pynbody SimArrays in simulation units (AU, scaled velocity, etc)

x1,x2 : SimArrays

Primary and secondary position arrays [AU]

v1,v2 : SimArrays

Primary and secondary velocity arrays [km/s]

Flag : bool

whether or not to convert to cgs. if false, already in cgs
Returns:

Omega : SimArray

longitude of the ascending node in degrees
diskpy.ICgen.AddBinary.calcMeanAnomaly(x1=1, x2=0, v1=1, v2=0, m1=1, m2=1, flag=True)[source]

Given pynbody arrays for positions and velocity of primary and secondary bodies and masses, calculates the Mean anomaly in the reduced two body system. Usage note: Intended for binary system, but pass x2 = v2 = 0 to use with any location in the disk.

Parameters:

Assumed as pynbody SimArrays [preferred units]

x1,x2 : SimArrays

Primary and secondary position arrays [AU]

v1,v2 : SimArrays

Primary and secondary velocity arrays [km/s]

m1, m2 : SimArrays

Primary and secondary masses [Msol]

Flag : bool

Whether or not to internally convert to cgs units
Returns:

M : float

Mean anomaly in degrees
diskpy.ICgen.AddBinary.calcMeanMotion(x1, x2, v1, v2, m1, m2, flag=True)[source]

Given pynbody arrays for positions and velocity of primary and secondary bodies and masses, calculates the Mean motion in the reduced two body system. Usage note: Intended for binary system, but pass x2 = v2 = 0 to use with any location in the disk.

Parameters:

Assumed as pynbody SimArrays [preferred units]

x1,x2 : SimArrays

Primary and secondary position arrays [AU]

v1,v2 : SimArrays

Primary and secondary velocity arrays [km/s]

m1, m2 : SimArrays

Primary and secondary masses [Msol]

Flag : bool

Whether or not to internally convert to cgs units
Returns:

n : SimArray

Mean motion in 1/s
diskpy.ICgen.AddBinary.calcOrbitalElements(x1, x2, v1, v2, m1, m2)[source]

Given as pynbody SimArrays the cental mass(es), the coodinate(s) and velocity(ies) of a CCW orbiting object, return the following orbital elements: eccentricity, semimajor axis, inclination, longitude of ascending node, argument of periapsis, and true anomaly. This function is designed to work for a binary star system but is general enough to also work for a ~massless gas particle orbiting around a central quasi-Keplerian mass.

Note: all params not required to be pynbody.array.SimArray, but it is strongly preferred. Returns floats to hopefully make this as general as possible

Parameters:

x1, x2 : SimArrays

Position arrays of primary and secondary x1, x2 (in AU)

v1, v2 : SimArrays

Velocity arrays of primary and secondary v1, v2 (in km/s)

m1, m2 : SimArrays

Masses of primary and secondary m1, m2 (in Msol)
Returns:

e : float

Eccentricity (unitless)

a : float

Semimajor Axis in Au

i : float

Inclination in degrees

Omega : float

Longitude of Ascending node in degrees

w : float

Argument of Periapsis in degrees

nu : float

True Anomaly in degrees
diskpy.ICgen.AddBinary.calcPositions(M=1, a=1, e=0, p=0.5)[source]

Given total mass of system M, semimajor axis a, and percentage of total mass contained in primary star p, calculate positions of binary components keep COM at origin (ex: If 1Msol system and MPri = Msec = 0.5 -> M =1 -> p = 0.5) Assume stars start a perihelion

Deprecated

Parameters:

M : float

Total mass of system (Msol)

a : float

Semimajor axis (au)

e : float

eccentricity

p : float

% of total mass contained in primary (m1)
Returns:

x1, x2: float

Semimajor axes of binary stars assuming that they start at perihelion.
diskpy.ICgen.AddBinary.calcRocheLobe(q, a)[source]

Given the mass ratio q = m1/m2 and the semimajor axis in AU, compute the radius of the Roche lobe around m1 given the Eggleton approximation http://en.wikipedia.org/wiki/Roche_lobe

Parameters:

q : float

Binary mass ratio q = m1/m2

a : float

Binary semimajor axis a (in AU)
Returns:

r : float

Radius of Roche lobe around m1 (in AU)
diskpy.ICgen.AddBinary.calcSemi(x1, x2, v1, v2, m1, m2, flag=True)[source]

Given pynbody arrays for positions and velocity of primary and secondary bodies and masses, calculates the binary’s semimajor axis.

Calculates a using the following: a = -mu/(2*e) where mu = G*(m1+m2) and e = (v^2)/2 - (mu/|r|) and e != eccentricity

Parameters:

x1,x2 : SimArrays

Primary and secondary position arrays [AU]

v1,v2 : SimArrays

Primary and secondary velocity arrays [km/s]

m1,m2 : SimArrays

Primary and secondary masses [Msol]

Flag : bool

Whether or not to internally convert to cgs units
Returns:

a : float

semimajor axis of binary orbit in AU
diskpy.ICgen.AddBinary.calcTrueAnomaly(x1, x2, v1, v2, m1, m2, flag=True)[source]

Given pynbody arrays for positions and velocity of primary and secondary bodies and masses, calculates the true anomaly in the reduced two body system. Usage note: Intended for binary system, but pass x2 = v2 = 0 to use with any location in the disk.

Parameters:

Assumed as pynbody SimArrays [preferred units]

x1,x2 : SimArrays

Primary and secondary position arrays [AU]

v1,v2 : SimArrays

Primary and secondary velocity arrays [km/s]

m1, m2 : SimArrays

Primary and secondary masses [Msol]

Flag: bool

Whether or not to internally convert to cgs units if flag is false, input assumed to be in cgs
Returns:

nu: float

True anomaly in degrees
diskpy.ICgen.AddBinary.calcV(m1=0.5, m2=0.5, a=1, e=0)[source]

Given total mass M, postions of stars at perihelion x1, x2, and eccentricity e, calculate the velocities of the stars assuming that they are located at the perihelion and rotate in same direction as disk (CCW)

Deprecated

Parameters:

m1, m2 : floats

masses of primary and secondary (Msol)

x1, x2 : floats

are semimajor axes of primary and secondary (au)

e : float

eccentricity
Returns:

v1, v2 : floats

velocities of m1, m2 in km/s oriented for CCW rotation (in xy plane)
diskpy.ICgen.AddBinary.computeLenAx(a)[source]

Given an numpy array, check to see if it’s length is 1 (aka it’s a float, or int or whatever) or otherwise and return it’s length and the axis over which calculations like normalization and cross product is to be taken. This function is useful when you expect data arrays of the form x = [n,3] with weird combinations of lists/numpy arrays

diskpy.ICgen.AddBinary.compute_E(a, m1, m2)[source]

Given an objects semimajor axis a and mass around which it is orbiting, M compute it’s total energy assuming a quasi-Keplerian orbit

E = -0.5*G*m1*m2/a

Parameters:

a : SimArray

semimajor axis

m1 : SimArray

mass of (primary) central body

m2 : SimArray

mass of secondary body
Returns:

E : SimArray

Keplerian total energy in cgs
diskpy.ICgen.AddBinary.compute_L(a, e, m1, m2)[source]

Given an objects semimajor axis a, eccentricity e, and mass around which it is orbiting, M compute it’s angular momentum assuming a quasi-Keplerian orbit

L = mu*sqrt(G*M*a*(1-e^2))

Parameters:

a : SimArray

semimajor axis

e : float (or array of float or SimArray with units = “1”)

eccentricity

m1 : SimArray

mass of (primary) central body

m2 : SimArray

mass of secondary body
Returns:

L : SimArray

Keplerian angular momentum in cgs
diskpy.ICgen.AddBinary.dotProduct(a, b)[source]

Given n x m numpy arrays, compute the dot product row-wise.

Parameters:

a,b : n x m numpy array
Returns:

dot product : numpy array

n rows containing dot products
diskpy.ICgen.AddBinary.initializeBinary(a, e, i, Omega, w, M, m1, m2, angleFlag=True, scaleFlag=True)[source]

Given the initial Kepler orbital parameters, compute the Cartesian positions and velocities for 2 mutually orbiting stars (binaries!).

Parameters:

Keplerian orbital elements : floats

in Au, degrees (see above for more details)

m1,m2 : floats

Masses of central objects (Msol)

angleFlag : bool

whether or not to convert from degrees->radians (True = convert)

scaleFlag : bool

whether or not to put v in sim units (True = do it/default option)
Returns:

x1, x2 : array

Positions of 2 objects about origin (Au)

v1, v2 : array

Velocities of 2 objects about origin (km/s in Sim Units)
diskpy.ICgen.AddBinary.keplerToCartesian(a, e, i, Omega, w, M, m1, m2, angleFlag=True, scaleFlag=True)[source]

TODO: make SimArray compatible Given the Keplerian orbital elements, compute the cartesian coordinates of the object orbiting in the reduced mass frame. Note: Requires all angles in degrees unless noted.

Note: A little redudant that I compute M when I typically already know the true anomaly nu, but most other schemes know M initially instead of nu so I’ll keep it for compatibility’s sake.

Parameters:

a : float

Semimajor axis (AU)

e : float

Eccentricity

i : float

inclination (degrees)

Omega : float

Longitude of Ascending Node (degrees)

w : float

Argument of Pericenter (degrees)

M : float

Mean Anomaly (degrees)

m1, m2 : float

Masses of central object(s) (Msol)

angleFlag, scaleFlag : bools

Tells code to convert degrees->rad and/or scale velocity to sim units, respectively
Returns:

x : numpy array

Position array of the object in reduced mass frame (AU)

v : numpy array

velocity array (km/s) with VEL_UNIT scaling factor optional
diskpy.ICgen.AddBinary.pToA(period=1, M=1)[source]

Converts period (in days) into semimajor axis (in au) given Kepler law

Parameters:

Period : float

[days]

M : float

COM mass [Msol]
Returns:

Semimajor axis a (au)
diskpy.ICgen.AddBinary.reduceToPhysical(r, v, m1, m2)[source]

Function converts from reduced mass coordinates to physical, origin-centered coords Works for arbitrary units, number of objects with coordinates as long as m1+m2 = central mass.

Parameters:

r : array

radius vector

v : array

velocity vector

m1, m2 : floats

mass of central object
Returns:

x1, x2 : arrays

position vectors of 2 mutually orbiting objects

v1, v2 : arrays

velocity vectors of 2 mutually orbiting objects
diskpy.ICgen.AddBinary.trueToMean(nu, e)[source]

Given the true anomaly nu in degrees and the eccentricity e, compute the mean anomaly M in degrees. Can take arbitary units.

Parameters:

nu : float

True anomaly (degrees)

e : float

eccentricity
Returns:

M: float

mean anomaly (degrees)

diskpy.ICgen.ICgen module

Created on Wed Mar 12 12:48:33 2014

@author: ibackus

class diskpy.ICgen.ICgen.IC(r=None, sigma=None, CDF=None, profile_kind=None, settings=None)[source]

Defines the IC class.

INITIALIZING NEW INITIAL CONDITIONS

# Initialize a blank IC object (with no surface density profile yet made)

IC = ICgen.IC()

# Initialize an IC object and surface density profile using default settings

IC = ICgen.IC(profile_kind=’powerlaw’) IC = ICgen.IC(profile_kind=’MQWS’)

# Initialize IC object from 1-D SimArrays r, sigma (surface density)

IC = ICgen.IC(r, sigma)

Optionally, the CDF for the surface density profile can be supplied to speed up generation of sigma. To do that:

IC = ICgen.IC(r, sigma, CDF)

Or, the input can be a the filename of a pickled dictionary containing ‘r’, ‘sigma’, and optionally ‘CDF’

Settings can also be entered manually if needed

settings = pickle.load(open(‘settings.p’,’r’)) IC = ICgen.IC(settings = settings)

GENERATING INITIAL CONDITIONS

Qest(r=None)[source]

Estimate Toomre Q at r (optional) for ICs, assuming omega=epicyclic frequency. Ignores disk self-gravity

generate(restart=False)[source]

Runs through all the steps to generate a set of initial conditions

IF restart=True, it picks up at the last completed step

diskpy.ICgen.ICgen.Qest(ICobj, r=None)[source]

Estimate Toomre Q at r (optional) for ICs, assuming omega=epicyclic frequency. Ignores disk self-gravity

class diskpy.ICgen.ICgen.add(ICobj)[source]

Contains modules to load data/parameters

rho(rho_dict)[source]

Generates a rho object and stores it in ICobj.rho

rho_dict should be a dictionary containing:
‘z’: 1D array of z values ‘r’: 1D array of r values ‘rho’: 2D array of rho evaluated at z,r

Exaple:

rho_dict = pickle.load(open(‘rhofile.p’, ‘rb’)) # Load up a rho dict ICobj.add.rho(rho_dict) # create ICobj.rho

diskpy.ICgen.ICgen.load(filename)[source]
class diskpy.ICgen.ICgen.maker(ICobj)[source]

A Wrapper containing various functions for generating initial conditions. Outputs of the functions are saved to the IC object. The IC object is referenced as self._parent. So to access temperature, simply call self._parent.T(r)

pos_gen(method=None)[source]

A wrapper for generating positions according to rho and sigma

Initializes a pos object (see pos_class.py) and saves it to ICobj.pos

IF called with method not set, the method used is:
ICobj.settings.pos_gen.method
rho_gen(**kwargs)[source]

A wrapper for calc_rho_zr.

Upon executing, generates rho and rho cdf inverse

Optional kwargs can be supplied and passed to the rootfinder used (scipy.optimize.newton_krylov)

sigma_gen(r=None, sigma=None, CDF=None)[source]

A Wrapper for make_sigma.sigma_gen See make_sigma.sigma_gen for documentation

Upon executing, generates sigma, pdf, and cdf_inv and saves to ICobj

USAGE:

# Generate sigma object from r, sigma arrays sigma = ICobj.maker.sigma_gen(r, sigma)

# Use pre-calculated CDF array sigma = ICobj.maker.sigma_gen(r, sigma, CDF)

# Generate using a profile defined in sigma_profile.py ICobj.settings.sigma.kind = ‘powerlaw’ sigma = ICobj.maker.sigma_gen()

r and sigma should be 1-D SimArrays. sigma is the surface density evaluated at r

snapshot_gen()[source]

A wrapper for generating a tipsy snapshot from the initial conditions

Uses make_snapshot.py

diskpy.ICgen.ICgen.save(ICobj, filename=None)[source]

diskpy.ICgen.ICgen_settings module

Defines the settings class for ICgen

USAGE:

Load settings from file:

import ICgen_settings settings = ICgen_settings.settings(‘filename’)

Load settings from defaults (Defined in ICgen_settings.py):

import ICgen_settings settings = ICgen_settings.settings()

AFTER LOADING: Access settings:

# Return number of particles: settings.pos_gen.nParticles

# Etc…

Echo settings:

# Echo ALL settings settings()

# Echo filenames settings.filenames()

# Etc…

Echo setting defaults (with descriptions)
settings.print_defaults() settings.sigma.print_defaults()

READ/WRITE

Write (save) settings to disk:

settings.write(‘filename’)

Load (read) settings from disk:

settings.load(‘filename’)

A note on implementation. There are 2 base settings classes: settingsBase and kindSettingsBase. The first is just for generic settings, the second is to allow users to have multiple layouts (see class sigma)

See the base classes for more info

class diskpy.ICgen.ICgen_settings.changa_run(**kwargs)[source]

Bases: diskpy.ICgen.ICgen_settings.settingsBase

Settings for running ChaNGa (needed when calculating velocities, eps, so on)

diskpy.ICgen.ICgen_settings.default_settings()[source]

Echos the current default settings with descriptions for ICgen settings

class diskpy.ICgen.ICgen_settings.filenames[source]

Bases: diskpy.ICgen.ICgen_settings.settingsBase

Filenames

To echo settings:
filenames()
To access settings:
filenames.settingname
class diskpy.ICgen.ICgen_settings.kindSettingsBase(descr=None)[source]

Bases: diskpy.ICgen.ICgen_settings.settingsBase

print_defaults()[source]
class diskpy.ICgen.ICgen_settings.physical(**kwargs)[source]

Bases: diskpy.ICgen.ICgen_settings.settingsBase

Defines default physical parameters

To echo settings:
physical()
To access settings:
physical.settingname
class diskpy.ICgen.ICgen_settings.pos_gen[source]

Bases: diskpy.ICgen.ICgen_settings.settingsBase

Settings for generating random positions [STEP 2]

To echo settings:
pos_gen()
To access settings:
pos_gen.settingname
diskpy.ICgen.ICgen_settings.print_settings(setting)[source]
diskpy.ICgen.ICgen_settings.repr_defaults(x, indent=0)[source]
Prints a defaults list. Defaults lists are a list of length 3 tuples:
[ (‘variable_name’, default_val, ‘Doc string for variable’),
(‘variable_name2’, default_val2, ‘Other doc string’) ]
diskpy.ICgen.ICgen_settings.repr_settings(setting)[source]
class diskpy.ICgen.ICgen_settings.rho_calc[source]

Bases: diskpy.ICgen.ICgen_settings.settingsBase

Settings for calculating rho(z,r)

To echo settings:
rho_calc()
To access settings:
rho_calc.settingname
class diskpy.ICgen.ICgen_settings.settings(settings_filename=None, kind=None)[source]

Bases: diskpy.ICgen.ICgen_settings.settingsBase

settings for ICgen

USAGE:

Load settings from file:

import ICgen_settings settings = ICgen_settings.settings(‘filename’)

Load settings from defaults (Defined in ICgen_settings.py):

import ICgen_settings settings = ICgen_settings.settings()

Load settings for a given surface density profile

import ICgen_settings settings = ICgen_settings.settings(‘powerlaw’) or try: settings = ICgen_settings.settings(‘MQWS’)

AFTER LOADING: Access settings:

# Return sigmaFileName: settings.filenames.sigmaFileName

# Return number of particles: settings.pos_gen.nParticles

# Etc…

Echo settings:

# Echo ALL settings settings()

# Echo filenames settings.filenames()

# Etc…

READ/WRITE

Write (save) settings to disk:

settings.write(‘filename’)

Load (read) settings from disk:

settings.load(‘filename’)
load(settings_filename)[source]

Load settings from settings_filename.

save(settings_filename=None)[source]

Save settings to settings_filename. If settings_filename is None (or is not set), settings_filename will be taken as the filename of the most recent save/load

class diskpy.ICgen.ICgen_settings.settingsBase(descr=None, **kwargs)[source]
print_defaults()[source]

Prints defaults defined in self._defaults. Defaults lists are a list of length 3 tuples: [ (‘variable_name’, default_val, ‘Doc string for variable’),

(‘variable_name2’, default_val2, ‘Other doc string’) ]

Also prints childrens defaults if they are settingsBase instances

class diskpy.ICgen.ICgen_settings.sigma(kind=’powerlaw’)[source]

Bases: diskpy.ICgen.ICgen_settings.kindSettingsBase

Settings for generating a surface density profile

To echo:
sigma()
To access settings:
sigma.settingname
class diskpy.ICgen.ICgen_settings.snapshot(nParticles=0)[source]

Bases: diskpy.ICgen.ICgen_settings.settingsBase

Settings for generating tipsy files (includes particle mass, temp, and vel.) [STEP 3]

diskpy.ICgen.ICgen_utils module

Created on Tue Aug 5 17:37:03 2014

@author: ibackus

diskpy.ICgen.ICgen_utils.Qeff(ICobj, bins=None)[source]
diskpy.ICgen.ICgen_utils.checkversion(obj)[source]

Checks the version of an object. Returns -1 if there is no version

diskpy.ICgen.ICgen_utils.delistify(in_list)[source]

Reverses listify. Takes an larray object (a listified array) and returns the original array

ARGUMENTS

in_list : larray
Listified array

RETURNS

array : array_like or SimArray
array from the original data

See Also : ICgen_utils.listify

diskpy.ICgen.ICgen_utils.est_eps(smoothlength_file, nstar=1)[source]

Estimate gravitational softening length (eps) from a ChaNGa output .smoothlength file. eps is estimated as 1/2 the mean smoothing length

ARGUENTS

smoothlength_file : str
Filename of the .smoothlength file
nstar : int
Number of star particles present

RETURNS

eps : float
Estimate of the gravitational softening length in simulation units
class diskpy.ICgen.ICgen_utils.larray(shape=None)[source]

Bases: list

A simple subclass of list for containing listified arrays (see ICgen_utils.listify) Should be instantiated by ICgen_utils.listify

USAGE:

Creating an larray object:

a = larray() # blank larray a = larray(shape) # sets a.shape = shape

Then one can append to an larray just as with a list

a.append(stuff)

To return return a normal array:

array = a.delistify()
delistify()[source]

Return an array made from self. See ICgen_utils.delistify

diskpy.ICgen.ICgen_utils.listify(array, max_element=10000000)[source]

Breaks up an array or SimArray into chunks and saves as an larray object (essentially, a list). Useful for pickling very large arrays which otherwise may throw an error with pickle

Whenever possible the array is NOT copied, rather a view is returned. This depends on the functionality of array.ravel() (see numpy.ravel)

ARGUMENTS

array : array_like or SimArray
Input array to listify
max_element : int
Maximimum number of elements per chunk (i.e., per list item)

RETURNS

list_array : larray
A listified array.

See Also : ICgen_utils.delistify

diskpy.ICgen.ICgen_utils.pickle_import(fname, moduledir=None)[source]

Performs a cPickle.load on fname If an ImportError is raised, moduledir is temporarily added to sys.path to try and find the module

Parameters:

fname : filename

filename of a file to open

moduledir : str or list of strings

The directory or directories to temporarily add to sys.path IF None, a normal cPickle.load is performed

diskpy.ICgen.binary module

author: @dflemin3 July 2015 Module for Binary star class. Holds the Cartesian position/velocity coordinates and Kepler orbital elements for the binary in the reduced mass frame. Can be initialized with orbital elements (preferred) or Cartesian position/velocity or a tipsy format snapshot that was read in using pynbody.load(“snapshot.name”). Just need to pass a string to let it know which input type it’s getting.

class diskpy.ICgen.binary.Binary(X, m1, m2, state)[source]

Bases: object

Defines the binary star class Binary. Binary star class used for storing positions/velocities and Keplerian Orbital elements for the system.

Input Units: Lengths = AU, Velocities = km/s (*29.8 == VEL_UNIT), angles in degrees Output Units: Same as input

Initializing:

With Cartesian:
Binary(X,m1,m2,”Cartesian”) where X is a list or numpy array of the form [x,y,z,vx,vy,vz] and m1, m2 are the masses of the primary and secondary (M_sol), respectively
With Kepler Orbital Elements:
Binary(X,m1,m2,”Kepler”) where X is e, a, i, Omega, w, nu which are the eccentricity, semimajor axis (AU), inclination (degrees), Longitude of Ascending Node (degrees), argument of periapsis (degrees), and the true anomaly (degrees). Note the units! m1, m2 are defined above.
With Snapshot:
X is a tipsy snapshot to be loaded. Positions and orbital elements will be calculated from it.
assignOrbElems(X)[source]

Given array of orbital elements, set them as class parameters. Intended as a shorthand function.

Input: X is e, a, i, Omega, w, nu which are the eccentricity, semimajor axis (AU), inclination (degrees), Longitude of Ascending Node (degrees), argument of periapsis (degrees), and the true anomaly (degrees). Convert all values to float for consistency and as a sanity check except for semimajor axis since that gets units we’re interested in, au.

Output: None

computeCartesian()[source]

Compute the Cartesian position and velocity in the reduced mass frame.

computeOrbElems()[source]

Compute the Kepler orbital elements.

Input: Self (must have r, v set and in proper units)

Output: sets and returns e,a,i,…

generateICs()[source]

From Kepler orbital elements, compute the position, velocities for two stars in ChaNGa-friendly units. Called “generateICs” since I’ll use this mainly to…generate initial conditions

reshapeData()[source]

Ensure data, specifically r and v arrays, are of the proper shape for further manipulation. Also ensure data has proper units This is useful because sometimes they come in as a list, a (1,3) numpy array or a (3,) numpy array. Much easier to clean up upon initialization then have many checks in later functions.

diskpy.ICgen.binaryUtils module

David Fleming Utilities to process/interact with binary star system in ChaNGa sims

Note on inputs: Most (if not all?) functions are designed to be used with SimArrays as inputs.

diskpy.ICgen.binaryUtils.calcCoMVsRadius(s, rBinEdges, starFlag=False)[source]

Calculates the system’s center of mass as a function of radius. At a given radius r, use the total enclosed mass of the star(s) and gas to compute the center of mass (CoM).

Parameters:

s: Tipsy-format snapshot readable by pynbody

rBinEdges: array

edges of the array of radii in xy plane

starFlag: bool

whether or not to consider stars in center of mass calculation
Returns:

com: array

Numpy array of len(r) * 3 containing location of CoM in Cartesian coordinates.
diskpy.ICgen.binaryUtils.calcDeDt(stars, tau)[source]

Calculates the change in binary orbital eccentricity over time at each radial bin due to the torque/mass from the surrounding CB disk.

Parameters:

stars: pynbody stars object (sim units)

tau: torque/mass on binary due to CB disk during a given snapshot (cgs)

tau is of the form tau[radii,(x,y,z) components]

Output:

(de/dt) at each radius (unitless/second)
diskpy.ICgen.binaryUtils.calcDiskRadialBins(s, r_in=0, r_out=0, bins=50)[source]

Cleanly partitions disk into radial bins and returns the bin edges and central bin values. Note, default ndim = 2 so all bins are in 2D plane (i.e. radius r is polar/cylindrical radius in xy plane which makes sense for thin disks)

Parameters:

s: Pynbody snapshot

r_in: float

Inner disk radius you’ll consider (AU)

r_out: float

Outer disk radius you’ll consider (AU)

bins: int

# of bins
Returns:

r: numpy array

central radial bin values (AU)

rBinEdges: numpy array

edges of radial bins
diskpy.ICgen.binaryUtils.calcNetTorque(stars, gas)[source]

Given pynbody snapshot (Tipsy format)arrays of the stars and gas of a binary surrounded by a CB disk, compute the net torque on the binary due to the CB disk. This function can be used to compute the net torque/mass due to any collection of gas (total disk, an annulus, etc) on the stars.

Parameters:

stars, gas: pynbody-readable Tipsy snapshot arrays

of binary + CB disk. Assumes units are in standard sim units (Msol,au…)
Returns:

net torque: numpy array

Net torque/mass vector (3D) acting on binary system in cgs.
diskpy.ICgen.binaryUtils.calcPoissonVsRadius(s, rBinEdges)[source]

Given a tipsy snapshot and radial bins, compute the Poisson noise, r/sqrt(N_particles), in each radial bin.

diskpy.ICgen.binaryUtils.calcQ(cs, kappa, sigma)[source]

Compute the Toomre Q parameter for a gaseous disk. Implimented here since pynbody calculates it for a stellar disk. Q = (c_s * kappa)/(pi*G*Sigma) > 1 -> axisymmetric stability

Parameters:

c_s: float

sound speed (cm/s)

Kappa: float

radially epicycle frequency (1/s). Can also be Omega for a quasi-Keplerian disk

Sigma: float

surface density at some radius (g/cm^2)

Output:

Q: unitless

Toomre parameter
diskpy.ICgen.binaryUtils.calcQVsRadius(s, a_in, a_out, bins)[source]

Given a tispy snapshot, compute the Toomre Q parameter at each radial point.

Parameters:

s: Tipsy snapshot

a_in, a_out: floats

Minimum and maximum radial points on which to calculate the profile [AU]

bins: int

number of radial bins
Returns:

r: array, AU

Q: array, unitless
diskpy.ICgen.binaryUtils.calcStableSigma(r, rd, Mstar, Mdisk, Q)[source]

Compute the surfance density sigma_0 such that if the quasi-keplerian disk had a surface density of sigma > sigma_0 at an OLR, the disk would be unstable to a m=1 mode. Condition comes from eqn. 110 in Shu 1990. Note: Assumes (b_n - c)^2 ~ 1 as authors did.

Parameters:

r: array

radii of OLR [AU]

rd: float

Maximum radially extent of the disk [AU]

Mstar, Mdisk: floats

Masses of the central star(s) and disk, respectively [Msol]

Q: float

Toomre Q stability parameter evalutated at rd
Returns:

sigma_0: array

critical surfance density [Msol/AU^2]
diskpy.ICgen.binaryUtils.calc_LB_resonance(s, m_min=1, m_max=3, l_min=1, l_max=3)[source]

Computes the locations of various Lindblad Resonances in the disk as a function of binary pattern speed.

Parameters:

s : Tipsy-format snapshot

m_min, l_min : ints

minimum orders of (m,l) LR

m_max,l_max : ints

maximum orders of (m,l) LR
Returns:

OLR, ILR, CR: numpy arrays

location in AU of (m,l)th order Lindblad resonances
diskpy.ICgen.binaryUtils.changaFloatSearch(name, simUnits=False)[source]

Given the name of a file containing line dumps for ChaNGa and outputs numpy arrays containing changa dumps line-by-line.

Genfromtxt will convert output into either an numpy array of strings or a numpy array of lists of strings. Determine what format is, use regex library to find all floats and output accordingly.

Default usage is searching for linear momentum dumps of the form mg,vx,vy,vz for gas Assume appropriate flag was used to grep data into input file

Input: name: Name of input file (something.txt) simUnits: whether or not to use simUnits (False -> convert to cgs)

Only for use for linear momentum values!

Output: Numpy array containing all floats from changa output

diskpy.ICgen.binaryUtils.computeCOM(stars, gas, cutoff=None, starFlag=True, gasFlag=True)[source]

Given pynbody star and gas arrays, compute the center of mass for the entire specified system.

Parameters:

stars: s.stars pynbody object

gas: s.gas pynbody object

cutoff: float

radius at which you only consider objects interior to it

starFlag: bool

whether or not to consider stars
diskpy.ICgen.binaryUtils.computeVelocityCOM(s, cutoff=None, starFlag=True, gasFlag=True)[source]

Given pynbody star and gas arrays, compute the center of mass velocity for the entire specified system.

Parameters:

s : pynbody snapshot

cutoff : float

radius at which you only consider objects interior to it [AU]

starFlag : bool

whether or not to consider stars

gasFlag : bool

whether or not to consider gas
Returns:

Center of mass velocity: SimArry

in AU for each vx,vy,vz component

Note: a lot of “strip_units” commands included to prevent throwing weird value errors. As long as all masses

are in solar masses and positions in AU before this is run, you won’t have any problems.
diskpy.ICgen.binaryUtils.diskAverage(s, r_out, bins=50, avgFlag=True)[source]

Computes the accretion disk mass-averaged for x via the following equation: integral of 2*pisigma*x*r*dr / integral of 2*pi*sigma*r*dr. Sigma, e,a… calculated on the fly to ensure that they are all evaluated at the same location.

Parameters:

s : tipsy snapshot

r_out : float

outer radii for averaging region. If none, use entire disk

bins : int

how many radial bins to calculate quantities over

avgFlag : bool

whether or not to average over all particles in a radial bin for orbital element calculation
Returns:

y : list

disk-averaged Keplerian orbital elements [e,a,i,Omega,w,nu] in AU, degrees (depending on unit)
diskpy.ICgen.binaryUtils.diskPrecession(s, radius)[source]

Computes the precession of the disk due to the binary quadrupole moment. The precessions considered are kappa_r and kappa_z corresponding to the precession of the argument of periapsis and longitude of th ascending node, respectively.

Precssion frequency: Omega_p = Omega - Kappa Omega = sqrt((G*mu/r^3)*(1 + 3*alpha/r^2)) == orbital frequency

Parameters:

s: Tipsy snapshot

r: numpy array

array of radial bins centers [AU]

Returns:

——-

Kappa array: numpy array

2 x len(rBinEdges)-1 array containing precession at each radial point in 1/s
diskpy.ICgen.binaryUtils.findCBResonances(s, r, r_min, r_max, m_max=4, l_max=4, bins=50)[source]
Given Tipsy snapshot, computes the resonances of disk on binary as a function of orbital angular frequency omega. Disk radius, in au, is convered to angular frequency which will then be used to compute corotation and inner/outer Lindblad resonances.

Note: r given MUST correspond to r over which de/dt was calculated. Otherwise, scale gets all messed up

!!! NOTE: This function is awful and deprecated — do NOT use it. Instead, use calc_LB_resonance !!!

Parameters:

s: Tipsy-format snapshot

r: array

radius array over which de/dt was calculated

r_min,r_max: floats

min/maximum disk radius for calculations (au)

bins: int

number of radial bins to calculate over

m_max,l_max: ints

maximum orders of (m,l) LR
Returns:

Orbital frequency: numpy array

for corotation and inner/outer resonances and radii as float and numpy arrays
diskpy.ICgen.binaryUtils.find_crit_radius(r, array, toFind, num=1000)[source]

Given an array as a function of radius, the array value to search for, the radial range to search over and the search resolution, find the corresponding radius. Assumed that array is calculated as a function of r.

Parameters:

r: array

array of radial points

array: array

array of quantity of interest (could be surface density?) as a function of r

toFind: float

array value you’re looking for

num: int

resolution of search
Returns:

critial_radius: float

radius at which array(critical_radius) = toFind (approximately)
diskpy.ICgen.binaryUtils.forcedEccentricity(binary_sys, r)[source]

Given a binary class object and an array of radial points in the disk, compute the forced eccentricity defined by Moriwaki et al. 2004 eqn 9 to first order. Extra factor of 2 to give e_pumped instead of e_forced. Note: This only applies when e_binary != 0 and when m2/(m1 + m2) != 0.5 (i.e. only applies for eccentric, non-equal mass binary)

Parameters:

binary_sys : binary.Binary class object

r : array

array of radii in AU
Returns:

e_forced : array

array of len(r)
diskpy.ICgen.binaryUtils.linearMomentumEffects(x1, x2, v1, v2, m1, m2, accretion)[source]

Given initial binary system parameters and an array tracking the accretion events, calculate the effects of accretion on the semimajor axis and eccentricity of the binary system.

Inputs: Assume all input arrays are in simulation units Masses of primary and secondary (m1, m2 in Msol) Position arrays of primary and secondary x1, x2 (in AU) Velocity arrays of primary and secondary v1, v2 (in km/s) Numpy array of accretion events of the form [m vx vy vz …] for each accreted gas particle at time of accretion

Output: Semimajor axis, eccentricity of binary system after accretion events

diskpy.ICgen.binaryUtils.orbElemsVsRadius(s, rBinEdges, average=False)[source]

Computes the orbital elements for disk particles about a binary system in given radial bins. Assumes center of mass has v ~ 0

Parameters:

s: Tipsy snapshot

rBinEdges: numpy array

Radial bin edges [AU] preferably calculated using binaryUtils.calcDiskRadialBins

average: bool

True -> average over all particles in bin, false -> randomly select 1 particle in bin
Returns:

orbElems: numpy array

6 x len(rBinEdges) - 1 containing orbital elements at each radial bin as e, a, i, Omega, w, nu
diskpy.ICgen.binaryUtils.torqueVsRadius(s, rBinEdges)[source]

Takes in pynbody snapshot s for a binary system with a CB disk returns torque per unit mass vs radius and the approximate radius of the bin where that torque was calculated. Note, I only care about the z component of the torque since this is a circumbinary disk system Note: This function is best for returning proper disk radial bins.

Parameters:

s: pynbody snapshot of binary + CB disk

Bins: int

number of radial bins to make the calculation

r_in, r_out: floats

Inner and outer radii over which torque is calculated (au)
Returns:

tau: numpy array

Torque per unit mass as function of radius (cgs vs au)

diskpy.ICgen.calc_temp module

Calculates the temperature at a given r for a power law:
T = T0(r/r0)^Tpower

If r_in has no units, au are assumed

Created on Wed Jan 29 13:20:04 2014

@author: ibackus

class diskpy.ICgen.calc_temp.T(ICobj)[source]

Calculates temperature as a function of radius. Updating the settings contained in the ICobject will automatically cause T to be calculated accordingly. ICobject should contain the attribute ‘settings’. Settings for temperature are contained in settings.physical

USAGE:

T = calc_temp.T(ICobject) # Generate function to calculate T

T(r) # Calculate T at r

There are multiple kinds of available temperature profiles. They are set in ICobject.settings.physical.kind. They are:

‘powerlaw’ : (default)

Follows a power law, T = T0 * (r/r0)**Tpower Settings needed:

T0 r0 Tpower Tmin (optional) Tmax (optional)
‘MQWS’
Mimics Mayer, Quinn et. all 2004 ICs. Settings needed:
T0 r0 (same as r_in in the paper) Tmin Tmax (optional)
T_adiabatic(r)[source]

Estimates the adiabatic temperature profile as a function of r

setup_interior must be run first

T_nocut(r)[source]

Calculates the temperature as a function of r ignoring any interior cutoff. The interior cutoff would affect T for adiabatic disks

Returns:

T : SimArray

temperature
setup_interior()[source]

Sets up the interior temperature profile.

For an adiabatic disk, this is an adiabatic profile inside of the cutoff. Otherwise, nothing is done.

Requires that sigma already be calculated

diskpy.ICgen.calc_velocity module

Created on Wed Apr 9 15:39:28 2014

@author: ibackus

diskpy.ICgen.calc_velocity.v_xy(f, param, changbin=None, nr=50, min_per_bin=100, changa_preset=None, max_particles=None, est_eps=True, changa_args=”)[source]

Attempts to calculate the circular velocities for particles in a thin (not flat) keplerian disk. Also estimates gravitational softening (eps) for the gas particles and a reasonable time step (dDelta)

Requires ChaNGa

Note that this will change the velocities IN f

Parameters:

f : tipsy snapshot

For a gaseous disk

param : dict

a dictionary containing params for changa. (see configparser)

changbin : str (OPTIONAL)

If set, should be the full path to the ChaNGa executable. If None, an attempt to find ChaNGa is made

nr : int (optional)

number of radial bins to use when averaging over accelerations

min_per_bin : int (optional)

The minimum number of particles to be in each bin. If there are too few particles in a bin, it is merged with an adjacent bin. Thus, actual number of radial bins may be less than nr.

changa_preset : str

Which ChaNGa execution preset to use (ie ‘mpi’, ‘local’, …). See ICgen_utils.changa_command

max_particles : int or None

Specifies the maximum number of particles to use for calculating accelerations and velocities. Setting a smaller number can speed up computation and save on memory but can yield noisier results. If None, max is unlimited.

est_eps : bool

Estimate eps (gravitational softening length). Default is True. If False, it is assumed eps has already been estimated

changa_args : str

Additional command line arguments to pass to changa
Returns:

dDelta : float

A reasonable time step for the simulation (in code units). Velocties and eps are updated in the snapshot.

diskpy.ICgen.iterativesolver module

Contains an iterative solver class to iteratively solve

Created on Thu Apr 20 10:36:48 2017

@author: ibackus

class diskpy.ICgen.iterativesolver.IterativeSolver(IC, R, z=None)[source]

Calculates the vertical density profile for a thin, extended gaseous disk orbiting a central star by solving vertical hydrostatic equilibrium using an iterative method.

The generator equation used is:

\[\begin{split}\rho_{n+1}(z) &= \rho_{n+1}(0)\exp\left[ -a \left( \frac{1}{R} - \frac{1}{\sqrt{R^2 + z^2}} \right) - 2 b II_{\rho}(z)\right] \\ a &\equiv \frac{G M m}{k_B T} \\ b &\equiv \frac{2\pi G m}{k_B T} \\ II_{\rho}(z) &\equiv \int_0^z\int_0^{z'} \rho(z'') dz'' dz'\end{split}\]

The normalization \(\rho_{n+1}(0)\) is calculated by the constraint that the integral of the density (numerically estimated) is equal to the surface density.

This generator typically converges for reasonable disk/star parameters. I’ve tried other generators, and they all failed horribly.

The initial guess is for a thin disk ignoring self gravity:

\[\begin{split}\rho_0(z) &= \frac{\Sigma}{\sqrt{2 \pi} h} \exp(-z^2/2 h^2) \\ h &\equiv \sqrt{R^3/a} = c_s/\Omega\end{split}\]

The 2nd integral of this is known which provides a basis for the first order solution:

\[\begin{split}II_{\rho,0}(z) &= z I_{\rho}(z) + h^2 \rho{0}(z) - \frac{\Sigma h}{\sqrt{2 \pi}} \\ I_{\rho,0}(z) &= \frac{\Sigma}{2}\mathrm{erf}(z/\sqrt{2}h)\end{split}\]

Which is then normalized numerically

fit(maxiter=50, ftol=1e-08)[source]

A convenience utility which iterates over self.iterate() and saves results to a results dict

This will also delete any guess except the first one

iterate()[source]

Perform one iteration. Guesses are appended as SimArrays to the list self.rho

rho0(z, h=None, sigma=None)[source]

Initial guess (a gaussian). this models a thin disk with no self-gravity. Generally, this will overestimate the disk scale height.

diskpy.ICgen.iterativesolver.estHeight2(rho, z)[source]

A simple way to estimate disk height from a density profile

diskpy.ICgen.iterativesolver.heightProfile(IC)[source]

Generates a disk height profile from the vertical density profile

Returns:

h, r : SimArray

Disk height and radius
diskpy.ICgen.iterativesolver.similarity(p1, p2, forcePositive=False)[source]

Attempts to quantify the similarity between 2 PDFs using a normalized and shifted version of the bhattacharyya distance, defined as:

1 - D(p1, p2)/D(p1, p1)

Where D(x, y) is the Bhattacharyya distance between x and y.

If force positive is set, negative values are set to zero.

A value << 1 indicates similar pdfs

diskpy.ICgen.make_sigma module

Calculates cubic spline interpolations for sigma(r) and probability(r) probability = 2*pi*r*sigma

Created on Mon Jan 27 13:00:52 2014

@author: ibackus

class diskpy.ICgen.make_sigma.sigma_gen(r_bins, sigmaBinned, CDF=None)[source]

A class to generate the surface density (sigma), probability density (pdf) and inverse cumulative distribution function (cdf_inv) as a function of r

USAGE:

# Generate sigma (settings is generated by ICgen_settings.py, see below)

import make_sigma sigma = make_sigma.sigma_gen(r, sigma)

# Calculate at various r positions:

sigma(r) # returns sigma evaluated at r sigma.sigma(r) # returns sigma evaluated at r pdf = sigma.pdf(r) # returns pdf evaluated at r cdf_inv = sigma.cdf_inv(m) # returns cdv_inv at m for 0 < m < 1

# Generate sigma with a precalulated CDF (as done before)

sigma = make_sigma.sigma_gen(r, sigma, CDF)

copy()[source]

Returns a copy of the sigma object

diskpy.ICgen.make_snapshot module

Contains functions for genrerating a simulation snapshot from an IC object which already has particle positions generated. Snapshot are made with snapshot_gen(IC)

Created on Fri Mar 21 15:11:31 2014

@author: ibackus

diskpy.ICgen.make_snapshot.init_director(IC, param=None)[source]

Creates a director dict for a PPD.

Parameters:

IC : ICobj

param : .param dict

If not supplied, param = IC.param
diskpy.ICgen.make_snapshot.init_param(IC, snapshot=None)[source]

Initializes a ChaNGa param dict (see also diskpy.pychanga.make_param) for an IC object.

Parameters:

IC : ICobject

Initial conditions object containing required settings for generating the param dict

snapshot : SimSnap

Snapshot to create param for. If None, IC.snapshot is used
Returns:

param : dict

Dict containing the param. Can be saved with diskpy.utils.configsave
diskpy.ICgen.make_snapshot.init_snapshot(IC)[source]

Initialize a snapshot for the IC object. Requires that positions have been created. Also sets:

  • pos
  • metals
  • temp
  • mass
  • star eps
Parameters:IC : ICobj
Returns:snapshot : SimSnap
diskpy.ICgen.make_snapshot.make_binary(IC, snapshot)[source]

Turns a snapshot for a single star into a snapshot of a binary system

Parameters:

IC : IC object

snapshot : SimSnap

Single star system to turn into a binary
Returns:

snapshotBinary : SimSnap

A binary version of the simulation snapshot
diskpy.ICgen.make_snapshot.setup_sinks(IC, snapshot, param)[source]

Sets up snapshot and param for stars that are sinks.

Parameters:

IC : IC obj

snapshot : SimSnap

param : param dict
Returns:

None
diskpy.ICgen.make_snapshot.sink_radius(IC)[source]

Determine a reasonable sink radius for the star particles depending on the star system type (e.g., single star, binary, etc…)

Parameters:

IC : IC object
Returns:

r_sink : SimArray

Sink radius for star particles
diskpy.ICgen.make_snapshot.snapshot_gen(IC)[source]

Generates a tipsy snapshot from the initial conditions object IC. Includes a routine to calculate the velocity.

Parameters:

IC : ICobj
Returns:

snapshot : SimSnap

Simulation snapshot with the velocity calculated

param : dict

dictionary containing info for a .param file

diskpy.ICgen.make_snapshotBinary module

Created on Fri Mar 21 15:11:31 2014

@author: ibackus @editor: dfleming -Note: indentation is 4 spaces in this file, not a tab!

diskpy.ICgen.make_snapshotBinary.snapshot_gen(ICobj)[source]

Generates a tipsy snapshot from the initial conditions object ICobj.

Returns snapshot, param

snapshot: tipsy snapshot param: dictionary containing info for a .param file

Note: Code has been edited (dflemin3) such that now it returns a snapshot for a circumbinary disk where initial conditions generated assuming star at origin of mass M. After gas initialized, replaced star at origin with binary system who’s center of mass lies at the origin and who’s mass m1 +m2 = M

diskpy.ICgen.make_snapshotSType module

Created on Fri Mar 21 15:11:31 2014

@author: ibackus @editor: dflemin3 -Note: indentation is 4 spaces in this file, not a tab!

This module initializes an S-type binary system in which the gas disk is around the primary, not both stars! Assumes a_bin >> r_disk such that the disk’s velocity is dominated by the influence of the primary.

diskpy.ICgen.make_snapshotSType.snapshot_gen(ICobj)[source]

Generates a tipsy snapshot from the initial conditions object ICobj.

Returns snapshot, param

snapshot: tipsy snapshot param: dictionary containing info for a .param file

Note: Code has been edited (dflemin3) such that now it returns a snapshot for a circumbinary disk where initial conditions generated assuming star at origin of mass M. After gas initialized, replaced star at origin with binary system who’s center of mass lies at the origin and who’s mass m1 +m2 = M

diskpy.ICgen.make_snapshotSTypeSupplied module

Created on Fri Mar 21 15:11:31 2014

@author: ibackus @editor: dflemin3 -Note: indentation is 4 spaces in this file, not a tab!

This module initializes an S-type binary system in which the gas disk is around the primary, not both stars! Assumes a_bin >> r_disk such that the disk’s velocity is dominated by the influence of the primary.

diskpy.ICgen.make_snapshotSTypeSupplied.snapshot_gen(ICobj)[source]

Generates a tipsy snapshot from the initial conditions object ICobj.

Returns snapshot, param

snapshot: tipsy snapshot param: dictionary containing info for a .param file

Note: Code has been edited (dflemin3) such that now it returns a snapshot for a circumbinary disk where initial conditions generated assuming star at origin of mass M. After gas initialized, replaced star at origin with binary system who’s center of mass lies at the origin and who’s mass m1 +m2 = M

diskpy.ICgen.pos_class module

Defines a function to randomly generate particle positions according to the desired surface density profile (sigma vs r) and the vertical profile (rho vs r,z).

Created on Mon Jan 27 18:48:04 2014

@author: ibackus

class diskpy.ICgen.pos_class.pos(ICobj, method=None, generate=True, seed=None)[source]

position class. Generates particle positions from rho and sigma

USAGE: # method = ‘grid’ or ‘random’ pos = pos_class.pos(ICobj, method)

ICobj should be an initial conditions object (ICgen.IC) with rho already calculated.

diskpy.ICgen.rhosolver module

Created on Mon Aug 3 16:28:37 2015

@author: ibackus

diskpy.ICgen.rhosolver.calc_rho(IC, r=None, **kwargs)[source]

Calculates the density for the initial conditions object IC by numerically solving hydrostatic equilibrium (see vertical_solver)

Parameters:

IC : IC object

Initial conditions object

r : SimArray

(optional) intial bin radii: not all will be used

**kwargs : keyword arguments to pass to the root finder used

(scipy.optimize.newton_krylov)
Returns:

R : 1D SimArray

Radial bins the density is calculated at

z : 2D SimArray

z points the density is calculated at. z[i,j] is the ith z position at R[j]

rho : 2D SimArray

Density as a function of (z,R). rho[i,j] is calculated at (z[i,j], R[j])
diskpy.ICgen.rhosolver.cdf_inv(z, r, rho)[source]

Calculates the inverse CDF as a function of r over the whole disk

Parameters:

z : SimArray or array

2D z-positions that rho is calculated at. z[i,j] is the ith z bin at r[j]

r : SimArray or array

Radial bins (1D array) the z and rho are calculated at

rho : SimArray or array

2D array of density values. rho[i,j] is rho at z[i,j], r[j]
Returns:

f : function

Inverse CDF. f(m, r) = z returns the z value for radius and 0<m<1. r and m are 1-D arrays of the same length, or numbers.
diskpy.ICgen.rhosolver.cdf_inv_z(z, rho)[source]

Calculates the inverse of the cumulative distribution function for probability as a function of z for a given r (ie gives you z as a function of the CDF)

Parameters:

z : array or SimArray

Positions to calculate over. 1D array

rho: array or SimArray

Density as a function of z. Treated as an un-normalized PDF. 1D array
Returns:

f : array

Normalized CDF

z : array or SimArray

z as a functin of the normalized CDF. Monotonically increasing

Notes

To ensure z, f are montonically increasing, some values are dropped. The remaining values of f are padded with 2, and the values of z are padded with z.max()

diskpy.ICgen.rhosolver.loadrho(IC, f)[source]

Loads a rhosolver object from f (a file or dictionary) and saves it to IC.rho

Parameters:

f : dict or str

Either a dictionary containing everything required to load rho or a filename pointing to a pickled dictionary
Returns:

rho : rhosolver

An intialied rho solver object
diskpy.ICgen.rhosolver.rhointerp(z, r, rho)[source]

Generates a callable interpolation of rho on the z, r points

Parameters:

z : 2D SimArray or array

z[i,j] is the ith z value at r[j]

r : 1D SimArray or array

Radial positions

rho : SimArray

density at points z[i,j], r[j]
Returns:

rhospline : function

An interpolation function for estimating rho as a function of z, r
class diskpy.ICgen.rhosolver.rhosolver(IC)[source]

Defines the rho class that allows the solving of vertical hydrostatic equilibrium over the disk, and generates callable methods for estimating density and the normalized inverse CDF over the disk.

Examples

Initialize rho, solve vertical equilibrium

>>> IC.maker.sigma_gen()
>>> rho = rhosolver(IC)
>>> rho.solve(maxiter=100)

Rho is now callable:

>>> rho(z, r)
>>> rho.cdf_inv(r, m)

Save rho

>>> rho_dict = {'z': rho.z_bins, 'r': rho.r_bins, 'rho': rho.rho_binned}
>>> pickle.dump(rho_dict, open('rhofile.p','w'))

Load rho

>>> rho.load('rhofile.p')
load(f)[source]

Load/initialize from a dictionary or a file

Parameters:

f : dict or str

Either a dictionary containing everything required to load rho or a filename pointing to a pickled dictionary
solve(**kwargs)[source]

Solves the hydrostatic vertical equilibrium for the disk to find the density. Also calculates the normalized inverse CDF

Parameters:

kwargs :

(optional) key word arguments to pass to the root finder [scipy.optimize.newton_krylov]
diskpy.ICgen.rhosolver.setup_r_bins(IC, r=None)[source]

diskpy.ICgen.sigma_profile module

Created on Mon Jun 23 10:17:53 2014

@author: ibackus

diskpy.ICgen.sigma_profile.MQWS(settings, T)[source]

Generates a surface density profile as the per method used in Mayer, Quinn, Wadsley, and Stadel 2004

** ARGUMENTS ** NOTE: if units are not supplied, assumed units are AU, Msol

settings : IC settings
settings like those contained in an IC object (see ICgen_settings.py)
T : callable
A function to calculate temperature as a function of radius

** RETURNS **

r : SimArray
Radii at which sigma is calculated
sigma : SimArray
Surface density profile as a function of R
diskpy.ICgen.sigma_profile.gaussian_ring(settings)[source]

Generates a gaussian ring surface density profile according to:

\[\Sigma = \Sigma_0 exp(-(R-R_d)^2/2a^2)\]
\[\Sigma_0 = M_d/(2\pi)^{3/2} a R_d\]

Here we call a the ringwidth.

The max radius is determined automatically

Parameters:

settings : IC settings

settings like those contained in an IC object (see ICgen_settings.py)
Returns:

R : SimArray

Radii at which sigma is calculated

sigma : SimArray

Surface density profile as a function of R
diskpy.ICgen.sigma_profile.make_profile(ICobj)[source]

A wrapper for generating surface density profiles according to the IC object.

Settings for the profile are defined in ICobj.settings. Which profile gets used is defined by ICobj.settings.sigma.kind

Currently available kinds are:

viscous powerlaw MQWS

RETURNS

r : SimArray
Radii at which sigma is calculated
sigma : SimArray
Surface density profile as a function of R
diskpy.ICgen.sigma_profile.powerlaw(settings, T=None)[source]

Generates a surface density profile according to a powerlaw sigma ~ r^p with a smooth interior cutoff and smooth exterior exponential cutoff.

ARGUMENTS

settings : IC settings
settings like those contained in an IC object (see ICgen_settings.py)
T : callable function
Function that returns temperature of the disk as a function of radius IF none, a powerlaw temperature is assumed

RETURNS

R : SimArray
Radii at which sigma is calculated
sigma : SimArray
Surface density profile as a function of R
diskpy.ICgen.sigma_profile.viscous(settings)[source]

Generates a surface density profile derived from a self-similarity solution for a viscous disk, according to:

sigma ~ r^-gamma exp(-r^(2-gamma))

Where r is a dimensionless radius and gamma is a constant less than 2. Rd (disk radius) is defined as the radius containing 95% of the disk mass

ARGUMENTS

settings : IC settings
settings like those contained in an IC object (see ICgen_settings.py)

RETURNS

R : SimArray
Radii at which sigma is calculated
sigma : SimArray
Surface density profile as a function of R

diskpy.ICgen.vertical_solver module

Defines the vertical_solver class for solving vertical hydrostatic equilibrium at a fixed radius.

Created on Mon Aug 3 16:01:44 2015

@author: ibackus

diskpy.ICgen.vertical_solver.dI(I, dz=1)[source]

Derivative of I with respect to z

diskpy.ICgen.vertical_solver.estHeight(sigma, T, M, m, R)[source]

Roughly estimates the height of the disk, in both the regime of disk self gravity dominating and star gravity dominating

diskpy.ICgen.vertical_solver.rho0(z, h=1)[source]

Approximate functional form for rho. All quantities are dimensionless

To first order, h should equal 1.

Parameters:

z : Array

z positions to estimate rho at

h : float

dimensionless scale height. Default=1
Returns:

rho : array

Dimensionless density
diskpy.ICgen.vertical_solver.setup(IC, R)[source]

Initialize values for solving vertical hydrostatic equlibrium.

Parameters:

IC : IC object

Initial conditions object

R : SimArray

Radius at which to solve hydrostatic equilibrium
Returns:

z : array

Dimensionless z bins, in units of h

r : float

Radius in units of h

c : float

Constant in residual equation

zscale : SimArray

Amount to scale z (all lengths) by to get dimensions back. ie, Z(actual) = z*zscale

rhoscale : SimArray

Scale for dimensionless rho
class diskpy.ICgen.vertical_solver.vertical_solver(IC, R, rescale_rho=True, rhoguess=None)[source]

Defines the class which solves vertical hydrostatic equilibrium for a thin Keplerian disk orbiting a central star

Examples

>>> IC = ICgen.load('IC.p')
>>> R = SimArray(0.5, 'au')
>>> solver = vertical_solver(IC, R)
>>> solver.fit(maxiter=100)
>>> z = solver.results['z']
>>> rho = solver.results['rho']

This example loads up an already created initial conditions object then solves for hydrostatic equilibrium at R=0.5 AU

A better solver is now implemented. see iterativesolver.IterativeSolver

fit(**kwargs)[source]

Attempts to find the root of the residual of the differential equation for rho at all z, given the initial guess self.rho0

Uses a Newton’s method with the Krylov approximation for the Jacobian (see scipy.optimize.newton_krylov)

Parameters:

**kwargs : keyword arguments

Additional arguments passed to the root finder [see scipy.optimize.newton_krylov]
Returns:

rho : array

Dimensionless density, saves it to self.rho
fitrobust(nh=10, scansize=5, **kwargs)[source]

Attempts to find the root of the residual of the differential equation for rho at all z, given the initial guess generated by different values of h.

This routine is slower than fit but can avoid local minima. If fit doesn’t seem to be converging to the right values, consider using this.

h is a dimensionless disk height, which is near 1. This algorithm scans over values of h to create initial guesses rho0. fit() is then called, and the results are stored. The fit value of rho which is globally closest to 0 is then used.

Uses a Newton’s method with the Krylov approximation for the Jacobian (see scipy.optimize.newton_krylov)

Parameters:

nh : int

(optional) Number of h values to scan over

scansize : float

(optional) Sets the range of h values to scan over. Scans from h0/scansize to h0*scansize

**kwargs : keyword arguments

Additional arguments passed to the root finder [see scipy.optimize.newton_krylov]
Returns:

rho : array

Dimensionless density, saves it to self.rho
residual(rho)[source]

Estimates the residual of the differential equation for rho for all z

Module contents

Created on Fri Aug 7 10:11:48 2015

@author: ibackus

class diskpy.ICgen.IC(r=None, sigma=None, CDF=None, profile_kind=None, settings=None)[source]

Defines the IC class.

INITIALIZING NEW INITIAL CONDITIONS

# Initialize a blank IC object (with no surface density profile yet made)

IC = ICgen.IC()

# Initialize an IC object and surface density profile using default settings

IC = ICgen.IC(profile_kind=’powerlaw’) IC = ICgen.IC(profile_kind=’MQWS’)

# Initialize IC object from 1-D SimArrays r, sigma (surface density)

IC = ICgen.IC(r, sigma)

Optionally, the CDF for the surface density profile can be supplied to speed up generation of sigma. To do that:

IC = ICgen.IC(r, sigma, CDF)

Or, the input can be a the filename of a pickled dictionary containing ‘r’, ‘sigma’, and optionally ‘CDF’

Settings can also be entered manually if needed

settings = pickle.load(open(‘settings.p’,’r’)) IC = ICgen.IC(settings = settings)

GENERATING INITIAL CONDITIONS

Qest(r=None)[source]

Estimate Toomre Q at r (optional) for ICs, assuming omega=epicyclic frequency. Ignores disk self-gravity

generate(restart=False)[source]

Runs through all the steps to generate a set of initial conditions

IF restart=True, it picks up at the last completed step

diskpy.ICgen.load(filename)[source]
class diskpy.ICgen.rhosolver(IC)[source]

Defines the rho class that allows the solving of vertical hydrostatic equilibrium over the disk, and generates callable methods for estimating density and the normalized inverse CDF over the disk.

Examples

Initialize rho, solve vertical equilibrium

>>> IC.maker.sigma_gen()
>>> rho = rhosolver(IC)
>>> rho.solve(maxiter=100)

Rho is now callable:

>>> rho(z, r)
>>> rho.cdf_inv(r, m)

Save rho

>>> rho_dict = {'z': rho.z_bins, 'r': rho.r_bins, 'rho': rho.rho_binned}
>>> pickle.dump(rho_dict, open('rhofile.p','w'))

Load rho

>>> rho.load('rhofile.p')
load(f)[source]

Load/initialize from a dictionary or a file

Parameters:

f : dict or str

Either a dictionary containing everything required to load rho or a filename pointing to a pickled dictionary
solve(**kwargs)[source]

Solves the hydrostatic vertical equilibrium for the disk to find the density. Also calculates the normalized inverse CDF

Parameters:

kwargs :

(optional) key word arguments to pass to the root finder [scipy.optimize.newton_krylov]
diskpy.ICgen.loadrho(IC, f)[source]

Loads a rhosolver object from f (a file or dictionary) and saves it to IC.rho

Parameters:

f : dict or str

Either a dictionary containing everything required to load rho or a filename pointing to a pickled dictionary
Returns:

rho : rhosolver

An intialied rho solver object