Rotating Asteroid¶

This notebook demonstrates how godot can be used to model a system of a spacecraft and a spinning asteroid.

In [1]:

from godot.core import tempo
from godot.core.autodif import bridge as br
from godot.core import autodif as ad
from godot.model import frames, param, common, prop, common
from godot import cosmos
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

# optionally avoid verbose logging messages
import godot.core.util as util
util.suppressLogger()

The Problem¶

We are evolving a system of a spinning asteroid and a spacecraft.

The spacecraft will fly by the asteroid and under the influence of its gravitational field. The asteroid is modeled as an ellipsoid and is rotating in a precessing motion. On the surface of the asteroid there are several landmarks distributed. The idea is that the spacecraft takes several pictures of the asteroid as it flies by it and can identify the landmarks of the asteroid.

We will try and determine at which time which landmarks are visible to the spaceraft and compute the vector between each visible landmark and the spacecraft.

Parameters of the problem¶

The system consists of a spacecraft and an asteroid. We will track the following parameters:

1.) initial orientation of asteroid: $\alpha$, $\delta$, $\omega$

2.) initial angular velocity of asteroid: $\dot{\alpha}$, $\dot{\delta}$, $\dot{\omega}$

3.) density of asteroid: $\rho$

4.) initial position of spacecraft: $x$, $y$, $z$

5.) initial velocity of spacecraft: $\dot{x}$, $\dot{y}$, $\dot{z}$

We set up the parameters using a parameter book as they enable us to conveniently store, access and track partial derivatives.

For the initial conditions, we put the spacecraft ca. 130 km from the asteroid and give it small initial velocity of 0.0006 km/s.

The asteroid is initially aligned with the ICRF system and is spinning with an angular velocity correspoding to one period per day.

In [2]:

uniConfig = cosmos.util.load_yaml('universe.yml')
uni = cosmos.Universe(uniConfig)


alpha = uni.parameters.add('alpha', np.pi/4)
delta = uni.parameters.add('delta', 0.)
omega = uni.parameters.add('omega', 0.)
vAngular = np.array([0., 0.5, 1.])
vAngular = 2 * np.pi* vAngular / (86400 * np.linalg.norm(vAngular))

alphaDot = uni.parameters.add("vAlpha", vAngular[0])
deltaDot = uni.parameters.add("vDelta", vAngular[1])
omegaDot = uni.parameters.add("vOmega", vAngular[2])

rho = uni.parameters.add("density", 1e12)

posx = uni.parameters.add("posx", -125.)
posy = uni.parameters.add("posy", 50.)
posz = uni.parameters.add("posz", 0.)

velx = uni.parameters.add("velx", 6e-4)
vely = uni.parameters.add("vely", 0.)
velz = uni.parameters.add("velz", 0.)

uni.parameters.track(['posx', 'posy', 'posz', 'velx', 'vely', 'velz', 'vOmega'])

Adding points and axes to the frames system¶

Next, we set up a frames system and add a point for the spacecraft and the asteroid, as well as ICRF axes and axes for the asteroid.

In [3]:

tscale = tempo.TimeScale.TDB

fra = uni.frames
icrf = fra.axesId('ICRF')
astAxes = fra.addAxes('astAxes', tscale)

astPoint = fra.addPoint('asteroid', tscale)
scPoint = fra.addPoint('spacecraft', tscale)

t0 = tempo.XEpoch()

Parameterizing initial conditions for the asteroid¶

We need to parameterize the initial orientation of the asteroid in terms of quaternion components. For the initial conditions and the accelerations/torques, the propagator works with time evaluable objects. We can create them with a class that inherits from the common module.

In [4]:

