Data to "check" Kepler’s first law

Question

I want to "check" Kepler’s first law by using real data of Mars. From the equation of the ellipse, I derived

$$\frac{1}{r}=\frac{a}{b^2}+\frac{a}{b^2}\cdot\epsilon\cdot\cos(\varphi),$$

where $a$ is the major semi-axis, $b$ is the minor semi-axis and $\epsilon$ is the eccentricity of the elliptic orbit. I'm looking for the following kind of data:

1. Mars' distance from Sun $r$
2. the angle $\varphi$ between Mars, Sun and the principal axis of the elliptic orbit.

Then, I want to check, whether $r$ and $\varphi$ fit the measured values of $a$, $b$ and $\epsilon$. If there is no such data (perpendicular view on Mars' orbital plane) available, how can I transform data given in other coordinate systems to the ones I need? On a NASA website (https://omniweb.gsfc.nasa.gov/coho/helios/heli.html) I found data in "Solar Ecliptic", "Heliographic" and "Heliographic Inertial" coordinates, but I don't know which come closest to my plan.

Update:

I tried it with uhoh's recommendations. Unfortunately I failed.

With the following python code, using the Horizons x, y, z data stored in an xlsx file,

from __future__ import division
import numpy as np
from statsmodels.regression.linear_model import OLS
from statsmodels.tools import add_constant
from statsmodels.tools.eval_measures import aicc
import pandas as pd
import matplotlib.pyplot as plt
horizons = pd.read_excel("horizons2.xlsx")
horizons = np.array(horizons)
horizonsxyz=horizons[:,2:5]
horizonsxyz=np.array(horizonsxyz, dtype=np.float64)
hx=horizonsxyz[:,0]
hy=horizonsxyz[:,1]
hz=horizonsxyz[:,2]

horizonsr=np.sqrt(hx**2+hy**2+hz**2)
horizonsr=horizonsr*6.68459*(10**(-9))

phi=np.arctan2(hy, hx) * 180 / np.pi
phi2=np.mod(phi+360, 360)
phia=np.mod(phi-286, 360)
phiganz=add_constant(phia)
horizonsdurchr=1/horizonsr




horizons_regr=OLS(horizonsdurchr, phiganz).fit()
print(horizons_regr.params)
print(horizons_regr.summary())
y_pred_horizons=np.dot(phiganz, horizons_regr.params)
print(horizons_regr.params)

I get a value of $7.1349\cdot10^{-1}$ for $\frac{a}{b^2}$. This is bad but at least in the right order of magnitude. However for $\frac{a}{b^2}\cdot\epsilon$ I get a really bad value of $-2.89228\cdot10^{-4}$. Deviding the two result yields an estimated excentricity of $0.00044$ which is really far away from the true $0.0934$.

I also tried another approach, using the heliographic data mentioned above. Here, I get closer, but only if I add 35 degrees to the angles, which doesnt make sense, since I should add 74 degrees or subtract 278 degrees, to get the angle relative to the perihelion.

Answer 1

Great project! and welcome to Stack Exchange. I'll post a short answer but I think someone can add a more detailed, thorough and insightful answer.

I think that website is not well suited, so I'll answer based on you switching to Horizons. If you like Python then it's more fun to use Skyfield.

If you want apply an equation based on a Kepler orbit model, you'll need to use data where the Sun stays in one place and Mars orbits around it. That would be Heliocentric with the Sun at (0, 0, 0).

That there are three zeros raises the issue of the number of dimensions; proper Kepler orbits are sort-of in 3D i.e. they have an orbital plane that can be tilted to a reference plane, but the orbits are planar. Two problems; your equation assumes 2D flat orbit because of the way $\varphi$ is defined. Ideally you'd like data in the plane of Mars' orbit and you may need to transform NASA/JPL Horizons data into Mars' orbital plane yourself because there are only two main "official" planes, no real planet remains perfectly in a plane.

So what you do depends on how far down the rabbit hole of pretending orbits are planes that you want to go.

Zeroth order approximation

Go to Horizons

Use this tutorial and set it up to match the following:

Current Settings
Ephemeris Type:      VECTORS
Target Body:         Mars [499]
Coordinate Origin:   Sun (body center) [500@10]
Time Span:           Start=2020-10-04, Stop=2020-10-05, Step=1 d
Table Settings:      quantities code=2; output units=KM-S; CSV format=YES
Display/Output:      default (formatted HTML)
             -- OR --
Display/Output:      download/save (plain text file)

Here's a sample line for Mars for today using the Sun as origin (I've truncated some decimal digits). You see right away that Mars is about 201 million km from the Sun, it is also about 4 million km below the J2000.0 ecliptic.

2459126.500, A.D. 2020-Oct-04 00:00:00.00,  2.036231544E+08,  5.355405115E+07, -3.872888712E+06...

From here you can approximate

$$r = \sqrt{x^2 + y^2 + z^2}$$

and

$$\varphi = \arctan2(y, x) - \text{286.502°}$$

Since you are going through all four quadrants it's better to use a computer's arctan2(y, x) or atan2(y, x) with two arguments, not $\arctan(y/x)$ which only works in two quadrants (i.e. 1/7 = -1/-7).

First order approximation

You see right away that Mars is about 201 million km from the Sun, it is also about 4 million km below the J2000.0 ecliptic.

If you want to correct for Mars' orbit's tilt with respect to the ecliptic, you can just find the best plane fit to one Martian year's of data and make your own Mars ecliptic.

But I recommend you do the zeroth order first and see how well or poorly it works, then you can decide if you want to tilt.

Answer 2

Kepler's first law is that a planet moves in an ellipse with the sun at one focus. Your equation is that of an ellipse about the focus, so, you have proven Kepler's first law. The $\varphi$ is what astronomers call true anomaly. To put your equation in the usual form, $a/b^2$ is $1/p$ so $$ r = \frac{p}{1+\epsilon cos\varphi}$$

With this equation, the ellipse can be plotted by choosing many values of the angle $\varphi$ and finding the corresponding r values, then plotting.

The p is called the parameter by astronomers and semi-latus rectum by mathematicians. As you can see, when $\varphi$ is 90 degrees, the value of r is p. Also, $p=a(1-\epsilon^2)$ which can be put in the equation above as an alternate form of the equation.

Keplers law gives no information about where the perihelion is in the orbit plane.