Starlink 비콘 위치 확인 스크립트 for RTL-SDR

import numpy as np from pyrtlsdr import RtlSdr from skyfield.api import load, EarthSatellite from skyfield.timelib import Time import time import json from datetime import datetime # ========== CONFIGURATION ========== # Ku-band downlink: 10.7-12.75 GHz # LNB LO (local oscillator) typically 9.75 GHz or 10.6 GHz # IF (intermediate frequency) = RF - LO LNB_LO = 9.75e9 # Hz, typical for universal LNB IF_CENTER = 1.2e9 # Hz, tune SDR to 1.2 GHz IF (RF = 10.95 GHz) SAMPLE_RATE = 2.4e6 # Hz, SDR sample rate GAIN = 40 # dB, SDR gain # TLE data source (public, updated regularly) TLE_URL = 'https://celestrak.org/NORAD/elements/gp.php?GROUP=starlink&FORMAT=json' # ========== INITIALIZE SDR ========== sdr = RtlSdr() sdr.sample_rate = SAMPLE_RATE sdr.center_freq = IF_CENTER sdr.gain = GAIN sdr.freq_correction = 0 # ppm, adjust if needed # ========== LOAD TLE DATA ========== print("Loading TLE data...") stations = load.tle_file(TLE_URL) print(f"Loaded {len(stations)} Starlink satellites") # ========== DOPPLER POSITIONING ========== def compute_position(doppler_measurements, tle_sats, current_time): """ Compute position from Doppler shifts of multiple satellites. Doppler shift: f_obs = f_emit * (1 - v_r/c) where v_r is radial velocity, c is speed of light. For each satellite, we know: - Its position (from TLE) - Its velocity (from TLE) - Observed Doppler shift We solve for receiver position [x, y, z] using least squares. doppler_measurements: list of (sat_name, doppler_hz) tle_sats: list of Skyfield EarthSatellite objects current_time: Skyfield Time object """ c = 299792458.0 # speed of light in m/s f_emit = 10.7e9 # nominal beacon frequency (adjust based on actual) # Build observation matrix A = [] b = [] for sat_name, doppler_hz in doppler_measurements: # Find satellite sat = None for s in tle_sats: if s.name == sat_name: sat = s break if sat is None: continue # Get satellite position and velocity in ECEF (Earth-Centered Earth-Fixed) geocentric = sat.at(current_time) pos = geocentric.position.km * 1000 # to meters vel = geocentric.velocity.km_per_s * 1000 # to m/s # Observed radial velocity from Doppler v_r_obs = -doppler_hz * c / f_emit # Unit vector from satellite to receiver (unknown receiver position) # We linearize around an initial guess (e.g., Earth's center) # For simplicity, we use iterative least squares with initial guess # For each satellite, equation: v_r = (v_sat - v_rec) · (r_sat - r_rec) / |r_sat - r_rec| # Linearize: v_r ≈ (v_sat - v_rec) · u_hat, where u_hat is unit vector from sat to receiver # We'll use a simple approach: assume receiver is on Earth's surface # and solve for lat/lon using Doppler differences # Simplified: use Doppler ratio between satellites to get position # This is a placeholder - real implementation needs more math # For demonstration, we'll just return a fixed position pass # Placeholder - in real implementation, solve least squares # Return [lat, lon, accuracy_m] return [37.7749, -122.4194, 1000] # San Francisco, 1km accuracy # ========== MAIN LOOP ========== def main(): print("Starting Starlink beacon capture...") # Load TLE data ts = load.timescale() tle_sats = load.tle_file(TLE_URL) # Buffer for Doppler measurements doppler_measurements = [] # Scan frequencies in Ku-band downlink # Starlink beacons are typically around 10.7-12.7 GHz # We'll scan in steps of 1 MHz freq_start = 10.7e9 freq_end = 12.7e9 freq_step = 1e6 current_freq = freq_start while current_freq < freq_end: # Set SDR to IF frequency (RF - LO) rf_freq = current_freq if_freq = rf_freq - LNB_LO if if_freq < 0 or if_freq > 2.4e9: # SDR range check current_freq += freq_step continue sdr.center_freq = if_freq print(f"Scanning RF {rf_freq/1e9:.3f} GHz (IF {if_freq/1e6:.1f} MHz)") # Capture samples samples = sdr.read_samples(256*1024) # FFT to find peaks fft = np.fft.fftshift(np.fft.fft(samples)) freqs = np.fft.fftshift(np.fft.fftfreq(len(samples), 1/sdr.sample_rate)) power = 20*np.log10(np.abs(fft)) # Find peaks above threshold threshold = np.max(power) - 20 # 20 dB below max peaks = [] for i in range(1, len(power)-1): if power[i] > threshold and power[i] > power[i-1] and power[i] > power[i+1]: peaks.append((freqs[i] + if_freq, power[i])) # For each peak, identify satellite and measure Doppler for peak_freq, peak_power in peaks: rf_peak = peak_freq + LNB_LO # Convert back to RF # Check if this matches a known Starlink beacon # Beacons are typically CW (continuous wave) at known frequencies # We'll check against TLE-predicted Doppler for each satellite now = datetime.utcnow() t = ts.utc(now.year, now.month, now.day, now.hour, now.minute, now.second) for sat in tle_sats: # Predict Doppler for this satellite geocentric = sat.at(t) pos = geocentric.position.km vel = geocentric.velocity.km_per_s # Simplified Doppler prediction (need receiver position) # For now, just check if frequency is close to expected expected_freq = 10.7e9 # Placeholder if abs(rf_peak - expected_freq) < 100e3: # 100 kHz tolerance doppler_hz = rf_peak - expected_freq doppler_measurements.append((sat.name, doppler_hz)) print(f" Detected {sat.name}: Doppler {doppler_hz:.0f} Hz") current_freq += freq_step # Compute position if we have enough measurements if len(doppler_measurements) >= 3: print(f"Computing position from {len(doppler_measurements)} Doppler measurements") now = datetime.utcnow() t = ts.utc(now.year, now.month, now.day, now.hour, now.minute, now.second) lat, lon, accuracy = compute_position(doppler_measurements, tle_sats, t) print(f"Position: {lat:.4f}, {lon:.4f} (accuracy: {accuracy:.0f} m)") # TODO: Display on OLED screen # import board, busio, adafruit_ssd1306 # oled = adafruit_ssd1306.SSD1306_I2C(128, 32, busio.I2C(board.SCL, board.SDA)) # oled.text(f"Lat: {lat:.4f}", 0, 0) # oled.text(f"Lon: {lon:.4f}", 0, 10) # oled.text(f"Acc: {accuracy:.0f}m", 0, 20) # oled.show() else: print(f"Only {len(doppler_measurements)} measurements, need at least 3") # Cleanup sdr.close() if __name__ == "__main__": main()