import datetime
import numpy as np

__authors__ = ["Ming Chen"]
__license__ = "Public Domain"
__maintainer__ = "Ming Chen"
__email__ = "ming.chen@noaa.gov"
__date__  = "March 2, 2026"


""" 
#####################################################################
# 
# leosp3.py
#
# Purpose:
# This python script is a lib with the function  
#        1) read BERNESE LEO orbit file and clock file, 
#        2) write the LEO orbit & clock in SP3
#
# History:  
#       Created by: Ming Chen 2026/03/02
#                   ming.chen@noaa.gov
# 
# Usage: python leosp3.py orbit_file clock_file sp3_file
#
#
#####################################################################

"""
# ------------------------------------------------------------
# Parse the clock file (RINEX CLK-like format)
# ------------------------------------------------------------

def parse_clock_file(clock_file):
    """
    Reads a clock file and returns a dictionary:
        clock[(sat, datetime)] = clock_value
    """
    clock_dict = {}

    with open(clock_file, "r") as f:
        for line in f:
            if not line.startswith("AR"):
                continue

            parts = line.split()

            # Satellite ID (e.g., G29)
            sat = parts[1].strip()

            # Date-time fields
            year  = int(parts[2])
            month = int(parts[3])
            day   = int(parts[4])
            hour  = int(parts[5])
            minute = int(parts[6])
            sec    = float(parts[7])

            # Build datetime object
            t = datetime.datetime(
                year, month, day, hour, minute,
                int(sec),
                int((sec % 1) * 1e6)
            )

            # Clock offset (seconds)
            clk = float(parts[9])
            clock_dict[t] = clk*1e6

    return sat, clock_dict



def leo_sp3(orbit_file, clock_file, sp3_file):
    """
    Reads an orbit file and rewrites it with updated clock values.
    Only "P" lines are modified.
    All other lines (header, V-lines, comments) remain unchanged.
    """
    from scipy.interpolate import interp1d, CubicSpline    
    
    # Load clock values
    sat, ts = parse_clock_file(clock_file)

    current_time = None
    output_lines = []
    leosats = {'L23': 'C2E1', 'L24': 'C2E2','L25': 'C2E3',\
               'L26': 'C2E4', 'L27': 'C2E5','L28': 'C2E6'}
    tref =  datetime.datetime(2000,1,1)
 
    epochs=[]
    with open(orbit_file, "r") as f:
        for line in f:

            # Epoch line: "* yyyy mm dd hh mm ss"
            if line.startswith("*"):
                parts = line.split()
                year, month, day = map(int, parts[1:4])
                hour, minute = map(int, parts[4:6])
                sec = float(parts[6])

                current_time = datetime.datetime(
                    year, month, day, hour, minute,
                    int(sec),
                    int((sec % 1) * 1e6)
                )

                epochs.append(current_time)
                continue
    
    xy = np.array([((t - tref).total_seconds(), ts[t]) for t in sorted(ts.keys())])
    x = xy[:,0] ; y = xy[:,1]
    tnew = np.array([(t - tref).total_seconds() for t in epochs])
   
    cs = CubicSpline(x,y)
    csclk =  cs(tnew)
    clk_new = dict(zip(epochs, csclk))

    with open(orbit_file, "r") as f:
        for line in f:

            # Epoch line: "* yyyy mm dd hh mm ss"
            if line.startswith("*"):
                parts = line.split()
                year, month, day = map(int, parts[1:4])
                hour, minute = map(int, parts[4:6])
                sec = float(parts[6])

                current_time = datetime.datetime(
                    year, month, day, hour, minute,
                    int(sec),
                    int((sec % 1) * 1e6)
                )

                output_lines.append(line)
                continue

            # Position line
            if line.startswith("P") or line.startswith("V") :
                sat = line[1:5].strip()   # e.g., "L23"

                prefix = line[:5]         # "PL23 "
                parts = line[5:].split()

                x, y, z = map(float, parts[:3])
                old_clock = float(parts[3])

                newline = (
                    f"{prefix}"
                    f"{x:13.6f}{y:14.6f}{z:14.6f}{old_clock:16.6f}\n"
                   )
                if line.startswith("P"): 
                   # Rebuild P-line with the same numeric formatting
                   # Lookup new clock (use exactly the same sat ID)
                   clk = clk_new[current_time]
                   newline = (
                      f"{prefix}"
                      f"{x:13.6f}{y:14.6f}{z:14.6f}{clk:16.6f}\n"
                     )
                     #f"{x:13.6f}{y:14.6f}{z:14.6f}{clk:16.12f}\n"

                output_lines.append(newline)
                continue

            # All other lines copied unchanged
            output_lines.append(line)

    # Write final output file
    with open(sp3_file, "w") as f:
        f.writelines(output_lines)

if __name__ == "__main__":
   import argparse
   import os, glob 
   parser = argparse.ArgumentParser(description='Main program')
   parser.add_argument("orbit_file", type=str,  help="leo orbit file")
   parser.add_argument("clock_file", type=str,  help="leo clock file")
   parser.add_argument("sp3_file", type=str,  help="outpu leo sp3")
   parser.add_argument("-v","--verbose", action='store_true',  help="print info for debug")
   args = parser.parse_args()
   if os.path.isfile(args.orbit_file) and \
           os.path.isfile(args.clock_file):
      leo_sp3(
        orbit_file = args.orbit_file,
        clock_file = args.clock_file,
        sp3_file   = args.sp3_file
      )
   else:
      print(f'Error: \n orbit_files not exist ...\n {args.orbit_file} or \n {args.clock_file}')

