Apophis: Applied Orbit Determination with godot¶

At first we import all the necessary modules.

In [1]:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse

In [2]:

from godot.core import tempo
from godot.model import common
from godot.core.autodif import bridge
from godot import cosmos
from godot.cosmos import orb
from godot.core import events, astro

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

Propagating the initial (perturbed) orbit¶

We are using the high level universe and trajectory modules of godot. All the information about the gravity model and the initial conditions of the trajectory are contained in the yml files used to set up the universe and the trajectory. For now, we use the perturbed initial state in "collisionTrajectory.yml".

In [3]:

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

#set up trajectory with perturbed data
traConfig = cosmos.util.load_yaml('collisionTrajectory.yml')
tra = cosmos.Trajectory(uni, traConfig)

In [4]:

#track partial derivatives and compute trajectory
uni.parameters.track(uni.parameters.contains())
tra.compute(partials = True)

We can test if Apophis will collide with earth by settig up a function which computes the distance between Apophis and earth minus the earth radius. When this function crosses 0 we have a collision. Godot allows us to compute the function everywhere using the events interface (it performs a root search).

In [5]:

earthRadius = 6357.

def collision(epoch):
    x = uni.frames.vector('Earth', 'ApophisPerturbed_center', 'ICRF', 0, epoch)
    return bridge.norm(x) - earthRadius

# Generate event interval set for observation periods
t0 = tempo.parseEpoch("2020-01-02T00:00:00 TDB")
t1 = tempo.parseEpoch("2030-01-01T00:00:00 TDB")
grid = tempo.EpochRange(t0, t1).createGrid(86400.)

collisions = events.generateEventSet(collision, 1.,  grid, 1.)

print(collisions)

2029-04-13T21:16:27.803839 TDB   2.046726750e-04 +++ {}
2029-04-13T21:31:08.122988 TDB  -2.530589783e-02 +++ {}

We see that, with the the current initial conditions, we expect a collision in April 2029. We can analyze this encounter further by mapping it on the B-plane. Using the partial derivatives and a guess for the initial covariance of the state we can also compute the covariance of the B-plane coordinates.

In [6]:

earthMu = uni.constants.getMu("Earth")

#select an epoch just before the encounter to map onto the B-plane
preEncounter = tempo.XEpoch(tempo.parseEpoch("2029-04-13T21:00:00 TDB"))

#get the position vector between earth and apophis
cartesian = uni.frames.vector6("Earth", "ApophisPerturbed_center", "ICRF", preEncounter)

#get the B-plane coordinates and resolve them
bplane = astro.vinfBplaneFromCart(cartesian, earthMu)
bplane = common.resolve(bplane)

#the B-plane coordinates
br = bplane[3]
bt = bplane[4]

#we guess an initial, large covariance
cov0 = np.identity(6)
cov0[0][0] = cov0[1][1] = cov0[2][2] = 1e2**2
cov0[3][3] = cov0[4][4] = cov0[5][5] = 0.0001**2

#get the partial derivatives with respect to the initial conditions
partials = np.append(
    [br.at(["start_ApophisPerturbed_center_pos_x", "start_ApophisPerturbed_center_pos_y", "start_ApophisPerturbed_center_pos_z",
    "start_ApophisPerturbed_center_vel_x", "start_ApophisPerturbed_center_vel_y", "start_ApophisPerturbed_center_vel_z"])],
    [bt.at(["start_ApophisPerturbed_center_pos_x", "start_ApophisPerturbed_center_pos_y", "start_ApophisPerturbed_center_pos_z",
    "start_ApophisPerturbed_center_vel_x", "start_ApophisPerturbed_center_vel_y", "start_ApophisPerturbed_center_vel_z"])], axis = 0)

#here we map the covariance of the initial conditions onto the B-plane
covBplane = np.matmul(partials, np.matmul(cov0, partials.T))

eigenval, eigenvec = np.linalg.eig(covBplane)

#the axes and rotation of the error ellipse
spread = np.sqrt(eigenval)
rot = np.rad2deg(np.arccos(eigenvec[0][0]))

#calculation of effective erth radius
vinf = np.linalg.norm(bplane[:3].value())
effEarthRadius = earthRadius * np.sqrt(1 + 2 * earthMu / (earthRadius * np.dot(vinf,vinf)))


#plotting of data
angles = np.linspace(0, 2 * np.pi, 10000)
fig = plt.figure(figsize = (8,8))
ax = fig.add_subplot()
ell = Ellipse( xy = (br.value(), bt.value()), width=spread[0]*2, height=spread[1]*2, 
              angle = rot, edgecolor='black', lw=1., facecolor='none')
ax.add_artist(ell)
ax.plot(earthRadius * np.cos(angles), earthRadius * np.sin(angles), label = "earth radius")
ax.plot(effEarthRadius * np.cos(angles), effEarthRadius * np.sin(angles), "--", label = "effective earth radius")
ax.scatter(br.value(), bt.value(), color = "blue")
plt.xlim(-1e5, 1e5)
plt.ylim(-1e5, 1e5)
plt.xlabel("km")
plt.ylabel("km")
plt.legend()

plt.show()

We see that initially the asteroid appears to be heading for earth with a large uncertainty. The uncertainty ellipse is stretched out along a line. We can try and reduce this uncertainty by generating observations.

Creating observations with catalogue data¶

So far we have been using the perturbed initial conditions. We can now use the catalogue value to create some initial conditions. We repeat the propagation steps with the catalogue data in "trajectory.yml".

In [7]:

#same as for perturbed trajectory earlier
traConfig2 = cosmos.util.load_yaml('trajectory.yml')
tra2 = cosmos.Trajectory(uni, traConfig2)
tra2.compute(partials = False)

No we create the observation file with the catalogue data. We compute an observation a day when the asteroid is visible to the ground station and the sun and moon are not.

In [8]:

# Generate event interval set for observation periods
t0 = tempo.parseEpoch("2017-01-01T00:00:00 TDB")
t1 = tempo.parseEpoch("2025-01-01T00:00:00 TDB")

grid = tempo.EpochRange(t0, t1).createGrid(100000.)


#get the elevations to compute visibility
def moonElev( epo ):
    x = uni.frames.vector3('New_Norcia', 'Moon', 'New_Norcia', epo)
    return astro.sphericalFromCart(x)[2]

def sunElev( epo ):
    x = uni.frames.vector3('New_Norcia', 'Sun', 'New_Norcia', epo)
    return astro.sphericalFromCart(x)[2]

def aphElev( epo ):
    x = uni.frames.vector3('New_Norcia', 'Apophis_center', 'New_Norcia', epo)
    return astro.sphericalFromCart(x)[2]

#computation of event intervals
moon_visible = events.generateEventIntervalSet( moonElev, 1e-7, grid, 1e-9 )
sun_visible = events.generateEventIntervalSet( sunElev, 1e-7, grid, 1e-9 )
aph_visible = events.generateEventIntervalSet( aphElev, 1e-7, grid, 1e-9 )
sun_or_moon_visible = events.merge( [ sun_visible , moon_visible ])
no_sun_or_moon_visible = sun_or_moon_visible.complement()
aph_visible_without_sun_and_moon = events.overlap( [ no_sun_or_moon_visible , aph_visible ])

/local/scratch/gvarga/godottut/env/lib64/python3.6/site-packages/ipykernel_launcher.py:25: DeprecationWarning: The merge(list()) function is deprecated and it will be soon removed, please use the "any" function instead
/local/scratch/gvarga/godottut/env/lib64/python3.6/site-packages/ipykernel_launcher.py:27: DeprecationWarning: The overlap(list()) function is deprecated and it will be soon removed, please use the "all" function instead

