# For table of nav data indices from GPS ICD: https://docs.google.com/spreadsheets/d/15VBzAbA_MOXZ_QZ6Ql8pBVfHuuysv2viKApt_ybzXv4/edit?usp=sharing 

import re
import glob
import matplotlib.pyplot as plt
import numpy as np
from numpy.lib import recfunctions as rfn

from mutils import *
import allantools
from bitarray import bitarray

import argparse
import os
import compare_nav_msg_data_bits as cmp_dbits


class Ephemeris:
    class Subframe1:
        def __init__(self, wn=None, tow=None, tgd=None, toc=None, af2=None, af1=None, af0=None, iodc=None):
            self.wn = wn
            self.tow = tow
            self.tgd = tgd
            self.toc = toc
            self.af2 = af2
            self.af1 = af1
            self.af0 = af0
            self.iodc = iodc

    class Subframe2:
        def __init__(self, iode=None, crs=None, delta_n=None, m0=None, cuc=None, ecc=None, cus=None, sqrt_a=None, toe=None, tow=None):
            self.tow = tow
            self.iode = iode
            self.crs = crs
            self.delta_n = delta_n
            self.m0 = m0
            self.cuc = cuc
            self.ecc = ecc
            self.cus = cus
            self.sqrt_a = sqrt_a
            self.toe = toe

    class Subframe3:
        def __init__(self, iode=None, cic=None, omega0=None, cis=None, i0=None, crc=None, omega=None, omega_dot=None, idot=None, tow=None):
            self. tow = tow
            self.iode = iode
            self.cic = cic
            self.omega0 = omega0
            self.cis = cis
            self.i0 = i0
            self.crc = crc
            self.omega = omega
            self.omega_dot = omega_dot
            self.idot = idot
    
    class Almanac:
        def __init__(self, toa=None, ecc=None, delta_i=None, omega_dot=None, sqrt_a=None, omega_o=None, i0=None, m0=None, af0=None, af1=None, tow=None):
            self.tow = tow
            self.toa = toa
            self.ecc = ecc
            self.delta_i = delta_i
            self.omega_dot = omega_dot
            self.sqrt_a = sqrt_a
            self.omega_o = omega_o
            self.i0 = i0
            self.m0 = m0
            self.af0 = af0
            self.af1 = af1
            
    def __init__(self):
        self.sf1 = Ephemeris.Subframe1()
        self.sf2 = Ephemeris.Subframe2()
        self.sf3 = Ephemeris.Subframe3()

def twos_complement(v, nb):
    v = int(v, 2)
    sign = v >> (nb-1)
    value = v
    if sign != 0:
        value = value - (1<<nb)
    return value

def extract_nav_data(cls_gps_dbits, prn):
        arr = cls_gps_dbits
        tow_col = arr.k['TOW']
        prn_col = arr.k['PRN']
        flags_col = arr.k['FLAGS']

        sf1 = []
        sf2 = []
        sf3 = []
        almanac = []

        for i, row in enumerate(arr):
            # Account for the potential of an extra flag byte (check if last bit is set)
            if row[flags_col] & 0b00000010:
                N1 = 9
            else:
                N1 = 10

            # Convert to bitarray
            bits = bitarray()
            bits.frombytes(row[N1:48].astype(np.uint8).tobytes())
            if row[flags_col] & 0b11 != 0:
                continue  # Not L1CA

            subf = int(bits[49:52].to01(), 2)
            TOW = int(bits[30:47].to01(), 2) * 6

            # Subframe 1
            if subf == 1 and row[prn_col] == prn:
                tgd = twos_complement(bits[196:204].to01(), 8) * pow(2, -31)
                toc = int(bits[218:234].to01(), 2) * pow(2, 4)
                af2 = twos_complement(bits[240:248].to01(), 16) * pow(2, -55)
                af1 = twos_complement(bits[248:264].to01(), 16) * pow(2, -43)
                af0 = twos_complement(bits[270:292].to01(), 22) * pow(2, -31)
                iodc = int((bits[82:84] + bits[210:218]).to01(), 2)
                sf1.append(Ephemeris.Subframe1(
                    wn=int(bits[17:27].to01(), 2),
                    tow=TOW, tgd=tgd, toc=toc, af2=af2, af1=af1, af0=af0, iodc=iodc
                ))

            # Subframe 2
            if subf == 2 and row[prn_col] == prn:
                iode = int(bits[60:68].to01(), 2)
                crs = twos_complement(bits[68:84].to01(), 16) * pow(2, -5)
                delta_n = twos_complement(bits[90:106].to01(), 16) * pow(2, -43)
                m0 = twos_complement((bits[106:114] + bits[120:144]).to01(), 32) * pow(2, -31)
                cuc = twos_complement(bits[150:166].to01(), 16) * pow(2, -29)
                cus = twos_complement(bits[210:226].to01(), 16) * pow(2, -29)
                ecc = int((bits[166:174] + bits[180:204]).to01(), 2) * pow(2, -33)
                sqrt_a = int((bits[226:234] + bits[240:264]).to01(), 2) * pow(2, -19)
                toe = int(bits[270:286].to01(), 2) * pow(2, 4)
                sf2.append(Ephemeris.Subframe2(
                    iode=iode, crs=crs, delta_n=delta_n, m0=m0, cuc=cuc,
                    ecc=ecc, cus=cus, sqrt_a=sqrt_a, toe=toe, tow=TOW
                ))

            # Subframe 3
            if subf == 3 and row[prn_col] == prn:
                cic = twos_complement(bits[60:76].to01(), 16) * pow(2, -29)
                omega0 = twos_complement((bits[76:84] + bits[90:114]).to01(), 32) * pow(2, -31)
                cis = twos_complement(bits[120:136].to01(), 16) * pow(2, -29)
                i0 = twos_complement((bits[136:144] + bits[150:174]).to01(), 32) * pow(2, -31)
                crc = twos_complement(bits[180:196].to01(), 16) * pow(2, -5)
                omega = twos_complement((bits[196:204] + bits[210:234]).to01(), 32) * pow(2, -31)
                omega_dot = twos_complement(bits[240:264].to01(), 24) * pow(2, -43)
                iode = int(bits[271:278].to01(), 2)
                idot = twos_complement(bits[278:292].to01(), 14) * pow(2, -43)
                sf3.append(Ephemeris.Subframe3(
                    iode=iode, cic=cic, omega0=omega0, cis=cis, i0=i0,
                    crc=crc, omega=omega, omega_dot=omega_dot, idot=idot, tow=TOW
                ))

            # Almanac
            if (subf == 4 or subf == 5) and row[prn_col] == prn:
                page = int(bits[63:68].to01(), 2)
                if subf == 5 and page > 24:
                    continue
                if page != prn:
                    continue
                ecc = int(bits[68:84].to01(), 2) * pow(2, -21)
                toa = int(bits[90:98].to01(), 2) * pow(2, 12)
                delta_i = twos_complement(bits[98:114].to01(), 16) * pow(2, -19)
                omega_dot = twos_complement(bits[120:136].to01(), 16) * pow(2, -38)
                sqrt_a = int(bits[150:174].to01(), 2) * pow(2, -11)
                omega_o = twos_complement(bits[180:204].to01(), 24) * pow(2, -23)
                i0 = twos_complement(bits[210:234].to01(), 24) * pow(2, -23)
                m0 = twos_complement(bits[240:264].to01(), 24) * pow(2, -23)
                af0 = twos_complement((bits[270:278] + bits[290:292]).to01(), 11) * pow(2, -20)
                af1 = twos_complement(bits[278:289].to01(), 11) * pow(2, -38)
                almanac.append(Ephemeris.Almanac(
                    ecc=ecc, toa=toa, delta_i=delta_i, omega_dot=omega_dot,
                    sqrt_a=sqrt_a, omega_o=omega_o, i0=i0, m0=m0, af0=af0, af1=af1, tow=TOW
                ))

        return sf1, sf2, sf3, almanac

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description="Process file paths for GPS data bits.")
    parser.add_argument('filepath', type=str, help='Path to the main test file')
    parser.add_argument('prn', type=int, help='PRN number to analyze')
    parser.add_argument('-filepath_ref', type=str, nargs='*', default=[], help='Paths to reference files')
    parser.add_argument('-sat_type', type=str, default='GPS', help='Satellite type')
    parser.add_argument('-s', type=str, default='', help='Start time for analysis (optional)')
    parser.add_argument('-e', type=str, default='', help='End time for analysis (optional)')


    args = parser.parse_args()

    args.sat_type = args.sat_type.upper()
    args.prn = int(args.prn)
    args.s = int(args.s) if args.s else None
    args.e = int(args.e) if args.e else None

    filepath = args.filepath

    # Load the databits in test file
    cls_gps_dbits = cmp_dbits.get_gps_dbits(filepath,args.sat_type)

    ref = False
    if len(args.filepath_ref) > 0:
        ref = True
        filepath_ref_list = args.filepath_ref

        combined_dbits = cmp_dbits.generate_dbits_reference(filepath_ref_list, "gps_dbits_reference.txt",args.sat_type)

        cls_gps_dbits_ref = cmp_dbits.dbits_arr_to_cls(combined_dbits)
        sf1_ref, sf2_ref, sf3_ref, almanac_ref = extract_nav_data(cls_gps_dbits_ref, args.prn)

    sf1, sf2, sf3, almanac = extract_nav_data(cls_gps_dbits, args.prn)
    
  
    if len(sf1) > 0:         
        fig, axs = plt.subplots(3, 3, figsize=(16, 12))
        fig.suptitle("GPS Subframe 1 Characteristics")

        if args.s is not None and args.e is not None:
            sf1_arr = np.array([s for s in sf1 if s.tow is not None and args.s <= s.tow <= args.e])
        else:
            sf1_arr = np.array(sf1)

        wn = [s.wn for s in sf1_arr]
        tow = [s.tow for s in sf1_arr]
        tgd = [s.tgd for s in sf1_arr]
        toc = [s.toc for s in sf1_arr]
        af2 = [s.af2 for s in sf1_arr]
        af1 = [s.af1 for s in sf1_arr]
        af0 = [s.af0 for s in sf1_arr]
        iodc = [s.iodc for s in sf1_arr]

        if ref:
            # Filter reference data by TOW range
            # Filter reference data by TOW range using list comprehension
            if args.s is not None and args.e is not None:
                sf1_ref_arr = [s for s in sf1_ref if s.tow is not None and args.s <= s.tow <= args.e]
            else:
                sf1_ref_arr = sf1_ref
            wn_ref = [s.wn for s in sf1_ref_arr]
            tow_ref = [s.tow for s in sf1_ref_arr]
            tgd_ref = [s.tgd for s in sf1_ref_arr]
            toc_ref = [s.toc for s in sf1_ref_arr]
            af2_ref = [s.af2 for s in sf1_ref_arr]
            af1_ref = [s.af1 for s in sf1_ref_arr]
            af0_ref = [s.af0 for s in sf1_ref_arr]
            iodc_ref = [s.iodc for s in sf1_ref_arr]
           
        
        axs[0, 0].plot(tow, wn, 'x', label="WN")
        if ref:
            axs[0, 0].plot(tow_ref, wn_ref, 'x', label="WN Ref")
            #axs[0, 0].legend(['Test Data', 'Reference Data'])

        axs[0, 0].set_xlabel("TOW")
        axs[0, 0].set_ylabel("WN")
        axs[0, 0].tick_params(axis='x', labelrotation=45)

        axs[0, 1].plot(tow, tgd, 'x', label="TGD")
        if ref:
            axs[0, 1].plot(tow_ref, tgd_ref, 'x',label ="TGD Ref")
        axs[0, 1].set_title("TGD")
        axs[0, 1].set_xlabel("TOW")
        axs[0, 1].set_ylabel("TGD")
        axs[0, 1].tick_params(axis='x', labelrotation=45)

        axs[0, 2].plot(tow, toc, 'x', label="TOC")
        if ref:
            axs[0, 2].plot(tow_ref, toc_ref, 'x', label="TOC Ref")
        axs[0, 2].set_title("TOC")
        axs[0, 2].set_xlabel("TOW")
        axs[0, 2].set_ylabel("TOC")
        axs[0, 2].tick_params(axis='x', labelrotation=45)

        axs[1, 0].plot(tow, af2, 'x', label="af2")
        if ref:
            axs[1, 0].plot(tow_ref, af2_ref, 'x', label="af2_ref")
        axs[1, 0].set_title("af2")
        axs[1, 0].set_xlabel("TOW")
        axs[1, 0].set_ylabel("af2")
        axs[1, 0].legend()
        axs[1, 0].tick_params(axis='x', labelrotation=45)

        axs[1, 1].plot(tow, af1, 'x', label="af1")
        if ref:
            axs[1, 1].plot(tow_ref, af1_ref, 'x', label="af1_ref")
        axs[1, 1].set_title("af1")
        axs[1, 1].set_xlabel("TOW")
        axs[1, 1].set_ylabel("af1")
        axs[1, 1].legend()
        axs[1, 1].tick_params(axis='x', labelrotation=45)

        axs[1, 2].plot(tow, af0, 'x', label="af0")
        if ref:
            axs[1, 2].plot(tow_ref, af0_ref, 'x', label="af0_ref")
        axs[1, 2].set_title("af0")
        axs[1, 2].set_xlabel("TOW")
        axs[1, 2].set_ylabel("af0")
        axs[1, 2].legend()
        axs[1, 2].tick_params(axis='x', labelrotation=45)

        axs[2, 0].plot(tow, iodc, 'x', label="iodc")
        if ref:
            axs[2, 0].plot(tow_ref, iodc_ref, 'x', label="iodc_ref")
        axs[2, 0].set_title("IODC")
        axs[2, 0].set_xlabel("TOW")
        axs[2, 0].set_ylabel("IODC")
        axs[2, 0].legend()
        axs[2, 0].tick_params(axis='x', labelrotation=45)

        # Hide unused subplots
        axs[2, 1].axis('off')
        axs[2, 2].axis('off')

        # Add PRN group title
        fig.suptitle(f"GPS Subframe 1 Nav Data (PRN {args.prn})", fontsize=16)
        test_filename = os.path.basename(filepath)
        ref_filename = os.path.basename(args.filepath_ref[0]) if ref and len(args.filepath_ref) > 0 else "None"
        fig.suptitle(
            f"GPS Subframe 1 Nav Data (PRN {args.prn})\nTest: {test_filename} | Ref: {ref_filename}",
            fontsize=16
        )
        plt.tight_layout(rect=[0, 0.03, 1, 0.95])
        plt.savefig("subframe1_characteristics.png")
        plt.close(fig)
    else:
        print("No subframe 1 data found.")
    
    if len(sf2) > 0:
        
        if args.s is not None and args.e is not None:
            sf2_arr = [s for s in sf2 if s.tow is not None and args.s <= s.tow <= args.e]
        else:
            sf2_arr = np.array(sf2)
        
        
        fig, axs = plt.subplots(3, 3, figsize=(16, 12))
        fig.suptitle("GPS Subframe 2 Characteristics")

        iode = [s.iode for s in sf2_arr]
        crs = [s.crs for s in sf2_arr]
        delta_n = [s.delta_n for s in sf2_arr]
        m0 = [s.m0 for s in sf2_arr]
        cuc = [s.cuc for s in sf2_arr]
        ecc = [s.ecc for s in sf2_arr]
        cus = [s.cus for s in sf2_arr]
        sqrt_a = [s.sqrt_a for s in sf2_arr]
        toe = [s.toe for s in sf2_arr]
        tow = [s.tow for s in sf2_arr]

        if ref:
            if args.s is not None and args.e is not None:
                sf2_ref_arr = [s for s in sf2_ref if s.tow is not None and args.s <= s.tow <= args.e]
            else:
                sf2_ref_arr = sf2_ref

            iode_ref = [s.iode for s in sf2_ref_arr]
            crs_ref = [s.crs for s in sf2_ref_arr]
            delta_n_ref = [s.delta_n for s in sf2_ref_arr]
            m0_ref = [s.m0 for s in sf2_ref_arr]
            cuc_ref = [s.cuc for s in sf2_ref_arr]
            ecc_ref = [s.ecc for s in sf2_ref_arr]
            cus_ref = [s.cus for s in sf2_ref_arr]
            sqrt_a_ref = [s.sqrt_a for s in sf2_ref_arr]
            toe_ref = [s.toe for s in sf2_ref_arr]
            tow_ref = [s.tow for s in sf2_ref_arr]

        axs[0, 0].plot(tow, iode, 'x', label="IODE")
        if ref:
            axs[0, 0].plot(tow_ref, iode_ref, 'x', label="IODE Ref")
        axs[0, 0].set_xlabel("TOW")
        axs[0, 0].set_ylabel("IODE")
        axs[0, 0].set_title("IODE")
        axs[0, 0].tick_params(axis='x', labelrotation=45)

        axs[0, 1].plot(tow, crs, 'x', label="CRS")
        if ref:
            axs[0, 1].plot(tow_ref, crs_ref, 'x', label="CRS Ref")
        axs[0, 1].set_xlabel("TOW")
        axs[0, 1].set_ylabel("CRS")
        axs[0, 1].set_title("CRS")
        axs[0, 1].tick_params(axis='x', labelrotation=45)

        axs[0, 2].plot(tow, delta_n, 'x', label="Delta N")
        if ref:
            axs[0, 2].plot(tow_ref, delta_n_ref, 'x', label="Delta N Ref")
        axs[0, 2].set_xlabel("TOW")
        axs[0, 2].set_ylabel("Delta N")
        axs[0, 2].set_title("Delta N")
        axs[0, 2].tick_params(axis='x', labelrotation=45)

        axs[1, 0].plot(tow, m0, 'x', label="M0")
        if ref:
            axs[1, 0].plot(tow_ref, m0_ref, 'x', label="M0 Ref")
        axs[1, 0].set_xlabel("TOW")
        axs[1, 0].set_ylabel("M0")
        axs[1, 0].set_title("M0")
        axs[1, 0].tick_params(axis='x', labelrotation=45)

        axs[1, 1].plot(tow, cuc, 'x', label="CUC")
        if ref:
            axs[1, 1].plot(tow_ref, cuc_ref, 'x', label="CUC Ref")
        axs[1, 1].set_xlabel("TOW")
        axs[1, 1].set_ylabel("CUC")
        axs[1, 1].set_title("CUC")
        axs[1, 1].tick_params(axis='x', labelrotation=45)

        axs[1, 2].plot(tow, ecc, 'x', label="Ecc")
        if ref:
            axs[1, 2].plot(tow_ref, ecc_ref, 'x', label="Ecc Ref")
        axs[1, 2].set_xlabel("TOW")
        axs[1, 2].set_ylabel("Eccentricity")
        axs[1, 2].set_title("Eccentricity")
        axs[1, 2].tick_params(axis='x', labelrotation=45)

        axs[2, 0].plot(tow, cus, 'x', label="CUS")
        if ref:
            axs[2, 0].plot(tow_ref, cus_ref, 'x', label="CUS Ref")
        axs[2, 0].set_xlabel("TOW")
        axs[2, 0].set_ylabel("CUS")
        axs[2, 0].set_title("CUS")
        axs[2, 0].tick_params(axis='x', labelrotation=45)

        axs[2, 1].plot(tow, sqrt_a, 'x', label="SQRT(A)")
        if ref:
            axs[2, 1].plot(tow_ref, sqrt_a_ref, 'x', label="SQRT(A) Ref")
        axs[2, 1].set_xlabel("TOW")
        axs[2, 1].set_ylabel("SQRT(A)")
        axs[2, 1].set_title("SQRT(A)")
        axs[2, 1].tick_params(axis='x', labelrotation=45)

        axs[2, 2].plot(tow, toe, 'x', label="TOE")
        if ref:
            axs[2, 2].plot(tow_ref, toe_ref, 'x', label="TOE Ref")
        axs[2, 2].set_xlabel("TOW")
        axs[2, 2].set_ylabel("TOE")
        axs[2, 2].set_title("TOE")
        axs[2, 2].tick_params(axis='x', labelrotation=45)

        fig.suptitle(f"GPS Subframe 2 Nav Data (PRN {args.prn})", fontsize=16)
        test_filename = os.path.basename(filepath)
        ref_filename = os.path.basename(args.filepath_ref[0]) if ref and len(args.filepath_ref) > 0 else "None"
        fig.suptitle(
            f"GPS Subframe 2 Nav Data (PRN {args.prn})\nTest: {test_filename} | Ref: {ref_filename}",
            fontsize=16
        )
        plt.tight_layout(rect=[0, 0.03, 1, 0.95])
        plt.savefig("subframe2_characteristics.png")
        plt.close(fig)
    
    if len(sf3) > 0:
        
        if args.s is not None and args.e is not None:
            sf3_arr = [s for s in sf3 if s.tow is not None and args.s <= s.tow <= args.e]
        else:
            sf3_arr = np.array(sf3)

        fig, axs = plt.subplots(3, 3, figsize=(16, 12))
        fig.suptitle("GPS Subframe 3 Characteristics")

        tow = [s.tow for s in sf3_arr]
        iode = [s.iode for s in sf3_arr]
        cic = [s.cic for s in sf3_arr]
        omega0 = [s.omega0 for s in sf3_arr]
        cis = [s.cis for s in sf3_arr]
        i0 = [s.i0 for s in sf3_arr]
        crc = [s.crc for s in sf3_arr]
        omega = [s.omega for s in sf3_arr]
        omega_dot = [s.omega_dot for s in sf3_arr]
        idot = [s.idot for s in sf3_arr]

        if ref:
            if args.s is not None and args.e is not None:
                sf3_ref_arr = [s for s in sf3_ref if s.tow is not None and args.s <= s.tow <= args.e]
            else:
                sf3_ref_arr = sf3_ref

            tow_ref = [s.tow for s in sf3_ref_arr]
            iode_ref = [s.iode for s in sf3_ref_arr]
            cic_ref = [s.cic for s in sf3_ref_arr]
            omega0_ref = [s.omega0 for s in sf3_ref_arr]
            cis_ref = [s.cis for s in sf3_ref_arr]
            i0_ref = [s.i0 for s in sf3_ref_arr]
            crc_ref = [s.crc for s in sf3_ref_arr]
            omega_ref = [s.omega for s in sf3_ref_arr]
            omega_dot_ref = [s.omega_dot for s in sf3_ref_arr]
            idot_ref = [s.idot for s in sf3_ref_arr]

        axs[0, 0].plot(tow, iode, 'x', label="IODE")
        if ref:
            axs[0, 0].plot(tow_ref, iode_ref, 'x', label="IODE Ref")
        axs[0, 0].set_xlabel("TOW")
        axs[0, 0].set_ylabel("IODE")
        axs[0, 0].set_title("IODE")
        axs[0, 0].tick_params(axis='x', labelrotation=45)

        axs[0, 1].plot(tow,cic, 'x', label="CIC")
        if ref:
            axs[0, 1].plot(tow_ref, cic_ref, 'x', label="CIC Ref")
        axs[0, 1].set_xlabel("TOW")
        axs[0, 1].set_ylabel("CIC")
        axs[0, 1].set_title("CIC")
        axs[0, 1].tick_params(axis='x', labelrotation=45)
        axs[0, 2].plot(tow, omega0, 'x', label="Omega0")
        if ref:
            axs[0, 2].plot(tow_ref, omega0_ref, 'x', label="Omega0 Ref")
        axs[0, 2].set_xlabel("TOW")
        axs[0, 2].set_ylabel("Omega0")
        axs[0, 2].set_title("Omega0")
        axs[0, 2].tick_params(axis='x', labelrotation=45)

        axs[1, 0].plot(tow, cis, 'x', label="CIS")
        if ref:
            axs[1, 0].plot(tow_ref, cis_ref, 'x', label="CIS Ref")
        axs[1, 0].set_xlabel("TOW")
        axs[1, 0].set_ylabel("CIS")
        axs[1, 0].set_title("CIS")
        axs[1, 0].tick_params(axis='x', labelrotation=45)

        axs[1, 1].plot(tow, i0, 'x', label="I0")
        if ref:
            axs[1, 1].plot(tow_ref, i0_ref, 'x', label="I0 Ref")
        axs[1, 1].set_xlabel("TOW")
        axs[1, 1].set_ylabel("I0")
        axs[1, 1].set_title("I0")
        axs[1, 1].tick_params(axis='x', labelrotation=45)

        axs[1, 2].plot(tow, crc, 'x', label="CRC")
        if ref:
            axs[1, 2].plot(tow_ref, crc_ref, 'x', label="CRC Ref")
        axs[1, 2].set_xlabel("TOW")
        axs[1, 2].set_ylabel("CRC")
        axs[1, 2].set_title("CRC")
        axs[1, 2].tick_params(axis='x', labelrotation=45)

        axs[2, 0].plot(tow, omega, 'x', label="Omega")
        if ref:
            axs[2, 0].plot(tow_ref, omega_ref, 'x', label="Omega Ref")
        axs[2, 0].set_xlabel("TOW")
        axs[2, 0].set_ylabel("Omega")
        axs[2, 0].set_title("Omega")
        axs[2, 0].tick_params(axis='x', labelrotation=45)

        axs[2, 1].plot(tow, omega_dot, 'x', label="Omega Dot")
        if ref:
            axs[2, 1].plot(tow_ref, omega_dot_ref, 'x', label="Omega Dot Ref")
        axs[2, 1].set_xlabel("TOW")
        axs[2, 1].set_ylabel("Omega Dot")
        axs[2, 1].set_title("Omega Dot")
        axs[2, 1].tick_params(axis='x', labelrotation=45)

        axs[2, 2].plot(tow, idot, 'x', label="IDOT")
        if ref:
            axs[2, 2].plot(tow_ref, idot_ref, 'x', label="IDOT Ref")
        axs[2, 2].set_xlabel("TOW")
        axs[2, 2].set_ylabel("IDOT")
        axs[2, 2].set_title("IDOT")
        axs[2, 2].tick_params(axis='x', labelrotation=45)

        fig.suptitle(f"GPS Subframe 3 Nav Data (PRN {args.prn})", fontsize=16)
        test_filename = os.path.basename(filepath)
        ref_filename = os.path.basename(args.filepath_ref[0]) if ref and len(args.filepath_ref) > 0 else "None"
        fig.suptitle(
            f"GPS Subframe 3 Nav Data (PRN {args.prn})\nTest: {test_filename} | Ref: {ref_filename}",
            fontsize=16
        )
        plt.tight_layout(rect=[0, 0.03, 1, 0.95])
        plt.savefig("subframe3_characteristics.png")
        plt.close(fig)
    
    if len(almanac) > 0:

        if args.s is not None and args.e is not None:
            almanac_arr = [s for s in almanac if s.tow is not None and args.s <= s.tow <= args.e]
        else:
            almanac_arr = np.array(almanac)

        fig, axs = plt.subplots(3, 4, figsize=(18, 12))
        fig.suptitle("GPS Almanac Characteristics")

        tow = [s.tow for s in almanac_arr]
        ecc = [s.ecc for s in almanac_arr]
        toa = [s.toa for s in almanac_arr]
        delta_i = [s.delta_i for s in almanac_arr]
        omega_dot = [s.omega_dot for s in almanac_arr]
        sqrt_a = [s.sqrt_a for s in almanac_arr]
        omega_o = [s.omega_o for s in almanac_arr]
        i0 = [s.i0 for s in almanac_arr]
        m0 = [s.m0 for s in almanac_arr]
        af0 = [s.af0 for s in almanac_arr]
        af1 = [s.af1 for s in almanac_arr]

        if ref:
            if args.s is not None and args.e is not None:
                almanac_ref_arr = [s for s in almanac_ref if s.tow is not None and args.s <= s.tow <= args.e]
            else:
                almanac_ref_arr = almanac_ref

            tow_ref = [s.tow for s in almanac_ref_arr]
            ecc_ref = [s.ecc for s in almanac_ref_arr]
            toa_ref = [s.toa for s in almanac_ref_arr]
            delta_i_ref = [s.delta_i for s in almanac_ref_arr]
            omega_dot_ref = [s.omega_dot for s in almanac_ref_arr]
            sqrt_a_ref = [s.sqrt_a for s in almanac_ref_arr]
            omega_o_ref = [s.omega_o for s in almanac_ref_arr]
            i0_ref = [s.i0 for s in almanac_ref_arr]
            m0_ref = [s.m0 for s in almanac_ref_arr]
            af0_ref = [s.af0 for s in almanac_ref_arr]
            af1_ref = [s.af1 for s in almanac_ref_arr]

        axs[0, 0].plot(tow, ecc, 'x', label="Eccentricity")
        if ref:
            axs[0, 0].plot(tow_ref, ecc_ref, 'x', label="Eccentricity Ref")
        axs[0, 0].set_xlabel("TOW")
        axs[0, 0].set_ylabel("Eccentricity")
        axs[0, 0].set_title("Eccentricity")
        axs[0, 0].tick_params(axis='x', labelrotation=45)

        axs[0, 1].plot(tow, toa, 'x', label="TOA")
        if ref:
            axs[0, 1].plot(tow_ref, toa_ref, 'x', label="TOA Ref")
        axs[0, 1].set_xlabel("TOW")
        axs[0, 1].set_ylabel("TOA")
        axs[0, 1].set_title("TOA")
        axs[0, 1].tick_params(axis='x', labelrotation=45)

        axs[0, 2].plot(tow, delta_i, 'x', label="Delta I")
        if ref:
            axs[0, 2].plot(tow_ref, delta_i_ref, 'x', label="Delta I Ref")
        axs[0, 2].set_xlabel("TOW")
        axs[0, 2].set_ylabel("Delta I")
        axs[0, 2].set_title("Delta I")
        axs[0, 2].tick_params(axis='x', labelrotation=45)

        axs[0, 3].plot(tow, omega_dot, 'x', label="Omega Dot")
        if ref:
            axs[0, 3].plot(tow_ref, omega_dot_ref, 'x', label="Omega Dot Ref")
        axs[0, 3].set_xlabel("TOW")
        axs[0, 3].set_ylabel("Omega Dot")
        axs[0, 3].set_title("Omega Dot")
        axs[0, 3].tick_params(axis='x', labelrotation=45)

        axs[1, 0].plot(tow, sqrt_a, 'x', label="SQRT(A)")
        if ref:
            axs[1, 0].plot(tow_ref, sqrt_a_ref, 'x', label="SQRT(A) Ref")
        axs[1, 0].set_xlabel("TOW")
        axs[1, 0].set_ylabel("SQRT(A)")
        axs[1, 0].set_title("SQRT(A)")
        axs[1, 0].tick_params(axis='x', labelrotation=45)

        axs[1, 1].plot(tow, omega_o, 'x', label="Omega_o")
        if ref:
            axs[1, 1].plot(tow_ref, omega_o_ref, 'x', label="Omega_o Ref")
        axs[1, 1].set_xlabel("TOW")
        axs[1, 1].set_ylabel("Omega_o")
        axs[1, 1].set_title("Omega_o")
        axs[1, 1].tick_params(axis='x', labelrotation=45)

        axs[1, 2].plot(tow, i0, 'x', label="I0")
        if ref:
            axs[1, 2].plot(tow_ref, i0_ref, 'x', label="I0 Ref")
        axs[1, 2].set_xlabel("TOW")
        axs[1, 2].set_ylabel("I0")
        axs[1, 2].set_title("I0")
        axs[1, 2].tick_params(axis='x', labelrotation=45)

        axs[1, 3].plot(tow, m0, 'x', label="M0")
        if ref:
            axs[1, 3].plot(tow_ref, m0_ref, 'x', label="M0 Ref")
        axs[1, 3].set_xlabel("TOW")
        axs[1, 3].set_ylabel("M0")
        axs[1, 3].set_title("M0")
        axs[1, 3].tick_params(axis='x', labelrotation=45)

        axs[2, 0].plot(tow, af0, 'x', label="af0")
        if ref:
            axs[2, 0].plot(tow_ref, af0_ref, 'x', label="af0 Ref")
        axs[2, 0].set_xlabel("TOW")
        axs[2, 0].set_ylabel("af0")
        axs[2, 0].set_title("af0")
        axs[2, 0].legend()
        axs[2, 0].tick_params(axis='x', labelrotation=45)

        axs[2, 1].plot(tow, af1, 'x', label="af1")
        if ref:
            axs[2, 1].plot(tow_ref, af1_ref, 'x', label="af1 Ref")
        axs[2, 1].set_xlabel("TOW")
        axs[2, 1].set_ylabel("af1")
        axs[2, 1].set_title("af1")
        axs[2, 1].legend()
        axs[2, 1].tick_params(axis='x', labelrotation=45)

        # Hide unused subplots
        axs[2, 2].axis('off')
        axs[2, 3].axis('off')

        fig.suptitle(f"GPS Almanac Nav Data (PRN {args.prn})", fontsize=16)
        test_filename = os.path.basename(filepath)
        ref_filename = os.path.basename(args.filepath_ref[0]) if ref and len(args.filepath_ref) > 0 else "None"
        fig.suptitle(
            f"GPS Almanac Nav Data (PRN {args.prn})\nTest: {test_filename} | Ref: {ref_filename}",
            fontsize=16
        )
        plt.tight_layout(rect=[0, 0.03, 1, 0.95])
        plt.savefig("almanac_characteristics.png")
        plt.close(fig)
    else:
        print("No almanac data found.")