From 4d20893b68c4c580318d3bb9ee24aa84600a4891 Mon Sep 17 00:00:00 2001 From: Ou Ku Date: Thu, 26 Sep 2024 09:12:36 +0200 Subject: [PATCH] separate orbit_fit to utils.py --- pydepsi/io.py | 93 ++---------------------------------------------- pydepsi/utils.py | 81 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 84 insertions(+), 90 deletions(-) create mode 100644 pydepsi/utils.py diff --git a/pydepsi/io.py b/pydepsi/io.py index d4b0843..042cda5 100644 --- a/pydepsi/io.py +++ b/pydepsi/io.py @@ -5,6 +5,8 @@ import numpy as np +from pydepsi.utils import _orbit_fit + # Define constants SC_N_PATTERN = r"\s+([-+]?\d+(?:\.\d+)?(?:[eE][-+]?\d+)?)" SPPED_OF_LIGHT = 299792458.0 @@ -88,7 +90,7 @@ def read_metadata(resfile, mode="raw", **kwargs): orbit[i][j] = float(table_rows[i][j]) # Generate the orbfit dictionary - orbfit = _orbit_fit(orbit, verbose=0, plot_flag=0) + orbfit = _orbit_fit(orbit, verbose=0) # ++++ 9 - range_spacing pattern = r"rangePixelSpacing:" + SC_N_PATTERN @@ -260,92 +262,3 @@ def read_metadata(resfile, mode="raw", **kwargs): } return datewise_metadata - - -def _orbit_fit(orbit, verbose=0, plot_flag=0, der=True): - """Return a orbit_fit dict. - - Satellite state vector interpolation using Chebyshev polynomials of - 7th order (according to DLR recommendations). Function returns Chebyshev - polynomial coefficients. Use these to evaluate orbit state at given time - via function 'orbitVal'. - - input: snappy 'orbit' object (as read by 'read_metadata' function) - - CHANGE LOG - - 30/6/2023: Modified to adapt the input to a np.array Nx4 (N number of timesamples) - - 22/09/23: add the flag for derivative or not - """ - # parse masterorb: - t = orbit[:, 0] - x = orbit[:, 1] - y = orbit[:, 2] - z = orbit[:, 3] - - # interpolate orbits using Chebyshev polynomials of 7th order: - t0 = (min(t) + max(t)) / 2 - px = t - t0 # time argument px (centered around mid interval) - cx = np.polynomial.chebyshev.chebfit(px, x, 7) # position - cy = np.polynomial.chebyshev.chebfit(px, y, 7) - cz = np.polynomial.chebyshev.chebfit(px, z, 7) - - if der: - cvx = np.polynomial.chebyshev.chebder(cx) # velocity - cvy = np.polynomial.chebyshev.chebder(cy) - cvz = np.polynomial.chebyshev.chebder(cz) - else: - x_vel = orbit[:, 4] - y_vel = orbit[:, 5] - z_vel = orbit[:, 6] - - cvx = np.polynomial.chebyshev.chebfit(px, x_vel, 7) # velocity - cvy = np.polynomial.chebyshev.chebfit(px, y_vel, 7) - cvz = np.polynomial.chebyshev.chebfit(px, z_vel, 7) - - cax = np.polynomial.chebyshev.chebder(cvx) # acceleration - cay = np.polynomial.chebyshev.chebder(cvy) - caz = np.polynomial.chebyshev.chebder(cvz) - - if verbose: - # position fit residuals: - x_res = np.polynomial.chebyshev.chebval(px, cx) - x - y_res = np.polynomial.chebyshev.chebval(px, cy) - y - z_res = np.polynomial.chebyshev.chebval(px, cz) - z - x_std = np.std(x_res) - y_std = np.std(y_res) - z_std = np.std(z_res) - print(f"Orbit fit position residuals: X {x_std:.4f} m, Y {y_std:.4f} m, Z {z_std:.4f} m. ") - # velocity residuals: - vx_res = np.polynomial.chebyshev.chebval(px, np.polynomial.chebyshev.chebder(cx)) - x_vel - vy_res = np.polynomial.chebyshev.chebval(px, np.polynomial.chebyshev.chebder(cy)) - y_vel - vz_res = np.polynomial.chebyshev.chebval(px, np.polynomial.chebyshev.chebder(cz)) - z_vel - vx_std = np.std(vx_res) - vy_std = np.std(vy_res) - vz_std = np.std(vz_res) - print(f"Orbit fit velocity residuals: vX {vx_std:.4f} m/s, vY {vy_std:.4f} m/s, vZ {vz_std:.4f} m/s. ") - ### fit using regular polynomials like this: - if plot_flag: - import matplotlib.pyplot as plt - - _, ax = plt.subplots() - ax.plot(px, x_res, label="x_res") - ax.plot(px, y_res, label="y_res") - ax.plot(px, z_res, label="z_res") - ax.set_ylabel("residuals [m]") - ax.set_xlabel("t [s]") - ax.set_title("Orbit fit position residuals") - ax.legend() - - orbit_fit = dict() - orbit_fit["t0"] = t0 - orbit_fit["cx"] = cx - orbit_fit["cy"] = cy - orbit_fit["cz"] = cz - orbit_fit["cvx"] = cvx - orbit_fit["cvy"] = cvy - orbit_fit["cvz"] = cvz - orbit_fit["cax"] = cax - orbit_fit["cay"] = cay - orbit_fit["caz"] = caz - - return orbit_fit diff --git a/pydepsi/utils.py b/pydepsi/utils.py new file mode 100644 index 0000000..995fc12 --- /dev/null +++ b/pydepsi/utils.py @@ -0,0 +1,81 @@ +import numpy as np + + +def _orbit_fit(orbit, verbose=0, der=True): + """Return a orbit_fit dict. + + Modified from the "orbitFit" function: + https://github.com/Pbaz98/Caroline-Radar-Coding-Toolbox/blob/main/gecoris/geoUtils.py#L325 + + Satellite state vector interpolation using Chebyshev polynomials of + 7th order (according to DLR recommendations). Function returns Chebyshev + polynomial coefficients. Use these to evaluate orbit state at given time + via function 'orbitVal'. + + input: snappy 'orbit' object (as read by 'read_metadata' function) + + CHANGE LOG + - 30/6/2023: Modified to adapt the input to a np.array Nx4 (N number of timesamples) + - 22/09/23: add the flag for derivative or not + """ + # parse masterorb: + t = orbit[:, 0] + x = orbit[:, 1] + y = orbit[:, 2] + z = orbit[:, 3] + + # interpolate orbits using Chebyshev polynomials of 7th order: + t0 = (min(t) + max(t)) / 2 + px = t - t0 # time argument px (centered around mid interval) + cx = np.polynomial.chebyshev.chebfit(px, x, 7) # position + cy = np.polynomial.chebyshev.chebfit(px, y, 7) + cz = np.polynomial.chebyshev.chebfit(px, z, 7) + + if der: + cvx = np.polynomial.chebyshev.chebder(cx) # velocity + cvy = np.polynomial.chebyshev.chebder(cy) + cvz = np.polynomial.chebyshev.chebder(cz) + else: + x_vel = orbit[:, 4] + y_vel = orbit[:, 5] + z_vel = orbit[:, 6] + + cvx = np.polynomial.chebyshev.chebfit(px, x_vel, 7) # velocity + cvy = np.polynomial.chebyshev.chebfit(px, y_vel, 7) + cvz = np.polynomial.chebyshev.chebfit(px, z_vel, 7) + + cax = np.polynomial.chebyshev.chebder(cvx) # acceleration + cay = np.polynomial.chebyshev.chebder(cvy) + caz = np.polynomial.chebyshev.chebder(cvz) + + if verbose: + # position fit residuals: + x_res = np.polynomial.chebyshev.chebval(px, cx) - x + y_res = np.polynomial.chebyshev.chebval(px, cy) - y + z_res = np.polynomial.chebyshev.chebval(px, cz) - z + x_std = np.std(x_res) + y_std = np.std(y_res) + z_std = np.std(z_res) + print(f"Orbit fit position residuals: X {x_std:.4f} m, Y {y_std:.4f} m, Z {z_std:.4f} m. ") + # velocity residuals: + vx_res = np.polynomial.chebyshev.chebval(px, np.polynomial.chebyshev.chebder(cx)) - x_vel + vy_res = np.polynomial.chebyshev.chebval(px, np.polynomial.chebyshev.chebder(cy)) - y_vel + vz_res = np.polynomial.chebyshev.chebval(px, np.polynomial.chebyshev.chebder(cz)) - z_vel + vx_std = np.std(vx_res) + vy_std = np.std(vy_res) + vz_std = np.std(vz_res) + print(f"Orbit fit velocity residuals: vX {vx_std:.4f} m/s, vY {vy_std:.4f} m/s, vZ {vz_std:.4f} m/s. ") + + orbit_fit = dict() + orbit_fit["t0"] = t0 + orbit_fit["cx"] = cx + orbit_fit["cy"] = cy + orbit_fit["cz"] = cz + orbit_fit["cvx"] = cvx + orbit_fit["cvy"] = cvy + orbit_fit["cvz"] = cvz + orbit_fit["cax"] = cax + orbit_fit["cay"] = cay + orbit_fit["caz"] = caz + + return orbit_fit