class Quaternion(common.Vector7TimeEvaluable):
    """
    Class that returns the initial conditions as expected by the propagator. 
    The eval function returns the 4 quaternion components as well as the angular velocities.
    """
    
    def __init__(self, alpha, delta, omega, alphaDot, deltaDot, omegaDot):
        
        #we have to manually call the constructor
        common.Vector7TimeEvaluable.__init__(self)
        self.alpha = alpha
        self.delta = delta
        self.omega = omega
        self.alphaDot = alphaDot
        self.deltaDot = deltaDot
        self.omegaDot = omegaDot
        
    def eval(self, epoch):
        
        al = self.alpha.eval(epoch)
        de = self.delta.eval(epoch)
        om = self.omega.eval(epoch)
        alDot = self.alphaDot.eval(epoch)
        deDot = self.deltaDot.eval(epoch)
        omDot = self.omegaDot.eval(epoch)
        
        #this converts Euler angles to quaternion components
        quart = br.concatenate([ad.cos(al/2) * ad.cos(de/2) * ad.cos(om/2) + ad.sin(al/2) * ad.sin(de/2) * ad.sin(om/2),
                                ad.sin(al/2) * ad.cos(de/2) * ad.cos(om/2) - ad.cos(al/2) * ad.sin(de/2) * ad.sin(om/2),
                                ad.cos(al/2) * ad.sin(de/2) * ad.cos(om/2) + ad.sin(al/2) * ad.cos(de/2) * ad.sin(om/2),
                                ad.cos(al/2) * ad.cos(de/2) * ad.sin(om/2) - ad.sin(al/2) * ad.sin(de/2) * ad.cos(om/2)])
        
        #the angular velocity
        anglesDot = br.concatenate([alDot, deDot, omDot])
        
        return br.concatenate([quart, anglesDot])

q0 = Quaternion(alpha, delta, omega, alphaDot, deltaDot, omegaDot)

We set up the asteroid as an ellipsoid with corresponding moment of inertia. The axes of the ellipsoid have length 3, 4 and 5 km, respectively.

In [5]:

class Ellipsoid():
    
    def __init__(self, a, b, c, dens):
        self.a = a
        self.b = b
        self.c = c
        self.dens = dens.eval(t0)  
        self.vol = 4/3 * np.pi * self.a *self.b*self.c
        self.M = self.vol * self.dens
        
    def inertia(self):
        inertia = np.identity(3)
        inertia[0,0] = self.b**2 + self.c**2
        inertia[1,1] = self.a**2 + self.c**2
        inertia[2,2] = self.a**2 + self.b**2
        
        return self.M.value()/5 * inertia
    
    def coords(self, theta, phi):
        x = self.a * np.sin(theta) * np.cos(phi)
        y = self.b * np.sin(theta) * np.sin(phi)
        z = self.c * np.cos(theta)
        
        return x,y,z
    
ast = Ellipsoid(3, 4, 5, rho)

Propagating the asteroid rotation¶

Now we can propagate the rotation of the asteroid. At the moment there is no torque acting on it.

In [6]:

#set up a propagator axes
astPAxes = prop.PropagatorAxes(fra, astAxes, icrf)

#we specify the torque on the asteroid to be zero.
torque = common.ConstantVector3TimeEvaluable(np.zeros(3))

#the inertia of the asteroid is now equal to the inertia of the corresponding ellipsoid
inertia = ast.inertia()

#the input to the propagator takes the object and its initial conditions
asAxInput = prop.AxesInput(astPAxes, q0, torque, inertia, np.full(7, 1e-9), np.full(7, 1e-9))

inpAxes = [asAxInput]

#the initial epoch for the propagator
epoch0 = tempo.Epoch(0.0, tscale)

#the final epoch for the propagator, 2 days later
epoch1 = tempo.Epoch(7.0, tscale)

#xepochs are used if we want to track partial derivatives
xepoch1 = tempo.XEpoch(epoch1)

#set up the propagator object and run the computation
proAxes = prop.Propagator("axes", epoch0, None, inpAxes)
proAxes.compute(xepoch1)

We can demonstrate that the axes are rotating with a plot of the deviation of the z axis from its original location. The plot shows that the asteroid is precessing.

In [7]:

#create a list of epochs with step size 1000 seconds
grid = tempo.EpochRange(epoch0, epoch1).createGrid(100.)

#get the julian day pair of each epoch (float representation of epoch)
t = [e.mjd() for e in grid]

#the unit z vector
uz = np.array([0., 0., 1.])

#the unit z vector at all points of time
xvecs = np.array([fra.rotate(uz, icrf, astAxes, e) for e in grid])

#calculating the deviation of the z axis at all points of time
radii = []
angles = []

for xvec in xvecs:
    #angle between new and old z-vector: radius in plot
    delta = np.arccos(np.dot(uz, xvec))
    radii.append(delta)

    p = xvec - np.dot(xvec, uz) * uz
                    
    #angle in the xy-plane: angle in plot
    alpha = np.arctan2(p[1], p[0])
    angles.append(alpha)

ax = plt.subplot(111, projection = 'polar')
ax.plot(angles, radii)
plt.show()

Propagating the spacecraft¶

