2a:T482,

This notebook builds on a previous one that compared Tycho Brahe's observations of Mercury, Venus, Mars, Jupiter, and Saturn with their estimated osculating orbital elements.

In this continuation, we employ joint optimization to determine the optimal osculating orbital elements for all the planets by minimizing the sum of squared differences between Tycho Brahe's observations and the predicted values. This approach includes bounded constraints on the osculating orbital elements to ensure the results align with underlying celestial mechanics. Through this process, the goal is to accurately reconstruct the planetary orbits as they unfolded at the close of the 16^th century.

2b:T86c,{ let dotplot = Plot.plot({ height: 600, width: 1100, grid: true, marginLeft: 50, marginRight: 30, marginTop: 45, marginBottom: 70, insetTop: 0, insetBottom: 0, style: { fontSize: "13px", color: "#000" }, x: { label: " Year → ", tickPadding: 10, labelOffset: 65, tickFormat: "f", tickSize: 9, ticks: 10, domain: [xMin, xMax], }, y: { label: "↑ Declination of "+planet.name+", δ", tickPadding: 6, labelOffset: 40, tickFormat: "f", ticks: 10, tickSize: 6, domain: [yMin, yMax] }, caption: html`
Figure 1: Comparison of Tycho Brahe's observations of ${planet.name} (${planet.symbol}) with joint optimization results.
⚪ Tycho Brahe's observations of the declination
– – – Derived from JPL DE440 orbital elements for year 1600
—— Joint optimization
`, marks: [ Plot.frame({stroke: "black"}), Plot.ruleY([0]), inputs_graph_display.show_initial_curve? Plot.line(createCurve(initialParams,planet.dataset), {x: "x", y: "y", curve: "catmull-rom", stroke: strokeColor1, strokeWidth: 0.8, strokeOpacity: 0.8, strokeDasharray: "4,4"}):null, inputs_graph_display.show_best_fit_curve? Plot.line(createCurve(optimizedParams,planet.dataset), {x: "x", y: "y", curve: "catmull-rom", stroke: strokeColor2, strokeWidth: 0.8, strokeOpacity: 0.8}):null, inputs_graph_display.show_data_points? Plot.dot(dataPoints .filter(point => point.dataset === planet.dataset) .filter(point => point.x >= xMin && point.x <= xMax), // Filter by x range {x: "x", y: "y", stroke:"#111", fill: dotColor, strokeOpacity: 0.4, r: dotSize}):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; }2c:T5fb,{ let dotplot = Plot.plot({ height: 600, width: 1100, grid: true, marginLeft: 50, marginRight: 40, marginTop: 40, marginBottom: 70, insetTop: 0, insetBottom: 0, style: { fontSize: "16px", background: "#fff", color: "#000" }, x: { label: " ← Right ascension of "+planet.name+", α ", tickPadding: 6, labelOffset: 64, transform: f => f/6, domain: [360/6,0], tickFormat: x => `${Math.round(6*x.toFixed(3))}°`, ticks: 10, tickSize: 9, tickRotate: 0 }, y: { label: "↑ Declination of "+planet.name+", δ", labelOffset: 35, tickFormat: x => x<=0 ?`${x}°`: `+${x}°`, ticks: 10, tickSize: 6, domain: [-30, 30] }, caption: html`
Figure 2: Theoretical astronomical coordinates of ${planet.name} (${planet.symbol}) from ${xMin} to ${xMax}.
The coordinates are derived from osculating orbital elements obtained through joint optimization.`, marks: [ Plot.frame({stroke: "black"}), Plot.ruleY([0]), Plot.line(createCurve_DEC_RA(initialParams,planet.dataset) .map(d => d.x <= 15 || d.x >= 355 ? NaN : d), {x: "x", y: "y", curve: "catmull-rom", stroke: strokeColor2, strokeWidth: 0.8, strokeOpacity: 0.8}) ] }) // Label in bold let select_labels = d3.select(dotplot) .selectAll("svg > g > text") .style("font-weight", "bold"); return dotplot; }2d:Td66,// Model to compute the declination of a planet as observed from Earth. function model(params, x, dataset) { const [ a_earth, e_earth, i_earth, n_earth, o_earth, w_earth, M_earth, a_mercury, e_mercury, i_mercury, n_mercury, o_mercury, w_mercury, M_mercury, a_venus, e_venus, i_venus, n_venus, o_venus, w_venus, M_venus, a_mars, e_mars, i_mars, n_mars, o_mars, w_mars, M_mars, a_jupiter, e_jupiter, i_jupiter, n_jupiter, o_jupiter, w_jupiter, M_jupiter, a_saturn, e_saturn, i_saturn, n_saturn, o_saturn, w_saturn, M_saturn, obliquity_ecliptic ] = params; let a_planet, e_planet, i_planet, n_planet, o_planet, w_planet, M_planet; switch (dataset) { case 'y0': // Mercury a_planet = a_mercury; e_planet = e_mercury; i_planet = i_mercury; n_planet = n_mercury; o_planet = o_mercury; w_planet = w_mercury; M_planet = M_mercury; break; case 'y1': // Venus a_planet = a_venus; e_planet = e_venus; i_planet = i_venus; n_planet = n_venus; o_planet = o_venus; w_planet = w_venus; M_planet = M_venus; break; case 'y2': // Mars a_planet = a_mars; e_planet = e_mars; i_planet = i_mars; n_planet = n_mars; o_planet = o_mars; w_planet = w_mars; M_planet = M_mars; break; case 'y3': // Jupiter a_planet = a_jupiter; e_planet = e_jupiter; i_planet = i_jupiter; n_planet = n_jupiter; o_planet = o_jupiter; w_planet = w_jupiter; M_planet = M_jupiter; break; case 'y4': // Saturn a_planet = a_saturn; e_planet = e_saturn; i_planet = i_saturn; n_planet = n_saturn; o_planet = o_saturn; w_planet = w_saturn; M_planet = M_saturn; break; default: throw new Error('Unknown dataset...'); } 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 - 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 = trueAnomaly(MA_earth, e_earth); const v_earth_rad = v_earth * Math.PI / 180; const v_planet = trueAnomaly(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 } = heliocentricCoordinates(r_earth, v_earth, o_earth, p_earth, i_earth); const { X: X_planet, Y: Y_planet, Z: Z_planet } = heliocentricCoordinates(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; var alpha = Math.atan(Y_q/X_q)*180/Math.PI; alpha = X_q <0? alpha + 180: alpha; // If X_q is negative then add 180 degrees to alpha alpha = X_q >0 && Y_q <0 ? alpha + 360: alpha; // If X_q is positive and Yq is negative then add 360 degrees to alpha return [DEC,alpha]; } 2e:T964,// 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; // var step = Opt.step; 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 }2f:T59f,paramsPlanet = function(params,dataset){ const i = planets.find(planet => planet.dataset === dataset).index; const result = [ { symbol: planets.find(planet => planet.dataset === dataset).symbol, description: planets.find(planet => planet.dataset === dataset).name, unit: "", value: "", }, { symbol: "n", description: "Mean motion", unit: "°/day", value: params[3+i*7], }, { symbol: "a", description: "Semi-major axis", unit: "AU", value: params[0+i*7], }, { symbol: "e", description: "Eccentricity", unit: "", value: params[1+i*7], }, { symbol: "i", description: "Inclination", unit: "°", value: params[2+i*7], }, { symbol: "Ω", description: "Longitude of ascending node", unit: "°", value: params[4+i*7], }, { symbol: "ω", description: "Argument of periapsis", unit: "°", value: params[5+i*7], }, { symbol: "M", description: "Mean anomaly", unit: "°", value: params[6+i*7], }, ]; // Provide Obliquity of the ecliptic only for Earth if (i === 0) { result.push({ symbol: "ε", description: "Obliquity of the ecliptic (Earth's tilt)", unit: "°", value: params[35], }); } return result; }30:T575,displayOrbitalElements = function(params) { // Define the keys and corresponding formatting functions const keys = ['planet', 'n', 'a', 'e', 'i', 'Ω', 'ω', 'M', 'ε']; const formats = [ x => x, // for 'planet' x => x.toFixed(4), // for 'n' x => x.toFixed(5), // for 'a' x => x.toFixed(5), // for 'e' x => x.toFixed(4), // for 'i' x => x.toFixed(2), // for 'Ω' x => x.toFixed(2), // for 'ω' x => x.toFixed(2), // for 'M' x => x.toFixed(3) // for 'ε' ]; const format = Object.fromEntries(keys.map((key, index) => [key, formats[index]])); // Apply the formatting to the params values const formattedParams = params.map(param => { if (format[param.symbol]) { return { ...param, value: format[param.symbol](param.value) }; } return param; }); return Inputs.table(formattedParams, { columns: [ "symbol", "description", "value", "unit", ], header: { symbol: "Symbol", description: "Description", value: "Value", unit: "Unit", }, align: { symbol: "center", description: "left", value: "right", unit: "left", }, width: { symbol: 60, description: 180, value: 80, unit: 40, }, maxWidth: 500, maxHeight: 250, height: 250 }); } 31:T6b3,// https://observablehq.com/@mootari/range-slider#interval function interval(range = [], options = {}) { const [min = 0, max = 1] = range; const { step = .001, label = null, value = [min, max], format = ([start, end]) => `${start} … ${end}`, color, width, theme, } = options; const __ns__ = DOM.uid('scope').id; const css = ` #${__ns__} { font: 13px/1.2 var(--sans-serif); display: flex; align-items: baseline; flex-wrap: wrap; max-width: 100%; width: auto; } @media only screen and (min-width: 30em) { #${__ns__} { flex-wrap: nowrap; width: 360px; } } #${__ns__} .label { width: 120px; padding: 5px 0 4px 0; margin-right: 6.5px; flex-shrink: 0; } #${__ns__} .form { display: flex; width: 100%; } #${__ns__} .range { flex-shrink: 1; width: 100%; } #${__ns__} .range-slider { width: 100%; } `; const $range = rangeInput({min, max, value: [value[0], value[1]], step, color, width, theme}); const $output = html``; const $view = html`

${label == null ? '' : html`

${label}`}

${$range}

${$output}

${html`