Class or Metaclass Design for Astrodynamics Engine

Question

Gurus out there:

The differential equations for modeling spacecraft motion can be described in terms of a collection of acceleration terms:

d2r/dt2 =  a0 + a1 + a2 + ... + an

Normally a0 is the point mass acceleration due to a body (a0 = -mu * r/r^3); the "higher order" terms can be due to other planets, solar radiation pressure, thrust, etc.

I'm implementing a collection of algorithms meant to work on this sort of system. I will start with Python for design and prototyping, then I will move on to C++ or Fortran 95.

I want to design a class (or metaclass) which will allow me to specify the different acceleration terms for a given instance, something along the lines of:

# please notice this is meant as "pseudo-code"
def some_acceleration(t):
    return (1*t, 2*t, 3*t)

def some_other_acceleration(t):
    return (4*t, 5*t, 6*t)

S = Spacecraft()
S.Acceleration += someacceleration + some_other_acceleration

In this case, the instance S would default to, say, two acceleration terms and I would add a the other two terms I want: some acceleration and some_other_acceleration; they return a vector (here represented as a triplet). Notice that in my "implementation" I've overloaded the + operator.

This way the algorithms will be designed for an abstract "spacecraft" and all the actual force fields will be provided on a case-by-case basis, allowing me to work with simplified models, compare modeling methodologies, etc.

How would you implement a class or metaclass for handling this?

I apologize for the rather verbose and not clearly explained question, but it is a bit fuzzy in my brain.

Thanks.

Answer 1

I would use instead some library which can work with vectors (in python, try numpy) and represent the acceleration as vector. Then, you are not reinventing the wheel, the + operator works like you wanted. Please correct me if I misunderstood your problem.

Answer 2

For those who would like to avoid numpy and do this in pure python, this may give you a few good ideas. I'm sure there are disadvantages and flaws to this little skit also. The "operator" module speeds up your math calculations as they are done with c functions:

from operator import sub, add, iadd, mul
import copy

class Acceleration(object):
   def __init__(self, x, y, z):
      super(Acceleration, self).__init__()
      self.accel = [x, y , z]
      self.dimensions = len(self.accel)

   @property
   def x(self):
      return self.accel[0]

   @x.setter
   def x(self, val):
      self.accel[0] = val


   @property
   def y(self):
      return self.accel[1]

   @y.setter
   def y(self, val):
      self.accel[1] = val

   @property
   def z(self):
      return self.accel[2]

   @z.setter
   def z(self, val):
      self.accel[2] = val

   def __iadd__(self, other):
      for x in xrange(self.dimensions):
         self.accel[x] = iadd(self.accel[x], other.accel[x])
      return self

   def __add__(self, other):
      newAccel = copy.deepcopy(self)
      newAccel += other
      return newAccel

   def __str__(self):
      return "Acceleration(%s, %s, %s)" % (self.accel[0], self.accel[1], self.accel[2])

   def getVelocity(self, deltaTime):
      return Velocity(mul(self.accel[0], deltaTime), mul(self.accel[1], deltaTime), mul(self.accel[2], deltaTime))

class Velocity(object):
   def __init__(self, x, y, z):
      super(Velocity, self).__init__()
      self.x = x
      self.y = y
      self.z = z

   def __str__(self):
      return "Velocity(%s, %s, %s)" % (self.x, self.y, self.z)

if __name__ == "__main__":
   accel = Acceleration(1.1234, 2.1234, 3.1234)
   accel += Acceleration(1, 1, 1)
   print accel

   accels = []
   for x in xrange(10):
      accel += Acceleration(1.1234, 2.1234, 3.1234)

   vel = accel.getVelocity(2)
   print "Velocity of object with acceleration %s after one second:" % (accel)
   print vel

prints the following:

Acceleration(2.1234, 3.1234, 4.1234)

Velocity of object with acceleration Acceleration(13.3574, 24.3574, 35.3574) after one second: Velocity(26.7148, 48.7148, 70.7148)

You can get fancy for faster calculations:

def getFancyVelocity(self, deltaTime):
   from itertools import repeat
   x, y, z = map(mul, self.accel, repeat(deltaTime, self.dimensions))
   return Velocity(x, y, z)

Answer 3

Are you asking how to store an arbitrary number of acceleration sources for a spacecraft class?

Can't you just use an array of functions? (Function pointers when you get to c++)

i.e:

#pseudo Python
class Spacecraft
    terms = []
    def accelerate(t):
       a = (0,0,0)
       for func in terms:
         a+= func(t)


s = Spacecraft
s.terms.append(some_acceleration)
s.terms.append(some_other_acceleration)
ac = s.accelerate(t)