Reference¶

GravHopper consists of the following submodules:

Simulation
IC
jbgrav
unitconverter
Exceptions

Simulation¶

General usage:

from gravhopper import Simulation
sim = Simulation(dt=<timestep>, eps=<softening>)
sim.add_IC(<initial conditions>)
sim.run(<N>)
sim.plot_particles()
# Other analysis using sim.positions, sim.velocities, sim.times

class gravhopper.Simulation(dt=<Quantity 1. Myr>, eps=<Quantity 100. pc>, algorithm='tree')[source]¶

Main class for N-body simulation.

Np¶

Number of particles

Type:: int

Nsnap¶

Number of snapshots

Type:: int

positions¶

numpy array of the positions of each particle at each snapshot

Type:: Quantity array of dimension (Nsnap,Np,3)

velocities¶

numpy array of the velocities of each particle at each snapshot

Type:: Quantity array of dimension (Nsnap,Np,3)

masses¶

numpy array of the masses of each particle

Type:: Quantity array of length Np

times¶

numpy array of the time of each snapshot

Type:: Quantity array of length Nsnap

lenunit¶

Default internal length unit

Type:: Unit

velunit¶

Default internal velocity unit

Type:: Unit

accelunit¶

Default internal acceleration unit

Type:: Unit

__init__(dt=<Quantity 1. Myr>, eps=<Quantity 100. pc>, algorithm='tree')[source]¶

Initializes Simulation object.

Parameters:

dt (Quantity) – Time step length. Default: 1 Myr
eps (Quantity) – Gravitational softening length. Default: 100 pc
algorithm (str ('tree' or 'direct')) – Gravitational force calculation algorithm. Default: ‘tree’

Example

To create a new simulation with a timestep of 1 year:

sim = Simulation(dt=1.0*u.yr)

add_IC(newIC)[source]¶

Adds particles to the initial conditions.

Parameters:

newIC (dict) –

Properties of new particles to add. Must have the following key-value pairs:

pos: an array of positions
vel: an array of velocities
mass: an array of masses

All should be astropy Quantities, with shape (Np,3) or just (3) if a single particle.

Example

Create a simulation object whose initial conditions consist of a particle of mass 10:sup:8 M:sub:sun at a position 10 kpc from the origin on the x-axis, and a velocity of 200 km/s in the positive y-direction:

sim = Simulation()
sim.add_IC( {'pos':np.array([10,0,0])*u.kpc, 'vel':np.array([0,200,0])*u.km/u.s,
    'mass':np.array([1e8])*u.Msun} )

See also

IC: Class containing static functions that generate initial conditions that can be added using add_IC().

add_external_force(fn, args=None, agama_units=None)[source]¶

Add an external position-dependent force to the simulation.

Forces can be in the form of:

A function that takes two arguments: an (Np,3) array of positions (must be Quantities) that contains the positions where the accelerations are to be calculated, and an additional argument that is a dictionary containing any extra parameters you need.
A galpy potential object.
A gala Potential object.
An Agama Potential object.

You may add as many external forces as you want - they will be summed together along with the N-body force.

Parameters:

fn (function or galpy potential.Potential object or gala potential.PotentialBase object) – External force function to add
args (dict) – Any extra parameters that should be passed to the function when it is called
agama_units (dict) – Agama potentials take unitless values, so assume that the Agama potential has been defined with this set of units. The dict should contain the keys ‘lenunit’ and ‘accelunit’, which should each by an Astropy Unit or Quantity. Required for Agama potentials, ignored otherwise.

Example

This function calculates an external force from a single point source:

def my_point_source_force(pos, args):
     # acceleration is G M / r^2
     GM = const.G * args['mass']
     d_pos = pos - args['pos']
     r2 = (d_pos**2).sum(axis=1)
     rhat = d_pos / np.sqrt(r2)
     return rhat * GM / r2

To add the force from a particle at 10kpc on the x-axis with a mass of 1e8 Msun:

mysimulation = Simulation()
mysimulation.add_external_force(my_point_source_force, {'mass':1e8*u.Msun,
  'pos':np.array([10,0,0])*u.kpc})

