"""Generators to simulate orbit system dynamics. The viewer classes can be found in QtOrbits."""
################################################################
# Written in 2003-2019 by Garth Zeglin <garthz@cmu.edu>
# To the extent possible under law, the author has dedicated all copyright
# and related and neighboring rights to this software to the public domain
# worldwide. This software is distributed without any warranty.
# You should have received a copy of the CC0 Public Domain Dedication along with this software.
# If not, see <http://creativecommons.org/publicdomain/zero/1.0/>.
################################################################
import math, logging
import numpy as np
# set up logger for module
log = logging.getLogger('orbits')
# filter out most logging; the default is NOTSET which passes along everything
# log.setLevel(logging.WARNING)
################################################################
[docs]class Orbits(object):
"""Simulation of a set of particles orbiting around fixed attractors."""
def __init__(self, channels=20):
# scalar constants
self.P = channels # Number of particles.
self.B = 3 # Number of fixed bodies.
self.g = 0.1 # Gravitational constant.
self.b = 0.01 # Viscous damping.
self.kC = 1000.0 # Particle collision stiffness.
self.bC = 0.1 # Particle collision damping.
self.limit = np.sqrt(3) # Radius of the spherical constraint limit.
self.Ft = 0.0 # Global applied tangential force field magnitude.
self.Fr = 0.0 # Global applied radial force field magnitude.
# q i s a (P,3) matrix of particle positions, initially uniformly spaced
# around a circle in the XY plane with periodic Z variation
iP = np.arange(self.P)
dtheta = 2*math.pi/self.P
self.q = np.stack((np.cos(dtheta*iP), np.sin(dtheta*iP), 0.2*np.sin(4*dtheta*iP))).T
# Rp is a (P,) vector of particle radii
self.Rp = 0.05 * np.ones(self.P)
# qd is a (P,3) matrix of particle velocities, initially directed tangentially around the circle
rot = np.array([[0, -1, 0], [1, 0, 0], [0, 0, 1]])
self.qd = 0.5 * np.dot(self.q, rot)
# qB is a (B, 3) matrix of fixed body attractor locations, uniformly spaced around the circle
iB = np.arange(self.B)
dtheta = 2*math.pi/self.B
self.qB = 0.5*np.stack((np.cos(dtheta*iB), np.sin(dtheta*iB), 0*iB)).T
# Rp is a (B,) vector of body radii
self.Rb = 0.2 * np.ones(self.B)
# Boundary collision flags are used to detect events.
self.b_colliding = np.zeros((self.P), dtype=bool)
self.b_collision_listener = None
# elapsed integrator time
self.t = 0.0
return
[docs] def set_applied_force(self, tangential, radial):
self.Ft = tangential
self.Fr = radial
[docs] def set_damping(self, b):
self.b = b
[docs] def set_gravity(self, g):
self.b = g
[docs] def connect_boundary_collision_listener(self, callback):
self.b_collision_listener = callback
[docs] def get_particles_and_bodies(self):
return self.q.copy(), self.qB.copy()
[docs] def update_for_interval(self, interval):
"""Run the simulator for the given interval, which may include one or more integration steps."""
while interval > 0.0:
dt = min(interval, 0.001)
interval -= dt
self._step(dt)
def _step(self, dt):
# ---- Calculate the inverse-square law gravity forces -------
# Compute the pairwise vectors from each particle q[:] to each body
# qB[:] using numpy matrix broadcasting with shapes as follows:
# pairs(B, P, 3) = qB(B,1,3) - q(1,P,3)
pairs = self.qB.reshape((self.B,1,3)) - self.q.reshape((1,self.P,3))
# Compute the pairwise center distances using shapes as follows:
# dist (B, P) = sqrt(sum((pairs (B, P, 3) * pairs.(B, P, 3), axis=2))
dist = np.sqrt(np.sum(pairs*pairs, axis=2))
# Normalize the pairwise displacement vectors to unit vectors from particles to bodies:
# unit (B, P, 3) = dist^-1 (B, P, 1) * pairs(B, P, 3)
unit = np.power(dist.reshape((self.B, self.P, 1)), -1) * pairs
# Calculate the pairwise gravity vectors, assuming unit masses for simplicity:
# F = G / R^2
# = G * R^(-2) * unit(pair)
# F (B, P, 3) = G * dist**(-2)(B, P, 1) * unit (B,P,3)
Fgrav = self.g * np.power(dist, -2).reshape((self.B, self.P, 1)) * unit
# --- Calculate the particle-body collision forces, if any ------
# Calculate the pairwise distance thresholds.
# min_dist (B, P) = Rb (B, 1) + Rp(1, P)
min_dist = self.Rb.reshape((self.B, 1)) + self.Rp.reshape((1, self.P))
# Calculate the body-particle interferences as a matrix of non-negative values.
# inter (B, P)
inter = np.maximum(min_dist - dist, 0)
# Apply a Hooke's Law proportional spring force (F = k x) along the
# body-particle normals if interfering.
# Fcoll (B, P, 3) = -kB * inter (B, P, 1) * unit
Fcoll = (-self.kC) * inter.reshape((self.B, self.P, 1)) * unit
# Apply pairwise collision damping.
# F = -b * (qd . unit)*unit for pairs in collision
# F = -b * (qdR*unit)
# qdR is a ((B, P)) matrix of the radial velocity component of each particle along the axis toward a body
qdR = np.sum(self.qd.reshape((1, self.P, 3)) * unit, axis=2)
Fcoll -= (self.bC * (inter > 0.0)).reshape((self.B, self.P, 1)) * (qdR.reshape((self.B, self.P, 1)) * unit)
# --- Calculate the per-particle forces --------------
# General viscous drag.
Fdamping = (-self.b) * self.qd
# Apply an elastic spherical boundary.
# Calculate the radial position of each particle.
# R (P) = sqrt(sum(q(P, 3) * q(P,3)), axis=2)
R = np.sqrt(np.sum(self.q * self.q, axis=1))
# Calculate the outward normals.
Qunit = np.power(R.reshape((self.P, 1)), -1) * self.q
# Calculate the particle-boundary interferences as a non-zero value.
bounds = np.maximum((R + self.Rp) - self.limit, 0)
# Apply a proportional spring force inward.
Fbound = (-self.kC) * bounds.reshape((self.P, 1)) * Qunit
# Detect collision events
colliding = bounds > 0.0
new_collisions = np.logical_and(colliding, np.logical_not(self.b_colliding))
self.b_colliding = colliding
if self.b_collision_listener is not None:
# iterate over the set of indicies of true elements
for i in np.nonzero(new_collisions)[0]:
self.b_collision_listener(i, np.linalg.norm(self.qd[i]))
# Apply global fields. These can be used to steer the whole system externally.
# Calculate the tangential and radial directions for each particle.
unit_r = Qunit
unit_t = np.dot(unit_r, np.array([[0, -1, 0], [1, 0, 0], [0, 0, 1]]))
Ffield = self.Fr * unit_r + self.Ft * unit_t
# ---- Calculate the particle dynamics ------------
# Sum all pairwise forces to a net force on each particle:
# forces(P, 3) = np.sum(Fgrav+Fcoll, axis=0) (B, P, 3)
forces = np.sum(Fgrav + Fcoll, axis=0)
# Add in per-particle forces
forces += Fdamping + Fbound + Ffield
# integrate one time step
self.q += self.qd * dt
self.qd += forces * dt
self.t += dt
return
################################################################
if __name__ == "__main__":
obj = Orbits()
obj.update_for_interval(0.040)