Orbital calculation headaches

[Maltman]

I'm trying to get more into the math of orbiting Kerbin. 

One concept I am struggling with is how to come up with all the orbital information when you don't know the periapsis or apoapsis or semi major axis or anything. 

I want to know how to calculate the periapsis, apoapsis ( and all that the game always gives you) using only the velocity vector, distance from center of planet, and the planets mass. (AKA using only your position, your velocity vector and the planets gravitational force).

 

Please help me

And thank you. 

[Freshmeat]

Check this

I cannot give a full answer, but on the top of my head: v squared gives kinetic energy, distance potential energy. Your semi major axis is dependent on only mechanical energy. From the velocity vector and distance to center you get angular momentum which is conserved. Thus you should be able to get ap and pe by solving for a velocity and distance that satisfies both angular momentum and mechanical energy. I'm off to bed, and will be off for a day or two, but I think it is doable using only algebra.

[Maltman]

I understand how to get the velocity of a spacecraft using the range and semi major axis.  These are just for the magnitude of velocity and distance though and I want to figure out how to get the semi major axis using velocity and distance vector. 

 

I am having a hell of a time finding any information about this on the internet oddly enough. I would think it would be one of the most important concepts.

 

I should have payed more attention in school!

[bewing]

Square both sides. Divide both sides by GM. Subtract 2/r from both sides. Then negate and invert both sides. That should give you "a" in terms of V and R.

[Maltman]

3 hours ago, bewing said:

Square both sides. Divide both sides by GM. Subtract 2/r from both sides. Then negate and invert both sides. That should give you "a" in terms of V and R.

I got this equation from that method . I can find the length of the semi major axis through this and then the distances for the periapsis and apoapsis.  Is there any way I can determine the coordinates of where the semi major axis, periapsis, and apoasis will be?

[bewing]

OK, assuming you got the math right so far, then yes. It's not too hard to describe in words.

Your ship is on an ellipse. So you make up a coordinate system with the x axis along the Pe-Ap line, in the plane containing your ship, with the origin at the center of the ellipse. The basic equation of your ellipse is: (x + (Pe - Ap)/2))^2 / a^2 + y^2 / b^2 = 1. In the new coordinate system that is.

Eccentricity is probably a good thing to calculate:

From that, you should be able to easily calculate b (the semi-minor axis). And you basically know x and y for your ship's current position. However, it's in the wrong coordinate system.

So you also need a "transform" for the coordinate system. If you aren't familiar with those -- they are a bit annoyingly complicated, but look it up on wikipedia. To calculate it:

The cross product of your velocity vector with the Kerbin radius vector always gives you a vector perpendicular to your new coordinate system's XY plane -- ie. the new Z-axis vector. But there is still one more rotation angle. And you need to solve backwards for that -- using the ellipse equation, and solving for x and y simultaneously.

Once you have calculated the coordinate system transform, then you just plug in (Ap - (Ap - Pe)/2, 0, 0) and (-Pe + (Ap - Pe)/2, 0, 0) and solve backwards to find their locations in your original coordinate system.

 

Edited June 10, 2017 by bewing

[Cunjo Carl]

Hmmm.... This is a thing I've wanted to do for ages!  My goal is to provide the math for @Maltman's request with high school level math/physics concepts (the kind I work best with). It'll be long, but hopefully intuitive and informative.

 

Spoiler

