warehack.ing/skywalker-1
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Add Phase 1 experimenter tools: MCP server, H21cm, beacon logger, arc survey

Browse source 
Four new tools transforming the SkyWalker-1 from satellite TV receiver into
a general-purpose RF observatory:

- skywalker-mcp: FastMCP server exposing 20 tools, 4 resources, 2 prompts.
  Thread-safe DeviceBridge with motor safety (continuous drive opt-in),
  input validation on all frequency/symbol rate/step parameters,
  try/finally on TS capture, path traversal sanitization, and reduced
  lock scope so emergency motor halt isn't blocked during long surveys.

- h21cm.py: Hydrogen 21 cm drift-scan radiometer at 1420.405 MHz with
  Doppler velocity calculation, control band comparison, and CSV output.

- beacon_logger.py: Long-term Ku-band beacon SNR/AGC logger with auto-relock,
  dual CSV/JSONL output, signal handlers, and systemd unit generation.

- arc_survey.py: Multi-satellite orbital arc census with USALS motor control,
  per-slot catalog persistence, resume support, and defensive motor halt
  on all error/interrupt paths.

Documentation: experimenter's roadmap guide + 4 tool reference pages (48 pages total).
This commit is contained in:
Ryan Malloy 2026-02-17 14:45:02 -07:00
parent 6c00f941eb
commit a9dcf84c38
15 changed files with 4374 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

451
tools/arc_survey.py Normal file
View file