Now we can add the spacecraft to the system. We want the spacecraft to experience the gravitational acceleration of the asteroid. For this, we can set up a custom acceleration function by creating a class that inherits from an existing class. This creates a time evaluable obejct. We can then define the value of the acceleration at any time by overriding the eval method.

In [8]:

class Acceleration(common.Vector3TimeEvaluable):

    def __init__(self, frames, axes, p1, p2, M):
        
        # We have to explicitely call the inherited class constructor
        common.Vector3TimeEvaluable.__init__(self)

        # assign member variables
        self.fra = frames
        self.axes = axes
        self.p1 = p1
        self.p2 = p2
        self.M = M

    # the eval method returns the acceleration at the requested epoch
    def eval(self, epoch):
        
        # the method yields the vector between the spacecraft and the axes
        r = self.fra.vector3(self.p2, self.p1, self.axes, epoch)
        d = br.norm(r)
        acc = -self.M.value() *r / br.pow(d,3)
        return acc

For the initial conditions of the spacecraft we also create a class that inherits from common and returns the position and velocity of the spacecraft.

In [9]:

class SpacecraftState(common.Vector6TimeEvaluable):
    
    def __init__(self, posx, posy, posz, velx, vely, velz):
        
        common.Vector6TimeEvaluable.__init__(self)
        
        self.posx = posx
        self.posy = posy
        self.posz = posz
        self.velx = velx
        self.vely = vely
        self.velz = velz
        
    def eval(self, epoch):
        
        px = self.posx.eval(epoch)
        py = self.posy.eval(epoch)
        pz = self.posz.eval(epoch)
        vx = self.velx.eval(epoch)
        vy = self.vely.eval(epoch)
        vz = self.velz.eval(epoch)
        
        return br.concatenate([px, py, pz, vx, vy, vz])

Having set up the acceleration of the spacecraft we proceed exactly as before with the propagation. If we wanted to resolve the mutual gravitational attraction of the asteroid and the spacecraft, we would use the same propagator for both of them.

In [10]:

scPPoint = prop.PropagatorPoint(fra, icrf, scPoint)

#we add the asteroid mass to the parameter book as it is required for the acceleration of the spacecraft
G = 6.67408e-20
astMu = ast.M * G
scAcc = Acceleration(fra, icrf, scPoint, astPoint, astMu)

state0 = SpacecraftState(posx, posy, posz, velx, vely, velz)

scInput = prop.PointInput(scPPoint, astPoint, state0, scAcc, np.full(6, 1e-9), np.full(6, 1e-9))
inpPoints = [scInput]

scProp = prop.Propagator("points", epoch0, None, inpPoints)
scProp.compute(xepoch1)

Below, the orbit of the spacecraft is plotted.

In [11]:

scData = np.asarray([fra.vector3(astPoint, scPoint, icrf, e) for e in grid])

plt.figure(figsize = (8,8))
plt.scatter(0, 0, s = 10.)
plt.plot(scData[:, 0], scData[:, 1])
plt.xlabel('x')
plt.ylabel('y')
plt.show()

Adding landmarks to the asteroid's surface¶

Finally, we can add landmarks to the surface of the asteroid. They are randomly distributed at time $t = t_0$ and the propagation of the asteroid axes allows us to track their position.

The landmarks are added to the frames system as a fixed point so that they will always retain their orientation and position towards the asteroid.

The normals are also added as fixed points relative to the respective point.

In [12]:

#np.random.seed(12)

xyz = ['x', 'y', 'z']

#the landmarks
points = []

#the surface normals
normals = []

num_lm = 10

for i in range(num_lm):
    #first we generate random angles and map them to the surface of the asteroid
    theta = np.random.rand(1) * np.pi
    phi = np.random.rand(1) * 2 * np.pi
    coords = ast.coords(theta, phi)    
    
    #here the positions of the landmarks are added to the parameter book
    for j in range(3):
        uni.parameters.add("lm_{}{}".format(i, xyz[j]), coords[j])
    lmVec = uni.parameters.concat3(["lm_{}{}".format(i, c) for c in xyz])
    
    lmPoint = uni.frames.addPoint("lm_{}".format(i), tscale)
    normalPoint = uni.frames.addPoint("normal_{}".format(i), tscale)
    
    #the point to the frames system is fixed relative to the asteroid
    uni.frames.addFixedTranslation(lmPoint, astPoint, astAxes, lmVec.eval(t0).value())
    
    #the surface normal
    normal = np.array([2 * coords[0] / ast.a**2, 2 * coords[1] / ast.b**2, 2 * coords[2] / ast.c**2])
    
    #it is added as fixed to the corresponding point
    uni.frames.addFixedTranslation(normalPoint, astPoint, astAxes, normal)
    
    normals.append(normalPoint)
    points.append(lmPoint)