add_external_timedependent_force(fn, agama_units=None, args=None)[source]¶

Add an external time-dependent force to the simulation.

Forces can be in the form of:

A function that takes three arguments: an (Np,3) array of positions (must be Quantities) that contains the positions where the accelerations are to be calculated, a scalar time, and an additional argument that is a dictionary containing any extra parameters you need.
A galpy potential object.
A gala Potential object.
An Agama Potential object.

You may add as many external forces as you want - they will be summed together along with the N-body force.

Parameters:

fn (function or galpy potential.Potential object or gala potential.PotentialBase object) – External force function to add
args (dict) – Any extra parameters that should be passed to the function when it is called
agama_units (dict) – Agama potentials take unitless values, so assume that the Agama potential has been defined with this set of units. The dict should contain the keys ‘lenunit’ and ‘accelunit’, which should each by an Astropy Unit or Quantity. Required for Agama potentials, ignored otherwise.

Example

This function calculates an external force from a single point source whose mass oscilllates in time:

def my_oscillating_point_source_force(pos, time, args):
     # args has 3 parameters:
     #    args['pos']: position of mass
     #    args['massamplitude']: amplitude of oscillating mass (should have mass units)
     #    args['period']: period of oscillation of mass (should have time units)
     # acceleration is G M / r^2
     GM = const.G * args['massamplitude'] * (np.cos(2. * np.pi * time / args['period']) + 1.)
     d_pos = pos - args['pos']
     r2 = (d_pos**2).sum(axis=1)
     rhat = d_pos / np.sqrt(r2)
     return rhat * GM / r2

To add the force from a particle at 10kpc on the x-axis with a mass of 10:sup:8 M:sub:sun that oscillates every 100 Myr:

mysimulation = Simulation()
mysimulation.add_external_timedependent_force(my_point_source_force, {'massamplitude':1e8*u.Msun,
  'period':100*u.Myr, 'pos':np.array([10,0,0])*u.kpc})

add_external_velocitydependent_force(fn, args=None)[source]¶

Add an external velocity-dependent force to the simulation.

Forces can be in the form of:

A function that takes three arguments: an (Np,3) array of positions (must be Quantities) that contains the positions where the accelerations are to be calculated, an (Np,3) array of velocities (must be Quantities), and an additional argument that is a dictionary containing any extra parameters you need.
A galpy potential object (must derive from galpy.potential.DissipativeForce).

You may add as many external forces as you want - they will be summed together along with the N-body force.

Parameters:

fn (function or galpy.potential.DissipativeForcePotential object) – External force function to add
args (dict) – Any extra parameters that should be passed to the function when it is called

Example

This function adds on an external force that goes in the opposite directly of the current velocity of every particle with magnitude abs(velocity) / timescale given in args:

def my_friction_force(pos, vel, args):
   # args has 1 parameter:
   #    args['t0']:  timescale (should have time units)
   velmag = np.sqrt(np.sum(vel**2, axis=1))
   forcemag = velmag / args['t0']
   forcearray = -vel/velmag[:,np.newaxis] * forcemag[:,np.newaxis]
   return forcearray

Then to add a force that slows all particles down on a 100 Myr timescale:

mysimulation = Simulation()
mysimulation.add_external_velocitydependent_force(my_friction_force, {'t0':100*u.Myr})

calculate_acceleration(time=None)[source]¶

Calculate acceleration for particle positions at current_snap() due to gravitational N-body force, plus any external forces that have been added.

Parameters:: time (Quantity) – Time at which to calculate any time-dependent external forces. Default: None
Returns:: acceleration – An (Np,3) numpy array of the acceleration vector calculated for each particle
Return type:: array
Raises:: UnknownAlgorithmException – If the algorithm parameter is not ‘tree’ or ‘direct’

current_snap()[source]¶

Return the current snapshot.

Returns:

snap –

snap[‘pos’] is an (Np,3) array of positions
snap[‘vel’] is an (Np,3) array of velocities
snap[‘mass’] is a length-Np array of masses

