# Plot the ground path of a satellite

If you model a satellite as a free point orbiting a body, you can pretty easily see it has 6 degrees of freedom: three for the X, Y, and Z position, and three for the X, Y, and Z velocity. However, using those 6 numbers doesn't tell you much about the orbit (and are constantly changing over time for any one satellite), so instead we usually define the shape of an orbit with 6 variables ( + time for 7 total) called the 'orbital elements' which together describe the shape of the orbit:

Things you don't care about: The ascending node (☊) is the point where the body passes through the equatorial plane of its planet such that it is 'above'/'north' afterwards. The reference direction(♈︎) is some known vector on the equatorial plane, so may be defined differently on different bodies. Earth's reference direction always points towards where the sun would be during the vernal equinox.

The six elements (usually) used to define an orbit are:

• a: Semi-major axis. The size of the orbit.
• e: Eccentricity. How circular/stretched the orbit is.
• i: Inclination. The angle between the orbital plane and the equatorial plane.
• Ω: Longitude of the ascending node. The clockwise angle on the equatorial plane(as viewed from the positive Z axis / north pole) between the reference direction and the ascending node.
• ω: Argument of periapasis. The angle, on the orbital plane, from the ascending node to the periapsis (closest approach to the body).
• ν: True anomaly. The angle, on the orbital plane, from the periapsis to the satellite.

If you look closely, you can see that two of these (Ω and ν) are boring if you're plotting lat/long points, since they just shift the graph left or right. To rectify this, we are setting these two parameters to zero (eg line of nodes == reference direction && satellite starts at periapsis), replacing them with the standard gravitational parameter μ and the length of the sidereal day I. Finally, we will define our time t as how long it's been since the last solar noon at a latitude of 0° on the vernal equinox.

### Useful equations (or: What off Earth do I do with all these inputs?!?)

Radius: $$\ r = \frac {a (1-e^2)}{1+e \cos ν } \;\; {\rm{or}} \;\; \vec r = \frac {a (1-e^2)}{1+e \cos ν } \begin{bmatrix} \cos \nu \\ \sin \nu \end{bmatrix} \$$

Velocity: $$\ r' = \sqrt {\mu \left( \frac{2}{r} - \frac{1}{a} \right) } \;\; {\rm{or}} \;\; \vec r' = \sqrt{\frac{\omega}{a (1-e^2)}} \begin{bmatrix} \sin \nu \\ e + \cos \nu \end{bmatrix} \$$

Acceleration: $$\ r'' = \frac{\mu}{r^2} \;\; {\rm{or}} \;\; \vec r'' = - \frac{\mu}{r^3} \vec r \$$

Orbital period: $$\ T = 2 \pi \sqrt {\frac{a^3}{\mu}} \$$

Coord transfer from perifocal to equatorial (uses our assumption of Ω = 0):

\ \begin{align*} \begin{bmatrix} x_e \\ y_e \\ z_e \\ \end{bmatrix} &= \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos i & -\sin i \\ 0 & \sin i & \cos i \end{bmatrix} \, \begin{bmatrix} \cos ω & -\sin ω & 0 \\ \sin ω & \cos ω & 0 \\ 0 & 0 & 1 \end{bmatrix}\, \begin{bmatrix} x_p \\ y_p \\ 0 \\ \end{bmatrix} \\ {\rm or} \;\;\;\;&= \begin{bmatrix} \cosω & -\sin\omega & 0 \\ \cos i\sin\omega & \cos i\cos\omega & -\sin i \\ \sin i \sin\omega & \sin i \cos\omega & \cos i \\ \end{bmatrix} \begin{bmatrix} x_p \\ y_p \\ 0 \\ \end{bmatrix} \end{align*} \

From ECEF to lat/long: $$\ \lambda = \operatorname{atan2}(y_e,x_e), \phi = \frac{z_e}{\sqrt{x_e^2+y_e^2+z_e^2}} \$$

