Source code for irfpy.vexpvat.pvat2

''' More sophisticated pvat retriever.

.. warning::

    This module is deprecated. Use :mod:`irfpy.vexpvat.vexspice` module.

irfpy.vexpvat.pvat2 provides more sophisticated functinality
of :mod:`irfpy.vexpvat.pvat`. Compensations are the slow processing.

As the datafile in Peje's format is only 1 minutes resolution,
this may not enough to provide more detailed analysis.
This module provides the interplation.

For simple implementation, the module provides the orbit number
based functionality.
There is small data gap in the edges of one file,
but it should not be significant
because the data gap is very small and they are normally
in the apocenter where no spacecraft motion is expected much.

**Usage**

First, instance :class:`PvatFile`. Usually, one can use :func:`getPvatFile`
for this purpose.

>>> if not isdb():
...     pass    # doctest: +SKIP
...     pvat = getPvatFile(1500)   # Instance the PVAT file for 1500.
...     t0 = datetime.datetime(2010, 5, 29, 20, 15, 23)
...     pos = pvat.getInterpol(t0)
...     print('%.2f %.2f %.2f' % (pos[0], pos[1], pos[2]))
302.77 2579.72 -70981.57

Use :mod:`irfpy.vexpvat.orbnum` module for converting time and
the orbit number.
'''
import os
import logging
import gzip
import datetime

import numpy as np
import scipy as sp
from scipy import interpolate

import irfpy.util.irfpyrc
from irfpy.util.julday import Julday
from irfpy.util.utc import convert


class PvatFile:
    ''' PVAT file.

    Replacement of :class:`irfpy.vexpvat.pvat.PvatFile`, while
    only file name is allowed.  No uri can be specified.
    '''
    def __init__(self, filename):
        self.filename = filename
        fp = gzip.open(filename)

        self.tlist = []
        self.xlist = []
        self.ylist = []
        self.zlist = []

        for l in fp:
            lin = l.decode('latin-1')
            if lin.startswith('#'):
                continue
            elem = lin.split()
            if len(elem) != 20:
                continue

            t = Julday(float(elem[0])).juld()
            x = float(elem[1])
            y = float(elem[2])
            z = float(elem[3])

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

        self.xspl = interpolate.InterpolatedUnivariateSpline(
                                        self.tlist, self.xlist)
        self.yspl = interpolate.InterpolatedUnivariateSpline(
                                        self.tlist, self.ylist)
        self.zspl = interpolate.InterpolatedUnivariateSpline(
                                        self.tlist, self.zlist)

        fp.close()

    def getInterpol(self, t, extrapol=300):
        ''' Get the position of S/C in VSO frame.

        Interpolated position by spline.

        :param t: Time. Either iterative or non-iterative.
        :type t: Instances allowed in ``irfpy.util.utc.convert``.
        :keyword extrapol: Duration in seconds of the time
            allowed the extrapolation.
            Default is 300 seconds. If either of the given time exceed the
            time range of the data in the datafile by more than the value
            specified in extrapol, ValueError is thrown.
        :returns: Position.
        :rtype: np.array with shape of (3,) or (3, N), where N is len(t).
        '''
        if np.iterable(t):
            tt = t
        else:
            tt = [t, ]

        jd = [convert(ti, irfpy.util.julday.Julday).juld() for ti in tt]

        t0 = min(self.tlist) - extrapol / 86400.
        t1 = max(self.tlist) + extrapol / 86400.

        for jdi in jd:
            if jdi < t0:
                jdi_dt = irfpy.util.julday.Julday(jdi).getDatetime()
                t0_dt = irfpy.util.julday.Julday(t0).getDatetime()
                emsg = (('Extrapolation: Specified time %s is earlier ' +
                        'than the start data time of %s by %s') % (
                        jdi_dt, t0_dt, str(t0_dt - jdi_dt)))
                logging.warn(emsg)
                raise ValueError(emsg)
            elif jdi > t1:
                jdi_dt = irfpy.util.julday.Julday(jdi).getDatetime()
                t1_dt = irfpy.util.julday.Julday(t1).getDatetime()
                emsg = (('Extrapolation: Specified time %s is later ' +
                        'than the end data time of %s by %s') % (
                        jdi_dt, t1_dt, str(jdi_dt - t1_dt)))
                logging.warn(emsg)
                raise ValueError(emsg)

        xx = self.xspl(jd)
        yy = self.yspl(jd)
        zz = self.zspl(jd)

        return np.array([xx, yy, zz])


def getPvatFile(orbnr):
    logger = logging.getLogger('getPvatFile')
    rc = irfpy.util.irfpyrc.Rc()
    filebase = rc.get('vexpvat', 'pvatbase')
    if filebase is None:
        raise RuntimeError('No entry in RC.  vexpvat/pvatbase')
    file = os.path.join(filebase, 'venus_orbit_%04d_pos_sunframe.dat.gz'
                                                                % orbnr)
    logger.debug('File = %s' % file)

    return PvatFile(file)


def isdb():
    rc = irfpy.util.irfpyrc.Rc()
    filebase = rc.get('vexpvat', 'pvatbase')
    if filebase is None:
        return False
    if os.path.isdir(filebase):
        return True
    else:
        return False


import unittest
import doctest


def doctests():
    return unittest.TestSuite((
        doctest.DocTestSuite(),
        ))
if __name__ == '__main__':
    unittest.main(defaultTest='doctests')