Source code for rebound.particle

from ctypes import Structure, c_double, c_int, byref, memmove, sizeof, c_uint32, c_uint, c_ulong
from . import clibrebound
import math
import ctypes.util
import rebound
import sys
import random

__all__ = ["Particle"]

def notNone(a):
    """
    Returns True if array a contains at least one element that is not None. Returns False otherwise.
    """
    return a.count(None) != len(a)

[docs]class Particle(Structure): """ The main REBOUND particle data structure. This is an abstraction of the reb_particle structure in C. The Particle fields are set at the end of simulation.py to avoid circular references. Attributes ---------- x, y, z : float Particle positions vx, vy, vz : float Particle velocities ax, ay, az : float Particle accelerations m : float Particle mass r : float Particle radius lastcollision : float Last time the particle had a physical collision (if checking for collisions) c : c_void_p (C void pointer) Pointer to the cell the particle is currently in (if using tree code) hash : c_uint32 Particle hash (permanent identifier for the particle) ap : c_void_p (C void pointer) Pointer to additional parameters one might want to add to particles _sim : POINTER(rebound.Simulation) Internal pointer to the parent simulation (used in C version of REBOUND) a, e, inc, Omega, omega, f : float (Kepler Elements) Semi-major axis, eccentricity, inclination, longitude of the ascending node, argument of periapsis, and true anomaly respectively. The Keplerian Elements are in Jacobi coordinates (with mu = G*Minc, where Minc is the total mass from index 0 to the particle's index, inclusive). """ def __str__(self): """ Returns a string with the position and velocity of the particle. """ return "<rebound.Particle object, m=%s x=%s y=%s z=%s vx=%s vy=%s vz=%s>"%(self.m,self.x,self.y,self.z,self.vx,self.vy,self.vz) __repr__ = __str__
[docs] def __init__(self, simulation=None, particle=None, m=None, x=None, y=None, z=None, vx=None, vy=None, vz=None, primary=None, a=None, P=None, e=None, inc=None, Omega=None, omega=None, pomega=None, f=None, M=None, l=None, theta=None, T=None, r=None, date=None, variation=None, variation2=None, h=None, k=None, ix=None, iy=None, hash=0, jacobi_masses=False): """ Initializes a Particle structure. Rather than explicitly creating a Particle structure, users may use the ``add()`` member function of a Simulation instance, which will both create a Particle and then add it to the simulation with one function call. This function accepts either cartesian positions and velocities, classical orbital elements together with the reference Particle (the primary), as well as orbital parameters defined by Pal (2009). For convenience, optional keywords that are not passed default to zero (mass, cartesian and orbital elements). Whenever initializing a particle from orbital elements, one must specify either the semimajor axis or the period of the orbit. For classical orbital paramerers, one can specify the longitude of the ascending node by passing Omega, to specify the pericenter one can pass either omega or pomega (not both), and for the longitude/anomaly one can pass one of f, M, l or theta. See ipython_examples/OrbitalElements.ipynb for examples. See also Murray & Dermott Solar System Dynamics for formal definitions of angles in orbital mechanics. All angles should be specified in radians. Parameters ---------- simulation : Simulation Simulation instance associated with this particle (Required if passing orbital elements or setting up a variation). particle : Particle, optional If a particle is passed, a copy of that particle is returned. If a variational particle is initialized, then ``particle`` is original particle that will be varied. m : float Mass (Default: 0) x, y, z : float Positions in Cartesian coordinates (Default: 0) vx, vy, vz : float Velocities in Cartesian coordinates (Default: 0) primary : Particle Primary body for converting orbital elements to cartesian (Default: center of mass of the particles in the passed simulation, i.e., this will yield Jacobi coordinates as one progressively adds particles) a : float Semimajor axis (a or P required if passing orbital elements) P : float Orbital period (a or P required if passing orbital elements) e : float Eccentricity (Default: 0) inc : float Inclination (Default: 0) Omega : float Longitude of ascending node (Default: 0) omega : float Argument of pericenter (Default: 0) pomega : float Longitude of pericenter (Default: 0) f : float True anomaly (Default: 0) M : float Mean anomaly (Default: 0) l : float Mean longitude (Default: 0) theta : float True longitude (Default: 0) T : float Time of pericenter passage h : float h variable, see Pal (2009) for a definition (Default: 0) k : float k variable, see Pal (2009) for a definition (Default: 0) ix : float ix variable, see Pal (2009) for a definition (Default: 0) iy : float iy variable, see Pal (2009) for a definition (Default: 0) r : float Particle radius (only used for collisional simulations) date : string For consistency with adding particles through horizons. Not used here. variation : string (Default: None) Set this string to the name of an orbital parameter to initialize the particle as a variational particle. Can be one of the following: m, a, e, inc, omega, Omega, f, k, h, lambda, ix, iy. variation2 : string (Default: None) Set this string to the name of a second orbital parameter to initialize the particle as a second order variational particle. Only used for second order variational equations. Can be one of the following: m, a, e, inc, omega, Omega, f, k, h, lambda, ix, iy. hash : c_uint32 Unsigned integer identifier for particle. Can pass an integer directly, or a string that will be converted to a hash. User is responsible for assigning unique hashes. jacobi_masses: bool Whether to use jacobi primary mass in orbit initialization. Particle mass will still be set to physical value (Default: False) Examples -------- >>> sim = rebound.Simulation() >>> sim.add(m=1.) >>> p1 = rebound.Particle(simulation=sim, m=0.001, a=0.5, e=0.01) >>> p2 = rebound.Particle(simulation=sim, m=0.0, x=1., vy=1.) >>> p3 = rebound.Particle(simulation=sim, m=0.001, a=1.5, h=0.1, k=0.2, l=0.1) >>> p4 = rebound.Particle(simulation=sim, m=0.001, a=1.5, omega="uniform") # omega will be a random number between 0 and 2pi """ if Omega == "uniform": Omega = random.vonmisesvariate(0.,0.) if omega == "uniform": omega = random.vonmisesvariate(0.,0.) if pomega == "uniform": pomega = random.vonmisesvariate(0.,0.) if f == "uniform": f = random.vonmisesvariate(0.,0.) if M == "uniform": M = random.vonmisesvariate(0.,0.) if l == "uniform": l = random.vonmisesvariate(0.,0.) if theta == "uniform": theta = random.vonmisesvariate(0.,0.) self.hash = hash # set via the property, which checks for type if variation: if primary is None: primary = simulation.particles[0] # Find particle to differenciate lc = locals().copy() del lc["self"] del lc["variation"] del lc["variation2"] if particle is None: particle = Particle(**lc) # First or second order? if variation and variation2: variation_order = 2 else: variation_order = 1 # Shortcuts for variable names if variation == "l": variation = "lambda" if variation2 == "l": variation2 = "lambda" if variation == "i": variation = "inc" if variation2 == "i": variation2 = "inc" variationtypes = ["m","a","e","inc","omega","Omega","f","k","h","lambda","ix","iy"] if variation_order==1: if variation in variationtypes: method = getattr(clibrebound, 'reb_derivatives_'+variation) method.restype = Particle p = method(c_double(simulation.G), primary, particle) else: raise ValueError("Variational particles can only be initializes using the derivatives with respect to one of the following: %s."%", ".join(variationtypes)) elif variation_order==2: if variation in variationtypes and variation2 in variationtypes: # Swap variations if needed vi1 = variationtypes.index(variation) vi2 = variationtypes.index(variation2) if vi2 < vi1: variation, variation2 = variation2, variation method = getattr(clibrebound, 'reb_derivatives_'+variation+'_'+variation2) method.restype = Particle p = method(c_double(simulation.G), primary, particle) else: raise ValueError("Variational particles can only be initializes using the derivatives with respect to one of the following: %s."%", ".join(variationtypes)) else: raise ValueError("Variational equations beyond second order are not implemented.") self.m = p.m self.x = p.x self.y = p.y self.z = p.z self.vx = p.vx self.vy = p.vy self.vz = p.vz return if particle is not None: memmove(byref(self), byref(particle), sizeof(self)) return cart = [x,y,z,vx,vy,vz] orbi = [primary,a,P,e,inc,Omega,omega,pomega,f,M,l,theta,T] pal = [h,k,ix,iy] self.ax = 0. self.ay = 0. self.az = 0. if m is None: self.m = 0. else: self.m = m if r is None: self.r = 0. else: self.r = r self.lastcollision = 0. self.c = None self.ap = None if notNone([e,inc,omega,pomega,Omega,M,f,theta,T]) and notNone(pal): raise ValueError("You cannot mix Pal coordinates (h,k,ix,iy) with the following orbital elements: e,inc,Omega,omega,pomega,f,M,theta,T. If a longitude/anomaly is needed in Pal coordinates, use l.") if notNone(cart) and notNone(orbi): raise ValueError("You cannot pass cartesian coordinates and orbital elements (and/or primary) at the same time.") if notNone(orbi): if simulation is None: raise ValueError("Need to specify simulation when initializing particle with orbital elements.") if primary is None: clibrebound.reb_get_com.restype = Particle primary = clibrebound.reb_get_com(byref(simulation)) # this corresponds to adding in Jacobi coordinates if jacobi_masses is True: interior_mass = 0 for p in simulation.particles: interior_mass += p.m # orbit conversion uses mu=G*(p.m+primary.m) so set prim.m=Mjac-m so mu=G*Mjac primary.m = simulation.particles[0].m*(self.m + interior_mass)/interior_mass - self.m if a is None and P is None: raise ValueError("You need to pass either a semimajor axis or orbital period to initialize the particle using orbital elements.") if a is not None and P is not None: raise ValueError("You can pass either the semimajor axis or orbital period, but not both.") if a is None: a = (P**2*simulation.G*(primary.m + self.m)/(4.*math.pi**2))**(1./3.) if notNone(pal): # Pal orbital parameters if h is None: h = 0. if k is None: k = 0. if l is None: l = 0. if ix is None: ix = 0. if iy is None: iy = 0. if((ix*ix + iy*iy) > 4.0): raise ValueError("Passed (ix, iy) coordinates are not valid, squared sum exceeds 4.") clibrebound.reb_tools_pal_to_particle.restype = Particle p = clibrebound.reb_tools_pal_to_particle(c_double(simulation.G), primary, c_double(self.m), c_double(a), c_double(l), c_double(k), c_double(h), c_double(ix), c_double(iy)) else: # Normal orbital parameters if e is None: e = 0. if inc is None: inc = 0. if Omega is None: # we require that Omega be passed if you want to specify longitude of node Omega = 0. pericenters = [omega, pomega] # Need omega for C function. Can specify it either directly or through pomega indirectly. numNones = pericenters.count(None) if numNones == 0: raise ValueError("Can't pass both omega and pomega") if numNones == 2: # Neither passed. Default to 0. omega = 0. if numNones == 1: if pomega is not None: # Only have to find omega is pomega was passed if math.cos(inc) > 0: # inc is in range [-pi/2, pi/2] (prograde), so pomega = Omega + omega omega = pomega - Omega else: omega = Omega - pomega # for retrograde orbits, pomega = Omega - omega longitudes = [f,M,l,theta,T] # can specify longitude through any of these four. Need f for C function. numNones = longitudes.count(None) if numNones < 4: raise ValueError("Can only pass one longitude/anomaly in the set [f, M, l, theta, T]") if numNones == 5: # none of them passed. Default to 0. f = 0. if numNones == 4: # Only one was passed. if f is None: # Only have to work if f wasn't passed. if theta is not None: # theta is next easiest if math.cos(inc) > 0: # for prograde orbits, theta = Omega + omega + f f = theta - Omega - omega else: f = Omega - omega - theta # for retrograde, theta = Omega - omega - f else: # Either M, l, or T was passed. Will need to find M first (if not passed) to find f if l is not None: if math.cos(inc) > 0: # for prograde orbits, l = Omega + omega + M M = l - Omega - omega else: M = Omega - omega - l # for retrograde, l = Omega - omega - M else: if T is not None: # works for both elliptical and hyperbolic orbits # TODO: has accuracy problems for M=n*(t-T) << 1 n = (simulation.G*(primary.m+self.m)/abs(a**3))**0.5 M = n*(simulation.t - T) clibrebound.reb_tools_M_to_f.restype = c_double f = clibrebound.reb_tools_M_to_f(c_double(e), c_double(M)) err = c_int() clibrebound.reb_tools_orbit_to_particle_err.restype = Particle p = clibrebound.reb_tools_orbit_to_particle_err(c_double(simulation.G), primary, c_double(self.m), c_double(a), c_double(e), c_double(inc), c_double(Omega), c_double(omega), c_double(f), byref(err)) if err.value == 1: raise ValueError("Can't set e exactly to 1.") if err.value == 2: raise ValueError("Eccentricity must be greater than or equal to zero.") if err.value == 3: raise ValueError("Bound orbit (a > 0) must have e < 1.") if err.value == 4: raise ValueError("Unbound orbit (a < 0) must have e > 1.") if err.value == 5: raise ValueError("Unbound orbit can't have f beyond the range allowed by the asymptotes set by the hyperbola.") if err.value == 6: raise ValueError("Primary has no mass.") self.x = p.x self.y = p.y self.z = p.z self.vx = p.vx self.vy = p.vy self.vz = p.vz else: if x is None: x = 0. if y is None: y = 0. if z is None: z = 0. if vx is None: vx = 0. if vy is None: vy = 0. if vz is None: vz = 0. self.x = x self.y = y self.z = z self.vx = vx self.vy = vy self.vz = vz
[docs] def copy(self): """ Returns a deep copy of the particle. The particle is not added to any simulation by default. """ np = Particle() memmove(byref(np), byref(self), sizeof(self)) return np
[docs] def calculate_orbit(self, primary=None, G=None): """ Returns a rebound.Orbit object with the keplerian orbital elements corresponding to the particle around the passed primary (rebound.Particle) If no primary is passed, defaults to Jacobi coordinates (with mu = G*Minc, where Minc is the total mass from index 0 to the particle's index, inclusive). Examples -------- >>> sim = rebound.Simulation() >>> sim.add(m=1.) >>> sim.add(x=1.,vy=1.) >>> orbit = sim.particles[1].calculate_orbit(sim.particles[0]) >>> print(orbit.e) # gives the eccentricity Parameters ---------- primary : rebound.Particle Central body (Optional. Default uses Jacobi coordinates) G : float Gravitational constant (Optional. Default takes G from simulation in which particle is in) Returns ------- A rebound.Orbit object """ if not self._sim: # Particle not in a simulation if primary is None: raise ValueError("Particle does not belong to any simulation and no primary given. Cannot calculate orbit.") if G is None: raise ValueError("Particle does not belong to any simulation and G not given. Cannot calculate orbit.") else: G = c_double(G) else: # First check whether this is particles[0] clibrebound.reb_get_particle_index.restype = c_int index = clibrebound.reb_get_particle_index(byref(self)) # first check this isn't particles[0] if index == 0 and primary is None: raise ValueError("Orbital elements for particle[0] not implemented unless primary is provided") if primary is None: # Use default, i.e., Jacobi coordinates clibrebound.reb_get_jacobi_com.restype = Particle # now return jacobi center of mass primary = clibrebound.reb_get_jacobi_com(byref(self)) G = c_double(self._sim.contents.G) err = c_int() clibrebound.reb_tools_particle_to_orbit_err.restype = rebound.Orbit o = clibrebound.reb_tools_particle_to_orbit_err(G, self, primary, byref(err)) if err.value == 1: raise ValueError("Primary has no mass.") if err.value == 2: raise ValueError("Particle and primary positions are the same.") return o
[docs] def sample_orbit(self, Npts=100, primary=None, trailing=True, timespan=None, useTrueAnomaly=None): """ Returns a nested list of xyz positions along the osculating orbit of the particle. If primary is not passed, returns xyz positions along the Jacobi osculating orbit (with mu = G*Minc, where Minc is the total mass from index 0 to the particle's index, inclusive). Parameters ---------- Npts : int, optional Number of points along the orbit to return (default: 100) primary : rebound.Particle, optional Primary to use for the osculating orbit (default: Jacobi center of mass) trailing: bool, optional Whether to return points stepping backwards in time (True) or forwards (False). (default: True) timespan: float, optional Return points (for the osculating orbit) from the current position to timespan (forwards or backwards in time depending on trailing keyword). Defaults to the orbital period for bound orbits, and to the rough time it takes the orbit to move by the current distance from the primary for a hyperbolic orbit. Implementation currently only supports this option if useTrueAnomaly=False. useTrueAnomaly: bool, optional Will sample equally spaced points in true anomaly if True, otherwise in mean anomaly. Latter might be better for hyperbolic orbits, where true anomaly can stay near the limiting value for a long time, and then switch abruptly at pericenter. (Default: True for bound orbits, False for unbound orbits) """ pts = [] if primary is None: primary = self.jacobi_com o = self.calculate_orbit(primary=primary) if timespan is None: if o.a < 0.: # hyperbolic orbit timespan = 2*math.pi*o.d/o.v # rough time to cross display box else: timespan = o.P lim_phase = abs(o.n)*timespan # n is negative for hyperbolic orbits if trailing is True: lim_phase *= -1 # sample phase backwards from current value phase = [lim_phase*i/(Npts-1) for i in range(Npts)] if useTrueAnomaly is None: if o.a <0.: useTrueAnomaly = False else: useTrueAnomaly = True for i,ph in enumerate(phase): if useTrueAnomaly is True: newp = Particle(a=o.a, f=o.f+ph, inc=o.inc, omega=o.omega, Omega=o.Omega, e=o.e, m=self.m, primary=primary, simulation=self._sim.contents) else: newp = Particle(a=o.a, M=o.M+ph, inc=o.inc, omega=o.omega, Omega=o.Omega, e=o.e, m=self.m, primary=primary, simulation=self._sim.contents) pts.append(newp.xyz) return pts
# Simple operators for particles. def __add__(self, other): if not isinstance(other, Particle): return NotImplemented c = self.copy() return c.__iadd__(other) def __iadd__(self, other): if not isinstance(other, Particle): return NotImplemented clibrebound.reb_particle_iadd(byref(self), byref(other)) return self def __sub__(self, other): if not isinstance(other, Particle): return NotImplemented c = self.copy() return c.__isub__(other) def __isub__(self, other): if not isinstance(other, Particle): return NotImplemented clibrebound.reb_particle_isub(byref(self), byref(other)) return self def __mul__(self, other): try: other = float(other) except: return NotImplemented c = self.copy() return c.__imul__(other) def __imul__(self, other): try: other = float(other) except: return NotImplemented clibrebound.reb_particle_imul(byref(self), c_double(other)) return self def __rmul__(self, other): try: other = float(other) except: return NotImplemented c = self.copy() return c.__imul__(other) def __div__(self, other): return self.__truediv__(other) def __idiv__(self, other): return self.__itruediv__(other) def __truediv__(self, other): try: other = float(other) except: return NotImplemented c = self.copy() if other==0.: raise ZeroDivisionError return c.__imul__(1./other) def __itruediv__(self, other): try: other = float(other) except: return NotImplemented if other==0.: raise ZeroDivisionError return self.__imul__(1./other) @property def index(self): clibrebound.reb_get_particle_index.restype = c_int return clibrebound.reb_get_particle_index(byref(self)) @property def xyz(self): """ Get or set the xyz position coordinates of the particle. """ return [self.x, self.y, self.z] @xyz.setter def xyz(self, value): if len(value)!=3: raise AttributeError("Can only set xyz positions to array with length 3") self.x = float(value[0]) self.y = float(value[1]) self.z = float(value[2]) @property def vxyz(self): """ Get or set the xyz velocity coordinates of the particle. """ return [self.vx, self.vy, self.vz] @vxyz.setter def vxyz(self, value): if len(value)!=3: raise AttributeError("Can only set xyz velocities to array with length 3") self.vx = float(value[0]) self.vy = float(value[1]) self.vz = float(value[2]) @property def d(self): return self.calculate_orbit().d @property def v(self): return self.calculate_orbit().v @property def h(self): return self.calculate_orbit().h @property def P(self): return self.calculate_orbit().P @property def n(self): return self.calculate_orbit().n @property def a(self): return self.calculate_orbit().a @property def rhill(self): return self.calculate_orbit().rhill @property def e(self): return self.calculate_orbit().e @property def inc(self): return self.calculate_orbit().inc @property def Omega(self): return self.calculate_orbit().Omega @property def omega(self): return self.calculate_orbit().omega @property def pomega(self): return self.calculate_orbit().pomega @property def f(self): return self.calculate_orbit().f @property def M(self): return self.calculate_orbit().M @property def l(self): return self.calculate_orbit().l @property def theta(self): return self.calculate_orbit().theta @property def T(self): return self.calculate_orbit().T @property def orbit(self): return self.calculate_orbit() @property def jacobi_com(self): clibrebound.reb_get_jacobi_com.restype = Particle return clibrebound.reb_get_jacobi_com(byref(self)) @property def hash(self): """ Get or set the particle's hash. If set to a string, the corresponding integer hash is calculated. """ return c_uint32(self._hash) @hash.setter def hash(self, value): PY3 = sys.version_info[0] == 3 hash_types = c_uint32, c_uint, c_ulong if PY3: string_types = str, int_types = int, else: string_types = basestring, int_types = int, long, if isinstance(value, hash_types): self._hash = value.value elif isinstance(value, string_types): self._hash = rebound.hash(value).value elif isinstance(value, int_types): self._hash = value else: raise AttributeError("Hash must be set to an integer, a ctypes.c_uint32 or a string. See UniquelyIdentifyingParticles.ipynb ipython_example.")