Return type:

dict

get_algorithm()[source]¶

Returns current simulation gravitational algorithm.

Returns:: algorithm – ‘direct’ for direct summation, ‘tree’ for Barnes-Hut tree
Return type:: str

See also

set_algorithm(): Sets the gravitational force algorithm.

get_dt()[source]¶

Returns simulation time step.

Returns:: dt – Time step length
Return type:: Quantity

See also

set_dt(): Sets simulation time step.

get_eps()[source]¶

Returns simulation gravitational softening length.

Returns:: eps – Gravitational softening length
Return type:: Quantity

See also

set_eps(): Sets gravitational softening length.

movie_particles(fname, fps=25, ax=None, skip=None, timeformat='{0:.1f}', *args, **kwargs)[source]¶

Create a movie of the particles using the plot_particles() function.

Parameters:

fname (str) – Movie output file name
fps (int) – Frames per second (default: 25)
ax (matplotlib Axis, optional) – Axis on which to plot. If missing or None, uses plt.subplot(111, aspect=1.0) to create a new Axis object and closes the figure after the movie is made.
skip (int, optional) – Skip every N frames (e.g. skip=5 only has 1/5th of the full number of frames).
xlim (array-like, optional) – x limits of the movie (default: encompasses the full extent any particle reaches)
ylim (array-like, optional) – y limits of the movie (default: encompasses the full extent any particle reaches)
*args (object) – Passed through to plot_particles()
**kwargs (dict) – Passed through to plot_particles()

plot_particles(parm='pos', coords='xy', snap='final', xlim=None, ylim=None, s=0.2, unit=None, ax=None, timeformat='{0:.1f}', nolabels=False, particle_range=None, **kwargs)[source]¶

Scatter plot of particle positions or velocities for a snapshot.

Parameters:

parms (str) – ‘pos’ or ‘vel’ to plot particle positions or velocities (default: ‘pos’)
coords (str) – ‘xy’, ‘xz’, or ‘yz’ to plot the different Cartesian projections (default ‘xy’)
snap (int or str) – Which snapshot to plot. Either a snapshot number, ‘final’ for the last snapshot, or ‘IC’ for the initial conditions (default: ‘final’)
xlim (array-like, optional) – x limits of plot
ylim (array-like, optional) – y limits of plot
s (float) – Scatter plot point size (default: 0.2)
unit (astropy Unit, optional) – Plot the quantities in the given unit system. The default is whatever units the positions or velocities array is in, which is kpc or km/s by default but can be changed.
ax (matplotlib Axis, optional) – Axis on which to plot. If missing or None, uses plt.subplot(111, aspect=1.0) to create a new Axis object.
timeformat (str or False) – Format string for snapshot time in title, or False for no title (default: ‘{0:.1f}’)
nolabels (bool, optional) – If True, do not label x and y axes (default: False)
particle_range (array-like, optional) – Only plot particles with indices in the slice particle_range[0]:particle_range[1]
**kwargs (dict) – All additional keyword arguments are passed onto plt.scatter()

Returns:

scatter – The scatter plot

Return type:

matplotlib PathCollection

prev_snap()[source]¶

Return the snapshot before the current snapshot.

Returns:

snap –

snap[‘pos’] is an (Np,3) array of positions
snap[‘vel’] is an (Np,3) array of velocities
snap[‘mass’] is a length-Np array of masses

Return type:

dict

pyn_snap(timestep=None)[source]¶

Return snapshot given by the timestep as a pynbody SimSnap.

Parameters:: timestep (int) – Snapshot number to return. Returns final snapshot if timestep is None (default: None)
Returns:: sim – Snapshot
Return type:: pynbody SimSnap

Notes

Returned snapshot will attempt to put the SimSnap in equivalent units as what the Simulation object is using internally.

Raises:: ExternalPackageException – If pynbody could not be imported.

reset()[source]¶: Reset a simulation to the state before init_run() was called. Initial conditions and external forces are preserved.

run(N=1)[source]¶