Let's make this derivation in terms of polar coordinates... And let's do it 2d! Finally, let's do it in terms of the fictitious rotational forces, because I think it'll make more intuitive sense than the other stuff I've tried to figure out (wiki / old texts). So... We'll solve for our object's motion in terms of Gravity, Centrifugal and Coriolis forces. The advantage of doing it this way is we won't (or at least I think we won't) need to deal with any confusing reference frame stuff this way. It'll just be F=ma with forces for all of its simplicity and occasional nastiness.

We can handle our position and velocity variables like this: we're at a certain distance (r) from the planet while traveling a certain velocity (v) and traveling at a certain angle compared to the 'horizon' on the planet ( φ ). So, for φ, if we're flying level with the planet (eg. circular orbit, or at apoapsis, periapsis) φ  = 0, and if we're flying straight up φ = pi/2 (aka 90 degrees). Finally there's θ , which will be our present angle in the orbit relative to our starting point. The initial values will be r₀, v₀,  φ₀ and θ₀. Present variables will always be bold, while constants will be normal type.

For convenience, we can re-package the information from v and φ by turning it into v_r and v_θ , meaning the velocities of our object going up/down and round/round. The conversion's the simple standard:
v_θ = v _*cos(φ)
v_r= v _*sin(φ)

Great! The stage is set. Let's plug in the forces F=ma style. Because it's the simplest part, we'll start with the rotational direction θ. The only force in the θ direction will be Coriolis force, however that one's confusing if you've never had to deal with it before, so let's talk about it. If you're on a ... spinny thing at the kid's park.... you push them and they spin round, and you get all sick? Anyways, let's say you're near the edge while it's spinning, and you start pulling yourself inwards. What happens? You're doing work to pull in (force*distance on centrifugal force) and in exchange, the spinny thing pushes you in the direction of rotation. Our velocity actually increases as we pull ourselves in, and we spin round faster as a result. Coriolis force!

Coriolis Force = m*v_θ*v_r/r in the direction of rotation. So, in order to feel a coriolis force we need to be spinning, and then it's proportional to the relative change in radius. I don't have my derivation with me anymore, but I'm pretty sure this was it. That said, I'll be using a workaround with an angular momentum balance instead until I'm sure I've got it right.

Centrifugal Force = m*v_θ²/r in the outwards direction

Gravity = G _* M _* m / r² in the inwards direction. G is the gravitational constant and M is the mass of our planet.

All of our forces are proportional to our object's mass, so we can write things immediately in terms of accelerations where: a = F/m

a_θ = -2v_θ*v_r/r                                  (a = coriolis)  (probably... sorry, I'll doublecheck this one soon! I use another approach momentarily)
a_r = v_θ²/r   -   G_*M/ r²                   (a = centrifugal out - gravity in)

Now the next step would be to find v_θ as a function of r. Why? We'd like to avoid solving this as a parametric (simultaneous) system, so if we can solve for that v_θ hanging out in a_r 's equation, we can solve for a_r much more simply. You can think of this as: centrifugal force should only depend on the initial conditions and the present radius. That makes sense! I want to do with an energy balance on coriolis force eventually, but for the moment, we happen to know the answer from conservation of angular momentum is:

v_θ = v_θ,0 * r₀ / r    (because L = constant = m _*v_{θ *}r at all times. This means v_{θ *}r = v_θ,0 * r₀ , and we get our equation from here   I will replace this later - It's just Coriolis force in disguise!)

Plug it in~!

a_r = v_θ,0² * r₀² / r³   -   G_*M/ r²                        This is our acceleration away from kerbin knowing only our present radius and initial starting conditions

And this is great- it relates exactly the things we want to know in a simple way without fuss... and it happens to be completely unsolvable. Thanks math. Even if we can't solve it directly, it does make a great vantage point for us to stop and take a look to make sure we're on the right track. Let's see. This makes sense! When you get too far in, centrifugal force quickly pushes you out (inverse r cubed) . When you get too far out, gravity eventually pulls you back in (inverse r squared). When orbiting, you tend to swing past the periapsis quickly and the apoapsis slowly, so I'll claim we're on the right track. I'm lucky I had to go make pancakes at this point, because that's a nasty integral. The next couple steps took me the rest of the day to puzzle through, and I think it's important to realize that these next steps don't gain us ANYTHING conceptually. They're purely mathematical manipulations to cram our formula into one of the solvable forms for a second order ordinary difeq.

have    a_r = ~~~/r³ - ~~~/r²       Right now the formula looks like this, which is nasty unsolvable analytically
want     a_u =      -u + ~~~      If it instead looked like this, it would be solvable as a standard harmonic oscillator (sin and cos) Let's get it into this form.

We would adore some variable substitution to pop an r³ in there, but sadly there are none. Because of this, we have to make a concession and solve for just the orbit's shape, not its motion as a function of time. It's pretty yuck, but it's a hallmark of any orbit derivation. Over the centuries, people have made hundreds of increasingly complicated ways to partially undo the damage done in this next tragic-but-necessary-step. We'll use Leibnitz notation to convert our equation from time domain, where it tells us exactly what our craft does over time, into theta domain where it will only tell us the shape of its path. We'll also introduce a helper variable u = 1/r just to make the math work nicely.

have   a_r = d²r/(dt)²                    This is just the definition, showing off Leibnitz notation.
want    a_u = d²u/(dθ)²

Ok, so to get from where we are to where we want to be we'll need to convert from t to θ, and from r to u. First we'll make a conversion factor of (dt/dθ)² to convert the domain (t->θ) Leibnitz style, and then a du/dr factor to convert to our helper (r->u).

dθ/dt * r = v_θ = v_θ,0 * r₀ / r         That dθ/dt will help in our conversion. dθ/dt is the rate you're spinning in circles expressed in radians, and if you multiply it by your radius, you'll get your radial velocity. Hey we have a handy relationship for radial velocity from before thanks to Coriolis force / angular momentum.
dθ/dt = v_θ,0 * r₀ / r²                    In this step we solve for dθ/dt. In the next we invert and square it to make the conversion factor.
(dt/dθ)² = r⁴/ (v_θ,0 * r₀)²             [1]  This factor will let us convert to theta domain. We just need to find the conversion factor for the u = 1/r helper variable now.

u = 1/r
du/dr = -1/r²                              [2]  Solving for du/dr was much easier.

All set! Let's multiply our a_r equations with the conversion factors [1] and [2], and replace any remaining 1/r with u. Finally, move the u across to the left, and we'll have the classic harmonic equation:
a_r = d²r/(dt)²=  v_θ,0² * r₀² / r³   -   G_*M/ r²                           reminder of the ar equation. Next we plug in the equations [1] and [2].

d²r/(dt)² * (dt/dθ)² * du/dr  =   (v_θ,0² * r₀² / r³   -   G_*M/ r²) * (r⁴/ (v_θ,0 * r₀)²) * (-1/r²)        Plugged in equations... and then it simplifies to!

d²u/(dθ)² + u  =  G_*M/( v_θ,0² * r₀² )                                   This is the happy form, which shows up for all ideal harmonic oscillators. We can now get our solution from an integral table:

u(θ) = G_*M/( v_θ,0² * r₀² )  + C₁*sin(θ) + C₂*cos(θ)             Ok, next we need to solve for those integration constants... Let's make θ₀ = 0. at t = 0. then we can plug in r₀ to get C₂:

 

u₀ = 1/r₀ = G_*M/( v_θ,0² * r₀² ) + C₂                                                    Solving for C₂

C₂ = 1/r₀ - G_*M/( v_θ,0² * r₀² )                                                                                   Gotcha! No time like the present, plug it in.

u(θ) = G_*M/( v_θ,0² * r₀² )  + C₁*sin(θ) + ( 1/r₀ - G_*M/( v_θ,0² * r₀² ) ) *cos(θ)            Working equation...

 

Alright, the next one's gonna be the trick. We need to bring in our starting angle φ₀, which will interact with the first derivative of u(θ) to generate C₁.  When we take the first derivative, the constant drops out, the cosine turns into a -sine (which is 0 at θ=0) and the sine with C₁ into a cosine (which is 1 at θ=0)

du/dθ (evaluated at θ = 0) = C₁                                                            We're gonna need 3 derivatives for piecing together du/dθ, Fortunately, all of which we've determined already.

1. du/dr = -1/r²           2.   dr/dt = v_r= v _*sin(φ)        3.  dθ/dt = v_θ/r = v _*cos(φ)/r   .....   from these    Eq1 * Eq2 / Eq3   gives du/dθ
du/dr * dr/dt / (dθ/dt ) = (-1/r²) * v _*sin(φ) / (v _*cos(φ)/r)                   Simplify...
du/dθ = -tan(φ)/r                                                                                 And finally, evaluating at t = 0...
C₁ = -tan(φ₀)/r₀                                                                                    This was easier than I feel like it should have been. I'll take it! Wrong units. Huh. Good to go.

 

u(θ) = G _* M/( v_θ,0² * r₀² )  + (-tan(φ₀)/r₀)*sin(θ) + ( 1/r₀ - G _* M/( v_θ,0² * r₀²) )*cos(θ)   Last step, convert back to r = 1/u

r(θ) = 1/[  G _* M ( v_θ,0² * r₀² )  + (-tan(φ₀)/r₀)*sin(θ) + ( 1/r₀ - G _* M/( v_θ,0² * r₀²) )*cos(θ) ]  Hee! Does it work? .... bug fixing.... bug fixing.... Now it does. Er, I mean, I always get it right first try!

 

1/  ( (3.5*10^12)/(2275^2*675000^2) +  (-tan(3.14/16)/675000)*sin(t)   (1/675000 - (3.5*10^12)/(2275^2*675000^2))*cos(t) )  Checks out. I'll need to double check exact values later, but it seems to fit the bill.

r(θ) = r₀²/[  G_*M/ v_θ,0²  - r₀*tan(φ₀)*sin(θ) + ( r₀ - G_*M/ v_θ,0²)*cos(θ) ]           One simplification form.... I think this form is my favorite, but lastly:

-------------------------------------------------------------------------------------------------------------------------
|     r(θ) = r₀/[  G_*M/(r_0*v₀²_*cos²(φ₀))  - tan(φ₀)*sin(θ)  +  (1 - G_*M/ (r_0*v₀²_*cos²(φ₀)))*cos(θ) ]       |                  Yatta!
-------------------------------------------------------------------------------------------------------------------------

=============================================================
||   r(θ) = r₀/[ cos(θ) + (G_*M/(r_0*v₀²_*cos²(φ₀)))*(1-cos(θ))  - tan(φ₀)*sin(θ)  ]    ||                                                (another alternative form, simplified for the constants, easier to play with)
=============================================================              

 

So we finally have our full orbit equation in terms of only initial conditions. At some point, if the equations don't spaghetti  I'll complete the derivation of orbital parameters from initial conditions here.

Day 2:

OK, in a lot of ways the double boxed equation is super ideal! If you type it in to a calculator you can quickly go from your present conditions to an equation for your entire orbit. The question was how to get specific orbital values though, so in this next part we'll shoe-horn this equation into the standard format so we can find periapsis, argument of periapsis, eccentricity, and so on.

First off, what form are we trying to fit?

have r(θ) = r₀/[    P  + Q*sin(θ)  +  (1 - P)*cos(θ)    ]                       (the single boxed formula we just derived with convenience constants P and Q)
want r(θ) = a*(1-e²) / ( 1 - e*cos(θ + α ) )                                       (the standard formula with a phase shift so our variables will fit. Has convenience constants a=semimajor and e=eccentricity)

   We had to add the  α  into the standard form because our theta is 0 at t=0, rather than 0 at periapsis. The first step is going to be cramming the sin and cos terms together into one. Trigonometry textbook formula page go!

a*sin(x) + b*cos(x) = R*cos(x-α)   ,   R=(a² + b²)^1/2  ,   α=arctan(a/b)  

Now we just fill the form:
R=(Q² + (1-P)²  )^1/2       where P = (G_*M/(r_0*v₀²_*cos²(φ₀)))  and Q = -tan(φ₀)           We're just plugging our convenience variables in.
α = arctan(Q / (1-P) )                                                                                               This has a tan in an arctan... I want to fix it but I can't

r(θ) = r₀/[   P  + R*cos(θ-α)  ]              We're plugging in our new simplified form with convenience variables composed of convenience variables! How deep will the rabbit hole go? We'll get apoapsis and periapsis from this equation. The cosine has a max of +1 and a min of -1, which correspond to periapsis and apoapsis respectively.

r(θ) = (r₀/P)/[ 1 - (-R/P)*cos(θ-α)  ]         Here we have our equation shoe-horned in to the standard form. From this we can get eccentricity, true anomaly and semi-major axis. Argument of periapsis could also be determined if we had a set reference angle, but we never declared one. I suppose we could now, but eh?

 

Periapsis = r₀ / (P+R) =   r₀ / (P+(Q + (1-P)²  )^1/2)
Periapsis =   r₀ / [(G_*M/(r_0*v₀²_*cos²(φ₀))) +(tan²(φ₀) + (1-(G_*M/(r_0*v₀²_*cos²(φ₀)))²  )^1/2]    It's honestly not too bad, but it's a bit spaghetti.

Apoapsis = r₀ / (P-R) =   r₀ / (P-(Q + (1-P)²  )^1/2)
Apoapsis =   r₀ / [(G_*M/(r_0*v₀²_*cos²(φ₀))) - (tan²(φ₀) + (1-(G_*M/(r_0*v₀²_*cos²(φ₀)))²  )^1/2]    There's other forms, but I haven't found any simpler

eccentricity = -R/P  = -(Q² + (1-P)²  )^1/2 / P 
eccentricity =  -(tan²(φ₀) + (1-(G_*M/(r_0*v₀²_*cos²(φ₀)))²  )^1/2 /  (G_*M/(r_0*v₀²_*cos²(φ₀)))        These are all easiest to just calculate P, Q and R first.

Semi-major axis = r₀ / (P*(1-(R/P)²))                                                                                   This one in particular!

True anomaly = -arctan(Q / (1-P) )
True anomaly = arctan( tan(φ₀) / (1-(G_*M/(r_0*v₀²_*cos²(φ₀)))) )

All this said, I've had way worse as far as spaghetti goes. Hey! We set out to find our orbit's critical values in terms of present position/velocity and here we have it! We built a boat and found ourselves on the far shore. What to do now we're here? I think I'll just enjoy the view. Hopefully for you this is enough to take the next trek farther!

Day 3:  (WIP)

*cough cough*

Speaking of enjoying the view, I couldn't help but consider what we could learn from the equations we generated. In particular, the double-boxed equation is surprisingly concise. It wraps up the entire orbit using our starting variables only once or twice each, and can be invoked from any point in the orbit. Given the complexity of orbit and of orbital equations, that's surprising!

Let's look at that thing I dubbed 'P', which is central in the double-boxed equation
P = G_*M/(r_0*v₀²_*cos²(φ₀))
Q = - tan(φ₀)

I actually wrote a section about how it's an energy ratio between escape energy and rotational energy  E_grav = GMm/r  and   E_rot = 1/2 m v_θ²   ... and while it is (mind the factor of 2) there's actually a better explanation.... It's a force balance!! I may be the only person for whom this seems insane- I wish I could convey it. There's a force balance in an equation that holds for all time. Forces are present, all time is not. This isn't a thing that I've ever seen happen for a system so complex. Energy sure, but force?? I think I'll leave the mistaken section at the bottom because I think it's funny. Also, the reason it isn't clear what these things mean is that they started as integration constants, which need to be divined from starting conditions through mathematical tomfoolery, they don't start as physical relationships. Ok, enough gabbing. Time to copy paste.

Gravity = G _* M _* m / r²
Centrifugal Force = m*v_θ²/r
P = Gravity/Centrifugal

Then we notice the only difference between centrifugal force and coriolis force is that a v_r and v_θ are swapped. Suddenly, the tan(φ) makes sense!
Coriolis Force = -2m*v_θ*v_r/r
v_θ = v _*cos(φ)   ,   v_r= v _*sin(φ)    ,    v_r/v_θ = tan(φ)
Q = Coriolis/Centrifugal    (well ¹/₂Q the way I defined it)

And then, in the greatest logic leap I've made in a while, I'll claim that:
1 = Centrifugal/Centrifugal

So plugging it all back in, we see:

r(θ) = r₀  *    F_{Centrifugal,0} /[ F_{Centrifugal,0} * cos(θ) + F_Gravityl,0*(1-cos(θ))  - F_Coriolis,0*sin(θ)  ]     Then we remember that the only difference between these trig functions is the phase shift.

r(θ) = r₀  *    F_{Centrifugal,0} /[ F_Gravityl,0*(1+cos(θ-180deg))  - ¹/₂F_Coriolis,0*cos(θ-90deg) - F_{Centrifugal,0} * cos(θ-180deg) ]  

This is it! Values in the denominator increase your radius when they're small or negative, and increase your radius when they're big and positive. Because of that, I put minus signs on the forces that generally hold you out, so when that value's phase comes about is when it's doing its thing the most. We posited in the begining that orbits are caused by the interplay of these three forces, and now we can see the nature of that interplay very plainly. Hm.. maybe the best way to convey it would be a chart?

Total Coefficients of:
                        F_Gravity             F_Centrifugal                 F_Coriolis
Now                      0                      1                          0
90 deg                 1                       0                        -¹/₂
180 deg               2                      -1                         0
270 deg               1                       0                         ¹/₂
360 deg                0                      1                          0

Remember, negative numbers mean out! So, this is a measure of how much your initial radius changes relative to whatever imbalance in the forces there are.

too tired to go on. committing edit.

 

 

 

-----------------------------

If you look at it just right, you can see that it's actually an energy ratio, between the energy we have sunk into the gravity field versus the energy we have available as rotational kinetic energy.
E_grav = GMm/r              This is the energy required to escape the gravity field from this point
E_rot = 1/2 m v_θ²   =    1/2 m v² _*cos²(φ)         This is our kinetic energy in the rotational direction

E_grav,0/E_rot,0 = 2*G_*M/(r_0*v₀²_*cos²(φ₀)) = 2P     So 2P is this interaction between these two energies. Perhaps we might expect 2P to appear somewhere in our equation? Let's look.

r(180deg) = r₀/[ cos(180deg) + (G_*M/(r_0*v₀²_*cos²(φ₀)))*(1-cos(180deg))  - tan(φ₀)*sin(180deg)  ]         r(180deg) tells us our position when we're on the opposite side of the planet
r(180deg) = r₀/[ 2P - 1  ]                 We can see that if 2P > 1   (aka if E_rot > E_grav)   then we'll escape before we reach the far side (because the denominator will cross 0) 

That makes plenty of sense. Interestingly, our position on the far side of the planet doesn't depend on how fast up/down we're going, only on how fast round we're going. So if you burn normal or antinormal, your position on the far side of the planet shouldn't change.... neat.

Can we re-phrase the base equation to be more evocative of this?

r(180deg) = r₀ / (E_grav,0/E_rot,0  - 1)
r(180deg) = r₀ * E_rot,0 / (E_grav,0  - E_rot,0)     Well, there we go!
r(θ) = r_0*E_rot,0/[ E_rot,0 * cos(θ) + E_grav,0*(1-cos(θ))  - E_rot,0 tan(φ₀)*sin(θ) ]          I'm not immediately convinced.... that tan(phi-not) would need to be squared to turn Erot into Evertical ... So...

oh.

-------------------------------

 

 

 

Edited June 14, 2017 by Cunjo Carl

[xendelaar]

holy kraken @Cunjo Carl that was a very thorough derivation. must have taken you ages to write this paper. it's very clear and i hope @Maltman can use it to his advantage. 

[Cunjo Carl]

@xendelaar Thanks! In retrospect, the whole derivation from base principals may have been a little overkill ... Hopefully it makes for some nice food-for-thought! Thanks for checking through, by the way.

 

I did go back and finally solve for the standard orbital parameters (like the positions of apoapsis, periapsis) in terms of initial conditions, and just editted it in after the original derrivation of the orbital equation. Though I call it "a bit spaghetti" in the post, I'm actually surprised how clean the equations wound up now that I've started going back to look over other people's derivations (like wiki). Eccentricity is a fairly complicated little convenience variable, so an equation that can catch it from conditions on any point of the orbit will necessarily be a bit complicated. In the end it was all pretty fun! Please let me know if there's any questions, and I understand if it was a bit... beyond what's directly useful . Cheers!

 

Oh, @bewing, I'm burning with curiosity now. Do you happen to know how KSP goes from position/velocity to the orbit?

Edited June 12, 2017 by Cunjo Carl

[bewing]

9 minutes ago, Cunjo Carl said:

Oh, @bewing, I'm burning with curiosity now. Do you happen to know how KSP goes from position/velocity to the orbit?

Big linear algebra matrix math things that convert directly from parametrized keplerian orbital values to whatever local coordinate system is being used, and back.

 

[Cunjo Carl]

I think I finally have it boiled down to the base concepts. It seems like @Maltman is already set, so I'm a little hesitant to post it here, but without kicking a new thread (and where would it go?), I'm not sure where else to put it! Well, here we go:

Orbits for those finally tired of going round in circles

 

r(θ) = p/(1+e_*cos(θ-α)) 

This equation tells you the entire shape of your orbit in terms of 3 constants, which you can calculate from any point in your orbit. Those constants determine shape, size and direction:
e "eccentricity" which defines the shape
p"semilatus rectum" which defines the size
α "longitude of periapsis" which defines the direction the short end points.

Tasty elements of orbit like apoapsis and semimajor axis can all be found as simple relationships to these 3 constants, and in the next few paragraphs we'll see how these terms arise from from the fundamental concepts of an orbit. For a refresher: Orbits work by having you spin 'round a planet, and your centrifugal force holds you up even though gravity pulls you down. If these two forces happen to balance eachother your orbit will be a pleasant circle, but if they're out of balance, you'll bounce back and forth between them like a mass on a spring. At periapsis, extra centrifugal force becomes upwards momentum and you rise, then at apoapsis the higher relative gravity becomes downwards momentum and you fall. Then the cycle repeats, and because it happens the same way every time your craft orbits in the same shape every time. So the shape of the orbit must be based on the forces that make it happen. Let's see how.

Cheat Sheet:
 

Spoiler

 

r = radius from the planet
v_r = speed going away from the planet (dr/dt)
θ = present angle going 'round the planet
ω = rate your angle is changing (dθ/dt) "angular velocity"
v_θ = speed going around the planet (r*dθ/dt)

m = mass of your craft
M = mass of the planet
G = the gravitational constant

P_Up  (Upwards Momentum) = m _*v_r
F_G (Gravity) = G _* M _* m / r²
F_C (Centrifugal Force) = m*v_θ²/r

v = total velocity of your craft
φ = angle between velocity and horizon
v_r = v _* sin(φ) 
v_θ = v _* cos(φ)

v_θ = v_θ,0 * r₀ / r    Because angular momentum is conserved, we go 'round faster the further we go in. This mathematically describes the relationship. 0s mean initial, and bold means present.

 

 

1. Eccentricity: The Shape
e = sqrt[ (ω_*P_Up/F_G)² + (F_C/F_G - 1)² ]

The eccentricity measures the ovalness of an orbit, and to rephrase the picture painted in the intro, it captures this cycle: F_C >F_G ... we go up ... F_C <F_G ...  we go down... [repeat]. So, our constant 'e' needs to be able to catch the essence of this at any point in time. It handles this by adding these two squared terms:
(ω_*P_Up/F_G) , which is the upwards momentum relative to the force of gravity (then multiplied by the rate your angle is changing to make the units work out right... <waves hands>). Any amount of upwards/downwards momentum causes imbalance, but the amount is relative to the gravity it'll be pushing against.
(F_C/F_G - 1) When centrifugal force and gravity are equal to eachother, F_C/F_G becomes 1, and the term falls to 0. As either becomes greater than the other the term gets a greater and greater magnitude.
Since these two things are 'added in quadrature' (a² + b² = c² hypotenuse style), it doesn't matter if the terms are positive or negative, it only matters how big they are. And with this, all possibilities for imbalance are captured and brought together, and can be measured anywhere in the orbit!

2. Semilatus Rectum: The Size
p = r_*F_C/F_G

The semilatus rectum measures the size of an orbit. It tells us the radius your orbit will have in that single moment of the cycle when gravity and centrifugal force balance. We can measure what this magic radius should be from anywhere in the orbit by taking our current radius and multiplying it by a fudge factor: F_C/F_G . Why? By serendipity, F_C/F_G happens to be proportional to 1/r, so whenever our radius increases, the fudge factor decreases by exactly the same amount. The product is the same as our radius when F_C/F_G = 1, or in other words, when the forces balance. In this way we can measure this special radius from any point in our orbit.
 

3. Longitude of Periapsis: The Direction
α = θ - arctan( (ω_*P_Up/F_G) / (F_C/F_G - 1) )

The Argument of Periapsis measures the direction (angle) of our orbit's periapsis relative to some pre-chosen galactic standard. Notice the (ω_*P_Up/F_G) and (F_C/F_G - 1) terms are just the measurements from before for how far out of balance we are through either up/down momentum or imbalanced forces. Their ratio tells us where we are in the cycle, and the arctan turns that information into an angle. That whole piece tells us where the periapsis is relative to us, so when we add it to θ (which is where we are relative to the standard) we get the longitude of periapsis: the angle our periapsis points relative to the standard. Because it's always relative to where we are, we can figure out the direction of the orbit from anywhere along it.

 

And that's it! Now we know the shape, size and direction of our orbit we can plot it or gather other information through simple relationships. I hope it's brought some insight to the nature of these seemingly ineffable constants!
 

Edited June 15, 2017 by Cunjo Carl

[Maltman]

Hey thanks for all these quality posts! I especially appreciate yours Cunjo! 

 

I am struggling to understand it all. I am also making a 2D orbit game and I want to code some of this math into the game so that the player can get information on his orbit. I find all of this shockingly complex sometimes and rather overwhelming. So far I have the eccentricity down as well as the angular momentum. I'd like to code an entire orbital read out with prediction like you get in KSP or in Kerbal Engineer.

 

I've just gotten over the cold. I have to devote some quality time to understanding these concepts. I almost need someone to explain it over a discord or IRC for me because I am derp.

Edited June 15, 2017 by Maltman

[Cunjo Carl]

15 hours ago, Maltman said:

 I almost need someone to explain it over a discord or IRC...

I'm in! My free time is pretty sparse, sadly. Do you have availability tonight 9-10pm PDT (4-5am GMT)? If not, I can try to swing a time this weekend.

Just like anything, it's tricky stuff when you first get going, but gets easier as you play with it. Hopefully I can get you started! What language program, btw?

 

[Maltman]

On 6/15/2017 at 2:21 PM, Cunjo Carl said:

I'm in! My free time is pretty sparse, sadly. Do you have availability tonight 9-10pm PDT (4-5am GMT)? If not, I can try to swing a time this weekend.

Just like anything, it's tricky stuff when you first get going, but gets easier as you play with it. Hopefully I can get you started! What language program, btw?

 

 

My free time is infinite thankfully! So we should be able to talk some time at night or during the weekends. I'll send you a PM with like my gmail and you can contact me on hangouts is that OK?

Yeah I struggle with the feeling of being overwhelmed a lot, but I'm learning to ask others for help! I'm making it in Unity with C#.