#!/usr/bin/env python

#######################################################
# Copyright Trimble Inc 2017
#######################################################

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

from numpy import *
from pylab import *
import mutils as m
import os;
import sys;
import argparse

# Override pylab find() to avoid deprecation warning
# pylint: disable=function-redefined
def find(x):
    return where(x)[0]

#######################################################
# Start of script
#######################################################
if __name__ == "__main__":

  ######################################################################
  # Parse arguments
  parser = argparse.ArgumentParser(description='Parse positions from spoof test')
  parser.add_argument('filename', help='T02/T04 filename')
  parser.add_argument('sched', help='CSV file providing the spoof schedule')
  parser.add_argument('spoofDir', help='subdirectory where we\'ll save the results')
  parser.add_argument('-t','--truth_lla', help='Truth lat,lon,hgt, e.g. -t 37.3848990528,-122.0054260639,-6.260')
  parser.add_argument('-r','--record35_16', help='Analyze the record 35:16 instead of the 29 (T02 files) or 35:2 (T04)', action="store_true")
  args = parser.parse_args()
  ######################################################################

  filename = args.filename
  sched    = args.sched
  spoofDir = args.spoofDir
  
  cwd = os.getcwd()
  print(cwd)

  if not os.path.exists(spoofDir):
    os.makedirs(spoofDir)

  # ToDo: update viewdat so we can automatically convert record 29's to 35:2's, then
  # all we will need is something like the following
  #(ref,k)=m.vd2arr(filename,rec='-d35:2',opt=['--dec=1000 --translate_rec29_to_rec35_sub2'])

  # First remove all the records except position and obs records. For each
  # segment we will call t012kml so by removing the non-position/obs records
  # once we should be able to complete the processing faster. We need
  # the obs record to fill in the tracking table in the KMZ
  # placemarkers. What we'll remove:
  #
  # - FFT data
  # - CPU load data
  # - FFT based ACF data
  # - ephemeris
  # - nav/SBAS bits
  # - diagnostic records
  # - error / warnings
  # - received corrections
  # - inertial records
  # - key/value pair data
  # - ...
  name, file_extension = os.path.splitext(filename)
  if(file_extension == ".T02" or file_extension == ".t02"):
    os.system("t0x2t0x -i29,35:2,35:16,27,35:19 " + filename + " " + name + ".t04")
    filename = name + ".t04"
    datafile = filename
  else:
    if(args.record35_16):
      os.system("t0x2t0x -i35:16,27,35:19 " + filename + " " + name + "-tmp.t04")
    else:
      os.system("t0x2t0x -i29,35:2,27,35:19 " + filename + " " + name + "-tmp.t04")
    datafile = name + "-tmp.t04"

  if(args.record35_16):
    print("Extracting 35:16")
    (ref,k)=m.vd2arr(datafile,rec='-d35:16')
  else:
    print("Extracting 35:2")
    (ref,k)=m.vd2arr(datafile,rec='-d35:2')

  Time    = ref[:,k.TIME]
  used    = ref[:,k.NUSED]
  tracked = ref[:,k.NTRK]

  GPSUsd   = ref[:,k.GPSUsd]
  SBASUsd  = ref[:,k.SBASUsd]
  GLOUsd   = ref[:,k.GLOUsd]
  GALUsd   = ref[:,k.GALUsd]
  QZSSUsd  = ref[:,k.QZSSUsd]
  BDSUsd   = ref[:,k.BDSUsd]
  IRNSSUsd = ref[:,k.IRNSSUsd]

  ClkOffset = ref[:,k.CLKBIAS]
  ClkDrift  = ref[:,k.CLKDRIFT]
  GPSClkDelta = ref[:,k.GPSCLKDELTA]
  GLNClkDelta = ref[:,k.GLNCLKDELTA]
  GALClkDelta = ref[:,k.GALCLKDELTA]
  SBASClkDelta = ref[:,k.SBASCLKDELTA]
  BDSClkDelta = ref[:,k.BDSCLKDELTA]
  
  HorzEstErr = sqrt(ref[:,k.SigE]**2 + ref[:,k.SigN]**2)

  LLH = ref[:,k.LAT:k.LAT+3]
  
  dPP_TIME=0
  dPP_LAT=1

  if args.truth_lla:
    lla_str = args.truth_lla
    LLH_Ref = array([float(x) for x in lla_str.split(',')])
    TimeSecs = Time
    print(LLH_Ref)
  else:
    # No truth - use the median as it protects to a first order
    # against outliers
    LLH_Ref = array([median(LLH[:,0]), median(LLH[:,1]), median(LLH[:,2])])
    TimeSecs = Time
    print(LLH_Ref)

  (dE, dN, dU) = m.llh2enu( LLH, LLH_Ref, False, False)
  err_2d = sqrt( dE**2 + dN**2 )
  err_3d = sqrt( dE**2 + dN**2 + dU**2 )
  colArray = ['b','r','g','m','c','k','#FF8C00','#FF99FF','#COCOCO','#CCFFCC','#FFCCCC'];
  DayNum = floor(TimeSecs[0]/86400)


  timeVsHour = 0

  if(timeVsHour == 1):
    Time = Time/3600.0 - DayNum*24

  statsfile = open(filename + '-stats.csv', 'w')
  htmlfile  = open(filename + '-stats.html','w')
  htmlfile.write("<html><head>\n")
  htmlfile.write("<title>%s</title>\n" % filename)
  htmlfile.write("<link rel='stylesheet' href='../../js/style.css' type='text/css'/>\n")
  htmlfile.write("<script type='text/javascript' src='../../js/jquery-latest.js'></script>\n")
  htmlfile.write("<script type='text/javascript' src='../../js/jquery.tablesorter.js'></script>\n")
  htmlfile.write("<!-- Needs some work -->")
  htmlfile.write("<!--\n")
  htmlfile.write("<script type='text/javascript'>")
  htmlfile.write("$(document).ready(function() {\n")
  htmlfile.write("  $.tablesorter.defaults.sortList = [[0,0]];\n")
  htmlfile.write("  $('table').tablesorter({\n")
  htmlfile.write("    textExtraction: function(node) {\n")
  htmlfile.write("      return node.childNodes[0].innerHTML;\n")
  htmlfile.write("} }); } );\n")
  htmlfile.write("</script>")
  htmlfile.write("-->\n")
  htmlfile.write("</head>\n")

  htmlfile.write("<body>\n")
  if("BX940" in cwd):
    tokens = cwd.split("/")
    htmlfile.write("<a href='Freq1583-%s-Freq2D-120s.png'>L1 Spectrum History</a>\n" % tokens[-2])

  htmlfile.write("<table class='tablesorter'>\n")
  htmlfile.write("<thead>\n")
  htmlfile.write("<tr>\n")
  htmlfile.write("<th class =\"{sorter: 'custom'}\" rowspan='2'>Segment</th>")
  if("Loci" in filename):
    htmlfile.write("<th rowspan='2'>History</th>\n")
  htmlfile.write("<th rowspan='2'>Start</th>\n")
  htmlfile.write("<th rowspan='2'>Stop</th>\n")
  htmlfile.write("<th rowspan='2'>Epochs</th>\n")
  htmlfile.write("<th rowspan='2'>Yield [%]</th>\n")
  htmlfile.write("<th rowspan='2'>2D HMI* [%]</th>\n")
  htmlfile.write("<th rowspan='2'>Alarm** [%]</th>\n")
  if("BX940" in cwd):
    htmlfile.write("<th rowspan='2'>ACF</th>\n")
  htmlfile.write("<th rowspan='2'>KMZ</th>\n")
  htmlfile.write("<th colspan='3'>Clock</th>\n")
  htmlfile.write("<th rowspan='2'>ENU</th>\n")
  if(args.truth_lla):
    htmlfile.write("<th rowspan='2'>ErrVsEst</th>\n")
  htmlfile.write("<th colspan=7>Mean SVs Used</th>\n")
  htmlfile.write("<th colspan=5>2D Percentiles [m]</th>\n")
  htmlfile.write("<th colspan=5>abs(Hgt) Percentiles [m]</th>\n")
  htmlfile.write("</tr>\n");
  
  htmlfile.write("<tr>\n")
  htmlfile.write("<th>Bias</th>\n");
  htmlfile.write("<th>Drift</th>\n");
  htmlfile.write("<th>Sys</th>\n");
  htmlfile.write("<th>Total</th>\n");
  htmlfile.write("<th>GPS</th>\n");
  htmlfile.write("<th>SBAS</th>\n");
  htmlfile.write("<th>GLN</th>\n");
  htmlfile.write("<th>GAL</th>\n");
  htmlfile.write("<th>QZSS</th>\n");
  htmlfile.write("<th>BDS</th>\n");
  
  htmlfile.write("<th>50%</th>\n");
  htmlfile.write("<th>68%</th>\n");
  htmlfile.write("<th>90%</th>\n");
  htmlfile.write("<th>95%</th>\n");
  htmlfile.write("<th>100%</th>\n");
  
  htmlfile.write("<th>50%</th>\n");
  htmlfile.write("<th>68%</th>\n");
  htmlfile.write("<th>90%</th>\n");
  htmlfile.write("<th>95%</th>\n");
  htmlfile.write("<th>100%</th>\n");
  htmlfile.write("</tr>\n");
  htmlfile.write("</thead>\n")
  
  htmlfile.write("<tbody>\n")

  with open(sched) as fid:
    for line in fid:
      tokens = line.rstrip().split(",")

      segment = tokens[0]
      week    = int(tokens[1])
      start   = float(tokens[2])
      stop    = float(tokens[3])

      # As we may be processing SNPC data that does not include the
      # string "BX992" (which was 10Hz) check the full path as BX992
      # is in the path where the data is located and we're operating

      #if(filename.startswith('BX992')):
      if("BX992" in cwd):
        # 10Hz
        expected = (stop - start) * 10
      else:
        # 5Hz
        expected = (stop - start) * 5

      print(segment,start,stop)

      i = find( (TimeSecs >= start) & (TimeSecs < stop) )


      HourStart = remainder(floor(start / 3600), 24)
      HourStop  = remainder(floor(stop  / 3600), 24)
      MinStart  = remainder(floor(start / 60), 60)
      MinStop   = remainder(floor(stop  / 60), 60)
      SecStart  = remainder(start, 60)
      SecStop   = remainder(stop, 60)

      statsfile.write("%s,%d,%d," % (segment,start, stop))
      statsfile.write("%d:%02d:%02d," % (HourStart,MinStart,SecStart))
      statsfile.write("%d:%02d:%02d," % (HourStop,MinStop,SecStop))

      if(len(i) == 0):
        statsfile.write("0,0.0,NaN,NaN,")
        statsfile.write("NaN,NaN,NaN,NaN,NaN,")
        statsfile.write("NaN,NaN,NaN,NaN,NaN,")
        statsfile.write("0,0,0,0,0,0,\n")
      
        #######
        ##  HTML File
        #######

        htmlfile.write("<tr><th>%s</th>" % (segment))
        if("Loci" in filename):
          htmlfile.write("<td><a href='Trend/Trend-%s-%s.html'>X</a></td>" % (segment,filename))
      
        htmlfile.write("<td>%.0f</td><td>%.0f</td><td>%d</td><td>0.0</td>" % (start, stop, len(i)))

        htmlfile.write("<td>NaN</td><td>NaN</td>")
        if("BX940" in cwd):
          htmlfile.write("<td><a href='%s.mp4'>X</a></td>" % segment)

        htmlfile.write("<td>&nbsp;</td>")
        htmlfile.write("<td>&nbsp;</td>")
        htmlfile.write("<td>&nbsp;</td>")
        htmlfile.write("<td>&nbsp;</td>")
        htmlfile.write("<td>&nbsp;</td>")
        htmlfile.write("<td>&nbsp;</td>")
      
        # Mean num SVs in the solution
        htmlfile.write("<td>0</td>")
        htmlfile.write("<td>0</td>")
        htmlfile.write("<td>0</td>")
        htmlfile.write("<td>0</td>")
        htmlfile.write("<td>0</td>")
        htmlfile.write("<td>0</td>")
        htmlfile.write("<td>0</td>")
        
        # 2D Stats
        htmlfile.write("<td>NaN</td>")
        htmlfile.write("<td>NaN</td><td>NaN</td><td>NaN</td><td>NaN</td>")
        # Hgt stats
        htmlfile.write("<td>NaN</td>")
        htmlfile.write("<td>NaN</td><td>NaN</td><td>NaN</td><td>NaN</td>")

      else:
        posYield = 100 * (len(i) / expected)

        cx_2d,cy_2d = m.docdf( err_2d[i] )
        i_50_2d = find( cy_2d >= .5 )[0]
        i_68_2d = find( cy_2d >= .68 )[0]
        i_90_2d = find( cy_2d >= .9 )[0]
        i_95_2d = find( cy_2d >= .95 )[0]
        i_98_2d = find( cy_2d >= .98 )[0]
        print('Error:       50%    68%    90%    95%   100%')
        print('2D      : %6.2f %6.2f %6.2f %6.2f %6.2f' % (cx_2d[i_50_2d],cx_2d[i_68_2d],cx_2d[i_90_2d],cx_2d[i_95_2d],cx_2d[-1]))


        cx_hgt,cy_hgt = m.docdf( abs(dU[i]) )
        i_50_hgt = find( cy_hgt >= .5 )[0]
        i_68_hgt = find( cy_hgt >= .68 )[0]
        i_90_hgt = find( cy_hgt >= .9 )[0]
        i_95_hgt = find( cy_hgt >= .95 )[0]
        i_98_hgt = find( cy_hgt >= .98 )[0]
        print('abs(Hgt): %6.2f %6.2f %6.2f %6.2f %6.2f' % (cx_hgt[i_50_hgt],cx_hgt[i_68_hgt],cx_hgt[i_90_hgt],cx_hgt[i_95_hgt],cx_hgt[-1]))

        cx_3d,cy_3d = m.docdf( err_3d[i] )
        i_50_3d = find( cy_3d >= .5 )[0]
        i_68_3d = find( cy_3d >= .68 )[0]
        i_90_3d = find( cy_3d >= .9 )[0]
        i_95_3d = find( cy_3d >= .95 )[0]
        print('3D      : %6.2f %6.2f %6.2f %6.2f %6.2f' % (cx_3d[i_50_3d],cx_3d[i_68_3d],cx_3d[i_90_3d],cx_3d[i_95_3d],cx_3d[-1]))

        figure()
        plot( Time[i], used[i],     label='Total Used',color=colArray[0])
        plot( Time[i], GPSUsd[i],   label='GPS',color=colArray[1])
        plot( Time[i], SBASUsd[i],  label='SBAS',color=colArray[2])
        plot( Time[i], GLOUsd[i],   label='GLONASS',color=colArray[3])
        plot( Time[i], GALUsd[i],   label='GALILEO',color=colArray[4])
        plot( Time[i], QZSSUsd[i],  label='QZSS',color=colArray[5])
        plot( Time[i], BDSUsd[i],   label='BDS',color=colArray[6])
        plot( Time[i], IRNSSUsd[i], label='IRNSS',color=colArray[7])

        xlabel('Time [Hour]')
        ylabel('Num Satellites')
        grid(True)
        legend(loc='best')
        tight_layout()
        show()
        # Save the data as a PNG file
        savefig(spoofDir + '/' + segment + '-' + filename + ".SVsUsed.png",dpi=150)
        close()

        # Set the ENU and 2D to approx the 95% horizontal
        Lim2D = cx_2d[i_95_2d]
        if(Lim2D < 0.9):
          Lim2D = ceil(Lim2D*10)/10
        else:
          Lim2D = ceil(Lim2D)
      
        LimENU = cx_2d[i_98_2d];
        if(LimENU < cx_hgt[i_98_hgt]):
          LimENU = cx_hgt[i_98_hgt]

        if(LimENU < 0.9):
          LimENU = ceil(LimENU*10)/10
        else:
          LimENU = ceil(LimENU)

  
        # Set the height to approx 95%
        LimHgt = cx_hgt[i_95_hgt]
        if(LimHgt < 0.9):
          LimHgt = ceil(LimHgt*10)/10
        else:
          LimHgt = ceil(LimHgt)
  

        figure()
        plot( Time[i], err_2d[i], label='2D Err',   color=colArray[0])
        plot( Time[i], HorzEstErr[i], label='Est Err',   color=colArray[1])
        tight_layout()
        show()
        savefig(spoofDir + '/' + segment + '-' + filename + ".ErrVsEst.png",dpi=150)
        close()

        Alarm = 20; # 20m alarm limit
        # For GPS a failure is if the error in the satellite is greater
        # than 4.42 times the URA. Let's use the same metric
        ErrorMult = 4.42 # Same as the GPS SPS Performance Standard (4th Ed, 2008)
        j = find( (HorzEstErr[i] < Alarm) & (err_2d[i] > (ErrorMult*HorzEstErr[i])) )
        if(len(j) > 0):
          figure()
          plot( Time[i[j]], err_2d[i[j]], label='2D Err',   color=colArray[0])
          plot( Time[i[j]], HorzEstErr[i[j]], label='Est Err',   color=colArray[1])
          tight_layout()
          show()
          savefig(spoofDir + '/' + segment + '-' + filename + ".ErrVsEst4Sig.png",dpi=150)
          close()
          HMI = 100.0 * float(len(j)) / float(len(i))
          print("Num Err = %d | Total = %d | HMI = %.3f" % (len(j), len(i), HMI))
        else:
          HMI = 0.0
  
        k = find(HorzEstErr[i] >= Alarm)
        if(len(k) > 0):
          AlarmRate = 100.0 * float(len(k)) / float(len(i))
          print("Num Alarm = %d | Total = %d | Alarm = %.3f" % (len(k), len(i), AlarmRate))
        else:
          AlarmRate = 0.0
      
        figure()
        plot( Time[i], dE[i], label=r'East $\sigma$=%.3fm'%std(dE[i]),   color=colArray[0])
        plot( Time[i], dN[i], label=r'North $\sigma$=%.3fm'%std(dN[i]),  color=colArray[1])
        plot( Time[i], dU[i], label=r'Height $\sigma$=%.3fm'%std(dU[i]), color=colArray[2])
        ylim([-LimENU,LimENU])

        xlabel('Time [Hour]')
        ylabel('Error [m]')

        if args.truth_lla:
          title('Relative to Truth : 2D 95%% = %.3fm' % cx_2d[i_95_2d])
        else:
          title('Relative to Median');

        grid(True)
        legend(loc='best')
        tight_layout()
        show()
        # Save the data as a PNG file
        savefig(spoofDir + '/' + segment + '-' + filename + ".ENUErr.png",dpi=150)
        close()

        figure()
        plot( Time[i], ClkOffset[i], label='ClockOffset')
        xlabel('Time [Hour]')
        ylabel('Receiver Clock Offset [ms]')
        grid(True)
        legend(loc='best')
        tight_layout()
        show()
        # Save the data as a PNG file
        savefig(spoofDir + '/' + segment + '-' + filename + ".ClkOffset.png",dpi=150)
        close()

        figure()
        plot( Time[i], ClkDrift[i], label='ClockDrift')
        xlabel('Time [Hour]')
        ylabel('Receiver Clock Drift [ppm]')
        grid(True)
        legend(loc='best')
        tight_layout()
        show()
        # Save the data as a PNG file
        savefig(spoofDir + '/' + segment + '-' + filename + ".ClkDrift.png",dpi=150)
        close()


        clks = concatenate((GPSClkDelta[i],GLNClkDelta[i], GALClkDelta[i], SBASClkDelta[i], BDSClkDelta[i]), axis=0)
        j = find(isfinite(clks))
        if(len(j) > 0):
          maxDelta = 10 * ceil(max(clks[j])/10)
          minDelta = 10 * floor(min(clks[j])/10)

          plot( Time[i], GPSClkDelta[i],  label='GPS')
          plot( Time[i], GLNClkDelta[i],  label='GLN')
          plot( Time[i], GALClkDelta[i],  label='GAL')
          plot( Time[i], SBASClkDelta[i], label='SBAS')
          plot( Time[i], BDSClkDelta[i],  label='BDS')
          ylim([minDelta,maxDelta])
          xlabel('Time [Hour]')
          ylabel('Reference Sys - Sys Clock Offset [ns]')
          grid(True)
          legend(loc='best')
          tight_layout()
          show()
          # Save the data as a PNG file
          savefig(spoofDir + '/' + segment + '-' + filename + ".RefMinusSysClkDelta.png",dpi=150)
          close()

        figure()
        plot( cx_2d, cy_2d, label='2D Error%c68%%=%.3fm%c95%%=%.3fm%c100%%=%.3fm'%(10,cx_2d[i_68_2d],10,cx_2d[i_95_2d],10,cx_2d[-1]),   color=colArray[0])
        xlim([0,Lim2D])
        ylim([0,1])
        xlabel('Error [m]')
        ylabel('CDF')
        if args.truth_lla:
          title('Relative to Truth');
        else:
          title('Relative to Median');
        grid(True)
        legend(loc='best')
        tight_layout()
        show()
        # Save the data as a PNG file
        savefig(spoofDir + '/' + segment + '-' + filename + ".2DCDF.png",dpi=150)
        close()

        figure()
        plot( cx_hgt, cy_hgt, label='abs(Hgt) Error%c68%%=%.3fm%c95%%=%.3fm%c100%%=%.3fm'%(10,cx_hgt[i_68_hgt],10,cx_hgt[i_95_hgt],10,cx_hgt[-1]), color=colArray[0])
        xlim([0,LimHgt])
        ylim([0,1])
        xlabel('Error [m]')
        ylabel('CDF')
        if args.truth_lla:
          title('Relative to Truth');
        else:
          title('Relative to Median');
        grid(True)
        legend(loc='best')
        tight_layout()
        show()
        # Save the data as a PNG file
        savefig(spoofDir + '/' + segment + '-' + filename + ".HgtCDF.png",dpi=150)
        close()

        statsfile.write("%d,%.3f,%.3f,%.3f," % ( len(i), posYield, HMI, AlarmRate))
      
        # 2D Stats
        statsfile.write("%.2f,%.2f,%.2f,%.2f,%.2f," % (cx_2d[i_50_2d],cx_2d[i_68_2d],cx_2d[i_90_2d],cx_2d[i_95_2d],cx_2d[-1]))
        # Hgt stats
        statsfile.write("%.2f,%.2f,%.2f,%.2f,%.2f," % (cx_hgt[i_50_hgt],cx_hgt[i_68_hgt],cx_hgt[i_90_hgt],cx_hgt[i_95_hgt],cx_hgt[-1]))
      
        # Mean num SVs in the solution
        statsfile.write("%.2f,%.2f,%.2f,%.2f,%.2f,%.2f,%.2f\n" % (mean(used[i]), 
                                                                  mean(GPSUsd[i]), 
                                                                  mean(SBASUsd[i]),
                                                                  mean(GLOUsd[i]),
                                                                  mean(GALUsd[i]),
                                                                  mean(QZSSUsd[i]),
                                                                  mean(BDSUsd[i])))

      
        #######
        ##  HTML File
        #######

        htmlfile.write("<tr><th>%s</th>" % (segment))
        if("Loci" in filename):
          htmlfile.write("<td><a href='Trend/Trend-%s-%s.html'>X</a></td>" % (segment,filename))
      
        htmlfile.write("<td>%.0f</td><td>%.0f</td><td>%d</td><td>%.2f</td>" % (start, stop, len(i), posYield))

        htmlfile.write("<td>%.2f</td><td>%.2f</td>" % (HMI,AlarmRate))
        #if(filename.startswith('BX940')):
        if("BX940" in cwd):
          htmlfile.write("<td><a href='%s.mp4'>X</a></td>" % segment)

        htmlfile.write("<td><a href='%s/%s-%s.kmz'>X</a></td>" % (spoofDir,segment,filename))
        htmlfile.write("<td><a href='%s/%s'>X</a></td>" % (spoofDir, segment + '-' + filename + '.ClkOffset.png'))
        htmlfile.write("<td><a href='%s/%s'>X</a></td>" % (spoofDir, segment + '-' + filename + '.ClkDrift.png'))
        htmlfile.write("<td><a href='%s/%s'>X</a></td>" % (spoofDir, segment + '-' + filename + '.RefMinusSysClkDelta.png'))
        htmlfile.write("<td><a href='%s/%s'>X</a></td>" % (spoofDir, segment + '-' + filename + '.ENUErr.png'))
      
        if(args.truth_lla):
          if("SPS461" in cwd):
            htmlfile.write("<td><a href='%s/%s'>X</a></td>" % (spoofDir, segment + '-' + name + '.T02-errPlot.png'))
          else:
            htmlfile.write("<td><a href='%s/%s'>X</a></td>" % (spoofDir, segment + '-' + filename + '-errPlot.png'))

        # Mean num SVs in the solution
        htmlfile.write("<td><a href='%s/%s'>%.2f</a></td>" % (spoofDir, segment + '-' + filename + '.SVsUsed.png', mean(used[i])) )
        #htmlfile.write("<td>%.2f</td><td>%.2f</td><td>%.2f</td><td>%.2f</td><td>%.2f</td><td>%.2f</td>" % 
        #    (mean(GPSUsd[i]),mean(SBASUsd[i]),mean(GLOUsd[i]),mean(GALUsd[i]),mean(QZSSUsd[i]),mean(BDSUsd[i])))

        htmlfile.write("<td><a href='%s/%s'>%.2f</a></td>" % (spoofDir, segment + '-' + name + '.T04-GPS-RAIM.png',  mean(GPSUsd[i])))
        htmlfile.write("<td><a href='%s/%s'>%.2f</a></td>" % (spoofDir, segment + '-' + name + '.T04-SBAS-RAIM.png', mean(SBASUsd[i])))
        htmlfile.write("<td><a href='%s/%s'>%.2f</a></td>" % (spoofDir, segment + '-' + name + '.T04-GLN-RAIM.png',  mean(GLOUsd[i])))
        htmlfile.write("<td><a href='%s/%s'>%.2f</a></td>" % (spoofDir, segment + '-' + name + '.T04-GAL-RAIM.png',  mean(GALUsd[i])))
        htmlfile.write("<td><a href='%s/%s'>%.2f</a></td>" % (spoofDir, segment + '-' + name + '.T04-QZSS-RAIM.png', mean(QZSSUsd[i])))
        htmlfile.write("<td><a href='%s/%s'>%.2f</a></td>" % (spoofDir, segment + '-' + name + '.T04-BDS-RAIM.png',  mean(BDSUsd[i])))
          
        # 2D Stats
        htmlfile.write("<td><a href='%s/%s'>%.2f</a></td>" % 
            (spoofDir, segment + '-' + filename + '.2DCDF.png',cx_2d[i_50_2d]))
        htmlfile.write("<td>%.2f</td><td>%.2f</td><td>%.2f</td><td>%.2f</td>" % 
            (cx_2d[i_68_2d],cx_2d[i_90_2d],cx_2d[i_95_2d],cx_2d[-1]))

        # Hgt stats
        htmlfile.write("<td><a href='%s/%s'>%.2f</a></td>" % 
            (spoofDir, segment + '-' + filename + '.HgtCDF.png',cx_hgt[i_50_hgt]))
        htmlfile.write("<td>%.2f</td><td>%.2f</td><td>%.2f</td><td>%.2f</td></tr>\n" % 
            (cx_hgt[i_68_hgt],cx_hgt[i_90_hgt],cx_hgt[i_95_hgt],cx_hgt[-1]))

        if(args.record35_16):
          os.system("t012kml --ignore_err -o%s-%s --sub16 -is%d,%.0f -es%d,%.0f %s %s/%s-%s.kmz" % (segment,filename,week,start,week,stop,datafile,spoofDir,segment,filename))
        else:
          os.system("t012kml --ignore_err -o%s-%s -is%d,%.0f -es%d,%.0f %s %s/%s-%s.kmz" % (segment,filename,week,start,week,stop,datafile,spoofDir,segment,filename))

  # Close the stats file
  statsfile.close()
  htmlfile.write("</tbody>\n")
  htmlfile.write("</table>\n")
  htmlfile.write("<p>*HMI 2D - Hazadous Misleading Information - percentage of time that the estimated 2D error is less than 20m (alarm limit) *and* the 2D error is more than 4.42 times the reported 2D estimated sigma</p>\n")
  htmlfile.write("<p>**Alarm - Percentage of time the estimated error is greater than the alarm threshold of 20m</p>\n")
  htmlfile.write("</body>\n")
  htmlfile.write("</html>\n")
  htmlfile.close()

  # Remove the temporary file we created
  os.remove(datafile)