''' Trajectory related class.
'''
import numpy as np
import logging
logging.basicConfig()
from pyana.util import keplernumeric as kepler
from . import observer
from . import parameters
class Trajectory:
''' Trajectory class.
First, :class:`observer.VirtualObserver` object should be generated.
>>> sc = observer.VirtualSpacecraft([1e6, 0, 0], [0., 10, 0])
>>> snsr = observer.VirtualSensor([-1., 0, -.1], 100.)
>>> obsr = observer.VirtualObserver(sc, snsr)
``sc`` is the spacecraft object located at 1000000 km in ``x``,
moving with 10km/s in y direction.
``snsr`` is the sensor object looking toward -x with velocity of 100 km/s
(namely observing vx~-100 km/s particle, if ram velocity is neglected),
and ``obsr`` is the observer object with above measurement.
>>> traj = Trajectory(obsr)
will create the trajectory instance, and :meth:`calc_coordinates` will
return the trajectory coordinates at a given time.
>>> tlist = np.array([0, -1000, -2000, -3000, -4000, -5000, -6000, -7000, -8000, -9000])
>>> coord = traj.get_coordinates(tlist)
>>> print(coord.shape)
(6, 10)
>>> length = traj.get_length(tlist)[0]
>>> for i in range(len(length)): print('%.2f' % length[i])
0.00
100565.93
201287.13
302204.13
403375.49
504889.91
606890.24
709624.25
813556.04
919520.19
'''
logger = logging.getLogger('Trajectory')
# logger.setLevel(logging.DEBUG)
def __init__(self, observer):
self.observer = observer
def get_coordinates(self, tlist, mj=parameters.mj):
''' Return 6-D coordinates of the trajectory of particles.
According to the given time list,
particle trajectories are calculated.
:param tlist: List of time. Unit in seconds.
The first element (``tlist[0]``) will be the origin.
:type tlist: ``np.array``
:returns: The trajectory coordinates. (6, N) shaped ``np.array``.
The first dimension (6) is x, y, z, vx, vy, vz.
Units are km or km/s.
The second dimension, ``N``, is time, and ``N=len(tlist)``.
'''
logger = self.logger # Alias.
### S/C pos, i.e. ENA pos at t=t0, MKSA.
enapos_t0 = self.observer.spacecraft.pos * 1e3
### ENA vel at t=t0, considering RAM velocity, MKSA.
enavel_t0 = self.observer.get_observe_velocity() * 1e3
logger.debug(enapos_t0)
logger.debug(enavel_t0)
### Kepler motion.
kep = kepler.KeplerRungeKutta(mj, 1., enapos_t0, enavel_t0)
posvel = kep.get_posvel(tlist).T # (6, N) after transpose.
logger.debug(posvel.shape)
return posvel / 1e3 # in km or km/s
def get_length(self, tlist, mj=parameters.mj):
''' Return the accumulate length in km along the trajectory.
:returns: (l, coord) where l is the distsnce list (N,) shape
and coord is the same output as :meth:`get_coordinates`.
'''
posvel = self.get_coordinates(tlist, mj=mj)
coord = posvel[:3, :]
l = np.zeros([len(coord[0, :])])
l[1:] = np.sqrt(((coord[:, 1:] - coord[:, :-1]) ** 2).sum(0))
# for i in range(1, len(l)):
# l[i] = l[i - 1] + l[i]
l = np.cumsum(l)
return l, posvel
import unittest
import doctest
def doctests():
return unittest.TestSuite((
doctest.DocTestSuite(),
))
if __name__ == '__main__':
unittest.main(defaultTest='doctests')