For each day where apophis is visible and the sun and the moon are not, we now measure the right ascension and declination to the asteroid and add some Gaussian noise to it. All the observations are stored in the json format.

In [9]:

import json
jsonObj = []
for ei in aph_visible_without_sun_and_moon:
    for time in tempo.EpochRange(ei.start().value(),ei.end().value() ).createGrid(864000):
        obs = uni.frames.vector('New_Norcia', 'Apophis_center', 'ICRF', 0, time)
        obs /= np.linalg.norm(obs)
        ra = np.arctan2(obs[1], obs[0]) + np.random.randn(1)[0] * 1e-6
        de = np.arctan2(obs[2], np.linalg.norm(obs[:2])) + np.random.randn(1)[0] * 1e-6
        d = {"target": "start_ApophisPerturbed_center", "right ascension": ra, "declination": de, "noise": 1e-4}
        jsonObj.append({"epoch": str(time), "camera": "New_Norcia", "axes": "ICRF", "vectors": [d] })
with open("cosmosapophisObservations.json", "w") as f:
    json.dump(jsonObj, f, indent=4)

Fitting the observations to the perturbed model¶

We now have a perturbed model and some observations of the actual data. Therefore we can perform orbit determination and correct the perturbed orbit so that the residuals between observations and model are minimized. We do this with the module godot.problem and godot.Solver. At first, we set up the configuration file with the information on the observables and the free and consider parameters.

In [10]:

config =    {
                "parameters" : 
                {
                    "free" : ["start_ApophisPerturbed_center_pos_x", "start_ApophisPerturbed_center_pos_y", "start_ApophisPerturbed_center_pos_z",
                              "start_ApophisPerturbed_center_vel_x", "start_ApophisPerturbed_center_vel_y", "start_ApophisPerturbed_center_vel_z"],
#                    "consider" : []
                },
                "equations":
                [
                    {
                        "type" : "direction",
                        "name" : "ground observations",
                        "config":
                        {
                            "file" : "./cosmosapophisObservations.json"
                        }
                    }
                ]
            }


prob = orb.Problem(uni, config)
sol = orb.Solver(prob)

#calculate the residuals
pre_residuals = sol.processProblemEquations()
sol.solve()

#update the model and recompute
sol.update()
tra.compute(partials = True)

#calculate residuals with new model
prob = orb.Problem(uni, config)
sol = orb.Solver(prob)
post_residuals = np.array(sol.processProblemEquations())
sol.solve()
sol.update()

#get the filter covariance
con = sol.considerSolution()

# get the bplane coordinates
cartesian = uni.frames.vector6("Earth", "ApophisPerturbed_center", "ICRF", preEncounter)
bplane = astro.vinfBplaneFromCart(cartesian, earthMu)
bplane = common.resolve(bplane)
#the B-plane coordinates
br1 = bplane[3]
bt1 = bplane[4]


#compute the covariance of the bplane coordinates
bplaneCov = con.covariance(bplane)

#calculate the new error ellipse
eigenval, eigenvec = np.linalg.eig(bplaneCov)
spread1 = np.sqrt(eigenval)

rot1 = np.rad2deg(np.arccos(eigenvec[0][0]))

/local/scratch/gvarga/godottut/env/lib64/python3.6/site-packages/ipykernel_launcher.py:57: RuntimeWarning: invalid value encountered in sqrt

In [ ]:

We can estimate the quality of the fit by comparing the pre-fit and post-fit residuals.

In [11]:

# ####
# # WARNING RESIDUALS WILL NOT WORK FOR NOW - SKIP TO NEXT
# ####

# data = cosmos.util.load_json("cosmosapophisObservations.json")

# pre_residuals = np.vstack(pre_residuals)
# post_residuals = np.vstack(post_residuals)

# #compute the dates for plotting
# dates = [str(d["epoch"]).split("T")[0] for d in data[:len(pre_residuals)]]
# dates = np.array(dates, dtype=np.datetime64)

