#!/usr/bin/env python

# Copyright (c) 2018-2019 Trimble Inc. ###############################
# $Id: sdp4_predict.py,v 1.3 2019/01/09 20:14:11 acartmel Exp $
######################################################################

usage="""\
Example:
  ./sdp4_predict.py --week 2028 --tow 86400 --file bds34tle
or for all command line options
  ./sdp4_predict.py -h

Plot the predicted elevation for a given SDP4/TLE almanac
"""

##################################
# Imports

from pylab import *
import ephem # pip install pyephem if missing
from datetime import datetime, timedelta
from sdp4_check import read_tle, compute_elev
import argparse


##################################
# Default configuration parameters

def_step = 60
def_dur  = 86400

def_lat = 37.384896686
def_lon = -122.005463678
def_hgt = -6.225

tleline1 = ""#1 43647U 18078A   18317.01550391 -.00000098  00000-0  10000-3 0  9994"
tleline2 = ""#2 43647  55.0057  35.6772 0006416 316.1729 276.9015  1.86232309   525"


##################################
# Command line parsing

parser = argparse.ArgumentParser()
parser.add_argument('--week', default=-1, help='Start GPS week #', type=int)
parser.add_argument('--tow', default=0, help='Start GPS ToW [s]', type=int)
parser.add_argument('--step', default=def_step, help='Time step size [s] (default to '+str(def_step)+'s)',type=int)
parser.add_argument('--dur', default=def_dur, help='Duration [s] (default to '+str(def_dur)+'s)',type=int)
parser.add_argument('--lat', default=def_lat, help='User latitude [deg.] (default to '+str(def_lat)+'deg.)',type=float)
parser.add_argument('--lon', default=def_lon, help='User longitude [deg.] (default to '+str(def_lon)+'deg.)',type=float)
parser.add_argument('--hgt', default=def_hgt, help='User height [m] (default to '+str(def_hgt)+'m)',type=float)
parser.add_argument('--min_el', default=0, help='Min. elevation to display [deg.] (default to 0deg.)', type=int)
parser.add_argument('--tle1', default="", help='TLE line 1 (enclose in double quotes "")')
parser.add_argument('--tle2', default="", help='TLE line 2 (enclose in double quotes "")')
parser.add_argument('--file', default="", help='File containing TLE (overrides --tle)')
parser.add_argument('--sv', default="", help='SV name to display on plot (defaults to empty string)')
parser.add_argument('--save', default="", help='Save the plot as PNG with this file name, leave empty to not save (defaults to no save)')
parser.add_argument('--plot_sec', help='Show seconds on plot instead of hours', action="store_true")

args = parser.parse_args()

# Take the TLE lines from a file by first choice
# Otherwise assumed supplied from --tle1 "" and --tle2 ""
if (args.file != ""):
  tlecnt = 1
  with open(args.file) as fp:
    for line in fp:
      if (tlecnt == 1 and line.find("1 ") == 0):
        tleline1 = line
        tlecnt = 2
      elif (tlecnt == 2 and line.find("2 ") == 0):
        tleline2 = line
        tlecnt = -1
else:
  if (args.tle1 != ""):
    tleline1 = args.tle1
  if (args.tle2 != ""):
    tleline2 = args.tle2


##################################
# Compute the predicted elevation

svelev = []
taxis  = []

svtle = read_tle( [tleline1, tleline2] )
gps_start = datetime(1980, 1, 6)

if args.week >= 0:
  gps_sec = args.week*7*24*60*60 + args.tow
  start_time = gps_start + timedelta(seconds = gps_sec)
else:
  start_time = datetime.utcnow()
prev_elev = -1
max_elev = -1
rise_time = None
for pred_offset in range(0, args.dur, args.step):
  pred_time = start_time + timedelta(seconds = pred_offset)
  elev = compute_elev(svtle, pred_time, lat=args.lat, lon=args.lon, hgt=args.hgt)
  if (elev < args.min_el):
    elev = None
  if elev is not None:
    if prev_elev < 0:
      rise_time = pred_time
      print("Start of pass UTC:",pred_time)
    if elev < prev_elev and max_elev < 0:
      max_elev = elev
    prev_elev = elev
  elif prev_elev >= 0:
    print("  End max elev %.1f, length"%max_elev,(pred_time-rise_time))
    max_elev = -1
    prev_elev = -1
  svelev.append(elev)
  taxis.append(pred_time)

if args.plot_sec:
  secs_in_week = 7*86400
  taxis = [(x - gps_start).total_seconds()%secs_in_week for x in taxis]

##################################
# Plot the predicted elevation

figure(figsize=(12,6))
plot(taxis, svelev)

tmp = 'Elevation Plot'
if (args.sv != ""):
  tmp = args.sv+' '+tmp
title(tmp)
xlabel('UTC time')
ylabel('Elevation [deg]')

if (args.min_el<0):
  ymin = args.min_el
else:
  ymin = 0
axis([taxis[0], taxis[-1], ymin, 90])
grid()

if (args.save != ""):
  savefig(args.save, format='png')
  print('Plot saved as '+args.save+', no display')
else:
  show()