There are many ways you could go about implementing this challenge, but for now I'll outline the two ways I think would be easiest. You do not need to use either of these methods if you don't want to; but if you're feeling stuck they're here to help. If you think of some other way, but get stuck bc you need an equation I haven't added yet, let me know in a comment and I'll try to add it. Or if you spot an error in my equations, please let me know.

##### Method the First: The Long March

Uses vector form of acceleration equation + initial conditions. Start by picking some small Δt to be your step size. You add velocity * Δt to the position, and then add acceleration * Δt to velocity. If you keep going back and forth like this you get a list of points evenly spaced in time. With this method it doesn't matter if you convert from perifocal to equatorial before or after the march. Once you have the points in the equatorial coordinate system, you can rotate each point around the z axis by angular velocity * the time corresponding to that point to put it into ECEF. With your points in ECEF, you can use that last equation to end up in lat/long.

##### Method the Second: Bendy Curves

Define a curve in perifocal coords by replacing r and ν in the 'radius' equation with sqrt(x^2+y^2) and atan2(y,x). Convert the curve from perifocal to equatorial coords. Starting from the periapsis, walk along th curve in Δt*(velocity from the equation) sized steps, rotating each point by angular velocity * time corresponding to that point. This puts your point into ECEF, so you can then use the last equation to end up in lat/long.

### I/O requirements

##### Input

7 total inputs: eccentricity, semi-major axis, inclination, argument of periapasis, standard gravitational parameter, length of the sidereal day, and initial time.

You may accept them in any order you want, and in whatever form you want. You may even swap individual elements with equivalent constructs, eg taking a true anomaly or mean anomaly as input instead of initial time. (do ask before doing anything too crazy though)

##### Output

This is a challenge, so you should output a graphical plot of the equatorial (not geodetic) latitude vs longitude of the ground path of the satellite. (so, technically, it's declination rather than latitude, but it's easier to just say lat/long)

If your language has no concept of graphics, like bf or sed; first of all, why would you choose that language for this challenge; and secondly, outputting a list of points work okay. (such that if someone else did plot those points as on line graph, it would look correct)

Your plot should cover at least a time period equal to the lcm of orbital period and sidereal day length. This means you are plotting the full path of the orbit relative to the ground, as after that it (should) repeat. If you're march along the curve, please try to keep error low enough that the end looks like it basically connects to the start.

##### Edge cases

If a value is undefined, you can treat it as if it were. For example, if e is 0 (meaning ω is undefined), ω will be input as 0, and you can treat it as if it was 0. (eg the rocket starts at the ascending node)

The semi-major axis, standard gravitational parameter, and length of the sidereal day are guaranteed to be positive; the eccentricity, argument of periapsis, inclination, and initial time are guaranteed to be a non-negative number less than one, 360, 180, and ∞ respectively. Or 1, 2π, π, and ∞ if you take angles as radians.

This is , so lowest bytes in each language wins

### Test cases

Inputs:

e = .737
a = 26553
i = 63.4
ω = 270
μ = 398600
I = 86164
t = 19147


Output (molniya orbit):

Inputs:

e = .1
a = 42164
i = 10
ω = 20
μ = 398600
I = 86164
t = 0


Output (geosynchronous orbit):

Inputs:

e = 0
a = 1
i = 90
ω = 0
μ = 39.5 (= (2pi)^2)
I = 2
t = 0


Output:

Inputs:

e = 0
a = 1
i = 0
ω = 0
μ = 39.5
I = 1
t = 0


Output ([not geo, but] -stationary orbit):

Inputs:

e = 0.5
a = 1
i = 50
ω = 0
μ = 39.5
I = 1
t = 0


Output:

Inputs:

e = 0.75
a = 2
i = 130
ω = 180
μ = 79
I = 1
t = 0.5


Output:

Inputs:

e = 0.5
a = 2
i = 40
ω = 180
μ = 79
I = 5
t = 1.25


Output:

Inputs:

e = 0.5
a = 40
i = 60
ω = 180
μ = 787
I = 20
t = 0


Output:

• sandbox Commented Jul 5 at 12:22