apps130226_orbit_finder.sara_wake_earth_facing

Look for the time when SARA in wake and Earth facing.

It is to understand if the auroral ENA affect the observation or not.

''' Look for the time when SARA in wake and Earth facing.

It is to understand if the auroral ENA affect the observation or not.
'''
import time
import datetime
import logging
logging.basicConfig()
logging.getLogger().setLevel(logging.INFO)
import pickle as pickle

import numpy as np
import matplotlib.dates as mpldates

from irfpy.util.timeseries import dtrange
from irfpy.cy1orb.pvat2 import getmepos, getlsepos, getdatarange

from irfpy.util import timeinterval

def main_pole():
    ''' Return the list of time interval when the s/c is at pole
    '''

    # Here is a template...
    dt = datetime.timedelta(minutes=1)
    trng = getdatarange()
    t0 = mpldates.num2date(min(trng.interval.boundary.args)).replace(tzinfo=None)
    t1 = mpldates.num2date(max(trng.interval.boundary.args)).replace(tzinfo=None)
    print(trng, trng.interval.boundary.args)
    print(t0, t1)
#    t0 = datetime.datetime(2009, 2, 6)
#    t1 = datetime.datetime(2009, 2, 6, 10)
    dth = dt / 2
    tlist = dtrange(t0, t1, dt)
    xyzlist = np.zeros([len(tlist), 3])

    t0 = time.time()
    for i, t in enumerate(tlist):
        if i % 1000 == 0:
            logging.info('Reading data.   %d / %d (%.2f s ellapsed)' % (i, len(tlist), time.time() - t0))

        pos = getmepos(t)    # You may also use getlsepos for LSE condition.
                             # Hybrid could be used with a small effort.
        xyzlist[i, 0] = pos[0]
        xyzlist[i, 1] = pos[1]
        xyzlist[i, 2] = pos[2]

    ## After reading data, you may use the interval set for storing the OK time.
    tint = timeinterval.empty_interval()

    t0 = time.time()
    for i, txyz in enumerate(list(zip(tlist, xyzlist))):
        if i % 1000 == 0:
            logging.info('Judging data.   %d / %d (%.2f s ellapsed)' % (i, len(tlist), time.time() - t0))

        t, (lon, lat, hei) = txyz
        
#        print t, lon, lat, hei
        # Add condition here.
        if lat > 60 or lat < -60:
            tint = tint + timeinterval.timeinterval(t - dth, t + dth)

    print(tint)
    pickle.dump(tint, open('time_interval_pole.pickle.txt', 'wb'))

    return tint

def main_wake():
    '''Main script'''

    dt = datetime.timedelta(minutes=1)

    trng = getdatarange()
    t0 = mpldates.num2date(min(trng.interval.boundary.args)).replace(tzinfo=None)
    t1 = mpldates.num2date(max(trng.interval.boundary.args)).replace(tzinfo=None)
    dth = dt / 2

    tlist = dtrange(t0, t1, dt)
    print(len(tlist))

    xyzlist = np.zeros([len(tlist), 3])

    for i, t in enumerate(tlist):
        if i % 1000 == 0: print(t)

        pos = getlsepos(t)
        xyzlist[i, 0] = pos[0]
        xyzlist[i, 1] = pos[1]
        xyzlist[i, 2] = pos[2]

    xneg = timeinterval.empty_interval()
    rin = timeinterval.empty_interval()

    for i, txyz in enumerate(list(zip(tlist, xyzlist))):
        t, (x, y, z) = txyz
        print(i, t, x, y, z)

        ### Where x-negative
        if x < 0:
            xneg = xneg + timeinterval.timeinterval(t - dth, t + dth)

    for i, txyz in enumerate(list(zip(tlist, xyzlist))):
        t, (x, y, z) = txyz
        print(i, t, x, y, z)

        ### Where r < 1738
        if y * y + z * z < 1738 * 1738.:
            rin = rin + timeinterval.timeinterval(t - dth, t + dth)



    print('X negative', xneg)
    print('R in radius', rin)

    print('Wake', xneg.intersection(rin))

    wake = xneg.intersection(rin)
        
    pickle.dump(wake, open('time_interval_wake.pickle.txt', 'w'))
    

if __name__ == "__main__":
    main_pole()
    main_wake()