Landmarks visible to snapshots by spacecraft¶

We assume that the spacecraft takes 5 snapshots of the asteroid. At each snapshot we check which of the landmarks are visible. The condition for this is that the angle spanned by the vector connecting the spacecraft and a landmark and the surface normal vector is less than 90 degrees. If they are visible we compute the vector between the spacecraft and the asteroid.

In [13]:

num_snaps = 3

num_secs = epoch1 - epoch0

snaps = [epoch0 + i * num_secs/ (num_snaps -1)  for i in range(num_snaps)]

#a list of all vectors between the spacecraft and visible points
visVecs = []
for i, e in enumerate(snaps):
    vec = []
    print("Points visible at time {}:".format(e))
    for j, p in enumerate(points):
        
        #vector between the spacecraft and the landmark
        sc_p = uni.frames.vector3(scPoint, p, icrf, tempo.XEpoch(e))
        
        #the resolve method resolves the partial derivatives to the initial parameters
        sc_p = common.resolve(sc_p)
        
        #the normal vector
        normal = uni.frames.vector3(p, normals[j], icrf, tempo.XEpoch(e))
            
        #if the angle between the vectors is greater than 90 deg it is visible
        if br.dot(sc_p, normal).value() > 0.:
            var = sc_p / br.norm(sc_p)
            vec.append((j, var, br.norm(sc_p)))
            print(j)

    visVecs.append(vec)
    print("\n")

Points visible at time 2000-01-01T00:00:00.000000 TDB:
0
1
4
6
8


Points visible at time 2000-01-04T12:00:00.000000 TDB:
0
3
4
7
8


Points visible at time 2000-01-08T00:00:00.000000 TDB:
1
2
5
6
9

The unit vectors also carry the partial derivatives w.r.t. all the parameters which we "tracked" using pb.track() at the beginning of the code.

In [14]:

import json

jsonObj = []

for i, snap in enumerate(visVecs):
    if len(snap) == 0:
        continue
    lmList = []
    for lm in snap:
        u = lm[1].value()
        ra = np.arctan2(u[1], u[0])
        de = np.arctan2(u[2], np.linalg.norm(u[:2]))
        d = {"target": "lm_{}".format(lm[0]), "right ascension": ra, "declination": de, "noise": 1e-4}
        lmList.append(d)
    jsonObj.append({"epoch": str(snaps[i]), "camera": "spacecraft", "axes": "ICRF", "vectors": lmList })

with open("../cosmoslandmarks.json", "w") as f:
    json.dump(jsonObj, f, indent=4)

The plots below show the position of the spacecraft and the visible landmarks at each snapshot projected onto the xy-plane.

In [15]:

scPos = np.empty((num_snaps,3))

scData = np.asarray([fra.vector3(astPoint, scPoint, icrf, e) for e in grid])

for i,s in enumerate(snaps):
    plt.figure(figsize = (8,8))
    plt.plot(scData[:, 0], scData[:, 1])

    scPos[i] = fra.vector3(astPoint, scPoint, icrf, s)
    
    for p in visVecs[i]:
        ast_p = p[1] * p[2] + scPos[i]
        plt.scatter(ast_p.value()[0], ast_p.value()[1], s = 10.)
        
    plt.scatter(scPos[i,0], scPos[i,1], label = "snapshot", color = 'blue')
    plt.scatter(0,0, label = 'asteroid', s = 15., color = 'black')
    plt.legend()
    plt.show()

In [16]:

## BUG, consider cannot be null vector. 
# from godot.cosmos import orb
# config =    {
#                "parameters" : 
#                {
#                    "free" : ["posx", "posy", "posz", "velx", "vely", "velz"],
#                    "consider" : []
#                },
#                "equations":
#                [
#                    {
#                        "type" : "direction" , 
#                        "file" : "../cosmoslandmarks.json"
#                    }
#                ]
#             }
#
#
#np.set_printoptions(linewidth=150)
#prob = orb.Problem(uni, config)
#sol = orb.Solver(prob)
#sol.process()
#sol.solve()
#print(sol.filter().covariance())
#
#solution = sol.filter().solution()
#cov = sol.filter().covariance()
#
#for i,c in enumerate(cov):
#    print(solution[i], "+/-", np.sqrt(c[i]))
#sol.update()