import mutils as m
import sys
from pylab import *

if len(sys.argv) != 5:
    print("""
Show SNR and cycle slips for highest N satellites with L2C.
Used mainly to compare a 0-baseline receiver to SoftGNSS.

Usage: From shell, run "ipython --pylab" then:
  %run diag_snr.py N do_L2 filename1.T0x filename2.T0x
where:
  N = maximum # of plots to generate
  do_L2 = if non-zero, then analyze L2 data as well as L1
  filename1.T0x is the reference 0-baseline file
  filename2.T0x is the test SoftGNSS file
Example (note - the '%' character below is needed):
  %run diag_snr.py 2 receiver.T04 test.T02
""")
    sys.exit(1)


print("Loading data")
N_max = int(sys.argv[1])
do_L2 = int(sys.argv[2])
(ref,k)=m.vd2arr(sys.argv[3],rec='-d27',opt=['--translate_rec35_sub9_to_rec27'])
(tst,k)=m.vd2arr(sys.argv[4],rec='-d27',opt=['--translate_rec35_sub9_to_rec27'])

iGPS=find(ref[:,k.SAT_TYPE]==0)
sv_list = unique( ref[iGPS,k.SV] )
el = []
for sv in sv_list:
    irL1 = find( (ref[:,k.SV]==sv) & (ref[:,k.SAT_TYPE]==0) & (ref[:,k.FREQ]==0) )
    irL2 = find( (ref[:,k.SV]==sv) & (ref[:,k.SAT_TYPE]==0) & (ref[:,k.FREQ]==1) & (ref[:,k.TRACK]>=3) )
    if len(irL1) > 10:
        L1mean = mean( ref[irL1,k.EL] )
        if do_L2:
            if len(irL2) > 0:
                L2C = 1
                el.append( [L1mean, L2C, sv] )
        else:
                el.append( [L1mean, 0, sv] )

el = array(el)
el2 = el[el[:,0].argsort(),:]

sv_list = el2[-1::-1,2]

print("Plotting data")
n = 0
for sv in sv_list:
    n += 1
    if n > N_max:
        break

    irL1 = find( (ref[:,k.SV]==sv) & (ref[:,k.SAT_TYPE]==0) & (ref[:,k.FREQ]==0) )
    itL1 = find( (tst[:,k.SV]==sv) & (tst[:,k.SAT_TYPE]==0) & (tst[:,k.FREQ]==0) )
    [tmp, ir, it] = m.intersect( ref[irL1,k.TIME], tst[itL1,k.TIME] )
    irL1 = irL1[ir]
    itL1 = itL1[it]

    if do_L2:
        irL2 = find( (ref[:,k.SV]==sv) & (ref[:,k.SAT_TYPE]==0) & (ref[:,k.FREQ]==1) & (ref[:,k.TRACK]==5) )
        itL2 = find( (tst[:,k.SV]==sv) & (tst[:,k.SAT_TYPE]==0) & (tst[:,k.FREQ]==1) & (tst[:,k.TRACK]==5) )
        [tmp, ir, it] = m.intersect( ref[irL2,k.TIME], tst[itL2,k.TIME] )
        irL2 = irL2[ir]
        itL2 = itL2[it]

    figure()
    if do_L2:
        ax = subplot(211)
    ir = find(ref[irL1,k.SLIP_FLG])
    it = find(tst[itL1,k.SLIP_FLG])
    plot(ref[irL1,k.TIME], ref[irL1,k.CNO], 'x', ref[irL1[ir],k.TIME], ref[irL1[ir],k.CNO], 'ro', \
         tst[itL1,k.TIME], tst[itL1,k.CNO], 'x', tst[itL1[it],k.TIME], tst[itL1[it],k.CNO], 'mo')
    title('L1: sv %d elev %.1f. slips ref=%d tst=%d' % (sv, mean(ref[irL1,k.EL]),len(ir),len(it)))
    ylabel('SNR')

    if do_L2:
        subplot(212,sharex=ax)
        ir = find(ref[irL2,k.SLIP_FLG])
        it = find(tst[itL2,k.SLIP_FLG])
        plot(ref[irL2,k.TIME], ref[irL2,k.CNO], 'x', ref[irL2[ir],k.TIME], ref[irL2[ir],k.CNO], 'ro', \
             tst[itL2,k.TIME], tst[itL2,k.CNO], 'x', tst[itL2[it],k.TIME], tst[itL2[it],k.CNO], 'mo')
        xlabel('Time [s]')
        ylabel('SNR')
        title('L2.  slips ref=%d tst=%d' % (len(ir),len(it)))

show()