@ -0,0 +1,451 @@
#!/usr/bin/env python3
"""
Multi-satellite arc survey for the Genpix SkyWalker-1.
Automated "satellite census": points the dish motor to each known GEO
longitude, runs a full-band carrier survey at each position, and aggregates
results into a comprehensive sky map. The diff capability tracks changes
between survey runs.
Usage:
python arc_survey.py --observer-lon -96.8 --slots "97W,99W,101W,103W"
python arc_survey.py --observer-lon -96.8 --file slots.json
python arc_survey.py --observer-lon -96.8 --arc -120 -60 --step 3
python arc_survey.py --resume arc-survey-2026-02-17.json
The tool saves progress after each orbital slot, so interrupted surveys
can be resumed. Each slot's catalog is saved individually, and a summary
report covers the entire arc.
"""
import sys
import os
import argparse
import time
import json
from datetime import datetime, timezone
from pathlib import Path
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
from skywalker_lib import SkyWalker1, usals_angle
from survey_engine import SurveyEngine
from carrier_catalog import CarrierCatalog, CATALOG_DIR
# Common North American GEO orbital slots
NA_ORBITAL_SLOTS = {
"129W": -129.0, "125W": -125.0, "123W": -123.0, "121W": -121.0,
"119W": -119.0, "118.7W": -118.7, "116.8W": -116.8, "114.9W": -114.9,
"113W": -113.0, "111.1W": -111.1, "110W": -110.0, "107.3W": -107.3,
"105W": -105.0, "103W": -103.0, "101W": -101.0, "99W": -99.0,
"97W": -97.0, "95W": -95.0, "93W": -93.0, "91W": -91.0,
"89W": -89.0, "87W": -87.0, "85W": -85.0, "83W": -83.0,
"82W": -82.0, "79W": -79.0, "77W": -77.0, "75W": -75.0,
"72.7W": -72.7, "70W": -70.0, "67W": -67.0, "65W": -65.0,
"63W": -63.0, "61.5W": -61.5, "58W": -58.0, "55.5W": -55.5,
}
ARC_SURVEY_DIR = CATALOG_DIR.parent / "arc-surveys"
class ArcSurvey:
"""Multi-position orbital arc survey with persistence and resume."""
def __init__(self, sw: SkyWalker1, observer_lon: float,
observer_lat: float = 0.0, settle_time: float = 15.0):
self.sw = sw
self.observer_lon = observer_lon
self.observer_lat = observer_lat
self.settle_time = settle_time
def survey_slot(self, name: str, sat_lon: float,
coarse_step: float = 5.0,
band: str = "", pol: str = "",
callback=None) -> CarrierCatalog:
"""Survey a single orbital slot: move dish, wait, run survey."""
# Calculate motor angle
angle = usals_angle(self.observer_lon, sat_lon, self.observer_lat)
direction = "west" if angle < 0 else "east"
if callback:
callback("moving", 0,
f"Moving to {name} ({sat_lon:.1f}), "
f"angle {abs(angle):.1f} deg {direction}")
# Command the motor
self.sw.motor_goto_x(self.observer_lon, sat_lon)
# Wait for motor to settle (larger angles need more time)
settle = max(self.settle_time, abs(angle) * 0.3)
if callback:
callback("settling", 20, f"Settling {settle:.0f}s...")
time.sleep(settle)
# Verify we have signal (check AGC for any RF energy)
sig = self.sw.signal_monitor()
if callback:
callback("signal_check", 30,
f"AGC1={sig['agc1']}, power={sig['power_db']:.1f} dB")
# Run the six-stage survey
def survey_cb(stage, pct, msg):
overall_pct = 30 + int(pct * 0.7)
if callback:
callback(stage, overall_pct, msg)
engine = SurveyEngine(self.sw, callback=survey_cb)
catalog = engine.run_full_scan(
coarse_step=coarse_step,
ts_capture_secs=2.0,
)
catalog.name = f"{name} ({sat_lon:.1f})"
catalog.band = band
catalog.pol = pol
catalog.notes = (f"Arc survey position: {name}, "
f"observer: {self.observer_lon:.2f} lon, "
f"motor angle: {angle:.2f} deg")
if callback:
callback("complete", 100,
f"{name}: {len(catalog.carriers)} carriers, "
f"{sum(1 for c in catalog.carriers if c.locked)} locked")
return catalog
def run_arc(self, slots: list[tuple[str, float]],
coarse_step: float = 5.0,
band: str = "", pol: str = "",
save_individual: bool = True,
resume_state: dict | None = None) -> dict:
"""Survey an entire arc of orbital slots.
slots: list of (name, sat_lon) tuples
resume_state: previous arc survey state dict for resuming
Returns a complete arc survey result dict.
"""
ARC_SURVEY_DIR.mkdir(parents=True, exist_ok=True)
date_str = datetime.now().strftime("%Y-%m-%d")
# Initialize or resume state
if resume_state:
state = resume_state
completed_names = set(state.get("completed_slots", {}).keys())
else:
state = {
"started": datetime.now(timezone.utc).isoformat(),
"observer_lon": self.observer_lon,
"observer_lat": self.observer_lat,
"total_slots": len(slots),
"completed_slots": {},
"skipped_slots": {},
"summary": {
"total_carriers": 0,
"total_locked": 0,
"total_services": 0,
},
}
completed_names = set()
state_path = ARC_SURVEY_DIR / f"arc-survey-{date_str}.json"
for i, (name, sat_lon) in enumerate(slots):
if name in completed_names:
print(f" [{i+1}/{len(slots)}] Skipping {name} (already surveyed)")
continue
print(f"\n [{i+1}/{len(slots)}] Surveying {name} ({sat_lon:.1f} lon)")
def progress_cb(stage, pct, msg):
print(f" [{pct:3d}%] {stage}: {msg}")
try:
catalog = self.survey_slot(
name, sat_lon,
coarse_step=coarse_step,
band=band, pol=pol,
callback=progress_cb,
)
# Save individual catalog
if save_individual:
slot_filename = f"arc-{date_str}-{name.replace('.', '_')}.json"
cat_path = catalog.save(slot_filename)
print(f" Saved: {cat_path}")
# Update state
carrier_count = len(catalog.carriers)
locked_count = sum(1 for c in catalog.carriers if c.locked)
service_count = sum(len(c.services) for c in catalog.carriers)
state["completed_slots"][name] = {
"sat_lon": sat_lon,
"completed": datetime.now(timezone.utc).isoformat(),
"carriers": carrier_count,
"locked": locked_count,
"services": service_count,
"catalog_file": slot_filename if save_individual else None,
}
state["summary"]["total_carriers"] += carrier_count
state["summary"]["total_locked"] += locked_count
state["summary"]["total_services"] += service_count
except KeyboardInterrupt:
print(f"\n Survey interrupted at {name}")
try:
self.sw.motor_halt()
except Exception:
pass
state["interrupted_at"] = name
_save_state(state, state_path)
print(f" Motor halted. Progress saved to {state_path}")
print(f" Resume with: python arc_survey.py --resume {state_path}")
return state
except Exception as e:
print(f" Error at {name}: {e}")
try:
self.sw.motor_halt()
except Exception:
pass
state["skipped_slots"][name] = {
"sat_lon": sat_lon,
"error": str(e),
}
# Save state after each slot for resume capability
_save_state(state, state_path)
# Final summary
state["completed"] = datetime.now(timezone.utc).isoformat()
_save_state(state, state_path)
return state
def _save_state(state: dict, path: Path) -> None:
"""Save arc survey state to JSON."""
with open(path, 'w') as f:
json.dump(state, f, indent=2)
def parse_slot_string(slot_str: str) -> list[tuple[str, float]]:
"""Parse a comma-separated slot string like '97W,99W,101W'.
Accepts formats: '97W', '97.5W', '3E', '-97', '-97.5'
"""
slots = []
for part in slot_str.split(','):
part = part.strip()
if not part:
continue
if part in NA_ORBITAL_SLOTS:
slots.append((part, NA_ORBITAL_SLOTS[part]))
elif part.upper().endswith('W'):
lon = -float(part[:-1])
slots.append((part.upper(), lon))
elif part.upper().endswith('E'):
lon = float(part[:-1])
slots.append((part.upper(), lon))
else:
lon = float(part)
name = f"{abs(lon):.1f}{'W' if lon < 0 else 'E'}"
slots.append((name, lon))
return slots
def generate_arc_range(start_lon: float, stop_lon: float,
step: float) -> list[tuple[str, float]]:
"""Generate orbital slots at regular intervals across an arc."""
slots = []
lon = start_lon
while lon <= stop_lon:
name = f"{abs(lon):.1f}{'W' if lon < 0 else 'E'}"
slots.append((name, lon))
lon += step
return slots
def print_summary(state: dict) -> None:
"""Print a human-readable arc survey summary."""
print(f"\n Arc Survey Summary")
print(f" ==================")
print(f" Observer: {state['observer_lon']:.2f} lon")
print(f" Slots surveyed: {len(state['completed_slots'])} / {state['total_slots']}")
print(f" Total carriers: {state['summary']['total_carriers']}")
print(f" Total locked: {state['summary']['total_locked']}")
print(f" Total services: {state['summary']['total_services']}")
if state.get("skipped_slots"):
print(f" Skipped: {len(state['skipped_slots'])}")
print(f"\n Per-slot results:")
for name, info in sorted(state["completed_slots"].items(),
key=lambda x: x[1]["sat_lon"]):
lock_str = f"{info['locked']}/{info['carriers']}"
svc_str = f"{info['services']} svc" if info['services'] else ""
print(f" {name:>8s} ({info['sat_lon']:+7.1f}): "
f"{lock_str:>7s} locked {svc_str}")
def build_parser() -> argparse.ArgumentParser:
parser = argparse.ArgumentParser(
prog="arc_survey.py",
description="Multi-satellite arc survey for SkyWalker-1",
formatter_class=argparse.RawDescriptionHelpFormatter,
epilog="""\
examples:
# Survey specific slots (North American arc)
%(prog)s --observer-lon -96.8 --slots "97W,99W,101W,103W"
# Survey an arc range at 3-degree intervals
%(prog)s --observer-lon -96.8 --arc -120 -60 --step 3
# Load slots from a JSON file
%(prog)s --observer-lon -96.8 --file my-slots.json
# Resume an interrupted survey
%(prog)s --resume ~/.skywalker1/arc-surveys/arc-survey-2026-02-17.json
# List common North American orbital slots
%(prog)s --list-slots
slot file format (JSON):
[
{"name": "97W", "lon": -97.0},
{"name": "99W", "lon": -99.0}
]
notes:
- Motor settle time scales with angle (min 15s, + 0.3s per degree)
- Each slot takes 5-15 minutes depending on carrier density
- Progress is saved after each slot; Ctrl-C to pause safely
- Individual catalogs saved to ~/.skywalker1/surveys/
- Arc survey state saved to ~/.skywalker1/arc-surveys/
""",
)
parser.add_argument('-v', '--verbose', action='store_true')
parser.add_argument('--observer-lon', type=float,
help="Observer longitude (negative=west, e.g. -96.8)")
parser.add_argument('--observer-lat', type=float, default=0.0,
help="Observer latitude (default: 0.0)")
source = parser.add_mutually_exclusive_group()
source.add_argument('--slots', type=str,
help="Comma-separated slot list (e.g. '97W,99W,101W')")
source.add_argument('--file', type=str,
help="JSON file with slot definitions")
source.add_argument('--arc', nargs=2, type=float, metavar=('START', 'STOP'),
help="Arc range in degrees longitude")
source.add_argument('--resume', type=str,
help="Resume from a saved arc survey state file")
source.add_argument('--list-slots', action='store_true',
help="List common NA orbital slots and exit")
parser.add_argument('--step', type=float, default=3.0,
help="Step size for --arc mode (default: 3.0 degrees)")
parser.add_argument('--coarse-step', type=float, default=5.0,
help="Coarse sweep step in MHz (default: 5.0)")
parser.add_argument('--settle-time', type=float, default=15.0,
help="Minimum motor settle time in seconds (default: 15)")
parser.add_argument('--pol', type=str, default="",
help="Polarization label (H/V, for catalog metadata)")
parser.add_argument('--band', type=str, default="",
help="Band label (low/high, for catalog metadata)")
return parser
def main():
parser = build_parser()
args = parser.parse_args()
if args.list_slots:
print("Common North American GEO orbital slots:")
for name in sorted(NA_ORBITAL_SLOTS, key=lambda n: NA_ORBITAL_SLOTS[n]):
lon = NA_ORBITAL_SLOTS[name]
print(f" {name:>8s} {lon:+7.1f}")
return
# Determine slot list
resume_state = None
if args.resume:
with open(args.resume) as f:
resume_state = json.load(f)
observer_lon = resume_state["observer_lon"]
observer_lat = resume_state.get("observer_lat", 0.0)
# Reconstruct slots from state
all_slot_names = (
list(resume_state.get("completed_slots", {}).keys()) +
list(resume_state.get("skipped_slots", {}).keys())
)
# We need the original slot list — reconstruct from completed + remaining
slots = []
for name, info in resume_state.get("completed_slots", {}).items():
slots.append((name, info["sat_lon"]))
for name, info in resume_state.get("skipped_slots", {}).items():
slots.append((name, info["sat_lon"]))
# Sort by longitude
slots.sort(key=lambda x: x[1])
print(f"Resuming arc survey: {len(resume_state.get('completed_slots', {}))} "
f"of {len(slots)} slots completed")
else:
if not args.observer_lon and args.observer_lon != 0:
parser.error("--observer-lon is required (or use --resume)")
observer_lon = args.observer_lon
observer_lat = args.observer_lat
if args.slots:
slots = parse_slot_string(args.slots)
elif args.file:
with open(args.file) as f:
data = json.load(f)
slots = [(d["name"], d["lon"]) for d in data]
elif args.arc:
start, stop = sorted(args.arc)
slots = generate_arc_range(start, stop, args.step)
else:
parser.error("Specify --slots, --file, --arc, or --resume")
if not slots:
print("No orbital slots to survey", file=sys.stderr)
sys.exit(1)
print(f"Arc Survey")
print(f" Observer: {observer_lon:.2f} lon, {observer_lat:.2f} lat")
print(f" Orbital slots: {len(slots)}")
for name, lon in slots:
angle = usals_angle(observer_lon, lon, observer_lat)
direction = "W" if angle < 0 else "E"
print(f" {name:>8s} {lon:+7.1f} (motor: {abs(angle):.1f} deg {direction})")
print()
with SkyWalker1(verbose=args.verbose) as sw:
sw.ensure_booted()
survey = ArcSurvey(
sw, observer_lon, observer_lat,
settle_time=args.settle_time,
)
state = survey.run_arc(
slots,
coarse_step=args.coarse_step,
band=args.band, pol=args.pol,
resume_state=resume_state,
)
print_summary(state)
if __name__ == "__main__":
main()

