32:Tba3,

During the late 16^th century, spanning over two decades, Tycho Brahe meticulously observed Mars (♂) and other planets. His precise measurements, unparalleled in his time, played a pivotal role in Johannes Kepler's establishment of the heliocentric and elliptical nature of planetary motion.

To accurately predict the orbit of a planet, a set of six parameters, known as osculating orbital elements, are required. These parameters are readily accessible through the NASA Horizons system, which extensively relies on data from the Jet Propulsion Laboratory's Development Ephemeris (JPL DE).

By combining Tycho Brahe's observations of Mars with the osculating orbital elements provided by JPL, we can effectively reconstruct the orbit of Mars as it unfolded at the close of the 16^th century.

This notebook aims to use nonlinear regression to derive the best-fitting osculating orbital elements for Mars during the period of 1582-1600, leveraging the observational data recorded by Tycho Brahe. The JavaScript code fits the astrophysical model to the observed data points and estimates the orbital elements using the method of least squares optimization. Additionally, bounded limits are imposed on the parameters to ensure an optimal fit while maintaining consistency with the underlying celestial mechanics.

We utilize the JPL DE440 parameters estimated for the year 1600 as an initial approximation. Subsequently, nonlinear regression is conducted employing least squares optimization with bounded parameters to refine the estimation. The results are visualized in the following graph, juxtaposing the JPL DE440 trajectory with the optimized fitting curve derived in this notebook. The number of iterations can be adjusted; it is currently set to ${max_iterations} (note that larger values increase the computation time).

