diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index bfba5883..c852e18f 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -319,6 +319,52 @@ jobs: output/ppp_static_products_comparison.png if-no-files-found: ignore + madoca-parity: + needs: [changes, hygiene] + if: needs.changes.outputs.run_heavy == 'true' + runs-on: ubuntu-latest + timeout-minutes: 30 + steps: + - uses: actions/checkout@v4 + + - name: Checkout MADOCALIB + uses: actions/checkout@v4 + with: + repository: QZSS-Strategy-Office/madocalib + ref: 0089f7dc97e8e2ba283a40be2edf4b73a140df6c + path: external/madocalib + + - name: Install dependencies + run: | + sudo apt-get update + sudo apt-get install -y libeigen3-dev libgtest-dev cmake build-essential python3-dev pybind11-dev + + - name: Configure CMake with MADOCALIB parity link + run: | + cmake -B build-madoca-parity \ + -DCMAKE_BUILD_TYPE=Release \ + -DMADOCALIB_PARITY_LINK=ON \ + -DMADOCALIB_ROOT_DIR="${GITHUB_WORKSPACE}/external/madocalib" + + - name: Build MadocaParity tests + run: cmake --build build-madoca-parity --config Release --parallel --target run_tests + + - name: Run MadocaParity + working-directory: build-madoca-parity + run: ./tests/run_tests --gtest_filter='*MadocaParity*' + + - name: Upload MadocaParity diagnostics + if: failure() + uses: actions/upload-artifact@v4 + with: + name: madoca-parity-diagnostics + path: | + build-madoca-parity/CMakeCache.txt + build-madoca-parity/Testing/Temporary/LastTest.log + build-madoca-parity/Testing/Temporary/LastTestsFailed.log + build-madoca-parity/CMakeFiles/CMakeConfigureLog.yaml + if-no-files-found: ignore + static-analysis: needs: [changes, hygiene] if: needs.changes.outputs.run_heavy == 'true' diff --git a/.gitignore b/.gitignore index 749d964e..536aa9ac 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,7 @@ # Build artifacts build/ build_local/ +build_madoca/ bin/ *.a *.o diff --git a/CMakeLists.txt b/CMakeLists.txt index bcba7635..658572b5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -22,11 +22,17 @@ find_package(nav_msgs QUIET) find_package(std_msgs QUIET) find_package(OpenCV QUIET COMPONENTS core imgproc imgcodecs) +option(MADOCALIB_PARITY_LINK "Link external MADOCALIB C sources for opt-in parity oracle tests" OFF) +set(MADOCALIB_ROOT_DIR + "$ENV{MADOCALIB_ROOT_DIR}" + CACHE PATH "Path to the MADOCALIB checkout root") + # List source files for the main library set(GNSS_SOURCES src/gnss.cpp src/io/ntrip.cpp src/io/rinex.cpp + src/io/rinex4.cpp src/io/rtcm.cpp src/io/rtcm_stream.cpp src/io/ubx.cpp @@ -40,6 +46,10 @@ set(GNSS_SOURCES src/algorithms/ppp_osr.cpp src/algorithms/ppp_wlnl.cpp src/algorithms/ppp_clas_epoch.cpp + src/algorithms/madoca_core.cpp + src/algorithms/madoca_parity.cpp + src/algorithms/madocalib_oracle.cpp + src/algorithms/madocalib_bridge.cpp src/algorithms/ppp.cpp ) @@ -66,6 +76,78 @@ target_include_directories(gnss_lib PUBLIC ) target_link_libraries(gnss_lib PUBLIC Eigen3::Eigen) +if(MADOCALIB_PARITY_LINK) + set(MADOCALIB_SRC_DIR "${MADOCALIB_ROOT_DIR}/src") + if(NOT EXISTS "${MADOCALIB_SRC_DIR}/rtklib.h") + message(FATAL_ERROR "MADOCALIB source tree not found: ${MADOCALIB_SRC_DIR}") + endif() + + add_library(madocalib STATIC + ${MADOCALIB_SRC_DIR}/rtkcmn.c + ${MADOCALIB_SRC_DIR}/rinex.c + ${MADOCALIB_SRC_DIR}/rtkpos.c + ${MADOCALIB_SRC_DIR}/postpos.c + ${MADOCALIB_SRC_DIR}/solution.c + ${MADOCALIB_SRC_DIR}/lambda.c + ${MADOCALIB_SRC_DIR}/geoid.c + ${MADOCALIB_SRC_DIR}/sbas.c + ${MADOCALIB_SRC_DIR}/preceph.c + ${MADOCALIB_SRC_DIR}/pntpos.c + ${MADOCALIB_SRC_DIR}/ephemeris.c + ${MADOCALIB_SRC_DIR}/options.c + ${MADOCALIB_SRC_DIR}/ppp.c + ${MADOCALIB_SRC_DIR}/ppp_ar.c + ${MADOCALIB_SRC_DIR}/ppp_iono.c + ${MADOCALIB_SRC_DIR}/rtcm.c + ${MADOCALIB_SRC_DIR}/rtcm2.c + ${MADOCALIB_SRC_DIR}/rtcm3.c + ${MADOCALIB_SRC_DIR}/rtcm3e.c + ${MADOCALIB_SRC_DIR}/ionex.c + ${MADOCALIB_SRC_DIR}/tides.c + ${MADOCALIB_SRC_DIR}/mdccssr.c + ${MADOCALIB_SRC_DIR}/mdciono.c + ) + target_include_directories(madocalib PRIVATE ${MADOCALIB_SRC_DIR}) + target_compile_definitions(madocalib + PRIVATE + TRACE + ENAGLO + ENAQZS + ENAGAL + ENACMP + ENAIRN + NFREQ=5 + NEXOBS=5 + ) + target_link_libraries(madocalib PRIVATE m Threads::Threads) + if(UNIX AND NOT APPLE) + target_link_libraries(madocalib PRIVATE rt) + endif() + + target_include_directories(gnss_lib PRIVATE ${MADOCALIB_SRC_DIR}) + target_compile_definitions(gnss_lib + PUBLIC + GNSSPP_HAS_MADOCALIB_BRIDGE=1 + GNSSPP_HAS_MADOCALIB_ORACLE=1 + PRIVATE + TRACE + ENAGLO + ENAQZS + ENAGAL + ENACMP + ENAIRN + NFREQ=5 + NEXOBS=5 + "GNSSPP_MADOCALIB_ROOT=\"${MADOCALIB_ROOT_DIR}\"" + ) + target_link_libraries(gnss_lib PUBLIC madocalib) +else() + target_compile_definitions(gnss_lib + PUBLIC + GNSSPP_HAS_MADOCALIB_BRIDGE=0 + GNSSPP_HAS_MADOCALIB_ORACLE=0) +endif() + # Create the no-optimization library add_library(gnss_lib_noopt STATIC ${GNSS_SOURCES_NO_OPTIMIZATION}) target_compile_options(gnss_lib_noopt PRIVATE -O0) @@ -75,8 +157,13 @@ target_include_directories(gnss_lib_noopt PUBLIC ) target_link_libraries(gnss_lib_noopt PUBLIC Eigen3::Eigen) +set(GNSS_INSTALL_TARGETS gnss_lib gnss_lib_noopt) +if(MADOCALIB_PARITY_LINK) + list(APPEND GNSS_INSTALL_TARGETS madocalib) +endif() + install( - TARGETS gnss_lib gnss_lib_noopt + TARGETS ${GNSS_INSTALL_TARGETS} EXPORT libgnssppTargets ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} diff --git a/apps/gnss_clas_ppp.py b/apps/gnss_clas_ppp.py index 974a21a2..0345fe54 100644 --- a/apps/gnss_clas_ppp.py +++ b/apps/gnss_clas_ppp.py @@ -22,10 +22,12 @@ ROOT_DIR = Path(__file__).resolve().parent.parent PPP_STATUSES = {5, 6} DEFAULT_SERIAL_BAUD = 115200 +DEFAULT_CLAS_FILTER_ITERATIONS = 1 +DEFAULT_QZSS_WINDOW_MARGIN_SECONDS = 90.0 COMPACT_FILE_COLUMNS = ( "week,tow,system,prn,dx,dy,dz,dclock_m" "[,high_rate_clock_m][,ura_sigma_m=][,cbias:=...][,pbias:=...]" - "[,bias_network_id=][,atmos_=...]" + "[,pdi:=...][,bias_network_id=][,atmos_=...]" ) @@ -39,6 +41,15 @@ def parse_args() -> argparse.Namespace: ) parser.add_argument("--obs", type=Path, required=True, help="Rover observation RINEX file.") parser.add_argument("--nav", type=Path, required=True, help="Broadcast navigation RINEX file.") + parser.add_argument( + "--navsys", + type=int, + default=None, + help=( + "RTKLIB navsys mask passed to PPP observation admission. " + "Default for --profile clas is 25 (GPS+Galileo+QZSS); use 0 for all observed systems." + ), + ) parser.add_argument( "--ssr-rtcm", default=None, @@ -57,12 +68,30 @@ def parse_args() -> argparse.Namespace: default=None, help="Direct QZSS L6 frame source: local file or serial:// URI.", ) + parser.add_argument( + "--expanded-ssr", + type=Path, + default=None, + help="Already-expanded sampled SSR CSV passed directly to native PPP.", + ) parser.add_argument( "--qzss-gps-week", type=int, default=None, help="Optional GPS week override used when decoding raw QZSS L6 corrections.", ) + parser.add_argument( + "--qzss-expanded-cache", + type=Path, + default=None, + help="Optional expanded SSR CSV cache path for --qzss-l6; reused when it already exists.", + ) + parser.add_argument( + "--qzss-window-margin-seconds", + type=float, + default=DEFAULT_QZSS_WINDOW_MARGIN_SECONDS, + help="Margin around the observation TOW window when expanding raw QZSS L6 for --max-epochs.", + ) parser.add_argument( "--compact-flush-policy", choices=qzss_l6_info.COMPACT_SSR_FLUSH_POLICIES, @@ -131,8 +160,18 @@ def parse_args() -> argparse.Namespace: ) parser.add_argument("--out", type=Path, required=True, help="Output PPP .pos file.") parser.add_argument("--summary-json", type=Path, default=None, help="Optional summary JSON path.") + parser.add_argument("--ppp-filter-log", type=Path, default=None, help="Optional native PPP filter iteration CSV path.") + parser.add_argument("--ppp-residual-log", type=Path, default=None, help="Optional native PPP residual CSV path.") + parser.add_argument("--ppp-state-log", type=Path, default=None, help="Optional native PPP state log path.") + parser.add_argument("--ppp-correction-log", type=Path, default=None, help="Optional native PPP correction diagnostic CSV path.") parser.add_argument("--sp3", type=Path, default=None, help="Optional precise SP3 file.") parser.add_argument("--clk", type=Path, default=None, help="Optional precise CLK file.") + parser.add_argument("--antex", type=Path, default=None, help="Optional ANTEX file for PPP.") + parser.add_argument( + "--receiver-antenna-type", + default=None, + help="Override the RINEX receiver antenna type for ANTEX lookup.", + ) parser.add_argument("--kml", type=Path, default=None, help="Optional KML output path.") parser.add_argument("--max-epochs", type=int, default=0, help="Stop after N epochs.") parser.add_argument( @@ -141,6 +180,15 @@ def parse_args() -> argparse.Namespace: default=None, help="Minimum epochs before PPP convergence/AR checks.", ) + parser.add_argument( + "--filter-iterations", + type=int, + default=None, + help=( + "Override native PPP measurement update iterations per epoch. " + "Default for --profile clas is 1; MADOCA leaves native PPP default." + ), + ) parser.add_argument( "--ssr-step-seconds", type=float, @@ -157,6 +205,157 @@ def parse_args() -> argparse.Namespace: default=3.0, help="AR ratio threshold when ambiguity fixing is enabled.", ) + parser.add_argument( + "--wlnl-wl-max-fractional", + type=float, + default=None, + help="Max MW wide-lane fractional cycle before rejecting WL fix.", + ) + parser.add_argument( + "--native-pntpos-parity-seed", + dest="native_pntpos_parity_seed", + action="store_true", + default=None, + help="Use native SPP settings that mirror pntpos behavior for PPP seed.", + ) + parser.add_argument( + "--no-native-pntpos-parity-seed", + dest="native_pntpos_parity_seed", + action="store_false", + help="Disable native pntpos-parity PPP seeding.", + ) + parser.add_argument( + "--clas-ppc-profile", + action="store_true", + help=( + "Use the native PPC kinematic CLAS profile " + "(CLAS OSR, residual ionosphere, 1 m anchor, 0.25 m^2 iono prior, " + "2x code variance, kinematic reseed, " + "pntpos-parity seed)." + ), + ) + parser.add_argument( + "--clas-anchor-sigma", + type=float, + default=None, + help="Override the native CLAS OSR SPP anchor sigma in meters.", + ) + parser.add_argument( + "--clas-code-variance-scale", + type=float, + default=None, + help="Override the native CLAS OSR code variance multiplier.", + ) + parser.add_argument( + "--clas-phase-variance", + type=float, + default=None, + help="Override the native CLAS OSR carrier-phase base variance in m^2.", + ) + parser.add_argument( + "--clas-phase-continuity", + choices=[ + "full-repair", + "sis-continuity-only", + "repair-only", + "raw-phase-bias", + "no-phase-bias", + ], + default=None, + help="Override the native CLAS phase-bias continuity policy.", + ) + parser.add_argument( + "--clas-phase-bias-values", + choices=["full", "phase-bias-only", "compensation-only"], + default=None, + help="Override which native CLAS phase-bias correction values are applied.", + ) + parser.add_argument( + "--clas-phase-bias-reference-time", + choices=["phase-bias-reference", "clock-reference", "observation-epoch"], + default=None, + help="Override the native CLAS phase-bias continuity reference epoch.", + ) + parser.add_argument( + "--clas-iono-prior-variance", + type=float, + default=None, + help="Override the native CLAS residual-ionosphere prior variance in m^2.", + ) + parser.add_argument( + "--clas-code-outlier-sigma-scale", + type=float, + default=None, + help="Override the native CLAS code-only outlier sigma gate; 0 disables.", + ) + parser.add_argument( + "--clas-code-outlier-min-residual", + type=float, + default=None, + help="Override the native CLAS code-only outlier absolute residual floor in meters.", + ) + parser.add_argument( + "--clas-phase-outlier-sigma-scale", + type=float, + default=None, + help="Override the native CLAS phase-only outlier sigma gate; 0 disables.", + ) + parser.add_argument( + "--clas-phase-outlier-inflated-variance", + type=float, + default=None, + help="Override native CLAS phase-only inflated measurement variance in m^2.", + ) + parser.add_argument( + "--clas-prior-outlier-sigma-scale", + type=float, + default=None, + help="Override the native CLAS prior-constraint outlier sigma gate; 0 disables.", + ) + parser.add_argument( + "--clas-prior-outlier-min-residual", + type=float, + default=None, + help="Override the native CLAS prior-constraint outlier absolute residual floor in meters.", + ) + parser.add_argument( + "--clas-phase-outlier-min-residual", + type=float, + default=None, + help="Override the native CLAS phase-only outlier absolute residual floor in meters.", + ) + parser.add_argument( + "--clas-reset-phase-ambiguity-on-outlier-inflation", + action="store_true", + help="Reset native CLAS phase ambiguities touched by inflated phase rows.", + ) + parser.add_argument( + "--clas-reset-phase-ambiguity-outlier-min-rows", + type=int, + default=None, + help="Minimum inflated CLAS phase rows before ambiguity reset is applied.", + ) + parser.add_argument( + "--clas-reset-phase-ambiguity-outlier-min-residual", + type=float, + default=None, + help="Minimum inflated CLAS phase residual for ambiguity reset membership.", + ) + parser.add_argument( + "--clas-wlnl-fixed-position-max-shift", + type=float, + default=None, + help=( + "Reject native CLAS WLNL fixed-position replacement when it shifts " + "farther than this from the float solution; 0 disables." + ), + ) + parser.add_argument( + "--claslib-pntpos-seed", + dest="native_pntpos_parity_seed", + action="store_true", + help="Deprecated alias for --native-pntpos-parity-seed.", + ) parser.add_argument( "--no-estimate-troposphere", dest="estimate_troposphere", @@ -169,7 +368,7 @@ def parse_args() -> argparse.Namespace: action="store_true", help="Enable zenith troposphere estimation.", ) - parser.set_defaults(estimate_troposphere=True) + parser.set_defaults(estimate_troposphere=None) parser.add_argument( "--require-valid-epochs-min", type=int, @@ -206,6 +405,76 @@ def parse_args() -> argparse.Namespace: default=None, help="Fail if PPP-applied atmospheric ionosphere corrections are below this count.", ) + parser.add_argument( + "--enable-wlnl-par", + dest="enable_wlnl_par", + action="store_true", + default=None, + help="Enable WLNL Partial AR (greedy worst-|frac| exclusion).", + ) + parser.add_argument( + "--no-enable-wlnl-par", + dest="enable_wlnl_par", + action="store_false", + help="Disable WLNL Partial AR.", + ) + parser.add_argument( + "--clas-spp-clock-overwrite", + dest="clas_spp_clock_overwrite", + action="store_true", + default=None, + help="Overwrite KF receiver-clock state with SPP every kinematic CLAS-OSR epoch (default on).", + ) + parser.add_argument( + "--no-clas-spp-clock-overwrite", + dest="clas_spp_clock_overwrite", + action="store_false", + help="Skip the SPP clock overwrite on the CLAS-OSR path; let the KF track the receiver clock.", + ) + parser.add_argument( + "--clas-kf-clock-seed-variance", + type=float, + default=None, + help="Initial KF clock variance when --no-clas-spp-clock-overwrite is set.", + ) + parser.add_argument( + "--ppp-process-noise-clock", + type=float, + default=None, + help="Override KF receiver-clock process noise (m^2/s).", + ) + parser.add_argument( + "--process-noise-troposphere", + type=float, + default=None, + help="Override troposphere process noise.", + ) + parser.add_argument( + "--reset-clock-to-spp", + dest="reset_clock_to_spp", + action="store_true", + default=None, + help="Force KF receiver-clock state to SPP every epoch (default on).", + ) + parser.add_argument( + "--no-reset-clock-to-spp", + dest="reset_clock_to_spp", + action="store_false", + help="Skip SPP overwrite of KF receiver-clock state on the standard PPP path.", + ) + parser.add_argument( + "--spp-seed-iono-free-code", + dest="spp_seed_iono_free_code", + action="store_true", + default=None, + help="Use ionosphere-free L1+L2 code combination in the SPP seed.", + ) + parser.add_argument( + "--no-spp-seed-iono-free-code", + dest="spp_seed_iono_free_code", + action="store_false", + help="Force single-frequency L1 code in the SPP seed (default).", + ) return parser.parse_args() @@ -223,17 +492,45 @@ def classify_transport(source: str) -> str: def classify_encoding(args: argparse.Namespace) -> str: + if args.expanded_ssr: + return "expanded" if args.qzss_l6: return "qzss_l6" return "compact" if args.compact_ssr else "rtcm" def selected_correction_source(args: argparse.Namespace) -> str: + if args.expanded_ssr: + return str(args.expanded_ssr) if args.qzss_l6: return args.qzss_l6 return args.compact_ssr if args.compact_ssr else args.ssr_rtcm +def effective_navsys_mask(args: argparse.Namespace) -> int: + if args.navsys is not None: + return args.navsys + return 25 if args.profile == "clas" or args.clas_ppc_profile else 0 + + +def effective_estimate_troposphere(args: argparse.Namespace) -> bool: + if args.estimate_troposphere is not None: + return bool(args.estimate_troposphere) + return not bool(args.clas_ppc_profile) + + +def effective_native_pntpos_parity_seed(args: argparse.Namespace) -> bool: + if args.native_pntpos_parity_seed is not None: + return bool(args.native_pntpos_parity_seed) + return bool(args.clas_ppc_profile) + + +def effective_filter_iterations(args: argparse.Namespace) -> int | None: + if args.filter_iterations is not None: + return args.filter_iterations + return DEFAULT_CLAS_FILTER_ITERATIONS if args.profile == "clas" else None + + def parse_serial_path(path: str) -> tuple[str, int]: raw = path[len("serial://") :] if path.startswith("serial://") else path if "?baud=" in raw: @@ -349,6 +646,7 @@ def expand_compact_ssr_text(text: str, output_path: Path) -> dict[str, object]: ura_sigma_token = None code_bias_tokens: list[str] = [] phase_bias_tokens: list[str] = [] + phase_discontinuity_tokens: list[str] = [] bias_network_tokens: list[str] = [] atmos_tokens: list[str] = [] extras = columns[8:] @@ -365,6 +663,9 @@ def expand_compact_ssr_text(text: str, output_path: Path) -> dict[str, object]: if token.startswith("pbias:") and "=" in token: phase_bias_tokens.append(token) continue + if token.startswith("pdi:") and "=" in token: + phase_discontinuity_tokens.append(token) + continue if token.startswith("bias_network_id="): bias_network_tokens.append(token) continue @@ -391,6 +692,7 @@ def expand_compact_ssr_text(text: str, output_path: Path) -> dict[str, object]: output_tokens.append(ura_sigma_token) output_tokens.extend(code_bias_tokens) output_tokens.extend(phase_bias_tokens) + output_tokens.extend(phase_discontinuity_tokens) output_tokens.extend(bias_network_tokens) output_tokens.extend(atmos_tokens) handle.write(",".join(output_tokens) + "\n") @@ -404,8 +706,34 @@ def expand_compact_ssr_text(text: str, output_path: Path) -> dict[str, object]: } -def infer_gps_week_from_obs(obs_path: Path) -> int: +def gps_week_tow_from_datetime( + year: int, + month: int, + day: int, + hour: int, + minute: int, + second: float, +) -> tuple[int, float]: + whole_seconds = int(second) + microseconds = int(round((second - whole_seconds) * 1_000_000)) + epoch = dt.datetime( + year, + month, + day, + hour, + minute, + whole_seconds, + microseconds, + tzinfo=dt.timezone.utc, + ) gps_epoch = dt.datetime(1980, 1, 6, tzinfo=dt.timezone.utc) + delta_seconds = (epoch - gps_epoch).total_seconds() + week = int(delta_seconds // (7 * 24 * 60 * 60)) + tow = round(delta_seconds - week * 7 * 24 * 60 * 60, 7) + return week, tow + + +def infer_gps_week_from_obs(obs_path: Path) -> int: with obs_path.open(encoding="ascii", errors="ignore") as handle: for line in handle: if not line.startswith(">"): @@ -413,34 +741,138 @@ def infer_gps_week_from_obs(obs_path: Path) -> int: parts = line[1:].split() if len(parts) < 6: continue - year = int(parts[0]) - month = int(parts[1]) - day = int(parts[2]) - hour = int(parts[3]) - minute = int(parts[4]) - second = float(parts[5]) - whole_seconds = int(second) - microseconds = int(round((second - whole_seconds) * 1_000_000)) - epoch = dt.datetime( - year, - month, - day, - hour, - minute, - whole_seconds, - microseconds, - tzinfo=dt.timezone.utc, + week, _tow = gps_week_tow_from_datetime( + int(parts[0]), + int(parts[1]), + int(parts[2]), + int(parts[3]), + int(parts[4]), + float(parts[5]), ) - delta = epoch - gps_epoch - return delta.days // 7 + return week raise SystemExit(f"Could not infer GPS week from observation file: {obs_path}") +def observation_tow_window(obs_path: Path, max_epochs: int) -> tuple[int, float, float] | None: + if max_epochs <= 0: + return None + epoch_tows: list[tuple[int, float]] = [] + in_header = True + with obs_path.open(encoding="ascii", errors="ignore") as handle: + for line in handle: + if in_header: + if "END OF HEADER" in line: + in_header = False + continue + if not line.startswith(">"): + continue + parts = line[1:].split() + if len(parts) < 6: + continue + epoch_tows.append( + gps_week_tow_from_datetime( + int(parts[0]), + int(parts[1]), + int(parts[2]), + int(parts[3]), + int(parts[4]), + float(parts[5]), + ) + ) + if len(epoch_tows) >= max_epochs: + break + if not epoch_tows: + return None + weeks = {week for week, _tow in epoch_tows} + if len(weeks) != 1: + return None + week = epoch_tows[0][0] + tows = [tow for _week, tow in epoch_tows] + return week, min(tows), max(tows) + + +def correction_in_window( + correction: qzss_l6_info.CompactSSRCorrection, + window: tuple[int, float, float] | None, + margin_seconds: float, +) -> bool: + if window is None: + return True + week, start_tow, end_tow = window + return ( + correction.week == week + and correction.tow >= max(0.0, start_tow - margin_seconds) + and correction.tow <= min(604800.0, end_tow + margin_seconds) + ) + + +def expanded_correction_tokens(correction: qzss_l6_info.CompactSSRCorrection) -> list[str]: + system = normalize_system_token(correction.system) + tokens = [ + str(correction.week), + f"{correction.tow:.3f}", + f"{system}{correction.prn:02d}", + f"{correction.dx:.6f}", + f"{correction.dy:.6f}", + f"{correction.dz:.6f}", + f"{correction.dclock_m + correction.high_rate_clock_m:.6f}", + ] + if correction.ura_sigma_m is not None: + tokens.append(f"ura_sigma_m={correction.ura_sigma_m:.6f}") + code_bias_m = correction.code_bias_m or {} + phase_bias_m = correction.phase_bias_m or {} + phase_discontinuity = correction.phase_discontinuity or {} + for signal_id in sorted(code_bias_m): + tokens.append(f"cbias:{signal_id}={code_bias_m[signal_id]:.6f}") + for signal_id in sorted(phase_bias_m): + tokens.append(f"pbias:{signal_id}={phase_bias_m[signal_id]:.6f}") + for signal_id in sorted(phase_discontinuity): + tokens.append(f"pdi:{signal_id}={phase_discontinuity[signal_id]}") + if correction.bias_network_id is not None: + tokens.append(f"bias_network_id={correction.bias_network_id}") + if correction.atmos_network_id is not None: + tokens.append(f"atmos_network_id={correction.atmos_network_id}") + if correction.atmos_trop_avail is not None: + tokens.append(f"atmos_trop_avail={correction.atmos_trop_avail}") + if correction.atmos_stec_avail is not None: + tokens.append(f"atmos_stec_avail={correction.atmos_stec_avail}") + if correction.atmos_grid_count is not None: + tokens.append(f"atmos_grid_count={correction.atmos_grid_count}") + if correction.atmos_selected_satellites is not None: + tokens.append(f"atmos_selected_satellites={correction.atmos_selected_satellites}") + atmos_tokens = correction.atmos_tokens or {} + for key in sorted(atmos_tokens): + if key in { + "atmos_network_id", + "atmos_trop_avail", + "atmos_stec_avail", + "atmos_grid_count", + "atmos_selected_satellites", + }: + continue + tokens.append(f"{key}={atmos_tokens[key]}") + return tokens + + +def write_expanded_corrections( + path: Path, + corrections: list[qzss_l6_info.CompactSSRCorrection], +) -> int: + path.parent.mkdir(parents=True, exist_ok=True) + with path.open("w", encoding="ascii") as handle: + handle.write("# week,tow,sat,dx,dy,dz,dclock_m\n") + for correction in corrections: + handle.write(",".join(expanded_correction_tokens(correction)) + "\n") + return len(corrections) + + def expand_qzss_l6_source( source: str, gps_week: int, output_path: Path, *, + correction_window: tuple[int, float, float] | None = None, + correction_window_margin_seconds: float = DEFAULT_QZSS_WINDOW_MARGIN_SECONDS, compact_flush_policy: str = qzss_l6_info.COMPACT_SSR_FLUSH_POLICY_LAG_TOLERANT, compact_atmos_merge_policy: str = qzss_l6_info.COMPACT_ATMOS_MERGE_POLICY_STEC_COEFF_CARRY, compact_atmos_subtype_merge_policy: str = qzss_l6_info.COMPACT_ATMOS_SUBTYPE_MERGE_POLICY_UNION, @@ -469,9 +901,13 @@ def expand_qzss_l6_source( bias_row_materialization_policy=compact_bias_row_materialization, row_construction_policy=compact_row_construction_policy, ) - compact_source = output_path.parent / "qzss_l6_compact.csv" - qzss_l6_info.write_compact_corrections(compact_source, corrections) - expand_compact_ssr_text(compact_source.read_text(encoding="ascii"), output_path) + decoded_rows = len(corrections) + corrections = [ + correction + for correction in corrections + if correction_in_window(correction, correction_window, correction_window_margin_seconds) + ] + write_expanded_corrections(output_path, corrections) atmos_messages = sum( 1 for message in messages @@ -491,10 +927,21 @@ def expand_qzss_l6_source( "subframes": len(subframes), "messages": len(messages), "rows_written": len(corrections), + "rows_decoded": decoded_rows, "atmos_messages": atmos_messages, "atmos_rows": atmos_rows, - "compact_csv": str(compact_source), + "compact_csv": None, "expanded_csv": str(output_path), + "correction_window": ( + None + if correction_window is None + else { + "week": correction_window[0], + "start_tow": correction_window[1], + "end_tow": correction_window[2], + "margin_seconds": correction_window_margin_seconds, + } + ), "compact_atmos_merge_policy": compact_atmos_merge_policy, "compact_atmos_subtype_merge_policy": compact_atmos_subtype_merge_policy, "compact_phase_bias_merge_policy": compact_phase_bias_merge_policy, @@ -542,6 +989,12 @@ def rounded(value: float) -> float: return round(value, 6) +def optional_int(value: object) -> int | None: + if value is None: + return None + return int(value) + + def _parse_ppp_summary_counts(ppp_stdout: str) -> dict[str, int | float]: parsed = { "ppp_atmospheric_trop_corrections": 0, @@ -586,12 +1039,20 @@ def build_summary_payload( "correction_encoding": classify_encoding(args), "obs": str(args.obs), "nav": str(args.nav), + "navsys_mask": effective_navsys_mask(args), + "navsys_all_observed": effective_navsys_mask(args) == 0, + "estimate_troposphere": effective_estimate_troposphere(args), + "native_pntpos_parity_seed": effective_native_pntpos_parity_seed(args), + "filter_iterations": effective_filter_iterations(args), "sp3": str(args.sp3) if args.sp3 is not None else None, "clk": str(args.clk) if args.clk is not None else None, "ssr_rtcm": args.ssr_rtcm, "compact_ssr": args.compact_ssr, "qzss_l6": args.qzss_l6, + "expanded_ssr": str(args.expanded_ssr) if args.expanded_ssr is not None else None, + "qzss_expanded_cache": str(args.qzss_expanded_cache) if args.qzss_expanded_cache is not None else None, "solution_pos": str(args.out), + "ppp_correction_log": str(args.ppp_correction_log) if args.ppp_correction_log is not None else None, "epochs": len(records), "ppp_float_epochs": ppp_float_epochs, "ppp_fixed_epochs": ppp_fixed_epochs, @@ -600,6 +1061,12 @@ def build_summary_payload( "mean_satellites": rounded(mean_satellites), "ambiguity_resolution_enabled": bool(args.enable_ar), "ar_ratio_threshold": rounded(args.ar_ratio_threshold), + "clas_phase_continuity": args.clas_phase_continuity, + "clas_phase_bias_values": args.clas_phase_bias_values, + "clas_phase_bias_reference_time": args.clas_phase_bias_reference_time, + "wlnl_wl_max_fractional": rounded(args.wlnl_wl_max_fractional) + if args.wlnl_wl_max_fractional is not None + else None, "mode": "kinematic" if args.kinematic else "static", "compact_atmos_merge_policy": args.compact_atmos_merge_policy, "compact_atmos_subtype_merge_policy": args.compact_atmos_subtype_merge_policy, @@ -619,8 +1086,20 @@ def build_summary_payload( "ppp_atmospheric_ionosphere_meters": 0.0, } if compact_summary is not None: - payload["atmos_messages"] = int(compact_summary.get("atmos_messages", 0)) - payload["atmos_rows"] = int(compact_summary.get("atmos_rows", 0)) + if not compact_summary.get("cache_hit", False): + payload["atmos_messages"] = int(compact_summary.get("atmos_messages", 0)) + payload["atmos_rows"] = int(compact_summary.get("atmos_rows", 0)) + expanded_csv = compact_summary.get("expanded_csv") + payload["qzss_expanded_csv"] = str(expanded_csv) if expanded_csv is not None else None + payload["qzss_expanded_csv_size_bytes"] = ( + Path(str(expanded_csv)).stat().st_size + if expanded_csv is not None and Path(str(expanded_csv)).exists() + else None + ) + payload["qzss_rows_written"] = optional_int(compact_summary.get("rows_written")) + payload["qzss_rows_decoded"] = optional_int(compact_summary.get("rows_decoded")) + payload["qzss_cache_hit"] = bool(compact_summary.get("cache_hit", False)) + payload["qzss_correction_window"] = compact_summary.get("correction_window") if ppp_summary_counts is not None: payload["ppp_atmospheric_trop_corrections"] = int( ppp_summary_counts.get("ppp_atmospheric_trop_corrections", 0) @@ -635,7 +1114,6 @@ def build_summary_payload( float(ppp_summary_counts.get("ppp_atmospheric_ionosphere_meters", 0.0)) ) if args.summary_json is not None: - args.summary_json.parent.mkdir(parents=True, exist_ok=True) args.summary_json.write_text(json.dumps(payload, indent=2, sort_keys=True) + "\n", encoding="utf-8") return payload @@ -699,13 +1177,29 @@ def main() -> int: ensure_exists(args.nav, "navigation file") ensure_exists(args.sp3, "SP3 file") ensure_exists(args.clk, "CLK file") - selected_sources = [bool(args.ssr_rtcm), bool(args.compact_ssr), bool(args.qzss_l6)] + ensure_exists(args.antex, "ANTEX file") + ensure_exists(args.expanded_ssr, "expanded SSR CSV") + selected_sources = [bool(args.ssr_rtcm), bool(args.compact_ssr), bool(args.qzss_l6), bool(args.expanded_ssr)] if sum(1 for selected in selected_sources if selected) != 1: - raise SystemExit("Specify exactly one of --ssr-rtcm, --compact-ssr, or --qzss-l6") + raise SystemExit("Specify exactly one of --ssr-rtcm, --compact-ssr, --qzss-l6, or --expanded-ssr") + navsys_mask = effective_navsys_mask(args) + if navsys_mask < 0 or navsys_mask > 127: + raise SystemExit("--navsys must be a RTKLIB navsys mask from 0 to 127") + filter_iterations = effective_filter_iterations(args) + if filter_iterations is not None and filter_iterations < 0: + raise SystemExit("--filter-iterations must be non-negative") + if args.qzss_window_margin_seconds < 0.0: + raise SystemExit("--qzss-window-margin-seconds must be non-negative") + if args.qzss_expanded_cache is not None and args.qzss_l6 is None: + raise SystemExit("--qzss-expanded-cache requires --qzss-l6") args.out.parent.mkdir(parents=True, exist_ok=True) if args.summary_json is not None: args.summary_json.parent.mkdir(parents=True, exist_ok=True) + if args.ppp_residual_log is not None: + args.ppp_residual_log.parent.mkdir(parents=True, exist_ok=True) + if args.ppp_state_log is not None: + args.ppp_state_log.parent.mkdir(parents=True, exist_ok=True) with tempfile.TemporaryDirectory(prefix="gnss_clas_ppp_") as temp_dir: compact_summary: dict[str, object] | None = None @@ -719,35 +1213,64 @@ def main() -> int: "--out", str(args.out), ] + if navsys_mask != 0: + command.extend(["--navsys", str(navsys_mask)]) if args.qzss_l6 is not None: - compact_csv = Path(temp_dir) / "qzss_l6_expanded.csv" + compact_csv = args.qzss_expanded_cache or (Path(temp_dir) / "qzss_l6_expanded.csv") gps_week = args.qzss_gps_week if args.qzss_gps_week is not None else infer_gps_week_from_obs(args.obs) - compact_summary = expand_qzss_l6_source( - args.qzss_l6, - gps_week, - compact_csv, - compact_flush_policy=args.compact_flush_policy, - compact_atmos_merge_policy=args.compact_atmos_merge_policy, - compact_atmos_subtype_merge_policy=args.compact_atmos_subtype_merge_policy, - compact_phase_bias_merge_policy=args.compact_phase_bias_merge_policy, - compact_phase_bias_source_policy=args.compact_phase_bias_source_policy, - compact_code_bias_composition_policy=args.compact_code_bias_composition_policy, - compact_code_bias_bank_policy=args.compact_code_bias_bank_policy, - compact_phase_bias_composition_policy=args.compact_phase_bias_composition_policy, - compact_phase_bias_bank_policy=args.compact_phase_bias_bank_policy, - compact_bias_row_materialization=args.compact_bias_row_materialization, - compact_row_construction_policy=args.compact_row_construction_policy, - ) - print( - "decoded qzss l6 corrections:", - f"frames={compact_summary['frames']}", - f"subframes={compact_summary['subframes']}", - f"messages={compact_summary['messages']}", - f"rows={compact_summary['rows_written']}", - f"atmos_messages={compact_summary['atmos_messages']}", - f"atmos_rows={compact_summary['atmos_rows']}", - f"csv={compact_summary['expanded_csv']}", - ) + correction_window = observation_tow_window(args.obs, args.max_epochs) + if args.qzss_expanded_cache is not None and args.qzss_expanded_cache.exists(): + compact_summary = { + "frames": 0, + "subframes": 0, + "messages": 0, + "rows_written": None, + "rows_decoded": None, + "atmos_messages": None, + "atmos_rows": None, + "expanded_csv": str(args.qzss_expanded_cache), + "correction_window": ( + None + if correction_window is None + else { + "week": correction_window[0], + "start_tow": correction_window[1], + "end_tow": correction_window[2], + "margin_seconds": args.qzss_window_margin_seconds, + } + ), + "cache_hit": True, + } + print("reusing qzss l6 expanded cache:", f"csv={args.qzss_expanded_cache}") + else: + compact_summary = expand_qzss_l6_source( + args.qzss_l6, + gps_week, + compact_csv, + correction_window=correction_window, + correction_window_margin_seconds=args.qzss_window_margin_seconds, + compact_flush_policy=args.compact_flush_policy, + compact_atmos_merge_policy=args.compact_atmos_merge_policy, + compact_atmos_subtype_merge_policy=args.compact_atmos_subtype_merge_policy, + compact_phase_bias_merge_policy=args.compact_phase_bias_merge_policy, + compact_phase_bias_source_policy=args.compact_phase_bias_source_policy, + compact_code_bias_composition_policy=args.compact_code_bias_composition_policy, + compact_code_bias_bank_policy=args.compact_code_bias_bank_policy, + compact_phase_bias_composition_policy=args.compact_phase_bias_composition_policy, + compact_phase_bias_bank_policy=args.compact_phase_bias_bank_policy, + compact_bias_row_materialization=args.compact_bias_row_materialization, + compact_row_construction_policy=args.compact_row_construction_policy, + ) + print( + "decoded qzss l6 corrections:", + f"frames={compact_summary['frames']}", + f"subframes={compact_summary['subframes']}", + f"messages={compact_summary['messages']}", + f"rows={compact_summary['rows_written']}/{compact_summary['rows_decoded']}", + f"atmos_messages={compact_summary['atmos_messages']}", + f"atmos_rows={compact_summary['atmos_rows']}", + f"csv={compact_summary['expanded_csv']}", + ) command.extend(["--ssr", str(compact_csv)]) elif args.compact_ssr is not None: compact_csv = Path(temp_dir) / "compact_expanded.csv" @@ -759,11 +1282,18 @@ def main() -> int: f"csv={compact_summary['expanded_csv']}", ) command.extend(["--ssr", str(compact_csv)]) + elif args.expanded_ssr is not None: + command.extend(["--ssr", str(args.expanded_ssr)]) else: command.extend(["--ssr-rtcm", args.ssr_rtcm, "--ssr-step-seconds", str(args.ssr_step_seconds)]) command.append("--kinematic" if args.kinematic else "--static") - command.append("--estimate-troposphere" if args.estimate_troposphere else "--no-estimate-troposphere") + if args.estimate_troposphere is not None: + command.append( + "--estimate-troposphere" + if args.estimate_troposphere + else "--no-estimate-troposphere" + ) # When atmospheric corrections are available from CLAS, disable IFLC and # enable per-satellite ionosphere estimation with STEC constraints. command.append("--no-ionosphere-free") @@ -772,12 +1302,137 @@ def main() -> int: command.extend(["--sp3", str(args.sp3)]) if args.clk is not None: command.extend(["--clk", str(args.clk)]) + if args.antex is not None: + command.extend(["--antex", str(args.antex)]) + if args.receiver_antenna_type: + command.extend(["--receiver-antenna-type", args.receiver_antenna_type]) if args.kml is not None: command.extend(["--kml", str(args.kml)]) if args.max_epochs > 0: command.extend(["--max-epochs", str(args.max_epochs)]) if args.convergence_min_epochs is not None: command.extend(["--convergence-min-epochs", str(args.convergence_min_epochs)]) + if filter_iterations is not None: + command.extend(["--filter-iterations", str(filter_iterations)]) + if args.ppp_filter_log is not None: + command.extend(["--ppp-filter-log", str(args.ppp_filter_log)]) + if args.ppp_residual_log is not None: + command.extend(["--ppp-residual-log", str(args.ppp_residual_log)]) + if args.ppp_state_log is not None: + command.extend(["--ppp-state-log", str(args.ppp_state_log)]) + if args.ppp_correction_log is not None: + command.extend(["--ppp-correction-log", str(args.ppp_correction_log)]) + if effective_native_pntpos_parity_seed(args): + command.append("--native-pntpos-parity-seed") + elif args.native_pntpos_parity_seed is False: + command.append("--no-native-pntpos-parity-seed") + if args.clas_ppc_profile: + command.append("--clas-ppc-profile") + if args.clas_anchor_sigma is not None: + command.extend(["--clas-anchor-sigma", str(args.clas_anchor_sigma)]) + if args.clas_code_variance_scale is not None: + command.extend(["--clas-code-variance-scale", str(args.clas_code_variance_scale)]) + if args.clas_phase_variance is not None: + command.extend(["--clas-phase-variance", str(args.clas_phase_variance)]) + if args.clas_phase_continuity is not None: + command.extend(["--clas-phase-continuity", args.clas_phase_continuity]) + if args.clas_phase_bias_values is not None: + command.extend(["--clas-phase-bias-values", args.clas_phase_bias_values]) + if args.clas_phase_bias_reference_time is not None: + command.extend([ + "--clas-phase-bias-reference-time", + args.clas_phase_bias_reference_time, + ]) + if args.clas_iono_prior_variance is not None: + command.extend(["--clas-iono-prior-variance", str(args.clas_iono_prior_variance)]) + if args.clas_code_outlier_sigma_scale is not None: + command.extend([ + "--clas-code-outlier-sigma-scale", + str(args.clas_code_outlier_sigma_scale), + ]) + if args.clas_code_outlier_min_residual is not None: + command.extend([ + "--clas-code-outlier-min-residual", + str(args.clas_code_outlier_min_residual), + ]) + if args.clas_phase_outlier_sigma_scale is not None: + command.extend([ + "--clas-phase-outlier-sigma-scale", + str(args.clas_phase_outlier_sigma_scale), + ]) + if args.clas_phase_outlier_inflated_variance is not None: + command.extend([ + "--clas-phase-outlier-inflated-variance", + str(args.clas_phase_outlier_inflated_variance), + ]) + if args.clas_prior_outlier_sigma_scale is not None: + command.extend([ + "--clas-prior-outlier-sigma-scale", + str(args.clas_prior_outlier_sigma_scale), + ]) + if args.clas_prior_outlier_min_residual is not None: + command.extend([ + "--clas-prior-outlier-min-residual", + str(args.clas_prior_outlier_min_residual), + ]) + if args.clas_phase_outlier_min_residual is not None: + command.extend([ + "--clas-phase-outlier-min-residual", + str(args.clas_phase_outlier_min_residual), + ]) + if args.clas_reset_phase_ambiguity_on_outlier_inflation: + command.append("--clas-reset-phase-ambiguity-on-outlier-inflation") + if args.clas_reset_phase_ambiguity_outlier_min_rows is not None: + command.extend([ + "--clas-reset-phase-ambiguity-outlier-min-rows", + str(args.clas_reset_phase_ambiguity_outlier_min_rows), + ]) + if args.clas_reset_phase_ambiguity_outlier_min_residual is not None: + command.extend([ + "--clas-reset-phase-ambiguity-outlier-min-residual", + str(args.clas_reset_phase_ambiguity_outlier_min_residual), + ]) + if args.clas_wlnl_fixed_position_max_shift is not None: + command.extend([ + "--clas-wlnl-fixed-position-max-shift", + str(args.clas_wlnl_fixed_position_max_shift), + ]) + if args.wlnl_wl_max_fractional is not None: + command.extend([ + "--wlnl-wl-max-fractional", + str(args.wlnl_wl_max_fractional), + ]) + if args.enable_wlnl_par is True: + command.append("--enable-wlnl-par") + elif args.enable_wlnl_par is False: + command.append("--no-wlnl-par") + if args.clas_spp_clock_overwrite is True: + command.append("--clas-spp-clock-overwrite") + elif args.clas_spp_clock_overwrite is False: + command.append("--no-clas-spp-clock-overwrite") + if args.clas_kf_clock_seed_variance is not None: + command.extend([ + "--clas-kf-clock-seed-variance", + str(args.clas_kf_clock_seed_variance), + ]) + if args.ppp_process_noise_clock is not None: + command.extend([ + "--ppp-process-noise-clock", + str(args.ppp_process_noise_clock), + ]) + if args.process_noise_troposphere is not None: + command.extend([ + "--process-noise-troposphere", + str(args.process_noise_troposphere), + ]) + if args.reset_clock_to_spp is True: + command.append("--reset-clock-to-spp") + elif args.reset_clock_to_spp is False: + command.append("--no-reset-clock-to-spp") + if args.spp_seed_iono_free_code is True: + command.append("--spp-seed-iono-free-code") + elif args.spp_seed_iono_free_code is False: + command.append("--no-spp-seed-iono-free-code") if args.enable_ar: command.extend(["--enable-ar", "--ar-ratio-threshold", str(args.ar_ratio_threshold)]) @@ -787,6 +1442,11 @@ def main() -> int: print("Finished CLAS/MADOCA PPP run.") print(f" profile: {args.profile}") + print(f" navsys: {effective_navsys_mask(args)}") + print( + " estimate troposphere: " + f"{'on' if effective_estimate_troposphere(args) else 'off'}" + ) print(f" transport: {payload['ssr_transport']}") print(f" encoding: {payload['correction_encoding']}") print(f" solution: {args.out}") diff --git a/apps/gnss_ppp.cpp b/apps/gnss_ppp.cpp index 6424c79d..e0902e0c 100644 --- a/apps/gnss_ppp.cpp +++ b/apps/gnss_ppp.cpp @@ -1,22 +1,32 @@ #include +#include +#include #include +#include #include #include #include +#include #include #include +#include #include +#include #include +#include +#include +#include #include #include +#include #include namespace { struct Options { std::string obs_path; - std::string nav_path; + std::vector nav_paths; std::string sp3_path; std::string clk_path; std::string ssr_path; @@ -24,21 +34,55 @@ struct Options { std::string ionex_path; std::string dcb_path; std::string antex_path; + std::string receiver_antenna_type; std::string blq_path; std::string ocean_loading_station_name; std::string out_path; std::string summary_json_path; + std::string ppp_correction_log_path; + std::string ppp_filter_log_path; + std::string ppp_residual_log_path; + std::string ppp_state_log_path; std::string kml_path; int max_epochs = 0; + double start_tow = -1.0; // negative disables the gate int convergence_min_epochs = 20; + bool convergence_min_epochs_set = false; + int phase_measurement_min_lock_count = 1; + bool enable_initial_phase_admission_warm_start = false; + bool enable_all_frequency_initial_phase_admission_warm_start = false; + int initial_phase_admission_warm_start_navsys_mask = 0; + std::string initial_phase_admission_warm_start_satellites_arg; + std::string initial_phase_admission_warm_start_frequency_indexes_arg; + std::string initial_phase_admission_warm_start_satellite_frequency_pairs_arg; + std::string phase_admission_excluded_satellite_frequency_pairs_arg; + std::string phase_admission_excluded_before_satellite_frequency_pairs_arg; + std::string phase_admission_residual_floor_satellite_frequency_pairs_arg; + std::string exclude_known_bad_madoca_sats_preset; + bool reset_phase_ambiguity_on_before_exclusion = false; + double elevation_mask_deg = 15.0; double ssr_step_seconds = 1.0; + bool enforce_ssr_orbit_iode = false; + bool enforce_ssr_orbit_iode_admission_only = false; + int ssr_orbit_iode_admission_gate_warmup_epochs = 0; + std::string ar_method_arg; // Empty means profile default/auto. bool estimate_troposphere = true; + bool estimate_troposphere_set = false; bool estimate_ionosphere = false; bool use_ionosphere_free = true; + bool enable_per_frequency_phase_bias_states = false; + bool enable_extra_band_observations = false; + bool initialize_phase_ambiguity_with_ionosphere_state = false; + bool enable_ppp_outlier_detection = true; + bool claslib_parity = false; + bool clas_ppc_profile = false; + bool native_pntpos_parity_seed = false; + bool native_pntpos_parity_seed_set = false; bool use_clas_osr_filter = false; std::string clas_epoch_policy = "strict-osr"; std::string clas_osr_application = "full-osr"; std::string clas_phase_continuity = "full-repair"; + bool clas_phase_continuity_set = false; std::string clas_phase_bias_values = "full"; std::string clas_phase_bias_reference_time = "phase-bias-reference"; std::string clas_ssr_timing = "lag-tolerant"; @@ -47,10 +91,88 @@ struct Options { std::string clas_residual_sampling = "indexed-or-mean"; std::string clas_atmos_selection = "grid-first"; double clas_atmos_stale_after_seconds = 15.0; + double clas_anchor_sigma = -1.0; + double clas_code_variance_scale = -1.0; + double clas_phase_variance = -1.0; + double clas_iono_prior_variance = -1.0; + double clas_trop_prior_variance = -1.0; + double clas_initial_ambiguity_std_cycles = -1.0; + bool clas_initial_ambiguity_std_cycles_set = false; + double clas_outlier_sigma_scale = -1.0; + double clas_code_outlier_sigma_scale = -1.0; + double clas_phase_outlier_sigma_scale = -1.0; + double clas_prior_outlier_sigma_scale = -1.0; + double clas_phase_outlier_inflated_variance_m2 = -1.0; + double clas_code_outlier_min_residual_m = -1.0; + double clas_phase_outlier_min_residual_m = -1.0; + double clas_prior_outlier_min_residual_m = -1.0; + bool clas_reset_phase_ambiguity_on_outlier_inflation = false; + bool clas_reset_phase_ambiguity_on_outlier_inflation_set = false; + int clas_reset_phase_ambiguity_outlier_min_rows = -1; + double clas_reset_phase_ambiguity_outlier_min_residual_m = -1.0; + bool clas_kinematic_position_reseed = false; + bool clas_kinematic_position_reseed_set = false; + double clas_kinematic_position_reseed_variance = -1.0; + double clas_kinematic_position_reseed_max_residual_rms_m = -1.0; + // SPP-clock overwrite policy. Default true preserves historical CLAS + // behaviour. Set false (via --no-clas-spp-clock-overwrite) to let KF + // track the receiver clock — see clas_per_sat_residual_root_cause_2026_05_08. + bool clas_spp_clock_overwrite = true; + bool clas_spp_clock_overwrite_set = false; + double clas_kf_clock_seed_variance = -1.0; + // KF clock process noise override (m^2/s). Default <0 keeps internal + // mode-dependent fallback (0 broadcast / 100 SSR). Used to investigate + // standard PPP clock-axis residual systematic — see + // clas_per_sat_residual_root_cause_2026_05_08. + double ppp_process_noise_clock = -1.0; + // SPP-clock seed reset gate. ppp_config defaults to true. Set false to + // disable per-epoch SPP overwrite of KF clock state — used to investigate + // whether the +9m sat-uniform residual systematic comes from KF being + // re-anchored to SPP each epoch. + bool reset_clock_to_spp_each_epoch = true; + bool reset_clock_to_spp_each_epoch_set = false; + // Targeted ionosphere-free SPP seed knob. Bypasses + // --native-pntpos-parity-seed's other side effects (zero-init position, + // preserve-default-clock). Used to test whether single-frequency SPP's + // implicit absorption of ionosphere into receiver-clock state explains + // the +9m residual systematic. + bool spp_seed_iono_free_code = false; + bool spp_seed_iono_free_code_set = false; + double clas_wlnl_fixed_position_max_shift_m = -1.0; + bool enable_wlnl_par = false; + bool enable_wlnl_par_set = false; + int wlnl_par_max_exclusions = -1; + double wlnl_par_exclude_frac_threshold = -1.0; + double wlnl_wl_max_fractional_cycles = -1.0; + bool madocalib_bridge = false; + std::string madocalib_config_path; + std::vector madoca_l6e_paths; + std::vector madoca_l6d_paths; + int navsys_mask = 0; + bool navsys_mask_set = false; + int madoca_navsys_mask = 0; + std::vector madocalib_l6_paths; + std::vector madocalib_mdciono_paths; + std::string madocalib_start_time; + std::string madocalib_end_time; + double madocalib_time_interval_seconds = 0.0; + int madocalib_trace_level = 0; bool kinematic_mode = false; bool low_dynamics_mode = false; bool enable_ar = false; double ar_ratio_threshold = 3.0; + bool ar_ratio_threshold_set = false; + bool enable_ppp_holdamb = false; + double ppp_holdamb_innovation_gate_m = 0.0; // 0 = no gate + double kinematic_preconvergence_phase_residual_floor_m = 200.0; + int filter_iterations = 0; // 0 keeps PPPConfig default + double initial_ionosphere_variance = -1.0; // negative keeps default + double initial_troposphere_variance = -1.0; + double process_noise_ionosphere = -1.0; // negative keeps default + double process_noise_troposphere = -1.0; // negative keeps default + double code_phase_error_ratio_l1 = -1.0; + double code_phase_error_ratio_l2 = -1.0; + bool apply_solid_earth_tides = true; bool quiet = false; }; @@ -58,30 +180,129 @@ void printUsage(const char* program_name) { std::cout << "Usage: " << program_name << " --obs [options]\n" << "Options:\n" - << " --nav Optional broadcast navigation file\n" + << " --nav Optional broadcast navigation file; repeat to merge multiple files\n" << " --sp3 Precise orbit file\n" << " --clk Precise clock file\n" << " --ssr Simple SSR orbit/clock corrections CSV\n" << " --ssr-rtcm \n" << " RTCM SSR source converted/read for PPP use\n" + << " --madoca-l6e Native MADOCA L6E correction file; repeat for multiple channels\n" + << " --madoca-l6d Native MADOCA L6D ionosphere file; repeat for multiple channels\n" + << " Requires --madoca-l6e and an approximate receiver position in RINEX header\n" + << " --navsys \n" + << " RTKLIB navsys mask for PPP observation admission (0=all, 25=GPS+Galileo+QZSS)\n" + << " --madoca-navsys \n" + << " Native MADOCA RTKLIB navsys mask (0=all, 29 excludes BDS, 61 includes BDS)\n" << " --ionex Optional IONEX TEC map product\n" << " --dcb Optional DCB / Bias-SINEX product\n" - << " --antex Optional ANTEX file for receiver antenna PCO\n" + << " --antex Optional ANTEX file for receiver PCO/PCV and satellite PCO\n" + << " --receiver-antenna-type \n" + << " Override the RINEX receiver antenna type for ANTEX lookup\n" << " --blq Optional BLQ ocean loading coefficient file\n" << " --ocean-loading-station \n" << " Station name to select from the BLQ file\n" << " --ssr-step-seconds \n" << " Sampling step for RTCM SSR conversion (default: 1.0)\n" + << " --enforce-ssr-orbit-iode\n" + << " Require SSR orbit IODE to match a broadcast ephemeris\n" + << " --enforce-ssr-orbit-iode-admission-only\n" + << " Gate admission on SSR orbit IODE match without\n" + << " swapping broadcast-state selection for range modeling\n" + << " --ssr-orbit-iode-admission-gate-warmup-epochs \n" + << " Skip the admission-only IODE gate for the first N\n" + << " epochs so stale-CSSR-IODE observations can seed the\n" + << " filter before the gate activates (default: 0)\n" + << " --ar-method \n" + << " PPP ambiguity resolution method (default: dd-iflc).\n" + << " dd-madoca-cascaded adds MADOCALIB-style Kalman\n" + << " WL state-constraint injection before N1 LAMBDA\n" << " --out Output position file (required)\n" << " --summary-json \n" << " Optional machine-readable run summary\n" + << " --ppp-correction-log \n" + << " Optional per-epoch PPP correction diagnostics CSV\n" + << " --ppp-filter-log \n" + << " Optional per-iteration PPP filter diagnostics CSV\n" + << " --ppp-residual-log \n" + << " Optional per-row PPP residual diagnostics CSV\n" + << " --ppp-state-log \n" + << " Optional per-epoch KF state log in MADOCALIB .stat format\n" << " --kml Optional KML output\n" << " --max-epochs Limit processed epochs (default: all)\n" + << " --start-tow Drop observation epochs with TOW < this value\n" << " --convergence-min-epochs \n" << " Minimum epochs before PPP convergence/AR checks (default: 20)\n" + << " --phase-measurement-min-lock-count \n" + << " Minimum lock count before SSR/MADOCA phase rows are admitted (default: 1; 0 admits immediately)\n" + << " --enable-initial-phase-admission-warm-start\n" + << " Admit low-residual initial primary SSR/MADOCA phase rows before the lock gate\n" + << " --disable-initial-phase-admission-warm-start\n" + << " Disable initial phase admission warm start (default)\n" + << " --enable-all-frequency-initial-phase-admission-warm-start\n" + << " Admit low-residual initial SSR/MADOCA phase rows on all frequencies before the lock gate\n" + << " --disable-all-frequency-initial-phase-admission-warm-start\n" + << " Disable all-frequency initial phase admission warm start (default)\n" + << " --initial-phase-admission-warm-start-navsys \n" + << " Limit initial phase warm-start rows to a RTKLIB navsys mask (0=all; default)\n" + << " --initial-phase-admission-warm-start-sats \n" + << " Limit initial phase warm-start rows to satellite IDs (e.g. R09,C30; default all)\n" + << " --initial-phase-admission-warm-start-frequency-indexes \n" + << " Limit initial phase warm-start rows to frequency indexes (e.g. 0,1; default all)\n" + << " --initial-phase-admission-warm-start-sat-frequency-pairs \n" + << " Limit initial phase warm-start rows to satellite:frequency pairs (e.g. R09:0,C30:1; default all)\n" + << " --phase-admission-exclude-sat-frequency-pairs \n" + << " Exclude diagnostic carrier-phase rows by satellite:frequency pairs (e.g. E31:0,C29:1; default none)\n" + << " --phase-admission-exclude-sat-frequency-pairs-before \n" + << " Exclude diagnostic carrier-phase rows before GPS week/TOW (e.g. E31:0:2360:173610; default none)\n" + << " --phase-admission-residual-floor-sat-frequency-pairs \n" + << " Override phase residual floor by satellite:frequency:meters (e.g. C30:1:60; default none)\n" + << " --exclude-known-bad-madoca-sats \n" + << " Expand a named preset to --phase-admission-exclude-sat-frequency-pairs.\n" + << " Presets target satellites whose MADOCA SSR phase bias does not capture\n" + << " the per-sat hardware bias (see docs). Supported: mizu-2025-04.\n" + << " --reset-phase-ambiguity-on-before-exclusion\n" + << " Reset ambiguity state while before-time phase exclusion is active (default: keep state)\n" + << " --elevation-mask Satellite elevation mask in degrees (default: 15)\n" << " --no-estimate-troposphere\n" << " Disable zenith troposphere estimation\n" << " --estimate-troposphere Enable zenith troposphere estimation (default)\n" + << " --no-ionosphere-free\n" + << " Use per-frequency observations instead of the ionosphere-free combination\n" + << " --ionosphere-free Use the ionosphere-free dual-frequency combination (default)\n" + << " --estimate-ionosphere Estimate per-satellite ionosphere states\n" + << " --no-estimate-ionosphere\n" + << " Disable per-satellite ionosphere state estimation (default)\n" + << " --enable-extra-band-observations\n" + << " Emit RINEX observations for bands beyond the\n" + << " primary/secondary pair (enables 3rd-freq states).\n" + << " Experimental; SSR bias path for extra bands may\n" + << " not be fully wired for all product types.\n" + << " --enable-per-frequency-phase-bias-states\n" + << " Add guarded non-IFLC carrier-phase ambiguity states per frequency\n" + << " --disable-per-frequency-phase-bias-states\n" + << " Disable per-frequency phase-bias states (default)\n" + << " --enable-ionosphere-aware-phase-ambiguity-init\n" + << " Seed phase ambiguities with the current estimated ionosphere state\n" + << " --disable-ionosphere-aware-phase-ambiguity-init\n" + << " Disable ionosphere-aware phase ambiguity seeding (default)\n" + << " --disable-ppp-outlier-detection\n" + << " Disable PPP residual outlier rejection for parity experiments\n" + << " --enable-ppp-outlier-detection\n" + << " Enable PPP residual outlier rejection (default)\n" + << " --claslib-parity Use the native CLAS parity profile\n" + << " (CLAS OSR + per-frequency ionosphere estimation)\n" + << " --ported-clasnat Alias for --claslib-parity\n" + << " --clas-ppc-profile Use the PPC kinematic CLAS profile\n" + << " (CLAS OSR, residual ionosphere, 1 m anchor,\n" + << " 0.25 m^2 iono prior,\n" + << " 2x code variance, CLASLIB phase variance, kinematic reseed)\n" + << " --native-pntpos-parity-seed\n" + << " Use native SPP settings that mirror pntpos behavior for PPP seed\n" + << " --no-native-pntpos-parity-seed\n" + << " Use the default PPP SPP seed settings\n" + << " --claslib-pntpos-seed Deprecated alias for --native-pntpos-parity-seed\n" + << " --no-claslib-pntpos-seed\n" + << " Deprecated alias for --no-native-pntpos-parity-seed\n" << " --clas-osr Use the CLAS OSR-based PPP-RTK filter path\n" << " --clas-epoch-policy \n" << " CLAS epoch boundary policy (default: strict-osr)\n" @@ -105,14 +326,105 @@ void printUsage(const char* program_name) { << " CLAS atmosphere token selection policy (default: grid-first)\n" << " --clas-atmos-stale-after-seconds \n" << " Balanced policy stale threshold (default: 15.0)\n" + << " --clas-anchor-sigma \n" + << " CLAS OSR SPP anchor sigma (default: 5.0)\n" + << " --clas-code-variance-scale \n" + << " CLAS OSR code variance multiplier (default: 8.0)\n" + << " --clas-phase-variance \n" + << " CLAS OSR carrier-phase base variance (default: 0.01)\n" + << " --clas-iono-prior-variance \n" + << " CLAS OSR residual-ionosphere prior variance (default: 0.25)\n" + << " --clas-trop-prior-variance \n" + << " CLAS OSR troposphere prior variance (default: 0.0001)\n" + << " --clas-initial-ambiguity-std-cycles \n" + << " CLAS ambiguity init std per frequency; 0 uses the PPP default variance\n" + << " --clas-outlier-sigma-scale \n" + << " CLAS measurement outlier sigma gate; 0 disables inflation\n" + << " --clas-code-outlier-sigma-scale \n" + << " CLAS code-only outlier sigma gate; 0 disables, default follows --clas-outlier-sigma-scale\n" + << " --clas-code-outlier-min-residual \n" + << " CLAS code-only outlier gate absolute residual floor; 0 disables\n" + << " --clas-phase-outlier-sigma-scale \n" + << " CLAS phase-only outlier sigma gate; 0 disables, default follows --clas-outlier-sigma-scale\n" + << " --clas-phase-outlier-inflated-variance \n" + << " CLAS phase-only inflated measurement variance; default uses 1e10\n" + << " --clas-prior-outlier-sigma-scale \n" + << " CLAS prior-constraint outlier sigma gate; 0 disables, default follows --clas-outlier-sigma-scale\n" + << " --clas-prior-outlier-min-residual \n" + << " CLAS prior-constraint outlier gate absolute residual floor; 0 disables\n" + << " --clas-phase-outlier-min-residual \n" + << " CLAS phase-only outlier gate absolute residual floor; 0 disables\n" + << " --clas-reset-phase-ambiguity-on-outlier-inflation\n" + << " Reset CLAS phase ambiguity states touched by inflated phase rows\n" + << " --no-clas-reset-phase-ambiguity-on-outlier-inflation\n" + << " Disable reset of inflated CLAS phase ambiguity states (default)\n" + << " --clas-reset-phase-ambiguity-outlier-min-rows \n" + << " Minimum inflated phase rows in an epoch before reset is applied\n" + << " --clas-reset-phase-ambiguity-outlier-min-residual \n" + << " Minimum inflated phase residual for ambiguity reset membership\n" + << " --clas-kinematic-reseed-position\n" + << " Re-init position state from SPP every kinematic epoch (default off)\n" + << " --no-clas-kinematic-reseed-position\n" + << " Disable kinematic position re-seed\n" + << " --clas-kinematic-reseed-position-variance \n" + << " Variance reset on kinematic re-seed (default: 10000.0)\n" + << " --clas-kinematic-reseed-position-max-residual-rms \n" + << " Skip reseed when SPP residual_rms exceeds this (m); guards against urban-canyon SPP pin-to-bad. <=0 disables (default).\n" + << " --clas-spp-clock-overwrite\n" + << " Overwrite KF receiver-clock state with SPP every kinematic epoch (default on; historical behaviour).\n" + << " --no-clas-spp-clock-overwrite\n" + << " Skip the SPP clock overwrite and let the KF track receiver clock. Targets the +11m sat-uniform residual systematic documented in clas_per_sat_residual_root_cause_2026_05_08.\n" + << " --clas-kf-clock-seed-variance \n" + << " Initial KF clock variance when --no-clas-spp-clock-overwrite is set (default: 10000.0). Ignored when SPP overwrite is on.\n" + << " --ppp-process-noise-clock \n" + << " Override KF receiver-clock process noise (default: 0 broadcast / 100 SSR). Set >0 to use a single value for both modes — used to investigate the +11m clock-axis residual systematic in standard PPP path.\n" + << " --reset-clock-to-spp Force KF receiver-clock state to SPP every epoch (default on; historical behaviour for broadcast-rtklib path).\n" + << " --no-reset-clock-to-spp Skip SPP overwrite of KF receiver-clock state — used to investigate whether the +9m sat-uniform residual systematic comes from per-epoch SPP re-anchoring.\n" + << " --spp-seed-iono-free-code\n" + << " Use ionosphere-free L1+L2 code combination in the SPP seed. Targeted alternative to --native-pntpos-parity-seed without its zero-init / preserve-default-clock side effects. Used to test whether single-frequency SPP absorbing iono into receiver-clock state explains the +9m residual systematic.\n" + << " --no-spp-seed-iono-free-code\n" + << " Force single-frequency L1 code in the SPP seed (default).\n" + << " --enable-wlnl-par Enable WLNL Partial AR (greedy worst-|frac| exclusion when ratio fails)\n" + << " --no-wlnl-par Disable WLNL Partial AR\n" + << " --wlnl-par-max-exclusions \n" + << " Cap PAR exclusion iterations (default: 4)\n" + << " --wlnl-par-exclude-frac-threshold \n" + << " Skip excluding pairs whose |frac| is below this (default: 0.20)\n" + << " --wlnl-wl-max-fractional \n" + << " Max MW wide-lane fractional cycle before rejecting WL fix (default: 0.25)\n" + << " --madocalib-bridge Delegate this run to linked MADOCALIB postpos()\n" + << " (requires CMake -DMADOCALIB_PARITY_LINK=ON)\n" + << " --madocalib-l6 Extra MADOCA L6 input file; repeat for two-channel L6E\n" + << " --madocalib-mdciono \n" + << " Extra MADOCA L6D ionosphere file; repeat up to three\n" + << " --madocalib-config \n" + << " MADOCALIB rnx2rtkp config (default: sample.conf)\n" + << " --madocalib-start