====== Calculating orbiting body's position at a given time T ======
{{tag>orbital_elements kepler_orbit orbit}}
=== Creating an orbit for a planet and place the planet on it ===
To create a planetary orbit and place a planet on it so that one can calculate precise position of the planet at a given time can be done by using //Orbital Elements//. You can read detailed information about them from this [[https://en.wikipedia.org/wiki/Orbital_elements|wikipedia page]]. In this guide, I will show detailed steps on how to do the calculation and implement it in Blender or Unity.

There are six parameters in orbital elements that are closely related to the 3-D positions and 3-D velocity of a planet. 

<grid>

<col md="4">
<thumbnail>
{{ proposals:megagrant:tech:orbit1.svg.png?400&direct }}
<caption>
In this diagram, the orbital plane (yellow) intersects a reference plane (gray). For Earth-orbiting satellites, the reference plane is usually the Earth's equatorial plane, and for satellites in solar orbits it is the ecliptic plane. The intersection is called the line of nodes, as it connects the center of mass with the ascending and descending nodes. The reference plane, together with the vernal point (♈︎), establishes a reference frame [from Wikipedia].
</caption>
</thumbnail>
</col>

<col md="3">
  * semi-major axis ($a$): determines the size of an orbit. For a given mass of the central object (e.g., Sun), $a$ is directly related to an orbital period ($P$) through Kepler's 3rd law.
  * eccentricity ($e$): shape of the ellipse. In an elliptical orbit, Sun is located at one of two foci.
  * inclination ($i$): the tilt angle of the orbital plane with respect to observer's line of sight.
  * longitude of the ascending node ($\Omega$): In an inclined orbit (i.e., non-zero $i$), there are two points where a planetary orbit crosses the plane of reference. We choose one of these points as $\Omega$ where a planet moves from below the reference planet to above the plane. 
  * argument of periapsis ($\omega$): determines the orientation of the ellipse in the orbital plane. From the focus where the Sun is located at, the direction toward the periapsis can be anywhere on the plane. An angle measured from the ascending node to the periapsis can fix the orientation of the ellipse.
  * true anomaly ($\nu$): the position of the planet measured in angle from the periapsis for a given time (or epoch).
</col>

<col md="5">
<thumbnail>
{{ topics:orbitalelements_demo.gif?700 }}
<caption>
Simulation showing each orbital element (from pyorb site)
</caption>
</thumbnail>
</col>
</grid>

----

=== Calculation ===
Check the Wikipedia for "Orbit Modeling" for useful information.

From a given set of orbital elements, let's see how we can calculate the planet's position at a given time $T$. At any arbitrary epoch (say $T=0$), orbital elements can be expressed as $(a_0,e_0,i_0,\Omega_0,\omega_0,\nu_0)$. Since the first five parameters do not change over time (without a perturbation like additional object), at a later time ($T+\delta T$), we only need to find the new $\nu$; i.e., at $T=T+\delta T \equiv new$, we have $(a_0,e_0,i_0,\Omega_0,\omega_0,\nu_{new})$ where $\nu_{new}=\nu_0 + \delta\nu$.
 
In this calculation, the reference plane is denoted as $(\hat{I},\hat J,\hat K)$ where $\hat I$ is the direction toward the vernal equinox and the orbital plane is denoted as $(\hat x,\hat y,\hat z)$ where $\hat x$ is toward the periapsis and $\hat z$ is toward the north pole. In the Solar System, we choose the ecliptic plane as the reference plane.

Then, the transformation between the reference plane and the orbital plane is defined by three Euler rotations. 

{{ proposals:megagrant:tech:eulerrotations.svg?direct }} where the new $xyz$ values are {{ proposals:megagrant:tech:new_xyz.svg?direct }}

=== Step by step calculation tutorial ===
See this python script for calculating positions of inner Solar System planet.
++++ Click here to see |
<code python>
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse, Circle

#Set axes aspect to equal as orbits are almost circular; hence square is needed
ax = plt.figure(0).add_subplot(111, aspect='equal')

#Setting the title, axis labels, axis values and introducing a grid underlay
#Variable used so title can indicate user inputed date
plt.title('Inner Planetary Orbits at[user input date]')
plt.ylabel('x10^6 km')
plt.xlabel('x10^6 km')
ax.set_xlim(-300, 300)
ax.set_ylim(-300, 300)
plt.grid()

