#!/usr/bin/env python

# simple script to plot scint_level from Rec35:19 for all satellites/signals in a T0x file

import argparse
parser = argparse.ArgumentParser(description='Plot Legacy Amp Scint Levels for all SVs')
parser.add_argument('filename', help='T04 filename')
parser.add_argument('-p','--png',
                    help='Generate PNG files in headless operation',action="store_true")
parser.add_argument('--sat_types',
                    help='List to sat types to process (e.g., 1,10)',
                    default=None)
parser.add_argument('--prns',
                    help='List to satellites to process (e.g., 1,24)',
                    default=None)
parser.add_argument('--viewdat',
                    help='Additional arguments to pass to viewdat (e.g. --dec=1000)',
                    default="")
parser.add_argument('--elev',
                    help='Minimum elevation [default 0]',
                    default=0., type=float)
parser.add_argument('--min_pts',
                    help='Minimum # of epochs in a segment [default 1000]',
                    default=1000, type=int)
args = parser.parse_args()
if args.min_pts < 60:
    parser.error("--min_pts must be >= 60")

do_sat_types = None
if args.sat_types is not None:
    do_sat_types = [int(x) for x in args.sat_types.split(',')]

do_prns = None
if args.prns is not None:
    do_prns = [int(x) for x in args.prns.split(',')]

print("Elevation mask %.1f deg"%args.elev)

if args.png:
    import matplotlib
    # Allow running headless from the command line
    matplotlib.use("agg")

from mutils import *
from collections import defaultdict

if not 'd' in globals():
    if args.viewdat != "":
        print("Extra viewdat args:" + args.viewdat)
    d = vd2cls(args.filename,rec='-d35:19 '+args.viewdat)

all_figs = {}
sv_list = get_sv_list(d,d.k)
for sv,sat_type in sv_list:
    if do_sat_types is not None and not sat_type in do_sat_types:
        continue
    if do_prns is not None and not sv in do_prns:
        continue

    i = find( (d.SV==sv) & (d.SAT_TYPE==sat_type) )
    d0 = d[i]
    info = sv_to_Lx_info( d0, d0.k, sv, sat_type )
    for freq,track,desc in info:
        i1 = get_Lx_idx( d0, d0.kf, sv, sat_type,
                         freq, track,
                         args.elev, min_seg=args.min_pts )
        if len(i1) == 0:
            print('%s SV %d - too few points'%
                  (desc,sv))
            continue
        if desc in all_figs:
            fig1,ax1 = all_figs[desc]
        else:
            fig1,ax1 = subplots()
            title(desc)
            all_figs[desc] = (fig1,ax1)
            ax1.grid();
        ax1.plot( d0.TIME[i1], (d0.MEAS[i1].astype(int)>>13)&3, label='SV %d'%(sv) )
        ax1.set_xlabel('Time (GPS Sec)')
        ax1.set_ylabel('Scint Level')

def finalize_fig(prefix,fig_info):
    for desc,(fig,ax) in fig_info.items():
        figure(fig.number)
        figlegend()
        make_legend_interactive(fig)
        if args.png:
            desc = desc.replace(' ','_')
            desc = desc.replace('-','_')
            savefig('%s_%s.png'%(prefix,desc))

finalize_fig('scint', all_figs)

if not args.png:
    show()