from datetime import datetime, timedelta class DeferOccultation(Exception): """STAR 2026-11-05 (exph_skip_dup = 'defer'): the orbit after this occultation is not in the tars received so far (POD gap not yet fittable, or flag-2 stretch running past the end of the tar); step4 keeps the occultation and processes it again with the next tar of the same LEO.""" pass import numpy as np def log_message(msg): with open(f"log/progress_rt__{datetime.now():%Y%m%d}.txt", 'a') as fid: fid.write(f"{datetime.now():%H:%M:%S} {msg}\n") def to_datetime(dtlist): # TC: updated 2026-10-27 (1 line, fix the error for 0.9999995 sec case) return datetime(int(dtlist[0]),int(dtlist[1]),int(dtlist[2]),int(dtlist[3]),int(dtlist[4]))+timedelta(seconds=float(dtlist[5])) def dt2float(dt,dt0): return ((dt - dt0) / timedelta(seconds=1)).astype('float') def get_sat_freq(sat): GPRN = int(sat[1:]) if sat[0] == 'G': #GPS FL1 = 1.57542E9 FL2 = 1.22760E9 elif sat[0] == 'R': #GLONASS Kfrq = [1, -4, 5, 6, 1, -4, 5, 6, -2, -7, 0, -1, -2, -7, 0, -1, 4, -3, 3, 2, 4, -3, 3, 2] baseFL1 = 1.602E9 baseFL2 = 1.246E9 DFRQ1_GLO = 5.625E5 DFRQ2_GLO = 4.375E5 FL1 = baseFL1 + Kfrq[GPRN - 1] * DFRQ1_GLO FL2 = baseFL2 + Kfrq[GPRN - 1] * DFRQ2_GLO GPRN = GPRN + 100 elif sat[0] == 'E': #GALILEO FL1 = 1.57542E9 FL2 = 1.20714E9 GPRN = GPRN + 200 elif sat[0] == 'C': # BEIDOU FL1 = 1.57542E9 FL2 = 1.17645E9 GPRN = GPRN + 400 else: # ERROR!! FL1 = 0 FL2 = 0 GPRN = 0 return FL1, FL2, GPRN def load_gnss(file): with np.load(file, allow_pickle=True) as fid: Gdata = fid['Gdata'] Gdt = fid['Gdt'] Gsats = fid['Gsats'] return Gdata, Gdt, Gsats def read_gnss(dt,prev=False): import os Gdata1 = np.array([]) if prev: for x in ['NRT','RT']: gnssfile = f'../gnss/gnssOrb/GNSS_{dt - timedelta(days=1):%Y.%j}_{x}.npz' if os.path.isfile(gnssfile): if len(Gdata1) == 0: Gdata1, Gdt1, Gsats1 = load_gnss(gnssfile) else: Gdata0, Gdt0, Gsats0 = load_gnss(gnssfile) if Gdt1.shape == Gdt0.shape and np.all(Gdt1 == Gdt0): idx = ~np.isin(Gsats0, Gsats1) Gsats1 = np.hstack([Gsats1, Gsats0[idx]]) Gdata1 = np.hstack([Gdata1, Gdata0[:, idx, :]]) del Gdata0, Gdt0, Gsats0 Gdata2 = np.array([]) for x in ['NRT','RT']: gnssfile = f'../gnss/gnssOrb/GNSS_{dt:%Y.%j}_{x}.npz' if os.path.isfile(gnssfile): if len(Gdata2) == 0: Gdata2, Gdt2, Gsats2 = load_gnss(gnssfile) else: Gdata0, Gdt0, Gsats0 = load_gnss(gnssfile) if Gdt2.shape == Gdt0.shape and np.all(Gdt2 == Gdt0): idx = ~np.isin(Gsats0, Gsats2) Gsats2 = np.hstack([Gsats2, Gsats0[idx]]) Gdata2 = np.hstack([Gdata2, Gdata0[:, idx, :]]) del Gdata0, Gdt0, Gsats0 Gsats = Gsats2 if prev: idx = Gdt1dt_rng[1]) if np.any(tidx): gfiles = gfiles[tidx] tstart = tstart[tidx] tidx = tstart == max(tstart) gfiles = gfiles[tidx] tstart = tstart[tidx] # tend = tend[tidx] Gdt = None Gsats = [] Gdata = [] for ss in src_types: sidx = np.where([ss in gf for gf in gfiles])[0] if sidx.size == 0: continue gfile = gfiles[sidx[0]] Gdata0, Gdt0, Gsats0 = load_gnss(gfile) tidx = np.logical_and(Gdt0>=dt_rng[0]-timedelta(minutes=10),Gdt0 dt_rng[0]) gfiles = gfiles[tidx] tstart = tstart[tidx] Gdt = None Gsats = [] Gdata = [] for gfile in gfiles: Gdata0, Gdt0, Gsats0 = load_gnss(gfile) tidx = np.logical_and(Gdt0>=dt_rng[0]-timedelta(minutes=10),Gdt0dt_rng[1]) gfiles = gfiles[tidx] tstart = tstart[tidx] # tend = tend[tidx] Gdt = None for gsat in gsats: gfile = None for ss in src_types: sidx = np.where([ss in gf for gf in gfiles])[0] if sidx.size==0: continue sfile = gfiles[sidx[np.argmin(np.abs(dt_rng[0] - tstart[sidx] - timedelta(hours=6)))]] with np.load(sfile, allow_pickle=True) as fid: if gsat in fid['Gsats']: gfile = sfile break if gfile is None: raise FileNotFoundError(f'Could not find a GNSS Orb file containing {gsat}.') Gdata0, Gdt0, Gsats0 = load_gnss(gfile) tidx = np.logical_and(Gdt0>=dt_rng[0]-timedelta(minutes=10),Gdt0 nX) idxL = max(min(idx - order, nX - 2 * order), 0) idxR = min(max(idx + order, 2 * order), nX) J = np.arange(idxL, idxR) pX = x[ii] - X[J] # method 1, quick as it multiplies all numerators and all denominators before division pX = np.prod(pX) / pX R = np.prod(T[J][:,J], axis=1) y[ii] = (pX / R) @ Y[J] # pX = np.tile(x[ii] - X[J], (len(J), 1)) # method 2, reduce error as it pairs numerators and denominators # np.fill_diagonal(pX, 1) # R = T[J][:,J] # y[ii] = np.prod(pX / R, axis=1) @ Y[J] return y def inv_quat_rotate(q, v): t = np.array([q[:,3]*v[0]+q[:,1]*v[2]-q[:,2]*v[1], q[:,3]*v[1]+q[:,2]*v[0]-q[:,0]*v[2], q[:,3]*v[2]+q[:,0]*v[1]-q[:,1]*v[0], -q[:,0]*v[0]-q[:,1]*v[1]-q[:,2]*v[2]]) v = np.array([-t[3]*q[:,0]+t[0]*q[:,3]-t[1]*q[:,2]+t[2]*q[:,1], -t[3]*q[:,1]+t[1]*q[:,3]-t[2]*q[:,0]+t[0]*q[:,2], -t[3]*q[:,2]+t[2]*q[:,3]-t[0]*q[:,1]+t[1]*q[:,0]]).transpose() return v def pod_quality(sec, pos_km, vel_dms, t0, t1, gap_min=60.0, pv_win=30.0): """STAR 2026-11-01: quality of the 1 s LEO POD orbit around an occultation [t0, t1] (s, same scale as sec). Returns dict: gap_after_s : seconds from the end of the nearest preceding POD gap (> gap_min s) to the occultation start (time since the orbit resumed; inf if no such gap before the occultation in the file) gap_before_s : seconds from the occultation end to the start of the next POD gap (inf if none) pv_max : max over [t0 - pv_win, t1 + pv_win] of the pv_win-s mean of |dP/dt - V| (m/s); dP/dt by centred differences within continuous 1 s segments. Normally ~0.002 m/s; after a POD outage the PlanetiQ orbit converges over a few minutes (0.31 m/s in the first 30 s after the 13.5 min gap on 8 Oct).""" sec = np.asarray(sec, float) P = np.asarray(pos_km, float)[:, :3] * 1e3; V = np.asarray(vel_dms, float)[:, :3] * 0.1 # m, m/s (dm/s in file) gi = np.where(np.diff(sec) > gap_min)[0] ends = sec[gi + 1]; starts = sec[gi] # gap end (data resumes) / gap start (data stops) a = ends[ends <= t0]; b = starts[starts >= t1] out = dict(gap_after_s=float(t0 - a.max()) if a.size else np.inf, gap_before_s=float(b.min() - t1) if b.size else np.inf) k = (sec >= t0 - pv_win - 1) & (sec <= t1 + pv_win + 1) pv = 0.0 if k.sum() > 3: ii = np.where(k)[0] ok = (np.abs(sec[ii[2:]] - sec[ii[:-2]] - 2.0) < 1e-3) # centred differences inside a 1 s segment tc = sec[ii[1:-1]][ok] r = ((P[ii[2:]] - P[ii[:-2]]) / 2.0 - V[ii[1:-1]])[ok] for c0 in np.arange(t0 - pv_win, t1 + 1e-6, pv_win / 2): q = (tc >= c0) & (tc < c0 + pv_win) if q.sum() >= pv_win / 2: pv = max(pv, float(np.linalg.norm(r[q].mean(0)))) out['pv_max'] = pv return out def smooth_1s(sec, pos, vel, hw=15): """STAR 2026-10-31 (exph_leo_interp = 'smooth'): local least-squares cubic over +-hw s (Savitzky-Golay) of the 1 s LEO positions (x, y, z; clock column unchanged) and velocities, applied to each continuous 1 s segment.""" from scipy.signal import savgol_filter pos = pos.copy(); vel = vel.copy() brk = np.where(np.abs(np.diff(sec) - 1.0) > 1e-3)[0] + 1 for seg in np.split(np.arange(len(sec)), brk): if len(seg) < 2 * hw + 1: continue pos[seg, :3] = savgol_filter(pos[seg, :3], 2 * hw + 1, 3, axis=0) vel[seg, :3] = savgol_filter(vel[seg, :3], 2 * hw + 1, 3, axis=0) return pos, vel def phase_correction(l1, pos_gnss, vel_gnss, sec_gnss, pos_leo, vel_leo, att_leo, sec_leo, sec_ro, antoffset_leo, calc_exph, leo_order=11): C = 299792458 G = 6.67408E-11 Mearth = 5.9722E24 ocd = 0 const = -6.9693E-10 GMoCC = 2*G*Mearth/C/C # Schwarzschild radius pos_gnss_intp = lagrange(sec_gnss - pos_gnss[:, 3] / 1E6, pos_gnss[:, :3], sec_ro, 11) # vel_gnss_intp = lagrange(sec_gnss - pos_gnss[:, 3] / 1E6, vel_gnss[:, :3], sec_ro, 11) pos_leo_intp = lagrange(sec_leo, pos_leo[:, :3], sec_ro, leo_order) # STAR 2026-10-31: order 4 (8 nodes) on 1 s nodes vel_leo_intp = lagrange(sec_leo, vel_leo[:, :3], sec_ro, leo_order) # for exph_leo_interp = 'full1s' / 'smooth' att_leo_intp = lagrange(sec_leo, att_leo, sec_ro, 4) att_leo_intp /= np.sum(att_leo_intp ** 2, axis=1)[:,None] ** 0.5 antoffset_leo_intp = inv_quat_rotate(att_leo_intp, antoffset_leo) pos_leo_intp += antoffset_leo_intp antvel_leo_intp = np.full(antoffset_leo_intp.shape, np.nan) antvel_leo_intp[1:-1] = (antoffset_leo_intp[2:] - antoffset_leo_intp[:-2]) / (sec_ro[2:,None] - sec_ro[:-2,None]) antvel_leo_intp[0] = antvel_leo_intp[1] antvel_leo_intp[-1] = antvel_leo_intp[-2] vel_leo_intp += antvel_leo_intp obj_leo = dict() obj_leo['pos'] = pos_leo_intp obj_leo['vel'] = vel_leo_intp obj_leo['clk_s'] = lagrange(sec_leo, pos_leo[:, 3], sec_ro, 6) / 1E6 obj_leo['clk_m'] = obj_leo['clk_s'] * C obj_gnss = dict() obj_gnss['clk_s'] = lagrange(sec_gnss, pos_gnss[:, 3], sec_ro, 6) / 1E6 obj_gnss['clk_m'] = obj_gnss['clk_s'] * C tao_new = np.array([1E5]) for _ in range(50): range_lg = np.sum((pos_leo_intp - pos_gnss_intp) ** 2, axis=1) ** 0.5 if calc_exph: exph = l1 - range_lg vec_ol = pos_leo_intp[[0,-1]] vec_og = pos_gnss_intp[[0,-1]] vec_lg = vec_og - vec_ol proj_lg2 = np.sum(np.cross(vec_ol, vec_og) ** 2, axis=1) / np.sum(vec_lg ** 2, axis=1) ocd = proj_lg2[0] > proj_lg2[1] exph -= exph[0 if ocd else -1] else: exph = np.zeros_like(l1) range_leo = np.sum(pos_leo_intp ** 2, axis=1) ** 0.5 range_gnss = np.sum(pos_gnss_intp ** 2, axis=1) ** 0.5 delta_rho = GMoCC * np.log((range_gnss + range_leo + range_lg) / (range_gnss + range_leo - range_lg)) tao_old = tao_new.copy() tao_new = (range_lg + exph + delta_rho) / C pos_gnss_intp = lagrange(sec_gnss - pos_gnss[:, 3] / 1E6, pos_gnss[:, :3], sec_ro - tao_new, 9) if not np.any(np.abs(tao_new - tao_old) > 1E-12): break vel_gnss_intp = lagrange(sec_gnss - pos_gnss[:, 3] / 1E6, vel_gnss[:, :3], sec_ro - tao_new, 9) range_dT = np.sum((pos_leo_intp - pos_gnss_intp) ** 2, axis=1) ** 0.5 # from scipy.interpolate import make_interp_spline # speed_dT = make_interp_spline(sec_ro, range_dT, k=2).derivative()(sec_ro) obj_gnss['pos'] = pos_gnss_intp obj_gnss['vel'] = vel_gnss_intp coef_leo = 1 - (GMoCC / range_leo - np.sum(vel_leo_intp ** 2, axis=1) / C / C) / 2 - const propertimedelay_gnss = (coef_leo - 1) * obj_gnss['clk_s'] * C delt_t_rel_gnss = -2 * np.sum(pos_gnss_intp * vel_gnss_intp, axis=1) / C / C delt_t_rel_leo = -2 * np.sum(pos_leo_intp * vel_leo_intp, axis=1) / C / C obj_gnss['range_diff_T'] = range_dT # obj_gnss['range'] = range_lg obj_gnss['GL_term'] = -delta_rho + C * delt_t_rel_gnss + propertimedelay_gnss + obj_gnss['clk_s'] * C obj_leo['delt_t_rel_r'] = delt_t_rel_leo obj_gnss['delt_t_rel_s'] = delt_t_rel_gnss obj_leo['ocd'] = ocd return obj_gnss, obj_leo def remove_cycleslip(ph, ph_m, freq, sat_l, ocd): c0 = 299792458 if not ocd: ph = ph[::-1] ph_m = ph_m[::-1] ph_new = (ph - ph_m) * (2 * freq / c0) if sat_l == 'G1': idx = np.unique(np.where(np.round(np.diff(ph_new[:100])) % 2 == 1)[0] % 2) if idx.size != 1: ph_new -= np.cumsum(np.concatenate(([0], np.round(np.diff(ph_new))))) else: for ii in range(idx[0] + 1, len(ph_new), 2): ph_new[ii:ii + 2] -= np.round(ph_new[ii] - ph_new[ii - 1]) ph_new -= np.cumsum(np.concatenate(([0], np.round(np.diff(ph_new) / 2) * 2))) elif sat_l[0] == 'R': ph_new -= np.cumsum(np.concatenate(([0], np.round(np.diff(ph_new))))) else: ph_new -= np.cumsum(np.concatenate(([0], np.round(np.diff(ph_new) / 2) * 2))) ph_new = ph_new * (c0 / (2 * freq)) + ph_m if not ocd: ph_new = ph_new[::-1] return ph_new # TC: updated 2026-10-27 (added 93 lines, new function to account for navigation bits) def remove_cycleslip_with_navbit(ph, ph_m, freq, gnss_ph, ocd, dtr_ph, sec_ph): c0 = 299792458 ph_new = (ph - ph_m) * (2 * freq / c0) with_navbit = False if gnss_ph[0] in 'GR': try: rate = 100 smp = 1 vld = 170 if gnss_ph[0]=='R' else 200 match_len = 1000 m_s = 0 if ocd else ph_new.shape[0] - match_len m_e = m_s + match_len h1, h2 = int((sec_ph[0] - 1) // 3600), int(sec_ph[-1] // 3600) bad = np.abs((np.diff(ph_new[m_s:m_e, 0], prepend=ph_new[m_s:m_s+1, 0]) + 0.5) % 1 - 0.5) > 0.15 nbit_ref = (np.cumsum(np.round(np.diff(ph_new[m_s:m_e, 0], prepend=ph_new[m_s:m_s+1, 0]))) % 2).astype('int8') nbit_sec = sec_ph[m_s] bad = np.convolve(bad, np.ones(11,), mode='same')>0 glen = np.sum(~bad) if glen<500: raise ValueError('Go to except') gnss_dir = '../gnss/gnssNavBit/' with np.load(gnss_dir + f'{dtr_ph+timedelta(hours=h1):%Y.%j}/navBit_{gnss_ph}_{dtr_ph+timedelta(hours=h1):%Y.%j.%H}.npz') as fid: navbit = np.unpackbits(fid['navbit']) navbitmask = np.unpackbits(fid['navbitmask']) if h2>h1: with np.load(gnss_dir + f'{dtr_ph+timedelta(hours=h2):%Y.%j}/navBit_{gnss_ph}_{dtr_ph+timedelta(hours=h2):%Y.%j.%H}.npz') as fid: navbit = np.hstack([navbit,np.unpackbits(fid['navbit'])]) navbitmask = np.hstack([navbitmask,np.unpackbits(fid['navbitmask'])]) if gnss_ph[0] == 'G': navbit = np.repeat(navbit, 2) navbitmask = np.repeat(navbitmask, 2) else: navbit = np.column_stack([navbit, 1 - navbit]).ravel() navbit = navbit.reshape((-1, 200)) navbit[:, 170:] = np.array([int(t) for t in '111110001101110101000010010110'], dtype='uint8') navbit = navbit.ravel() navbitmask = np.repeat(navbitmask, 2) navbit[navbitmask == 0] = 0 idx0 = int((nbit_sec - h1 * 3600) * rate) nbit_check = np.full((20,2), np.nan) for ii in range(20): if not np.all(navbitmask[idx0 - smp*ii:idx0 - smp*ii + smp*match_len:smp]): continue idxn = np.arange(idx0 - smp * ii, idx0 - smp * ii + smp * match_len, smp) % 200 >= vld nbit = navbit[idx0 - smp*ii:idx0 - smp*ii + smp*match_len:smp].copy() for jj in [0,1]: if jj == 1: nbit[idxn] = 1 - nbit[idxn] temp = (nbit + nbit_ref) % 2 nbit_check[ii,jj] = np.min([np.sum(temp[~bad]==0),np.sum(temp[~bad]==1)])/glen ii = np.nanargmin(np.min(nbit_check,axis=1)) # print(f'{ii:d}, {nbit_check[ii][0]:.2f}, {nbit_check[ii][1]:.2f}') if min(nbit_check[ii])<0.01: idx0 = int((nbit_sec - h1 * 3600) * rate) - m_s idx = np.arange(idx0 - smp * ii, idx0 - smp * ii + smp * ph_new.shape[0], smp) if not np.all(navbitmask[idx]): raise ValueError('Go to except') if nbit_check[ii][0] > nbit_check[ii][1]: navbit[idx % 200 > vld] = 1 - navbit[idx % 200 > vld] ph_new[:,0] += navbit[idx] if gnss_ph[0] == 'R': ph_new[:, 1] += navbit[idx] ph_new -= np.cumsum(np.round(np.diff(ph_new, axis=0, prepend=ph_new[:1]) / 2) * 2, axis=0) else: raise ValueError('Go to except') with_navbit = True except: if gnss_ph[0] == 'G': match_len = 100 m_s = 0 if ocd else ph_new.shape[0] - match_len m_e = m_s + match_len idx = np.unique(np.where(np.round(np.diff(ph_new[m_s:m_e, 0])) % 2 == 1)[0] % 2) if idx.size != 1: ph_new[:, 0] -= np.cumsum(np.round(np.diff(ph_new[:, 0], axis=0, prepend=ph_new[:1, 0])), axis=0) else: for ii in range(idx[0] + 1, len(ph_new), 2): ph_new[ii:ii + 2, 0] -= np.round(ph_new[ii, 0] - ph_new[ii - 1, 0]) ph_new[:, 0] -= np.cumsum(np.round(np.diff(ph_new[:, 0], prepend=ph_new[:1, 0]) / 2) * 2) ph_new[:, 1] -= np.cumsum(np.round(np.diff(ph_new[:, 1], axis=0, prepend=ph_new[:1, 1]) / 2) * 2, axis=0) if gnss_ph[0] == 'R': for ii in range(1, ph_new.shape[0]): ph_diff = ph_new[ii] - ph_new[ii - 1] if np.all(np.round(ph_diff) % 2 == 1) or np.all( (*(np.round(ph_diff) % 2 == [1, 0]), *(np.diff(np.abs(ph_diff - np.round(ph_diff))) > 0))) or np.all( (*(np.round(ph_diff) % 2 == [0, 1]), *(np.diff(np.abs(ph_diff - np.round(ph_diff))) < 0))): ph_new[ii] += 1 ph_new -= np.cumsum(np.round(np.diff(ph_new, axis=0, prepend=ph_new[:1]) / 2) * 2, axis=0) else: ph_new -= np.cumsum(np.round(np.diff(ph_new, axis=0, prepend=ph_new[:1]) / 2) * 2, axis=0) if not ocd: ph_new -= np.round(ph_new[-1] / 2) * 2 ph_new = ph_new * (c0 / (2 * freq)) + ph_m return with_navbit, ph_new.T # ====================================================================== # STAR 2026-10-30: excess-phase processing options (see README_exph_options.md) # exph_clock : 'single_diff' (legacy) LEO clock removed with the POD reference link # 'zero_diff' (proposed) LEO clock + relativity from the POD orbit file # exph_clk_interp (zero_diff only): # 'linear' (default) linear interpolation of the full-rate (1 s) POD clock, # as CDAAC conPhs ("clkinterp=linear") # 'lagrange' 6th-order Lagrange on the decimated (<= 20 point) orbit grid # (first test, 2026-10-30) # 'poly' polynomial fit to the 1 s POD clock over exph_clk_poly_dur s centred # on the occultation (degree exph_clk_poly_deg), as CDAAC genXphase52 # GTLCL0 for zero differencing (idur = 300, ndeg = 2) # exph_connect : 'star' (legacy) nav bits + STAR half-cycle rules + Fortran cyclesliperemoval # 'starc' STARC: STAR implementation of the CDAAC newroam rawphs.f algorithm: down-conversion with the # CIRA+Q climatological Doppler model (ray1/dop1/BENMOD, Python port in # planetiq/starc_dopmodel.py), nav bits # (STAR bit files) or phasor-flip rule, phase connection; no Fortran step # exph_cycleslip : (exph_connect = 'star' only) 'python' (default) port of the cyclesliperemoval Fortran # (planetiq/ropp_openloop.py, MSIS model phase from exph_msis_file), or 'fortran' (PhaseModel/cyclesliperemoval) # exph_cycleslip_pilot : (python only) True = no internal nav-bit search on dataless pilot L1 signals # (Galileo E1-C, BeiDou B1C-P); False (default) = as the Fortran (20 ms bits on every 1575.42 MHz signal) # exph_starc_nobit : (starc) GPS/GLONASS without a bit-file match: 'star' (default) STAR fallback on the # residual to the CIRA+Q model, or 'flip' (CDAAC phasor-flip rule) # exph_connect_gnss : constellations that use exph_connect = 'starc' (default 'GREC'); others keep 'star' # exph_leo_interp : LEO orbit/attitude interpolation. Default: 'smooth' with zero_diff, 'decimated' with single_diff. # 'decimated' (delivered) 1 s POD orbit decimated to >= 20 nodes (~30 s), 22-node Lagrange: result depends # on the window; 'full1s' all 1 s epochs as nodes, 8-node Lagrange: passes the 1 mm rounding of the POD # positions into the phase (+0.27 mm at 20-60 km vs CDAAC); 'smooth' as full1s after a local least-squares # cubic over +-exph_leo_smooth_hw (15) s: window independent, same agreement with CDAAC as 'decimated' # exph_out_date : output date folder of s04_atmPhs and of the in-line ROPP/BUFR outputs: 'observation' (default, # delivered: day of the occultation) or 'input' (name of the input folder YYYY-MM-DD the tar came from, # i.e. the arrival date; needs a driver that passes it: process_excess_phase_rt.py does) # exph_out_subdir : sub-folder of s04_atmPhs (and of the in-line ROPP outputs) before the date folder ('' default: # s04_atmPhs/); set per model by apply_exph_model when reprocessing several models (/) # exph_inline_ropp : True (default, delivered) runs the in-line L1A -> L1C and BUFR when ropp_path exists; False skips it # exph_ropp_cf, exph_ropp_timeout : in-line L1A -> L1C (when ropp_path exists): ropp_pp_occ_tool configuration file in the # working directory and timeout in s (default 'cosmic_pp_rfsi_planetiq.cf', 20 s as delivered) # exph_ropp_pp : ropp_pp folder of the in-line tool under ropp-10.0_ ('ropp_pp' default; 'ropp_pp_y' = Y_PM_GO) # exph_pod_qc : LEO POD quality check: 'flag' (default) writes pod_qc / pod_gap_after_s / pod_gap_before_s / pod_pv_max # and sets the L1A quality flag bad = "1" (errstr = "POD ") when the POD is still bad after # the fixes; 'reject' (STAR 2026-11-04) does not write those occultations; 'off' as 'flag'. # pod_qc is the record of what was found and done; bad = "1" says the POD is still bad: # gap_unfitted (bit 1 without bit 4), pv_after (|dP/dt - V| on the orbit used still > exph_pod_pv_max), # fit_rms (a fit residual > exph_pod_fit_rms_max, 0.5 m). # pod_qc bit 1: the orbit resumed after a gap (> 60 s) less than exph_pod_gap_win (120) s before the # occultation; bit 2: max 30 s mean |dP/dt - V| of the 1 s POD > exph_pod_pv_max (0.02 m/s; normal ~0.002) # (8 Oct: all occultations damaged by the POD transient started <= 88 s after a gap, pv >= 0.024 m/s; # those from 129 s on, pv <= 0.014 m/s, agree with CDAAC as elsewhere) # exph_pod_fit : STAR 2026-11-02: seconds after a POD gap (> 60 s) whose orbit is replaced by a dynamic fit (0 / missing = # off, as delivered; 300 recommended). The fit (planetiq/pod_fit.py: EGM2008 degree 40 via satkit, 3 empirical # accelerations) uses the POD positions of [gap end + exph_pod_fit, + exph_pod_fit_len (600 s)] and is # integrated backward; it needs exph_pod_fit_min_after (900 s) of POD after the gap in the orbit merged up to # the current tar, otherwise only the flag is set. The POD clock is kept. Replaced orbits set pod_qc bit 3 # (value 4) and attributes pod_fit ('fit', 'nodata', 'nosatkit', 'failed', 'none', 'off'), pod_fit_rms_m, # pod_fit_maxdiff_m. One fit per gap, cached in s01_leoOrb/podfit/. # exph_pod_fit_pv : STAR 2026-11-04: True = DMF for flag-2 stretches away from gaps (POD velocity inconsistent with the # positions). The stretch (30 s mean |dP/dt - V| > 0.005 m/s, joined over < 120 s, peak > exph_pod_pv_max, # +- exph_pod_fit_pv_margin 30 s) is replaced by one dynamic orbit fitted to 600 s of POD before AND after # it in the orbit merged up to the current tar (exph_pod_fit_pv_mode 'bridge'), blended into the POD over # exph_pod_fit_pv_taper (120 s) at both edges. Sets pod_qc bit value 8, pod_method = DMF, attributes # pod_fit_pv ('bridge', 'nodata' = no POD after the stretch yet, 'gap', 'long', 'none', 'off', ...), # pod_fit_pv_span_s, _rms_m, _maxdiff_m and pod_pv_max_after. Test 8 Oct, 47 occultations: RMS O-B # 40-50 km 10.3 -> 6.3 % (CDAAC 9.2 %, unflagged 4.0 %). # exph_zd_noref : (zero_diff only) True (default) = no reference-link search: the occultation is processed # without a paired POD reference profile and the orbit windows are taken from the occultation # (+- exph_zd_win, default 300 s); # False = as legacy (a reference profile is still required, though not used) # Settings come from configuration.py; missing keys give the legacy behaviour. # ====================================================================== # STAR 2026-11-01: named settings ("models") for reprocessing several settings in one run (process_excess_phase_rt.py # --models); each sets these configuration keys, all other exph_* keys stay as in configuration.py EXPH_MODELS = { 'STAR_Legacy': dict(exph_clock='single_diff', exph_connect='star'), 'STAR_Updated': dict(exph_clock='zero_diff', exph_clk_interp='linear', exph_connect='star'), 'STARC': dict(exph_clock='zero_diff', exph_clk_interp='linear', exph_connect='starc'), } def apply_exph_model(name, subdir=None): """Set the configuration keys of a named setting in the loaded configuration module (before cal_exph, so the worker processes inherit them); outputs go to s04_atmPhs// (subdir defaults to the name).""" import configuration as cfg if name not in EXPH_MODELS: raise ValueError(f'unknown model {name}; known: {", ".join(EXPH_MODELS)}') for k, v in EXPH_MODELS[name].items(): setattr(cfg, k, v) cfg.exph_out_subdir = name if subdir is None else subdir def exph_options(): import configuration as cfg opt = dict(clock=getattr(cfg, 'exph_clock', 'single_diff'), clk_interp=getattr(cfg, 'exph_clk_interp', 'linear'), clk_poly_dur=getattr(cfg, 'exph_clk_poly_dur', 300.0), clk_poly_deg=getattr(cfg, 'exph_clk_poly_deg', 2), zd_noref=getattr(cfg, 'exph_zd_noref', True), zd_win=float(getattr(cfg, 'exph_zd_win', 300.0)), leo_interp=getattr(cfg, 'exph_leo_interp', None), pod_qc=getattr(cfg, 'exph_pod_qc', 'flag'), out_date=getattr(cfg, 'exph_out_date', 'observation'), out_subdir=str(getattr(cfg, 'exph_out_subdir', '') or ''), inline_ropp=bool(getattr(cfg, 'exph_inline_ropp', True)), ropp_cf=getattr(cfg, 'exph_ropp_cf', 'cosmic_pp_rfsi_planetiq.cf'), ropp_timeout=float(getattr(cfg, 'exph_ropp_timeout', 20)), ropp_pp=str(getattr(cfg, 'exph_ropp_pp', 'ropp_pp') or 'ropp_pp'), pod_gap_win=float(getattr(cfg, 'exph_pod_gap_win', 120.0)), pod_pv_max=float(getattr(cfg, 'exph_pod_pv_max', 0.02)), pod_fit=float(getattr(cfg, 'exph_pod_fit', 0.0) or 0.0), pod_fit_len=float(getattr(cfg, 'exph_pod_fit_len', 600.0)), pod_fit_min_after=float(getattr(cfg, 'exph_pod_fit_min_after', 900.0)), pod_fit_pv=bool(getattr(cfg, 'exph_pod_fit_pv', False)), pod_fit_pv_margin=float(getattr(cfg, 'exph_pod_fit_pv_margin', 30.0)), pod_fit_pv_side_max=float(getattr(cfg, 'exph_pod_fit_pv_side_max', 420.0)), pod_fit_pv_taper=float(getattr(cfg, 'exph_pod_fit_pv_taper', 120.0)), pod_fit_pv_mode=str(getattr(cfg, 'exph_pod_fit_pv_mode', 'bridge')), pod_fit_rms_max=float(getattr(cfg, 'exph_pod_fit_rms_max', 0.5)), skip_dup=getattr(cfg, 'exph_skip_dup', True), defer_now=bool(getattr(cfg, 'exph_defer_now', False)), defer_max_h=float(getattr(cfg, 'exph_defer_max_h', 3.0)), leo_smooth_hw=int(getattr(cfg, 'exph_leo_smooth_hw', 15)), connect=getattr(cfg, 'exph_connect', 'star'), connect_gnss=getattr(cfg, 'exph_connect_gnss', 'GREC'), starc_nobit=getattr(cfg, 'exph_starc_nobit', 'star'), cycleslip=getattr(cfg, 'exph_cycleslip', 'python'), cycleslip_pilot=getattr(cfg, 'exph_cycleslip_pilot', False), msis_file=getattr(cfg, 'exph_msis_file', 'PhaseModel/MSIS_coeff.nc'), benmod=getattr(cfg, 'starc_benmod', 'PhaseModel/BENMOD')) assert opt['clock'] in ('single_diff', 'zero_diff'), opt['clock'] if opt['leo_interp'] is None: # STAR 2026-11-01: 'smooth' for zero difference, delivered 'decimated' for single difference opt['leo_interp'] = 'smooth' if opt['clock'] == 'zero_diff' else 'decimated' assert opt['leo_interp'] in ('decimated', 'full1s', 'smooth'), opt['leo_interp'] assert opt['pod_qc'] in ('flag', 'reject', 'off'), opt['pod_qc'] # STAR 2026-11-05: exph_skip_dup True = 'first' (keep the copy written first), 'replace_bad' (default in configuration.py: # a re-sent copy replaces a first copy whose L1A bad = "1"), 'replace_flagged' (also when the first copy has pod_qc > 0), # False = 'off' (delivered check, which never matched) sd = opt['skip_dup'] opt['skip_dup'] = {True: 'first', False: 'off', None: 'off'}.get(sd, sd) if not isinstance(sd, str) else sd.lower() assert opt['skip_dup'] in ('first', 'replace_bad', 'replace_flagged', 'defer', 'off'), opt['skip_dup'] assert opt['out_date'] in ('observation', 'input'), opt['out_date'] assert opt['clk_interp'] in ('linear', 'lagrange', 'poly'), opt['clk_interp'] if opt['connect'] == 'cdaac': opt['connect'] = 'starc' # earlier name assert opt['connect'] in ('star', 'starc'), opt['connect'] assert opt['cycleslip'] in ('python', 'fortran'), opt['cycleslip'] assert opt['starc_nobit'] in ('star', 'flip'), opt['starc_nobit'] return opt def navbit_aligned(ph, ph_m, freq, gnss_ph, ocd, dtr_ph, sec_ph): """Nav-bit sequence aligned to the samples (0/1 per sample, in half cycles), using the same bit files and matching as remove_cycleslip_with_navbit; None when unavailable.""" c0 = 299792458 if gnss_ph[0] not in 'GR': return None ph_new = (ph - ph_m) * (2 * freq / c0) try: rate = 100 smp = 1 vld = 170 if gnss_ph[0] == 'R' else 200 match_len = 1000 m_s = 0 if ocd else ph_new.shape[0] - match_len m_e = m_s + match_len h1, h2 = int((sec_ph[0] - 1) // 3600), int(sec_ph[-1] // 3600) bad = np.abs((np.diff(ph_new[m_s:m_e, 0], prepend=ph_new[m_s:m_s+1, 0]) + 0.5) % 1 - 0.5) > 0.15 nbit_ref = (np.cumsum(np.round(np.diff(ph_new[m_s:m_e, 0], prepend=ph_new[m_s:m_s+1, 0]))) % 2).astype('int8') nbit_sec = sec_ph[m_s] bad = np.convolve(bad, np.ones(11,), mode='same') > 0 glen = np.sum(~bad) if glen < 500: return None gnss_dir = '../gnss/gnssNavBit/' with np.load(gnss_dir + f'{dtr_ph+timedelta(hours=h1):%Y.%j}/navBit_{gnss_ph}_{dtr_ph+timedelta(hours=h1):%Y.%j.%H}.npz') as fid: navbit = np.unpackbits(fid['navbit']) navbitmask = np.unpackbits(fid['navbitmask']) if h2 > h1: with np.load(gnss_dir + f'{dtr_ph+timedelta(hours=h2):%Y.%j}/navBit_{gnss_ph}_{dtr_ph+timedelta(hours=h2):%Y.%j.%H}.npz') as fid: navbit = np.hstack([navbit, np.unpackbits(fid['navbit'])]) navbitmask = np.hstack([navbitmask, np.unpackbits(fid['navbitmask'])]) if gnss_ph[0] == 'G': navbit = np.repeat(navbit, 2) navbitmask = np.repeat(navbitmask, 2) else: navbit = np.column_stack([navbit, 1 - navbit]).ravel() navbit = navbit.reshape((-1, 200)) navbit[:, 170:] = np.array([int(t) for t in '111110001101110101000010010110'], dtype='uint8') navbit = navbit.ravel() navbitmask = np.repeat(navbitmask, 2) navbit[navbitmask == 0] = 0 idx0 = int((nbit_sec - h1 * 3600) * rate) nbit_check = np.full((20, 2), np.nan) for ii in range(20): if not np.all(navbitmask[idx0 - smp*ii:idx0 - smp*ii + smp*match_len:smp]): continue idxn = np.arange(idx0 - smp * ii, idx0 - smp * ii + smp * match_len, smp) % 200 >= vld nbit = navbit[idx0 - smp*ii:idx0 - smp*ii + smp*match_len:smp].copy() for jj in [0, 1]: if jj == 1: nbit[idxn] = 1 - nbit[idxn] temp = (nbit + nbit_ref) % 2 nbit_check[ii, jj] = np.min([np.sum(temp[~bad] == 0), np.sum(temp[~bad] == 1)]) / glen ii = np.nanargmin(np.min(nbit_check, axis=1)) if min(nbit_check[ii]) >= 0.01: return None idx0 = int((nbit_sec - h1 * 3600) * rate) - m_s idx = np.arange(idx0 - smp * ii, idx0 - smp * ii + smp * ph_new.shape[0], smp) if not np.all(navbitmask[idx]): return None if nbit_check[ii][0] > nbit_check[ii][1]: navbit[idx % 200 > vld] = 1 - navbit[idx % 200 > vld] return navbit[idx].astype(float) except Exception: return None def starc_doppler_model(pos_g, vel_g, pos_l, vel_l, sec, mon, fr1, ocd, benmod): """CDAAC newroam L1 excess-Doppler model dfrm (Hz) from the CIRA+Q bending-angle climatology (Python port of ray1.f/dop1.f/benmod.f, planetiq/starc_dopmodel.py). Positions/velocities in ECI (m or km) at the RO samples; returns dfrm in the original time order and the model excess phase (m), integrated from the top of the occultation as in rawphs.f.""" from planetiq.starc_dopmodel import dopmodel sc = 1e-3 if np.nanmax(np.abs(pos_l)) > 1e5 else 1.0 # -> km, km/s o = slice(None) if ocd else slice(None, None, -1) # top-down order sg = 1.0 if ocd else -1.0 # rising: velocities negated (newroam) dfrm_td, _, ierr = dopmodel(pos_g[o] * sc, pos_l[o] * sc, sg * vel_g[o] * sc, sg * vel_l[o] * sc, mon, fr1, benmod) if ierr: raise ValueError(f'STARC Doppler model (CDAAC ray1) failed at sample {ierr}') t_td = sec[o] dt = np.abs(np.diff(t_td, prepend=t_td[0])) phmod_td = -np.cumsum(dfrm_td * dt) * (299792458.0 / fr1) # rawphs.f: wrk0 = -sum(dfrm*dt) cycles return dfrm_td[o], phmod_td[o] # RINEX-3 signal attributes of dataless pilot channels (no nav-data bits): the phase needs no bit removal, # only full-cycle connection. GPS L2C(L)/L5Q/L1C-P, Galileo E1-C/E5a-Q/E5b-Q/E6-C, BeiDou B1C-P/B2a-P/B2b-P/B3-Q. PILOT_ATTR = {'G': 'LQX', 'R': '', 'E': 'CQX', 'C': 'PQZX'} def is_pilot(gnss, obs_code): """True for a dataless pilot signal, from the RINEX observation code (e.g. 'L1C', 'L2L', 'L5P').""" return len(obs_code) >= 3 and obs_code[2] in PILOT_ATTR.get(gnss, '') def star_nobit_fallback(ph_new, gnss_ph, ocd): """STAR no-bit-file fallback of remove_cycleslip_with_navbit (its except branch), on a residual ph_new (n,2: L1, L2) in half cycles; returns the connected residual. GPS: the 20 ms bit-boundary parity is taken from the top 1 s, and half cycles are removed only at those boundaries (if the parity is ambiguous, any half-cycle jump is removed); L2 (pilot) full cycles. GLONASS: joint L1/L2 half-cycle test.""" ph_new = ph_new.copy() if gnss_ph[0] == 'G': match_len = 100 m_s = 0 if ocd else ph_new.shape[0] - match_len m_e = m_s + match_len idx = np.unique(np.where(np.round(np.diff(ph_new[m_s:m_e, 0])) % 2 == 1)[0] % 2) if idx.size != 1: ph_new[:, 0] -= np.cumsum(np.round(np.diff(ph_new[:, 0], axis=0, prepend=ph_new[:1, 0])), axis=0) else: for ii in range(idx[0] + 1, len(ph_new), 2): ph_new[ii:ii + 2, 0] -= np.round(ph_new[ii, 0] - ph_new[ii - 1, 0]) ph_new[:, 0] -= np.cumsum(np.round(np.diff(ph_new[:, 0], prepend=ph_new[:1, 0]) / 2) * 2) ph_new[:, 1] -= np.cumsum(np.round(np.diff(ph_new[:, 1], axis=0, prepend=ph_new[:1, 1]) / 2) * 2, axis=0) elif gnss_ph[0] == 'R': for ii in range(1, ph_new.shape[0]): ph_diff = ph_new[ii] - ph_new[ii - 1] if np.all(np.round(ph_diff) % 2 == 1) or np.all( (*(np.round(ph_diff) % 2 == [1, 0]), *(np.diff(np.abs(ph_diff - np.round(ph_diff))) > 0))) or np.all( (*(np.round(ph_diff) % 2 == [0, 1]), *(np.diff(np.abs(ph_diff - np.round(ph_diff))) < 0))): ph_new[ii] += 1 ph_new -= np.cumsum(np.round(np.diff(ph_new, axis=0, prepend=ph_new[:1]) / 2) * 2, axis=0) else: ph_new -= np.cumsum(np.round(np.diff(ph_new, axis=0, prepend=ph_new[:1]) / 2) * 2, axis=0) return ph_new def connect_phase_starc(ph, ph_mod, freq, gnss_ph, ocd, bits, codes=('', ''), nobit='star'): """STARC phase connection: STAR implementation of the CDAAC newroam rawphs.f algorithm, in STAR units. ph: (n,2) excess phase L1, L2 (m) before cycle-slip handling; ph_mod: L1 model excess phase (m) from starc_doppler_model; bits: aligned nav bits (half cycles) or None; codes: RINEX observation codes of the two signals (e.g. 'L1C', 'L2L'), used to tell data from pilot signals. For each signal, residual to the model (half cycles), then: - data signal with bits (GPS L1 C/A; GLONASS L1, L2): bits removed, connection to +-pi (pha2dop); - data signal without bits: nobit = 'star' (default): STAR fallback (star_nobit_fallback) on the residual to the model; nobit = 'flip': CDAAC phasor-flip rule (|dphi| > pi/2 flipped), i.e. connection to +-pi/2 (rawphs.f); - pilot signal (no data bits; GPS L2C-L, Galileo E1-C/E5b-Q, BeiDou B1C-P/B2a-P): connection to +-pi. The L2 model is the L1 model in metres (rawphs.f: L1 Doppler model scaled by f2/f1). rawphs.f connects L2 with pha2dop1 (removes only multiple-cycle jumps), which assumes an already connected L2; PlanetiQ open-loop L2 still carries single-cycle wraps, so the rules above are used for L2 as well. Processing runs top-down (rawphs.f reverses rising occultations).""" c0 = 299792458.0 if bits is None and nobit == 'star' and gnss_ph[0] in 'GR': # no bit file for a data signal: STAR's fallback on the residual to the CIRA+Q model # (better than the CDAAC phasor-flip rule for PlanetiQ: GPS bit-boundary parity, GLONASS L1/L2 test) r = (ph - ph_mod[:, None]) * (2 * freq / c0) r = star_nobit_fallback(r, gnss_ph, ocd) return r * (c0 / (2 * freq)) + ph_mod[:, None] o = slice(None) if ocd else slice(None, None, -1) out = np.empty_like(ph) for k in (0, 1): m = ph_mod r = ((ph[:, k] - m) * (2 * freq[0, k] / c0))[o] # half cycles, top-down d = np.diff(r, prepend=r[:1]) pilot = is_pilot(gnss_ph[0], codes[k]) bits_k = bits is not None and not pilot and (k == 0 or gnss_ph[0] == 'R') if bits_k: r = r + bits[o] d = np.diff(r, prepend=r[:1]) d = d - 2 * np.round(d / 2) # pha2dop (+-pi) elif pilot: d = d - 2 * np.round(d / 2) # pha2dop (+-pi) else: d = d - np.round(d) # phasor-flip rule (+-pi/2) out[:, k] = (np.cumsum(d)[o]) * (c0 / (2 * freq[0, k])) + m return out def calculate_exph(rofile, reffiles, tarsource='', out_date=None): from re import sub from scipy.interpolate import CubicSpline from planetiq.frame_conversions import ro_ecef2eci, ro_eci2ecef #Jun modified from subprocess import run, TimeoutExpired, CalledProcessError import subprocess # TC: updated on 2026/04/08: add zoneinfo module for output UTC time from zoneinfo import ZoneInfo from configuration import dev_dir #Jun modified exopt = exph_options() # STAR 2026-10-30: clock / phase-connection options ## Xinjia modified, for RFSI/ROPP import os from configuration import version, rfsi_dir, ropp_path, suffix ucar2ropp = f"{ropp_path}/ropp/ropp-10.0_{version}/ropp_io/tools/ucar2ropp" # STAR 2026-11-04: exph_ropp_pp = 'ropp_pp_y' runs Y_PM_GO (Y. Chen's ropp_pp build) instead of the STAR/delivered ropp_pp RFSI_occ_tool = f"{ropp_path}/ropp/ropp-10.0_{version}/{exopt['ropp_pp']}/tools/ropp_pp_occ_tool" ropp2bufr = f"{ropp_path}/ropp/ropp-10.0_{version}/ropp_io/tools/ropp2bufr_eccodes" gtsheaders_bufr = f"{ropp_path}/ropp/ropp-10.0_{version}/ropp_io/tools/gtsheaders_bufr.pl" bufr_table = f"{ropp_path}/bufr_decoder/tables" env = os.environ.copy() env["BUFR_TABLES"] = f"{ropp_path}/ropp/ropp-10.0_{version}/data/bufr/" env["ROPP_THIN"] = f"{ropp_path}/ropp/ropp-10.0_{version}/ropp_io/data/ropp_thin_asglog-247.dat" ## Xinjia modified, for RFSI/ROPP and BUFR # Constants version = '1.0' dt0_gps = datetime(1980,1,6) antennas = ['POD0', 'POD1', 'OCC2', 'OCC3', 'OCC4', 'OCC5', 'OCC6', 'OCC7'] # Choose reference file dt_ro = [datetime.strptime(x,'%Y%m%d%H%M%S') for x in sub(r'.*(\d{14}\.\d{14}).*',r'\1',rofile).split('.')] dt_rr = datetime.strptime(sub(r'.*\.(\d{4}\.\d{3})\..*',r'\1',rofile),'%Y.%j') sec_rr = (dt_rr - dt0_gps) / timedelta(seconds=1) good = False # STAR 2026-10-31: zero differencing does not use the reference link, so no reference profile is searched for: # one pass with reffile = None, and the orbit windows are taken from the occultation itself, widened by # exph_zd_win (default 300 s) on each side. The LEO orbit is decimated to >= 20 nodes over the window and # interpolated with an 11-point Lagrange; a window only 60 s wider than the occultation leaves asymmetric # stencils near its ends (phase drift vs CDAAC worse than with the reference window); 300 s is better than both. zero_ref = exopt['clock'] == 'zero_diff' and exopt['zd_noref'] for reffile in ([None] if zero_ref else reffiles): if reffile is not None: with np.load(reffile, allow_pickle=True) as fid: data_ref = fid['data'] idx = np.logical_and(data_ref[0, 0] + data_ref[:, 1] > (dt_ro[0] - dt0_gps) / timedelta(seconds=1) - 30, data_ref[0, 0] + data_ref[:, 1] < (dt_ro[1] - dt0_gps) / timedelta(seconds=1) + 30) yy = data_ref[idx, 3] - data_ref[idx, 2] if max(yy)-min(yy)>1000: continue # Extract info from file names leoName = sub(r'.*occTab_(....)\..*',r'\1',rofile) ant_leo_ro = int(sub(r'.*\.A(\d\d)\..*',r'\1',rofile)) gnss_ro = sub(r'.*\.(...).npz',r'\1',rofile) if reffile is None: # zero differencing: windows from the occultation, the RO link stands in for the reference ant_leo_ref = ant_leo_ro dt_ref = [dt_ro[0] - timedelta(seconds=exopt['zd_win']), dt_ro[1] + timedelta(seconds=exopt['zd_win'])] gnss_ref = gnss_ro else: ant_leo_ref = int(sub(r'.*\.A(\d\d)\..*',r'\1',reffile)) dt_ref = [datetime.strptime(x,'%Y%m%d%H%M%S') for x in sub(r'.*(\d{14}\.\d{14}).*',r'\1',reffile).split('.')] gnss_ref = sub(r'.*\.(...).npz', r'\1', reffile) f1_ro, f2_ro, gprn_ro = get_sat_freq(gnss_ro) c2_ro = 1 / ((f1_ro / f2_ro) ** 2 - 1) c1_ro = c2_ro + 1 f1_ref, f2_ref, gprn_ref = get_sat_freq(gnss_ref) c2_ref = 1 / ((f1_ref / f2_ref) ** 2 - 1) c1_ref = c2_ref + 1 # ant_gps = np.loadtxt("GPS_Antenna_Infos.Final") # ant_gps_start_time = [datetime(*map(int, x[1:7])) for x in ant_gps] # ant_gps_end_time = [datetime(*map(int, x[7:13])) for x in ant_gps] with np.load('../gnss/GPS_Antenna_Infos.npz', allow_pickle=True) as fid: antoffset_gnss_ro = fid[gnss_ro][()] antoffset_gnss_ref = fid[gnss_ref][()] antoffset_gnss_ro = antoffset_gnss_ro['antenna'][np.where(np.logical_and( antoffset_gnss_ro['dt_rng'][:, 0] < dt_rr, antoffset_gnss_ro['dt_rng'][:, 1] > dt_rr))[0][0]][:3] / 1E3 antoffset_gnss_ref = antoffset_gnss_ref['antenna'][np.where(np.logical_and( antoffset_gnss_ref['dt_rng'][:, 0] < dt_rr, antoffset_gnss_ref['dt_rng'][:, 1] > dt_rr))[0][0]][:3] / 1E3 # Read in gnss position/velocity file accordingly # Gdata, Gdt, Gsats = read_gnss(dt_ref[1] + timedelta(minutes=10), # (dt_ref[0] - timedelta(minutes=10)).day != (dt_ref[1] + timedelta(minutes=10)).day) # if gnss_ro not in Gsats: # raise LookupError(f'RO satellite {gnss_ro} not in Gsats') # if gnss_ref not in Gsats: # raise LookupError(f'REF satellite {gnss_ref} not in Gsats') # gidx = np.hstack([np.where(Gsats == gnss_ro), np.where(Gsats == gnss_ref)])[0] # Gdata = Gdata[:,gidx,:] Gdata, Gdt = read_gnss_v2(dt_ref, [gnss_ro, gnss_ref]) tidx = np.logical_and(Gdt > dt_ref[0] - timedelta(minutes=10), Gdt < dt_ref[1] + timedelta(minutes=10)) dts_gnss = Gdt[tidx] sec_gnss = dt2float(dts_gnss, dt_rr) pos_gnss = Gdata[tidx,:,:4] vel_gnss = Gdata[tidx,:,4:] del Gdata, Gdt # Read in leo position/velocity/attitude files accordingly # TC: updated 2026-10-27 (added 1 line, new attribute) pod_method = 'PlanetiQ' with np.load(f"s01_leoOrb/leoOrb_{leoName}.{dt_rr:%Y.%j}.npz", allow_pickle=True) as fid: Ldata = fid['data'] Ldts = fid['dts'] # STAR 2026-11-01: POD quality around the occultation (gap and position/velocity consistency) podq = pod_quality(dt2float(Ldts, dt_rr), Ldata[:, :3], Ldata[:, 4:7], dt2float(np.array([dt_ro[0]]), dt_rr)[0], dt2float(np.array([dt_ro[1]]), dt_rr)[0], pv_win=30.0) podq['qc'] = int(podq['gap_after_s'] < exopt['pod_gap_win']) + 2 * int(podq['pv_max'] > exopt['pod_pv_max']) # STAR 2026-11-02: replace the converging POD orbit after a gap by a dynamic fit (exph_pod_fit s; bit 3 = replaced) podfit = dict(status='off', rms=-1.0, maxdiff=-1.0) if exopt['pod_fit'] > 0: from planetiq.pod_fit import replace_after_gap Ldata, podfit = replace_after_gap(Ldts, Ldata, dt_ro[0], dt_ro[1], exopt['pod_fit'], exopt['pod_fit_len'], exopt['pod_fit_min_after'], leo=leoName) if podfit['status'] == 'fit': podq['qc'] |= 4 pod_method = 'DMF' # STAR 2026-11-03: orbit after the gap from Dynamic Model Fitting (else 'PlanetiQ') elif podfit['status'] in ('nosatkit', 'failed'): log_message(f"-- POD fit {podfit['status']} for {os.path.basename(rofile)}: POD orbit kept") # STAR 2026-11-04: DMF for a flag-2 stretch (POD velocity inconsistent with the positions, away from a gap): # replace the stretch by a dynamic fit bridging the POD before and after it in the current tar, blended into the # POD over exph_pod_fit_pv_taper s at both edges (pod_qc bit value 8). No POD after the stretch -> flag 2 kept. pvfit = dict(status='off', rms=-1.0, maxdiff=-1.0, span_s=0.0) if exopt['pod_fit_pv']: from planetiq.pod_fit import replace_pv_segment Ldata, pvfit = replace_pv_segment(Ldts, Ldata, dt_ro[0], dt_ro[1], pv_max=exopt['pod_pv_max'], margin=exopt['pod_fit_pv_margin'], fit_len=exopt['pod_fit_len'], mode=exopt['pod_fit_pv_mode'], max_side_span=exopt['pod_fit_pv_side_max'], taper=exopt['pod_fit_pv_taper'], leo=leoName) if pvfit['status'] in ('bridge', 'before', 'after'): podq['qc'] |= 8 pod_method = 'DMF' podq['pv_after'] = pod_quality(dt2float(Ldts, dt_rr), Ldata[:, :3], Ldata[:, 4:7], dt2float(np.array([dt_ro[0]]), dt_rr)[0], dt2float(np.array([dt_ro[1]]), dt_rr)[0], pv_win=30.0)['pv_max'] # STAR 2026-11-04: L1A quality flag bad = "1" (existing attribute, errstr = reason) when the POD is still bad # after the fixes (pod_qc records what was found and done): # gap_unfitted : orbit resumed after a gap < exph_pod_gap_win s before the occultation (bit 1) and not replaced # pv_after : 30 s mean |dP/dt - V| around the occultation still > exph_pod_pv_max on the orbit used # fit_rms : a fit was applied but its residual exceeds exph_pod_fit_rms_max (m) reasons = [] if podq['qc'] & 1 and not podq['qc'] & 4: reasons.append('gap_unfitted') if podq['pv_after'] > exopt['pod_pv_max']: reasons.append('pv_after') if (podfit['status'] == 'fit' and podfit['rms'] > exopt['pod_fit_rms_max']) or \ (pvfit['status'] in ('bridge', 'before', 'after') and pvfit['rms'] > exopt['pod_fit_rms_max']): reasons.append('fit_rms') podq['bad'] = int(bool(reasons)); podq['bad_reason'] = ','.join(reasons) if reasons else 'none' # STAR 2026-11-05: exph_skip_dup = 'defer' -- the POD needed to fix the orbit is not there yet (gap without 900 s of # POD after it, or a flag-2 stretch reaching the end of the POD): do not write; step4 retries with the next tar if exopt['defer_now'] and podq['bad'] and (('gap_unfitted' in reasons and podfit['status'] == 'nodata') or ('pv_after' in reasons and pvfit['status'] == 'nodata')): raise DeferOccultation(podq['bad_reason']) if exopt['pod_qc'] == 'reject' and podq['bad']: # STAR 2026-11-04: reject when still bad after the fixes raise ValueError(f"POD QC: {podq['bad_reason']} (|dP/dt - V| {podq['pv_after']:.3f} m/s after the fixes, " f"orbit resumed {podq['gap_after_s']:.0f} s before the occultation)") tidx = np.logical_and(Ldts > dt_ref[0] - timedelta(minutes=1), Ldts < dt_ref[1] + timedelta(minutes=1)) dts_leoorb = Ldts[tidx] pos_leo = Ldata[tidx,:4] vel_leo = Ldata[tidx,4:] del Ldata, Ldts with np.load(f"s01_leoAtt/leoAtt_{leoName}.{dt_rr:%Y.%j}.npz", allow_pickle=True) as fid: Ldata = fid['data'] Ldts = fid['dts'] Lvar = fid['vars'] tidx = np.logical_and(Ldts > dt_ref[0] - timedelta(minutes=1), Ldts < dt_ref[1] + timedelta(minutes=1)) dts_leoatt = Ldts[tidx] att_leo = Ldata[tidx,:][:,['att' in x for x in Lvar]] aridx = np.cumsum(np.hstack(([False], np.any(np.abs(np.diff(att_leo, axis=0)) > 0.1, axis=1)))) % 2 == 1 att_leo[aridx] *= -1 del Ldata, Ldts, Lvar if dts_leoatt.size < 2 or dts_leoorb.size < 2: raise ValueError(f'No LEO orbit/attitude data between {dt_ref[0]} and {dt_ref[1]}') tidx = np.logical_and(dts_leoorb >= dts_leoatt[0], dts_leoorb <= dts_leoatt[-1]) dts_leo = dts_leoorb[tidx] sec_leo = dt2float(dts_leo, dt_rr) pos_leo = pos_leo[tidx] vel_leo = vel_leo[tidx] att_leo = CubicSpline(dt2float(dts_leoatt,dts_leo[0]), att_leo)(dt2float(dts_leo,dts_leo[0])) att_leo /= np.sum(att_leo ** 2, axis=1)[:,None] ** 0.5 del dts_leoorb, dts_leoatt sec_leo_full, clk_leo_full = sec_leo.copy(), pos_leo[:, 3].copy() # STAR 2026-10-30: 1 s POD clock (us) tsmp = 30 tspc = np.round(np.linspace(0,len(sec_leo)-1,min(max(-((len(sec_leo)-1)//-tsmp)+1,20),len(sec_leo)))).astype(int) leo_order = 11 if exopt['leo_interp'] != 'decimated': # STAR 2026-10-31: keep the 1 s POD orbit/attitude as nodes if exopt['leo_interp'] == 'smooth': # local least-squares cubic over +-15 s (Savitzky-Golay), per continuous segment pos_leo, vel_leo = smooth_1s(sec_leo, pos_leo, vel_leo, exopt['leo_smooth_hw']) tspc = np.arange(len(sec_leo)) leo_order = 4 # 8-node local Lagrange dts_leo = dts_leo[tspc] sec_leo = sec_leo[tspc] pos_leo = pos_leo[tspc] vel_leo = vel_leo[tspc] att_leo = att_leo[tspc] # Load ro/ref files with np.load(rofile, allow_pickle=True) as fid: data_ro = fid['data'] header_ro = fid['header'] if reffile is not None: with np.load(reffile, allow_pickle=True) as fid: data_ref = fid['data'] header_ref = fid['header'] # data_ro[:, 2] = data_ro[:, 2] + data_ro[:, 3] # data_ro[:, 4] = data_ro[:, 4] + data_ro[:, 5] # TC: updated 2026-10-27 (added 12 lines, trim boundary nan data and interp inbetween) for ii in range(data_ro.shape[0]): if np.all(np.isfinite(data_ro[ii])): break data_ro = data_ro[ii:] for ii in range(data_ro.shape[0]-1,-1,-1): if np.all(np.isfinite(data_ro[ii])): break data_ro = data_ro[:ii+1] for ii in range(2, data_ro.shape[1]): idx = np.isfinite(data_ro[:, ii]) if not np.all(idx): data_ro[~idx, ii] = np.interp(np.where(~idx)[0],np.where(idx)[0],data_ro[idx, ii]) l_ro = np.where([x[0]=='L' and '(M)' not in x for x in header_ro])[0] if reffile is not None: l_ref = np.where([x[0]=='L' and '(M)' not in x for x in header_ref])[0] m_ro = np.where(['(M)' in x for x in header_ro])[0] s_ro = np.where([x[0]=='S' for x in header_ro])[0] sec_ro = data_ro[:, 0] - sec_rr + data_ro[:, 1] if reffile is not None: s_ref = np.where([x[0]=='S' for x in header_ref])[0] sec_ref = data_ref[:, 0] - sec_rr + data_ref[:, 1] tidx = np.logical_and(sec_ref > sec_ro[0] - 30, sec_ref < sec_ro[-1] + 30) data_ref = data_ref[tidx] sec_ref = sec_ref[tidx] if sec_leo[0]>sec_ro[0] or sec_leo[-1]> {sec_leo[-1]:.1f}) could not cover RO event ({sec_ro[0]:.1f} >> {sec_ro[-1]:.1f})!') sec_ro -= CubicSpline(sec_leo, pos_leo[:,-1])(sec_ro)/1E6 if reffile is not None: if sec_leo[0]>sec_ref[0] or sec_leo[-1]> {sec_leo[-1]:.1f}) could not cover REF event ({sec_ref[0]:.1f} >> {sec_ref[-1]:.1f})!') sec_ref -= CubicSpline(sec_leo, pos_leo[:,-1])(sec_ref)/1E6 tidx = np.logical_and(sec_gnss > sec_ro[0] - 905, sec_gnss < sec_ro[-1] + 905) dts_gnss = dts_gnss[tidx] sec_gnss = sec_gnss[tidx] pos_gnss = pos_gnss[tidx] vel_gnss = vel_gnss[tidx] if ~np.all(np.isfinite(pos_gnss[:, 0, 3])): raise ValueError(f'{np.sum(~np.isfinite(pos_gnss[:, 0, 3]))}/{sec_gnss.size} clock bias value for {gnss_ro} is not finite.') tidx = np.logical_and(sec_leo > sec_ro[0] - 905, sec_leo < sec_ro[-1] + 905) dts_leo = dts_leo[tidx] sec_leo = sec_leo[tidx] pos_leo = pos_leo[tidx] vel_leo = vel_leo[tidx] att_leo = att_leo[tidx] ant_ro_ecef = conv_gnss_ant(pos_gnss[:, 0, :3], antoffset_gnss_ro, dts_gnss) ant_ref_ecef = conv_gnss_ant(pos_gnss[:, 1, :3], antoffset_gnss_ref, dts_gnss) # TC: updated 2026-11-09 (1 line, pass argument without clock) pos_leo_eci, vel_leo_eci = ro_ecef2eci(pos_leo[:, :3], vel_leo, dts_leo) #pos_leo_eci, vel_leo_eci = ro_ecef2eci(pos_leo, vel_leo, dts_leo) pos_ro_eci, vel_ro_eci = ro_ecef2eci(pos_gnss[:, 0, :3] + ant_ro_ecef, vel_gnss[:, 0, :], dts_gnss) pos_ref_eci, vel_ref_eci = ro_ecef2eci(pos_gnss[:, 1, :3] + ant_ref_ecef, vel_gnss[:, 1, :], dts_gnss) pos_leo_eci = np.hstack([pos_leo_eci, pos_leo[:,[-1]]]) vel_leo_eci = np.hstack([vel_leo_eci, pos_leo[:, [-1]]]) pos_ro_eci = np.hstack([pos_ro_eci, pos_gnss[:, 0, [3]]]) vel_ro_eci = np.hstack([vel_ro_eci, vel_gnss[:, 0, [3]]]) pos_ref_eci = np.hstack([pos_ref_eci, pos_gnss[:, 1, [3]]]) vel_ref_eci = np.hstack([vel_ref_eci, vel_gnss[:, 1, [3]]]) with np.load('antenna_offset.npz') as fid: for a_ in (ant_leo_ro, ant_leo_ref): # STAR 2026-11-01: clear message for a satellite without offsets if f'ANT_{leoName}_{antennas[a_]}' not in fid.files: raise ValueError(f"no antenna offset ANT_{leoName}_{antennas[a_]} in antenna_offset.npz") antoffset_leo_ro = fid[f'ANT_{leoName}_{antennas[ant_leo_ro]}'][12:15] antoffset_leo_ref = fid[f'ANT_{leoName}_{antennas[ant_leo_ref]}'][12:15] obj_gnss_ro, obj_leo_ro = phase_correction(data_ro[:, l_ro[0]], pos_ro_eci, vel_ro_eci, sec_gnss, pos_leo_eci, vel_leo_eci, att_leo, sec_leo, sec_ro, antoffset_leo_ro, True, leo_order) phL1_ro = data_ro[:, l_ro[0]] - obj_gnss_ro['range_diff_T'] + obj_gnss_ro['GL_term'] phL2_ro = data_ro[:, l_ro[1]] - obj_gnss_ro['range_diff_T'] + obj_gnss_ro['GL_term'] if reffile is not None: obj_gnss_ref, obj_leo_ref = phase_correction(data_ref[:, l_ref[0]], pos_ref_eci, vel_ref_eci, sec_gnss, pos_leo_eci, vel_leo_eci, att_leo, sec_leo, sec_ref, antoffset_leo_ref, False, leo_order) phL1_ref = data_ref[:, l_ref[0]] - obj_gnss_ref['range_diff_T'] + obj_gnss_ref['GL_term'] coef = np.polyfit(sec_ref, phL1_ref, 1) if leoName.startswith('YM') and (np.abs(coef[0]) > 1 or np.max(np.abs(phL1_ref - np.polyval(coef, sec_ref))) > 10): continue phL2_ref = data_ref[:, l_ref[1]] - obj_gnss_ref['range_diff_T'] + obj_gnss_ref['GL_term'] dph_ref = phL1_ref - phL2_ref dph_ref_smth = np.hstack( [np.cumsum(dph_ref[:49])[::2] / np.arange(1, 51, 2), np.convolve(dph_ref, np.ones((51,)), 'valid') / 51, (np.cumsum(dph_ref[:-50:-1])[::2] / np.arange(1, 51, 2))[::-1]]) termL1_ref = phL1_ref + c2_ref * dph_ref_smth termL2_ref = phL2_ref + c1_ref * dph_ref_smth if exopt['clock'] == 'single_diff': # legacy: LEO clock from the POD reference link termL1_ro = CubicSpline(sec_ref, termL1_ref)(sec_ro) termL2_ro = CubicSpline(sec_ref, termL2_ref)(sec_ro) else: # proposed: zero differencing with the POD LEO clock if exopt['clk_interp'] == 'linear': clk_m_ro = np.interp(sec_ro, sec_leo_full, clk_leo_full) / 1E6 * 299792458.0 elif exopt['clk_interp'] == 'poly': # CDAAC genXphase52 GTLCL0 tmid = 0.5 * (sec_ro[0] + sec_ro[-1]) kk = np.abs(sec_leo_full - tmid) <= exopt['clk_poly_dur'] / 2 if kk.sum() <= exopt['clk_poly_deg'] + 1: raise ValueError('Too few LEO clock values for the polynomial clock fit') cpf = np.polyfit(sec_leo_full[kk] - tmid, clk_leo_full[kk], exopt['clk_poly_deg']) clk_m_ro = np.polyval(cpf, sec_ro - tmid) / 1E6 * 299792458.0 else: clk_m_ro = obj_leo_ro['clk_m'] termL1_ro = clk_m_ro + obj_leo_ro['delt_t_rel_r'] * 299792458.0 termL2_ro = termL1_ro.copy() # termL1_ro = obj_leo_ro['clk_m'] + obj_leo_ro['delt_t_rel_r'] * c0 # Alternative method # termL2_ro = obj_leo_ro['clk_m'] + obj_leo_ro['delt_t_rel_r'] * c0 exL1 = phL1_ro - termL1_ro exL2 = phL2_ro - termL2_ro exL1 -= exL1[0 if obj_leo_ro['ocd'] else -1] exL2 -= exL2[0 if obj_leo_ro['ocd'] else -1] dt_ros = [dt_rr + timedelta(seconds=x) for x in sec_ro] #pos_leo_ecef, vel_leo_ecef = ro_eci2ecef(obj_leo_ro['pos'], obj_leo_ro['vel'], dt_ros) #pos_gnss_ecef, vel_gnss_ecef = ro_eci2ecef(obj_gnss_ro['pos'], obj_gnss_ro['vel'], dt_ros) ## Jan 8, 2026 at 10:51 AM pos_leo_ecef, vel_leo_ecef, pos_leo_tod, vel_leo_tod = ro_eci2ecef(obj_leo_ro['pos'], obj_leo_ro['vel'], dt_ros, tod=True) pos_gnss_ecef, vel_gnss_ecef, pos_gnss_tod, vel_gnss_tod = ro_eci2ecef(obj_gnss_ro['pos'], obj_gnss_ro['vel'], dt_ros, tod=True) exL1_LM = data_ro[:,m_ro[0]] - (data_ro[:,l_ro[0]] - exL1) # TC: updated 2026-10-27 (3 lines, adapt new remove_cycleslip method) connect_ro = exopt['connect'] if gnss_ro[0] in exopt['connect_gnss'] else 'star' if connect_ro == 'star': # legacy with_navbit, (exL1, exL2) = remove_cycleslip_with_navbit(np.column_stack([exL1, exL2]), exL1_LM[:, None], np.array([[f1_ro, f2_ro]]), gnss_ro, obj_leo_ro['ocd'], dt_rr, sec_ro) else: # STARC (CDAAC newroam rawphs.f algorithm) bits = navbit_aligned(np.column_stack([exL1, exL2]), exL1_LM[:, None], np.array([[f1_ro, f2_ro]]), gnss_ro, obj_leo_ro['ocd'], dt_rr, sec_ro) with_navbit = bits is not None _, phmod_ro = starc_doppler_model(obj_gnss_ro['pos'], obj_gnss_ro['vel'], obj_leo_ro['pos'], obj_leo_ro['vel'], sec_ro, dt_rr.month, f1_ro, obj_leo_ro['ocd'], exopt['benmod']) exL1, exL2 = connect_phase_starc(np.column_stack([exL1, exL2]), phmod_ro, np.array([[f1_ro, f2_ro]]), gnss_ro, obj_leo_ro['ocd'], bits, (str(header_ro[l_ro[0]]), str(header_ro[l_ro[1]])), exopt['starc_nobit']).T # exL1 = remove_cycleslip(exL1, exL1_LM, f1_ro, gnss_ro[0] + '1', obj_leo_ro['ocd']) # exL2 = remove_cycleslip(exL2, exL1_LM, f2_ro, gnss_ro[0] + '2', obj_leo_ro['ocd']) rofile_tkn = sub(r'.*_(.*).npz',r'\1',rofile) #fortran_input = f'/dev/shm/jzhou128/in_{rofile_tkn}.txt' #fortran_output = f'/dev/shm/jzhou128/out_{rofile_tkn}.txt' fortran_input = f'{dev_dir}/in_{rofile_tkn}.txt' #Jun modified fortran_output = f'{dev_dir}/out_{rofile_tkn}.txt' #Jun modified from os import remove, makedirs, path if connect_ro == 'star' and exopt['cycleslip'] == 'python': # legacy cycle-slip step, Python port (planetiq/ropp_openloop.py) from planetiq.ropp_openloop import cyclesliperemoval as _csr try: f_l1, f_l2, f_lm = _csr(sec_ro, exL1.copy(), exL2.copy(), pos_leo_ecef, pos_gnss_ecef, 1 if obj_leo_ro['ocd'] else -1, float(int(f1_ro)), float(int(f2_ro)), dt_rr.month, exopt['msis_file'], pilot_l1=exopt['cycleslip_pilot'] and is_pilot(gnss_ro[0], str(header_ro[l_ro[0]]))) if np.max(np.abs(f_l1 - exL2)) < np.max(np.abs(exL1 - exL2)): exL1 = f_l1 exL1_LM = f_lm except Exception as e: # the Fortran program failing is not fatal either log_message(f'cycle-slip step (Python) skipped for {rofile}: {type(e).__name__} {e}') elif connect_ro == 'star': # legacy Fortran cycle-slip step (not part of the CDAAC algorithm) with open(fortran_input, 'w') as fid: fid.write(f'{len(sec_ro)} {1 if obj_leo_ro["ocd"] else -1} {int(f1_ro)} {int(f2_ro)}\n') for n in range(len(sec_ro)): fid.write(f"{sec_ro[n]:15.8f}{exL1[n]:20.8f}{exL2[n]:20.8f}{data_ro[n,m_ro[0]]:20.8f}" f"{pos_leo_ecef[n,0]:15.5f}{pos_leo_ecef[n,1]:15.5f}{pos_leo_ecef[n,2]:15.5f}" f"{pos_gnss_ecef[n,0]:15.5f}{pos_gnss_ecef[n,1]:15.5f}{pos_gnss_ecef[n,2]:15.5f}3{dt_rr.month:2d}\n") run(['PhaseModel/cyclesliperemoval', '-i', fortran_input, '-o', fortran_output], stdout=-3, stderr=-3, timeout=5) # TC: updated on 2026/04/07 (7 -> 8 lines) remove(fortran_input) if path.exists(fortran_output): fortran_data = np.loadtxt(fortran_output) remove(fortran_output) if np.max(np.abs(fortran_data[:, 0] - exL2)) < np.max(np.abs(exL1 - exL2)): exL1 = fortran_data[:, 0] exL1_LM = fortran_data[:, 2] idx = np.where(np.isfinite(exL1))[0] exL1 -= exL1[idx[0 if obj_leo_ro['ocd'] else -1]] exL2 -= exL2[idx[0 if obj_leo_ro['ocd'] else -1]] exLC = exL1 + c2_ro * (exL1 - exL2) vec_lg = pos_gnss_ecef - pos_leo_ecef proj_lg = pos_leo_ecef - vec_lg * (np.sum(pos_leo_ecef * vec_lg, axis=1) / np.sum(vec_lg * vec_lg, axis=1))[:, None] from astropy.coordinates import EarthLocation alt_km = EarthLocation.from_geocentric(*proj_lg.T, unit='m').geodetic.height.value / 1E3 if np.logical_or(np.all(alt_km < 90), np.all(alt_km > 0)): raise ValueError(f'Non-RO direct height range: {np.min(alt_km):.3f} >> {np.max(alt_km):.3f} (km)') # TC: updated 2026/07/14 (remove 69 lines, where I kept as comments below) # xx = sec_ro[::1 if obj_leo_ro['ocd'] else -1] # yy = exLC[::1 if obj_leo_ro['ocd'] else -1] # idx = np.isfinite(yy) # xx = xx[idx] # yy = yy[idx] # len_old = np.where(np.abs(yy)>0.1)[0][0] # coef = np.polyfit(xx[:len_old], yy[:len_old], 1) # len_new = np.where(np.abs(np.polyval(coef, xx) - yy) > 0.1)[0][0] # while len_new > len_old: # len_old = len_new # coef = np.polyfit(xx[:len_old], yy[:len_old], 1) # len_new = np.where(np.abs(np.polyval(coef, xx) - yy) > 0.1)[0][0] # trend_LC = np.polyval(coef, sec_ro) # # idx = np.where(np.logical_and(alt_km > max(90, max(alt_km) - 20), np.isfinite(exLC)))[0] # # idx = np.where(np.logical_and(np.abs(alt_km - 100) < 10, np.isfinite(exLC)))[0] # # trend_LC = np.polyval(np.polyfit(sec_ro[idx], exLC[idx], 1), sec_ro) # exL1 -= trend_LC # exL2 -= trend_LC # exLC -= trend_LC # # TC: updated on 2026/04/07 (36 lines, followed by the commented original block) # method = 'latbindn' # if method == 'uniform': # idx = np.logical_and(alt_km > 70, alt_km < 90) # coef = np.polyfit(sec_ro[idx], exLC[idx], 1) # if abs(coef[0]) < 0.001: # exL1 += (0.002 - abs(coef[0])) * sec_ro * (1 if obj_leo_ro['ocd'] else -1) # exL2 += (0.002 - abs(coef[0])) * sec_ro * (1 if obj_leo_ro['ocd'] else -1) # exLC += (0.002 - abs(coef[0])) * sec_ro * (1 if obj_leo_ro['ocd'] else -1) # elif method == 'latbin': # idx = np.logical_and(alt_km > 80, alt_km < 90) # coef = np.polyfit(sec_ro[idx], exLC[idx], 1) # idx = np.logical_and(alt_km > 89, alt_km < 91) # ecfLC = np.mean(proj_lg[idx], axis=0) # latLC = int(min((np.arcsin(ecfLC[2] / np.sum(ecfLC ** 2) ** 0.5) / np.pi * 180 + 90) // 20, 8)) # threshLC = [0.002, 0.0015, 0.001, 0.001, 0.00025, 0.001, 0.0005, 0.001, 0.0015][latLC] # slopeLC = [0.003, 0.0025, 0.0022, 0.0018, 0.0005, 0.0001, 0.001, 0.002, 0.0025][latLC] # if abs(coef[0]) < threshLC: # exL1 += (slopeLC - abs(coef[0])) * sec_ro * (1 if obj_leo_ro['ocd'] else -1) # exL2 += (slopeLC - abs(coef[0])) * sec_ro * (1 if obj_leo_ro['ocd'] else -1) # exLC += (slopeLC - abs(coef[0])) * sec_ro * (1 if obj_leo_ro['ocd'] else -1) # elif method == 'latbindn': # idx = np.logical_and(alt_km > 80, alt_km < 90) # coef = np.polyfit(sec_ro[idx], exLC[idx], 1) # idx = np.logical_and(alt_km > 89, alt_km < 91) # ecfLC = np.mean(proj_lg[idx], axis=0) # latLC = int(min((np.arcsin(ecfLC[2] / np.sum(ecfLC ** 2) ** 0.5) / np.pi * 180 + 90) // 20, 8)) # if (EarthLocation.from_geocentric(*ecfLC,unit='m').geodetic.lon.value/360+np.mean(sec_ro[idx])/86400-0.25)%1<0.5: # Daytime # threshLC = [0.002, 0.0015, 0.0015, 0.0015, 0.0005, 0.0012, 0.0005, 0.001, 0.0015][latLC] # slopeLC = [0.003, 0.0025, 0.0022, 0.002, 0.0015, 0.0001, 0.00075, 0.0018, 0.0025][latLC] # else: # threshLC = [0.002, 0.0015, 0.001, 0.001, 0.0005, 0.001, 0.0005, 0.001, 0.0015][latLC] # slopeLC = [0.003, 0.002, 0.002, 0.0018, 0.00075, 0.0001, 0.001, 0.002, 0.0025][latLC] # if abs(coef[0]) < threshLC: # exL1 += (slopeLC - abs(coef[0])) * sec_ro * (1 if obj_leo_ro['ocd'] else -1) # exL2 += (slopeLC - abs(coef[0])) * sec_ro * (1 if obj_leo_ro['ocd'] else -1) # exLC += (slopeLC - abs(coef[0])) * sec_ro * (1 if obj_leo_ro['ocd'] else -1) # # ## updated in 01/08/2026 # # idx = np.logical_and(alt_km > 70, alt_km < 90) # # coef = np.polyfit(sec_ro[idx], exLC[idx], 1) # # if abs(coef[0]) < 0.001: # # exL1 += (0.002 - abs(coef[0])) * sec_ro * (1 if obj_leo_ro['ocd'] else -1) # # exL2 += (0.002 - abs(coef[0])) * sec_ro * (1 if obj_leo_ro['ocd'] else -1) # # exLC += (0.002 - abs(coef[0])) * sec_ro * (1 if obj_leo_ro['ocd'] else -1) # idx = np.where(np.isfinite(exL1))[0] # exL1 -= exL1[idx[0 if obj_leo_ro['ocd'] else -1]] # exL2 -= exL2[idx[0 if obj_leo_ro['ocd'] else -1]] # exLC -= exLC[idx[0 if obj_leo_ro['ocd'] else -1]] # STAR 2026-11-01: output folder named after the observation day (delivered) or, with exph_out_date = 'input', # after the date of the input folder the tar came from (out_date, passed by the driver) dt_out = out_date if (out_date is not None and exopt['out_date'] == 'input') else dt_rr sub_out = exopt['out_subdir'] + '/' if exopt['out_subdir'] else '' # STAR 2026-11-01: / when reprocessing out_dir = f's04_atmPhs/{sub_out}{dt_out:%Y-%m-%d}/' makedirs(out_dir, exist_ok=True) # TC: updated file format on 2026/04/07 (2 lines) # TC: updated timezone on 2026/04/08 (1 line) # TC: updated 2026/05/07 (1 line, add %S second) filestamp = f'{leoName}_{dt_ro[0]:%Y%m%d%H%M%S}_{dt_ro[-1]:%Y%m%d%H%M%S}_{gnss_ro}_A{ant_leo_ro:02d}_c{datetime.now().astimezone(ZoneInfo("UTC")):%Y%m%d%H%M%S}' #f'{leoName}.{dt_ro[0 if obj_leo_ro["ocd"] else -1]:%Y.%j.%H.%M}.{gnss_ro}' out_file = out_dir + f'RONA_L1A_{filestamp}_nes_ops.nc' #f'atmPhs_STAR_{filestamp}_v{version}.nc' run(['ncgen','-b','./template.cdl','-o',out_file]) from netCDF4 import Dataset with Dataset(out_file,'r+') as fid: fid['time'][...] = data_ro[:, 1] fid['caL1Snr'][...] = data_ro[:, s_ro[0]] fid['pL1Snr'][...] = data_ro[:, s_ro[0]] fid['pL2Snr'][...] = data_ro[:, s_ro[1]] #fid['xLeo'][...] = obj_leo_ro['pos'][:, 0] / 1E3 #fid['yLeo'][...] = obj_leo_ro['pos'][:, 1] / 1E3 #fid['zLeo'][...] = obj_leo_ro['pos'][:, 2] / 1E3 #fid['xdLeo'][...] = obj_leo_ro['vel'][:, 0] / 1E3 #fid['ydLeo'][...] = obj_leo_ro['vel'][:, 1] / 1E3 #fid['zdLeo'][...] = obj_leo_ro['vel'][:, 2] / 1E3 #fid['xGps'][...] = obj_gnss_ro['pos'][:, 0] / 1E3 #fid['yGps'][...] = obj_gnss_ro['pos'][:, 1] / 1E3 #fid['zGps'][...] = obj_gnss_ro['pos'][:, 2] / 1E3 #fid['xdGps'][...] = obj_gnss_ro['vel'][:, 0] / 1E3 #fid['ydGps'][...] = obj_gnss_ro['vel'][:, 1] / 1E3 #fid['zdGps'][...] = obj_gnss_ro['vel'][:, 2] / 1E3 fid['xLeo'][...] = pos_leo_tod[:, 0] / 1E3 #obj_leo_ro['pos'][:, 0] / 1E3 fid['yLeo'][...] = pos_leo_tod[:, 1] / 1E3 #obj_leo_ro['pos'][:, 1] / 1E3 fid['zLeo'][...] = pos_leo_tod[:, 2] / 1E3 #obj_leo_ro['pos'][:, 2] / 1E3 fid['xdLeo'][...] = vel_leo_tod[:, 0] / 1E3 #obj_leo_ro['vel'][:, 0] / 1E3 fid['ydLeo'][...] = vel_leo_tod[:, 1] / 1E3 #obj_leo_ro['vel'][:, 1] / 1E3 fid['zdLeo'][...] = vel_leo_tod[:, 2] / 1E3 #obj_leo_ro['vel'][:, 2] / 1E3 fid['xGps'][...] = pos_gnss_tod[:, 0] / 1E3 #obj_gnss_ro['pos'][:, 0] / 1E3 fid['yGps'][...] = pos_gnss_tod[:, 1] / 1E3 #obj_gnss_ro['pos'][:, 1] / 1E3 fid['zGps'][...] = pos_gnss_tod[:, 2] / 1E3 #obj_gnss_ro['pos'][:, 2] / 1E3 fid['xdGps'][...] = vel_gnss_tod[:, 0] / 1E3 #obj_gnss_ro['vel'][:, 0] / 1E3 fid['ydGps'][...] = vel_gnss_tod[:, 1] / 1E3 #obj_gnss_ro['vel'][:, 1] / 1E3 fid['zdGps'][...] = vel_gnss_tod[:, 2] / 1E3 #obj_gnss_ro['vel'][:, 2] / 1E3 fid['exL1'][...] = exL1 fid['exL2'][...] = exL2 fid['exLC'][...] = exLC fid['xmdl'][...] = data_ro[:, m_ro[0]] - data_ro[0, m_ro[0]] fid.setncatts({ 'processingHistory': 'STAR_L1A_L1B_V1.0', 'timeOfProcessing': datetime.now().strftime('%d-%b-%Y %H:%M:%S'), 'missionId': 'planetiq', 'trkStatus': 0, 'leoID': int(sub(r'[A-Z]',r'',leoName)), 'occID': 1, 'conID': 'G', 'occfreq1': f1_ro, 'occfreq2': f2_ro, 'Obs_type_L1': header_ro[l_ro[0]], 'Obs_type_L2': header_ro[l_ro[1]], 'setting': int(obj_leo_ro['ocd']), 'refsatId': gnss_ref if reffile is not None else 'none', 'occsatId': gnss_ro, 'year': dt_ros[0].year, 'month': dt_ros[0].month, 'day': dt_ros[0].day, 'hour': dt_ros[0].hour, 'minute': dt_ros[0].minute, 'second': dt_ros[0].second, 'dayOfYear': int(f'{dt_ros[0]:%j}'), 'startTime': data_ro[0,0], 'stopTime': data_ro[-1,0]+data_ro[-1,1], # TC: updated 2026/05/07 (1 line, change organization) 'center': 'NESDIS', #'center': 'NOAA/STAR, ESSIC/CISESS/UMD', 'fileStamp': filestamp, 'VERSION': f"Version {version}: {'Single' if exopt['clock'] == 'single_diff' else 'Zero'} Differencing, Corrected ECI Coordinates", # TC: updated 2026-10-27 (4 lines, add new attributes) 'tarfile': tarsource, 'software_version': '1.3', 'pod_method': pod_method, # 'DMF' if the orbit was replaced by the dynamic fit (pod_qc bit 4), else 'PlanetiQ' 'navbit': int(with_navbit), # STAR 2026-10-30: excess-phase processing options 'exph_clock': exopt['clock'] + ('' if exopt['clock'] == 'single_diff' else '/' + exopt['clk_interp']), 'exph_connect': connect_ro, # STAR 2026-11-01: LEO POD quality record (pod_qc bit 1: orbit resumed after a gap < exph_pod_gap_win s # before the occultation; bit 2: max 30 s mean |dP/dt - V| > exph_pod_pv_max m/s; bit 3 (4): orbit after # the gap replaced by DMF, exph_pod_fit; bit 4 (8): flag-2 stretch replaced by DMF, exph_pod_fit_pv); # -1 s = no gap in the orbit file 'pod_qc': podq['qc'], 'pod_gap_after_s': podq['gap_after_s'] if np.isfinite(podq['gap_after_s']) else -1.0, 'pod_gap_before_s': podq['gap_before_s'] if np.isfinite(podq['gap_before_s']) else -1.0, 'pod_pv_max': round(podq['pv_max'], 5), 'pod_fit': podfit['status'], 'pod_fit_rms_m': round(podfit['rms'], 4), 'pod_fit_maxdiff_m': round(podfit['maxdiff'], 3), 'pod_fit_pv': pvfit['status'], 'pod_fit_pv_span_s': round(float(pvfit['span_s']), 1), 'pod_fit_pv_rms_m': round(float(pvfit['rms']), 4), 'pod_fit_pv_maxdiff_m': round(float(pvfit['maxdiff']), 3), 'pod_pv_max_after': round(float(podq['pv_after']), 5), # STAR 2026-11-04: L1A quality flag: bad = "1" when the POD is still bad after the fixes (reason in errstr: # POD gap_unfitted / pv_after / fit_rms); ucar2ropp reads it and ropp_pp keeps it with qc_l1a_bad = .true. 'bad': '1' if podq['bad'] else '0', 'errstr': ('POD ' + podq['bad_reason']) if podq['bad'] else 'OK' }) good = True if not exopt['inline_ropp']: # STAR 2026-11-01: exph_inline_ropp = False skips the in-line L1A -> L1C / BUFR break #--------------------------------------------------------- ## Xinjia modified, for RFSI/ROPP and BUFR atmPhs_ropp_dir = f"{rfsi_dir}/atmPhs_ropp{suffix}/{sub_out}{dt_out:%Y-%m-%d}/" atmPrf_dir = f"{rfsi_dir}/atmPrf{suffix}/{sub_out}{dt_out:%Y-%m-%d}/" bufr_dir = f"{rfsi_dir}/bfrPrf{suffix}/{sub_out}{dt_out:%Y-%m-%d}/" makedirs(atmPhs_ropp_dir, exist_ok=True) makedirs(atmPrf_dir, exist_ok=True) makedirs(bufr_dir, exist_ok=True) #out_file = out_dir + f'RONA_L1A_{filestamp}_nes_ops.nc' #f'atmPhs_STAR_{filestamp}_v{version}.nc' atmPhs_ropp_file = atmPhs_ropp_dir+ f'RONA_L1A_{filestamp}_nes_ops.nc' atmPrf_file = atmPrf_dir + f'RONA_L1C_{filestamp}_nes_ops.nc' bufr_file = bufr_dir + f'RONA_L1C_{filestamp}_nes_ops.bin' try: run([ucar2ropp,out_file,'-o',atmPhs_ropp_file]) except Exception as e: log_message(f'- ucar2ropp failed for {out_file}: {e}') break try: # STAR 2026-11-01: configuration file and timeout from configuration.py (delivered: cosmic_pp_rfsi_planetiq.cf, 20 s) run([RFSI_occ_tool,atmPhs_ropp_file,'-o',atmPrf_file,'-c',exopt['ropp_cf'],'-d','-full'],check=True,timeout=exopt['ropp_timeout']) except Exception as e: log_message(f'- RFSI_occ_tool failed for {out_file}: {e}') if path.exists(atmPrf_file) and path.getsize(atmPrf_file) == 0: remove(atmPrf_file) log_message(f"Removed empty atmPrf file: {atmPrf_file}") break #--------------------------------------------------------- # Remove empty atmPrf files #--------------------------------------------------------- try: if path.exists(atmPrf_file) and path.getsize(atmPrf_file) == 0: remove(atmPrf_file) log_message(f"Removed empty atmPrf file: {atmPrf_file}") break except Exception as e: log_message(f"WARNING: failed checking/removing empty atmPrf file {atmPrf_file}: {e}") #run([ropp2bufr,atmPrf_file,'-o',bufr_file,'-p',env["ROPP_THIN"]],env=env,check=True) import sys run([sys.executable, "../ropp_nc2bufr.py", atmPrf_file, "-o", bufr_file, "--heading", "nco", "--arh", "--corrected", "generic"], env=env, check=True) #~/miniconda3/bin/python ropp_nc2bufr.py ${rfsi_file} -o ${netcdf_path_rfsi}/bfrPrf${sub_version}/${thisDATE}/${bufr_file} --heading nco --arh --corrected generic #--------------------------------------------------------- # Remove empty BUFR files #--------------------------------------------------------- try: if path.exists(bufr_file) and path.getsize(bufr_file) == 0: remove(bufr_file) log_message(f"Removed empty BUFR file: {bufr_file}") break except Exception as e: log_message(f"WARNING: failed checking/removing empty BUFR file {bufr_file}: {e}") #--------------------------------------------------------- # Remove temporary ROPP atmPhs file (no longer needed) #--------------------------------------------------------- try: if path.exists(atmPhs_ropp_file): remove(atmPhs_ropp_file) except Exception as e: log_message(f"WARNING: failed removing {atmPhs_ropp_file}: {e}") ## Xinjia modified, for RFSI/ROPP and BUFR #--------------------------------------------------------- break if not good: raise LookupError('Could not find a good pod ref profile!')