How is the average distance between 2 objects orbiting around a third object calculated? I want to know the average distance between Io and Europa knowing that Io's semi-major axis (around Jupiter of course) is 421,800 km and Europa's semi-major axis is 671,100 km.
At first I thought it's the average of the closest approach (671,100 - 421,800 = 249,300 km) and the farthest approach (671,100 + 421,800 = 1,092,000 km), but the average of those is exactly 671,100 km, which is Europa's semi-major axis, and I find that hard to believe.
Assuming the orbits are co-planar and perfectly circular, which they more or less are, surely there must be a simple formula to calculate the average distance?
I asked this on another site and the answer was &quot;run a simulation!&quot;, which is a very good answer, and I would've already done that if I knew how lol.
If we set the radius of Io's orbit to 1, then Europa's is about $a=$1.591. Since they are in a 2:1 orbital resonance we would expect that number to the 3/2 power to be exactly 2. It's close (2.007) but there's enough of a difference to make it interesting; since Jupiter rotates rapidly I'm guessing that it's related to this.
Okay based on @JamesK's idea of keeping one fixed I tried to get the analytical integral
$$&lt;r_{12}&gt; = \frac{1}{2 \pi}\int_0^{2 \pi}\sqrt{(a-\cos \theta)^2 + \sin^2 \theta} \ d\theta$$
from Wolfram Alpha but I got a

Standard computation time exceeded...

message which I've never seen before, (screenshot) so I did it in Python.
While @JamesK's answer is &quot;the greater of the two&quot; which here would be 1.591 or there would be 671,100 kilometers, because this is a 2D and not a 1D problem and Pythagoras has something to say, I get a different value.
('ratio: ', 1.5910384068278804)
('d.mean(): ', 1.7524934914237922)
('James_K: ', 1.5910384068278804)

The Python script below returns 1.752 or about 739,200 kilometers (solid line) versus the mean projected 1 dimensional distance (dashed line).

import numpy as np
import matplotlib.pyplot as plt

a = 671100. / 421800
print('ratio: ', a)

theta = np.linspace(0, 2*np.pi, 100001)[::-1] # don't double-count the endpoints
d = np.sqrt((a - np.cos(theta))**2 + np.sin(theta)**2)
print('d.mean(): ', d.mean())
theta_degs = (180/np.pi) * theta

plt.plot(theta_degs, d)
plt.xlabel('phase angle (deg)')
plt.ylabel(&quot;distance normalized to Io's SMA&quot;)
plt.plot(theta_degs, d.mean() * np.ones_like(theta_degs), '-k')
plt.plot(theta_degs, a * np.ones_like(theta_degs), '--k')
plt.ylim(0, None)
plt.show()

James_K = ((a-1.) + (a+1.)) / 2.
print('James_K: ', James_K)