376
tools/beacon_logger.py Normal file
View file

@ -0,0 +1,376 @@
#!/usr/bin/env python3
"""
Long-term satellite beacon logger for the Genpix SkyWalker-1.
Locks onto a stable Ku-band transponder and logs SNR/AGC at configurable
intervals for hours, days, or weeks. Produces propagation datasets useful
for rain fade analysis, diurnal thermal drift measurement, antenna mount
stability assessment, and ITU propagation model validation.
Usage:
python beacon_logger.py --freq 12015 --sr 20000 # log to stdout
python beacon_logger.py --freq 12015 --sr 20000 -o log.csv # log to CSV
python beacon_logger.py --freq 12015 --sr 20000 --daemon # background mode
python beacon_logger.py --generate-systemd # print unit file
The tool automatically re-locks on signal loss and logs statistics per
reporting interval (min/max/mean/stddev of SNR over each window).
"""
import sys
import os
import argparse
import time
import csv
import math
import json
import signal
from datetime import datetime, timezone
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
from skywalker_lib import SkyWalker1, MODULATIONS, MOD_FEC_GROUP, FEC_RATES
def compute_stats(values: list[float]) -> dict:
"""Compute min/max/mean/stddev for a list of measurements."""
if not values:
return {"min": 0, "max": 0, "mean": 0, "stddev": 0, "count": 0}
n = len(values)
mean = sum(values) / n
variance = sum((v - mean) ** 2 for v in values) / n if n > 1 else 0
return {
"min": round(min(values), 3),
"max": round(max(values), 3),
"mean": round(mean, 3),
"stddev": round(math.sqrt(variance), 3),
"count": n,
}
class BeaconLogger:
"""Persistent signal logger with auto-relock and statistics."""
def __init__(self, sw: SkyWalker1, freq_khz: int, sr_sps: int,
mod_index: int = 0, fec_index: int = 5,
sample_interval: float = 1.0, report_interval: float = 60.0):
self.sw = sw
self.freq_khz = freq_khz
self.sr_sps = sr_sps
self.mod_index = mod_index
self.fec_index = fec_index
self.sample_interval = sample_interval
self.report_interval = report_interval
self._running = False
self._relock_count = 0
self._total_samples = 0
def tune_and_lock(self) -> bool:
"""Tune to the beacon frequency and check for lock."""
self.sw.tune(self.sr_sps, self.freq_khz, self.mod_index, self.fec_index)
time.sleep(0.5)
sig = self.sw.signal_monitor()
return sig.get("locked", False)
def run(self, duration_secs: float, csv_path: str | None = None,
json_path: str | None = None, quiet: bool = False) -> None:
"""Main logging loop.
Samples signal at sample_interval, computes statistics over
report_interval, outputs to CSV/JSON/stdout.
"""
self._running = True
# Register signal handlers for clean shutdown
def _stop(signum, frame):
self._running = False
signal.signal(signal.SIGTERM, _stop)
signal.signal(signal.SIGINT, _stop)
# Initial tune
locked = self.tune_and_lock()
if not locked:
print(f"Warning: no lock at {self.freq_khz} kHz, will keep trying",
file=sys.stderr)
# Open CSV
csv_file = None
csv_writer = None
if csv_path:
csv_file = open(csv_path, 'w', newline='')
csv_writer = csv.writer(csv_file)
csv_writer.writerow([
"timestamp", "elapsed_s", "snr_db", "agc1", "agc2",
"power_db", "locked", "relock_count",
])
# Open JSON log (append mode, one JSON object per report line)
json_file = None
if json_path:
json_file = open(json_path, 'a')
start_time = time.time()
last_report = start_time
window_snr = []
window_power = []
window_agc1 = []
lock_count_window = 0
sample_count_window = 0
try:
while self._running and (time.time() - start_time) < duration_secs:
now = time.time()
elapsed = now - start_time
# Sample
try:
sig = self.sw.signal_monitor()
except Exception as e:
if not quiet:
print(f" USB error: {e}", file=sys.stderr)
time.sleep(self.sample_interval)
continue
self._total_samples += 1
sample_count_window += 1
snr_db = sig["snr_db"]
agc1 = sig["agc1"]
agc2 = sig["agc2"]
power_db = sig["power_db"]
locked = sig["locked"]
if locked:
lock_count_window += 1
window_snr.append(snr_db)
window_power.append(power_db)
window_agc1.append(agc1)
# Write raw sample to CSV
if csv_writer:
csv_writer.writerow([
datetime.now(timezone.utc).isoformat(),
f"{elapsed:.1f}",
f"{snr_db:.3f}",
agc1, agc2,
f"{power_db:.3f}",
int(locked),
self._relock_count,
])
csv_file.flush()
# Auto-relock
if not locked:
if not quiet:
print(f" [{elapsed:.0f}s] Signal lost, attempting relock...",
file=sys.stderr)
if self.tune_and_lock():
self._relock_count += 1
if not quiet:
print(f" [{elapsed:.0f}s] Relocked (count: {self._relock_count})",
file=sys.stderr)
# Periodic report
if now - last_report >= self.report_interval:
report = {
"timestamp": datetime.now(timezone.utc).isoformat(),
"elapsed_s": round(elapsed, 1),
"samples": sample_count_window,
"lock_pct": round(100 * lock_count_window / max(sample_count_window, 1), 1),
"snr": compute_stats(window_snr),
"power": compute_stats(window_power),
"agc1": compute_stats(window_agc1),
"relock_count": self._relock_count,
}
if not quiet:
snr_s = report["snr"]
print(f" [{elapsed:7.0f}s] SNR {snr_s['mean']:5.1f} dB "
f"(min {snr_s['min']:.1f}, max {snr_s['max']:.1f}, "
f"std {snr_s['stddev']:.2f}) "
f"lock {report['lock_pct']:.0f}% "
f"relocks {self._relock_count}")
if json_file:
json_file.write(json.dumps(report) + "\n")
json_file.flush()
# Reset window
window_snr.clear()
window_power.clear()
window_agc1.clear()
lock_count_window = 0
sample_count_window = 0
last_report = now
time.sleep(self.sample_interval)
finally:
if csv_file:
csv_file.close()
if json_file:
json_file.close()
total_elapsed = time.time() - start_time
if not quiet:
print(f"\n Session complete: {self._total_samples} samples in "
f"{total_elapsed:.0f}s, {self._relock_count} relocks")
def generate_systemd_unit(args) -> str:
"""Generate a systemd unit file for daemon operation."""
cmd_parts = ["python3", os.path.abspath(__file__)]
cmd_parts.extend(["--freq", str(args.freq)])
cmd_parts.extend(["--sr", str(args.sr)])
if args.output:
cmd_parts.extend(["--output", os.path.abspath(args.output)])
if args.json_output:
cmd_parts.extend(["--json-output", os.path.abspath(args.json_output)])
cmd_parts.extend(["--duration", str(args.duration)])
cmd_parts.extend(["--sample-interval", str(args.sample_interval)])
cmd_parts.extend(["--report-interval", str(args.report_interval)])
cmd_parts.append("--quiet")
return f"""[Unit]
Description=SkyWalker-1 Beacon Logger ({args.freq} kHz)
After=network.target
[Service]
Type=simple
ExecStart={' '.join(cmd_parts)}
Restart=on-failure
RestartSec=30
StandardOutput=journal
StandardError=journal
[Install]
WantedBy=multi-user.target
"""
def build_parser() -> argparse.ArgumentParser:
parser = argparse.ArgumentParser(
prog="beacon_logger.py",
description="Long-term satellite beacon logger for SkyWalker-1",
formatter_class=argparse.RawDescriptionHelpFormatter,
epilog="""\
examples:
%(prog)s --freq 12015 --sr 20000 # Ku-band beacon, stdout
%(prog)s --freq 12015 --sr 20000 -o beacon.csv # log to CSV
%(prog)s --freq 12015 --sr 20000 --json-output beacon.jsonl # per-minute JSON
%(prog)s --freq 12015 --sr 20000 --duration 86400 # 24-hour log
%(prog)s --freq 12015 --sr 20000 --daemon # background
%(prog)s --generate-systemd --freq 12015 --sr 20000 # print unit file
The --freq is in kHz (IF frequency), not MHz. For Ku-band with a universal
LNB at LO 10750 MHz, a transponder at 12015 MHz has IF = 12015 - 10750 = 1265 MHz,
so you'd use --freq 1265000.
For IF frequencies, multiply MHz by 1000 (e.g., 1265 MHz = 1265000 kHz).
""",
)
parser.add_argument('-v', '--verbose', action='store_true')
parser.add_argument('--freq', type=int, required=True,
help="IF frequency in kHz (e.g., 1265000 for 1265 MHz)")
parser.add_argument('--sr', type=int, default=20000000,
help="Symbol rate in sps (default: 20000000)")
parser.add_argument('--mod', type=str, default="qpsk",
help="Modulation type (default: qpsk)")
parser.add_argument('--fec', type=str, default="auto",
help="FEC rate (default: auto)")
parser.add_argument('--output', '-o', type=str, default=None,
help="CSV output file (raw samples)")
parser.add_argument('--json-output', type=str, default=None,
help="JSONL output file (per-interval statistics)")
parser.add_argument('--duration', type=float, default=3600,
help="Logging duration in seconds (default: 3600)")
parser.add_argument('--sample-interval', type=float, default=1.0,
help="Seconds between samples (default: 1.0)")
parser.add_argument('--report-interval', type=float, default=60.0,
help="Seconds between summary reports (default: 60)")
parser.add_argument('--pol', type=str, default=None, choices=['H', 'V'],
help="LNB polarization (H=18V, V=13V)")
parser.add_argument('--band', type=str, default=None, choices=['low', 'high'],
help="LNB band (low=no tone, high=22kHz)")
parser.add_argument('--daemon', action='store_true',
help="Run as daemon (suppress stdout)")
parser.add_argument('--quiet', action='store_true',
help="Suppress progress output to stderr")
parser.add_argument('--generate-systemd', action='store_true',
help="Print a systemd unit file and exit")
return parser
def main():
parser = build_parser()
args = parser.parse_args()
if args.generate_systemd:
print(generate_systemd_unit(args))
return
# Resolve modulation/FEC indices
mod_entry = MODULATIONS.get(args.mod)
if mod_entry is None:
print(f"Unknown modulation '{args.mod}'. Valid: {list(MODULATIONS.keys())}",
file=sys.stderr)
sys.exit(1)
mod_idx = mod_entry[0]
fec_group = MOD_FEC_GROUP.get(args.mod, "dvbs")
fec_table = FEC_RATES.get(fec_group, {})
fec_idx = fec_table.get(args.fec, fec_table.get("auto", 0))
quiet = args.daemon or args.quiet
with SkyWalker1(verbose=args.verbose) as sw:
sw.ensure_booted()
# Configure LNB
if args.pol:
sw.set_lnb_voltage(args.pol.upper() in ("H", "L"))
if args.band:
sw.set_22khz_tone(args.band == "high")
freq_mhz = args.freq / 1000.0
sr_msps = args.sr / 1e6
if not quiet:
print(f"Beacon Logger")
print(f" Frequency: {freq_mhz:.3f} MHz IF ({args.freq} kHz)")
print(f" Symbol rate: {sr_msps:.3f} Msps")
print(f" Modulation: {args.mod}, FEC: {args.fec}")
print(f" Sample interval: {args.sample_interval}s")
print(f" Report interval: {args.report_interval}s")
print(f" Duration: {args.duration}s ({args.duration/3600:.1f}h)")
if args.output:
print(f" CSV output: {args.output}")
if args.json_output:
print(f" JSON output: {args.json_output}")
print()
logger = BeaconLogger(
sw, args.freq, args.sr,
mod_index=mod_idx, fec_index=fec_idx,
sample_interval=args.sample_interval,
report_interval=args.report_interval,
)
logger.run(
duration_secs=args.duration,
csv_path=args.output,
json_path=args.json_output,
quiet=quiet,
)
if __name__ == "__main__":
main()

