محاسبه مختصات ماهواره از فایل SP3مدارات دقیق IGS

نمونه فایل SP3 (هدر و اولین اپوک):

#dP2024  1  7  0  0  0.00000000      96 ORBIT IGb14 HLM  IGS
## 2296 518400.00000000   900.00000000 60316 0.0000000000000
+   32   G01G02G03G04G05G06G07G08G09G10G11G12G13G14G15G16G17
+        G18G19G20G21G22G23G24G25G26G27G28G29G30G31G32  0  0
+         0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
++         2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
++         2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  0  0
%c G  cc GPS ccc cccc cccc cccc cccc ccccc ccccc ccccc ccccc
%f  1.2500000  1.025000000  0.00000000000  0.000000000000000
%i    0    0    0    0      0      0      0      0         0
/*  IGS FINAL ORBIT COMBINATION FROM IGNAAC
/*  IGS = International GNSS Service
/*  Computed using 12 Analysis Center solutions
/*
*  2024  1  7  0  0  0.00000000
PG01  -11044.578312  -11717.015832   21620.942143    260.887523
PG02   15501.434825   -4523.175039   21190.685038    481.610574
PG03   12052.762653   22403.839972    7804.261313    -68.379814
PG12  -14867.283145    4238.112308   21568.198472     39.498721
...
*  2024  1  7  0 15  0.00000000
PG01  -12261.345127  -10429.876543   21142.567891    260.891234
...

الگوریتم خواندن:

def read_sp3(filename, target_prn):
    """Read SP3 file and extract positions for a satellite."""
    epochs = []
    positions = []  # list of (x_km, y_km, z_km, clk_us)

    current_epoch = None

    with open(filename, 'r') as f:
        for line in f:
            # Skip header lines (start with #, +, %, /)
            if line[0] in ('#', '+', '%', '/'):
                continue

            # Epoch line
            if line[0] == '*':
                # *  2024  1  7  0  0  0.00000000
                parts = line.split()
                year  = int(parts[1])
                month = int(parts[2])
                day   = int(parts[3])
                hour  = int(parts[4])
                minute = int(parts[5])
                sec   = float(parts[6])
                current_epoch = datetime(year, month, day,
                                        hour, minute, int(sec))
                continue

            # Position line
            if line[0] == 'P':
                prn = line[1:4].strip()  # e.g. "G12"
                if prn != target_prn:
                    continue

                x_km  = float(line[4:18])
                y_km  = float(line[18:32])
                z_km  = float(line[32:46])
                clk_us = float(line[46:60])

                # Check for missing data
                if abs(x_km) > 999999 or abs(clk_us) > 999999:
                    continue

                epochs.append(current_epoch)
                positions.append((x_km, y_km, z_km, clk_us))

    return epochs, positions

فرمول درون‌یابی لاگرانژ:

هر مولفه به‌طور مستقل درون‌یابی می‌شود.

پیاده‌سازی در شبه‌کد:

def lagrange_interpolation(t_points, y_points, t_target):
    """
    Lagrange polynomial interpolation.
    t_points: list of known times (length n)
    y_points: list of known values (length n)
    t_target: desired time
    Returns: interpolated value at t_target
    """
    n = len(t_points)
    result = 0.0

    for i in range(n):
        # Compute Lagrange basis polynomial L_i(t)
        Li = 1.0
        for j in range(n):
            if j != i:
                Li *= (t_target - t_points[j]) /                       (t_points[i] - t_points[j])

        result += y_points[i] * Li

    return result


def interpolate_sp3(epochs, positions, target_time, n_points=9):
    """
    Interpolate SP3 positions using Lagrange polynomial.
    n_points: number of points to use (9 recommended)
    """
    # Convert epochs to seconds from first epoch
    t0 = epochs[0]
    t_sec = [(e - t0).total_seconds() for e in epochs]
    t_target = (target_time - t0).total_seconds()

    # Find closest n_points centered on target_time
    center_idx = min(range(len(t_sec)),
                     key=lambda i: abs(t_sec[i] - t_target))
    half = n_points // 2
    start = max(0, center_idx - half)
    end = min(len(t_sec), start + n_points)
    start = max(0, end - n_points)

    t_sub = t_sec[start:end]
    pos_sub = positions[start:end]

    # Interpolate each coordinate independently
    x = lagrange_interpolation(t_sub,
            [p[0] for p in pos_sub], t_target)
    y = lagrange_interpolation(t_sub,
            [p[1] for p in pos_sub], t_target)
    z = lagrange_interpolation(t_sub,
            [p[2] for p in pos_sub], t_target)
    clk = lagrange_interpolation(t_sub,
            [p[3] for p in pos_sub], t_target)

    return (x, y, z, clk)

داده‌های SP3 برای G12 (۹ اپوک مجاور):

Epoch (UTC)      X (km)         Y (km)         Z (km)        CLK (μs)
23:00:00      -15421.283145   3012.456308   21998.342472    39.412341
23:15:00      -15189.876234   3421.789012   21876.123456    39.435678
23:30:00      -14945.234567   3828.345678   21741.567890    39.458912
23:45:00      -14688.567890   4231.012345   21594.890123    39.482145
00:00:00      -14867.283145   4238.112308   21568.198472    39.498721
00:15:00      -14142.890123   5023.567890   21267.234567    39.528234
00:30:00      -13855.123456   5416.890123   21094.567890    39.551567
00:45:00      -13555.456789   5806.123456   20909.890123    39.574890
01:00:00      -13244.789012   6191.345678   20713.123456    39.598123

Target time: 00:07:30 (t = 450 s from 00:00:00)

محاسبه درون‌یابی لاگرانژ (فقط مولفه X نشان داده شده):

# 9 known times (seconds from 23:00:00):
t = [0, 900, 1800, 2700, 3600, 4500, 5400, 6300, 7200]

# X coordinates (km):
x = [-15421.28, -15189.88, -14945.23, -14688.57,
     -14867.28, -14142.89, -13855.12, -13555.46, -13244.79]

# Target time: 00:07:30 = 3600 + 450 = 4050 s from 23:00

t_target = 4050

# Lagrange interpolation:
# L_0(4050) = (4050-900)(4050-1800)...(4050-7200) /
#             (0-900)(0-1800)...(0-7200)
# ... (9 basis polynomials computed)

# Result:
X_interp ≈ -14508.156 km
Y_interp ≈ 4632.845 km
Z_interp ≈ 21419.712 km
CLK_interp ≈ 39.513 μs

نام‌گذاری فایل‌ها (فرمت جدید IGS):

# فرمت جدید (long filename):
# AAA0OPSTYP_YYYYDDDHHMM_LEN_INT_CNT.SP3.gz
# مثال:
IGS0OPSFIN_20240070000_01D_15M_ORB.SP3.gz
│   │  │    │         │   │   │
│   │  │    │         │   │   └─ نوع: ORB=orbit
│   │  │    │         │   └─ بازه: 15M=15 دقیقه
│   │  │    │         └─ مدت: 01D=1 روز
│   │  │    └─ تاریخ: 2024, day 007, 00:00
│   │  └─ نوع: FIN=Final, RAP=Rapid, ULT=Ultra-rapid
│   └─ Operation center
└─ Agency: IGS, COD, GFZ, JPL, ...

# محصولات مهم مراکز تحلیل:
COD0OPSFIN  → CODE (Switzerland)
GFZ0OPSFIN  → GFZ Potsdam (Germany)
JPL0OPSFIN  → JPL/NASA (USA)
ESA0OPSFIN  → ESA/ESOC (Europe)
IGS0OPSFIN  → IGS Combined (بهترین دقت)

انتخاب محصول مناسب: