#!/usr/bin/env python

from datetime import datetime,  timezone
import re, os


def custom_scientific_notation(x, precision=6):
    if x == 0:
        return f"{0:.{precision}f}E+00"

    import math

    sign = "-" if x < 0 else " "
    x_abs = abs(x)

    # Compute exponent so that coefficient < 1
    exponent = math.floor(math.log10(x_abs)) + 1
    coefficient = x_abs / (10 ** exponent)

    # Format the coefficient with the specified precision
    coefficient_str = f"{sign}{coefficient:.{precision}f}"

    return f"{coefficient_str}E{exponent:03d}"

# === Helper Functions ===
def parse_epoch_line(line):
    """Extract datetime from epoch header line starting with *"""
    parts = line.strip().split()
    year, month, day, hour, minute = map(int, parts[1:6])
    second = float(parts[6])
    return datetime(year, month, day, hour, minute, int(second)), second

def parse_sat_line(line):
    """Parse a satellite line starting with P and return PRN and clock bias"""
    prn = line[1:4]
    #clock_bias_str = line[-13:].strip()
    clock_bias_str = line.strip().split()[4]
    try:
        clock_bias = float(clock_bias_str)
    except ValueError:
        clock_bias = None  # Invalid entry like 999999.999999
    return prn, clock_bias


# === Main Parsing  ===

def read_records(input_file):
   epochs = []
   clock_records = {}  # {prn: [(datetime, sec_of_day, bias), ...]}

   with open(input_file, "r") as f:
      current_epoch = None
      current_sec = None
      for line in f:
        if line.startswith('*'):
            current_epoch, current_sec = parse_epoch_line(line)
            epochs.append(current_epoch)
        elif line.startswith('P'):
            prn, clock_bias = parse_sat_line(line)
            if clock_bias is not None:
                if prn not in clock_records:
                    clock_records[prn] = []
                clock_records[prn].append((current_epoch, current_sec, clock_bias*1.0e-6))
   return(clock_records)



def code_headers(dd, sats, n_sats):
   #dd = '20250629T00:00:00'
   d = datetime.strptime(dd,'%Y%m%dT%H:%M:%S')
   date_str = d.strftime("%d-%b-%y %H:%M").upper()
   prns_sorted = sorted(sats)
   code_header_lines = [
        "3.04                 C                    M                      RINEX VERSION / TYPE",
        f"CCRNXC V5.5          AIUB                 {date_str}        PGM / RUN BY / DATE ",
        "CODE RAPID PRECISE CLOCK INFORMATION                             COMMENT             ",
        "BASED ON CODE AND PHASE MEASUREMENTS                             COMMENT             ",
        "MERGED CLOCK INFORMATION FROM 5-MINUTE AND                       COMMENT             ",
        "30-SECOND CLOCK CORRECTIONS                                      COMMENT             ",
        "   GPS                                                           TIME SYSTEM ID      ",
        "    18                                                           LEAP SECONDS GNSS   ",
        "E  GPSEST V5.5        IGS20                                      SYS / PCVS APPLIED  ",
        "G  GPSEST V5.5        IGS20                                      SYS / PCVS APPLIED  ",
        "R  GPSEST V5.5        IGS20                                      SYS / PCVS APPLIED  ",
        "E  GPSEST V5.5        CODE.BIA @ ftp.aiub.unibe.ch/CODE/         SYS / DCBS APPLIED  ",
        "G  GPSEST V5.5        CODE.BIA @ ftp.aiub.unibe.ch/CODE/         SYS / DCBS APPLIED  ",
        "R  GPSEST V5.5        CODE.BIA @ ftp.aiub.unibe.ch/CODE/         SYS / DCBS APPLIED  ",
        "     2    AR    AS                                               # / TYPES OF DATA   ",
        "COD  Center for Orbit Determination in Europe                    ANALYSIS CENTER     ",
        "     1                                                           # OF CLK REF        ",
 
        "BREW00USA 40473M001                           0.000000000000E+00 ANALYSIS CLK REF    ",   
        "     1    IGS20                                                  # OF SOLN STA / TRF    ",
        "BREW00USA 40473M001           -2112007408 -3705351826  4726827014SOLN STA NAME / NUM ", 

       f"{n_sats:6d}                                                           # OF SOLN SATS"
    ]

   header_lines = code_header_lines
   # Add PRN list lines (exactly 15 IDs per line, pad with "00" if needed)
   prn_lines = []

   nid=16 ; padstr="  " ; arbias=0.0
   for i in range(0, len(prns_sorted), nid):
      chunk = prns_sorted[i:i+nid]
      # Pad chunk if less than nid
      if len(chunk) < nid:
        chunk += [padstr] * (nid - len(chunk))
      # Format line with right-aligned 3-char wide fields
      prn_line = ' '.join(f"{prn:>3}" for prn in chunk) + "  PRN LIST"
      prn_lines.append(prn_line)
   header_lines.extend(prn_lines)
   # Final COMMENT line
   header_lines.append("                                                                 END OF HEADER")
   return(header_lines)


# === Write Output File ===
def write_clk_codefmt(output_file, header_lines, clock_records):
  # Extract the directory path
  file_path = os.path.dirname(output_file)

  # Create the directory if it doesn't exist
  if file_path:  # Check if directory is not an empty string (e.g., if file_path is just a filename)
       os.makedirs(file_path, exist_ok=True)

  with open(output_file, "w") as out:
    # Write header
    for line in header_lines:
        out.write(line.rstrip() + "\n")

    # Write clock bias records
    prns_sorted = sorted(clock_records.keys())
    nrecords = len(clock_records[prns_sorted[0]])
    for i in range(nrecords):
      for p in range(len(prns_sorted)):
            prn = prns_sorted[p]
            dt, sec_of_day, bias = clock_records[prn][i]
            clock_bias = custom_scientific_notation(bias, precision=12)
            arbias = custom_scientific_notation(0.0, precision=12)
            #out.write(f"AS {prn}       {dt.year:4d} {dt.month:02d} {dt.day:02d} {dt.hour:02d} {dt.minute:02d} {dt.second:010.7f}  1 {bias:21.12E}\n")
            #out.write(f"AS {prn}       {dt.year:4d} {dt.month:02d} {dt.day:02d} {dt.hour:02d} {dt.minute:02d} {dt.second:9.6f}  2   {clock_bias}  {arbias}\n")
            out.write(f"AS {prn}       {dt.year:4d} {dt.month:02d} {dt.day:02d} {dt.hour:02d} {dt.minute:02d} {dt.second:9.6f}  1   {clock_bias}\n")

def extract_clk(input_file, output_file, datestr):

   dd = datestr
   if not os.path.exists(input_file):
      exit
   clock_records = read_records(input_file)
   n_sats = len(clock_records)
   sats = sorted(clock_records.keys())

   header_lines = code_headers(dd, sats, n_sats)
   write_clk_codefmt(output_file, header_lines, clock_records)
   print(f" RINEX clock bias file written to: {output_file}")

if __name__ == "__main__":
   import argparse

   parser = argparse.ArgumentParser(description='Main program')
   parser.add_argument("input", type=str,  help="input file")
   parser.add_argument("output", type=str,  help="output file")
   parser.add_argument("datestr", type=str,  help="date string in YYYYMMDDT00:00:00")
   parser.add_argument("-v","--verbose", action='store_true',  help="print info for debug")
 

   args = parser.parse_args()
   input_file= args.input ; output_file = args.output
   datestr = args.datestr
   if not os.path.exists(input_file):
      exit
   extract_clk(input_file, output_file, datestr)
