<p>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 <span class="math-container">$\varphi$</span> is what astronomers call true anomaly.  To put your equation in the usual form, <span class="math-container">$a/b^2$</span> is <span class="math-container">$1/p$</span> so
<span class="math-container">$$ r = \frac{p}{1+\epsilon cos\varphi}$$</span></p>
<p>With this equation, the ellipse can be plotted by choosing many values of the angle <span class="math-container">$\varphi$</span> and finding the corresponding r values, then plotting.</p>
<p>The p is called the parameter by astronomers and semi-latus rectum by mathematicians. As you can see, when <span class="math-container">$\varphi$</span> is 90 degrees, the value of r is p. Also, <span class="math-container">$p=a(1-\epsilon^2)$</span> which can be put in the equation above as an alternate form of the equation.</p>
<p>Keplers law gives no information about where the perihelion is in the orbit plane.</p>

CC BY-SA 4.0
2020-10-04T19:03:51.400Z
2020-10-04T19:03:51.400Z
W H G
<p>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.</p>
<p>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 <a href="https://rhodesmill.org/skyfield/" rel="noreferrer">Skyfield</a>.</p>
<p>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 <em>Heliocentric</em> with the Sun at (0, 0, 0).</p>
<p>That there are <em>three zeros</em> 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 <em>assumes</em> 2D flat orbit because of the way <span class="math-container">$\varphi$</span> 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 &quot;official&quot; planes, no real planet remains perfectly in a plane.</p>
<p>So what you do depends on how far down the rabbit hole of pretending orbits are planes that you want to go.</p>
<h2>Zeroth order approximation</h2>
<p>Go to <a href="https://ssd.jpl.nasa.gov/horizons.cgi" rel="noreferrer">Horizons</a></p>
<p>Use <a href="https://space.stackexchange.com/a/25357/12102">this tutorial</a> and set it up to match the following:</p>
<pre><code>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)
</code></pre>
<p><a href="https://i.sstatic.net/Q388H.png" rel="noreferrer"><img src="https://i.sstatic.net/Q388H.png" alt="JPL Horizons setup" /></a></p>
<p>Here's a sample line for Mars for today <strong>using the Sun as origin</strong> (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.</p>
<pre><code>2459126.500, A.D. 2020-Oct-04 00:00:00.00,  2.036231544E+08,  5.355405115E+07, -3.872888712E+06...
</code></pre>
<p>From here you can approximate</p>
<p><span class="math-container">$$r = \sqrt{x^2 + y^2 + z^2}$$</span></p>
<p>and</p>
<p><span class="math-container">$$\varphi = \arctan2(y, x) - \text{286.502°}$$</span></p>
<p>Since you are going through all four quadrants it's better to use a computer's <code>arctan2(y, x)</code> or <code>atan2(y, x)</code> with two arguments, not <span class="math-container">$\arctan(y/x)$</span> which only works in two quadrants (i.e. 1/7 = -1/-7).</p>
<h2>First order approximation</h2>
<p>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.</p>
<p>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.</p>
<p>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.</p>

