moon_oplot

Orbit plotter for moon

USAGE: moon_oplot.py [-h] [-v] [-b baseuri] [-o output] yyyy-mm-ddThh-mm-ss

OPTS: -h help

-v

verbose

-b basename

Specify baseuri of the orbit data.

-o output

Specify the output file name. If not specified (default), X window output.

#!/usr/bin/python

'''Orbit plotter for moon

USAGE:  moon_oplot.py [-h] [-v] [-b baseuri] [-o output] yyyy-mm-ddThh-mm-ss

OPTS:     -h      help
          -v      verbose
          -b basename   Specify baseuri of the orbit data.
          -o output     Specify the output file name.  If not specified (default), X window output.
'''

from irfpy.util.julday import Julday
from irfpy.util import utc 
from irfpy.earth.bowshock import BsModelFairfield as BsModel
from irfpy.earth.magnetosphere import MsModelShue as MsModel

from pylab import *
from matplotlib.patches import Arc
from numpy.ma import masked_where

import os,sys
import getopt
import datetime
import dateutil.parser

import logging

import subprocess

import re

from irfpy.cy1orb.MoonPos import MoonGse

def oplot(args, p_figure=None, p_subplot='111'):
    opts, arg = getopt.getopt(args, 'hvb:o:')
    basedir=None
    output=None

    for o,oa in opts:
        if o in ('-h'):
            usage()
            return
        if o in ('-v'):
            logging.basicConfig(level=logging.DEBUG)
        if o in ('-b'):
            logging.debug('Basedir specified.')
            basedir=oa
        if o in ('-o'):
            logging.debug('Output specified.')
            output=oa

    logging.info('Basedir=%s'%basedir)
    logging.info('Output=%s'%output)

    if len(arg) != 1:
        usage()
        return
    else:
        plottime = dateutil.parser.parse(arg[0])

    logging.info('Plotting time is %s'%str(plottime))
    t0 = plottime - datetime.timedelta(15, 0)   # 15 days
    t1 = plottime + datetime.timedelta(15, 0)   # 15 days

    logging.info('Retrieving data from %s'%t0)
    logging.info('Retrieving data to %s'%t1)

    if basedir==None:
        moon=MoonGse()
    else:
        moon=MoonGse(baseurl=basedir)

    moon.setFromDefaultUri()

    pos = moon.getPosition(start=t0, stop=t1)
    rm=6378.
    pos[1:4,:]=pos[1:4,:]/rm   # in Rm.

    print(pos.shape)

    def getNeighborPos(poslist,dt):
      '''Obtain the closest time of poslist to dt.'''
      jd=Julday(dt.year, dt.month, dt.day, dt.hour, dt.minute, dt.second+dt.microsecond/1e6)
      djd = array(abs(poslist[0,:] - jd.juld()))
      minidx = djd.argmin()
      return array(poslist[:,minidx])

    if p_figure == None:
      p_figure=figure()

    # Orbit

    ax=p_figure.add_subplot(p_subplot, aspect='equal', title='Moon-Earth [GSE]', xlabel='Xgse [Rm]', ylabel='Ygse [Rm]')

    ax.plot(pos[1,:], pos[2,:], 'b-')

    posnow = getNeighborPos(pos, plottime)
    ax.plot(posnow[1], posnow[2], 'ro')
    timenow = Julday(posnow[0]).getDatetime()
    ax.text(posnow[1], posnow[2], timenow.strftime('%FT%T'))
    

    # Ticks

    tic = datetime.datetime(t0.year, t0.month, t0.day, 0, 0, 0)+datetime.timedelta(1)
    tic1 = datetime.datetime(t1.year, t1.month, t1.day, 0, 0, 0)
    while tic <= tic1:
      postic = getNeighborPos(pos, tic)
      ax.plot(postic[1],postic[2], 'k+')
      timetic = Julday(postic[0]).getDatetime()
      if timetic.day % 5 == 0:
        ax.text(postic[1], postic[2], timetic.strftime('%FT%T'))
      
      tic = tic+datetime.timedelta(1)



    # Earth

    earth = Arc((0,0), 1, 1, fill=False)
    ax.add_patch(earth)

    # Bowshock/Magnetopause
    bs = BsModel()
    bnd=array(bs.getBounds2D())
    ax.plot(bnd[:,0], bnd[:,1], 'g-.')
    ax.plot(bnd[:,0], -bnd[:,1], 'g-.')

    ms = MsModel()
    bnd=array(ms.getBounds2D())
    ax.plot(bnd[:,0], bnd[:,1], 'g-.')
    ax.plot(bnd[:,0], -bnd[:,1], 'g-.')


    ax.set_xlim(-80,80)
    ax.set_ylim(-80,80)

    return (p_figure, (ax))
      

#    fig, axes = plot_data(timelist, poslist, ticks=[timelist_t, poslist_t], output=output, posttitle='@%04d'%orbnr)

def usage(fd=sys.stderr, exit=False):
    print(__doc__, file=fd)


if __name__ == '__main__':
    if len(sys.argv)==1:
        usage()
        sys.exit(-1)
    f,a=oplot(sys.argv[1:])

    plt.show()