from datetime import datetime, timedelta, timezone from skyfield.api import load from skyfield.framelib import itrs import math import numpy as np from numpy.polynomial.chebyshev import chebfit ## topocentric Moon vector # START_TIME = datetime.now(timezone.utc) # START_TIME = datetime(2026, 8, 18, 0, 0, tzinfo=timezone.utc) START_TIME = datetime(2025, 8, 18, 0, 0, tzinfo=timezone.utc) END_TIME_YEARS = 5 # PLANET = "mercury barycenter" PLANET = "moon" """fitting""" CHEB_DEGREE = 5 ## polynomial approximation window ## bigger means less accurate but less memory SEGMENT_DAYS = 12 ## all tests were done with days = time_array_segments[0] # mercury | CHEB_DEGREE = 5, SEGMENT_DAYS = 32, mean error = 0.0xx degrees, std 0.0xx degrees # venus | CHEB_DEGREE = 5, SEGMENT_DAYS = 64, mean error = 0.012 degrees, std 0.007 degrees # moon | CHEB_DEGREE = 5, SEGMENT_DAYS = 12, mean error = 0.0xx degrees, std 0.0xx degrees # mars | CHEB_DEGREE = 4, SEGMENT_DAYS = 160, mean error = 0.016 degrees, std 0.008 degrees # jupiter | CHEB_DEGREE = 4, SEGMENT_DAYS = 160, mean error = 0.018 degrees, std 0.011 degrees # saturn | CHEB_DEGREE = 4, SEGMENT_DAYS = 160, mean error = 0.012 degrees, std 0.007 degrees # uranus | CHEB_DEGREE = 3, SEGMENT_DAYS = 160, mean error = 0.018 degrees, std 0.011 degrees # neptune | CHEB_DEGREE = 3, SEGMENT_DAYS = 160, mean error = 0.014 degrees, std 0.008 degrees # pluto | CHEB_DEGREE = 3, SEGMENT_DAYS = 160, mean error = 0.014 degrees, std 0.008 degrees ## dense samples/segment for better chebyshev polynomials lsq fit SAMPLES_PER_SEGMENT = 256 TEST_RUNS = 200000 RUN_TEST = False CREATE_EPHEMERIS = True ## build time, should not need to edit code below this line ########################################################################################## ## build time array that has the next END_TIME_YEARS years in SEGMENT_DAYS segments with samples = SAMPLES_PER_SEGMENT END_TIME = START_TIME + timedelta(days=END_TIME_YEARS * 365.25) time_array_segments = [] segment_start = START_TIME # print(f"created {math.ceil(END_TIME_YEARS * 365.25 / SEGMENT_DAYS)} chebyshev segments for {PLANET}") while segment_start < END_TIME: segment_end = segment_start + timedelta(days=SEGMENT_DAYS) ## the end data point will go over by some value between 0 and SEGMENT_DAYS dt = (segment_end - segment_start) / (SAMPLES_PER_SEGMENT - 1) segment = [segment_start + i * dt for i in range(SAMPLES_PER_SEGMENT)] time_array_segments.append(segment) segment_start = segment_end ts = load.timescale() ## get earth's rotation for ECEF calculations ## R0_full is earth rotation matrix that would take a ECI vector ## and rotate it to its ECEF spot on earth ## But we cant store this many rotation matrices ## SIDEREAL_DAY_SECONDS = 86164.0905 ## 1) Get total seconds, say dt_seconds and use OMEGA = 2.0 * np.pi / SIDEREAL_DAY_SECONDS ## 2) Use theta0 to get theta = theta0 + OMEGA * dt_seconds, then c = cos(theta), s = sin(theta) ## 2.5) theta0 is how far has the Earth-fixed X axis rotated around the Earth’s ECI spin axis ## 3) Build a Rz matrix, Rz = [[c, s, 0.0],[-s, c, 0.0],[0.0, 0.0, 1.0]] ## 4) Tilt the ECI vector by R_slow, then spin fast fast with Rz ## 4.5) ecef_approx = Rz @ R_slow @ eci t0 = ts.from_datetime(START_TIME) R0_full = itrs.rotation_at(t0) theta0 = np.arctan2(R0_full[0, 1], R0_full[0, 0]) c = np.cos(theta0) s = np.sin(theta0) Rz0 = np.array([ [c, s, 0.0], [-s, c, 0.0], [0.0, 0.0, 1.0]]) R_slow = Rz0.T @ R0_full ## https://rhodesmill.org/skyfield/planets.html eph = load("de440s.bsp") earth = eph["earth"] target = eph[PLANET] ## test earth rotation # test_datetime = START_TIME + timedelta(days=4 * 365.25) # t_test = ts.from_datetime(test_datetime) # observation = earth.at(t_test).observe(target) # eci = np.asarray(observation.position.au, dtype=np.float64) # eci /= np.linalg.norm(eci) # ecef_truth = np.asarray(observation.frame_xyz(itrs).au, dtype=np.float64) # ecef_truth /= np.linalg.norm(ecef_truth) # SIDEREAL_DAY_SECONDS = 86164.0905 # OMEGA = 2.0 * np.pi / SIDEREAL_DAY_SECONDS # dt_seconds = (test_datetime - START_TIME).total_seconds() # theta = theta0 + OMEGA * dt_seconds # theta = theta % (2.0 * np.pi) # c = np.cos(theta) # s = np.sin(theta) # Rz = np.array([ # [c, s, 0.0], # [-s, c, 0.0], # [0.0, 0.0, 1.0] # ]) # ecef_approx = Rz @ R_slow @ eci # ecef_approx /= np.linalg.norm(ecef_approx) # print("date:") # print(test_datetime) # print("\nECI:") # print(eci) # print("\nECEF Skyfield:") # print(ecef_truth) # print("\nECEF approximation:") # print(ecef_approx) # dot = np.dot(ecef_truth, ecef_approx) # dot = np.clip(dot, -1.0, 1.0) # error_deg = np.degrees(np.arccos(dot)) # print("\nAngular error:") # print(error_deg, "degrees") def skyfield_target_positions(datetimes): """ normalized geocentric ICRF AU vectors to target body. this is the vector pointing from earth the the target, accounting for the speed of light (it takes 4 hours for neptunes light to reach us, so yes speed of light needs to be factored) this reference frame is inertial and does not account earth spin """ times = ts.from_datetimes(datetimes) observation = earth.at(times).observe(target) xyz = np.asarray(observation.position.au, dtype=np.float64) lengths = np.sqrt(np.sum(xyz * xyz, axis=0)) xyz /= lengths return xyz # ## saves data as representing data in int16, must remember this when we reverse calculate # days = time_array_segments[0] # xyz = skyfield_target_positions(days) # # # build ephemeris # u = np.linspace(-1.0, 1.0, len(days), dtype=np.float64) # coeff_x = chebfit(u, xyz[0], CHEB_DEGREE) # coeff_y = chebfit(u, xyz[1], CHEB_DEGREE) # coeff_z = chebfit(u, xyz[2], CHEB_DEGREE) # # ## I noticed chebyshev coefficient were all less than one, so i mapped out the data to int 16 which ranges # ## +/- 32767 to save memory. Then I found the moon at 16 day segments broke my assumption... # ## if your max |chebyshev coefficient| is over one, adjust until it works. this unfortunately is a bad solution # ## but when it works (99%) of the time saves a lot of memory # max_chebyshev_coefficient = np.max(np.abs(np.array([coeff_x, coeff_y, coeff_z]))) # # print(f"largest chebyshev coefficient is {max_chebyshev_coefficient}") # if max_chebyshev_coefficient > 1.0: # raise ValueError(f"Chebyshev coefficient exceeds is too large, adjust fitting params ̄\_(ツ)_/ ̄") # # # to get back to chebyshev coefficient for positions in AU divide coeffs by 32767.0 ## wait, everything above is still true but the moon consistently gives vlaues over 1, so increase the range ## by coeff_range > 1.0 coeff_range = 1.15 int_16_scale = 32767.0 / coeff_range # int_16_scale = 32767.0 # coeff_x_int16 = np.round(coeff_x * int_16_scale).astype(np.int16) # coeff_y_int16 = np.round(coeff_y * int_16_scale).astype(np.int16) # coeff_z_int16 = np.round(coeff_z * int_16_scale).astype(np.int16) ## coefficient array and a value between -1 and +1 def chebyshev_eval(coefficients, x): b1 = 0.0 b2 = 0.0 for i in range(len(coefficients) - 1, 0, -1): c = float(coefficients[i]) / int_16_scale b0 = (2.0 * x * b1 - b2 + c) b2 = b1 b1 = b0 c0 = float(coefficients[0]) / int_16_scale return x * b1 - b2 + c0 ## test if RUN_TEST: from tqdm import tqdm for idx, days in tqdm(enumerate(time_array_segments)): xyz = skyfield_target_positions(days) u = np.linspace(-1.0, 1.0, len(days), dtype=np.float64) coeff_x = chebfit(u, xyz[0], CHEB_DEGREE) coeff_y = chebfit(u, xyz[1], CHEB_DEGREE) coeff_z = chebfit(u, xyz[2], CHEB_DEGREE) max_chebyshev_coefficient = np.max(np.abs(np.array([coeff_x, coeff_y, coeff_z]))) if max_chebyshev_coefficient > coeff_range: raise ValueError(f"Chebyshev coefficient exceeds is too large") coeff_x_int16 = np.round(coeff_x * int_16_scale).astype(np.int16) coeff_y_int16 = np.round(coeff_y * int_16_scale).astype(np.int16) coeff_z_int16 = np.round(coeff_z * int_16_scale).astype(np.int16) fractions = np.random.random(TEST_RUNS) u_test = 2.0 * fractions - 1.0 segment_start = time_array_segments[idx][0] segment_end = time_array_segments[idx][-1] segment_duration = segment_end - segment_start test_times = [segment_start + segment_duration * f for f in fractions] ## use numpy for fast fast test coeff_x_float = coeff_x_int16.astype(np.float64) / int_16_scale coeff_y_float = coeff_y_int16.astype(np.float64) / int_16_scale coeff_z_float = coeff_z_int16.astype(np.float64) / int_16_scale approx = np.vstack([ np.polynomial.chebyshev.chebval(u_test, coeff_x_float), np.polynomial.chebyshev.chebval(u_test, coeff_y_float), np.polynomial.chebyshev.chebval(u_test, coeff_z_float), ]) ## should I re-normalize each vector? I think so because of int_16_scale scaling and independent coeff axes...? approx /= np.linalg.norm(approx, axis=0, keepdims=True) ## skyfield truth truth = skyfield_target_positions(test_times) # # error (this one doesnt mean much since the arrays have been normalized) # vector_error = np.linalg.norm(approx - truth, axis=0) # print(np.max(vector_error)) # print(np.std(vector_error)) # print(np.mean(vector_error)) ## angular error (angle between approx to truth) dot = np.sum(approx * truth, axis=0) ## same as dot product dot = np.clip(dot, -1.0, 1.0) ## be safe angular_error_rad = np.arccos(dot) angular_error_deg = np.degrees(angular_error_rad) print() print(f"{PLANET} angular errors results:") print(f"mean: {np.mean(angular_error_deg):.3f} degrees") print(f"std: {np.std(angular_error_deg):.3f} degrees") print(f"max: {np.max(angular_error_deg):.3f} degrees") ## 95% of tests have error ≤ what ever value it prints° print(f"95th pct: {np.percentile(angular_error_deg, 95):.3f} degrees") print(f"99th pct: {np.percentile(angular_error_deg, 99):.3f} degrees") # import matplotlib.pyplot as plt # plt.hist(angular_error_deg, bins=100) # plt.xlabel("Angular Error (degrees)") # plt.ylabel("Count") # plt.title(f"{PLANET} Chebyshev Angular Error") # plt.show() ## print the N largest errors will show worst errors are at Chebyshev boundaries where u~1 ## intuitively, error at boundaries makes sense, why at 1 and not -1, i have no idea ## probably some really interesting Chebyshev math # TOP_N = 50 # worst_indices = np.argsort(angular_error_deg)[-TOP_N:][::-1] # print() # print(f"Top {TOP_N} largest errors") # for i in worst_indices: # print( # f"{test_times[i]} " # f"u={u_test[i]: .6f} " # f"error={angular_error_deg[i]:.6f} deg") # ephemeris = [] # for days in time_array_segments: # xyz = skyfield_target_positions(days) # u = np.linspace(-1.0, 1.0, len(days), dtype=np.float64) # coeff_x = chebfit(u, xyz[0], CHEB_DEGREE) # coeff_y = chebfit(u, xyz[1], CHEB_DEGREE) # coeff_z = chebfit(u, xyz[2], CHEB_DEGREE) # max_coeff = np.max(np.abs([coeff_x, coeff_y, coeff_z])) # # if max_coeff > coeff_range: # raise ValueError( # f"coefficient too large: {max_coeff}" # ) # # coeff_x_int16 = np.round(coeff_x * int_16_scale).astype(np.int16) # coeff_y_int16 = np.round(coeff_y * int_16_scale).astype(np.int16) # coeff_z_int16 = np.round(coeff_z * int_16_scale).astype(np.int16) # # ephemeris.append(np.vstack([coeff_x_int16, coeff_y_int16, coeff_z_int16])) # # ephemeris = np.asarray(ephemeris, dtype=np.int16) # import json # # ephemeris_json = { # "planet": PLANET, # "start_time_utc": START_TIME.isoformat(), # "end_time_utc": END_TIME.isoformat(), # "segment_days": SEGMENT_DAYS, # "cheb_degree": CHEB_DEGREE, # "coeff_range": coeff_range, # "coeff_dtype": "int16", # "earth_theta0_rad": float(theta0), # "earth_r_slow": R_slow.tolist(), # "segments": [] # } # # for segment in ephemeris: # ephemeris_json["segments"].append({ # "x": segment[0].tolist(), # "y": segment[1].tolist(), # "z": segment[2].tolist(), # }) # # with open(f"{PLANET} ephemeris.json", "w") as f: # json.dump(ephemeris_json, f, indent=2) # # print(f"saved {PLANET} ephemeris.json") def master_ephemeris_json(filename="master_ephemeris.json"): global PLANET, target configs = { "mercury": ("mercury", 5, 32), "venus": ("venus", 5, 64), "moon": ("moon", 5, 12), "mars": ("mars barycenter", 4, 160), "jupiter": ("jupiter barycenter", 4, 160), "saturn": ("saturn barycenter", 4, 160), "uranus": ("uranus barycenter", 3, 160), "neptune": ("neptune barycenter", 3, 160), "pluto": ("pluto barycenter", 3, 160), } master = { "start_time_utc": START_TIME.isoformat(), "end_time_utc": END_TIME.isoformat(), "coeff_range": coeff_range, "coeff_dtype": "int16", "earth_theta0_rad": float(theta0), "earth_r_slow": R_slow.tolist(), "bodies": {} } for name, (skyfield_name, cheb_degree, segment_days) in configs.items(): print(f"building {name}...") PLANET = skyfield_name target = eph[skyfield_name] segments = [] segment_start = START_TIME while segment_start < END_TIME: segment_end = segment_start + timedelta(days=segment_days) dt = ((segment_end - segment_start) / (SAMPLES_PER_SEGMENT - 1)) days = [segment_start + i * dt for i in range(SAMPLES_PER_SEGMENT)] xyz = skyfield_target_positions(days) u = np.linspace(-1.0, 1.0, len(days), dtype=np.float64) coeff_x = chebfit(u, xyz[0], cheb_degree) coeff_y = chebfit(u, xyz[1], cheb_degree) coeff_z = chebfit(u, xyz[2], cheb_degree) max_coeff = np.max(np.abs([coeff_x, coeff_y, coeff_z])) if max_coeff > coeff_range: raise ValueError(f"{name}: coefficient too large: {max_coeff}") coeff_x_int16 = np.round(coeff_x * int_16_scale).astype(np.int16) coeff_y_int16 = np.round(coeff_y * int_16_scale).astype(np.int16) coeff_z_int16 = np.round(coeff_z * int_16_scale).astype(np.int16) segments.append({ "x": coeff_x_int16.tolist(), "y": coeff_y_int16.tolist(), "z": coeff_z_int16.tolist(), }) segment_start = segment_end master["bodies"][name] = { "skyfield_name": skyfield_name, "segment_days": segment_days, "cheb_degree": cheb_degree, "segment_count": len(segments), "segments": segments, } with open(filename, "w") as f: json.dump(master, f, indent=2) print(f"saved {filename}") return master def master_ephemeris_c(filename="master_ephemeris.c"): global PLANET, target configs = { "mercury": ("mercury", 5, 32), "venus": ("venus", 5, 64), "moon": ("moon", 5, 12), "mars": ("mars barycenter", 4, 160), "jupiter": ("jupiter barycenter", 4, 160), "saturn": ("saturn barycenter", 4, 160), "uranus": ("uranus barycenter", 3, 160), "neptune": ("neptune barycenter", 3, 160), "pluto": ("pluto barycenter", 3, 160), } bodies = {} for name, (skyfield_name, cheb_degree, segment_days) in configs.items(): print(f"building {name}...") PLANET = skyfield_name target = eph[skyfield_name] body_segments = [] segment_start = START_TIME while segment_start < END_TIME: segment_end = segment_start + timedelta(days=segment_days) dt = ((segment_end - segment_start) / (SAMPLES_PER_SEGMENT - 1)) days = [segment_start + i * dt for i in range(SAMPLES_PER_SEGMENT)] xyz = skyfield_target_positions(days) u = np.linspace(-1.0, 1.0, len(days), dtype=np.float64) coeff_x = chebfit(u, xyz[0], cheb_degree) coeff_y = chebfit(u, xyz[1], cheb_degree) coeff_z = chebfit(u, xyz[2], cheb_degree) max_coeff = np.max(np.abs([coeff_x, coeff_y, coeff_z])) if max_coeff > coeff_range: raise ValueError(f"{name}: coefficient too large: {max_coeff}") coeffs = np.vstack([ np.round(coeff_x * int_16_scale).astype(np.int16), np.round(coeff_y * int_16_scale).astype(np.int16), np.round(coeff_z * int_16_scale).astype(np.int16), ]) body_segments.append(coeffs) segment_start = segment_end bodies[name] = { "degree": cheb_degree, "segment_days": segment_days, "segments": body_segments, } ## number of days from START_TIME to the first day of each month ## lets the STM32 convert year/month/day into ephemeris time in sec start_year = START_TIME.year end_year = END_TIME.year month_offset_rows = [] for year in range(start_year, end_year + 1): row = [] for month in range(1, 13): month_start = datetime( year, month, 1, tzinfo=timezone.utc ) days_from_start = (month_start - START_TIME).days row.append(days_from_start) month_offset_rows.append(row) with open(filename, "w") as f: f.write("#include \n") f.write("#include \n") f.write("#include \n\n") # calendar data f.write(f"#define EPHEMERIS_START_YEAR {start_year}\n") f.write(f"#define EPHEMERIS_END_YEAR {end_year}\n") f.write( f"#define EPHEMERIS_YEAR_COUNT " f"{len(month_offset_rows)}\n\n") f.write( "static const int16_t " "ephemeris_month_start_days" "[EPHEMERIS_YEAR_COUNT][12] = {\n") for year, row in zip(range(start_year, end_year + 1), month_offset_rows): values = ", ".join(str(v) for v in row) f.write(f" {{{values}}}, // {year}\n") f.write("};\n\n") f.write(f"#define EPHEMERIS_COEFF_RANGE {coeff_range:.17g}\n") f.write( "#define EPHEMERIS_INT16_SCALE " "(32767.0 / EPHEMERIS_COEFF_RANGE)\n\n" ) f.write( f"static const float earth_theta0_rad = " f"{float(theta0):.17g};\n\n" ) f.write("static const float earth_r_slow[3][3] = {\n") for row in R_slow: f.write( " {" + ", ".join(f"{float(v):.17g}" for v in row) + "},\n" ) f.write("};\n\n") for name, body in bodies.items(): degree = body["degree"] coeff_count = degree + 1 segments = body["segments"] f.write( f"// {name}: degree={degree}, " f"segment_days={body['segment_days']}\n" ) f.write( f"static const int16_t {name}_ephemeris" f"[{len(segments)}][3][{coeff_count}] = {{\n" ) for segment in segments: f.write(" {\n") for axis in range(3): values = ", ".join( str(int(v)) for v in segment[axis] ) f.write(f" {{{values}}},\n") f.write(" },\n") f.write("};\n\n") f.write( "typedef struct {\n" " const void *coefficients;\n" " uint32_t segment_count;\n" " uint16_t segment_days;\n" " uint8_t cheb_degree;\n" "} ephemeris_body_t;\n\n" ) f.write("static const ephemeris_body_t ephemeris_bodies[] = {\n") for name, body in bodies.items(): f.write( f' {{ {name}_ephemeris, ' f'{len(body["segments"])}, ' f'{body["segment_days"]}, ' f'{body["degree"]} }}, // {name}\n' ) f.write("};\n") print(f"saved {filename}") return bodies if CREATE_EPHEMERIS: master_ephemeris_json() master_ephemeris_c()