CC BY-SA 4.0
2020-10-04T08:54:04.243Z
That's really helpful. However, I dont understand where the 286.502 deg comes from.
Joe_base
CC BY-SA 4.0
2020-10-04T09:01:50.103Z
@Joe_base Wikipedia's [Mars](https://en.wikipedia.org/wiki/Mars) gives that angle as the [argument of perihelion](https://en.wikipedia.org/wiki/Argument_of_periapsis) The $x, y, z$ data in Horizons is in a standard coordinate system called J2000.0 (the details are indeed quite detailed but let's just say on a certain day in the year 2000 the ecliptic and Earth's axis and equator were "frozen" and that's *roughly speaking* the coordinate system everyone uses. The +x direction is like 0 degrees, the +y direction is like 90 degrees.
uhoh
CC BY-SA 4.0
2020-10-04T09:02:31.877Z
So if you get an angle based on $x, y$ you need to subtract 286.502 to get the angle $\varphi$  relative to the perihelion.
uhoh
CC BY-SA 4.0
2020-10-13T16:59:54.670Z
@Joe_base okay I see. Usually it's bad to change a question once answers are posted, but this is a special case and you're basically refining the same question. I don't have much time in the next few days to look at this so instead I'll add a bounty to attract attention and maybe someone else can chime in. Thanks!
uhoh
CC BY-SA 4.0
2020-10-19T16:49:50.083Z
https://ssd.jpl.nasa.gov/txt/aprx_pos_planets.pdf, maybe you need complement with this docuement
Adrian R
CC BY-SA 4.0
2020-10-19T17:03:49.767Z
@AdrianR That is a great idea, nice! I'd seen those years ago and forgotten, but that is a much easier solution. I don't have time to do go further on this, so instead I'd added the bounty which expires 2020-10-20 17:00 UTC, plus has as 24 hour grace period. If you write up an short answer based on that it would be helpful to the OP and competitive towards the bounty, maybe I'll even add a second one.
uhoh
CC BY-SA 4.0
2020-11-09T06:16:23.637Z
Just to keep you updated: I finally found out, that I had to take the longitude of the ascending node into account. Thus, in the J2000 coordinate system (Argument of Perihelion + Longitude of Ascending Node) has to be deducted to get the right values. Doind that it worked perfectly fine. By doing that, both the heliographic coordinates I started with as well as uhoh's x-y-z account work. I did not further pursue the Keplerian elements.
Joe_base
CC BY-SA 4.0
2020-10-04T00:51:20.237Z
2020-10-04T00:51:20.237Z
uhoh
<p>I want to &quot;check&quot; Kepler’s first law by using real data of Mars. From the equation of the ellipse, I derived</p>
<p><span class="math-container">$$\frac{1}{r}=\frac{a}{b^2}+\frac{a}{b^2}\cdot\epsilon\cdot\cos(\varphi),$$</span></p>
<p>where <span class="math-container">$a$</span> is the major semi-axis, <span class="math-container">$b$</span> is the minor semi-axis and  <span class="math-container">$\epsilon$</span> is the eccentricity of the elliptic orbit. I'm looking for the following kind of data:</p>
<ol>
<li>Mars' distance from Sun <span class="math-container">$r$</span></li>
<li>the angle <span class="math-container">$\varphi$</span> between Mars, Sun and the principal axis of the elliptic orbit.</li>
</ol>
<p>Then, I want to check, whether <span class="math-container">$r$</span> and <span class="math-container">$\varphi$</span> fit the measured values of <span class="math-container">$a$</span>, <span class="math-container">$b$</span> and <span class="math-container">$\epsilon$</span>.
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 (<a href="https://omniweb.gsfc.nasa.gov/coho/helios/heli.html" rel="nofollow noreferrer">https://omniweb.gsfc.nasa.gov/coho/helios/heli.html</a>) I found data in &quot;Solar Ecliptic&quot;, &quot;Heliographic&quot; and &quot;Heliographic Inertial&quot; coordinates, but I don't know which come closest to my plan.</p>
<p>Update:</p>
<p>I tried it with uhoh's recommendations. Unfortunately I failed.</p>
<p>With the following python code, using the Horizons x, y, z data stored in an xlsx file,</p>
<pre><code>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(&quot;horizons2.xlsx&quot;)
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)
</code></pre>
<p>I get a value of <span class="math-container">$7.1349\cdot10^{-1}$</span> for <span class="math-container">$\frac{a}{b^2}$</span>. This is bad but at least in the right order of magnitude. However for <span class="math-container">$\frac{a}{b^2}\cdot\epsilon$</span> I get a really bad value of <span class="math-container">$-2.89228\cdot10^{-4}$</span>. Deviding the two result yields an estimated excentricity of <span class="math-container">$0.00044$</span> which is really far away from the true <span class="math-container">$0.0934$</span>.</p>
<p>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.</p>

<coordinate><mathematics><kepler><raw-data><space-geometry>
Data to "check" Kepler’s first law
Joe_base
# Data to "check" Kepler’s first law

(Asked by: username_0 on 2020-10-03)

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:

Mars' distance from Sun $r$
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 username_2'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.

**Tags:** coordinate, mathematics, kepler, raw-data, space-geometry


---

## 2 answers

### Answer by username_2 on 2020-10-04. Score: 4 | Accepted answer

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.

> **Comments:**
> **username_0**: That's really helpful. However, I dont understand where the 286.502 deg comes from.  
> **username_2**: @username_0 Wikipedia's [Mars](https://en.wikipedia.org/wiki/Mars) gives that angle as the [argument of perihelion](https://en.wikipedia.org/wiki/Argument_of_periapsis) The $x, y, z$ data in Horizons is in a standard coordinate system called J2000.0 (the details are indeed quite detailed but let's just say on a certain day in the year 2000 the ecliptic and Earth's axis and equator were "frozen" and that's *roughly speaking* the coordinate system everyone uses. The +x direction is like 0 degrees, the +y direction is like 90 degrees.  
> **username_2**: So if you get an angle based on $x, y$ you need to subtract 286.502 to get the angle $\varphi$  relative to the perihelion.  
> **username_2**: @username_0 okay I see. Usually it's bad to change a question once answers are posted, but this is a special case and you're basically refining the same question. I don't have much time in the next few days to look at this so instead I'll add a bounty to attract attention and maybe someone else can chime in. Thanks!  
> **username_3**: https://ssd.jpl.nasa.gov/txt/aprx_pos_planets.pdf, maybe you need complement with this docuement  
> **username_2**: @AdrianR That is a great idea, nice! I'd seen those years ago and forgotten, but that is a much easier solution. I don't have time to do go further on this, so instead I'd added the bounty which expires 2020-10-20 17:00 UTC, plus has as 24 hour grace period. If you write up an short answer based on that it would be helpful to the OP and competitive towards the bounty, maybe I'll even add a second one.  
> **username_0**: Just to keep you updated: I finally found out, that I had to take the longitude of the ascending node into account. Thus, in the J2000 coordinate system (Argument of Perihelion + Longitude of Ascending Node) has to be deducted to get the right values. Doind that it worked perfectly fine. By doing that, both the heliographic coordinates I started with as well as username_2's x-y-z account work. I did not further pursue the Keplerian elements.  

---

### Answer by username_1 on 2020-10-04. Score: 3

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.


---

