Source code for irfpy.pep.galileo_orbit

''' Galileo orbital information based on PDS.

The data is downloaded from
http://ppi.pds.nasa.gov/search/view/?f=yes&id=pds://PPI/GOMW_5001/DATA/TRAJECT

.. warning::

    Note that this file contains a browse data, thus, the orbit information
    may have worse time resolution especially during the flyby.

    In such a case you should use spice data with :mod:`galileo_spice` module.

Then the downloaded files are manually merged and processed.

.. todo::

    Data file in the data repositry.


Raw data contains 12 columns.

..

    Column 1  | UT
    Column 2  | R in Rj (71492 km)
    Column 3  | Latitude in System III
    Column 4  | Longitude in System III
    Column 5  | JSE_X
    Column 6  | JSE_Y
    Column 7  | JSE_Z
    Column 8  | JSM_X
    Column 9  | JSM_Y
    Column 10 | JSM_Z
    Column 11 | Local hour
    Column 12 | Local magnetic hour

We only need the column 1, 5-7.

.. code-block:: sh

    % cat /Users/futaana/Downloads/GOMW_5001/SURVEY/*.TAB > galileo_orbit_raw.txt
    % awk '{print $1, $5, $6, $7}' < galileo_orbit_raw.txt > galileo_orbit_sample.txt
    % gzip galileo_orbit_sample.txt
    % vi .irfpyrc
    [galileo]
    galileo_orbit_browse = <path_to:galileo_orbit_sample.txt.gz>

'''
import os
import gzip
import io
import logging
logging.basicConfig()
from pkg_resources import resource_filename
import dateutil.parser

import numpy as np
import scipy as sp
from scipy import interpolate
import matplotlib.dates

from irfpy.util import irfpyrc
from irfpy.util import exception


[docs]class GalileoOrbit:
    ''' Galileo orbit from discretely sampled points

    You may get the default instance by

    .. code-block:: python

        go = GalileoOrbit.get_default_instance()

    Reading and parsing takes time.

    >>> go = GalileoOrbit.get_default_instance()
    >>> print(go.t0(), go.t1())    # Data coverage   # doctest: +SKIP
    1995-10-09 00:00:00 1997-05-03 23:50:00

    >>> import datetime
    >>> print(go.get_position(datetime.datetime(1995, 10, 25, 0, 0, 0)))
    [   8.769 -394.688  -36.523]
    '''

    __default_instance = None
    __reduced_instance = None


[docs]    @classmethod
    def get_default_instance(cls):
        if cls.__default_instance == None:
            cls.__default_instance = GalileoOrbit()
        return cls.__default_instance



[docs]    @classmethod
    def get_reduced_instance(cls):
        import warnings
        warnings.warn('Use of get_reduced_instance is not supported now.')
        warnings.warn('Just return the full instance.')
        
        return cls.get_default_instance()


    def __init__(self, filename=None):
        self.logger = logging.getLogger(self.__class__.__name__)
        self.logger.setLevel(logging.INFO)
        if filename:
            fn = filename
        else:
            rc = irfpyrc.Rc()
            fn = rc.get('galileo', 'galileo_orbit_browse')
            if fn == None:
                raise exception.PyanaException('No entry found: [galileo]-galileo_orbit_browse')
        self.logger.info('File name: %s' % fn)

        if not os.path.exists(fn):
            raise IOError('File %s not found.' % fn)
        
#        self.dat = np.loadtxt(fn, converters={0: dateutil.parser.parse})
        self.tlist = []  # Stores datetime instnce
        self.xlist = []
        self.ylist = []
        self.zlist = []

### Works only for Python3.3
#        f = gzip.open(fn, 'rt')

### Better for Python3.
        f = io.TextIOWrapper(gzip.open(fn))
        for l in f:
            tok = l.split()
            t = dateutil.parser.parse(tok[0])
#            p = np.array([float(tok[1]), float(tok[2]), float(tok[3])])
#            self.posdict[t] = p
            x = float(tok[1])
            y = float(tok[2])
            z = float(tok[3])

            self.tlist.append(t)
            self.xlist.append(x)
            self.ylist.append(y)
            self.zlist.append(z)
            
        f.close()

        self.tnlist = matplotlib.dates.date2num(self.tlist)   # Stores float instance.

        # Interpolation.
        self.xfunc = interpolate.InterpolatedUnivariateSpline(self.tnlist, self.xlist)
        self.yfunc = interpolate.InterpolatedUnivariateSpline(self.tnlist, self.ylist)
        self.zfunc = interpolate.InterpolatedUnivariateSpline(self.tnlist, self.zlist)


[docs]    def t0(self):
        return min(self.tlist)



[docs]    def t1(self):
        return max(self.tlist)



[docs]    def get_position(self, t):
        ''' Get the position of Galileo spacecraft at the time of *t*.

        Interpolation by each component (x, y, z) is used to calculate.
        The position is in the unit of Rj.
        '''
        tn = matplotlib.dates.date2num(t)
        if tn < self.tnlist.min() or tn > self.tnlist.max():
            raise ValueError('Given time %s is not in the range of data [%s %s]' %
                            (t, min(self.tlist), max(self.tlist)))
        x = self.xfunc(tn)
        y = self.yfunc(tn)
        z = self.zfunc(tn)

        return np.array([x, y, z])

        


[docs]    def get_positions(self, t):
        ''' Same as get_position but t can be an array.
        '''
        if np.iterable(t):
            return np.array([self.get_position(ti) for ti in t])
        else:
            return self.get_position(t)







import unittest
import doctest



[docs]def doctests():
    return unittest.TestSuite((
        doctest.DocTestSuite(),
        ))

if __name__ == '__main__':
    unittest.main(defaultTest='doctests')