# ## Copyright (c) 2026 NOAA/NESDIS/STAR ## All rights reserved # #/**---------------------------------------------------------------------- # ionPhs_to_ionPrf_COSMIC_final.py # # This program takes as input L1 and L2 excess phase for an occultation # and computes an ionospheric electron density profile at the ray path # perigee points. # It also computes calibrated occultation TEC below the LEO orbit. # # Input: -- Level1b file (ionPhs) containing ionospheric excess phase # # Output: Vectors: electron density and calibrated TEC # Scalars: F0F2 (edmax) and HmF2 (critfreq) (plus many others) # # author: Yong Chen, Yong.Chen@noaa.gov # initial verion 04/12/2026 #-----------------------------------------------------------------------*/ import numpy as np from netCDF4 import Dataset from pyproj import Transformer #from datetime import datetime, UTC from datetime import datetime, timezone from frame_conversions import ro_ecef2eci, ro_eci2ecef from pygeodesy import GeoidEGM96 # load EGM96 coefficient for undulation geoid = GeoidEGM96("./geoids/WW15MGH.GRD") # consant K_ION = 40.30819294518825 #40.3 TECU = 1e16 FILL = -999.0 OMEGA_E = 7.2921150e-5 # rad/s c0_speed = 299792458.0 # m/s #User-configurable parameters from dataclasses import dataclass @dataclass class Parms: # ---- altitude coverage (meters) ---- tp_min: float = 150e3 # 150 km h_sat_margin: float = 10e3 # 10 km below satellite altitude fitlim: float = 5e3 # 5 km below satellite altitude dhmin: float = 1e3 # 1 km minimu difference between points in vertical grid edlim: float = 1e7 # maximum check limit for peak density (1/cm3) icutoff : int = 1 # Cutoff the first few seconds to get rid of acquisition errors # ---- calibration ---- use_aux_calibration: bool = False # use aux for calibration only mode 1 only #----- sampling rate --- sampling_rate: float = 1.0 # 1Hz and 1.0 second sigmatec: float = 0.01 # assumed standard deviation for the TEC near the top (TECU) eta: float = 7.5 # a priori scale height parameter = sqrt(porb/hscale) (unit-less) sigmaeta: float = 0.05 # assumed standard deviation of scale height parameter (unit-less) #local ENU unit vectors def enu_basis(lat, lon): """ lat, lon: arrays or scalars (degrees) returns: east, north, up with shape (..., 3) """ latr = np.deg2rad(lat) lonr = np.deg2rad(lon) east = np.stack(( -np.sin(lonr), np.cos(lonr), np.zeros_like(lonr) ), axis=-1) north = np.stack(( -np.sin(latr)*np.cos(lonr), -np.sin(latr)*np.sin(lonr), np.cos(latr) ), axis=-1) up = np.stack(( np.cos(latr)*np.cos(lonr), np.cos(latr)*np.sin(lonr), np.sin(latr) ), axis=-1) return east, north, up def unit(v): return v / np.linalg.norm(v, axis=-1, keepdims=True) def vector_angle(a, b, ref): """ Signed angle between a and b, using ref for sign. """ a_u = unit(a) b_u = unit(b) cross_ab = np.cross(a_u, b_u) sinang = np.linalg.norm(cross_ab, axis=-1) cosang = np.sum(a_u * b_u, axis=-1) angle = np.arctan2(sinang, cosang) sign = np.sign(np.sum(ref * cross_ab, axis=-1)) return angle * sign def occ_side_index(r_gps, r_leo, r_tp, eps=1e-6): """ Determine occultation geometry index: -1 : GPS and LEO on opposite sides of tangent point +1 : GPS and LEO on same side 0 : indeterminate / near-orthogonal """ v_gps = r_gps - r_tp v_leo = r_leo - r_tp dot = np.einsum("ij,ij->i", v_gps, v_leo) idx = np.zeros_like(dot, dtype=int) idx[dot > eps] = +1 idx[dot < -eps] = -1 return idx import numpy as np from scipy.signal import savgol_filter def compute_sat_velocity(r_eci, t, window=9, poly=3): """ Compute satellite velocity from ECI positions using Savitzky-Golay derivative. Parameters ---------- r_eci : ndarray (N,3) Satellite position in ECI coordinates [meters] t : ndarray (N,) Time array [seconds] (e.g., GPS seconds) window : int Savitzky-Golay window length (must be odd) poly : int Polynomial order for filter Returns ------- v_eci : ndarray (N,3) Satellite velocity in ECI coordinates [m/s] speed : ndarray (N,) Satellite speed magnitude [m/s] """ dt = np.mean(np.diff(t)) vx = savgol_filter(r_eci[:, 0], window, poly, deriv=1, delta=dt) vy = savgol_filter(r_eci[:, 1], window, poly, deriv=1, delta=dt) vz = savgol_filter(r_eci[:, 2], window, poly, deriv=1, delta=dt) v_eci = np.column_stack((vx, vy, vz)) speed = np.linalg.norm(v_eci, axis=1) return v_eci, speed def compute_velocity(r, t): r = np.asarray(r) # Convert time safely if hasattr(t, "unix"): # astropy time t = (t - t[0]).sec else: t = np.asarray(t) v = np.zeros_like(r) dt = t[2:] - t[:-2] v[1:-1] = (r[2:] - r[:-2]) / dt[:, None] v[0] = (r[1] - r[0]) / (t[1] - t[0]) v[-1] = (r[-1] - r[-2]) / (t[-1] - t[-2]) return v # elev=arcsin(s*rleo/(|s|*|rleo|)) def elevation_angle(r_leo, r_gns): """ Compute GNSS elevation angle at LEO (radians) Parameters ---------- r_leo : (N,3) ndarray LEO satellite positions (ECEF, meters) r_gns : (N,3) ndarray GNSS satellite positions (ECEF, meters) Returns ------- elev : ndarray Elevation angle (radians). Negative = occultation side """ s = r_gns - r_leo # normalize s_norm = np.linalg.norm(s, axis=1) r_norm = np.linalg.norm(r_leo, axis=1) # cosine of angle between r_leo and s cos_theta = np.einsum('ij,ij->i', s, r_leo) / (s_norm * r_norm) # elevation relative to local horizon elev = np.degrees(np.arcsin(cos_theta)) # degree return elev def tangent_point_vectorized(r_leo, r_gns, r_coc): """ Python equivalent of ROM SAF tangent_point.f90 Straight-line (no bending), vectorized """ pa = np.array([0.0, 0.0, 1.0]) # polar axis npts = r_leo.shape[0] # Shift to center-of-curvature frame r_leo_new = r_leo - r_coc r_gns_new = r_gns - r_coc # Line-of-sight vectors s = r_gns_new - r_leo_new s_hat = unit(s) # Impact parameter (straight-line) ip = np.linalg.norm(np.cross(r_leo_new, s_hat), axis=1) ro = np.linalg.norm(r_leo_new, axis=1) alpha = np.arccos(ip / ro) # Rotation axis axis = unit(np.cross(r_leo_new, r_gns_new)) # Rodrigues rotation (vectorized) v = r_leo_new perigee = ( v * np.cos(alpha)[:, None] + np.cross(axis, v) * np.sin(alpha)[:, None] + axis * np.sum(axis * v, axis=1)[:, None] * (1 - np.cos(alpha))[:, None] ) * (ip / ro)[:, None] occ_index = occ_side_index(r_gns, r_leo, perigee) elev = elevation_angle(r_leo, r_gns) # ---- Lat / Lon ---- rnorm = np.linalg.norm(perigee, axis=1) lat_tp = 90.0 - np.degrees(np.arccos(perigee[:, 2] / rnorm)) lon_tp = np.degrees(np.arctan2(perigee[:, 1], perigee[:, 0])) # ---- Azimuth (profile!) ---- v1 = np.cross(perigee, pa) v2 = np.cross(r_gns - r_coc, r_leo - r_coc) theta = vector_angle(v1, v2, -perigee) azimuth_tp = np.degrees(theta) azimuth_tp[azimuth_tp < 0] += 360.0 return lat_tp, lon_tp, azimuth_tp, ip, perigee, occ_index, elev # Use a representative mean ionospheric tangent point for r_coc, roc, and undulation # because ionospheric signal energy is concentrated near 300?400 km, # not near the neutral atmosphere as RO profile (occ point) # Using mean of tangent points above ~150 km def representative_point_vectorized(dtime, r_leo, r_gns, vleo): """ Representative Earth geometry for ionospheric RO """ # Earth constants (WGS-84) Re = 6378137.0 #meter Rp = 6356752.314245 # -------------------------------------------------- # 1. Straight-line tangent points (vectorized) # -------------------------------------------------- s = r_gns - r_leo s_hat = s / np.linalg.norm(s, axis=1)[:, None] p = np.linalg.norm(np.cross(r_leo, s_hat), axis=1) ro = np.linalg.norm(r_leo, axis=1) alpha = np.arccos(p / ro) axis = np.cross(r_leo, r_gns) axis /= np.linalg.norm(axis, axis=1)[:, None] # Rodrigues rotation (vectorized) v = r_leo v_rot = ( v * np.cos(alpha)[:, None] + np.cross(axis, v) * np.sin(alpha)[:, None] + axis * np.sum(axis * v, axis=1)[:, None] * (1 - np.cos(alpha))[:, None] ) perigee = v_rot * (p / ro)[:, None] # -------------------------------------------------- # 2. Find representative ionospheric tangent point perigee # -------------------------------------------------- # perigee: (N, 3) array in ECEF meters lat_per, lon_per, height_per = cart2geod(perigee) iono_mask = height_per > 150e3 lat_rep = np.mean(lat_per[iono_mask]) lon_rep = np.mean(lon_per[iono_mask]) perigee_rep = np.mean(perigee[iono_mask], axis=0) #iocc = np.argmin(np.abs(height_per)) #lat_occ = lat_per[iocc] # degrees #lon_occ = lon_per[iocc] # degrees #time_offset = dtime[iocc] # -------------------------------------------------- # 3. First lat/lon estimate # -------------------------------------------------- # -------------------------------------------------- # 4. Center of curvature refinement # -------------------------------------------------- latr = np.deg2rad(lat_rep) lonr = np.deg2rad(lon_rep) print("perigee lat, lon:", lat_rep, lon_rep, latr,lonr) n = np.array([ np.cos(latr) * np.cos(lonr), np.cos(latr) * np.sin(lonr), np.sin(latr) ]) r_coc = perigee_rep - Re * n print("perigee_rep:", perigee_rep ) print("r_coc norm (km):", np.linalg.norm(r_coc) / 1000) # perigee relative to CoC perigee_coc = perigee_rep - r_coc lat = 90.0 - np.degrees( np.arccos(perigee_coc[2] / np.linalg.norm(perigee_coc)) ) lon = np.degrees(np.arctan2(perigee_coc[1], perigee_coc[0])) # -------------------------------------------------- # 5. Azimuth (ROM SAF / UCAR definition) # -------------------------------------------------- pa = np.array([0.0, 0.0, 1.0]) n1 = np.cross(perigee_coc, pa) n2_all = np.cross( r_gns[iono_mask] - r_coc, r_leo[iono_mask] - r_coc ) # normalize each n2_all /= np.linalg.norm(n2_all, axis=1)[:, None] # representative plane normal n2 = np.mean(n2_all, axis=0) n1 /= np.linalg.norm(n1) n2 /= np.linalg.norm(n2) azimuth = np.degrees(np.arctan2( np.dot(perigee_coc, np.cross(n1, n2)), np.dot(n1, n2) )) if azimuth < 0: azimuth += 360.0 # -------------------------------------------------- # 6. Radius of curvature (ROM SAF exact) # -------------------------------------------------- theta = np.deg2rad(azimuth) r_coc, roc = curvature(lat, lon, theta) # undulation at the occ point undulation = geoid.height(lat, lon) print('undulation (m):', undulation) # -------------------------------------------------- # 7. Mean LEO altitude # -------------------------------------------------- sat_alt = np.linalg.norm(r_leo - r_coc, axis=1) - roc r_alt = np.mean(sat_alt) # -------------------------------------------------- # 8. View angle # -------------------------------------------------- vec = np.mean(r_gns[iono_mask] - r_leo[iono_mask], axis=0) vleo_mean = np.mean(vleo[iono_mask], axis=0) viewang = np.degrees(np.arccos( np.dot(vec, vleo_mean) / (np.linalg.norm(vec) * np.linalg.norm(vleo_mean)) )) return (lat, lon, r_coc, roc, undulation, r_alt, azimuth, viewang) def curvature(lat, lon, theta): """ Python equivalent of ROM SAF curvature.f90 Parameters ---------- lat : float Latitude [deg] lon : float Longitude [deg] theta : float Azimuth angle [rad] Returns ------- r_coc : ndarray (3,) Center of curvature in ECEF [m] roc : float Radius of curvature [m] """ # WGS-84 Re = 6378137.0 Rp = 6356752.314245 latr = np.deg2rad(lat) lonr = np.deg2rad(lon) cosl = np.cos(latr) sinl = np.sin(latr) # ---- surface point (h = 0) ---- x = geod2cart(lat, lon, 0.0) # ---- curvature radii ---- atoc2 = (Re / Rp) ** 2 xy2 = x[0]**2 + x[1]**2 z2 = (atoc2 * x[2])**2 cphi2 = xy2 / (xy2 + z2) sphi2 = z2 / (xy2 + z2) r_EW = 1.0 / np.sqrt(cphi2 + sphi2 / atoc2) r_NS = r_EW**3 / atoc2 roc = Re / (np.cos(theta)**2 / r_NS + np.sin(theta)**2 / r_EW) # ---- surface normal ---- n = np.array([ cosl * np.cos(lonr), cosl * np.sin(lonr), sinl ]) r_coc = x - roc * n return r_coc, roc # WGS-84 constants (must match ROM SAF) Re = 6378137.0 Rp = 6356752.314245 e = np.sqrt(1.0 - (Rp / Re)**2) rad2deg = 180.0 / np.pi deg2rad = np.pi / 180.0 f = 1.0 / 298.2572235630 def cart2geod_scalar(cart): """ Python equivalent of cart2geod_0d.f90 Parameters ---------- cart : array_like (3,) ECEF position vector [m] Returns ------- lat : float Geodetic latitude [deg] lon : float Longitude [deg] height : float Ellipsoidal height [m] """ x, y, z = cart Rxy = np.sqrt(x*x + y*y) Theta = np.arctan2(z * Re, Rxy * Rp) lon = np.arctan2(y, x) * rad2deg sinT = np.sin(Theta) cosT = np.cos(Theta) rlat = np.arctan2( z + e * Rp * sinT**3, Rxy - e * Re * cosT**3 ) lat = rlat * rad2deg cosPhi = np.cos(rlat) sinPhi = np.sin(rlat) height = -( -cosPhi * Rp**2 * Rxy - Re**2 * sinPhi * z + Re * Rp * np.sqrt( cosPhi**2 * Rp**2 + Re**2 * sinPhi**2 - sinPhi**2 * Rxy**2 + 2 * cosPhi * sinPhi * Rxy * z - cosPhi**2 * z**2 ) ) / (cosPhi**2 * Rp**2 + Re**2 * sinPhi**2) return lat, lon, height def cart2geod(cart): """ Vectorized Python equivalent of cart2geod.f90 Parameters ---------- cart : ndarray Shape (3,) or (N,3), ECEF coordinates [m] Returns ------- lat : float or ndarray Latitude [deg] lon : float or ndarray Longitude [deg] height : float or ndarray Ellipsoidal height [m] """ cart = np.asarray(cart, dtype=float) # ---- single point ---- if cart.ndim == 1: if cart.size != 3: raise ValueError("cart must have 3 elements") return cart2geod_scalar(cart) # ---- vectorized ---- if cart.ndim != 2 or cart.shape[1] != 3: raise ValueError("cart must have shape (N,3)") x = cart[:, 0] y = cart[:, 1] z = cart[:, 2] Rxy = np.sqrt(x*x + y*y) Theta = np.arctan2(z * Re, Rxy * Rp) lon = np.arctan2(y, x) * rad2deg sinT = np.sin(Theta) cosT = np.cos(Theta) rlat = np.arctan2( z + e * Rp * sinT**3, Rxy - e * Re * cosT**3 ) lat = rlat * rad2deg cosPhi = np.cos(rlat) sinPhi = np.sin(rlat) height = -( -cosPhi * Rp**2 * Rxy - Re**2 * sinPhi * z + Re * Rp * np.sqrt( cosPhi**2 * Rp**2 + Re**2 * sinPhi**2 - sinPhi**2 * Rxy**2 + 2 * cosPhi * sinPhi * Rxy * z - cosPhi**2 * z**2 ) ) / (cosPhi**2 * Rp**2 + Re**2 * sinPhi**2) return lat, lon, height def geod2cart(lat, lon, height): """ Vectorized Python equivalent of geod2cart.f90 Parameters ---------- lat : float or ndarray Latitude [deg] lon : float or ndarray Longitude [deg] height : float or ndarray Height above ellipsoid [m] Returns ------- cart : ndarray Shape (3,) or (N,3), ECEF coordinates [m] """ lat = np.asarray(lat, dtype=float) lon = np.asarray(lon, dtype=float) height = np.asarray(height, dtype=float) rlat = lat * deg2rad rlon = lon * deg2rad flatfn = (2.0 - f) * f funsq = (1.0 - f)**2 gd = Re / np.sqrt(1.0 - flatfn * np.sin(rlat)**2) x = np.cos(rlat) * np.cos(rlon) * (gd + height) y = np.cos(rlat) * np.sin(rlon) * (gd + height) z = np.sin(rlat) * (gd * funsq + height) cart = np.stack((x, y, z), axis=-1) return cart import numpy as np import warnings def reorder_occultation_arcs( p, impact_height, elev, time, ex_ph, *arrays, min_samples=3 ): """ Reorder RO samples into occultation and auxiliary arcs using p_max, with elevation-angle consistency checks and degenerate-arc handling. Parameters ---------- p : 1D array Impact parameter impact_height : 1D array Impact height (m), must be > 0 elev : 1D array Elevation angle (rad), negative = occultation elev : 1D array Elevation angle (rad), negative = occultation ex_ph : 1D array excess phase (L1 - L2) array *arrays : additional arrays to reorder min_samples : int Minimum number of samples required for a valid arc Returns ------- dict with keys: 'occultation' : tuple or None 'auxiliary' : tuple or None """ # ------------------------------------------------------------------ # 1. Convert & validate # ------------------------------------------------------------------ p = np.asarray(p) impact_height = np.asarray(impact_height) elev = np.asarray(elev) time = np.asarray(time) ex_ph = np.asarray(ex_ph) n = len(p) if not all(len(a) == n for a in (impact_height, elev, time, ex_ph, *arrays)): raise ValueError("All input arrays must have the same length") # ------------------------------------------------------------------ # 2. Filter by impact height > 0, and excess phase (l1-L2) is not far away form median # ------------------------------------------------------------------ median = np.median(ex_ph) valid = (impact_height > 0) & (abs(ex_ph - median) < 1000.0) n_valid = np.sum(valid) if n_valid < min_samples: warnings.warn( f"Only {n_valid} samples with impact_height > 0; " "cannot form valid arcs." ) return {'occultation': None, 'auxiliary': None} p = p[valid] impact_height = impact_height[valid] elev = elev[valid] time = time[valid] ex_ph = ex_ph[valid] arrays = [a[valid] for a in arrays] # ------------------------------------------------------------------ # 3. Identify p_max # ------------------------------------------------------------------ i_pmax = np.argmax(p) # ------------------------------------------------------------------ # 4. Split around p_max (including p_max itself) # ------------------------------------------------------------------ idx_left = np.arange(0, i_pmax+1) idx_right = np.arange(i_pmax, len(p)) if len(idx_left) < min_samples and len(idx_right) < min_samples: warnings.warn( "Both sides around p_max are degenerate; " "unable to split occultation/auxiliary arcs." ) return {'occultation': None, 'auxiliary': None} # ------------------------------------------------------------------ # 5. Determine occultation side using elevation majority # ------------------------------------------------------------------ neg_left = np.sum(elev[idx_left] < 0) neg_right = np.sum(elev[idx_right] < 0) if neg_left >= neg_right: idx_occ = idx_left idx_aux = idx_right else: idx_occ = idx_right idx_aux = idx_left # ------------------------------------------------------------------ # 6. Elevation consistency warnings # ------------------------------------------------------------------ if len(idx_occ) > 0: frac_pos = np.mean(elev[idx_occ] > 0) if frac_pos > 0.1: warnings.warn( f"Occultation arc contains {frac_pos:.1%} " "positive elevation samples." ) if len(idx_aux) > 0: frac_neg = np.mean(elev[idx_aux] < 0) if frac_neg > 0.1: warnings.warn( f"Auxiliary arc contains {frac_neg:.1%} " "negative elevation samples." ) # ------------------------------------------------------------------ # 7. Handle degenerate arcs # ------------------------------------------------------------------ if len(idx_occ) < min_samples: warnings.warn( f"Occultation arc degenerate ({len(idx_occ)} samples); " "setting to None." ) idx_occ = None if len(idx_aux) < min_samples: warnings.warn( f"Auxiliary arc degenerate ({len(idx_aux)} samples); " "setting to None." ) idx_aux = None # ------------------------------------------------------------------ # 8. Sort by increasing impact parameter # ------------------------------------------------------------------ def extract(indices): if indices is None: return None indices = indices[np.argsort(p[indices])] return ( p[indices], impact_height[indices], elev[indices], time[indices], ex_ph[indices], *[a[indices] for a in arrays] ) return { 'occultation': extract(idx_occ), 'auxiliary': extract(idx_aux) } # calibration mode determination def determine_calibration_mode(p_occ, p_aux): """ Determine calibration mode. Mode 2: auxiliary-side coverage available, but high within 50 km of occ side Mode 1: auxiliary-side coverage available with full cover of occ side Mode 0: no auxiliary-side coverage Parameters ---------- p_occ : 1D array Impact parameters on occultation side p_aux : 1D array or None Impact parameters on auxiliary side Returns ------- int calibration_mode (0 or 1) """ if p_aux is None: return 0 p_occ = np.asarray(p_occ) p_aux = np.asarray(p_aux) if p_occ.size == 0 or p_aux.size == 0: return 0 if (np.min(p_aux) <= np.min(p_occ)) and \ (np.max(p_aux) >= np.max(p_occ)): return 1 #if (np.min(p_aux) - 5e4 <= np.min(p_occ)) and \ # (np.max(p_aux) >= np.max(p_occ)): # return 2 return 0 import numpy as np import warnings def check_occultation_altitude_coverage( p_occ, h_occ, p_aux=None, h_aux=None, sat_alt=800e3, h_min=150e3, h_sat_margin=1e3, use_aux_for_calibration=False ): """ Check altitude and impact-parameter coverage requirements. Returns ------- bool True if occultation passes QC, False otherwise """ # ------------------------------- # 1. Occultation-side altitude coverage # ------------------------------- h_occ = np.asarray(h_occ) if h_occ.size == 0: warnings.warn("Empty occultation-side altitude array.") return False if np.min(h_occ) > h_min: warnings.warn( f"Occultation does not reach minimum altitude " f"({np.min(h_occ)/1e3:.1f} km > {h_min/1e3:.1f} km)." ) return False if np.max(h_occ) < (sat_alt - h_sat_margin): warnings.warn( f"Occultation does not reach close to satellite altitude " f"({np.max(h_occ)/1e3:.1f} km < " f"{(sat_alt-h_sat_margin)/1e3:.1f} km)." ) return False # ------------------------------- # 2. Auxiliary-side coverage (optional) # ------------------------------- if use_aux_for_calibration: if p_aux is None or len(p_aux) == 0: warnings.warn( "Auxiliary-side data required for calibration but missing." ) return False p_occ = np.asarray(p_occ) p_aux = np.asarray(p_aux) if (np.min(p_aux) > np.min(p_occ)) or (np.max(p_aux) < np.max(p_occ)): warnings.warn( "Auxiliary-side impact parameter range does not " "fully cover occultation-side range." ) return False return True def determine_occultation_mode(p_occ, time_occ): """ setting: 1 rising: 0 """ slope = np.polyfit(time_occ, p_occ, 1)[0] #return "setting" if slope < 0 else "rising" return 1 if slope < 0 else 0 # Check for time gaps import numpy as np import warnings def check_time_gaps(time, fs=1.0, tol=0.5): """ Check for time gaps in a time series. Parameters ---------- time : 1D array Time array (seconds) fs : float Sampling rate (Hz) tol : float Allowed tolerance in seconds Returns ------- bool True if time gaps detected, False otherwise """ time = np.asarray(time) if time.size < 2: return False dt = np.diff(time) expected_dt = 1.0 / fs return np.any(dt > (expected_dt + tol)) def check_occultation_time_gaps( time_occ, time_aux=None, fs=1.0, calibration_mode=0 ): """ Check for time gaps on both occultation and auxiliary sides. Returns ------- bool True if occultation should be discarded """ # Occultation side (mandatory) if check_time_gaps(time_occ, fs=fs): warnings.warn("Time gap detected on occultation side.") return True # Auxiliary side (if present) if (calibration_mode==1): if time_aux is not None: if check_time_gaps(time_aux, fs=fs): warnings.warn("Time gap detected on auxiliary side.") return True return False import numpy as np def topest_core(p1, h1, t1, i1, m1): """ Estimate orbit-level impact parameter, altitude, and time. Parameters ---------- p1 : array_like Impact parameters (km) h1 : array_like Tangent altitudes (km) t1 : array_like Time since occultation start (s) i1 : int Index of lowest point used m1 : int Size of input arrays (last valid index = m1-1) maxlen : int, optional Not used (kept for compatibility) Returns ------- porb : float Impact parameter at orbit altitude horb : float Tangent altitude at orbit altitude torb : float Time at orbit altitude """ # Select last up to 20 points imin = max(i1, m1 - 20) t = t1[imin:m1] p = p1[imin:m1] h = h1[imin:m1] jmax = 5 t0 = t1[m1 - 1] dt = t - t0 G = np.vstack([dt**j for j in range(jmax)]).T a, *_ = np.linalg.lstsq(G, p, rcond=None) b, *_ = np.linalg.lstsq(G, h, rcond=None) torb = t0 for _ in range(5): dt0 = torb - t0 dpdt = sum(j * a[j] * dt0**(j - 1) for j in range(1, jmax)) d2pdt2 = sum(j * (j - 1) * a[j] * dt0**(j - 2) for j in range(2, jmax)) if abs(d2pdt2) < 1e-12: break torb -= dpdt / d2pdt2 dt0 = torb - t0 porb = sum(a[j] * dt0**j for j in range(jmax)) horb = sum(b[j] * dt0**j for j in range(jmax)) return porb, horb, torb def topest_driver(p1, h1, t1, p2=None, h2=None, t2=None, fitlim=5e3, iflg=1, starttime=0.0): """ Full logic including window selection + dual-side quadratic fallback. Parameters ---------- p1, h1, t1 : occ side arrays p2, h2, t2 : other side arrays (optional) fitlim : float Threshold for selecting fitting window iflg : int Direction flag (1 for setting or 0 for raising) starttime : float Returns ------- porb, horb, torb, edorbtime, m1, m2 """ p1 = np.asarray(p1) h1 = np.asarray(h1) t1 = np.asarray(t1) m1 = len(p1) if p2 is not None: p2 = np.asarray(p2) h2 = np.asarray(h2) t2 = np.asarray(t2) m2 = len(p2) # ------------------------------------------------- # 1. Determine i1 using fitlim # ------------------------------------------------- i = m1 - 1 while i > 0 and p1[i] > p1[m1 - 1] - fitlim: i -= 1 i1 = i if (m1 - i1) < 5: raise ValueError("Not enough points for fitting (need >=5)") # ------------------------------------------------- # 2. Two cases # ------------------------------------------------- if p2 is None or m2 == 1: # ---- Case 1: extrapolation using topest ---- porb, horb, torb = topest_core(p1, h1, t1, i1, m1) if porb < p1[m1 - 1]: raise RuntimeError("Extrapolated porb < last p1 value") else: # ---- Case 2: quadratic interpolation (3-point) ---- tA = t1[m1 - 2] tB = t1[m1 - 1] tC = t2[m2 - 2] pA = p1[m1 - 2] pB = p1[m1 - 1] pC = p2[m2 - 2] hA = h1[m1 - 2] hB = h1[m1 - 1] hC = h2[m2 - 2] sqdt1 = (tA - tC) * (tA - tB) sqdt0 = (tB - tA) * (tB - tC) sqdt_1 = (tC - tB) * (tC - tA) torb = 0.5 * ( pA * (tB + tC) / sqdt1 + pB * (tC + tA) / sqdt0 + pC * (tA + tB) / sqdt_1 ) / ( pA / sqdt1 + pB / sqdt0 + pC / sqdt_1 ) weight1 = (torb - tC) * (torb - tB) / sqdt1 weight0 = (torb - tA) * (torb - tC) / sqdt0 weight_1 = (torb - tB) * (torb - tA) / sqdt_1 porb = pA * weight1 + pB * weight0 + pC * weight_1 horb = hA * weight1 + hB * weight0 + hC * weight_1 # ------------------------------------------------- # 3. Update indices based on direction # ------------------------------------------------- if ((iflg == 1 and torb > t1[m1 - 1]) or (iflg == 0 and torb < t1[m1 - 1])): m1 -= 1 else: if p2 is not None: m2 -= 1 # ------------------------------------------------- # 4. Final time # ------------------------------------------------- edorbtime = starttime + torb return porb, horb, torb, edorbtime, m1, m2 import numpy as np def phase_normalization(ph1, ph2, t1, t2, m1, m2, torb, icalib): """ Apply phase calibration by removing phorb offset for calibration mode 1 Parameters ---------- ph1, ph2 : array_like Phase arrays (occ and other side) t1, t2 : array_like Time arrays m1, m2 : int Sizes (use first m1 / m2 elements) torb : float Orbit time icalib : int Calibration flag (1 = apply) Returns ------- ph1, ph2 : ndarray Updated phase arrays phorb : float or None """ ph1 = np.asarray(ph1).copy() ph2 = np.asarray(ph2).copy() t1 = np.asarray(t1) t2 = np.asarray(t2) phorb = None if icalib == 1: # --- Compute phorb (linear interpolation between two sides) --- denom = (t1[m1 - 1] - t2[m2 - 1]) if abs(denom) < 1e-12: raise ZeroDivisionError("t1[m1-1] and t2[m2-1] too close") phorb = ( ph1[m1 - 1] * (torb - t2[m2 - 1]) + ph2[m2 - 1] * (t1[m1 - 1] - torb) ) / denom # --- Subtract phorb from both arrays --- ph1[:m1] -= phorb ph2[:m2] -= phorb return ph1, ph2, phorb import numpy as np from scipy.interpolate import CubicSpline def spline_calibration(ph1, ph2, p1, p2, t1, t2, m1, m2, torb, porb, icalib): """ Perform spline interpolation and compute calibrated excess phase. Returns ------- dph1 : ndarray Excess phase (calibrated or uncalibrated) phc1 : ndarray or None Interpolated auxiliary phase on p1 grid """ ph1 = np.asarray(ph1) ph2 = np.asarray(ph2).copy() p1 = np.asarray(p1) p2 = np.asarray(p2).copy() t1 = np.asarray(t1) t2 = np.asarray(t2).copy() phc1 = None # ------------------------------------------------- # Spline interpolation (only if calibration enabled) # ------------------------------------------------- if icalib == 1: # --- Extend arrays with orbit point --- t2_ext = np.append(t2[:m2], torb) p2_ext = np.append(p2[:m2], porb) ph2_ext = np.append(ph2[:m2], 0.0) # --- Build transformed time coordinates --- dt2 = -np.abs(t2_ext - torb) # ensure ascending sqt2 = dt2**2 # enforce linear behavior near top # ---- Sort for spline (must be strictly increasing) ---- #idx = np.argsort(sqt2) idx = np.argsort(p2_ext) sqt2_sorted = sqt2[idx] p2_sorted = p2_ext[idx] dt2_sorted = dt2[idx] ph2_sorted = ph2_ext[idx] # --- Step 1: interpolate time^2 onto p1 grid --- spline1 = CubicSpline(p2_sorted, sqt2_sorted, bc_type="natural", extrapolate=False) sqtc1 = spline1(p1[:m1]) # --- Recover dt from squared time --- dtc1 = -np.sqrt(np.maximum(sqtc1, 0.0)) # --- Step 2: interpolate phase onto dt grid --- # dt2 must be increasing ? sort again idx2 = np.argsort(dt2_sorted) dt2_sorted2 = dt2_sorted[idx2] ph2_sorted2 = ph2_sorted[idx2] spline2 = CubicSpline(dt2_sorted2, ph2_sorted2, bc_type="natural", extrapolate=False) phc1 = spline2(dtc1) # ------------------------------------------------- # Compute excess phase # ------------------------------------------------- dph1 = np.zeros(m1) if icalib == 1: # calibrated dph1 = ph1[:m1] - phc1 else: # uncalibrated dph1 = ph1[:m1] - ph1[m1 - 1] return dph1, phc1 #TEC & calibration in TECU def phase_to_tec(phi, f1, f2): return (f1**2 * f2**2 / (K_ION * TECU * (f1**2 - f2**2))) * phi import numpy as np def tecfit1(p1, tec, m1, porb, fitlim=5e3): """ Fit a straight line to uppermost TEC values. Parameters ---------- p1 : array_like Impact parameters tec : array_like TEC values m1 : int Upper index (number of valid points) porb : float Orbit impact parameter fitlim : optional for lower index bound Returns ------- a : float Slope of fitted line topfit : float RMS fitting error """ p1 = np.asarray(p1) tec = np.asarray(tec) # ------------------------------------------------- # 0. Determine i1 using fitlim # ------------------------------------------------- i = m1 - 1 while i > 0 and p1[i] > p1[m1 - 1] - fitlim: i -= 1 i1 = i if (m1 - i1) < 5: raise ValueError("Not enough points for fitting (need >=5)") # ------------------------------------------------- # 1. Select up to last 100 points # ------------------------------------------------- imin = max(i1, m1 - 100) # Reverse order to match Fortran loop (m1 to imin) idx = np.arange(m1 - 1, imin - 1, -1) # ------------------------------------------------- # 2. Build x and y # ------------------------------------------------- x = np.sqrt(np.maximum(1.0 - p1[idx] / porb, 0.0)) y = tec[idx] imax = len(x) # ------------------------------------------------- # 3. Linear fit: y = a * x # (no intercept) # ------------------------------------------------- Syx = np.sum(y * x) Sxx = np.sum(x * x) if Sxx < 1e-12: raise ZeroDivisionError("Sxx too small, cannot compute slope") a = Syx / Sxx # ------------------------------------------------- # 4. RMS fitting error # ------------------------------------------------- residual = y - a * x topfit = np.sqrt(np.sum(residual**2) / imax) return a, topfit import numpy as np def tecfit0(p1, tec, m1, porb, sigmatec, sigmaeta, eta, a_init, fitlim=5e3, max_iter=100, tol=0.1): """ Nonlinear fit for TEC top parameters (a, b). Parameters ---------- p1 : array_like Impact parameters tec : array_like TEC observations m1 : int Number of valid points porb : float Orbit impact parameter sigmatec : float Std dev of TEC observations sigmaeta : float Std dev of scale height parameter eta : float A priori scale height parameter a_init : float Initial guess for 'a' fitlim : optional for lower index bound max_iter : int Iteration count (default 100) tol : float Convergence tolerance (same logic as Fortran) Returns ------- a : float b : float topfit : float """ p1 = np.asarray(p1) tec = np.asarray(tec) pi = np.pi # ------------------------------------------------- # 0. Determine i1 using fitlim # ------------------------------------------------- i = m1 - 1 while i > 0 and p1[i] > p1[m1 - 1] - fitlim: i -= 1 i1 = i if (m1 - i1) < 5: raise ValueError("Not enough points for fitting (need >=5)") # ------------------------------------------------- # 1. Select up to last 100 points # ------------------------------------------------- imin = max(i1, m1 - 100) idx = np.arange(m1 - 1, imin - 1, -1) # ------------------------------------------------- # 2. Build x, y, z # ------------------------------------------------- sqrt_term = np.sqrt(np.maximum(1.0 - p1 / porb, 0.0)) x = (sqrt_term[idx] - sqrt_term[m1 - 1]) / sigmatec y = tec[idx] / sigmatec z = (p1[m1 - 1] - p1[idx]) / porb / sigmatec imax = len(x) # Append constraint point x = np.append(x, 0.0) y = np.append(y, eta / sigmaeta) z = np.append(z, 0.0) # will be updated each iteration # ------------------------------------------------- # 3. Iterative nonlinear solve # ------------------------------------------------- a = a_init b = 0.0 for _ in range(max_iter): # Update nonlinear constraint term z[-1] = 2.0 / np.sqrt(pi) / a / sigmaeta Syx = np.sum(y * x) Syz = np.sum(y * z) Sxx = np.sum(x * x) Sxz = np.sum(x * z) Szz = np.sum(z * z) c = Sxz**2 - Sxx * Szz if abs(c) < 1e-14: raise ZeroDivisionError("Singular system in tecfit0") b = (Syx * Sxz - Syz * Sxx) / c a = (Syz * Sxz - Syx * Szz) / c # ------------------------------------------------- # 4. Compute residual (exclude constraint point) # ------------------------------------------------- residual = y[:-1] - a * x[:-1] - b * z[:-1] topfit = sigmatec * np.sqrt(np.sum(residual**2) / imax) # ------------------------------------------------- # 5. Convergence check (same logic as Fortran) # ------------------------------------------------- z_last_expected = 2.0 / np.sqrt(pi) / a / sigmaeta if abs(1.0 - z_last_expected / z[-1]) > tol: raise RuntimeError("ERROR IN tecfit0: solution did not converge") return a, b, topfit import numpy as np from scipy.special import erfc # derfc in Fortran def process_tec(p1, tec, m1, porb, icalib, sigmatec=None, sigmaeta=None, eta=None, fitlim=5e3): """ Process TEC using tecfit1 / tecfit0 logic. Returns ------- tec : ndarray (updated if icalib=0) acoef, bcoef : floats edorb : float hscale : float topfit : float """ p1 = np.asarray(p1) tec = np.asarray(tec).copy() # ------------------------------------------------- # 1. First fit (tecfit1) # ------------------------------------------------- acoef, topfit = tecfit1(p1, tec, m1, porb, fitlim) # ------------------------------------------------- # 2. Branch depending on calibration flag # ------------------------------------------------- if icalib == 1: edorb = 1e10 * acoef / np.sqrt(8.0) / porb hscale = -999.0 bcoef = None else: # ---- Prepare initial guess ---- acoef = max(1.0, acoef) # ---- Nonlinear fit ---- acoef, bcoef, topfit = tecfit0( p1, tec, m1, porb, sigmatec, sigmaeta, eta, acoef, fitlim) # ---- Check solution validity ---- if acoef < 0.0 or bcoef < 0.0: raise RuntimeError("Invalid solution: acoef or bcoef < 0") # ---- Derived parameters ---- edorb = 1e10 * acoef / np.sqrt(2.0) / porb hscale = porb / ((4.0 / np.pi) * (bcoef / acoef)**2) # ---- TEC correction ---- eps0 = (porb - p1[m1 - 1]) / porb for i in range(m1): beta = (porb - p1[i]) / hscale tec[i] = ( tec[i] + acoef * np.sqrt(eps0) + bcoef * eps0 + bcoef * hscale / porb * (1.0 - np.exp(beta) * erfc(np.sqrt(beta))) ) hscale = hscale/1000.0 return tec, acoef, bcoef, edorb, hscale, topfit # Prevents oversampling and reduces correlated noise import numpy as np def filter_by_height(horb, m1, t1, h1, p1, tec, *arrays, dhmin=1e3): """ Remove points too close in height (within dhmin). Parameters ---------- horb : float Top (orbit) height m1 : int Number of valid points t1, h1, p1, tec : ndarray *arrays : additional aligned arrays dhmin : float Minimum vertical spacing Returns ------- m1_new, t1, h1, p1, tec, *arrays """ if any(len(a) < m1 for a in [t1, h1, p1, tec, *arrays]): raise ValueError("Some arrays are shorter than m1") # Convert to numpy and truncate to valid size t1 = np.asarray(t1)[:m1] h1 = np.asarray(h1)[:m1] p1 = np.asarray(p1)[:m1] tec = np.asarray(tec)[:m1] arrays = [np.asarray(a)[:m1] for a in arrays] keep = np.ones(m1, dtype=bool) hlast = horb # Traverse from top to bottom for i in range(m1 - 1, -1, -1): if h1[i] > hlast - dhmin: keep[i] = False else: hlast = h1[i] # Apply mask t1 = t1[keep] h1 = h1[keep] p1 = p1[keep] tec = tec[keep] arrays = [a[keep] for a in arrays] m1_new = len(h1) return (m1_new, t1, h1, p1, tec, *arrays) import numpy as np def invtec(p, tec, porb, edorb): """ Invert TEC to electron density profile. using straight signal path assumption, spherical symmetry assumption, and circular LEO orbit assumption. Parameters ---------- p : array_like Impact parameters (size n) tec : array_like TEC values (size n) porb : float Orbit impact parameter edorb : float Electron density at orbit Returns ------- ed : ndarray Electron density profile (size n) """ p = np.asarray(p) tec = np.asarray(tec) n = len(p) ed = np.zeros(n) c0 = np.zeros(n) c = np.zeros((n, n)) # ------------------------------------------------- # 1. Compute coefficients c0 and c # ------------------------------------------------- for i in range(n): # eps array (size up to n-i) eps = np.zeros(n - i ) for k in range(n - i -1 ): eps[k] = (p[i + k +1 ] - p[i]) / p[i] eps[n - i -1] = (porb - p[i]) / p[i] # Initialize ep0 = 1.0 ep1 = 1.0 + eps[0] sq0 = 0.0 sq1 = np.sqrt(max(eps[0] * (1.0 + ep1), 0.0)) # c0 if abs(eps[0]) < 1e-12: c0[i] = 1.0 else: c0[i] = (ep1 * sq1 - np.log(ep1 + sq1)) / eps[0] # main loop for k in range(n - i -1): ep_1 = ep0 ep0 = ep1 ep1 = 1.0 + eps[k + 1] sq_1 = sq0 sq0 = sq1 sq1 = np.sqrt(max(eps[k + 1] * (1.0 + ep1), 0.0)) denom1 = ep1 - ep0 denom2 = ep0 - ep_1 term1 = 0.0 term2 = 0.0 if abs(denom1) > 1e-12: term1 = ( ep1 * (sq1 - sq0) - np.log((ep1 + sq1) / (ep0 + sq0)) ) / denom1 if abs(denom2) > 1e-12: term2 = ( ep_1 * (sq0 - sq_1) - np.log((ep0 + sq0) / (ep_1 + sq_1)) ) / denom2 c[k, i] = term1 - term2 # last coefficient k = n - i denom = ep1 - ep0 if abs(denom) > 1e-12: c[k-1, i] = -( (2 * ep0 - ep1) * sq1 - ep0 * sq0 - np.log((ep1 + sq1) / (ep0 + sq0)) ) / denom else: c[k-1, i] = 0.0 # ------------------------------------------------- # 2. Back-substitution for electron density # ------------------------------------------------- for i in range(n - 1, -1, -1): ed[i] = 1e10 * tec[i] / p[i] for k in range(n - i - 1): ed[i] -= c[k, i] * ed[i + k + 1] ed[i] = (ed[i] - c[n - i - 1, i] * edorb) / c0[i] # force the boundary condtion #ed[n-1] = edorb return ed from geographiclib.geodesic import Geodesic def horizontal_smear(lat_bot, lon_bot, lat_top, lon_top): """ UCAR definition of horizontal smear: geodesic distance between bottom and top tangent points projected onto Earth's surface. Returns: meters """ g = Geodesic.WGS84 r = g.Inverse(lat_bot, lon_bot, lat_top, lon_top) return r["s12"] import numpy as np def find_f_layer_peak( alt_msl, # m Ne, # cm^-3 lat, # deg lon, # deg az, # deg time, # time alt_min=80e3, edlim=1e7 ): """ UCAR-style F-layer peak detection. Returns ------- edmax : max electron density (cm^-3) edmaxalt : altitude of peak (km) edmaxlat : latitude at peak (deg) edmaxlon : longitude at peak (deg) edmaxaz : azimuth angle at peak (deg) time : time at peak (deg) critfreq : critical frequency (MHz) """ alt_msl = np.asarray(alt_msl) Ne = np.asarray(Ne) # ------------------------------------------- # 1. Mask above 80 km # ------------------------------------------- mask = ( (alt_msl >= alt_min) & np.isfinite(Ne) & (Ne > 0) ) if np.count_nonzero(mask) < 3: raise RuntimeError("Insufficient valid Ne samples above 80 km") Ne_use = Ne[mask] alt_use = alt_msl[mask] lat_use = lat[mask] lon_use = lon[mask] az_use = az[mask] time_use = time[mask] # ------------------------------------------- # 2. Simple maximum search # ------------------------------------------- i = np.argmax(Ne_use) edmax = Ne_use[i] if edmax > edlim: raise RuntimeError("Insufficient valid Ne samples above 80 km") edmaxalt = alt_use[i] edmaxlat = lat_use[i] edmaxlon = lon_use[i] edmaxaz = az_use[i] edmaxtime = time_use[i] # ------------------------------------------- # 3. Critical frequency (MHz) # ------------------------------------------- critfreq = 8.98e-3 * np.sqrt(edmax) return ( edmax, edmaxalt, edmaxlat, edmaxlon, edmaxaz, edmaxtime, critfreq ) def vertical_tec(h1, ed, hmin=80e3): """ Compute TEC between 80 km and top. Parameters ---------- h1 : array_like Heights (m) ed : array_like Electron density hmin : optional input minimum height Returns ------- tec0 : float TEC (in TECU) """ h1 = np.asarray(h1) ed = np.asarray(ed) m1 = len(h1) # ------------------------------------------------- # 1. Find first index above 80 km # ------------------------------------------------- i = 0 while i < m1 and h1[i] < hmin: i += 1 if i >= m1: return 0.0 # no valid region n80 = i # ------------------------------------------------- # 2. Trapezoidal integration # ------------------------------------------------- tec0 = 0.0 for i in range(n80 + 1, m1): delta = (h1[i] - h1[i - 1]) * (ed[i] + ed[i - 1]) / 2.0 if delta < 0.0: delta = 0.0 tec0 += delta # ------------------------------------------------- # 3. Unit conversion # ------------------------------------------------- tec0 = 1e6 * tec0 # el/m^2 tec0 = 1e-16 * tec0 # TECU return tec0 import numpy as np def extrapolate_tec(h1, p1, ed, hmax): """ Extrapolate electron density upward and compute TEC contribution. Parameters ---------- h1 : array_like Heights (m) p1: array_like impact height (m) ed : array_like Electron density hmax: height at max electron density Returns ------- tec1 : float Extrapolated TEC (TECU), or -999 if invalid """ h1 = np.asarray(h1) p1 = np.asarray(p1) ed = np.asarray(ed) m1 = len(h1) # ------------------------------------------------- # 1. Find n_100 (point ~100 km below top) # ------------------------------------------------- i = m1 - 1 while i >= 0 and (h1[m1 - 1] - h1[i]) < 100.0e3: i -= 1 if i < 0: return -999.0 n_100 = i # ------------------------------------------------- # 2. Check height condition # ------------------------------------------------- if h1[n_100] < hmax + 50.0e3: return -999.0 # ------------------------------------------------- # 3. Number of regression points # ------------------------------------------------- nreg = m1 - n_100 # ------------------------------------------------- # 4. Check for negative densities # ------------------------------------------------- if np.any(ed[n_100:m1] < 0.0): return -999.0 # ------------------------------------------------- # 5. Build regression variables # ------------------------------------------------- xx = p1[m1 - 1] - p1[n_100:m1] yy = np.log(ed[n_100:m1]) # ------------------------------------------------- # 6. Linear regression: yy = c0 + c1 * xx # ------------------------------------------------- # Equivalent to Fortran lreg A = np.vstack([np.ones(nreg), xx]).T coeffs, residuals, _, _ = np.linalg.lstsq(A, yy, rcond=None) c0, c1 = coeffs # Optional: RMS (not used further, but matches Fortran) if len(residuals) > 0: rms = np.sqrt(residuals[0] / nreg) else: rms = 0.0 # ------------------------------------------------- # 7. Check slope # ------------------------------------------------- if c1 <= 0.0: return -999.0 # ------------------------------------------------- # 8. Compute extrapolated TEC # ------------------------------------------------- # zz(nreg) = fitted value at last point zz_last = c0 + c1 * xx[-1] tec1 = 1e6 * np.exp(zz_last) / c1 # el/m^2 tec1 = 1e-16 * tec1 # TECU return tec1 import numpy as np from scipy.interpolate import CubicSpline def interp_to_ucar_height( alt_occ, tec_occ, alt_ucar ): """ Interpolate TEC from this work onto UCAR height grid using natural cubic spline. All altitudes in km. """ alt_occ = np.asarray(alt_occ) tec_occ = np.asarray(tec_occ) alt_ucar = np.asarray(alt_ucar) # Remove NaNs mask = np.isfinite(alt_occ) & np.isfinite(tec_occ) alt_occ = alt_occ[mask] tec_occ = tec_occ[mask] # Ensure increasing altitude idx = np.argsort(alt_occ) alt_occ = alt_occ[idx] tec_occ = tec_occ[idx] # Ensure UCAR grid within range valid = ( (alt_ucar >= alt_occ.min()) & (alt_ucar <= alt_occ.max()) ) tec_interp = np.full_like(alt_ucar, np.nan, dtype=float) spline = CubicSpline( alt_occ, tec_occ, bc_type="natural", extrapolate=False ) tec_interp[valid] = spline(alt_ucar[valid]) return tec_interp import numpy as np from astropy.time import Time import astropy.units as u # convert gps time to local time def gps_to_lct(gps_time, lon_deg): """ Convert GPS time and longitude to Local Solar Time (LCT). Parameters ---------- gps_time : float or array-like GPS seconds lon_deg : float or array-like Longitude in degrees (east positive) Returns ------- lct : float or ndarray Local solar time in hours [0, 24) """ gps_time = np.asarray(gps_time) lon_deg = np.asarray(lon_deg) # GPS UTC (handles leap seconds correctly) t = Time(gps_time, format="gps").utc # UTC hour in decimal hours utc_hour = ( t.datetime.hour + t.datetime.minute / 60.0 + ( t.datetime.second + t.datetime.microsecond * 1e-6 ) / 3600.0 ) # Local solar time lct = utc_hour + lon_deg / 15.0 # wrap to [0, 24) return np.mod(lct, 24.0) import numpy as np def apply_cutoff(arrays, icutoff=5): """ Apply cutoff to a list of NumPy arrays. Parameters ---------- arrays : list of np.ndarray List of arrays with the same length (nmax). icutoff : int Number of elements to remove from the beginning. Returns ------- new_arrays : list of np.ndarray Sliced arrays after cutoff. new_nmax : int Updated length after cutoff. """ if icutoff <= 0: return arrays, len(arrays[0]) # Safety check nmax = len(arrays[0]) if any(len(arr) != nmax for arr in arrays): raise ValueError("All arrays must have the same length") if icutoff >= nmax: raise ValueError("icutoff must be smaller than array length") # Apply slicing new_arrays = [arr[icutoff:] for arr in arrays] new_nmax = nmax - icutoff return new_arrays, new_nmax import os import re from datetime import datetime, timedelta def parse_cdaac_filename(filename): """ General UCAR/CDAAC RO filename parser. Examples: ionPhs_C006.2009.001.23.55.G28_2021.0390_nc ionPrf_C2E4.2026.001.14.18.R15_0001.0001_nc """ fname = os.path.basename(filename) pattern = ( r'(?P\w+)_' r'(?P[A-Z0-9]+)\.' r'(?P\d{4})\.' r'(?P\d{3})\.' r'(?P\d{2})\.' r'(?P\d{2})\.' r'(?P[GRCE]\d{1,2})_' r'(?P\d+)\.' r'(?P\d+)' r'_nc' ) m = re.match(pattern, fname) if not m: raise ValueError(f"Unrecognized CDAAC filename: {fname}") d = m.groupdict() fileStamp = (f"{d['mission']}."f"{d['year']}."f"{d['doy']}."f"{d['hour']}."f"{d['minute']}."f"{d['gnss']}") # ---------------------------- # GNSS # ---------------------------- gnss_code = d["gnss"][0] gnss_prn = int(d["gnss"][1:]) gnss_map = { "G": "GPS", "R": "GLONASS", "E": "GALILEO", "C": "BEIDOU", } gnss_system = gnss_map.get(gnss_code, "UNKNOWN") # ---------------------------- # Time (UTC) # ---------------------------- year = int(d["year"]) doy = int(d["doy"]) hour = int(d["hour"]) minute = int(d["minute"]) start_datetime = ( datetime(year, 1, 1) + timedelta(days=doy - 1, hours=hour, minutes=minute) ) return { "fileStamp": fileStamp, "product": d["product"], "mission": d["mission"], "year": year, "doy": doy, "hour": hour, "minute": minute, "datetime_utc": start_datetime, "gnss_code": gnss_code, "gnss_system": gnss_system, "gnss_prn": gnss_prn, "occultation_id": int(d["occid"]), "arc_id": int(d["arcid"]), } import os import re from datetime import datetime def parse_STAR_podTec_filename(filename): """Parse STAR ROSW_TEC filenames used as phase-input containers. Example: ROSW_TEC_C2E2_20260921041001_20260921050340_G21_A01_c20260921183901_nes_ops.nc The EDP retrieval still uses exph_L1/exph_L2 as its science input; this parser only adapts the filename metadata to the TEC container naming. """ fname = os.path.basename(filename) pattern = ( r'ROSW_TEC_' r'(?P[A-Za-z0-9]+)_' r'(?P\d{14})_' r'(?P\d{14})_' r'(?P[GRCE]\d{1,2})_' r'(?PA\d{2})_' r'c(?P\d{14})_' r'(?P
[A-Za-z0-9]+)_ops\.nc$' ) m = re.match(pattern, fname) if not m: raise ValueError(f"Unrecognized ROSW_TEC filename: {fname}") d = m.groupdict() start_datetime = datetime.strptime(d["start"], "%Y%m%d%H%M%S") end_datetime = datetime.strptime(d["end"], "%Y%m%d%H%M%S") year = start_datetime.year doy = int(start_datetime.strftime("%j")) gnss_code = d["gnss"][0] gnss_prn = int(d["gnss"][1:]) gnss_map = {"G": "GPS", "R": "GLONASS", "E": "GALILEO", "C": "BEIDOU"} fileStamp = f"{d['sat']}.{year}.{doy:03d}.{start_datetime.hour:02d}.{start_datetime.minute:02d}.{d['gnss']}" return { "fileStamp": fileStamp, "product": "ROSW_TEC", "mission": d["sat"], "year": year, "doy": doy, "hour": start_datetime.hour, "minute": start_datetime.minute, "datetime_utc": start_datetime, "end_datetime_utc": end_datetime, "gnss_code": gnss_code, "gnss_system": gnss_map.get(gnss_code, "UNKNOWN"), "gnss_prn": gnss_prn, "arc_id": d["arcid"], } def get_gnss_frequencies(gnss_code, gnss_prn): """ Return GNSS L1/L2 frequencies (Hz) based on GNSS system code. Parameters ---------- gnss_code : str GNSS system code: G,R,E,C Returns ------- f1, f2 : float Carrier frequencies (Hz) """ k = [1, -4, 5, 6, 1, -4, 5, 6, -2, -7, 0, -1, -2, -7, 0, -1, 4, -3, 3, 2, 4, -3, 3, 2] freq = { "G": (1575.42e6, 1227.60e6), # GPS L1, L2 "R": ((1602 + k[gnss_prn-1] * 0.5625) * 1e6, (1246 + k[gnss_prn-1] * 0.4375) * 1e6), # GLONASS nominal "E": (1575.42e6, 1207.140e6), # Galileo L1C, L7Q "C": (1575.42e6, 1176.450e6) # BeiDou L1P, L5P } # For GLONASS, te frequences actually depend on the frequency channel k #f1 = (1602 + k * 0.5625) * 1e6 #f2 = (1246 + k * 0.4375) * 1e6 if gnss_code not in freq: raise ValueError(f"Unknown GNSS system: {gnss_code}") return freq[gnss_code] #INLINE ionPrf netCDF WRITER def write_ionPrf( filename, alt, lat, lon, az, tec_cal, Ne, meta ): n = len(alt) nc = Dataset(filename, "w", format="NETCDF4") # Dimension nc.createDimension("MSL_alt", n) def v(name): return nc.createVariable(name, "f4", ("MSL_alt",), fill_value=FILL) v_alt = v("MSL_alt") v_lat = v("GEO_lat") v_lon = v("GEO_lon") v_azi = v("OCC_azi") v_tec = v("TEC_cal") v_ed = v("ELEC_dens") v_alt.units = "km" v_alt.valid_range = np.array([0.0, 9999.0], dtype="f4") v_alt.long_name = "Mean sea level altitude of perigee point" v_lat.units = "degrees_north" v_lat.valid_range = np.array([-90.0, 90.0], dtype="f4") v_lat.long_name = "Geographical latitude of perigee point" v_lon.units = "degrees_east" v_lon.valid_range = np.array([-180.0, 180.0], dtype="f4") v_lon.long_name = "Geographical longitude of perigee point" v_azi.units = "deg" v_azi.valid_range = np.array([-180.0, 180.0], dtype="f4") v_azi.long_name = "Azimuth angle of occ. plane with respect to north" v_tec.units = "TECU" v_tec.valid_range = np.array([-1.e+08, 1.e+08], dtype="f4") v_tec.long_name = "Calibrated occultation TEC below LEO orbit" v_ed.units = "el/cm3" v_ed.valid_range = np.array([-1.e+08, 1.e+08], dtype="f4") v_ed.long_name = "Electron density" v_alt[:] = alt v_lat[:] = lat v_lon[:] = lon v_azi[:] = az v_tec[:] = tec_cal v_ed[:] = Ne # cm^-3 for k, val in meta.items(): setattr(nc, k, val) UTC = timezone.utc nc.creation_time = datetime.now(UTC).strftime("%d-%b-%y %H:%M").upper() nc.close() """ run_gmrion.py Operational driver for GMR ionospheric profile processing. Example: python run_gmrion.py \ --ionPhs input_ionPhs.nc \ --out output_ionPrf_directory Batch example: python run_gmrion.py --list filelist.txt """ import argparse import logging import sys import os #from datetime import datetime, UTC # --------------------------------------------------------- # Logging setup # --------------------------------------------------------- def setup_logging(level="INFO", logfile=None): log_format = "%(asctime)s | %(levelname)s | %(message)s" logging.basicConfig( level=getattr(logging, level.upper()), format=log_format, handlers=[logging.StreamHandler(sys.stdout)] ) if logfile: file_handler = logging.FileHandler(logfile) file_handler.setFormatter(logging.Formatter(log_format)) logging.getLogger().addHandler(file_handler) ##MASTER DRIVER (FULL ALGORITHM, END-TO-END) #def run_gmrion(ionPhs_file, ionPrf_file_UCAR, ionPrf_file_out): def run_gmrion(tec_file, edp_file_out_dir=None): ionPhs_file = tec_file # TEC netCDF is the container; exph_L1/exph_L2 remain the retrieval inputs ionPrf_file_out_dir = edp_file_out_dir from astropy.time import Time import astropy.units as u import os parms = Parms() logging.info(f"Processing TEC phase container: {ionPhs_file}") if not os.path.exists(ionPhs_file): raise FileNotFoundError(f"Input TEC file not found: {ionPhs_file}") logging.info("Running ionospheric electron density profile retrieval algorithm...") info = parse_STAR_podTec_filename(ionPhs_file) fileStamp = info["fileStamp"] print("fileStamp:", fileStamp) gnss_prn = info["gnss_prn"] gnss_code = info["gnss_code"] nc = Dataset(ionPhs_file) # Disable automatic masking to prevent the time # time:valid_range = 0., 360. issue # bug fix on 01/26/2026 nc.set_auto_mask(False) #time = nc["time"][:] exL1 = nc["exph_L1"][:] exL2 = nc["exph_L2"][:] # --- Time --- #startTime = nc.startTime # scalar, GPS seconds # the time already adding "add_offset" which is startTime for planetiq data startTime = 0.0 # scalar, GPS seconds # the time already adding "add_offset" which is startTime for planetiq data time_rel = nc['time'][:] # seconds since startTime print("time_rel shape:", time_rel.shape, time_rel.min(), time_rel.max()) # Absolute GPS seconds gps_seconds = startTime + time_rel # Build Astropy Time object # GPS time does not have leap seconds # GPS time = TAI - 19 s (constant offset since 1980) # GPS time = UTC + 15 seconds (leap seconds 15 s in 2009) t = Time(gps_seconds, format='gps', scale='tai') print(t[0].tai.iso) print(t[0].utc.iso) # need UTC (e.g., for Earth rotation) # Astropy automatically applies the correct leap second table for the year 2009 t_utc = t.utc # ECEF positions, convert ECEF km to meters r_leo_ecef = 1e3 * np.column_stack((nc['x_LEO'][:], nc['y_LEO'][:], nc['z_LEO'][:])) r_gps_ecef = 1e3 * np.column_stack((nc['x_GNSS'][:], nc['y_GNSS'][:], nc['z_GNSS'][:])) # If ECEF velocity is not provided # Compute velocities #v_leo_ecef, speed_leo = compute_sat_velocity(r_leo_eci, gps_seconds) #v_gps_ecef, speed_gps = compute_sat_velocity(r_gps_eci, gps_seconds) #print("LEO mean speed:", speed_leo.mean()/1000, "km/s") #print("GPS mean speed:", speed_gps.mean()/1000, "km/s") # --- ECI to ECEF --- #r_leo_ecef, v_leo_ecef = ro_eci2ecef(r_leo_eci, v_leo_eci, t.utc) #r_gps_ecef, v_gps_ecef = ro_eci2ecef(r_gps_eci, v_gps_eci, t.utc) v_leo_ecef = compute_velocity(r_leo_ecef, t.utc) v_gps_ecef = compute_velocity(r_gps_ecef, t.utc) #f1 = nc.occfreq1 #f2 = nc.occfreq2 f1, f2 = get_gnss_frequencies(gnss_code, gnss_prn) print(f1, f2) arrays = [time_rel, r_leo_ecef, r_gps_ecef, v_leo_ecef, v_gps_ecef, exL1, exL2] # Cutoff the first few seconds to get rid of acquisition errors arrays, nmax = apply_cutoff(arrays, icutoff=parms.icutoff) # Unpack back (time_rel, r_leo_ecef, r_gps_ecef, v_leo_ecef, v_gps_ecef, exL1, exL2) = arrays # Calculation of the excess phases ex_ph = exL1 - exL2 (occ_lat, occ_lon, r_coc, roc, undulation, mean_alt, occ_azi, viewang ) = \ representative_point_vectorized(time_rel, r_leo_ecef, r_gps_ecef, v_leo_ecef ) lat_tp, lon_tp, az_tp, p, r_per, occ_index, elev = tangent_point_vectorized( r_leo_ecef, r_gps_ecef, r_coc ) plotfig = 0 if plotfig: import matplotlib.pyplot as plt print("time_rel shape:", time_rel.shape, time_rel.min(), time_rel.max()) print("impact height shape:", p.shape ) impact_height = (p-roc-undulation)/1000.0 print("impact height min and max:", impact_height.min(), impact_height.max() ) plt.figure() plt.plot(time_rel, impact_height, linestyle='-') plt.xlabel('Time') plt.ylabel('Impact Height (km)') plt.title('Impact Height vs Time') plt.grid(True) plt.xlim(time_rel.min(), time_rel.max()) plt.ylim(impact_height.min(), impact_height.max() ) #plt.show() plt.savefig('piq_p_vs_time_v6.png', dpi=300, bbox_inches='tight') plt.close() plt.figure() plt.plot(time_rel, elev, linestyle='-') plt.xlabel('Time') plt.ylabel('Elevation Angle') plt.title('Elevation Angle vs Time') plt.grid(True) plt.xlim(time_rel.min(), time_rel.max()) plt.ylim(elev.min(), elev.max() ) #plt.show() plt.savefig('piq_elev_vs_time_v6.png', dpi=300, bbox_inches='tight') plt.close() # Occulation plane normal n_occ = np.cross(r_leo_ecef, r_gps_ecef) n_occ /= np.linalg.norm(n_occ) # Project into horizontal plane east, north, up = enu_basis(lat_tp, lon_tp) # project occultation-plane normal into horizontal plane n_h = n_occ - np.sum(n_occ * up, axis=1, keepdims=True) * up n_h /= np.linalg.norm(n_h, axis=1, keepdims=True) az = np.degrees(np.arctan2( np.sum(n_h * east, axis=1), np.sum(n_h * north, axis=1) )) #az = az % 360.0 az_tp = (az + 90.0) % 360.0 # convert the azimuth angle from 0-360 to -180 to 180 az_tp = ((az_tp + 180.0) % 360.0) - 180.0 print("mean_alt:", mean_alt) print("impact height range:", p.min()- roc-undulation, p.max()- roc-undulation) i0 = np.argmax(p) i1 = np.argmin(p) print("min and max impact height index:", i1, i0) #print("occ_index shape:", occ_index.shape, occ_index[1000]) print("occ_index shape:", occ_index.shape) # ------------------------------------------------------------------ # Step 3. Re-ordering of impact parameter and phase arrays at # both sides of the maximum impact parameter (for TEC # calibration) # ------------------------------------------------------------------ impact_height = p-roc-undulation reordered = reorder_occultation_arcs( p, impact_height, elev, time_rel, ex_ph, lat_tp, lon_tp, az_tp, r_per, r_leo_ecef, r_gps_ecef, exL1, exL2 ) # Occultation side occ = reordered['occultation'] if occ is None: raise RuntimeError("No occultation arc") else: p_occ, h_occ, elev_occ, time_occ, ex_ph_occ, lat_occ, lon_occ, az_occ, r_per_occ, \ r_leo_occ, r_gps_occ, exL1_occ, exL2_occ = occ #print("Occultation p:", p_occ.shape, p_occ.min()-roc-undulation, p_occ.max()-roc-undulation, p_occ-roc-undulation) print("Occultation p:", p_occ.shape, p_occ.min()-roc-undulation, p_occ.max()-roc-undulation) # Auxiliary side aux = reordered['auxiliary'] if aux is not None: p_aux, h_aux, elev_aux, time_aux, ex_ph_aux, lat_aux, lon_aux, az_aux, \ r_per_aux, r_leo_aux, r_gps_aux, exL1_aux, exL2_aux = aux print("Auxiliary p:", p_aux.shape, p_aux.min()-roc-undulation, p_aux.max()-roc-undulation) else: p_aux = None h_aux = None calibration_mode = determine_calibration_mode(p_occ, p_aux) print("calibration_mode:", calibration_mode) occ_mode = determine_occultation_mode(p_occ, time_occ) if calibration_mode == 0 and p_aux is not None: print( "Auxiliary-side coverage insufficient: " f"p_occ=[{p_occ.min()-roc:.1f}, {p_occ.max()-roc:.1f}], " f"p_aux=[{p_aux.min()-roc:.1f}, {p_aux.max()-roc:.1f}]" ) # ------------------------------------------------------------------ # Step 4. Check for sufficient altitude range # ------------------------------------------------------------------ ok = check_occultation_altitude_coverage( p_occ, h_occ, p_aux, h_aux, sat_alt=mean_alt, h_min=parms.tp_min, h_sat_margin=parms.h_sat_margin, use_aux_for_calibration=parms.use_aux_calibration ) if not ok: raise RuntimeError("Insufficient altitude range") # ------------------------------------------------------------------ # Step 5. Check for time gaps # ------------------------------------------------------------------ discard = check_occultation_time_gaps(time_occ, time_aux, fs=parms.sampling_rate, calibration_mode=calibration_mode) if discard: raise RuntimeError("Time gaps detected") # ------------------------------------------------------------------ # Step 6a. Estimate time and impact parameter at the top of the occulation # Notice that array sizes for OCC and AUX have been changed to m1_occ and m2_aux # ------------------------------------------------------------------ porb, horb, torb, edorbtime_rel, m1_occ, m2_aux = topest_driver(p_occ, h_occ, time_occ, \ p_aux, h_aux,time_aux, fitlim=parms.fitlim, iflg=occ_mode, starttime=0.0) # ------------------------------------------------------------------ # Step 6b. Estimation of the phase at torb by linear interpolation (for calibration # mode 1) and normalization of the phases (to be zero at the top) # ------------------------------------------------------------------ if calibration_mode == 1: ph_occ, ph_aux, phorb = phase_normalization(ex_ph_occ, ex_ph_aux, time_occ, time_aux, m1_occ, m2_aux, torb, calibration_mode) # ------------------------------------------------------------------ # Step 6c. Spline interpolation of the array ph2 onto the p1 grid # output the Calibrated/uncalibrated phases # ------------------------------------------------------------------ if calibration_mode == 1: dph1, phc1 = spline_calibration(ph_occ, ph_aux, p_occ, p_aux, time_occ, time_aux, m1_occ, m2_aux, torb, porb, calibration_mode) else: # uncalibrated dph1 = ex_ph_occ[:m1_occ] - ex_ph_occ[m1_occ - 1] print("excess phase:", dph1.shape, dph1.min(), dph1.max()) # tec in TECU unit with dimension of m1_occ tec = phase_to_tec(dph1, f1, f2) #print("tec:", tec.shape, tec.min(), tec.max(), tec) # ------------------------------------------------------------------ # Step 7. Estimation of the electron density at the top of the # occultation (orbit altitude) # and updated tec if icalib=0 for mode 0 # ------------------------------------------------------------------ tec, acoef, bcoef, edorb, hscale, topfit = process_tec(p_occ, tec, m1_occ, porb, calibration_mode, sigmatec=parms.sigmatec, sigmaeta=parms.sigmaeta, eta=parms.eta, fitlim=parms.fitlim) # ------------------------------------------------------------------ # Step 8. Inversion of TEC to electron density # ------------------------------------------------------------------ (m1_new, time_occ, h_occ, p_occ, tec, lat_occ, lon_occ, az_occ) = filter_by_height \ (horb, m1_occ, time_occ, h_occ, p_occ, tec, \ lat_occ, lon_occ, az_occ, dhmin=parms.dhmin) ed = invtec(p_occ, tec, porb, edorb) #print("ed:", ed.shape, ed.min(), ed.max(), ed) # ------------------------------------------------------------------ # Step 8.1 check the calibrated TEC results with UCAR # ------------------------------------------------------------------ plotfig_comp = 0 if plotfig_comp: import matplotlib.pyplot as plt nc_ucar = Dataset(ionPrf_file_UCAR) nc_ucar.set_auto_mask(False) MSL_alt_UCAR = nc_ucar["MSL_alt"][:] TEC_cal_UCAR = nc_ucar["TEC_cal"][:] ELEC_dens_UCAR = nc_ucar["ELEC_dens"][:] MSL_alt_occ = h_occ/1000.0 TEC_occ = tec #tec if MSL_alt_occ[0] > MSL_alt_occ[-1]: MSL_alt_occ = MSL_alt_occ[::-1] TEC_occ = TEC_occ[::-1] plt.figure(figsize=(5, 7)) plt.plot(TEC_cal_UCAR, MSL_alt_UCAR, linestyle='-', linewidth=2, label='UCAR') plt.plot(TEC_occ, MSL_alt_occ, '--', linewidth=2, label='STAR') plt.xlabel('TEC (TECU)') plt.ylabel('MSL altitude (km)') plt.title('TEC profile comparison') plt.grid(True) plt.legend() plt.tight_layout() plt.savefig('piq_TEC_vs_MSL_comparison_v6.png', dpi=300) plt.close() plt.figure(figsize=(5, 7)) plt.plot(ELEC_dens_UCAR, MSL_alt_UCAR, color='blue', label='UCAR') plt.plot(ed, MSL_alt_occ, color='red', label='STAR') #plt.plot(Ne_use * TECU/1e6, (r_use-roc)/1000.0, color='black', label='STAR Truncated') #plt.plot(Ne_use * TECU/1e6, (r_use-roc-undulation)/1000.0, color='black', label='STAR Truncated') plt.xlabel('ELEC Dens (el/cm3)') plt.ylabel('MSL altitude (km)') plt.title('ELEC Dens Profile Comparison') plt.grid(True) plt.legend() plt.tight_layout() plt.savefig('piq_ELEC_vs_MSL_comparison_v6.png', dpi=300) plt.close() TEC_this_on_UCAR = interp_to_ucar_height( MSL_alt_occ, TEC_occ, MSL_alt_UCAR) TEC_diff = TEC_this_on_UCAR - TEC_cal_UCAR mask = np.isfinite(TEC_diff) & np.isfinite(MSL_alt_UCAR) plt.figure(figsize=(5, 7)) plt.plot(TEC_diff[mask], MSL_alt_UCAR[mask], linewidth=2, label='STAR-UCAR') plt.axvline(0, color='k', linestyle='--', linewidth=1) plt.xlabel(r'$\Delta$TEC (TECU)') plt.ylabel('MSL altitude (km)') plt.title('TEC difference (STAR vs UCAR)') plt.grid(True) plt.legend() plt.tight_layout() plt.savefig('piq_TEC_difference_vs_MSL_v6.png', dpi=300) plt.close() # bottom/top bottime = float(int(startTime + time_occ[0])) toptime = float(int(startTime + time_occ[-1])) edorbtime= float(startTime + edorbtime_rel) edorbalt = float(horb/1000.0) # the geodesic distance between the first and last tangent-point locations, # projected onto the Earth's surface (meters). smear = horizontal_smear( lat_occ[0], lon_occ[0], lat_occ[-1], lon_occ[-1]) (edmax, edmaxalt, edmaxlat, edmaxlon, edmaxaz, edmaxtime, critfreq) = \ find_f_layer_peak(h_occ, ed, lat_occ, lon_occ, az_occ, time_occ) edmaxtime1 = float(int(startTime + edmaxtime )) # edmaxtime calendar time t0 = Time(startTime+edmaxtime, format='gps') t0_utc = t0.utc year = t0_utc.datetime.year month = t0_utc.datetime.month day = t0_utc.datetime.day hour = t0_utc.datetime.hour minute = t0_utc.datetime.minute second = float(int( t0_utc.datetime.second + t0_utc.datetime.microsecond * 1e-6)) print("Time:", year,month,day,hour,minute,second, bottime, toptime) # Electron density peak tec0 = vertical_tec(h_occ, ed) tec1 = extrapolate_tec(h_occ, p_occ, ed, edmaxalt) # local time (LCT) botlct = gps_to_lct(startTime + time_occ[0], lon_occ[0]) toplct = gps_to_lct(startTime + time_occ[-1], lon_occ[-1]) edmaxlct = gps_to_lct(startTime+edmaxtime, edmaxlon) #fileStamp = fileStamp.format(hour=hour, minute=minute) parts = fileStamp.split(".") parts[3] = f"{hour:02d}" parts[4] = f"{minute:02d}" fileStamp = ".".join(parts) print("fileStamp:", fileStamp) meta = { "occ_id" : np.int32(0), #"fiducial_id" : " " , "reference_sat_id": np.int32(-999) , "occulting_sat_id": np.int32(gnss_prn), "year" : np.int32(year), "month" : np.int32(month), "day" : np.int32(day), "hour" : np.int32(hour), "minute" : np.int32(minute), "second" : second, "offset" : np.int32(0), "shortlen": np.int32(m1_new), "setting" : np.int32(occ_mode), "icalib" : np.int32(calibration_mode), "botnum" : np.int32(1), "bottime" : bottime, "botlct" : float(botlct), "botalt" : float(h_occ[0]/1000.0), "botlat" : float(lat_occ[0]), "botlon" : float(lon_occ[0]), "botaz" : float(az_occ[0]), "topnum" : np.int32(m1_new), "toptime" : toptime, "toplct" : float(toplct), "topalt" : float(h_occ[-1]/1000.0), "toplat" : float(lat_occ[-1]), "toplon" : float(lon_occ[-1]), "topaz" : float(az_occ[-1]), "edmaxtime" : float(edmaxtime1), "edmaxlct" : float(edmaxlct), "edmaxalt" : float(edmaxalt/1000.0), "edmaxlat": float(edmaxlat), "edmaxlon" : float(edmaxlon), "edmaxaz" : float(edmaxaz), "smear" : float(smear/1000.0), # convert from meter to km "edmax" : float(edmax), "critfreq" : float(critfreq), "tec0" : float(tec0), "tec1" : float(tec1), "edorbtime": float(edorbtime), "edorbalt" : float(edorbalt), "edorb" : float(edorb), "hscale" : float(hscale), # only valid for icalib=0 "topfit" : float(topfit), # Post-fit residual near the top "nmax" : np.int32(m1_new), "fileStamp" : fileStamp, "inverter" : 'STAR_roion', #"parmsfile" : "center": "NOAA/STAR", "mission": "GNOMES" } nc.close() # Output filename: preserve the TEC filename and replace product token TEC -> EDP. # If no output directory is supplied, write beside the input TEC file. if ionPrf_file_out_dir is None: ionPrf_file_out_dir = os.path.dirname(os.path.abspath(ionPhs_file)) os.makedirs(ionPrf_file_out_dir, exist_ok=True) in_base = os.path.basename(ionPhs_file) if "_TEC_" not in in_base: raise ValueError(f"Input filename does not contain _TEC_: {in_base}") # Modified by Yong Chen on 09/29/2026 # to construct the filename based on the STAR filename convension bot_t = Time(bottime, format='gps') bot_t_utc = bot_t.utc bot_timestamp = bot_t_utc.datetime.strftime('%Y%m%d%H%M%S') print(bot_timestamp) top_t = Time(toptime, format='gps') top_t_utc = top_t.utc top_timestamp = top_t_utc.datetime.strftime('%Y%m%d%H%M%S') print(top_timestamp) start_timestamp = bot_timestamp end_timestamp = top_timestamp if bottime > toptime: start_timestamp = top_timestamp end_timestamp = bot_timestamp UTC = timezone.utc create_time = datetime.now(UTC) create_string = create_time.strftime("%Y%m%d%H%M%S") parts = in_base.split("_") edp_base = ( f"{parts[0]}_EDP_{parts[2]}_" f"{start_timestamp}_{end_timestamp}_" f"{parts[5]}_{parts[6]}_" f"c{create_string}_{parts[8]}_{parts[9]}" ) #edp_base = in_base.replace("_TEC_", "_EDP_", 1) ionPrf_file_out = os.path.join(ionPrf_file_out_dir, edp_base) logging.info(f"Output written to: {ionPrf_file_out}") write_ionPrf( ionPrf_file_out, h_occ/1000.0, lat_occ, lon_occ, az_occ, tec, ed, meta) # main program #run_gmrion('/data2/ychen131/TEC_EDP/cosmic/level1b/ionPhs/ionPhs_C001.2009.001.00.12.G26_2021.0390_nc', \ # '/data2/ychen131/TEC_EDP/cosmic/level2/ionPrf_C001.2009.001.00.12.G26_2021.0390_nc', \ # '/data2/ychen131/TEC_EDP/cosmic_STAR/level2/ionPrf_C001.2009.001.00.12.G26_2021.0390_nc') #run_gmrion('/data2/ychen131/TEC_EDP/cosmic/level1b/ionPhs/ionPhs_C004.2009.001.08.20.G12_2021.0390_nc', \ # '/data2/ychen131/TEC_EDP/cosmic/level2/ionPrf_C004.2009.001.08.20.G12_2021.0390_nc', \ # '/data2/ychen131/TEC_EDP/cosmic_STAR/level2/ionPrf_C004.2009.001.08.20.G12_2021.0390_nc') #ionPrf_C004.2009.001.08.20.G12_2021.0390_nc calibration mode=0 #ionPrf_C006.2009.001.23.45.G03_2021.0390_nc calibration mode=0 def process_tec_file(tec_file, out_dir=None): """Process one ROSW_TEC netCDF into ROSW_EDP using the original phase-based retrieval. exph_L1 and exph_L2 from the TEC file are passed through the unchanged ionPhs->ionPrf processing chain. The TEC variable itself is not substituted for the original phase input. """ return run_gmrion(tec_file, out_dir) # --------------------------------------------------------- # Batch Processing # --------------------------------------------------------- def process_batch(list_file): logging.info(f"Batch mode using file list: {list_file}") if not os.path.exists(list_file): raise FileNotFoundError(f"List file not found: {list_file}") with open(list_file) as f: lines = [l.strip() for l in f if l.strip()] success = 0 fail = 0 for line in lines: try: ionPhs_file, ionPrf_file_out_dir = line.split() run_gmrion( ionPhs_file, ionPrf_file_out_dir ) success += 1 except Exception as e: logging.error(f"Failed processing line: {line}") logging.error(str(e)) fail += 1 logging.info(f"Batch finished: success={success}, failed={fail}") # --------------------------------------------------------- # Argument Parsing # --------------------------------------------------------- def parse_args(): parser = argparse.ArgumentParser( description="GMR ionosphere EDP processing driver" ) parser.add_argument("--TEC", "--ionPhs", dest="tec_file", help="Input ROSW_TEC netCDF file containing exph_L1/exph_L2") parser.add_argument("--out", help="Output EDP directory (default: input file directory)") parser.add_argument("--list", help="Batch file list") parser.add_argument("--log", default=None, help="Optional log file") parser.add_argument("--loglevel", default="INFO", help="Logging level") return parser.parse_args() # --------------------------------------------------------- # Main Driver # --------------------------------------------------------- def main(): class Tee: def __init__(self, *files): self.files = files def write(self, obj): for f in self.files: f.write(obj) f.flush() def flush(self): for f in self.files: f.flush() logfile = open("TEC_to_EDP_output.txt", "w") sys.stdout = Tee(sys.stdout, logfile) args = parse_args() setup_logging(args.loglevel, args.log) UTC = timezone.utc start = datetime.now(UTC) logging.info("===== GMRION PROCESS START =====") try: if args.list: process_batch(args.list) else: if not args.tec_file: raise ValueError("Missing required --TEC input file") process_tec_file( args.tec_file, args.out ) except Exception as e: logging.exception("Fatal error occurred") sys.exit(1) end = datetime.now(UTC) logging.info(f"Processing completed in {(end-start).total_seconds():.1f} seconds") if __name__ == "__main__": main() # /data3/tcliu/gps_pod/python/planetiq/p05_podTec/2026-01-01/podTec_GN04.2026.001.20251231235955.20260101003411.A00.G03.nc # /data2/ychen131/TEC_EDP/planetiq_STAR/level2