#Creating the point to represent the sun at the origin (not to scale),
ax.scatter(0,0,s=200,color='y')
plt.annotate('Sun', xy=(0,-30))

#Implementing ellipse equations to generate the values needed to plot an ellipse
#Using only the planet's min (m) and max (M) distances from the sun
#Equations return '2a' (the ellipses width) and '2b' (the ellipses height)
def OrbitLength(M, m):
    a=(M+m)/2
    c=a-m
    e=c/a
    b=a*(1-e**2)**0.5
    print(a)
    print(b)
    return 2*a, 2*b

#This function uses the returned 2a and 2b for the ellipse function's variables
#Also generating the orbit offset (putting the sun at a focal point) using M and m
def PlanetOrbit(Name, M, m):
    w, h = OrbitLength(M, m)
    Xoffset= ((M+m)/2)-m
    Name = Ellipse(xy=((Xoffset),0), width=w, height=h, angle=0, linewidth=1, fill=False)
    ax.add_artist(Name)



from math import *

EPSILON = 1e-12
def solve_bisection(fn, xmin,xmax,epsilon = EPSILON):
  while True:
      xmid = (xmin + xmax) * 0.5
      if (xmax-xmin < epsilon):
        return xmid
      fn_mid = fn(xmid)
      fn_min = fn(xmin)
      if fn_min*fn_mid < 0:
          xmax = xmid
      else:
          xmin = xmid

'''
Found something similar at this gamedev question:
https://gamedev.stackexchange.com/questions/11116/kepler-orbit-get-position-on-the-orbit-over-time?newreg=e895c2a71651407d8e18915c38024d50

Equations taken from:
https://en.wikipedia.org/wiki/Kepler%27s_laws_of_planetary_motion#Position_as_a_function_of_time
'''
def SolveOrbit(rmax, rmin, t):
  # calculation precision
  epsilon = EPSILON
  # mass of the sun [kg]
  Msun = 1.9891e30
  # Newton's gravitational constant [N*m**2/kg**2]
  G = 6.6740831e-11
  # standard gravitational parameter
  mu = G*Msun
  # eccentricity
  eps = (rmax - rmin) / (rmax + rmin)
  # semi-latus rectum
  p = rmin * (1 + eps)
  # semi/half major axis
  a = p / (1 - eps**2)
  # period
  P = sqrt(a**3 / mu)
  # mean anomaly
  M = (t / P) % (2*pi)
  # eccentric anomaly
  def fn_E(E):
    return M - (E-eps*sin(E))
  E = solve_bisection(fn_E, 0, 2*pi)
  # true anomaly
  # TODO: what if E == pi?
  theta = 2*atan(sqrt((((1+eps)*tan(E/2)**2)/(1-eps))))
  # if we are at the second half of the orbit
  if (E > pi):
    theta = 2*pi - theta
  # heliocentric distance
  r = a * (1 - eps * cos(E))
  return theta, r

def DrawPlanet(name, rmax, rmin, t):
  SCALE = 1e9
  theta, r = SolveOrbit(rmax * SCALE, rmin * SCALE, t)
  x = -r * cos(theta) / SCALE
  y = r * sin(theta) / SCALE
  planet = Circle((x, y), 8)
  ax.add_artist(planet)

#These are the arguments taken from hyperphysics.phy-astr.gsu.edu/hbase/solar/soldata2.html
#They are the planet names, max and min distances, and their longitudinal angle
#Also included is Halley's Comet, used to show different scale  and eccentricity
PlanetOrbit('Mercury', 69.8, 46.0)
PlanetOrbit('Venus', 108.9, 107.5)
PlanetOrbit('Earth', 152.1, 147.1)
PlanetOrbit('Mars', 249.1, 206.7)
PlanetOrbit("Halley's Comet",45900,88)
for i in range(0, 52):
  DrawPlanet('Earth', 152.1, 147.1, i/52 * 365.25 *60*60*24)
for i in range(-2, 30):
  DrawPlanet("Halley's Comet",45900,88, 7*i *60*60*24)
print(60*60*24*365)

plt.xlim(-1000,1000)
plt.ylim(-1000,1000)
plt.show()
</code>
++++

See {{ :tech:tutorials:orbitcalculation.pdf |this PDF file}} for the step-by-step calculations.