# #plot the residuals
# fig, ax  = plt.subplots(nrows = 3, ncols = 2, sharex = True, figsize = (12,10))

# ax[0,0].scatter(dates , pre_residuals[:,0], color = 'C1', s = 1.)
# ax[1,0].scatter(dates , pre_residuals[:,1], color = 'C1', s = 1.)

# ax[0,1].scatter(dates , post_residuals[:,0], color = 'C1', s = 1.)
# ax[1,1].scatter(dates , post_residuals[:,1], color = 'C1', s = 1.)


# t0 = tempo.parseEpoch(data[0]["epoch"])
# t1 = tempo.parseEpoch(data[-1]["epoch"])
# grid = tempo.EpochRange(t0, t1).createGrid(86400)

# from godot.core.autodif import bridge
# dists = [bridge.norm(uni.frames.vector("Earth", "Apophis_center", "ICRF", 0, e)) for e in grid]
# plotGrid = np.array([str(d).split("T")[0] for d in grid], dtype = np.datetime64)
# ax[2,0].plot(plotGrid, np.array(dists) / 149598073 , color = 'C1')
# ax[2,1].plot(plotGrid, np.array(dists) / 149598073 , color = 'C1')

# for a in ax[0]:
#     a.set_title("Right Ascension", fontsize = 14)
#     a.set_ylabel("Residual", fontsize = 14)

    
# for a in ax[1]:
#     a.set_title("Declination", fontsize = 14)
#     a.set_ylabel("Residual", fontsize = 14)

# for a in ax[2]:
#     a.set_title("Distance to Earth", fontsize = 14)
#     a.set_ylabel(r"Distance [AU]",  fontsize = 14)
#     a.set_xlabel("Epoch", fontsize = 14)
    
# for row in ax:
#     for a in row:
#         a.grid()


# plt.tight_layout()
# plt.show()

We see that the fit has reduced the initial resiudals to the Gaussian noise we added to the observations. Now we can see how the fit looks like in the B plane:

In [12]:

# plot solution, again

angles = np.linspace(0, 2 * np.pi, 10000)
fig = plt.figure(figsize = (8,8))
ax = fig.add_subplot()

ell = Ellipse( xy = (br.value(), bt.value()), width=spread[0]*2, height=spread[1]*2, 
              angle = rot, edgecolor='b', lw=1., facecolor='none', alpha = 0.5)
ax.add_artist(ell)

ell1 = Ellipse( xy = (br1.value(), bt1.value()), width=spread1[0]*2, height=spread1[1]*2, 
              angle = rot1, edgecolor='black', lw=1., facecolor='none')
ax.add_artist(ell1)

ell2 = Ellipse( xy = (br1.value(), bt1.value()), width=spread1[0]*2 * 2, height=spread1[1]*2 * 2, 
              angle = rot1, edgecolor='black', lw=1., facecolor='none')
ax.add_artist(ell2)

ell3 = Ellipse( xy = (br1.value(), bt1.value()), width=spread1[0]*2 * 3, height=spread1[1] * 2 * 3, 
              angle = rot1, edgecolor='black', lw=1., facecolor='none')
ax.add_artist(ell3)


ax.plot(earthRadius * np.cos(angles), earthRadius * np.sin(angles))
ax.plot(effEarthRadius * np.cos(angles), effEarthRadius * np.sin(angles), "--")
ax.scatter(br.value(), bt.value(), label = "before observations")
ax.scatter(br1.value(), bt1.value(), color = "black", label = "after observations")
plt.legend()
plt.xlim(-1e5, 1e5)
plt.ylim(-1e5, 1e5)

plt.xlabel("km")
plt.ylabel("km")

plt.show()

The solution has converged nicely onto the catalogue value. The error ellipses represent 1, 2 and 3 standard deviations, respectively.