357
tools/h21cm.py Normal file
View file

@ -0,0 +1,357 @@
#!/usr/bin/env python3
"""
Hydrogen 21 cm drift-scan radiometer for the Genpix SkyWalker-1.
Detects neutral hydrogen emission at 1420.405 MHz — directly in the IF range
with no LNB required. Connect an L-band antenna (patch, helical, or horn)
directly to the F-connector.
The Milky Way's spiral arms create a velocity-dispersed emission profile
detectable even with the BCM4500's ~346 kHz RBW. Earth's rotation provides
a natural drift-scan across the sky.
Usage:
python h21cm.py # single sweep, print spectrum
python h21cm.py --drift --duration 3600 # 1-hour drift scan
python h21cm.py --drift --motor-step 5 # step motor between sweeps
python h21cm.py --output data.csv # log to CSV
The c in 21 cm stands for centimeters. The frequency (1420.405 MHz) comes from
the hyperfine transition in neutral hydrogen — when the electron's spin flips
relative to the proton. This is the most fundamental spectral line in radio
astronomy, and you can detect it with a $30 DVB-S dongle.
"""
import sys
import os
import argparse
import time
import csv
import math
from datetime import datetime, timezone
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
from skywalker_lib import SkyWalker1, agc_to_power_db
# Physical constants
H1_FREQ_MHZ = 1420.405751 # Hydrogen 21 cm rest frequency
C_KM_S = 299792.458 # Speed of light
def freq_to_velocity(freq_mhz: float) -> float:
"""Convert observed frequency to radial velocity via Doppler shift.
v = c * (f_rest - f_obs) / f_rest
Positive velocity = receding (redshifted, lower frequency).
Negative velocity = approaching (blueshifted, higher frequency).
"""
return C_KM_S * (H1_FREQ_MHZ - freq_mhz) / H1_FREQ_MHZ
def sweep_h1_band(sw: SkyWalker1, center_mhz: float = H1_FREQ_MHZ,
span_mhz: float = 4.0, step_mhz: float = 0.5,
dwell_ms: int = 50, averages: int = 1) -> dict:
"""Sweep the hydrogen line band and return power measurements.
Higher dwell_ms and multiple averages improve SNR for this weak signal.
Default 50ms dwell is 5x longer than typical satellite sweeps.
Returns dict with frequencies, powers, velocities, and statistics.
"""
start = center_mhz - span_mhz / 2
stop = center_mhz + span_mhz / 2
# Accumulate multiple sweeps for averaging
all_powers = None
for avg in range(averages):
freqs, powers, raw = sw.sweep_spectrum(
start, stop, step_mhz=step_mhz, dwell_ms=dwell_ms,
sr_ksps=1000, mod_index=0, fec_index=5,
)
if all_powers is None:
all_powers = [0.0] * len(powers)
for i in range(len(powers)):
all_powers[i] += powers[i]
# Average
avg_powers = [p / averages for p in all_powers]
# Calculate velocities
velocities = [freq_to_velocity(f) for f in freqs]
# Baseline: edges of the band should be "empty" (no hydrogen)
edge_count = max(2, len(avg_powers) // 5)
baseline = (sum(avg_powers[:edge_count]) + sum(avg_powers[-edge_count:])) / (2 * edge_count)
# Excess power above baseline
excess = [p - baseline for p in avg_powers]
# Find peak excess (the hydrogen line center)
peak_idx = max(range(len(excess)), key=lambda i: excess[i])
peak_freq = freqs[peak_idx]
peak_excess = excess[peak_idx]
peak_velocity = velocities[peak_idx]
return {
"timestamp": datetime.now(timezone.utc).isoformat(),
"freqs_mhz": freqs,
"powers_db": avg_powers,
"velocities_km_s": velocities,
"excess_db": excess,
"baseline_db": baseline,
"peak_freq_mhz": peak_freq,
"peak_excess_db": peak_excess,
"peak_velocity_km_s": peak_velocity,
"averages": averages,
"dwell_ms": dwell_ms,
}
def sweep_control_band(sw: SkyWalker1, step_mhz: float = 0.5,
dwell_ms: int = 50) -> dict:
"""Sweep a control band (1430-1434 MHz) where no hydrogen is expected.
Comparing the control band to the hydrogen band reveals whether a
detected power bump is real emission or just system noise variation.
"""
freqs, powers, _ = sw.sweep_spectrum(
1430.0, 1434.0, step_mhz=step_mhz, dwell_ms=dwell_ms,
sr_ksps=1000, mod_index=0, fec_index=5,
)
mean_power = sum(powers) / len(powers) if powers else 0
return {
"control_freqs_mhz": freqs,
"control_powers_db": powers,
"control_mean_db": mean_power,
}
def print_spectrum(result: dict, show_velocity: bool = True) -> None:
"""Print an ASCII spectrum of the hydrogen band."""
freqs = result["freqs_mhz"]
excess = result["excess_db"]
velocities = result["velocities_km_s"]
baseline = result["baseline_db"]
# Scale for display
max_excess = max(excess) if excess else 1.0
min_excess = min(excess)
span = max(max_excess - min_excess, 0.5)
print(f"\n Hydrogen 21 cm Spectrum")
print(f" Baseline: {baseline:.2f} dB | Peak excess: {result['peak_excess_db']:.2f} dB")
print(f" Peak at {result['peak_freq_mhz']:.3f} MHz ({result['peak_velocity_km_s']:+.1f} km/s)")
print()
bar_width = 50
for i in range(len(freqs)):
f = freqs[i]
e = excess[i]
v = velocities[i]
# Normalize to bar width
filled = int((e - min_excess) / span * bar_width)
filled = max(0, min(filled, bar_width))
bar = '#' * filled + '-' * (bar_width - filled)
# Mark the hydrogen rest frequency
marker = " *" if abs(f - H1_FREQ_MHZ) < 0.3 else " "
if show_velocity:
print(f" {f:8.3f} MHz {v:+7.1f} km/s [{bar}] {e:+.2f} dB{marker}")
else:
print(f" {f:8.3f} MHz [{bar}] {e:+.2f} dB{marker}")
print()
print(" * = hydrogen rest frequency (1420.405 MHz)")
def drift_scan(sw: SkyWalker1, duration_secs: float, interval_secs: float,
step_mhz: float, dwell_ms: int, averages: int,
motor_step: int, output_path: str | None) -> None:
"""Run a drift scan: repeated sweeps over time.
Earth's rotation naturally scans the sky. Each sweep captures the
hydrogen profile at the current sky position. Over hours, you trace
out the galactic plane.
"""
csv_writer = None
csv_file = None
header_written = False
if output_path:
csv_file = open(output_path, 'w', newline='')
csv_writer = csv.writer(csv_file)
start_time = time.time()
scan_num = 0
try:
while time.time() - start_time < duration_secs:
scan_num += 1
elapsed = time.time() - start_time
remaining = duration_secs - elapsed
print(f"\n--- Scan #{scan_num} (elapsed {elapsed:.0f}s, "
f"remaining {remaining:.0f}s) ---")
# Motor step between scans (for declination scanning)
if motor_step and scan_num > 1:
print(f" Stepping motor {motor_step} steps east...")
sw.motor_drive_east(motor_step)
time.sleep(1.0)
result = sweep_h1_band(sw, step_mhz=step_mhz,
dwell_ms=dwell_ms, averages=averages)
print_spectrum(result, show_velocity=True)
# Write CSV
if csv_writer:
if not header_written:
csv_writer.writerow([
"timestamp", "scan_num", "freq_mhz", "power_db",
"excess_db", "velocity_km_s", "baseline_db",
])
header_written = True
for i in range(len(result["freqs_mhz"])):
csv_writer.writerow([
result["timestamp"],
scan_num,
f"{result['freqs_mhz'][i]:.3f}",
f"{result['powers_db'][i]:.3f}",
f"{result['excess_db'][i]:.3f}",
f"{result['velocities_km_s'][i]:.1f}",
f"{result['baseline_db']:.3f}",
])
csv_file.flush()
# Wait for next scan
if remaining > interval_secs:
print(f" Next scan in {interval_secs:.0f}s...")
time.sleep(interval_secs)
except KeyboardInterrupt:
print("\n Drift scan interrupted")
finally:
if csv_file:
csv_file.close()
print(f" Data saved to {output_path}")
def build_parser() -> argparse.ArgumentParser:
parser = argparse.ArgumentParser(
prog="h21cm.py",
description="Hydrogen 21 cm line radiometer for SkyWalker-1",
formatter_class=argparse.RawDescriptionHelpFormatter,
epilog="""\
examples:
%(prog)s # single sweep, print spectrum
%(prog)s --averages 8 # 8x averaging for better SNR
%(prog)s --drift --duration 3600 # 1-hour drift scan
%(prog)s --drift --motor-step 5 # step motor between sweeps
%(prog)s --output h21cm-data.csv # log to CSV
%(prog)s --control # include control band comparison
notes:
- Connect an L-band antenna directly to the F-connector (no LNB)
- LNB power is disabled automatically for direct input
- Hydrogen emission is weak; use --averages 4-16 for best results
- The --dwell option increases per-step integration time (default 50ms)
- Earth rotation provides natural sky drift at ~15 deg/hour
""",
)
parser.add_argument('-v', '--verbose', action='store_true',
help="Show raw USB traffic")
parser.add_argument('--center', type=float, default=H1_FREQ_MHZ,
help=f"Center frequency in MHz (default: {H1_FREQ_MHZ})")
parser.add_argument('--span', type=float, default=4.0,
help="Frequency span in MHz (default: 4.0)")
parser.add_argument('--step', type=float, default=0.5,
help="Frequency step in MHz (default: 0.5)")
parser.add_argument('--dwell', type=int, default=50,
help="Dwell time per step in ms (default: 50)")
parser.add_argument('--averages', type=int, default=1,
help="Number of sweeps to average (default: 1)")
parser.add_argument('--output', '-o', type=str, default=None,
help="CSV output file path")
parser.add_argument('--control', action='store_true',
help="Include control band (1430-1434 MHz) for comparison")
parser.add_argument('--no-velocity', action='store_true',
help="Don't show velocity axis in spectrum display")
drift_group = parser.add_argument_group('drift scan')
drift_group.add_argument('--drift', action='store_true',
help="Enable drift scan mode (repeated sweeps)")
drift_group.add_argument('--duration', type=float, default=3600,
help="Drift scan duration in seconds (default: 3600)")
drift_group.add_argument('--interval', type=float, default=60,
help="Seconds between sweeps (default: 60)")
drift_group.add_argument('--motor-step', type=int, default=0,
help="Motor steps between sweeps (0=no motor, default: 0)")
return parser
def main():
parser = build_parser()
args = parser.parse_args()
with SkyWalker1(verbose=args.verbose) as sw:
sw.ensure_booted()
# Disable LNB power for direct input
sw.start_intersil(on=False)
print("LNB power disabled (direct L-band input mode)")
if args.drift:
drift_scan(sw, duration_secs=args.duration,
interval_secs=args.interval,
step_mhz=args.step, dwell_ms=args.dwell,
averages=args.averages, motor_step=args.motor_step,
output_path=args.output)
else:
# Single sweep
print(f"\nSweeping {args.center - args.span/2:.1f} - "
f"{args.center + args.span/2:.1f} MHz "
f"(step={args.step} MHz, dwell={args.dwell}ms, "
f"avg={args.averages}x)")
result = sweep_h1_band(sw, center_mhz=args.center,
span_mhz=args.span, step_mhz=args.step,
dwell_ms=args.dwell, averages=args.averages)
print_spectrum(result, show_velocity=not args.no_velocity)
if args.control:
print(" Control band (1430-1434 MHz, no hydrogen expected):")
ctrl = sweep_control_band(sw, step_mhz=args.step, dwell_ms=args.dwell)
print(f" Control mean: {ctrl['control_mean_db']:.2f} dB")
print(f" H1 baseline: {result['baseline_db']:.2f} dB")
diff = result["peak_excess_db"]
print(f" H1 peak excess above baseline: {diff:+.2f} dB")
# Write single sweep to CSV if requested
if args.output:
with open(args.output, 'w', newline='') as f:
writer = csv.writer(f)
writer.writerow([
"freq_mhz", "power_db", "excess_db",
"velocity_km_s", "baseline_db",
])
for i in range(len(result["freqs_mhz"])):
writer.writerow([
f"{result['freqs_mhz'][i]:.3f}",
f"{result['powers_db'][i]:.3f}",
f"{result['excess_db'][i]:.3f}",
f"{result['velocities_km_s'][i]:.1f}",
f"{result['baseline_db']:.3f}",
])
print(f" Data saved to {args.output}")
if __name__ == "__main__":
main()