Run N timesteps of the simulation. Will either initialize a simulation that has not yet been run, or continue from the last snapshot if it has.

Parameters:: N (int) – Number of time steps to perform

set_algorithm(algorithm)[source]¶

Sets the simulation gravitational force algorithm.

Parameters:: algorithm (str ('tree' or 'direct')) – Direct summation or Barnes-Hut tree algorithm
Raises:: ValueError – If algorithm is not ‘tree’ or ‘direct’.

See also

get_algorithm(): Returns current gravitational force algorithm.

set_dt(dt)[source]¶

Sets the simulation time step.

Parameters:: dt (Quantity with dimensions of time) – Time step length
Raises:: ValueError – If dt does not have dimensions of time.

See also

get_dt(): Returns simulation time step.

set_eps(eps)[source]¶

Sets the simulation gravitational softening length.

Parameters:: eps (Quantity with dimensions of length) – Gravitational softening length
Raises:: ValueError – If eps does not have dimensions of length.

See also

get_eps(): Returns gravitational softening length.

snap(step)[source]¶

Return the given snapshot.

Parameters:: step (int) – Time step to return
Returns:: snap – snap[‘pos’] is an (Np,3) array of positions snap[‘vel’] is an (Np,3) array of velocities snap[‘mass’] is a length-Np array of masses
Return type:: dict

IC¶

General usage:

from gravhopper import IC
IC_array = IC.Plummer(totmass=<mass>, a=<length>, N=<N>)
sim.add_IC(IC_array)

Because random sampling is not guaranteed to have the center of mass or mean velocity of exactly zero, all functions that involve random sampling take an argument force_origin (default: True) that shifts the final particle positions and velocities to be at the origin, or to a pre-set other location or velocity using the center_pos and center_vel arguments.

All functions where GravHopper does the random sampling allow a random seed to be set using the seed argument to facilitate reproducibility. This seed is fed directly into numpy.random.default_rng().

class gravhopper.IC[source]¶

Namespace for holding static methods that create a variety of initial conditions.

static Hernquist(N=None, a=None, totmass=None, cutoff=10.0, center_pos=None, center_vel=None, force_origin=True, seed=None)[source]¶

Generate the initial conditions for an isotropic Hernquist model (Hernquist 1990).

Parameters:

