Function for current distance between in light-seconds between Mars and Earth

Question

I'm looking for a function that given a time (for example Unix time) gives back the current distance between Mars and Earth in light-seconds.

It can be a mathematical function or programming function (preferably in Python or pseudocode), but must not rely on an Internet connection.

This question is similar to a question I asked erlier, but again I need to get the distance without an Internet connection.

Thanks in advance!

Sidenote: I'm planning on learning orbital mechanics, but I need this function earlier. :)

Answer 1

There is no simple function that is accurate. Both Earth and Mars are orbiting in ellipses (Kepler's first law). They are travelling at varying speeds (Kepler's second). To find the current distance between them you have to calculate their positions in three dimensions.

In other words you have to do orbital kinematics. How do I calculate the positions of objects in orbit?

If accuracy is not important. You can approximate the orbits with circles in the same plane. The Earths position is then $$(\cos(2\pi t + \alpha), \sin(2\pi t +\alpha))$$ and Mars is $$(1.5\cos(2\pi t/1.882 +\beta), 1.5\sin(2\pi t/1.882 + \beta))$$

The distance between them is calculated by pythagoras

The distances are in AU and the time is in years (converting to lightseconds is then easy). $\alpha$ and $\beta$ are constants chosen to fit the positions of the planets at time 0. The values 1.5 and 1.882 are chosen because Mars's orbit is about 1.5 times bigger than Earth, and a Mars year is 1.882 times longer than Earth.

The circular approximation however is poor for Mars, and this ignores the fact that Mars isn't in exactly the same orbital plane as the earth, so the values will only be rather approximate.

Answer 2

note: This is an "extended comment" type answer, addressing two potential solutions that may not be exactly what you need, but are close.

A thing which gives you positions of things in the solar system given a specified time is called an Ephemeris.

Originally they were tables of predictions (past and future) based on calculations. By the time of Newton, they became quite good, since Newton developed a method to iteratively solve Kepler's equation (e.g. Newton's Method). See the question How did Newton and Kepler (actually) do it? and the answers there for a little more background.

1. The Python package PyEphem has been around and well supported, and is the pythonic reincarnation of XEphem. I haven't used it, but I believe it keeps enough information about orbital parameters at certain epochs to generate an ephemeris, including some gravitational perturbations. In other words, it's much more than planets moving on fixed elliptical orbits around a fixed sun. So I believe it runs without internet connection.

2. I never used it because I was recommended to look at Skyfield and it's exactly what I needed. It downloads a standard JPL ephemeris that you choose, and then just uses it from your hard drive after that. However, in order to deal with leap seconds and other time related effects, it occasionally needs to check the internet for leap second information updates, since these are arbitrary.

I don't know if Skyfield has a mode to avoid that. Actually that's a good question. If you work with a timescale that doesn't have leap seconds, I am not sure if it will run in its current version.

Both Skyfield and PyEphem Python packages have been written and are maintained by @BrandonRhodes.

If you can allow an occasional connection to the internet (say every month or few months) then Skyfield is extremely easy and pythonic to use. For example here is a script I used in this answer. If you want to convert time formats from system time, you can search my questions in stackoverflow and space exploration. If you want it for your local computers system time now, you just use t = ts.now() directly instead of utc.

import numpy as np
import matplotlib.pyplot as plt

from skyfield.api import load

data   = load('de421.bsp')
ts     = load.timescale()
t      = ts.utc(2016, 7, 5, 9, 50, 0)

jupiter, earth  = data['Jupiter barycenter'], data['Earth']
jpos, epos      = jupiter.at(t).position.km, earth.at(t).position.km
d_instantaneous = np.sqrt(((jpos - epos)**2).sum())

d_light = earth.at(t).observe(jupiter).distance().km  # where WAS Jupiter 48 minutes ago?

clight = 299792.458  # km/s

print "d_instantaneous / c = ", d_instantaneous/clight
print "d_light / c =         ", d_light/clight