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{\mu}{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 = \sin^{-1}\left(\frac{z_e}{\sqrt{x_e^2+y_e^2+z_e^2}}\right) \$
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 graphical-output 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 code-golf, 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:
αinstead ofμfor the gravitational parameter and introduce a symbolp=a(1-e**2)for the semi-latus rectum. There's an error in your cartesian (x,y) velocity formula - it should sayμ(orα) notωsee equations 46 & 47 in the linked page. Also for ECEF to lat/long, the quantity on the right hand of the latitude formula issin ϕnotϕAn additionalasinis needed to convert. \$\endgroup\$ℓ=a(1-e^2)since it's just another symbol to use; I was trying to keep it as simple as possible. thanks for the corrections and wiki page, no idea how I never spotted the kepler orbit page when I was linking everything. \$\endgroup\$pso you could check my correction against the wikipedia formula. I'm going to try to provide an answer this week. \$\endgroup\$H=sqrt(αp)=constant=magnitude of cross product of displacementrand velocityr dotequations (26 and 3.) This reduces to scalar radius * tangential velocity, with the latter being simplyH/r(equation 19) giving the result dTheta/dt=H/r**2(equation 3) which is consistent with Kepler's observation that the orbit sweeps out equal area in equal time. I intend to use bendy curves method with dt/dTheta=r**2/HI tried the other method but this is shorter, more stable & gives the same result \$\endgroup\$