From 32294044520e58fb8706b07f24825760396b2a59 Mon Sep 17 00:00:00 2001 From: Colin Orion Chandler Date: Fri, 6 Nov 2026 12:39:27 -0800 Subject: [PATCH 1/3] Use MPC Halley-Marcus photometry by default for comets --- .github/workflows/comet-photometry.yml | 29 ++ README.md | 53 +++ pyproject.toml | 1 + src/ponder/__main__.py | 14 + src/ponder/comet_photometry.py | 94 ++++++ src/ponder/runner.py | 182 +++++++++-- src/ponder/utils.py | 81 +++-- tests/data/comet_external.ini | 28 ++ tests/ponder/test_comet_photometry.py | 429 +++++++++++++++++++++++++ 9 files changed, 859 insertions(+), 52 deletions(-) create mode 100644 .github/workflows/comet-photometry.yml create mode 100644 src/ponder/comet_photometry.py create mode 100644 tests/data/comet_external.ini create mode 100644 tests/ponder/test_comet_photometry.py diff --git a/.github/workflows/comet-photometry.yml b/.github/workflows/comet-photometry.yml new file mode 100644 index 0000000..935eb9a --- /dev/null +++ b/.github/workflows/comet-photometry.yml @@ -0,0 +1,29 @@ +name: Comet photometry integration + +on: + pull_request: + branches: [main] + workflow_dispatch: + +permissions: + contents: read + +jobs: + comet-photometry: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v6 + - uses: actions/setup-python@v6 + with: + python-version: '3.13' + - uses: astral-sh/setup-uv@v7 + - name: Install companion add-on and Ponder + # Temporary immutable dependency for the draft PR. Replace with the + # released sorcha-addons minimum version before merging the default. + run: | + uv pip install --system 'sorcha-addons @ git+https://github.com/ColinOrionChandler/sorcha-addons.git@7cfc0e9f173cc194379991f75d468c4cc482cd13' + uv pip install --system -e '.[dev]' + - name: Verify add-on availability and run tests + run: | + python -c 'from sorcha_addons.activity.mpc_comet import MPCCometActivity' + python -m pytest tests -q diff --git a/README.md b/README.md index 750f609..34a24b6 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,59 @@ # ponder +### Default MPC comet photometry + +Comet runs (`ponder --comet` or `run_ponder(..., comet=True)`) use the MPC comet +brightness exponent and the Schleicher–Marcus (Halley–Marcus) dust phase function +by default. Asteroid runs retain their existing photometry. Select +`--comet-photometry legacy` explicitly to reproduce the old distance-only law. + +Ponder requires sorcha-addons. Until the new model is released, install **the +sorcha-addons checkout containing `activity/mpc_comet`** into the same Python +environment as Ponder; an older add-ons release is insufficient. Ponder checks +for the model before starting a comet run and reports an error if it is missing. +For a local checkout: + +```bash +python -m pip install -e /path/to/sorcha-addons +ponder --comet \ + --config /path/to/sorcha.ini --orbits /path/to/comets.json \ + --db /path/to/pointings.db --work-dir /path/to/work --results-dir /path/to/results +``` + +The add-on computes +`m = H_filter + 5 log10(delta) + 2.5 G log10(r_h) - 2.5 log10(Phi_HM(alpha))` +inside Sorcha, before observational noise and selection. Ponder treats MPC `H` as `H_r`, preserves `g-r = 0.5`, and supplies +fixed colors for the other Rubin bands (`u-r = 1.8`, `i-r = -0.2`, +`z-r = -0.3`, `y-r = -0.4`). The catalog `G` is the comet +heliocentric exponent, **not** asteroid `GS`. Missing/nonfinite `G` becomes 4; +`comet_G_assumed` flags this in each `*_physparams_complex.csv` alongside the +actual `comet_G`. Missing H stays unknown. Those columns are also available +in Sorcha's `output_columns = all` output. + +Ponder derives an effective configuration with `comet_activity = mpc_hm`, +`phase_function = none`, and `brute_force = True` (no early faint-object cull). +Conflicting phase, activity, or rotational-lightcurve settings are rejected. +Comet runs bypass the saturation cut. Subsequent detection +filters remain as configured; disable magnitude/SNR, fading, and linking cuts +and observational effects to retain unknown magnitudes in geometry-only output. + +New-model work and results are isolated below each supplied root in +`photometry/mpc-hm-/`. The signature includes the add-on implementation, +model version, input configuration, color assumptions, and numerical dependency versions. The result directory contains +`photometry.json` with those details and the effective-config path. State files live in the model-specific work directory and +are separately keyed, so the first default comet run processes full history without reusing legacy +completion state. Changes to H/G also trigger historical +reruns. Use the model-specific results directory when consuming or publishing results. + +This first-order total-magnitude model assumes dust-dominated light and treats H +as a zero-phase normalization. Catalog H/G may already absorb phase or aperture +effects. It does not recalibrate bands, model outbursts/gas emission, or convert +integrated coma brightness to point-source photometry. See the add-on's +`docs/community_modules/activity/mpc_comet.rst` for references and examples. + + + [![Template](https://img.shields.io/badge/Template-LINCC%20Frameworks%20Python%20Project%20Template-brightgreen)](https://lincc-ppt.readthedocs.io/en/latest/) diff --git a/pyproject.toml b/pyproject.toml index 0032807..fd37f67 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -20,6 +20,7 @@ dependencies = [ "astropy", "pandas < 3.0", "sorcha", + "sorcha-addons", "tqdm", "pyarrow", ] diff --git a/src/ponder/__main__.py b/src/ponder/__main__.py index a3c38e7..2da8439 100644 --- a/src/ponder/__main__.py +++ b/src/ponder/__main__.py @@ -2,6 +2,7 @@ from datetime import datetime, timezone from pathlib import Path +from .comet_photometry import COMET_PHOTOMETRY_CHOICES, resolve_model from .runner import DEFAULT_DEBUG_FAILED_CHUNK_SIZE, DEFAULT_ISOLATE_FAILING_ROWS, run_ponder from .utils import get_current_orbits @@ -12,6 +13,15 @@ def main(): parser.add_argument("--orbits", type=Path, help="Path to the orbits JSON file") parser.add_argument("--config", type=Path, help="Path to the configuration file") parser.add_argument("--comet", action="store_true", help="Whether to run comet analysis") + parser.add_argument( + "--comet-photometry", + choices=COMET_PHOTOMETRY_CHOICES, + default=None, + help="Comet magnitude model (default for --comet: mpc-hm; legacy restores the distance-only law).", + ) + parser.add_argument( + "--results-dir", type=Path, default=Path.cwd() / "results", help="Directory to store run results" + ) parser.add_argument( "--work-dir", type=Path, default=Path.cwd() / "work", help="Directory to store intermediate files" ) @@ -87,6 +97,7 @@ def main(): ) args = parser.parse_args() + args.comet_photometry = resolve_model(args.comet_photometry, comet=args.comet) orbit_path = args.orbits if args.download_orbits: date_str = datetime.now(timezone.utc).strftime("%Y-%m-%d") @@ -109,6 +120,9 @@ def main(): orbit_path, args.config, comet=args.comet, + comet_photometry=args.comet_photometry, + work_dir=args.work_dir, + results_dir=args.results_dir, filter_orbits=not args.no_filter_orbits, chunk_size=args.chunk_size, sorcha_workers=args.sorcha_workers, diff --git a/src/ponder/comet_photometry.py b/src/ponder/comet_photometry.py new file mode 100644 index 0000000..4188246 --- /dev/null +++ b/src/ponder/comet_photometry.py @@ -0,0 +1,94 @@ +"""Catalog and configuration plumbing for Sorcha comet photometry. + +The magnitude calculation lives in sorcha-addons, not in Ponder. +""" + +import configparser +import hashlib +import inspect +import json +from importlib.metadata import version +from pathlib import Path + +import numpy as np +import pandas as pd + +COMET_PHOTOMETRY_CHOICES = ("legacy", "mpc-hm") +COMET_COMPLEX_COLUMNS = ("comet_G", "comet_G_assumed") + + +def validate_model(model, comet=True): + if model not in COMET_PHOTOMETRY_CHOICES: + raise ValueError(f"Unknown comet photometry model: {model}") + if model != "legacy" and not comet: + raise ValueError("--comet-photometry mpc-hm requires --comet") + + +def resolve_model(model, comet): + """Default comet runs to MPC/Halley–Marcus, leaving asteroid runs unchanged.""" + if model is None: + model = "mpc-hm" if comet else "legacy" + validate_model(model, comet) + return model + + +def comet_slopes(catalog): + """Return catalog G and an explicit flag for the G=4 fallback.""" + slope = pd.to_numeric(catalog.reindex(columns=["G"])["G"], errors="coerce") + assumed = ~np.isfinite(slope) + return slope.mask(assumed, 4.0), assumed + + +def model_metadata(): + """Fail early if the selected add-on is absent; fingerprint its implementation.""" + try: + from sorcha_addons.activity.mpc_comet import mpc_comet_activity as model + except ImportError as exc: + raise RuntimeError( + "mpc-hm requires sorcha-addons with the mpc_comet activity model. " + "Install the matching sorcha-addons checkout in this environment, or explicitly " + "select --comet-photometry legacy to use the old distance-only law." + ) from exc + from sorcha_addons import __version__ as addons_version + + from .utils import SORCHA_DEFAULT_COLOR_OFFSETS + + return { + "ponder_photometry_version": 1, + "sorcha_addons_version": addons_version, + "color_offsets": dict(SORCHA_DEFAULT_COLOR_OFFSETS), + "model": "mpc-hm", + "model_version": model.MODEL_VERSION, + "model_sha256": hashlib.sha256(inspect.getsource(model).encode()).hexdigest(), + "sorcha_version": version("sorcha"), + "sbpy_version": version("sbpy"), + "scipy_version": version("scipy"), + "missing_G": 4.0, + "H_band_assumption": "MPC H treated as H_r; existing fixed colors", + } + + +def metadata_signature(metadata): + return hashlib.sha256(json.dumps(metadata, sort_keys=True).encode()).hexdigest()[:12] + + +def configured_model(path): + parser = configparser.ConfigParser(interpolation=None) + if not parser.read(path): + raise FileNotFoundError(path) + activity = parser.get("ACTIVITY", "comet_activity", fallback="none").strip().lower() + return "mpc-hm" if activity == "mpc_hm" else "legacy" + + +def complex_parameters_path(physical_path): + path = Path(physical_path) + return path.with_name(path.stem + "_complex.csv") + + +def write_physical_inputs(phys, path): + """Split paired physical/complex inputs after each chunk has been sliced.""" + if "comet_G" in phys: + columns = ["ObjID", *COMET_COMPLEX_COLUMNS] + phys[columns].to_csv(complex_parameters_path(path), index=False) + phys = phys.drop(columns=list(COMET_COMPLEX_COLUMNS)) + phys.to_csv(path, index=False) diff --git a/src/ponder/runner.py b/src/ponder/runner.py index 528647a..cd98d29 100644 --- a/src/ponder/runner.py +++ b/src/ponder/runner.py @@ -1,18 +1,30 @@ +import configparser import gzip import hashlib +import io import json import os import signal import subprocess +import sys from concurrent.futures import ThreadPoolExecutor, as_completed from dataclasses import dataclass from datetime import datetime, timezone +from functools import wraps from pathlib import Path from time import perf_counter import pandas as pd from tqdm.auto import tqdm +from .comet_photometry import ( + complex_parameters_path, + configured_model, + metadata_signature, + model_metadata, + resolve_model, + write_physical_inputs, +) from .utils import * # TODO: should be some sort of configurable @@ -26,6 +38,7 @@ DEFAULT_DEBUG_FAILED_CHUNK_SIZE = 250 DEFAULT_ISOLATE_FAILING_ROWS = True CATALOG_ROW_COLUMN = "ponder_catalog_row" +COMET_CONFIG_SUFFIX = "_comet_no_bright_limit" # Sorcha output schemas have shifted between releases, so audits look for the # first usable object/time columns instead of assuming one fixed CSV shape. @@ -51,8 +64,8 @@ ] -def object_hashes(objects): - return {obj_id(obj): hash_orbit(obj) for obj in objects} +def object_hashes(objects, include_photometry=False): + return {obj_id(obj): hash_orbit(obj, include_photometry=include_photometry) for obj in objects} def _state_path_token(path): @@ -62,11 +75,14 @@ def _state_path_token(path): return safe[:80] or "pointings" -def state_files_for_run(db_path, comet): +def state_files_for_run(db_path, comet, photometry_signature=None): """Scope incremental state to the pointing DB and object mode.""" mode = "comets" if comet else "asteroids" db_path = Path(db_path).expanduser().resolve() - digest = hashlib.sha256(f"{db_path}|{mode}".encode()).hexdigest()[:12] + identity = f"{db_path}|{mode}" + if photometry_signature: + identity += f"|photometry={photometry_signature}" + digest = hashlib.sha256(identity.encode()).hexdigest()[:12] stem = f"{_state_path_token(db_path)}_{mode}_{digest}" return Path(f"state_{stem}.json"), Path(f"element_hashes_{stem}.json") @@ -95,6 +111,66 @@ def run_context_digest(db_path, config_path, comet, **metadata): return hashlib.sha256("\0".join(parts).encode()).hexdigest() +def prepare_sorcha_config_for_mode(config_path, comet, work_dir, comet_photometry=None): + """Configure comet photometry and disable saturation; preserve asteroid configs.""" + comet_photometry = resolve_model(comet_photometry, comet) + config_path = Path(config_path) + if not comet: + return config_path + + lines = config_path.read_text().splitlines(keepends=True) + output_lines = [] + in_saturation = False + removed = 0 + for line in lines: + stripped = line.strip() + if stripped.startswith("[") and stripped.endswith("]"): + in_saturation = stripped[1:-1].strip().casefold() == "saturation" + if ( + in_saturation + and stripped + and not stripped.startswith(("#", ";")) + and "=" in line + and line.split("=", 1)[0].strip().casefold() == "bright_limit" + ): + output_lines.append("# bright_limit omitted by Ponder for comet mode (no saturation cut)\n") + removed += 1 + continue + output_lines.append(line) + + suffix = config_path.suffix or ".ini" + effective_path = Path(work_dir) / (f"{config_path.stem}{COMET_CONFIG_SUFFIX}{suffix}") + content = "".join(output_lines) + if comet_photometry == "mpc-hm": + parser = configparser.ConfigParser(interpolation=None) + parser.read_string(content) + for section, key, allowed in ( + ("PHASECURVES", "phase_function", {"none"}), + ("LIGHTCURVE", "lc_model", {"none"}), + ("ACTIVITY", "comet_activity", {"none", "mpc_hm"}), + ): + value = parser.get(section, key, fallback="none").strip().lower() + if value not in allowed: + raise ValueError(f"mpc-hm conflicts with [{section}] {key}={value}") + for section in ("PHASECURVES", "ACTIVITY", "EXPERT"): + if not parser.has_section(section): + parser.add_section(section) + parser.set("PHASECURVES", "phase_function", "none") + parser.set("ACTIVITY", "comet_activity", "mpc_hm") + parser.set("EXPERT", "brute_force", "True") + stream = io.StringIO() + parser.write(stream) + content = stream.getvalue() + effective_path = effective_path.with_name(effective_path.stem + "_mpc_hm" + suffix) + if not effective_path.exists() or effective_path.read_text() != content: + effective_path.write_text(content) + print( + f" Comet Sorcha config — photometry: {comet_photometry}; saturation disabled; " + f"removed bright_limit entries: {removed}; effective config: {effective_path}" + ) + return effective_path + + @dataclass(frozen=True) class SorchaChunk: """A row range plus the concrete input, output, and marker files it owns.""" @@ -160,8 +236,9 @@ def _terminate_process_group(proc): def run_sorcha(orbits, physparams, output, db, config, timeout=None): print(output.with_suffix("")) command = [ - "sorcha", - "run", + sys.executable, + "-m", + "sorcha_cmdline.run", "-c", str(config), "--ob", @@ -178,6 +255,11 @@ def run_sorcha(orbits, physparams, output, db, config, timeout=None): str(output.stem), "-f", ] + if configured_model(config) == "mpc-hm": + complex_path = complex_parameters_path(physparams) + if not complex_path.is_file(): + raise FileNotFoundError(f"Missing MPC comet parameters: {complex_path}") + command.extend(["--cp", str(complex_path)]) proc = subprocess.Popen(command, start_new_session=True) try: returncode = proc.wait(timeout=timeout) @@ -189,13 +271,14 @@ def run_sorcha(orbits, physparams, output, db, config, timeout=None): raise subprocess.CalledProcessError(returncode, command) -def build_id_set_inputs(objects, ids, job_name, comet): +def build_id_set_inputs(objects, ids, job_name, comet, comet_photometry=None): + comet_photometry = resolve_model(comet_photometry, comet) if len(ids) == 0: print(f" [skipping {job_name} — no objects to process]") return if comet: - orbs, phys = comets_to_sorcha_inputs(objects, ids) + orbs, phys = comets_to_sorcha_inputs(objects, ids, comet_photometry=comet_photometry) else: orbs, phys = mpcorb_to_sorcha_inputs(objects, ids) @@ -212,7 +295,11 @@ def build_id_set_inputs(objects, ids, job_name, comet): f"{len(catalog_rows)} != {len(orbs)}" ) - return orbs.reset_index(drop=True), phys.reset_index(drop=True), catalog_rows.reset_index(drop=True) + return ( + orbs.reset_index(drop=True), + phys.reset_index(drop=True), + catalog_rows.reset_index(drop=True), + ) def catalog_rows_for_sorcha_inputs(objects, ids, comet): @@ -226,18 +313,7 @@ def catalog_rows_for_sorcha_inputs(objects, ids, comet): ) df = df[id_col.isin(ids)].copy() - required_epoch_columns = [ - "Epoch_year", - "Epoch_month", - "Epoch_day", - "Year_of_perihelion", - "Month_of_perihelion", - "Day_of_perihelion", - ] - bad = [] - for row_index, row in df.iterrows(): - if any(pd.isna(row.get(column)) for column in required_epoch_columns): - bad.append(row_index) + bad = [index for index, row in df.iterrows() if not valid_comet_sorcha_row(row)] return df.drop(bad).reset_index(drop=True) return df[df["Principal_desig"].isin(ids)].copy().reset_index(drop=True) @@ -445,7 +521,7 @@ def write_chunk_inputs(chunk, orbs, phys): chunk.orbits_path.parent.mkdir(parents=True, exist_ok=True) chunk.physparams_path.parent.mkdir(parents=True, exist_ok=True) orbs.iloc[chunk.row_start : chunk.row_end].to_csv(chunk.orbits_path, index=False) - phys.iloc[chunk.row_start : chunk.row_end].to_csv(chunk.physparams_path, index=False) + write_physical_inputs(phys.iloc[chunk.row_start : chunk.row_end], chunk.physparams_path) return perf_counter() - started @@ -1569,9 +1645,9 @@ def run_sorcha_chunks( physpar_path = WORK_DIR / f"job_{job_name}_physparams.csv" out_path = RESULTS_DIR / f"{time[:10]}_job_{job_name}.csv" orbs.to_csv(orbits_path, index=False) - phys.to_csv(physpar_path, index=False) + write_physical_inputs(phys, physpar_path) run_sorcha(orbits_path, physpar_path, out_path, db, config, timeout=timeout) - return + return True digest = dataframe_digest(orbs, phys) if context_digest: @@ -1799,7 +1875,8 @@ def run_id_set( force_debug_chunking=False, isolate_failing_rows=DEFAULT_ISOLATE_FAILING_ROWS, ): - inputs = build_id_set_inputs(objects, ids, job_name, comet) + photometry = configured_model(config) if comet else "legacy" + inputs = build_id_set_inputs(objects, ids, job_name, comet, comet_photometry=photometry) if not inputs: return True @@ -1825,6 +1902,20 @@ def run_id_set( ) +def _restore_run_directories(function): + @wraps(function) + def wrapped(*args, **kwargs): + global WORK_DIR, RESULTS_DIR + previous = WORK_DIR, RESULTS_DIR + try: + return function(*args, **kwargs) + finally: + WORK_DIR, RESULTS_DIR = previous + + return wrapped + + +@_restore_run_directories def run_ponder( db_path, object_path, @@ -1841,8 +1932,26 @@ def run_ponder( debug_failed_chunk_size=DEFAULT_DEBUG_FAILED_CHUNK_SIZE, force_debug_chunking=False, isolate_failing_rows=DEFAULT_ISOLATE_FAILING_ROWS, + work_dir=None, + results_dir=None, + comet_photometry=None, ): - """Run Ponder on the given configs.""" + """Run Ponder, using MPC/Halley–Marcus photometry by default for comets.""" + global WORK_DIR, RESULTS_DIR + comet_photometry = resolve_model(comet_photometry, comet) + photometry = model_metadata() if comet_photometry == "mpc-hm" else None + if photometry: + photometry["input_config_sha256"] = hashlib.sha256(Path(config_path).read_bytes()).hexdigest() + signature = metadata_signature(photometry) if photometry else None + if work_dir is not None: + WORK_DIR = Path(work_dir) + if results_dir is not None: + RESULTS_DIR = Path(results_dir) + if signature: + WORK_DIR = WORK_DIR / "photometry" / f"mpc-hm-{signature}" + RESULTS_DIR = RESULTS_DIR / "photometry" / f"mpc-hm-{signature}" + if configured_model(config_path) == "mpc-hm" and comet_photometry != "mpc-hm": + raise ValueError("Config selects mpc_hm; use --comet with mpc-hm photometry") db_path = Path(db_path) object_path = Path(object_path) config_path = Path(config_path) @@ -1855,8 +1964,12 @@ def run_ponder( if isolate_failing_rows and debug_failed_chunk_size <= 0: raise ValueError("--isolate-failing-rows requires --debug-failed-chunk-size > 0") - WORK_DIR.mkdir(exist_ok=True) - RESULTS_DIR.mkdir(exist_ok=True) + WORK_DIR.mkdir(parents=True, exist_ok=True) + RESULTS_DIR.mkdir(parents=True, exist_ok=True) + config_path = prepare_sorcha_config_for_mode(config_path, comet, WORK_DIR, comet_photometry) + if photometry: + photometry["effective_config"] = str(config_path.resolve()) + (RESULTS_DIR / "photometry.json").write_text(json.dumps(photometry, indent=2) + "\n") ts = datetime.now(timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ") catalog_snapshot_path = write_mpc_catalog_snapshot(object_path, RESULTS_DIR, ts) @@ -1864,7 +1977,10 @@ def run_ponder( # Persisted state is scoped by DB path and object mode so runs against a # different pointing database do not reuse old "last_mjd" or orbit hashes. - state_file, hashes_file = state_files_for_run(db_path, comet) + state_file, hashes_file = state_files_for_run(db_path, comet, photometry_signature=signature) + if photometry: + state_file = WORK_DIR / state_file + hashes_file = WORK_DIR / hashes_file state = json.loads(state_file.read_text()) if state_file.exists() else {"last_mjd": 0.0} prev_hashes = json.loads(hashes_file.read_text()) if hashes_file.exists() else {} @@ -1882,7 +1998,9 @@ def run_ponder( print(f"\n[{ts}] Starting cycle") # -- load and diff inputs -- - new_ids, updated_ids, unchanged_ids = diff_objects(objects, prev_hashes) + new_ids, updated_ids, unchanged_ids = diff_objects( + objects, prev_hashes, include_photometry=photometry is not None + ) print( f" Objects — new: {len(new_ids)} updated: {len(updated_ids)} " f"unchanged: {len(unchanged_ids)}" ) @@ -1896,6 +2014,7 @@ def run_ponder( db_path, config_path, comet, + **({"photometry_signature": signature} if signature else {}), pointing_scope="full", db_max_mjd=db_last_mjd, db_row_count=total_pts, @@ -1906,6 +2025,7 @@ def run_ponder( db_path, config_path, comet, + **({"photometry_signature": signature} if signature else {}), pointing_scope="new", previous_last_mjd=state["last_mjd"], db_max_mjd=db_last_mjd, @@ -1984,7 +2104,7 @@ def run_ponder( # -- persist state for next run -- state["last_mjd"] = db_last_mjd state_file.write_text(json.dumps(state)) - hashes_file.write_text(json.dumps(object_hashes(objects))) + hashes_file.write_text(json.dumps(object_hashes(objects, include_photometry=photometry is not None))) def main(): diff --git a/src/ponder/utils.py b/src/ponder/utils.py index a404604..be679fb 100644 --- a/src/ponder/utils.py +++ b/src/ponder/utils.py @@ -1,5 +1,6 @@ import gzip import hashlib +import math import re import sqlite3 @@ -8,6 +9,8 @@ import requests from astropy.time import Time +from .comet_photometry import comet_slopes, validate_model + # orbital element fields - a change in any of these triggers a full-history rerun ORBIT_FIELDS = [ "e", @@ -30,6 +33,24 @@ MPCORB_OBSERVATION_ARC_LIMIT_DAYS = 3.0 +SORCHA_DEFAULT_GS = 0.15 +SORCHA_DEFAULT_COLOR_OFFSETS = { + "u-r": 1.8, + "g-r": 0.5, + "i-r": -0.2, + "z-r": -0.3, + "y-r": -0.4, +} +COMET_REQUIRED_INPUT_COLUMNS = ( + "Epoch_year", + "Epoch_month", + "Epoch_day", + "Year_of_perihelion", + "Month_of_perihelion", + "Day_of_perihelion", +) + + # -- comet diff helpers -- @@ -41,12 +62,13 @@ def obj_id(comet): ) -def hash_orbit(comet): - key = "|".join(str(comet.get(f, "")) for f in ORBIT_FIELDS) +def hash_orbit(comet, include_photometry=False): + fields = ORBIT_FIELDS + (["H", "G"] if include_photometry else []) + key = "|".join(str(comet.get(f, "")) for f in fields) return hashlib.md5(key.encode()).hexdigest() -def diff_objects(objects, prev_hashes): +def diff_objects(objects, prev_hashes, include_photometry=False): """ Given a list of objects and a dict of previous comet hashes, return lists of new, updated, and unchanged comets. @@ -58,7 +80,7 @@ def diff_objects(objects, prev_hashes): new, updated, unchanged = [], [], [] for c in objects: oid = obj_id(c) - h = hash_orbit(c) + h = hash_orbit(c, include_photometry=include_photometry) if oid not in prev_hashes: new.append(oid) elif prev_hashes[oid] != h: @@ -205,7 +227,21 @@ def mpcorb_to_sorcha_inputs(mpcorb_json, ids): # -- comet -> sorcha format converters -- -def comets_to_sorcha_inputs(comets_json, ids): +def valid_comet_sorcha_row(obj): + """Return whether an MPC comet row has every value required by Sorcha.""" + return not any(_is_missing(obj.get(column)) for column in COMET_REQUIRED_INPUT_COLUMNS) + + +def add_default_sorcha_physical_columns(phys): + """Add generic colour columns needed for all Rubin filters in Sorcha.""" + phys["GS"] = SORCHA_DEFAULT_GS + for column, value in SORCHA_DEFAULT_COLOR_OFFSETS.items(): + phys[column] = value + return phys + + +def comets_to_sorcha_inputs(comets_json, ids, comet_photometry="mpc-hm"): + validate_model(comet_photometry) # TODO: potentially refactor so that it only does the conversion once df = pd.DataFrame(comets_json) @@ -215,20 +251,11 @@ def comets_to_sorcha_inputs(comets_json, ids): ) df = df[id_col.isin(ids)].copy() - # build epoch and perihelion MJD strings, dropping rows with missing epochs + # Build epoch and perihelion MJD strings, dropping rows that cannot form a + # complete Sorcha orbit/physical-parameter pair. epoch_strings, tp_strings, bad = [], [], [] for r, row in df.iterrows(): - if any( - pd.isna(row.get(c)) - for c in [ - "Epoch_year", - "Epoch_month", - "Epoch_day", - "Year_of_perihelion", - "Month_of_perihelion", - "Day_of_perihelion", - ] - ): + if not valid_comet_sorcha_row(row): # print(f" [warn] dropping bad epoch: {row.get('Designation_and_name')}") bad.append(r) continue @@ -243,8 +270,12 @@ def comets_to_sorcha_inputs(comets_json, ids): df = df.drop(bad).reset_index(drop=True) - df["epochMJD_TDB"] = Time(epoch_strings).mjd - df["t_p_MJD_TDB"] = Time([s for s, _ in tp_strings], format="iso").mjd + [f for _, f in tp_strings] + if df.empty: + df["epochMJD_TDB"] = pd.Series(dtype=float) + df["t_p_MJD_TDB"] = pd.Series(dtype=float) + else: + df["epochMJD_TDB"] = Time(epoch_strings).mjd + df["t_p_MJD_TDB"] = Time([s for s, _ in tp_strings], format="iso").mjd + [f for _, f in tp_strings] df["FORMAT"] = "COM" # -- orbits file -- @@ -275,10 +306,18 @@ def comets_to_sorcha_inputs(comets_json, ids): phys = pd.DataFrame( { "ObjID": orbs["ObjID"], - "H_r": df["H"].values, - "g-r": 0.5, + # Brightness is not an orbit requirement. Preserve unknown H as + # NaN; geometry-only Sorcha configurations can retain these rows. + "H_r": pd.to_numeric(df.reindex(columns=["H"])["H"], errors="coerce").values, } ) + phys["g-r"] = 0.5 + if comet_photometry == "mpc-hm": + phys = add_default_sorcha_physical_columns(phys) + phys["H_r"] = phys["H_r"].where(phys["H_r"].map(math.isfinite)) + slope, assumed = comet_slopes(df) + phys["comet_G"] = slope.to_numpy() + phys["comet_G_assumed"] = assumed.to_numpy() # phys.to_csv(physparams_out, index=False) return orbs, phys diff --git a/tests/data/comet_external.ini b/tests/data/comet_external.ini new file mode 100644 index 0000000..688517b --- /dev/null +++ b/tests/data/comet_external.ini @@ -0,0 +1,28 @@ +[INPUT] +ephemerides_type = external +eph_format = csv +size_serial_chunk = 5000 +aux_format = csv +pointing_sql_query = SELECT observationId, observationStartMJD as observationStartMJD_TAI, visitTime, visitExposureTime, band as optFilter, seeingFwhmGeom as seeingFwhmGeom_arcsec, seeingFwhmEff as seeingFwhmEff_arcsec, fiveSigmaDepth as fieldFiveSigmaDepth_mag , fieldRA as fieldRA_deg, fieldDec as fieldDec_deg, rotSkyPos as fieldRotSkyPos_deg FROM observations order by observationId + +[SIMULATION] + +[FILTERS] +observing_filters = r,g,i,z,u,y + +[PHASECURVES] +phase_function = none + +[FOV] +camera_model = none + +[OUTPUT] +output_format = csv +output_columns = basic + +[LIGHTCURVE] +lc_model = none + +[ACTIVITY] +comet_activity = none + diff --git a/tests/ponder/test_comet_photometry.py b/tests/ponder/test_comet_photometry.py new file mode 100644 index 0000000..cf8d6bd --- /dev/null +++ b/tests/ponder/test_comet_photometry.py @@ -0,0 +1,429 @@ +import configparser +import json +from pathlib import Path + +import numpy as np +import pandas as pd +import pytest + +from ponder import runner +from ponder.comet_photometry import ( + complex_parameters_path, + configured_model, + metadata_signature, + model_metadata, + resolve_model, + write_physical_inputs, +) +from ponder.utils import comets_to_sorcha_inputs, diff_objects + + +def comet(name="A", **overrides): + obj = dict( + Designation_and_name=name, + Provisional_packed_desig=None, + Perihelion_dist=2.0, + e=0.5, + i=1.0, + Peri=2.0, + Node=3.0, + Epoch_year=2026, + Epoch_month=1, + Epoch_day=1, + Year_of_perihelion=2026, + Month_of_perihelion=2, + Day_of_perihelion=1.5, + H=12.0, + G=6.0, + ) + obj.update(overrides) + return obj + + +def config(tmp_path): + path = tmp_path / "input.ini" + path.write_text( + "[PHASECURVES]\nphase_function=none\n[ACTIVITY]\ncomet_activity=none\n" + "[LIGHTCURVE]\nlc_model=none\n[SATURATION]\nbright_limit=16\n" + ) + return path + + +def test_conversion_retains_alignment_and_flags_missing_g(): + objects = [ + comet("bad", Epoch_year=None), + comet("known"), + comet("missing", G=None), + comet("invalid", G="bad"), + comet("infinite", G=np.inf), + comet("unknown", H=None), + ] + orbits, phys = comets_to_sorcha_inputs(objects, [o["Designation_and_name"] for o in objects]) + assert ( + orbits.ObjID.tolist() == phys.ObjID.tolist() == ["known", "missing", "invalid", "infinite", "unknown"] + ) + assert phys.comet_G.tolist() == [6, 4, 4, 4, 6] + assert phys.comet_G_assumed.tolist() == [False, True, True, True, False] + assert np.isnan(phys.H_r.iloc[-1]) + _, legacy = comets_to_sorcha_inputs(objects, ["known"], "legacy") + assert "comet_G" not in legacy + + +def test_sidecars_for_parent_and_debug_chunks(tmp_path, monkeypatch): + monkeypatch.setattr(runner, "WORK_DIR", tmp_path / "work") + monkeypatch.setattr(runner, "RESULTS_DIR", tmp_path / "results") + objects = [comet(str(n), G=None if n == 2 else 6) for n in range(4)] + orbits, phys = comets_to_sorcha_inputs(objects, list("0123"), "mpc-hm") + parent = runner.plan_sorcha_chunks("new", "2026-10-31", 4, 2, "abcd")[1] + for chunk in [parent, *runner.plan_debug_subchunks(parent, 1)]: + runner.write_chunk_inputs(chunk, orbits, phys) + ordinary = pd.read_csv(chunk.physparams_path, dtype={"ObjID": str}) + extra = pd.read_csv(complex_parameters_path(chunk.physparams_path), dtype={"ObjID": str}) + assert ordinary.ObjID.tolist() == extra.ObjID.tolist() + assert "comet_G" not in ordinary + assert extra.comet_G.tolist() == phys.iloc[chunk.row_start : chunk.row_end].comet_G.tolist() + + +def test_effective_config_and_conflicts(tmp_path): + source = config(tmp_path) + original = source.read_text() + effective = runner.prepare_sorcha_config_for_mode(source, True, tmp_path) + parsed = configparser.ConfigParser() + parsed.read(effective) + assert configured_model(effective) == "mpc-hm" + assert parsed.getboolean("EXPERT", "brute_force") + assert not parsed.has_option("SATURATION", "bright_limit") + assert source.read_text() == original + for section, key, value in [ + ("PHASECURVES", "phase_function", "HG"), + ("ACTIVITY", "comet_activity", "lsst_comet"), + ("LIGHTCURVE", "lc_model", "sinusoidal"), + ]: + parser = configparser.ConfigParser() + parser.read_string(original) + parser.set(section, key, value) + with source.open("w") as handle: + parser.write(handle) + with pytest.raises(ValueError, match="conflicts"): + runner.prepare_sorcha_config_for_mode(source, True, tmp_path, "mpc-hm") + with pytest.raises(ValueError, match="requires --comet"): + runner.prepare_sorcha_config_for_mode(source, False, tmp_path, "mpc-hm") + + +def test_photometry_hashes_and_state_do_not_reuse_legacy(tmp_path): + objects = [comet()] + baseline = runner.object_hashes(objects, include_photometry=True) + assert diff_objects(objects, baseline, include_photometry=True) == ([], [], ["A"]) + for field in ["H", "G"]: + changed = [comet(**{field: 5.0})] + assert diff_objects(changed, baseline, include_photometry=True) == ( + [], + ["A"], + [], + ) + assert runner.state_files_for_run(tmp_path / "db", True) != runner.state_files_for_run( + tmp_path / "db", True, photometry_signature="abc" + ) + + +def test_sorcha_command_supplies_complex_parameters(tmp_path, monkeypatch): + source = config(tmp_path) + effective = runner.prepare_sorcha_config_for_mode(source, True, tmp_path, "mpc-hm") + _, phys = comets_to_sorcha_inputs([comet()], ["A"], "mpc-hm") + physical = tmp_path / "phys.csv" + write_physical_inputs(phys, physical) + commands = [] + + class Process: + def __init__(self, command, **kwargs): + commands.append(command) + + def wait(self, **kwargs): + return 0 + + monkeypatch.setattr(runner.subprocess, "Popen", Process) + runner.run_sorcha( + tmp_path / "orbit.csv", + physical, + tmp_path / "out.csv", + tmp_path / "db", + effective, + ) + assert commands[0][-2:] == ["--cp", str(complex_parameters_path(physical))] + + +@pytest.mark.parametrize("chunk_size", [5000, 0]) +def test_runner_isolates_products_and_reprocesses_photometry_changes(tmp_path, monkeypatch, chunk_size): + pytest.importorskip("sorcha_addons.activity.mpc_comet") + # Real chunk writing, consolidation, and state cycle; only expensive Sorcha execution is substituted. + import sqlite3 + + monkeypatch.chdir(tmp_path) + source = config(tmp_path) + catalog = tmp_path / "catalog.json" + catalog.write_text(json.dumps([comet()])) + db = tmp_path / "pointings.db" + with sqlite3.connect(db) as con: + con.execute( + "CREATE TABLE observations (observationId INTEGER, observationStartMJD REAL, fieldRA REAL, fieldDec REAL)" + ) + con.execute("INSERT INTO observations VALUES (1, 61000., 10., -5.)") + calls = [] + + def fake_sorcha(orbits, physparams, output, db, config, timeout=None): + calls.append(output) + assert pd.read_csv(complex_parameters_path(physparams)).comet_G.iloc[0] in ( + 6, + 7, + ) + pd.DataFrame({"ObjID": ["A"], "fieldMJD_TAI": [61000.0]}).to_csv(output, index=False) + pd.DataFrame({"ObjID": ["A"], "fieldMJD_TAI": [61000.0]}).to_csv( + output.with_name(output.stem + "_ew.csv"), index=False + ) + + monkeypatch.setattr(runner, "run_sorcha", fake_sorcha) + base = dict( + db_path=db, + object_path=catalog, + config_path=source, + comet=True, + chunk_size=chunk_size, + work_dir=tmp_path / "work", + results_dir=tmp_path / "results", + ) + original_dirs = runner.WORK_DIR, runner.RESULTS_DIR + runner.run_ponder(**base) + assert len(calls) == 1 + assert "photometry" in calls[0].parts + assert (runner.WORK_DIR, runner.RESULTS_DIR) == original_dirs + runner.run_ponder(**base) + assert len(calls) == 1 + catalog.write_text(json.dumps([comet(G=7)])) + runner.run_ponder(**base) + assert len(calls) == 2 + assert "updated" in str(calls[-1]) + assert not list((tmp_path / "results").glob("*.parquet")) + metadata = list((tmp_path / "results").glob("photometry/*/photometry.json")) + assert len(metadata) == 1 + assert json.loads(metadata[0].read_text())["missing_G"] == 4 + + +def test_external_ephemeris_cli_comparison(tmp_path): + """Use real Sorcha subprocesses, readers, photometry, and detection filtering.""" + pytest.importorskip("sorcha_addons.activity.mpc_comet") + import os + import sqlite3 + import subprocess + import sys + + objects = [ + comet("A", H=20.0, G=4.0), + comet("B", H=12.0, G=6.0), + comet("C", H=None, G=None), + ] + orbits, phys = comets_to_sorcha_inputs(objects, ["A", "B", "C"], "mpc-hm") + orbits.to_csv(tmp_path / "orbits.csv", index=False) + write_physical_inputs(phys, tmp_path / "phys.csv") + # Deliberately reverse auxiliary order: the actual Sorcha reader must join by ID. + cp = complex_parameters_path(tmp_path / "phys.csv") + pd.read_csv(cp).iloc[::-1].to_csv(cp, index=False) + au = 149597870.7 + rows = [] + for obj in objects: + row = dict( + ObjID=obj["Designation_and_name"], + FieldID=1, + fieldMJD_TAI=61000.0, + Range_LTC_km=np.sqrt(3) * au, + RangeRate_LTC_km_s=0.0, + RA_deg=10.0, + RARateCosDec_deg_day=0.0, + Dec_deg=-5.0, + DecRate_deg_day=0.0, + phase_deg=30.0, + ) + for prefix, suffix in [ + ("Obj_Sun_", "_LTC_km"), + ("Obj_Sun_v", "_LTC_km_s"), + ("Obs_Sun_", "_km"), + ("Obs_Sun_v", "_km_s"), + ]: + row.update({prefix + axis + suffix: 0.0 for axis in "xyz"}) + row["Obj_Sun_x_LTC_km"] = 2 * au + row["Obs_Sun_x_km"], row["Obs_Sun_y_km"] = 0.5 * au, np.sqrt(3) / 2 * au + rows.append(row) + pd.DataFrame(rows).to_csv(tmp_path / "ephem.csv", index=False) + with sqlite3.connect(tmp_path / "pointings.db") as con: + pd.DataFrame( + [ + dict( + observationId=1, + observationStartMJD=61000.0, + visitTime=30.0, + visitExposureTime=30.0, + band="r", + seeingFwhmGeom=0.7, + seeingFwhmEff=0.8, + fiveSigmaDepth=24.0, + fieldRA=10.0, + fieldDec=-5.0, + rotSkyPos=0.0, + ) + ] + ).to_sql("observations", con, index=False) + parser = configparser.ConfigParser() + parser.read(Path(__file__).parents[1] / "data/comet_external.ini") + parser.set("INPUT", "ephemerides_type", "external") + parser["SIMULATION"] = {} + parser.set("FOV", "camera_model", "none") + parser.set("OUTPUT", "output_columns", "all") + for section in ["FADINGFUNCTION", "LINKINGFILTER", "SATURATION"]: + parser.remove_section(section) + parser["EXPERT"] = dict( + default_snr_cut="False", + randomization_on="False", + trailing_losses_on="False", + vignetting_on="False", + ) + baseline = tmp_path / "baseline.ini" + with baseline.open("w") as handle: + parser.write(handle) + active = runner.prepare_sorcha_config_for_mode(baseline, True, tmp_path, "mpc-hm") + parser.read(active) + parser.set("EXPERT", "mag_limit", "24") + cut = tmp_path / "cut-config.ini" + with cut.open("w") as handle: + parser.write(handle) + env = dict(os.environ, SORCHA_SEED="42") + # Keep source Ponder importable even if the parent test was invoked with a relative PYTHONPATH. + env["PYTHONPATH"] = str(Path(__file__).parents[2] / "src") + os.pathsep + env.get("PYTHONPATH", "") + outputs = {} + for name, cfg in [("legacy", baseline), ("mpc_hm", active), ("cut", cut)]: + command = [ + sys.executable, + "-m", + "sorcha_cmdline.run", + "-c", + str(cfg), + "--ob", + str(tmp_path / "orbits.csv"), + "-p", + str(tmp_path / "phys.csv"), + "--cp", + str(cp), + "--pd", + str(tmp_path / "pointings.db"), + "--er", + str(tmp_path / "ephem.csv"), + "-o", + str(tmp_path), + "-t", + name, + "-f", + ] + proc = subprocess.run(command, env=env, capture_output=True, text=True, timeout=60) + assert proc.returncode == 0, proc.stdout + proc.stderr + outputs[name] = pd.read_csv(tmp_path / f"{name}.csv").set_index("ObjID") + legacy, active = outputs["legacy"], outputs["mpc_hm"] + assert set(active.index) == {"A", "B", "C"} + assert np.isnan(active.loc["C", "trailedSourceMagTrue"]) + assert bool(active.loc["C", "comet_G_assumed"]) + assert active.loc["A", "trailedSourceMagTrue"] - legacy.loc["A", "trailedSourceMagTrue"] == pytest.approx( + 5 * np.log10(2) - 2.5 * np.log10(0.4142), abs=1e-6 + ) + assert set(outputs["cut"].index) == {"B"} + + +def test_cli_explicit_model_is_forwarded_and_requires_comet(tmp_path, monkeypatch): + import sys + + from ponder import __main__ as cli + + calls = [] + monkeypatch.setattr(cli, "run_ponder", lambda *args, **kwargs: calls.append(kwargs)) + args = [ + "ponder", + "--comet-photometry", + "mpc-hm", + "--orbits", + "orbits.json", + "--db", + "db.sqlite", + "--config", + "config.ini", + ] + monkeypatch.setattr(sys, "argv", args + ["--comet"]) + cli.main() + assert calls[0]["comet_photometry"] == "mpc-hm" + monkeypatch.setattr(sys, "argv", args) + with pytest.raises(ValueError, match="requires --comet"): + cli.main() + + +@pytest.mark.parametrize( + "population, options, expected", + [ + (["--comet"], [], "mpc-hm"), + ([], [], "legacy"), + (["--comet"], ["--comet-photometry", "legacy"], "legacy"), + ], +) +def test_cli_default_follows_population(monkeypatch, population, options, expected): + import sys + + from ponder import __main__ as cli + + calls = [] + monkeypatch.setattr(cli, "run_ponder", lambda *args, **kwargs: calls.append(kwargs)) + monkeypatch.setattr( + sys, + "argv", + ["ponder", "--orbits", "orbits.json", "--db", "db.sqlite", "--config", "config.ini"] + + population + + options, + ) + cli.main() + assert calls[0]["comet_photometry"] == expected + + +def test_python_model_defaults_and_explicit_override(): + assert resolve_model(None, comet=True) == "mpc-hm" + assert resolve_model(None, comet=False) == "legacy" + assert resolve_model("legacy", comet=True) == "legacy" + with pytest.raises(ValueError, match="requires --comet"): + resolve_model("mpc-hm", comet=False) + + +def test_missing_addon_fails_before_creating_run_outputs(tmp_path, monkeypatch): + import sys + + monkeypatch.setitem(sys.modules, "sorcha_addons.activity.mpc_comet", None) + with pytest.raises(RuntimeError, match="Install the matching sorcha-addons checkout"): + runner.run_ponder( + tmp_path / "db", + tmp_path / "catalog", + tmp_path / "config", + comet=True, + work_dir=tmp_path / "work", + results_dir=tmp_path / "results", + ) + assert not (tmp_path / "work").exists() + assert not (tmp_path / "results").exists() + + +@pytest.mark.parametrize("missing_epoch", [None, "", "nan"]) +def test_catalog_alignment_keeps_unknown_h_and_skips_invalid_epochs(missing_epoch): + objects = [comet("invalid", Epoch_year=missing_epoch), comet("unknown", H=None)] + orbits, physical, catalog = runner.build_id_set_inputs(objects, ["invalid", "unknown"], "new", True) + assert orbits.ObjID.tolist() == physical.ObjID.tolist() == ["unknown"] + assert catalog.Designation_and_name.tolist() == ["unknown"] + assert np.isnan(physical.H_r.iloc[0]) + + +def test_legacy_config_and_asteroids_remain_available(tmp_path): + source = config(tmp_path) + assert runner.prepare_sorcha_config_for_mode(source, False, tmp_path) == source + legacy = runner.prepare_sorcha_config_for_mode(source, True, tmp_path, "legacy") + assert configured_model(legacy) == "legacy" + _, physical = comets_to_sorcha_inputs([comet()], ["A"], "legacy") + assert physical.columns.tolist() == ["ObjID", "H_r", "g-r"] From 8a06f9ad8d1e43f2bc7d4af0732052a73d4c0708 Mon Sep 17 00:00:00 2001 From: Colin Orion Chandler Date: Fri, 6 Nov 2026 18:33:25 -0800 Subject: [PATCH 2/3] Pin the tested comet add-on until its upstream release --- .github/workflows/comet-photometry.yml | 8 ++------ README.md | 20 +++++++++++++------- pyproject.toml | 3 ++- 3 files changed, 17 insertions(+), 14 deletions(-) diff --git a/.github/workflows/comet-photometry.yml b/.github/workflows/comet-photometry.yml index 935eb9a..18dfe0e 100644 --- a/.github/workflows/comet-photometry.yml +++ b/.github/workflows/comet-photometry.yml @@ -17,12 +17,8 @@ jobs: with: python-version: '3.13' - uses: astral-sh/setup-uv@v7 - - name: Install companion add-on and Ponder - # Temporary immutable dependency for the draft PR. Replace with the - # released sorcha-addons minimum version before merging the default. - run: | - uv pip install --system 'sorcha-addons @ git+https://github.com/ColinOrionChandler/sorcha-addons.git@7cfc0e9f173cc194379991f75d468c4cc482cd13' - uv pip install --system -e '.[dev]' + - name: Install Ponder and its pinned add-on + run: uv pip install --system -e '.[dev]' - name: Verify add-on availability and run tests run: | python -c 'from sorcha_addons.activity.mpc_comet import MPCCometActivity' diff --git a/README.md b/README.md index 34a24b6..074e4ec 100644 --- a/README.md +++ b/README.md @@ -8,19 +8,21 @@ brightness exponent and the Schleicher–Marcus (Halley–Marcus) dust phase fun by default. Asteroid runs retain their existing photometry. Select `--comet-photometry legacy` explicitly to reproduce the old distance-only law. -Ponder requires sorcha-addons. Until the new model is released, install **the -sorcha-addons checkout containing `activity/mpc_comet`** into the same Python -environment as Ponder; an older add-ons release is insufficient. Ponder checks -for the model before starting a comet run and reports an error if it is missing. -For a local checkout: +Ponder temporarily pins sorcha-addons to the tested commit containing +`activity/mpc_comet`. Installing this Ponder branch installs that exact add-on +revision automatically (Git and network access are required): ```bash -python -m pip install -e /path/to/sorcha-addons +python -m pip install -e . ponder --comet \ --config /path/to/sorcha.ini --orbits /path/to/comets.json \ --db /path/to/pointings.db --work-dir /path/to/work --results-dir /path/to/results ``` +Replace the source pin with a released minimum version once the model is +available upstream, before publishing a Ponder package release. Ponder checks +for the model before starting a comet run and reports an error if it is missing. + The add-on computes `m = H_filter + 5 log10(delta) + 2.5 G log10(r_h) - 2.5 log10(Phi_HM(alpha))` inside Sorcha, before observational noise and selection. Ponder treats MPC `H` as `H_r`, preserves `g-r = 0.5`, and supplies @@ -49,7 +51,11 @@ reruns. Use the model-specific results directory when consuming or publishing re This first-order total-magnitude model assumes dust-dominated light and treats H as a zero-phase normalization. Catalog H/G may already absorb phase or aperture effects. It does not recalibrate bands, model outbursts/gas emission, or convert -integrated coma brightness to point-source photometry. See the add-on's +integrated coma brightness to point-source photometry. Steep catalog slopes +extrapolated to large heliocentric distances can predict extremely faint total +magnitudes; this model has no separate nucleus contribution. Numerical and +integration tests verify the calculation, but its accuracy against observed +comet photometry has not yet been established. See the add-on's `docs/community_modules/activity/mpc_comet.rst` for references and examples. diff --git a/pyproject.toml b/pyproject.toml index fd37f67..e94ff85 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -20,7 +20,8 @@ dependencies = [ "astropy", "pandas < 3.0", "sorcha", - "sorcha-addons", + # Temporary source pin until the MPC comet model is released upstream. + "sorcha-addons @ git+https://github.com/ColinOrionChandler/sorcha-addons.git@7cfc0e9f173cc194379991f75d468c4cc482cd13", "tqdm", "pyarrow", ] From 2eddee1cb0c05cbf306a6af961d60cf6db44748b Mon Sep 17 00:00:00 2001 From: Colin Orion Chandler Date: Sun, 8 Nov 2026 11:23:37 -0800 Subject: [PATCH 3/3] Document comet magnitude uncertainty and validation guidance --- README.md | 8 + docs/comet-photometry.md | 165 +++++++++++++ .../catalog_records.json | 233 ++++++++++++++++++ .../comet-photometry-20261108/comparison.csv | 13 + .../comet-photometry-20261108/summary.json | 31 +++ 5 files changed, 450 insertions(+) create mode 100644 docs/comet-photometry.md create mode 100644 docs/validation/comet-photometry-20261108/catalog_records.json create mode 100644 docs/validation/comet-photometry-20261108/comparison.csv create mode 100644 docs/validation/comet-photometry-20261108/summary.json diff --git a/README.md b/README.md index 074e4ec..9498870 100644 --- a/README.md +++ b/README.md @@ -40,6 +40,14 @@ Comet runs bypass the saturation cut. Subsequent detection filters remain as configured; disable magnitude/SNR, fading, and linking cuts and observational effects to retain unknown magnitudes in geometry-only output. +**For completeness-focused comet searches, consider retaining all geometrically +eligible comets without excluding any on predicted magnitude, including unknown +magnitudes.** Catalog H/G uncertainties can produce errors of many magnitudes. +The default comet model does not disable all downstream selection. See the +[comet photometry guidance](docs/comet-photometry.md) for configuration advice, +the Horizons/Kete spot-check results, and recommended future checks against +both independent models and observed photometry. + New-model work and results are isolated below each supplied root in `photometry/mpc-hm-/`. The signature includes the add-on implementation, model version, input configuration, color assumptions, and numerical dependency versions. The result directory contains diff --git a/docs/comet-photometry.md b/docs/comet-photometry.md new file mode 100644 index 0000000..a3c067e --- /dev/null +++ b/docs/comet-photometry.md @@ -0,0 +1,165 @@ +# Comet magnitudes: validation and selection guidance + +Ponder's default comet photometry uses MPC H/G and the Schleicher–Marcus +(Halley–Marcus, HM) dust phase function through sorcha-addons. The numerical +implementation has passed integration and independent model checks, but its +accuracy against observed comet photometry has not yet been established. +Predicted magnitudes should be treated as uncertain estimates. + +## Consider retaining all comets regardless of predicted magnitude + +For completeness-focused searches, archive searches, and cutout generation, +consider retaining **all geometrically eligible comets without any exclusion +based on predicted magnitude**, including missing/unknown magnitudes. Use the +predictions as annotations for inspection rather than as a guarantee that a +comet is too faint to be worth searching for. Retaining candidates increases +processing and storage costs and does not mean that they will be detectable. + +MPC H/G values may be uncertain, missing, or inappropriate for an epoch or +heliocentric distance. Catalog fits may absorb phase and aperture effects. +Activity changes, outbursts, gas emission, and an inactive nucleus are not +modeled. Ponder currently treats catalog H as zero-phase H_r and applies fixed +Rubin colors; these are assumptions, not measured comet-specific band +transformations. A total coma magnitude also does not directly predict a +point-source detection or the surface brightness in a particular aperture. + +There is no validated universal magnitude margin that makes a hard cutoff safe. +If a detection simulation needs selection cuts, retain a separate unfiltered +candidate baseline, document the cuts, and measure how the selected sample +changes under plausible photometric assumptions. Do not treat a model-based +rejection as an observational nondetection. + +### Configure selection explicitly + +For `mpc-hm`, Ponder disables the early faint-object cull with +`brute_force = True`; comet mode also bypasses the saturation cut. **Other Sorcha +detection filters still follow the supplied configuration.** `--comet` alone +does not retain every candidate, and `--no-filter-orbits` does not disable +photometric selection. + +For a baseline without photometric rejection, the following settings mirror the +Sorcha 1.2.0 external-ephemeris integration check: + +```ini +[EXPERT] +default_snr_cut = False +randomization_on = False +trailing_losses_on = False +vignetting_on = False +``` + +Edit the existing section rather than adding a duplicate. Remove `mag_limit` +and `snr_limit` from `[EXPERT]`; a large finite limit can still reject unknown +magnitudes. Omit the `[FADINGFUNCTION]` and `[LINKINGFILTER]` sections to disable +those filters in the tested version. The four settings above by themselves are +not a complete configuration or a switch that disables every selection effect. +Keep the input, simulation, and field-of-view configuration appropriate to the +search; disabling a camera footprint changes which geometrical candidates are +included. + +Inspect the effective configuration recorded by `photometry.json`, and audit +object/epoch retention on a small run containing faint and missing-H cases. +Repeat this check after configuration or Sorcha upgrades. The integration test +in [test_comet_photometry.py](../tests/ponder/test_comet_photometry.py) exercises +retention of unknown H with selection disabled and rejection after a magnitude +cut is enabled. These settings do not prevent unrelated orbit, geometry, or +execution failures from removing rows. + +## Initial model spot check: 2026-11-08 + +We checked 12 comet/epoch cases: eight saved Sorcha rows and four selected epochs +covering near-perihelion and strong forward-scattering geometries. This was a +small diagnostic sample, not a population-wide validation. The +[comparison CSV](validation/comet-photometry-20261108/comparison.csv), +[version and residual summary](validation/comet-photometry-20261108/summary.json), +and [selected original MPC records](validation/comet-photometry-20261108/catalog_records.json) +preserve the numerical findings and catalog provenance. Paths inside these +files identify the original analysis inputs; they are not portable dependencies. + +The check used the pinned sorcha-addons commit +`7cfc0e9f173cc194379991f75d468c4cc482cd13`, Sorcha 1.2.0, sbpy 0.6.0, and Kete +1.1.0. Horizons queries used Rubin observer center `X05@399`, UTC converted from +saved TAI timestamps, and quantities 9, 19, 20, and 43 (magnitudes, distances, +and true phase angle). + +| Check | Result | +|---|---| +| Same JPL parameters and geometry, phase correction disabled | All 12 match Horizons within its printed rounding; the 11 three-decimal outputs differ by at most 0.000429 mag. | +| Kete comet law with the same inputs | Agreement within 3.6e-15 mag for both JPL and MPC parameter sets. | +| Replay of eight saved Sorcha magnitudes | Agreement within 7.2e-15 mag. | +| Substitute Horizons geometry into those eight predictions | Magnitudes change by less than 0.001 mag. | +| HM correction at all 181 integer-degree phase angles | Agreement with Lowell's table within 0.000037 mag, explained by sbpy's five-significant-digit storage of phase factors. | + +The high-phase C/2023 A3 Horizons output was rounded to an integer (`6.`); the +unrounded phase-free calculation is 5.742284. Kete received vectors constructed +from Horizons geometry: this checked its independent photometry implementation, +not its orbit propagation. Numerical precision in these comparisons is **not** +an estimate of physical magnitude accuracy. + +### Catalog disagreements dominate some predictions + +These examples use identical Horizons geometry and omit HM to isolate catalog +parameter differences. They retain each catalog's H normalization and are not +calibrated band-matched comparisons: Horizons T-mag is visual, whereas Ponder +treats MPC H as H_r. + +| Comet | UTC date | MPC H / G | JPL M1 / K1 | MPC law, no HM (mag) | Horizons T-mag | +|---|---|---|---|---:|---:| +| P/2009 WX51 | 2026-05-25 | 1 / 4 (G assumed) | 20.6 / 10 | -0.043 | 19.557 | +| 38P | 2026-07-12 | 3.5 / 12 | 11.6 / 9.5 | 45.211 | 28.828 | +| C/2014 Q3 | 2025-06-01 | -0.5 / 20 | 8.8 / 22 | 73.980 | 45.317 | + +The source record for P/2009 WX51 contains H = 1 and no G; Ponder flags its G = 4 +fallback with `comet_G_assumed`. The 19.6-mag offset from JPL is already present +before the phase correction. C/2014 Q3's steep MPC G extrapolates to extremely +faint magnitudes at large heliocentric distance. The activity model has no +separate nucleus contribution to dominate when activity fades. Neither +catalog's prediction in this table is an observed measurement. + +The forward-scattering term also behaves as intended numerically: C/2023 A3 on +2024-10-09, at phase 172.088 degrees, receives a -6.100-mag HM correction. This +does not establish that every comet at that angle brightens by that amount. +[Lowell describes HM](https://asteroid.lowell.edu/comet/dustphase/) as a +first-order dust correction; applying it to total broadband light assumes that +dust dominates. + +## Repeat model and observational spot checks + +Repeat spot checks before relying on a new catalog, photometry implementation, +dependency version, or selection configuration. Include a mix of well-observed +periodic and long-period comets, near and far from perihelion, moderate and high +phase angles, extreme H/G values, missing parameters, and candidates near any +proposed selection boundary. + +1. **Match geometry and parameter conventions.** Use the same comet/apparition, + observer location, epoch/time scale, heliocentric and observer distances, and + phase definition. Record catalog versions and fitted parameters. + [Horizons](https://ssd.jpl.nasa.gov/horizons/manual.html) and + [Kete](https://caltech-ipac.github.io/kete/api/flux.html) use the standard total + law `M1 + 5 log10(delta) + K1 log10(r_h)`. To compare with Ponder's + `H + 5 log10(delta) + 2.5 G log10(r_h)`, use `M1 = H` and `K1 = 2.5 G`. + First compare identical inputs with the phase term off, then separate the + effects of differing catalogs, geometry, bands, and phase treatment. +2. **Check the phase correction separately.** Verify zero-phase normalization + and the sign of `-2.5 log10(Phi_HM)` against the + [published table](https://asteroid.lowell.edu/static/comet/dustphaseHM_table.txt). + The standard Horizons total law lacks HM; raw T-mag need not equal a + phase-corrected Ponder prediction. Record any special object-specific law. +3. **Compare with actual observations.** Use calibrated photometry with known + times, filters, magnitude systems, apertures, uncertainties, and quality + flags. Match total-coma, aperture, or nucleus measurements to the quantity + being predicted; do not silently compare integrated coma flux with a PSF + magnitude. Account for band transformations, seeing, trailing, saturation, + outbursts, and gas contamination where relevant. Reserve independent epochs + or objects for validation if observations are also used to fit H/G. +4. **Evaluate residuals and selection losses.** Examine predicted minus observed + magnitudes against distance, phase angle, epoch, and filter. Include measured + upper limits and inspect candidates that a proposed cutoff would reject, + rather than checking only surviving detections. Quantify the effect of + missing/assumed parameters and plausible alternative fits. Until justified + by these checks, consider retaining the full geometrical candidate sample. +5. **Preserve the evidence.** Save source records, queries/responses, software + versions, effective configuration, flags, observations and their provenance, + comparison tables, and the adopted selection policy. A successful + model-to-model comparison establishes numerical consistency; observational + checks are needed to assess predictive accuracy. diff --git a/docs/validation/comet-photometry-20261108/catalog_records.json b/docs/validation/comet-photometry-20261108/catalog_records.json new file mode 100644 index 0000000..658d8b0 --- /dev/null +++ b/docs/validation/comet-photometry-20261108/catalog_records.json @@ -0,0 +1,233 @@ +{ + "source": "/sdf/scratch/rubin/kbmod/ponder/comets/work/comet_orbits_24-08-2026.json", + "sha256": "ed514e27e7261e02ee5340e55434cc70ad22a3ade82524a38e08df0f200d6de0", + "records": [ + { + "Orbit_type": "P", + "Provisional_packed_desig": "K09W51X", + "Year_of_perihelion": 2026, + "Month_of_perihelion": 4, + "Day_of_perihelion": 15.0685, + "Perihelion_dist": 0.802315, + "e": 0.739734, + "Peri": 118.317, + "Node": 31.5002, + "i": 9.6029, + "Epoch_year": 2025, + "Epoch_month": 11, + "Epoch_day": 21, + "H": 1.0, + "Designation_and_name": "P/2009 WX51 (Catalina)", + "Ref": "MPC194119" + }, + { + "Orbit_type": "C", + "Provisional_packed_desig": "K14Q030", + "Year_of_perihelion": 2014, + "Month_of_perihelion": 11, + "Day_of_perihelion": 21.7503, + "Perihelion_dist": 1.639983, + "e": 0.941103, + "Peri": 47.2671, + "Node": 63.2591, + "i": 90.2697, + "Epoch_year": 2025, + "Epoch_month": 11, + "Epoch_day": 21, + "H": -0.5, + "G": 20.0, + "Designation_and_name": "C/2014 Q3 (Borisov)", + "Ref": "MPC194138" + }, + { + "Orbit_type": "C", + "Provisional_packed_desig": "K22E030", + "Year_of_perihelion": 2023, + "Month_of_perihelion": 1, + "Day_of_perihelion": 12.4173, + "Perihelion_dist": 1.115308, + "e": 1.000007, + "Peri": 145.8448, + "Node": 302.494, + "i": 109.1793, + "Epoch_year": 2025, + "Epoch_month": 11, + "Epoch_day": 21, + "H": 7.5, + "G": 4.0, + "Designation_and_name": "C/2022 E3 (ZTF)", + "Ref": "MPC194179" + }, + { + "Orbit_type": "C", + "Provisional_packed_desig": "K23A030", + "Year_of_perihelion": 2024, + "Month_of_perihelion": 9, + "Day_of_perihelion": 27.798, + "Perihelion_dist": 0.391476, + "e": 1.000172, + "Peri": 308.5665, + "Node": 21.6467, + "i": 139.1036, + "Epoch_year": 2026, + "Epoch_month": 6, + "Epoch_day": 9, + "H": 6.5, + "G": 3.2, + "Designation_and_name": "C/2023 A3 (Tsuchinshan-ATLAS)", + "Ref": "MPEC 2026" + }, + { + "Orbit_type": "C", + "Provisional_packed_desig": "K24G030", + "Year_of_perihelion": 2025, + "Month_of_perihelion": 1, + "Day_of_perihelion": 13.4234, + "Perihelion_dist": 0.093619, + "e": 1.000006, + "Peri": 108.1478, + "Node": 220.3771, + "i": 116.8459, + "Epoch_year": 2025, + "Epoch_month": 11, + "Epoch_day": 21, + "H": 9.0, + "G": 4.0, + "Designation_and_name": "C/2024 G3 (ATLAS)", + "Ref": "MPC194189" + }, + { + "Orbit_type": "C", + "Provisional_packed_desig": "K25A060", + "Year_of_perihelion": 2025, + "Month_of_perihelion": 11, + "Day_of_perihelion": 8.5355, + "Perihelion_dist": 0.529837, + "e": 0.995649, + "Peri": 132.959, + "Node": 108.0982, + "i": 143.6635, + "Epoch_year": 2026, + "Epoch_month": 6, + "Epoch_day": 9, + "H": 13.2, + "G": 4.0, + "Designation_and_name": "C/2025 A6 (Lemmon)", + "Ref": "MPEC 2026" + }, + { + "Comet_num": 10, + "Orbit_type": "P", + "Year_of_perihelion": 2026, + "Month_of_perihelion": 8, + "Day_of_perihelion": 2.1149, + "Perihelion_dist": 1.417737, + "e": 0.537442, + "Peri": 195.4679, + "Node": 117.7975, + "i": 12.0273, + "Epoch_year": 2026, + "Epoch_month": 6, + "Epoch_day": 9, + "H": 5.0, + "G": 10.0, + "Designation_and_name": "10P/Tempel", + "Ref": "MPEC 2026" + }, + { + "Comet_num": 24, + "Orbit_type": "P", + "Year_of_perihelion": 2026, + "Month_of_perihelion": 1, + "Day_of_perihelion": 8.3415, + "Perihelion_dist": 1.183996, + "e": 0.708198, + "Peri": 58.4866, + "Node": 78.2695, + "i": 11.502, + "Epoch_year": 2026, + "Epoch_month": 6, + "Epoch_day": 9, + "H": 6.5, + "G": 14.0, + "Designation_and_name": "24P/Schaumasse", + "Ref": "MPEC 2026" + }, + { + "Comet_num": 29, + "Orbit_type": "P", + "Year_of_perihelion": 2019, + "Month_of_perihelion": 4, + "Day_of_perihelion": 21.4415, + "Perihelion_dist": 5.798329, + "e": 0.042886, + "Peri": 51.5164, + "Node": 312.4003, + "i": 9.355, + "Epoch_year": 2026, + "Epoch_month": 6, + "Epoch_day": 9, + "H": 4.0, + "G": 4.0, + "Designation_and_name": "29P/Schwassmann-Wachmann", + "Ref": "MPEC 2026" + }, + { + "Comet_num": 38, + "Orbit_type": "P", + "Year_of_perihelion": 2018, + "Month_of_perihelion": 11, + "Day_of_perihelion": 14.0211, + "Perihelion_dist": 1.581095, + "e": 0.859702, + "Peri": 359.6089, + "Node": 77.8651, + "i": 18.4111, + "Epoch_year": 2025, + "Epoch_month": 11, + "Epoch_day": 21, + "H": 3.5, + "G": 12.0, + "Designation_and_name": "38P/Stephan-Oterma", + "Ref": "MPC191595" + }, + { + "Comet_num": 47, + "Orbit_type": "P", + "Year_of_perihelion": 2025, + "Month_of_perihelion": 10, + "Day_of_perihelion": 28.0049, + "Perihelion_dist": 2.807453, + "e": 0.318006, + "Peri": 357.9184, + "Node": 356.881, + "i": 13.0394, + "Epoch_year": 2025, + "Epoch_month": 11, + "Epoch_day": 21, + "H": 1.0, + "G": 11.2, + "Designation_and_name": "47P/Ashbrook-Jackson", + "Ref": "MPEC 2026" + }, + { + "Comet_num": 67, + "Orbit_type": "P", + "Year_of_perihelion": 2028, + "Month_of_perihelion": 4, + "Day_of_perihelion": 9.3091, + "Perihelion_dist": 1.211474, + "e": 0.649733, + "Peri": 22.2248, + "Node": 36.3054, + "i": 3.8675, + "Epoch_year": 2025, + "Epoch_month": 11, + "Epoch_day": 21, + "H": 11.0, + "G": 4.0, + "Designation_and_name": "67P/Churyumov-Gerasimenko", + "Ref": "MPC191598" + } + ] +} diff --git a/docs/validation/comet-photometry-20261108/comparison.csv b/docs/validation/comet-photometry-20261108/comparison.csv new file mode 100644 index 0000000..7ac594b --- /dev/null +++ b/docs/validation/comet-photometry-20261108/comparison.csv @@ -0,0 +1,13 @@ +name,kind,utc,MPC_H,MPC_G,MPC_G_assumed,horizons_Tmag_rounding,horizons_r,horizons_delta,horizons_phase,horizons_Tmag,horizons_Nmag,horizons_M1,horizons_K1,horizons_M2,horizons_K2,horizons_raw,ponder_mpc_no_phase_at_horizons,phase_correction_mag,ponder_mpc_hm_at_horizons,ponder_jpl_no_phase,jpl_same_parameters_residual,catalog_parameter_difference,kete_jpl_no_phase,kete_mpc_no_phase,kete_jpl_residual,kete_mpc_residual,saved_rh,saved_delta,saved_phase,saved_filter,saved_filter_offset,saved_mag,replay_residual,saved_Hband_hm,geometry_magnitude_difference +C/2014 Q3 (Borisov),saved_sorcha,2025-06-01T03:14:41.322,-0.5,20.0,False,0.0005,22.68945885763,21.7696682740869,1.101,45.317,28.701,8.8,22.0,15.2,5.0,raw/3505409b1b1e45d2.json,73.98046596201422,0.0492333768823513,74.02969933889658,45.31739009497673,0.00039009497672992666,28.66307586703749,45.31739009497673,73.98046596201422,0.0,0.0,22.689806876618615,21.770009823700377,1.100494674072499,r,0.0,74.03004418920972,0.0,74.03004418920972,0.00034485031314090975 +38P/Stephan-Oterma,saved_sorcha,2026-07-12T03:13:48.476,3.5,12.0,False,0.0005,15.64182860008,15.0154854304907,2.9934,28.828,25.344,11.6,9.5,13.4,5.0,raw/c01bdbf54b45cef9.json,45.211322564668606,0.13125711729313433,45.34257968196174,28.828428349953896,0.00042834995389640085,16.38289421471471,28.828428349953896,45.211322564668606,0.0,0.0,15.641825014193849,15.015477912686086,2.993993139915314,r,0.0,45.34260087455103,-7.105427357601002e-15,45.34260087455103,2.1192589287011288e-05 +C/2022 E3 (ZTF),saved_sorcha,2025-09-04T08:37:48.480,7.5,4.0,False,0.0005,9.758371810818,9.38256003656044,5.6051,20.114,,10.8,4.5,,,raw/42469c4b63106763.json,22.25538037476954,0.23929293923366263,22.494673314003204,20.113804886732648,-0.00019511326735255352,2.1415754880368922,20.113804886732648,22.25538037476954,0.0,0.0,9.758420255011035,9.382614222929698,5.6057220340933,r,0.0,22.4947324314006,0.0,22.4947324314006,5.9117397395169746e-05 +10P/Tempel,saved_sorcha,2026-06-12T08:09:22.719,5.0,10.0,False,0.0005,1.514289215306,0.63097202021544,29.6973,13.552,,13.2,7.5,,,raw/3a7fcb586026e069.json,8.505271240111718,0.9513079087908266,9.456579148902545,13.551616726734894,-0.000383273265105899,-5.046345486623176,13.551616726734894,8.505271240111718,0.0,0.0,1.5142887108305538,0.6309721695708272,29.692329272619546,g,0.5,9.956482509524866,0.0,9.456482509524866,-9.663937767889763e-05 +24P/Schaumasse,saved_sorcha,2024-10-26T02:38:32.238,6.5,14.0,False,0.0005,4.231437043649,3.46419096242987,9.4693,22.31,,14.6,8.0,,,raw/94d61a9763d6a1a4.json,31.125085051955992,0.38809699534807035,31.513182047304063,22.309912191515153,-8.78084848459082e-05,8.815172860440839,22.30991219151515,31.12508505195599,-3.552713678800501e-15,-3.552713678800501e-15,4.231429804507332,3.4641919180148038,9.471330214390708,r,0.0,31.513231066614274,0.0,31.513231066614274,4.901931021095152e-05 +C/2025 A6 (Lemmon),saved_sorcha,2026-05-11T23:25:50.397,13.2,4.0,False,0.0005,3.109876649004,2.47594593810824,16.2937,16.68,,10.4,8.75,,,raw/336b0591a099a083.json,20.096137422434275,0.618606862102336,20.71474428453661,16.680208468192255,0.00020846819225539548,3.4159289542420197,16.68020846819225,20.096137422434275,-3.552713678800501e-15,0.0,3.109913343919181,2.475984955560305,16.292205589398872,r,0.0,20.71478391855375,0.0,20.71478391855375,3.9634017138467925e-05 +47P/Ashbrook-Jackson,saved_sorcha,2025-11-17T02:01:01.689,1.0,11.2,False,0.0005,2.80985362783,2.38551196797633,19.8615,17.119,,13.1,4.75,,,raw/1046e9ce95d30756.json,15.451051515469768,0.7225940197337245,16.17364553520349,17.11915555954258,0.00015555954258061888,-1.6681040440728125,17.119155559542577,15.451051515469766,-3.552713678800501e-15,-1.7763568394002505e-15,2.809854447600896,2.3855134781485337,19.86528996601445,g,0.5,16.673754829433307,3.552713678800501e-15,16.173754829433307,0.00010929422981931225 +P/2009 WX51 (Catalina),saved_sorcha,2026-05-25T23:17:26.803,1.0,4.0,True,0.0005,1.033871217987,0.57870296559633,71.5886,19.557,,20.6,10.0,,,raw/919b85ed1a91d469.json,-0.04305700891862237,1.037992257861143,0.9949352489425205,19.556942991081378,-5.700891862048252e-05,-19.6,19.556942991081378,-0.043057008918622186,0.0,1.8735013540549517e-16,1.0345956733864574,0.5775337016386844,71.56236342000287,r,0.0,0.9939580950826606,8.881784197001252e-16,0.9939580950826606,-0.0009771538598589702 +67P/Churyumov-Gerasimenko,controlled_epoch,2021-11-02T00:00:00.000,11.0,4.0,False,0.0005,1.210612740199,0.42142838889537,49.8246,11.646,,12.9,7.5,,,raw/17b7f7783b789957.json,9.953671339408913,1.1717514792997508,11.125422818708664,11.646158239379256,0.00015823937925496523,-1.6924868999703424,11.646158239379256,9.953671339408913,0.0,0.0,,,,,,,,, +29P/Schwassmann-Wachmann,controlled_epoch,2026-07-01T00:00:00.000,4.0,4.0,False,0.0005,6.31650073586,6.58708447247047,8.7104,17.796,,10.1,4.5,,,raw/29e566539b6a4439.json,16.098231674953222,0.359976565913167,16.45820824086639,17.795610644065977,-0.0003893559340220065,-1.697378969112755,17.795610644065977,16.098231674953222,0.0,0.0,,,,,,,,, +C/2023 A3 (Tsuchinshan-ATLAS),controlled_epoch,2024-10-09T00:00:00.000,6.5,3.2,False,0.5,0.494784398781,0.50652296791479,172.0882,6.0,,8.9,5.5,,,raw/d372a672a5288b39.json,2.578323692808304,-6.100178885713249,-3.521855192904945,5.742283699601845,-0.25771630039815463,-3.163960006793541,5.742283699601845,2.578323692808305,0.0,8.881784197001252e-16,,,,,,,,, +C/2024 G3 (ATLAS),controlled_epoch,2025-01-13T00:00:00.000,9.0,4.0,False,0.0005,0.096652406277,0.9473576231874,109.2733,-5.202,10.487,7.6,12.5,12.4,5.0,raw/9b9f7e3da11e9674.json,-1.2653035214951702,0.06519131854430461,-1.2001122029508655,-5.20227184454871,-0.0002718445487097654,3.9369683230535397,-5.20227184454871,-1.26530352149517,0.0,2.220446049250313e-16,,,,,,,,, diff --git a/docs/validation/comet-photometry-20261108/summary.json b/docs/validation/comet-photometry-20261108/summary.json new file mode 100644 index 0000000..c0dc3c4 --- /dev/null +++ b/docs/validation/comet-photometry-20261108/summary.json @@ -0,0 +1,31 @@ +{ + "purpose": "Independent model cross-check, not observed-flux calibration", + "center": "X05@399", + "time_conversion": "saved fieldMJD_TAI -> UTC JD", + "phase_quantity": "43 (true phase angle)", + "input_catalog": "/sdf/scratch/rubin/kbmod/ponder/comets/work/comet_orbits_24-08-2026.json", + "saved_source": "/tmp/ponder-hm-historical-comparison/mpc_hm.csv", + "n_cases": 12, + "versions": { + "sorcha-addons": "0.1.dev76+g7cfc0e9f1", + "sorcha": "1.2.0", + "sbpy": "0.6.0", + "numpy": "2.3.5", + "astropy": "7.2.0", + "kete": "1.1.0" + }, + "max_horizons_formula_residual": 0.25771630039815463, + "max_replay_residual": 7.105427357601002e-15, + "max_geometry_magnitude_difference": 0.0009771538598589702, + "band_caveat": "MPC H is treated as H_r in Ponder; Horizons T-mag is visual. Removing Ponder fixed filter offsets does not establish a calibrated V-r conversion.", + "sources": [ + "https://ssd.jpl.nasa.gov/horizons/manual.html", + "https://ssd-api.jpl.nasa.gov/doc/horizons.html", + "https://caltech-ipac.github.io/kete/api/flux.html", + "https://asteroid.lowell.edu/static/comet/dustphaseHM_table.txt" + ], + "max_kete_residual": 3.552713678800501e-15, + "max_lowell_phase_table_residual": 3.652234273543442e-05, + "max_horizons_residual_for_three_decimal_outputs": 0.00042834995389640085, + "n_horizons_within_printed_rounding": 12 +}