33:T7ab,{ let dotplot = Plot.plot({ height: 700, width: 1200, grid: true, marginLeft: 50, marginRight: 30, marginTop: 50, marginBottom: 80, insetTop: 0, insetBottom: 0, style: { fontSize: "13px", color: "#000" }, x: { label: " Year → ", tickPadding: 10, labelOffset: 60, tickFormat: "f", tickSize: 9, ticks: 10, domain: [1582, 1601], tickRotate: -35 }, y: { label: "↑ Declination, δ", tickPadding: 6, labelOffset: 30, tickFormat: x => x<=0 ?`${x}°`: `+${x}°`, ticks: 10, tickSize: 6, domain: [-30, 30] }, caption: html`
Figure: Tycho Brahe's observations of the planet Mars (declination angle, period 1582 – 1600).
⚪ Tycho Brahe's observation.
– – – Declination of Mars derived from JPL DE440 orbital elements for year ${reference_date}.
—— Optimized fitting curve (least squares regression employed in this notebook).
`, marks: [ Plot.frame({stroke: "black"}), Plot.ruleY([0]), inputs_graph_display.show_initial_curve? Plot.line(initial_curve, {x: "Date", y: "Declination", curve: "catmull-rom", stroke: "red", strokeWidth: 0.8, strokeOpacity: 0.8, strokeDasharray: "4,4"}):null, inputs_graph_display.show_best_fit_curve? Plot.line(best_fit_curve, {x: "Date", y: "Declination", curve: "catmull-rom", stroke: "darkblue", strokeWidth: 0.8, strokeOpacity: 0.8}):null, inputs_graph_display.show_data_points? Plot.dot(Mars_data_Tycho_Brahe, {x: "Date", y: "Declination", stroke:"#111", fill: "orange", strokeOpacity: 0.4, r: dot_size}):null, ] }) // Label in bold/large font let select_labels = d3.select(dotplot) .selectAll("svg > g > text") .style("font-weight", "bold") .style("font-size", "16px"); return dotplot; }34:T97d,

Useful links

Planets: Reconstructing 16^th Century Orbits Using Tycho Brahe's Observations
Reconstructing 16^th Century Planetary Orbits Using Tycho Brahe's Observations and Joint Optimization

for all the planets

Visualizing Tycho Brahe's astronomical observations of Mars → Explore this earlier notebook showcasing Brahe's observations retrieved from Tychonis Brahe Dani Opera Omnia, the complete collection of Brahe's observation data.
Nonlinear regression with bounds in JavaScript → Detailed explanation of the underlying fitting algorithm used in this notebook is available here.
Mars ephemerides computation from osculating orbital elements (2000-2050) → This notebook focuses on computing the ephemerides of Mars for the contemporary period. While no fitting with observational data is conducted here, we can place trust in the accuracy of the osculating orbital elements used. The success of numerous NASA missions to Mars, among other evidence, underscores the reliability of these parameters.
Secular perturbation theory → Delve into perturbation theory in celestial mechanics with this notebook, providing insights into why the orbits of the planets are not immutable ellipses.

35:T6e2,// Model to compute the declination of planet Mars. model = function(x, p){return x.map(function(x_i){ const [a_earth, e_earth, i_earth, n_earth, o_earth, w_earth, M_earth, a_planet, e_planet, i_planet, n_planet, o_planet, w_planet, M_planet, obliquity_ecliptic] = p; const obliquity_rad = obliquity_ecliptic*Math.PI/180; const p_earth = o_earth + w_earth; const L_earth = p_earth + M_earth; const p_planet = o_planet + w_planet; const L_planet = p_planet + M_planet; const D = (x_i - reference_date) * 365.25; const MA_earth = n_earth * D + L_earth - p_earth; const MA_planet = n_planet * D + L_planet - p_planet; const v_earth = true_anomaly(MA_earth, e_earth); const v_earth_rad = v_earth*Math.PI/180; const v_planet = true_anomaly(MA_planet, e_planet); const v_planet_rad = v_planet*Math.PI/180; const r_earth = a_earth * (1 - Math.pow(e_earth,2)) / (1 + e_earth * Math.cos(v_earth_rad)); const r_planet = a_planet * (1 - Math.pow(e_planet,2)) / (1 + e_planet * Math.cos(v_planet_rad)); const { X: X_earth, Y: Y_earth, Z: Z_earth } = heliocentric_coordinates(r_earth, v_earth, o_earth, p_earth, i_earth); const { X: X_planet, Y: Y_planet, Z: Z_planet} = heliocentric_coordinates(r_planet, v_planet, o_planet, p_planet, i_planet); const X_prime = X_planet - X_earth; const Y_prime = Y_planet - Y_earth; const Z_prime = Z_planet - Z_earth; const X_q = X_prime; const Y_q = Y_prime * Math.cos(obliquity_rad) - Z_prime * Math.sin(obliquity_rad); const Z_q = Y_prime * Math.sin(obliquity_rad) + Z_prime * Math.cos(obliquity_rad); const DEC = Math.atan(Z_q/Math.sqrt(Math.pow(X_q,2)+Math.pow(Y_q,2)))*180/Math.PI; return DEC; })}36:T948,// Modified fminsearch function with bounds: Please note that this implementation is not optimized. // Adapted from: https://github.com/jonasalmeida/fminsearch // See associated notebook: https://observablehq.com/@christophe-yamahata/nonlinear-regression-with-bounds-in-javascript fminsearch = function(fun, Parm0, x, y, Opt) { if (!Opt) { Opt = {} }; if (!Opt.maxIter) { Opt.maxIter = 1000 }; if (!Opt.step) { // initial step is 1/100 of initial value (remember not to use zero in Parm0) Opt.step = Parm0.map(function(p) { return p / 100 }); Opt.step = Opt.step.map(function(si) { if (si == 0) { return 1 } else { return si } }); // convert null steps into 1's }; if (typeof (Opt.display) == 'undefined') { Opt.display = true }; if (!Opt.objFun) { Opt.objFun = function(y, yp) { return y.map(function(yi, i) { return Math.pow((yi - yp[i]), 2) }).reduce(function(a, b) { return a + b }) } } //SSD if (!Opt.bounds) { Opt.bounds = Parm0.map(function() { return [-Infinity, Infinity] }) }; // Default no bounds var cloneVector = function(V) { return V.map(function(v) { return v }) }; var ya, y0, yb, fP0, fP1; var P0 = cloneVector(Parm0), P1 = cloneVector(Parm0); var n = P0.length; let step = cloneVector(Opt.step); // Ensure step is reinitialized var funParm = function(P) { return Opt.objFun(y, fun(x, P)) } //function (of Parameters) to minimize // Apply bounds var applyBounds = function(P, bounds) { return P.map(function(p, i) { return Math.max(bounds[i][0], Math.min(bounds[i][1], p)); }); } // Silly multi-univariate screening for (var i = 0; i < Opt.maxIter; i++) { for (var j = 0; j < n; j++) { // take a step for each parameter P1 = cloneVector(P0); P1[j] += step[j]; P1 = applyBounds(P1, Opt.bounds); // Apply bounds after updating the parameter if (funParm(P1) < funParm(P0)) { // if parm value going in the right direction step[j] = 1.2 * step[j]; // then go a little faster P0 = cloneVector(P1); } else { step[j] = -(0.5 * step[j]); // otherwise reverse and go slower } } if (Opt.display) { if (i > (Opt.maxIter - 10)) { console.log(i + 1, funParm(P0), P0) } } } return P0 }37:T12be,python_script_1 = ` """ This Python script performs nonlinear regression to obtain the best-fitting osculating orbital elements for Mars for the period 1582-1600, using Tycho Brahe's observational data. The script fits a model to the observed data points and optimizes the orbital elements using least squares optimization from the SciPy library. For further details on Tycho Brahe's observations of Mars, visit: https://observablehq.com/@christophe-yamahata/visualizing-tycho-brahe-s-astronomical-observations-mars To use this script, download the required CSV file "Mars_data_1582-1600.csv" from the above URL. This file contains the necessary input data for the regression: dates as decimal years and the declination of Mars in degrees. In this version, the initial values for the osculating orbital elements are stored in the array "p_init". These values were obtained using the Skyfield astronomy package, configured with the JPL DE440 ephemeris (calculated for January 1, 1600). Author: C. Yamahata Date: May 2024 """ import numpy as np import pandas as pd from scipy.optimize import least_squares def true_anomaly(M, e): """Calculate true anomaly from mean anomaly and eccentricity.""" M_rad = np.radians(M) return M + (180 / np.pi) * ((2 * e - (1 / 4) * e**3) * np.sin(M_rad) + (5 / 4) * e**2 * np.sin(2 * M_rad) + (13 / 12) * e**3 * np.sin(3 * M_rad)) def heliocentric_coordinates(r, v, o, p, i): """Calculate heliocentric coordinates.""" v_rad = np.radians(v + p - o) o_rad = np.radians(o) i_rad = np.radians(i) X = r * (np.cos(o_rad) * np.cos(v_rad) - np.sin(o_rad) * np.sin(v_rad) * np.cos(i_rad)) Y = r * (np.sin(o_rad) * np.cos(v_rad) + np.cos(o_rad) * np.sin(v_rad) * np.cos(i_rad)) Z = r * np.sin(v_rad) * np.sin(i_rad) return X, Y, Z def model(x, p): """Model to compute the declination of Mars.""" (a_earth, e_earth, i_earth, n_earth, o_earth, w_earth, M_earth, a_mars, e_mars, i_mars, n_mars, o_mars, w_mars, M_mars, obliquity_ecliptic) = p p_earth = o_earth + w_earth L_earth = p_earth + M_earth p_mars = o_mars + w_mars L_mars = p_mars + M_mars D = (x - reference_date) * 365.25 MA_earth = n_earth * D + L_earth - p_earth MA_mars = n_mars * D + L_mars - p_mars v_earth = true_anomaly(MA_earth, e_earth) v_mars = true_anomaly(MA_mars, e_mars) r_earth = a_earth * (1 - e_earth**2) / (1 + e_earth * np.cos(np.radians(v_earth))) r_mars = a_mars * (1 - e_mars**2) / (1 + e_mars * np.cos(np.radians(v_mars))) X_earth, Y_earth, Z_earth = heliocentric_coordinates(r_earth, v_earth, o_earth, p_earth, i_earth) X_mars, Y_mars, Z_mars = heliocentric_coordinates(r_mars, v_mars, o_mars, p_mars, i_mars) X_prime = X_mars - X_earth Y_prime = Y_mars - Y_earth Z_prime = Z_mars - Z_earth obliquity_rad = np.radians(obliquity_ecliptic) Y_q = Y_prime * np.cos(obliquity_rad) - Z_prime * np.sin(obliquity_rad) Z_q = Y_prime * np.sin(obliquity_rad) + Z_prime * np.cos(obliquity_rad) DEC = np.degrees(np.arctan2(Z_q, np.sqrt(X_prime**2 + Y_q**2))) return DEC def objective_function(p, x, y): """Objective function for least squares optimization.""" return model(x, p) - y # Given data df = pd.read_csv("Mars_data_1582-1600.csv", delimiter=',', skiprows=1) x = df.iloc[:, 0].values y = df.iloc[:, 1].values reference_date = 1600 # Initial guess p_init = [1.000001, 0.016853, 0.052404, 0.985607, 355.777916, 105.93734, 3.725488, 1.523777, 0.092999, 1.88184, 0.523989, 50.727771, 283.539037, 141.282667, 23.491553] # Define lower and upper bounds for each parameter lower_bounds = [0.998 * p_init[0], 0.999 * p_init[1], 0.999 * p_init[2], 0.999 * p_init[3], 0, 0, 0, 0.998 * p_init[7], 0.999 * p_init[8], 0.999 * p_init[9], 0.999 * p_init[10], 0, 0, 0, min(23.4, p_init[14])] upper_bounds = [1.002 * p_init[0], 1.001 * p_init[1], 1.001 * p_init[2], 1.001 * p_init[3], 360, 360, 360, 1.002 * p_init[7], 1.001 * p_init[8], 1.001 * p_init[9], 1.001 * p_init[10], 360, 360, 360, max(23.5, p_init[14])] bounds = (lower_bounds, upper_bounds) # Optimize parameters using least squares with bounds result = least_squares(objective_function, p_init, bounds=bounds, args=(x, y)) # Variable names variable_names = ("a_earth", "e_earth", "i_earth", "n_earth", "o_earth", "w_earth", "M_earth", "a_mars", "e_mars", "i_mars", "n_mars", "o_mars", "w_mars", "M_mars", "obliquity_ecliptic") # Print best-fit parameters print("\\nBest-fit parameters:") for variable, value in zip(variable_names, result.x): print(f"{variable:<20} \\t {value:.6f}") `38:T1bff,python_script_2 = ` """ This Python script performs nonlinear regression to obtain the best-fitting osculating orbital elements for Mars for the period 1582-1600, using Tycho Brahe's observational data. The script fits a model to the observed data points and optimizes the orbital elements using least squares optimization from the SciPy library. For further details on Tycho Brahe's observations of Mars, visit: https://observablehq.com/@christophe-yamahata/visualizing-tycho-brahe-s-astronomical-observations-mars To use this script, download the required CSV file "Mars_data_1582-1600.csv" from the above URL. This file contains the necessary input data for the regression: dates as decimal years and the declination of Mars in degrees. In this version, the osculating orbital elements used for parameter initialization are directly loaded via the Skyfield astronomy package. Author: C. Yamahata Date: May 2024 """ import numpy as np import pandas as pd from scipy.optimize import least_squares from skyfield.api import load from skyfield.elementslib import osculating_elements_of from skyfield.data.spice import inertial_frames from datetime import datetime as dt def get_planet_elements(planets, planet, utc_time): """Get the osculating orbital elements of the specified planet.""" ecliptic = inertial_frames['ECLIPJ2000'] position = (planets[planet] - planets['sun']).at(utc_time) elements = osculating_elements_of(position, ecliptic) return elements def calculate_earth_obliquity(planets, utc_time): """Calculate the obliquity of the Earth.""" earth_position = (planets['earth barycenter'] - planets['sun']).at(utc_time) earth_elements = osculating_elements_of(earth_position) return earth_elements.inclination.degrees def true_anomaly(M, e): """Calculate true anomaly from mean anomaly and eccentricity.""" M_rad = np.radians(M) return M + (180 / np.pi) * ((2 * e - (1 / 4) * e**3) * np.sin(M_rad) + (5 / 4) * e**2 * np.sin(2 * M_rad) + (13 / 12) * e**3 * np.sin(3 * M_rad)) def heliocentric_coordinates(r, v, o, p, i): """Calculate heliocentric coordinates.""" v_rad = np.radians(v + p - o) o_rad = np.radians(o) i_rad = np.radians(i) X = r * (np.cos(o_rad) * np.cos(v_rad) - np.sin(o_rad) * np.sin(v_rad) * np.cos(i_rad)) Y = r * (np.sin(o_rad) * np.cos(v_rad) + np.cos(o_rad) * np.sin(v_rad) * np.cos(i_rad)) Z = r * np.sin(v_rad) * np.sin(i_rad) return X, Y, Z def model(x, p): """Model to compute the declination of Mars.""" (a_earth, e_earth, i_earth, n_earth, o_earth, w_earth, M_earth, a_mars, e_mars, i_mars, n_mars, o_mars, w_mars, M_mars, obliquity_ecliptic) = p p_earth = o_earth + w_earth L_earth = p_earth + M_earth p_mars = o_mars + w_mars L_mars = p_mars + M_mars D = (x - ref_date_decimal) * 365.25 MA_earth = n_earth * D + L_earth - p_earth MA_mars = n_mars * D + L_mars - p_mars v_earth = true_anomaly(MA_earth, e_earth) v_mars = true_anomaly(MA_mars, e_mars) r_earth = a_earth * (1 - e_earth**2) / (1 + e_earth * np.cos(np.radians(v_earth))) r_mars = a_mars * (1 - e_mars**2) / (1 + e_mars * np.cos(np.radians(v_mars))) X_earth, Y_earth, Z_earth = heliocentric_coordinates(r_earth, v_earth, o_earth, p_earth, i_earth) X_mars, Y_mars, Z_mars = heliocentric_coordinates(r_mars, v_mars, o_mars, p_mars, i_mars) X_prime = X_mars - X_earth Y_prime = Y_mars - Y_earth Z_prime = Z_mars - Z_earth obliquity_rad = np.radians(obliquity_ecliptic) Y_q = Y_prime * np.cos(obliquity_rad) - Z_prime * np.sin(obliquity_rad) Z_q = Y_prime * np.sin(obliquity_rad) + Z_prime * np.cos(obliquity_rad) DEC = np.degrees(np.arctan2(Z_q, np.sqrt(X_prime**2 + Y_q**2))) return DEC def objective_function(p, x, y): """Objective function for least squares optimization.""" return model(x, p) - y def to_year_fraction(date): """Convert a date to a decimal year.""" def days_in_year(year): return 365 if year % 4 != 0 or (year % 100 == 0 and year % 400 != 0) else 366 year = date.year start_of_this_year = dt(year=year, month=1, day=1) start_of_next_year = dt(year=year + 1, month=1, day=1) year_elapsed = (date - start_of_this_year).days year_duration = days_in_year(year) fraction = year_elapsed / year_duration return year + fraction # Reference date: January 1, 1600 ref_year = 1600 ref_month = 1 ref_day = 1 ref_date = dt(ref_year, ref_month, ref_day, 0, 0) ref_date_decimal = to_year_fraction(ref_date) timescale = load.timescale() ref_date_utc_time = timescale.utc(ref_year, ref_month, ref_day, 0, 0) # Load the JPL ephemeris DE440 planets = load('de440.bsp') # Get the osculating orbital elements for Earth and Mars (Used as initial guess) planet_names = ['earth barycenter', 'mars barycenter'] p_init = [] for planet_name in planet_names: elements = get_planet_elements(planets, planet_name, ref_date_utc_time) a = elements.semi_major_axis.au e = elements.eccentricity i = elements.inclination.degrees n = elements.mean_motion_per_day.degrees o = elements.longitude_of_ascending_node.degrees w = elements.argument_of_periapsis.degrees M = elements.mean_anomaly.degrees p_init.extend([round(a, 6), round(e, 6), round(i, 6), round(n, 6), round(o, 6), round(w, 6), round(M, 6)]) # Calculate and append Earth's obliquity to p_init earth_obliquity = calculate_earth_obliquity(planets, ref_date_utc_time) p_init.append(round(earth_obliquity, 6)) # Load Tycho Brahe's observation data from CSV df = pd.read_csv("Mars_data_1582-1600.csv", delimiter=',', skiprows=1) x = df.iloc[:, 0].values y = df.iloc[:, 1].values # Define lower and upper bounds for each parameter lower_bounds = [0.998 * p_init[0], 0.999 * p_init[1], 0.999 * p_init[2], 0.999 * p_init[3], 0, 0, 0, 0.998 * p_init[7], 0.999 * p_init[8], 0.999 * p_init[9], 0.999 * p_init[10], 0, 0, 0, min(23.4, p_init[14])] upper_bounds = [1.002 * p_init[0], 1.001 * p_init[1], 1.001 * p_init[2], 1.001 * p_init[3], 360, 360, 360, 1.002 * p_init[7], 1.001 * p_init[8], 1.001 * p_init[9], 1.001 * p_init[10], 360, 360, 360, max(23.5, p_init[14])] bounds = (lower_bounds, upper_bounds) # Optimize parameters using least squares with bounds result = least_squares(objective_function, p_init, bounds=bounds, args=(x, y)) # Variable names variable_names = ("a_earth", "e_earth", "i_earth", "n_earth", "o_earth", "w_earth", "M_earth", "a_mars", "e_mars", "i_mars", "n_mars", "o_mars", "w_mars", "M_mars", "obliquity_ecliptic") # Print best-fit parameters print("\\nVariable\\t\\tInitial guess\\tBest-fit") print('----------------------------------------------------') for variable, initial_guess, value in zip(variable_names, p_init, result.x): print(f"{variable:<20}\\t{initial_guess:.6f}\\t{value:.6f}") print(f'\\nReference date: {ref_date_utc_time.utc_datetime().strftime("%Y-%m-%d")}') print(f'TDB Julian date: {ref_date_utc_time.tdb:.2f}\\n') `39:[{"id":"1ae551704850c144a51fd14603f909f8b6b8b3026b15aaad01a735a81bfab4da6bb33218c13aca8153b002a0c884ecff1512ba2efa00587412b8a3b3d9e0964e","href":"https://static.observableusercontent.com/files/1ae551704850c144a51fd14603f909f8b6b8b3026b15aaad01a735a81bfab4da6bb33218c13aca8153b002a0c884ecff1512ba2efa00587412b8a3b3d9e0964e","name":"JPL_DE440_vs_least_squares_regression.mp4","type":"video/mp4","lastModified":1716985563387.866,"size":1322925},{"id":"e712dbd35aadce443d34977faa9f8488eea80b38e4b7e7912d9ab39b5f6bff9e288cd917e571700a857c7d317272d8f75724f31fffe839f3b096d018e303623a","href":"https://static.observableusercontent.com/files/e712dbd35aadce443d34977faa9f8488eea80b38e4b7e7912d9ab39b5f6bff9e288cd917e571700a857c7d317272d8f75724f31fffe839f3b096d018e303623a","name":"Mars_data_1582-1600.csv","type":"text/csv","lastModified":1716814041611.692,"size":17418}] 30:{"title":"Tycho Brahe's observations of Mars: osculating orbital elements for the period 1582-1600","cells":[{"id":0,"value":"# Tycho Brahe's observations of Mars: osculating orbital elements for the period 1582-1600","pinned":false,"hidden":false,"mode":"md"},{"id":521,"value":"

By combining Tycho Brahe's observations of Mars with the osculating orbital elements, we can effectively reconstruct the orbit of Mars as it unfolded at the close of the 16^th century.\n

","pinned":false,"hidden":false,"mode":"html"},{"id":372,"value":"$32","pinned":false,"hidden":false,"mode":"html"},{"id":256,"value":"viewof inputs_graph_display = columns({\n show_data_points : Inputs.toggle({label: \"Tycho Brahe's observations\", value: true}),\n show_initial_curve : Inputs.toggle({label: \"Declination of Mars derived from JPL DE440\", value: true}),\n show_best_fit_curve : Inputs.toggle({label: \"Optimized fitting curve (least squares regression)\", value: true}), \n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":547,"value":"display_dot_properties = {\n const html = inputs_graph_display.show_data_points? htl.html`${viewof dot_size}`: htl.html`

`; \n yield html\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":617,"value":"// As shown in the video below, the optimized fitting curve is achieved after about 200 iterations.\n// [see FileAttachment(\"JPL_DE440_vs_least_squares_regression.mp4\")]\n// Here, we randomly choose a number of iterations from the set [15, 20, 25, 30, 35, 40, 45, 50] \n// using the formula: 15 + 5 * Math.floor(Math.random() * 7)\nviewof iterations = Inputs.range([5, 200], {value: 15 + 5 * Math.floor(Math.random() * 7), step: 5, label: html`Num. of iterations
(least squares regression)

`})","pinned":false,"hidden":false,"mode":"ojs"},{"id":142,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":473,"value":"
\n

Python scripts

\nThe following Python scripts conduct nonlinear regression to derive the best-fitting osculating orbital elements for Mars \nduring the period 1582-1600, utilizing Tycho Brahe's observational data. Their implementation mirrors the JavaScript code used in this notebook.\n

\nThe scripts fit a model to the observed data points and optimize the orbital elements \nusing least squares optimization from the SciPy library.\n

\nTo use these script, download the required CSV file \"Mars_data_1582-1600.csv\". This file contains the necessary input data for the regression: \ndates as decimal years and the declination of Mars in degrees.\n

","pinned":false,"hidden":false,"mode":"html"},{"id":508,"value":"

\nScript #1 — In this shorter script, the osculating orbital elements' initial values are stored in the array \"p_init\". These values were acquired using the Skyfield astronomy package, configured with the JPL DE440 ephemeris (calculated for January 1, 1600).\n

","pinned":false,"hidden":false,"mode":"html"},{"id":469,"value":"Inputs.button(\"Copy script #1 to clipboard\", {value: null, reduce: () => navigator.clipboard.writeText(python_script_1)})","pinned":false,"hidden":false,"mode":"ojs"},{"id":510,"value":"

\nScript #2 — In this comprehensive version, the osculating orbital elements used for parameter initialization \nare directly loaded via the Skyfield astronomy package.\n

","pinned":false,"hidden":false,"mode":"html"},{"id":505,"value":"Inputs.button(\"Copy script #2 to clipboard\", {value: null, reduce: () => navigator.clipboard.writeText(python_script_2)})","pinned":false,"hidden":false,"mode":"ojs"},{"id":441,"value":"$34","pinned":false,"hidden":false,"mode":"html"},{"id":728,"value":"

Video: These screen captures illustrate Mars' orbit as measured by Tycho Brahe (period 1582-1600), compared with the orbits calculated using JPL DE440 and least squares regression over various iteration counts.

","pinned":false,"hidden":false,"mode":"html"},{"id":13,"value":"

Appendix

","pinned":false,"hidden":false,"mode":"html"},{"id":26,"value":"// Calculate true anomaly from mean anomaly and eccentricity.\n// tex`v = M + \\frac{180}{\\pi} \\times [ (2e - \\frac{1}{4}e^3) \\times \\sin(M) \n// + \\frac{5}{4}e^2 \\times \\sin(2M) + \\frac{13}{12} e^3 \\times \\sin(3M) ]`\ntrue_anomaly = function(M, e){\n var M_rad = M*Math.PI/180;\n return M + (180 / Math.PI) * ((2 * e - Math.pow(e,3)/4) * Math.sin(M_rad)\n + (5 / 4) * Math.pow(e,2) * Math.sin(2 * M_rad)\n + (13 / 12) * Math.pow(e,3) * Math.sin(3 * M_rad))\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":100,"value":"// Calculate heliocentric coordinates.\nheliocentric_coordinates = function(r, v, o, p, i){\n const v_rad = (v + p - o)*Math.PI/180;\n const o_rad = o*Math.PI/180;\n const i_rad = i*Math.PI/180;\n \n const X = r * (Math.cos(o_rad) * Math.cos(v_rad) - Math.sin(o_rad) * Math.sin(v_rad) * Math.cos(i_rad));\n const Y = r * (Math.sin(o_rad) * Math.cos(v_rad) + Math.cos(o_rad) * Math.sin(v_rad) * Math.cos(i_rad));\n const Z = r * Math.sin(v_rad) * Math.sin(i_rad);\n \n return {X, Y, Z};\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":81,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":195,"value":"create_fitting_curve = function(parameters){\n const data = [];\n const n = 1000;\n const year_min = 1582.5;\n const year_max = 1600.5;\n\n const years = d3.range(year_min, year_max, (year_max - year_min) / n).map(d => d.toFixed(4));\n const DEC = model(years,parameters);\n\n for (let i = 0; i < years.length; i++) {\n data.push({\n Date: parseFloat(years[i]),\n Declination: DEC[i]\n });\n }\n return data;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":76,"value":"// The initial values for the osculating orbital elements (stored in `p_init`)\n// were obtained for the reference date: January 1, 1600 (Julian Date, JD = 2305447.5)\nreference_date = 1600","pinned":false,"hidden":false,"mode":"ojs"},{"id":148,"value":"initial_curve = create_fitting_curve(p_init)","pinned":false,"hidden":false,"mode":"ojs"},{"id":199,"value":"best_fit_curve = create_fitting_curve(p_fitting_with_bounds)","pinned":false,"hidden":false,"mode":"ojs"},{"id":44,"value":"// Initial guess\n// p_init = [a_eath, e_earth, i_earth, n_earth, o_earth, w_earth, M_earth,\n// a_mars, e_mars, i_mars, n_mars, o_mars, w_mars, M_mars, \n// obliquity_ecliptic]\n// The initial values for the osculating orbital elements are stored in the array `p_init`. \n// These values were obtained using the Python Skyfield astronomy package, configured with the JPL DE440 ephemeris.\n// Reference date: January 1, 1600 (Julian Date, JD = 2305447.5)\np_init = [1.000001, 0.016853, 0.052404, 0.985607, 355.777916, 105.93734, 3.725488, \n 1.523777, 0.092999, 1.88184, 0.523989, 50.727771, 283.539037, 141.282667, \n 23.491553]","pinned":false,"hidden":false,"mode":"ojs"},{"id":131,"value":"p_fitting_with_bounds = fminsearch(model,p_init,x,y,Opt);","pinned":false,"hidden":false,"mode":"ojs"},{"id":49,"value":"// Lower bounds for each fitting parameter\nlower_bounds = [0.998 * p_init[0], 0.999 * p_init[1], 0.999 * p_init[2], 0.999 * p_init[3], 0, 0, 0,\n 0.998 * p_init[7], 0.999 * p_init[8], 0.999 * p_init[9], 0.999 * p_init[10], 0, 0, 0,\n Math.min(23.4, p_init[14])]","pinned":false,"hidden":false,"mode":"ojs"},{"id":53,"value":"// Upper bounds for each fitting parameter\nupper_bounds = [1.002 * p_init[0], 1.001 * p_init[1], 1.001 * p_init[2], 1.001 * p_init[3], 360, 360, 360,\n 1.002 * p_init[7], 1.001 * p_init[8], 1.001 * p_init[9], 1.001 * p_init[10], 360, 360, 360,\n Math.max(23.5, p_init[14])]","pinned":false,"hidden":false,"mode":"ojs"},{"id":273,"value":"// Maximum iterations for the nonlinear least squares optimization\nmax_iterations = iterations","pinned":false,"hidden":false,"mode":"ojs"},{"id":63,"value":"// Define the bounding parameters for nonlinear regression calculation using bounded parameters.\nOpt = {\n const bounds = lower_bounds.map((lower, index) => [lower, upper_bounds[index]]);\n return {\n maxIter: max_iterations, // Optional: Specify the maximum number of iterations\n bounds: bounds, // Bounds for the parameters\n }\n};","pinned":false,"hidden":false,"mode":"ojs"},{"id":17,"value":"// data: x = dates as decimal years\nx = {\n let data = [];\n for (let i = 0; i < Mars_data_Tycho_Brahe.length; i++) {\n data.push(Mars_data_Tycho_Brahe[i].Date);\n }\n return data;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":20,"value":"// data: y = declination of Mars in degrees\ny = {\n let data = [];\n for (let i = 0; i < Mars_data_Tycho_Brahe.length; i++) {\n data.push(Mars_data_Tycho_Brahe[i].Declination);\n }\n return data;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":11,"value":"// Tycho Brahe's observation of Mars (Declination) for the period 1582-1600\nMars_data_Tycho_Brahe = FileAttachment(\"Mars_data_1582-1600.csv\").csv({typed: true})","pinned":false,"hidden":false,"mode":"ojs"},{"id":709,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":258,"value":"// https://observablehq.com/@bcardiff/observable-columns\ncolumns = (args) => {\n const form = html``\n form.value = {}\n\n let cols = 0\n for (const key in args) {\n form.appendChild(args[key])\n cols++\n }\n \n form.style = `width:200px; display: grid; grid-gap: 10px 100px; grid-template-columns: repeat(${cols}, auto); grid-auto-flow: row;`\n \n form.oninput = () => {\n form.value = Object.keys(args).reduce((result, key) => {\n result[key] = args[key].value\n return result\n }, Array.isArray(args) ? [] : {})\n }\n form.oninput()\n \n return form\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":565,"value":"viewof dot_size = Inputs.range([0.4, 4], {value: 2.4, step: 0.2, label: \"Dot size:\"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":464,"value":"$37","pinned":false,"hidden":false,"mode":"ojs"},{"id":513,"value":"$38","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W39","stdlib":"1","titling":"auto"}