Class or Metaclass Design for Astrodynamics Engine

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.

From stackoverflow

• 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.

  Arrieta : It's not so much how to represent the force field (I would use a viable Vector package) but how to define the actual equations that model the force field.

• 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)
  

  Arrieta : Thanks for your solution, but it only provides means to add vectors and a method for approximating their integral. What I meant was a way to define the mathematical representation of the field (acceleration, in this case) but probably my question was not clear enough. I will rephrase it soon... Thanks again.

• 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)