N (int) – Number of particles
a (astropy Quantity of dimension length) – Scale radius
totmass (astropy Quantity of dimension mass) – Total mass
cutoff (float) – Cut off the distribution at cutoff times the scale radius
center_pos (3 element array-like Quantity, optional) – Force the center of mass of the IC to be at this position
center_vel (3 element array-like Quantity, optional) – Force the mean velocity of the IC to have this velocity
force_origin (bool) – Force the center of mass to be at the origin and the mean velocity to be zero; equivalent to setting center_pos=[0,0,0]*u.kpc and center_vel=[0,0,0]*u.km/u.s. Default is True unless center_pos and center_vel is set. If force_origin is True and only one of center_pos or center_vel is set, the other is set to zero.
seed ({None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional) – Seed to initialize random number generator to enable repeatable ICs.

Returns:

IC – Properties of new particles to add, which sample the given distribution function. Contains the following key/value pairs:

pos: an array of positions
vel: an array of velocities
mass: an array of masses

Each are astropy Quantities, with shape (Np,3).

Return type:

dict

Example

To create a Hernquist sphere with a total mass of 10:sup:10 solar masses, a scale radius of 1 kpc, sampled with 10,000 particles:

particles = IC.Hernquist(N=10000, a=1*u.kpc, totmass=1e10*u.Msun)

static Plummer(N=None, b=None, totmass=None, center_pos=None, center_vel=None, force_origin=True, seed=None)[source]¶

Generate the initial conditions for an isotropic Plummer model (BT eqs. 4.83 with n=5, 4.92, 2.44b).

Parameters:

N (int) – Number of particles
b (astropy Quantity of dimension length) – Scale radius
totmass (astropy Quantity of dimensions mass) – Total mass
center_pos (3 element array-like Quantity, optional) – Force the center of mass of the IC to be at this position
center_vel (3 element array-like Quantity, optional) – Force the mean velocity of the IC to have this velocity
force_origin (bool) – Force the center of mass to be at the origin and the mean velocity to be zero; equivalent to setting center_pos=[0,0,0]*u.kpc and center_vel=[0,0,0]*u.km/u.s. Default is True unless center_pos and center_vel is set. If force_origin is True and only one of center_pos or center_vel is set, the other is set to zero.
seed ({None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional) – Seed to initialize random number generator to enable repeatable ICs.

Returns:

IC – Properties of new particles to add, which sample the given distribution function. Contains the following key/value pairs:

pos: an array of positions
vel: an array of velocities
mass: an array of masses

Each are astropy Quantities, with shape (Np,3).

Return type:

dict

Example

To create a Plummer sphere with scale radius 1 pc and a total mass of 10:sup:6 M:sub:sun sampled with 10,000 particles:

particles = IC.Plummer(N=10000, b=1*u.pc, totmass=1e6*u.Msun)

static TSIS(N=None, maxrad=None, totmass=None, center_pos=None, center_vel=None, force_origin=True, seed=None)[source]¶

Generate the initial conditions for a Truncated Singular Isothermal Sphere (e.g. BT eqs 4.103 and 4.104 with a maximum radius imposed). Note that this is not a true equilibrium because of the truncation – it will be in apporoximate equilibrium in the inner regions, at least at first, but not in the outer regions.

Parameters:

N (int) – Number of particles
maxrad (astropy Quantity with length dimensions) – Truncation radius
totmass (astropy Quantity with mass dimensions) – Total mass
center_pos (3 element array-like Quantity, optional) – Force the center of mass of the IC to be at this position
center_vel (3 element array-like Quantity, optional) – Force the mean velocity of the IC to have this velocity
force_origin (bool) – Force the center of mass to be at the origin and the mean velocity to be zero; equivalent to setting center_pos=[0,0,0]*u.kpc and center_vel=[0,0,0]*u.km/u.s. Default is True unless center_pos and center_vel is set. If force_origin is True and only one of center_pos or center_vel is set, the other is set to zero.
seed ({None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional) – Seed to initialize random number generator to enable repeatable ICs.

Returns:

IC – Properties of new particles to add, which sample the given distribution function. Contains the following key/value pairs:

pos: an array of positions
vel: an array of velocities
mass: an array of masses

Each are astropy Quantities, with shape (Np,3).

Return type:

dict

Example

To create a truncated singular isothermal sphere with a total mass of 10:sup:11 solar masses, a truncation radius of 100 kpc, sampled with 10,000 particles:

particles = IC.TSiS(N=10000, maxrad=100*u.kpc, totmass=1e11*u.Msun)

static expdisk(sigma0=None, Rd=None, z0=None, sigmaR_Rd=None, external_rotcurve=None, N=None, center_pos=None, center_vel=None, force_origin=True, seed=None)[source]¶

Generates initial conditions of an exponential disk with a sech^2 vertical distribution that is in (very) approximate equilibrium: rho(R,z) = (sigma0 / 2 z0) exp(-R/Rd) sech^2(z/z0)

Parameters:

sigma0 (astropy Quantity with dimensions of surface density) – Central surface density
Rd (astropy Quantity with dimensions of length) – Radial exponential scale length
z0 (astropy Quantity with dimensions of length) – Vertical scale height
sigmaR_Rd (astropy Quantity with dimensions of velocity) – Radial velocity dispersion at R=Rd
external_rotcurve (function or None) – Function that returns the circular velocity of any external potential that contributes to the rotation curve aside from the disk itself. The function should accept input as an astropy Quantity of dimension length, and should return an astropy Quantity of dimension velocity.
N (int) – Number of particles
center_pos (3 element array-like Quantity, optional) – Force the center of mass of the IC to be at this position
center_vel (3 element array-like Quantity, optional) – Force the mean velocity of the IC to have this velocity
force_origin (bool) – Force the center of mass to be at the origin and the mean velocity to be zero; equivalent to setting center_pos=[0,0,0]*u.kpc and center_vel=[0,0,0]*u.km/u.s. Default is True unless center_pos and center_vel is set. If force_origin is True and only one of center_pos or center_vel is set, the other is set to zero.
seed ({None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional) – Seed to initialize random number generator to enable repeatable ICs.

Returns:

IC – Properties of new particles to add, which sample the given distribution function. Contains the following key/value pairs:

pos: an array of positions
vel: an array of velocities
mass: an array of masses

Each are astropy Quantities, with shape (Np,3).

Return type:

dict

Example

To create an exponential disk that is in a background logarithmic halo potential that generates a flat rotation curve of 200 km/s:

particles = IC.expdisk(N=10000, sigma0=200*u.Msun/u.pc**2, Rd=2*u.kpc,
    z0=0.5*u.kpc, sigmaR_Rd=10*u.km/u.s,
    external_rotcurve=lambda x: 200*u.km/u.s)

static from_galpy_df(df, N=None, totmass=None, center_pos=None, center_vel=None, force_origin=True)[source]¶

Sample a galpy sphericaldf distribution function object and return as an IC.

Parameters:

df (galpy df.sphericaldf object) – Distribution function to sample. Assumed to be in the ro=8, vo=220 unit system.
N (int) – Number of particles
totmass (astropy Quantity) – Total mass
center_pos (3 element array-like Quantity, optional) – Force the center of mass of the IC to be at this position
center_vel (3 element array-like Quantity, optional) – Force the mean velocity of the IC to have this velocity
force_origin (bool) – Force the center of mass to be at the origin and the mean velocity to be zero; equivalent to setting center_pos=[0,0,0]*u.kpc and center_vel=[0,0,0]*u.km/u.s. Default is True unless center_pos and center_vel is set. If force_origin is True and only one of center_pos or center_vel is set, the other is set to zero.

Returns:

IC – Properties of new particles to add, which sample the given distribution function. Contains the following key/value pairs:

pos: an array of positions
vel: an array of velocities
mass: an array of masses

Each are astropy Quantities, with shape (Np,3).

Return type:

dict

Note

Up to at least galpy v1.7, galpy df objects don’t fully incorporate astropy Quantity inputs and outputs. Therefore, it is recommended that you define any relevant potential using the default ro=8, vo=220 unit system, turn off physical output for the potential, create the df object from the potential, use from_galpy_df() to create the ICs, and then turn physical output back on again afterwards if needed (see Example).

Example

Sample an NFW halo with scale radius 20 kpc, scale amplitude 2x10:sup:11 solar masses, and a maximum radius of 1 Mpc with 10,000 particles:

from astropy import units as u
from galpy import potential, df

NFWamp = 2e11 * u.Msun
NFWrs = 20 * u.kpc
ro = 8.
vo = 220.
rmax = 1 * u.Mpc
rmax_over_ro = (rmax/(ro*u.kpc)).to(1).value
Nhalo = 10000
NFWpot = potential.NFWPotential(amp=NFWamp, a=NFWrs)
NFWmass = potential.mass(NFWpot, rmax)
potential.turn_physical_off(NFWpot)
NFWdf = df.isotropicNFWdf(pot=NFWpot, rmax=rmax_over_ro)

halo_IC = IC.from_galpy_df(NFWdf, N=Nhalo, totmass=NFWmass)

potential.turn_physical_on(NFWpot)   # optional if needed later

Raises:: ExternalPackageException – If galpy could not be imported.

static from_pyn_snap(pynsnap)[source]¶

Turn a pynbody SimSnap into a set of GravHopper initial conditions.

Parameters:

pynsnap (SimSnap) – pynbody snapshot to convert

Returns:

IC – Properties of new particles to add, which sample the given distribution function. Contains the following key/value pairs:

pos: an array of positions
vel: an array of velocities
mass: an array of masses

Each are astropy Quantities, with shape (Np,3).

Return type:

dict

Note

The IC will try to use equivalent units to what are in the SimSnap.

Raises:: ExternalPackageException – If pynbody could not be imported.

unitconverter¶

Converts between Astropy-style and Pynbody-style units, which are not directly compatible.

gravhopper.unitconverter.astropy_to_pynbody(apunit)[source]¶

Return a pynbody unit equivalent to the astropy unit input.

Parameters:: apunit (astropy.unit.Unit) – Astropy unit to convert
Returns:: pynunit – Equivalent Pynbody unit. If the input unit is a composite, the function will attempt to compose it in the same way, but sometimes that is not possible; however, it will always be equivalent numerically.
Return type:: pynbody.units.Unit
Raises:: pynbody.UnitsException – If it cannot successfully convert the unit

gravhopper.unitconverter.pynbody_to_astropy(pynunit)[source]¶

Return an astropy unit equivalent to the pynbody unit input.

Parameters:: pynunit (pynbody.units.Unit) – Pynbody unit to convert
Returns:: apunit – Equivalent astropy unit. If the input unit is a composite, the function will attempt to compose it in the same way, but sometimes that is not possible; however, it will always be equivalent numerically.
Return type:: astropy.unit.Unit
Raises:: pynbody.UnitsException – If it cannot successfully convert the unit

jbgrav¶

The jbgrav module contains the functions that calculate the N-body force within a simulation. They are used internally within Simulation, but can also be imported and called on their own, and can also calculate potentials:

from gravhopper import jbgrav
from astropy import units as u
from numpy import np
sun = {'pos':[0,0,0]*u.au, 'vel':[0,0,0]*u.km/u.s, 'mass':[1]*u.Msun}
earth = {'pos':[1,0,0]*u.au, 'vel':[0,29.8,0]*u.km/u.s, 'mass':[1]*u.Mearth}
solarsystem = {'pos':np.vstack((sun['pos'],earth['pos'])),
    'vel':np.vstack((sun['vel'],earth['vel'])),
    'mass':np.vstack((sun['mass'],earth['mass']))}
accelerations, potentials = jbgrav.direct_summation(solarsystem, 0.01*u.au, calc_potential=True)

Functions for calculating the N-body force from a set of particles. Calls the C extension functions to do the calculation.

gravhopper.jbgrav.direct_summation(snap, eps, calc_force=True, calc_potential=False)[source]¶

Calculate the gravitational acceleration on every particle in the simulation snapshot from every other particle in the snapshot using direct summation. Uses an external C module for speed.

Parameters:

snap (dict) – snap[‘pos’] must be an (Np,3) array of astropy Quantities for the positions of the particles, and snap[‘mass’] must be an (Np) array of astropy Quantities of the masses of each particle.
eps (Quantity) – Gravitational softening length.
calc_force (boolean) – True if the force should be calculated. Can be set to False to only calculate the potential. (default: True)
calc_potential (boolean) – True if the potential should be calculated. (default: False)

Returns:

acceleration_potential – Return type depends on the values of calc_force and calc_potential.

If calc_force=True then the accelerations are returned as an (Np,3) numpy array of the acceleration vector calculated for each particle.

If calc_potential=True then the potential energies are returned as an Np element numpy array of the potential energies calculated for each particle.

If both are set, then the return is a tuple (accelerations, potentials,).

If neither is set, the function does no work and returns None.

Return type:

array or tuple

gravhopper.jbgrav.direct_summation_position(snap, force_pos, eps, calc_force=True, calc_potential=False)[source]¶

Calculate the gravitational acceleration at the specified locations from every particle in the snapshot using direct summation. Uses an external C module for speed.

direct_summation(snap, eps) should be equivalent to but twice as slow as direct_summation_position(snap, snap[‘pos’], eps).

Parameters:

snap (dict) – snap[‘pos’] must be an (Np,3) array of astropy Quantities for the positions of the particles, and snap[‘mass’] must be an (Np) array of astropy Quantities of the masses of each particle.
force_pos (array) – An (N,3) array of astropy Quantities for the positions at which the force will be calculated. If any position exactly matches a particle position from the snapshot, the self-force from the particle at its own location is taken to be zero.
eps (Quantity) – Gravitational softening length.
calc_force (boolean) – True if the force should be calculated. Can be set to False to only calculate the potential. (default: True)
calc_potential (boolean) – True if the potential should be calculated. (default: False)

Returns:

acceleration_potential – Return type depends on the values of calc_force and calc_potential.

If calc_force=True then the accelerations are returned as an (N,3) numpy array of the acceleration vector calculated for each position.

If calc_potential=True then the potential energies are returned as an N element numpy array of the potential energies calculated for each position.

If both are set, then the return is a tuple (accelerations, potentials,).

If neither is set, the function does no work and returns None.

Return type:

array or tuple

gravhopper.jbgrav.tree_force(snap, eps, theta=0.7, calc_force=True, calc_potential=False)[source]¶

Calculate the gravitational acceleration on every particle in the simulation snapshot from every other particle in the snapshot using a Barnes-Hut tree. Uses an external C module for speed.

Parameters:

snap (dict) – snap[‘pos’] must be an (Np,3) array of astropy Quantities for the positions of the particles, and snap[‘mass’] must be an (Np) array of astropy Quantities of the masses of each particle.
eps (Quantity) – Gravitational softening length.
theta (float) – Opening angle in radians. (default: 0.7)
calc_force (boolean) – True if the force should be calculated. Can be set to False to only calculate the potential. (default: True)
calc_potential (boolean) – True if the potential should be calculated. (default: False)

Returns:

acceleration_potential – Return type depends on the values of calc_force and calc_potential.

If calc_force=True then the accelerations are returned as an (Np,3) numpy array of the acceleration vector calculated for each particle.

If calc_potential=True then the potential energies are returned as an Np element numpy array of the potential energies calculated for each particle.

If both are set, then the return is a tuple (accelerations, potentials,).

If neither is set, the function does no work and returns None.

Return type:

array or tuple

gravhopper.jbgrav.tree_force_position(snap, force_pos, eps, theta=0.7, calc_force=True, calc_potential=False)[source]¶

Calculate the gravitational acceleration at the specified locations from every particle in the snapshot using a Barnes-Hut tree. Uses an external C module for speed.

Parameters:

snap (dict) – snap[‘pos’] must be an (Np,3) array of astropy Quantities for the positions of the particles, and snap[‘mass’] must be an (Np) array of astropy Quantities of the masses of each particle.
force_pos (array) – An (N,3) array of astropy Quantities for the positions at which the force will be calculated. If any position exactly matches a particle position from the snapshot, the self-force from the particle at its own location is taken to be zero.
eps (Quantity) – Gravitational softening length.
theta (float) – Opening angle in radians. (default: 0.7)
calc_force (boolean) – True if the force should be calculated. Can be set to False to only calculate the potential. (default: True)
calc_potential (boolean) – True if the potential should be calculated. (default: False)

Returns:

acceleration_potential – Return type depends on the values of calc_force and calc_potential.

If calc_force=True then the accelerations are returned as an (N,3) numpy array of the acceleration vector calculated for each position.

If calc_potential=True then the potential energies are returned as an N element numpy array of the potential energies calculated for each position.

If both are set, then the return is a tuple (accelerations, potentials,).

If neither is set, the function does no work and returns None.

Return type:

array or tuple

Exceptions¶

exception gravhopper.gravhopper.GravHopperException[source]¶: Parent class for all error exceptions.

exception gravhopper.gravhopper.UninitializedSimulationException[source]¶: Exception for trying to run a simulation without any initial conditions.

exception gravhopper.gravhopper.ICException(msg)[source]¶: Exception for trying to add ICs that don’t make sense.

exception gravhopper.gravhopper.UnknownAlgorithmException[source]¶: Exception for using an unknown N-body algorithm name.

exception gravhopper.gravhopper.ExternalPackageException(msg)[source]¶: Exception for trying to call a pynbody/galpy/gala/agama function when not using them.

Reference¶

Simulation¶

IC¶

unitconverter¶

jbgrav¶

Exceptions¶

Table of Contents

Previous topic

This Page