This example notebook was originally developed by Lynne Jones (@rhiannonlynne) and subsequently updated.
This notebook is currently under revision and not working properly. Please refer to the examples provided in the documentation for now.
In [1]:
import os
import numpy as np
import pyoorb as oo
In [2]:
# Initialize oorb
oo.pyoorb.oorb_init()
Out[2]:
In [3]:
timeScales = {'UTC': 1, 'UT1': 2, 'TT': 3, 'TAI': 4}
elemType = {'CART': 1, 'COM': 2, 'KEP': 3, 'DEL': 4, 'EQX': 5}
In [4]:
# Set up some orbits
# orb is id, 6 elements, epoch_mjd, H, G, element type index
# keplerian appears to be element type index 3
# orbits = numpy.array([0.,1.,2.,3.,4.,5.,6.,5373.,1.,1.,3.])
o = np.loadtxt('test.des', dtype=([('objid', np.str_, 20), ('format', np.str_, 3),
('q', float), ('e', float), ('inc', float), ('Omega', float),
('argperi', float), ('tperi', float), ('H', float), ('epoch', float),
('index', int), ('npar', int), ('moid', float), ('compcode', np.str_, 10)]),
skiprows=1)
orbits = np.zeros([len(o), 12], dtype=np.double, order='F')
orbits[:, 0] = np.arange(0, len(o), 1, dtype=float)
orbits[:, 1] = o['q']
orbits[:, 2] = o['e']
orbits[:, 3] = o['inc']
orbits[:, 4] = o['Omega']
orbits[:, 5] = o['argperi']
orbits[:, 6] = o['tperi']
orbits[:, 7] = np.zeros(len(o), float) + elemType['COM']
orbits[:, 8] = o['epoch']
orbits[:, 9] = np.zeros(len(o), float) + timeScales['TT']
orbits[:, 10] = o['H']
orbits[:, 11] = np.zeros(len(o)) + 0.15
for i in range(3, 6):
orbits[:, i] = np.radians(orbits[:, i])
In [5]:
orbits
Out[5]:
In [6]:
offset = 30.
times = [orbits[0][8], orbits[0][8] + offset]
timescale = [timeScales['UTC']] * len(times)
ephem_dates = np.array(list(zip(times, timescale)), dtype=np.double, order='F')
ephs, err = oo.pyoorb.oorb_ephemeris_full(in_orbits=orbits,
in_obscode='I11',
in_date_ephems=ephem_dates,
in_dynmodel='N')
if err != 0:
print(err)
In [7]:
i = 0
start = ephs[i][0]
end = ephs[i][1]
print("@ %f" % (ephem_dates[0][0]), start)
print("@ %f" % (ephem_dates[1][0]), end)
In [8]:
type2 = 'CART'
print("Transforming from COM to %s" % (type2))
newElems, err = oo.pyoorb.oorb_element_transformation(in_orbits=orbits, in_element_type=elemType[type2])
if err != 0:
print("error", err)
print("Input and output orbital elements (COM / CART)")
print("in ", orbits[i])
print("out", newElems[i])
In [9]:
print("Transforming back from %s to COM" % (type2))
newElems2, err = oo.pyoorb.oorb_element_transformation(in_orbits=newElems, in_element_type=elemType['COM'])
if err != 0:
print("error", err)
print("Input and output orbital elements (CART / COM)")
print("in ", newElems[i])
print("out", newElems2[i])
In [13]:
print('Differences in the initial/final COM orb elements')
diffs = newElems2 - orbits
print(diffs)
print(newElems2[1,3])
print(orbits[1,3])
print(diffs[1,3])
There can be larger differences in the tPeri, as this is degenerate (previous or next orbit?).
In [14]:
print("Propagating orbits by %d days" % (offset))
epoch_orig = orbits[i][8]
epoch_new = epoch_orig + offset
newEpoch = np.array([epoch_new, timeScales['TT']], dtype='double', order='F')
# Note that currently you have to use CART format elements.
elems = newElems
newOorbElems, err = oo.pyoorb.oorb_propagation(in_orbits=elems, in_epoch=newEpoch, in_dynmodel='N')
if err != 0:
print("error", err)
print("Input and output orbital elements (0 and %d days)" % (offset))
print("in ", orbits[i])
print("out", newOorbElems[i])
In [15]:
print("Transforming new orbital elements to COM")
newElems3, err = oo.pyoorb.oorb_element_transformation(in_orbits=newOorbElems, in_element_type=elemType['COM'])
if err != 0:
print("error", err)
print("Input and output orbital elements (CART/COM)")
print("in ", orbits[i])
print("out", newElems3[i])
print('epoch', epoch_new, newElems3[i][8])
In [17]:
# Generate ephemerides with these new elements, to compare.
ephs2, err = oo.pyoorb.oorb_ephemeris_full(in_orbits=newElems3, in_obscode='I11',
in_date_ephems=ephem_dates, in_dynmodel='N')
for i in range(len(newElems3)):
print('Object %d' % i)
print(' Orbit: ', orbits[i] )
diffs = ephs[i] - ephs2[i]
#print("Difference in ephemerides from ephemeris generation vs propagation + ephemeris generation")
#print(diffs)
print("Just difference in mas in RA and then Dec")
deltaRA = diffs[:,1] * 3600. * 1000.
deltaDec = diffs[:,2] * 3600. * 1000.
print(' ', deltaRA, deltaDec)
Note that the difference can change depending on the length of time in 'offset'. Needs a bit further investigation, but I suspect this is due to the timescale of integration hard-coded into the different routines (?).
In [ ]: