Orbital mechanics for realistic space flight simulators

This guide covers the math required to calculate a wide range of orbital maneuvers, including:

1. AP and PE
   
2. Changing AP and PE
   
3. Changing planes to align with a target
   
4. Performing a phasing burn which aligns PE between two crafts
   
5. Interplanatary transfers

In addition, a craft file is available which contains Vizzy implementations of the math shown.

Introduction

There is a great deal of material available on the Internet about orbital mechanics, and most of it is even comprehendible by a mathematically oriented lay-man.

However, as far as I can tell, the “practical examples” that these resources provide are less then helpful if you are attempting to automate real-world common orbital mechanics tasks. Questions like “How do I setup phasing for a rendezvous” or even “How do I perform a rendezvous” are simply not covered. The reason is seems obvious to me – these sorts of problems are not something that could appear on an exam in college, simply because the math involved is too complex. Since the all the sources that I have found are aimed at college students who are taking a class in orbital mechanics there is a gap and this guide is intended to fill that gap.

Therefore, while this guide includes a list of formulas that will be used, I do not cover the derivation of those formulas unless I feel the existing sources are inadequate in some way. I do spend quite a bit of time describing algorithms to perform tasks that are useful for someone trying to automate orbital maneuvers in a realistic space flight simulator, such as Juno: New Origins, Kerbal Space Program (1 or 2), or Orbiter[orbit.medphys.ucl.ac.uk].

Status

This section will contain information on updates to the guide, assuming there are any.

Note that no changes are expected to occur, except as a result of feedback. While there are still things that could be added, I am very unlikely to have either the time nor the motivation to add them.

9/26/2023: Initial release

Pre-requisites

You should have a very good grounding in algebra, know the basics of trigonometric functions, know the basics of vectors, and know what a reference frame is and why they are useful (in the context of vectors), and the very, very basics of calculus.

I would expect most American high school seniors who have taken trig and introduction to calculus to be able to follow along with this guide. Someone who has taken high school physics may be OK – certainly AP Physics would be enough.

Audience

While this guide is published on Steam under “Juno: New Origins,” this guide has very, very little that is specific to Juno, and all Juno specific material is marked as such (almost always, by being placed in “quotation” boxes). In my opinion Juno is likely the best existing platform (at the time this guide was published) to implement this guide, by virtue of being the only realistic flight simulator that also includes an integrated programming environment, but this guide could also be implemented in kOS (a mod for KSP 1 that adds a programming environment) or as standalone mods for KSP 1 or 2 or Orbiter.

Still, if you are developing in Juno, I have a complimentary guide that discusses my specific implementation of everything covered in this guide:

https://steamcommunity.com/sharedfiles/filedetails/?id=3039843312
Which includes a downloadable craft file that implements everything.

References

• The vast majority of the math used in this guide came from this page[orbital-mechanics.space]. I cannot recommend this page enough if you want to understand the derivation of the formulas used in this guide.
  
• I referenced this[www.anilvrao.com] PDF file when I needed more detail than the first page provided. Most notably, I used this PDF to determine the unit vectors for the perifocal reference frame and the formula to determine the line on which the ascending / descending nodes reside to align planes.
  
• In addition to the above, I also referenced a number of other sources for specific details. These references will be called out as they come up in the guide.

On the subject of circular orbits

As far as this guide is concerned, they do not exist. 😊

Well, they could exist, but they are very rare whether you are talking about in a game or in real-life. Keep in mind that I am defining a circular orbit in the mathematical sense – one that has an exactly constant distance from the center of the planet throughout the entire orbit. If an orbit has an apoapsis of 4,000,000.01 and a periapsis of 4,000,000.0 it is not circular as far as this guide is concerned.

In Real Life™, the odds of lucking into an exactly circular orbit are nearly zero due to errors introduced by engine variability and sensing technology. Even if you did end up in such an orbit, the odds are good that you would not realize it, because the sensing instruments that you are using to establish the velocity and position of the craft are insufficiently precise.

In the context of a realistic space craft simulator, floating point errors make it pretty much impossible to achieve a perfectly circular orbit without resorting to debug commands – and, even if you used debug commands to achieve the desired orbit, as soon as you allowed time to advance errors would be introduced (see the next section for why) that would make it no longer circular.

Circular orbits have several unique problems, all related to the fact that there is no well-defined periapsis position. As a result, many of the equations and tasks described in this guide would require special case logic to handle circular orbits properly, and… Given how rare they are, and how complex this guide is already, I decided to omit them. Feel free to add support for circular orbits if you wish – I will not stop you. 😊

Parabolic orbits are even more rare, as few players will think to use debug commands to setup a parabolic orbit and second most orbits that look parabolic are in fact elliptical – they just exit the current sphere of influence before apoapsis. However, parabolic orbits need much less special handling than circular orbits (they only require a unique formula to convert from time to true anomaly and vice versa), so I do include support for them.

Accuracy

The formulas in this guide are 100% correct and therefore produce 100% correct results.

However…

Realistic space flight games such Kerbal Space Program, Orbiter, and Juno: New Origins update the positions and velocities of crafts using two distinctly different methods:

1. “Numerical Integration:” In this mode, the game:
   1. Sum up all the forces acting on an object, with the most common examples being gravity, thrust, drag, and lift,
      
   2. Divides the net force by the mass of the craft to produce an acceleration vector,
      
   3. Multiplies the acceleration vector by a small increment of time, which I will call a “physics frame tick,”
      
   4. Adds the vector from the previous step to the current velocity vector,
      
   5. Multiply the new velocity by the “physics frame tick” length,
      
   6. Add the vector from the previous step to the current position vector.
2. “Analytical Integration:” In this mode, the formulas used in this guide are used to produce the new velocity and position vectors after the passage of an arbitrary amount of time. Computationally, it costs the same amount to produce new velocity and position vectors 0.01 seconds in the future as it does to produce a prediction 1,000 years in the future. Certain operations take much longer (specifically, going from time to position or velocity is much more expensive than the reverse), but even then, the cost does not scale with the increment at all. This is commonly referred to as being “on rails,” and is how these games support high levels of time compression.

In Juno: New Origins, the numerical integrator is used at the 1x and 2x time compression, with all higher time compressions using the analytical integrator. Other programs may have different behavior, and some allow you to configure when the switch occurs.

Numerical integration is, fundamentally, inaccurate. It assumes that both the force (both magnitude and direction) acting on an object and the mass of the object will remain constant over some small period of time (in Juno, on my computer, typically ~0.0167 seconds). This is not the case for gravity (as the movement of the object will change the distance from the center of mass, plus the direction of the force will vary due to the object’s movement) nor even rocket engines (as the mass of the rocket drops as thrust is produced, so thrust performed now produces less acceleration than thrust produces in the future).

One might think “Hey, couldn’t you use an integral to fix this problem?” and… Yes, you can, and what you produce is the equations in this guide for simple (single force) situations, and there is no analytical solution for the general N-Force case. So… Yeah, simulations will be inaccurate when running in numerical integration mode.

How inaccurate? For that, refer to my other guide which contains detailed examples showing the magnitude of the problem. If you are on a different platform the problem may be worse or it may be better – but it will inevitably be present to some degree, because numerical integration simply cannot be as accurate as analytical integration.

What can you do about this? There is only one solution – do not allow the numeric integrator to run at all. You can achieve this by setting the time warp to a high enough value (in Juno, 10x) to force the use of the analytical integrator immediately after each burn, and leaving it at that setting (or higher) until the next burn. You do have to drop to 1x before the burn (to adjust the craft’s attitude to the burn’s attitude) and you need to leave it at 1x during the burn, but at the end of the burn reset the time warp before you calculate the new orbit. This will reduce the error to the maximum degree possible -- and it is significant enough to improve burn accuracy by a noticeable amount.

The downside of this is that you have less time real time to perform calculations between burns, but… I consider this a fair trade. Still, your mileage may vary.

Another, vastly inferior in my opinion, is to recalculate your maneuver continuously until it is time to begin the burn. I do not think it would be possible to split your “Calculate maneuver” code from your “Perform burn” code, in this case, but… It is not required that these be separate entities.

How is the problem addressed in the Real World™ (say, at NASA)? I do not know for sure, but I suspect that there are two ways that they deal with this issue:

1. First and foremost, their simulations do not need to run in real-time, so they can set the increment to a very, very small number (say, 0.0000000001). The smaller the increment the more accurate the result of the numerical integration, so this helps.
   
2. On the flip side, there are significant errors in the knowledge of the actual position and velocity of the craft due to a wide range of factors (including measurement error, changes in the center of mass due to fuel slosh, variability in the performance of engines, routine venting operations, and much more) so they simply expect their orbital predictions to be incorrect and plan for mid-course corrections to fix any errors that crop up.

Conventions

• In text, a vector will always be bolded (for example “v”).
  
• In an equation, a vector will be bold, in italics, and have an arrow over it. For example, (“”).
  
• In an equation, a unit vector (a vector whose length has been set to 1) will be bold, in italics and have a caret (also known as a “hat”) over it. For example, (“”).
  
• A dot above a variable in an equation, such as “,” represents the derivative of r.
  
• represents the operation “Find the magnitude of the vector v.”
  
• represents the operation “Normalize vector v.”
  
• When you see a scalar variable (in italics in equations) that shares the letter used for a vector (in bold and italics, with an arrow over it), the scalar variable is the magnitude of the matching vector. Thus, .

General definations

• are the three unit vectors for the planet centered inertial reference frame (PCI). In any sensible platform, these will be defined as (1, 0, 0), (0, 1, 0), and (0, 0, 1). By convention, the K unit vector starts at the center of the planet and points at the north pole on some specific date (the “start of the epoch”), the I vector starts at the center of the planet and points at the point where the equator crosses the prime meridian on the same time and date, and the J vector is the cross product of the other two unit vectors.
  

  The above is what I would expect from realistic space flight games. The actual definition is quite a bit more complex[en.wikipedia.org].

  For reasons that I cannot even speculate on, Juno defines the J vector as the one that points from the center of the planet to the north pole. This causes numerous equations to need to be adjusted to take this into account, which will be called out when it occurs.

• G = the universal gravitational constant (~6.6743 × 10-11)
  VERY, VERY, VERY IMPORTANT: G is a constant – but there is considerable debate as to exactly what the value is. When this guide was written, I found the following values:
  • Wikipedia[en.wikipedia.org] says 6.674 × 10-11,
    
  • Swinburne University[astronomy.swin.edu.au] says 6.67 × 10-11,
    
  • this[byjus.com] site (which is the first Google response as of this writing) says 6.67408 × 10-11,
    
  • The Britannica site[www.britannica.com] (who used to be well known as an encyclopedias manufacturer) says 6.6743 × 10-11, and
    
  • NASA[imagine.gsfc.nasa.gov] says 6.672 × 10-11.
  This is an unusually variable constant… 😊
  
  Using a value that is different than what your platform uses… “Bad.” Really, really bad, because it makes all of your future position and velocity predictions wrong, but only by a small amount. Do not do that.
  
  For the record, Juno: New Origins uses 6.67384 × 10-11, which does not even appear on the first page of the Google results. If you search for the number directly, this[cosmosmagazine.com] comes up – which discusses an experiment in 2012 that produced a different value, but mentions 6.67384 × 10-11 as the “official value.” Beyond this site, none of the Google results appear to be even slightly authoritative, so my guess is that the Juno developers set “G” back before 2012, and science has marched on further refining the number.
  
  Given the above, well designed programming platforms should make this value available as a constant.
  
• m = the mass of the active craft.
  
• Apoapsis = The maximum distance from the center of the planet being orbited and the orbiting object. In equations, this is commonly referred to as “AP.”
  
• Periapsis = The minimum distance between the center of the planet being orbited and the orbiting object. In equations, this is commonly referred to as “PE.”
  
• = the position of the object, in an inertial reference frame whose origin is at the center of the planet being orbited, as measured at some point in time (as noted in the subscript).
  
• = the velocity of the object, using the same reference frame as the position as measured at some point in time.
  
• State vector = both the position and velocity of an object measured at the same time, using the same reference frame. This is sufficient to fully determine the orbit of the object regardless of when you measure the two values.
  
• h = the angular momentum vector for an orbiting object. This fully defines the orbital plane of the object, and therefore corresponds directly to inclination.
  
• h = the magnitude of the angular momentum vector. This value can also be calculated in a number of other ways.
  
• e = the eccentricity vector. This vector points towards the PCI position of the craft at periapsis and defines where true anomaly = 0.
  
• e = the magnitude of the eccentricity vector, normally referred to as “eccentricity.” This value can be calculated from the apoapsis and periapsis of an orbit.
  
• a = the semi-major axis of an ellipse. This value can be calculated from the apoapsis and periapsis of an orbit. Note that if the orbit is circular (e=0), this will be the radius of the orbit.
  
• T = the period of an orbit. This is the amount of time it takes the object travel from periapsis back to periapsis in an elliptical or circular orbit. This value can be calculated from the apoapsis and periapsis of an orbit.
  
• θ = true anomaly. This is the value of the angle between two vectors:
  • The first pointing from the center of the object being orbited towards the position of the object at periapsis. This is the e vector.
    
  • The second pointing from the center of the object being towards the position of the object at some point of interest. This is the r (position) vector.
  Important: The standard symbol of this value is Greek letter Nu (“ν”). Yes, that glyph is different from a lower-case “v” – look very carefully. As you might imagine, using this glyph in the same context that you are referring to velocity can make it very difficult to determine which is which.
  

  After checking, the font used for Steam guides appears to render ν and v with the exact same glyph, so they really are identical. Try copying the characters into WordPad (in Windows) to see the difference.

  Some authors, recognizing the issue, use “f” for this value instead.
  
  Since this is my guide, and since I do not use θ for any other purpose, I decided to use θ to represent this angle.
  
• M = mean anomaly.
  
• Ε = Eccentric anomaly.
  This is the Greek letter Epsilon (“Ε”) rather than the letter “E.” Since E isn’t used for any other purpose in this guide, there is no ambiguity, so I left this definition alone.

If you want a full definition for mean anomaly and eccentric anomaly, please refer to this[orbital-mechanics.space] page. In the context of this guide, these can be defined as “intermediate values that you need to calculate to convert a true anomaly value into time since periapsis.”

Basic formulas

Formula for Mu (μ):


This term is just a convenient shorthand used in orbital mechanics to simplify the equations a bit. This value is constant for all orbits around the same central body.

Formula for the angular momentum vector (h):

The above formula is the “standard” formula for angular momentum, the one that you will find in all the available online sources. It assumes that the cross product is defined using the “right-hand rule.” Unity games (including Juno) use “left-handed” cross products. There are three options here, although two are very similar:

1. You can reverse the order in which r and v appear in the above equation – this will produce the correct angular momentum vector, but the formula itself will appear “wrong” if someone is checking against online sources.
   
2. You can negate the result of the above equation. This has the same effect as #1.
   
3. You can leave the formula as is, accepting the fact that the h vector will be pointing in the wrong direction (specifically, it will be pointing “down” when it should be pointing “up”).

Unless you are performing operations directly on the angular momentum vector, the direction does not matter. For myself, I choose to option #3 (leave the formula “correct” and accept the incorrect direction), and the rest of this guide will assume that you do the same.

Formula for the eccentricity vector of an orbit (e):

If you do not correct the order of the cross product for the angular momentum (the previous equation), the “left handedness” of this cross product cancels out the direction error introduced in the angular momentum equation, leaving an equation that both matches the standard form and produces the correct result.

If you reverse the operands for the cross product in the angular momentum equation, however, you will also need to reverse the operands here as well. If you do not, the eccentricity vector will point towards apoapsis instead of periapsis as is intended – and this will cause problems down the line.

If we know apoapsis and periapsis, we can calculate the magnitude of the eccentricity vector using these values instead:


Formula for the semi-major axis (a):


With a bit of algebra, we can solve this equation for h, providing an alternate formula for the magnitude of the angular momentum vector:


Due to the way in which h was defined earlier (as a magnitude of a vector), h will always be positive – the negative result from the square root can be ignored.

If you know apoapsis and periapsis, you can calculate the semi-major axis (a) directly:


The “orbit equation” (and yes, that is this formula's official name):


The more mathematically inclined readers will recognize this as an equation in polar coordinates.

When you solve this equation, there will be two inputs (θ) that produce the same r value. This is called “quadrant ambiguity” and is caused by the fact that all trig functions are cyclic with a period of π radians (although around which value they cycle varies with function – sine and tangent cycle between [-π/2…+π/2], while cosine cycles between [0...+π]. Quadrant ambiguity will cause problems throughout this guide.

The "orbit equation" can also be solved for θ, which is called the “true anomaly”:


When you use this formula, there will be two solutions, the one returned and another that is 2π-x, where x is the value returned from the equation. Which solution you are interested in will depend on context, and it requires more information than this formula has available to it. There is much more on this topic later, when we do have enough information to resolve the ambiguity.

The “parameter of the orbit” or "semi-latus rectum":

I swear, I am not making up these names. These really are the official names that can be found in reputable sources[orbital-mechanics.space].

This can be substituted back into the orbit equation, producing:


Since h is constant for a particular orbit, and μ is a constant for all orbits around a specific central body this both simplifies equations and saves CPU cycles.

Perifocal reference frame

I assume here that you are aware of what a reference frame is and how you convert vectors from one to another. If this is not true, please refer to my other guide:
https://steamcommunity.com/sharedfiles/filedetails/?id=2944674093
Especially the section entitled “Redefining axes for fun and profit.”

Unlike the reference frames defined in my previous guide, the perifocal reference frame is inertial – it does not change as the position of the active craft changes. Therefore, there is no need to recalculate the unit vectors for each physics tick, which is a very good thing. You must, of course, must recalculate the unit vectors if the orbit changes (due, for example, to a burn).
The perifocal reference frame looks like this (for an elliptical orbit):


It is important to note that the planet does not lie at the center of the ellipse but instead one of the two focal points (the “occupied focal point”). By design, the occupied focal point will always be the one on the right, placing periapsis on the right-hand side as well.
The unit vectors are defined as follows:

• +X = the direction from the center of the planet to the position of the object at periapsis. The unit vector for this direction is defined as:
  
  
• +Z (not shown) = one of the two vectors perpendicular (aka “normal”) to the orbital plane. The unit vector for this direction is defined as:
  
  

  If the cross product is defined using the "left-handed rule" and you did not fix the definition of h to account for this, you will need to invert this vector, like this:
  
  If you do not do this, the math will still work, but true anomaly (θ) will decrease with the passage of time and the craft will rotate in a clockwise direction. This will make it more difficult to follow the rest of this guide, so... Don't do that.

• +Y = the cross product between the +X and +Z vectors. This will be the vector that is tangent to the orbital path at periapsis. The unit vector for this direction is defined as:
  

The perifocal reference frame has a wide range of highly desirable characteristics that make its use all but mandatory when performing orbital mechanics:

• First and foremost, all orbital motion will occur within the perifocal X / Y plane – practically speaking, the Z component can be ignored altogether, turning a complex 3-dimensional problem into a much more tractable 2-dimensional problem.
  
• Second, regardless of whether the craft is in a retrograde or a prograde orbit, the craft always moves in a counter-clockwise direction in the perifocal reference frame.

Finally, the 2 dimensional Cartesian position of the craft can be very easily converted into a 2 dimensional polar form, which is even easier to work with. To do this, we describe the craft’s position in terms of an angle with the X axis (θ) and distance from the center of the planet (r). -- that is, as (θ, r). Geometrically, this looks like this:


Mathematically, we can calculate true anomaly like this:


However, this will only return values in the range [-π/2…+ π/2], leading to quadrant ambiguity. For many, but not all, platforms a function "atan2" exists which is specifically designed to address this problem. "atan2" takes two parameters (instead of 1), and those parameters are the numerator and denominator of the values that you would normally pass to arctangent. In terms of "atan2," we can define true anomaly this way:


However, atan2 returns values from [-π…+π] rather than the preferred [0…+2π], so it will be necessary to check the to see if the value returned is negative, and if it is add 2π

It is important to note that the definition of true anomaly does not care about the magnitude of the input vector, only its direction. This means that all points on the perifocal plane have true anomaly values associated with them, regardless of whether or not the points are actually part of the orbit. This insight makes certain math later much, much easier.

r is defined as the magnitude of the original vector, which is a simple formula:


Converting to a polar coordinate system has one additional benefit: we only need worry about one value in 99.99% of the cases, specifically true anomaly. Almost all of the equations from this point on will use true anomaly as a proxy for “The position of the craft.”

Note that the craft’s true anomaly value will always be 0 radians (0°) at periapsis, and π radians (180°) at apoapsis. Finally, true anomaly will always increase as the craft moves along its orbit, even if the craft’s orbit is retrograde.

We can now define use the orbital equation to define the position of the craft as a Cartesian vector in the perifocal reference frame:


We can find the derivative of this formula to find velocity vector of the craft as a Cartesian vector in the perifocal reference frame:

Conversion to and from the perifocal reference frame via unit vectors

The easiest way to convert a PCI vector to the perifocal reference frame is via the unit vectors defined above (p, q, and w). These three vectors are, themselves, PCI vectors and can be used with the dot product to convert an arbitrary vector v in PCI coordinates to perifocal coordinates. This looks like this:


The above formula can be used to convert the PCI unit vectors (I, J, and K) into perifocal unit vectors like this:


With the the PCI unit vectors converted to perifocal, we can now convert any perifocal vector into a PCI vector with this formula:


This should work. However, via experimentation in Juno, these conversions produce a small error -- not much, but enough to create significant errors when the predicted position of the craft is compared with the actual position of the craft.

Conversion to and from the perifocal reference frame via rotations (1/2)

An alternate way to convert vectors to and from the perifocal reference frame is via rotations. In order to do this, we need to define three angles of rotation. In orbital mechanics, all three have special names:

Inclination (i), which is the angle between the PCI axis that passes through the poles (normally, the PCI Z unit vector, (0, 0, 1), also known as K) and the angular momentum vector (h). The formula for this is:

While the result of this equation is restricted to the range [0… π] this is not only acceptable but required for inclination. Therefore, we do not need to deal with quadrant ambiguity.

The alert among you will have realized that this formula will not work properly in Juno – in this game, for unknown reasons, the Y axis is the one that passes through the poles. So, the correct equation in Juno is:

Due to the way that cosine works, this function returns the correct results without modification regardless of whether you calculated h using left-handed cross product or the right-handed cross product.

The “Right Ascension of the Ascending Node” (Ω), which is the angle between the PCI X axis and the point at which the orbit passes through the equatorial plane (the plane formed by the PCI unit vectors in the X and Y directions) when the Z component of the position transitions from a negative value to a positive value – that is, it is “ascending” from the perspective of a viewer who is viewing the equatorial plane “edge on.”

In order to calculate this value, we first need to find the line that is the intersection of the equatorial plane and the orbital plane. This formula will produce will work:


Next, we can use N to find desired angle:

When N is the zero vector, this will be undefined. By convention, Ω = 0 in this case (source[en.wikipedia.org]).

As with the previous use of arc-cosine, the range will be [0…π] and we do need values from [0…2π] for this parameter. This can be done with a piece-wise definition:

Again, this will not work properly in Juno, both because K is not the vector that points through the north pole and because Unity games use left-handed cross products instead of right-handed. Thankfully, if you use the correct equation for h (rather than modifying the equation to produce the correct result) then the error introduced by the change in handiness cancel out, and the only change to the definition of N is using the correct unit vector:


Where J is the unit vector (0, 1, 0). We also need to alter the piece-wise definition for Ω to the following:

The “Argument of Periapsis” (ω), which is angle from the ascending node (which we worked with in the previous formula) to the periapsis. Since both the ascending node and the periapsis are points on the orbit, this angle will be measured in the orbital plane, and an alternate definition would be “The true anomaly at the equatorial ascending node.” The formula for this is:


Alternatively, this could be defined this way:


Where “angle” refers to a function that returns the angle between two vectors.

This returns a value between [0…π], so we resort to a piece-wise definition again:

Yet again, we need to modify this for Juno due to the Y and Z axes being flipped. The correct equation for Juno is:

When N is the zero vector, ω is undefined. By convention, the following formula should be used if e is not zero (source[en.wikipedia.org]):


If the orbital is circular, e = 0 and again ω is undefined. In this case, the convention is to set ω = 0 (source[en.wikipedia.org]).

Conversion to and from the perifocal reference frame via rotations (2/2)

Using the three values given, we can now rotate a vector from the perifocal reference frame to the PCI reference frame as follows:

1. Rotate the input vector around the (0, 0, 1) axis by -ω.
   
2. Rotate the vector from the previous step around the (1, 0, 0) axis by -i.
   
3. Rotate the vector from the previous step around the (0, 0, 1) axis by -Ω.

As you might imagine, the above will not work if you are following along in this guide (and therefore “h” is pointing in the wrong direction, but you inverted the sign when defining the perifocal reference frame). This is an easy fix –the first (?) reference to (0, 0, 1) should be replaced by (0, 0, -1).

In addition, the above will not work in Juno (specifically), because the PCI Y and Z axis are swapped. This is much harder to fix, and I never worked out the proper solution. I did come up with a very hacky solution that involved calculating the rotations ignoring the issue altogether, then swapping the Y and Z components. After I implemented matrix rotations, I found a proper solution, but I am not sure how to translate it into three single axis rotations.

We can use Rodrigues' rotation formula[en.wikipedia.org] to rotate a vector around an arbitrary axis (it is not required to be a unit axis, but it does need to be normalized (magnitude = 1) ):


Where v is the vector to be rotated, k is the normalized vector to be rotated around, and ∅ is the amount of rotation desired (in radians, of course).

Performing reference frame conversions via this formula, however, will be both much slower and much less accurate than using the unit vectors. There is a very specific case where this is worthwhile -- namely, if you have an orbit defined in terms of apoapsis, periapsis, i, Ω, and ω, then you can use Rodrigues' rotation formula to convert the perifocal p and w unit vectors into normalized PCI vectors that will point in the correct direction for the vectors e and h. As such, orbits are routinely published in this format and this sort of conversion is necessary if you need to deal with the orbits in PCI reference frame.

In the context of realistic space flight simulators, it is very common for the orbits of planets and moons to be expressed in the "apoapsis, periapsis, i, Ω, ω" in the game definition files. As a result of this, using these parameters to generate the other orbital parameters may be produce better results than using the planet's or moon's state vector to do the same.

It does in Juno, at least.

We can dramatically improve performance by defining matrixes to perform the rotation, then pre-calculating the matrixes.

For perifocal to PCI conversion, the correct matrix is:


And to convert from PCI to perifocal, the following equation should be used:


It is not necessary to reverse the sign on the three angles for this matrix to work properly. The formulas swap cosine for sine as required to achieve the same effect.

For those who have forgotten, or never knew, this page[www.mathsisfun.com] describes how you perform matrix multiplication. The perifocal to PCI conversion matrix can be converted to this:


If that still does not seem “simpler” to you, then try this version on for size:


Ahaha – you see, all we need to do is multiply each component of the input vector by a scalar number, and we can re-use the rotation matrix as many times as we wish, so it is very fast even without explicit machine language support.

Ehhh… Your milage may vary.

The eventual reason that I ended up implementing rotations in addition to unit vectors is to compare the results – and it turns out that rotations produce better (more accurate) results. I believe that this is, in large part, because the game engine I was developing on (Juno) uses rotations internally, but it may be that the dot products between the unit vectors and the input vectors introduce significant error.

“Significant error” means “Several 100s of meters of difference in actual vs predicted position in a 3000 km x 100 km orbit at apoapsis. In orbital mechanics, this level of error is far, far too high to be accepted if you are trying to make accurate burns.

The above matrixes will not work in Juno due to the X and Y values needing to be swapped. These will, however:

The change here is swapping the second and third row, plus reversing the sign on the value that ended up in row 2, column 2.

Here, all that was required was swapping the second and third columns. No sign reversal this time, though.

Mean anomaly and time

If you wish to understand the derivation of these equations, I strongly recommend this page[orbital-mechanics.space].

On the other hand, it is not necessary to fully understand the derivation to make use of the formulas, so... Up to you.

If the orbit is elliptical (0 < e < 1), mean anomaly can be calculated from true anomaly (“θ”) and eccentricity (“e”) as follows:


Alternatively, you can define mean anomaly for elliptical orbits in terms of its period (“T”) and time since periapsis (“t”):

The above equation cannot be derived from the above equation – it has a separate derivation.

Where the period of the orbit (“T”) is defined as:


Which can be dramatically simplified by using the semi-major axis “a” instead of the angular momentum “h”:

Only elliptical and circular orbits have a period (“T”), thus the lack of a subscript for this term in the above equations.

If the orbit is parabolic (e=1), then mean anomaly in terms of true anomaly (“θ”) can be calculated with this formula (e will always be 1, of course, and therefore does not appear in the equation):


In terms of time since periapsis (“t”) you can use this equation:

Again, this equation is not derived from the previous equation.

If you know the mean anomaly (from the above equation), you can solve this equation for t:


If the orbit is hyperbolic (e>1), then mean anomaly can be calculated from true anomaly (“θ”) and eccentricity (“e”) as follows:


An alternate equation expresses mean anomaly in terms of time since periapsis (“t”):


This can again be dramatically simplified by using "a" instead of "h" and "e" as follows:

"a" will always be negative for hyperbolic orbits -- thus, the negative sign in the above equation.

The value of M increases to infinity as time increases and does so very quickly when a is small. This is expected and, while it makes the mean anomaly useless in isolation, it works for the purposes of the following formulas. Do not try to be clever and normalize the value of M (to fall between -𝜋 and 𝜋, for example). This will absolutely prevent the other equations from returning the correct values.

Eccentric anomaly and mean anomaly

For elliptical orbits (0 < e < 1), eccentric anomaly can be expressed in terms of eccentricity (“e”) and true anomaly (“θ”) with the following equation:


You can solve this equation for eccentric anomaly (“Ε”):


Or you can solve for true anomaly (“θ”):


Eccentric anomaly can be related to mean anomaly with the following equation (this is also known as “Kepler’s Equation”):


Note that this equation cannot be solved analytically (via algebra) to express “Ε” in terms of “M.” More on this in a bit.

As far as I can tell, nobody has bothered to even define eccentric anomaly for parabolic orbits (e = 1). So… Yeah. I am confident that the reason for this odd omission is that nobody is actually interested in eccentric anomaly in its own right – its value is as an intermediate stop between “time since periapsis” and “true anomaly.” As we will see later, this step is not required for parabolic orbits, so why define it?

For hyperbolic orbits (e > 1), eccentric anomaly can be calculated in terms of eccentricity (“e”) and true anomaly with the following equation:

By convention, this eccentric anomaly for a hyperbola is referred to as “F.” I have no idea why, so I am not going to do that. 😊

This can be solved for eccentric anomaly:


Or for true anomaly:

If your development environment does not implement hyperbolic trigonometric functions, such as “tanh,” never fear – these can be defined in a finite (and shockingly simple at that) set of standard algebra functions. See here[en.wikipedia.org] for more details.

It can also be expressed in terms of mean anomaly with the following equation (“Kepler’s Equation for the hyperbola”):


This also cannot be solved for eccentric anomaly analytically.

To convert time since periapsis (t) to true anomaly (θ)

For parabolas, this is easier as we can skip eccentric anomaly altogether. The process is:

1. Calculate the mean anomaly from t using this formula:
   
   
2. Evaluate the following formula:
   
   Where:
   

For elliptical orbits and hyperbolic orbits, the process is much more complex. The overview is simple enough:

1. Calculate the mean anomaly from t using the appropriate formula.
   For an ellipse, this formula will be:
   
   
   For a hyperbola, this formula will be:
   
   
2. Calculate the eccentric anomaly from the mean anomaly from mean anomaly using the appropriate formula.
   
3. Finally, calculate the true anomaly from the eccentric anomaly using the appropriate formula.
   For an ellipse, this formula will be:
   
   
   For a hyperbola, this formula will be:
   

No analytical solution exists for step 2, which means that no matter what you do, you cannot manipulate this equation algebraically to produce a single Εₑ term. Therefore, it is necessary to use an iterative process to calculate this value. One such method is known as Newton’s method for finding the roots of polynomials[en.wikipedia.org].

Implementing step #2 for an ellipse looks like this:

1. Put all terms of the equation to be solved on the same side:
   
   
2. Find the derivative of the equation in step 1:
   
   
3. Make an initial guess (x₀) as to the correct answer. The mean anomaly (M) works well here, and we just calculated it.
   
4. Evaluate:
   
   Until the error term
   
   is small enough.

For hyperbolas, step #2 looks like this:

1. Put all terms of the equation to be solved on the same side:
   
   
2. Find the derivative of the equation in step 1:
   
   
3. Make an initial guess (x₀) as to the correct answer. This equation works well enough:
   
   

   The “0 when M = 0” is absolutely required – logarithms are not defined at zero.

   If you wished, you could improve the initial estimate slightly by checking to see if M is between -1 and +1. In such cases, you should multiply the entire equation by -1. The goal of the sign function here is to ensure that the result of the logarithm has the same sign as M and the logarithm of a fractional number is negative. I suspect that the extra cycles spent implementing this logic would exceed the number of cycles saved by reducing the number of iterations required, however – the goal is to simply get close, and the extra error introduced by the incorrect sign in the initial guess is very, very small.

4. Evaluate:
   
   Until the error term
   
   is small enough.

How small is "small enough"? That is up to you, but if you want an error < 0.01°, the error term should be less than < ~ 0.0001745.

Thankfully, you only need to worry about this conversion when you are attempting to rendezvous two orbiting objects (and, even then, only when you are working out the orbital parameters to bring the two objects to the same point at the same time). If all you are doing is adjusting an orbit to specified values, you do not need to convert times into true anomaly values.

Converting from true anomaly (θ) to time since periapsis (t)

For parabolas, this is once again easier:

1. Calculate the mean anomaly from the true anomaly value using this formula:
   
   
2. Calculate time from this formula:
   

For elliptical and hyperbolic orbits, this process is a bit more complicated:

1. Calculate the eccentric anomaly from the true anomaly with the proper equation.
   The correct formula for elliptical orbits is:
   
   
   The correct formula for hyperbolic orbits is:
   
   
2. Calculate the mean anomaly from eccentric anomaly with the proper equation.
   The correct formula for elliptical orbits:
   
   
   The correct formula for hyperbolic orbits:
   
   
3. Calculate the time from periapsis with the proper equation.
   The correct formula for elliptical orbits:
   
   

   For this formula to return the correct results, M must be in the range [0...2π] but the previous equation will produce a result in the range [-π…+π]. Therefore, you must check the result and if it is negative add 2π before applying this equation.

   The correct formula for hyperbolic orbits:
   

Practical examples

The remainder of this guide will cover practical examples of using the equations described up to this point to calculate the required burns to perform a variety of common orbital mechanics tasks. There will be much less math, although there will still be some.

Performing burns (1/4)

For this section, we define, and in some cases redefine certain variables:

• F = The total thrust of all engines in the active stage
  
• m₀ = Initial mass of the spacecraft (including fuel)
  
• Iₛₚ = ISP of the spacecraft
  
• g₀ = the value of gravity that is used to calculate ISP – this is always ~ 9.80665, the gravity on Earth at sea level.
  
• e = Euler's number, the base of the natural logarithm function – this is always ~ 2.71828.

Propellant mass required by Delta-V and ISP

We can generate a formula for the change in mass as a function of force generated (aka “Propellant used per second of operation,” aka “The first derivative of mass with respect to time”):


We can then solve the rocket equation (here[en.wikipedia.org]) for the final mass of the craft (after the burn):


Next, we can calculate how much fuel must be burned to produce a specified change in velocity:


Finally, we can define the duration of the burn as “The amount of fuel that will be burned during the burn divided by the amount of fuel consumed per second”):


This particular formulation of this equation comes from here[space.stackexchange.com], and I find it delightfully simple and easy to understand.

To perform burns accurately, we need to know the mid-point of the burn – that is, the point at which half of the velocity change has occurred. You would think that you could just use the above formula with one half of the delta, but that is not correct (it is a reasonable approximation, however). The inaccuracy is due to a failure to account for the increased rate at which the engine can convert propellent into delta-v during the second half of the burn. This is because the craft will be lighter, of course, as fuel is burned off.

A kOS script that is reported to be working (here[github.com], lines 151 and 152) uses these two formulas:


Where:


These equations can be simplified considerably. Starting with the first equation:


This matches the equation that we derived earlier (once you substitute in the formula for dot-m).

The second equation simplifies like this:


This formula tells us how when the midpoint of the burn will be.

Performing Burns (2/4)

In general, burns will be scheduled to occur at a specified true anomaly (rather than a time) and the desired burn will be a vector. The basic process for performing a burn is:

1. Calculate the time from the craft’s current position to the true anomaly of the burn. If the result is negative, then add the period of the current orbit to make it positive.
   

   This will not work for parabolic or hyperbolic trajectories, as such orbits have no period to begin with. Furthermore, while this will always work for elliptical and circular orbits, if the purpose of the burn is to adjust phasing to achieve a rendezvous, adding the period of the orbit will totally invalidate the calculated intercept, so the code that generates true anomaly values for intercepts should always ensure the generated true anomaly value occurs in the future.

2. Calculated the midpoint of the burn and the duration of the burn using the appropriate formulas given above.
   
3. Calculate the burn start time by subtracting the time for the active craft to reach the burn true anomaly (as calculated in step #1) and the calculated duration to the burn midpoint (from step #2).
   
4. If the burn duration is very short (say, less than 5 seconds), reduce the starting throttle and try again – the minimum burn duration should be greater than 10 seconds.
   

   How important this is very platform dependent. If your platform has engines that respond instantly to the throttle then this can be skipped – otherwise, extending the duration of the burn will improve accuracy due to the reducing the amount of error introduced by the “spool up” and “spool down” periods (where the engine is commanded to produce a certain level of thrust, but has not yet reached the desired value). The unexpected low thrust at the beginning of the burn will mean the start time is slightly later than it should be, and the spool down thrust will result in the burn lasting longer than expected. While a lower throttle setting does not eliminate these issues completely, it does reduce the relative magnitude significantly.

5. Wait until you are “close” to the expected start time. How close depends on how maneuverable the craft is and the quality of the autopilot – you want to leave enough time to point the craft in the correct direction before the scheduled start time.
   
6. Align the craft with the desired delta-v vector.
   

   This is platform dependent. In Juno: New Origins, the autopilot expects a PCI vector as input, and the actual direction that the craft points for a given PCI vector relative to the prograde / normal / radial directions varies with the position of the craft in the orbit. You would think that this would produce an incorrect result (as start out pointing in one direction, then gradually change the attitude as the burn proceeds, only reaching the “correct” value once the craft arrives at the mid-point of the burn). You would be wrong – I worked out a way to point the craft in the “correct” direction at the start of the burn, and the errors were not only much larger but varied based on attribute (in particular, the e value was way off).

7. Start the burn at the scheduled time.
   
8. Stop the burn when some condition is met. The most obvious condition is “the expected duration has passed.”

Step #7 can be made much more complex – however, experimentation has shown that the simple “Do not worry about the spool up and spool down time for the engines, just stop the burn when the calculated burn duration has elapsed” produces the best results.

Performing Burns (3/4)

One method that seems like it would be more accurate is replacing step #7 with this:

• Start a loop that executes as frequently as possible (but not more than once / physics tick):
  1. Calculate the amount of acceleration that has occurred. The following formula will produce the right vector:
     
     
  2. Multiply the acceleration by the length of a “physics tick” to produce a velocity magnitude. This is likely somewhat variable, as it depends on a wide range of factors.
     

     In Juno, will return the correct result – if, and only if, you can complete the execution of all the code in your monitor loop in a single physics tick. It may also be incorrect if there are many active Vizzy threads, which could result in your burn monitoring code not completing in a single physics tick.
     
     My recommendation is to avoid "Frame Delta Time" and simply compare the "Time From Launch" timestamp from the previous iteration with the current one. You must use this method if your platform does not support direct reporting of the length of the physics ticks.

  3. Multiply the velocity magnitude by the vector pointing the same direction as the craft. If necessary, normalize this vector before performing the multiplication.
     

     In Juno, this Vizzy block contains the required vector , and it is already normalized.

     If you can avoid it, you should not use the direction you have requested the craft to point, as the craft may not be pointing in that direction (due to, for example, asymmetrical thrust, insufficient time allocated to maneuver to the burn attitude, or other factors).

  4. Add the velocity increment to the velocity accumulator. This is how much delta-v (and in what direction) you have added to the craft.
     
  5. Subtract the desired delta-v vector from the velocity accumulator. This is how much delta-v you need to add (and in what direction) to complete the burn.
     
  6. Order the autopilot to point in the direction of the vector from the previous step.
     

     As described earlier, this is absolutely required in Juno. It is also correct (doing it this way produces the most accurate burn).

     In Juno, is the right instruction. Be very aware that, despite the implication of “Lock heading,” it does not track the PCI vector specified.

  7. Set the throttle to something like the following expression:
     
     Where “clamp” is a function that returns the value of the first operand if the value lies between the second and third, but otherwise returns the value of the second operand if the first operands value is less than it, or the third operand if the first operands value is greater than it. This formula maintains a throttle setting of 100% (1.0) until the amount of velocity to gain is less than 15 m/s, then drops to zero. Ideally, you want the denominator to be as small as possible without producing over burn. Experimentation is highly recommended.
     
  8. Exit the loop when the magnitude of the velocity remaining to gain (as calculated in step #5 drops below a specified value. I used 0.01 -- too small of a value will result in "spinning", when the craft starts rotating as rapidly as it can trying to chase down a tiny amount of error while setting it too high will reduce the accuracy of the burn.

This method will be referenced in a bit as the “Moderately Complex” burn logic. The downside of this logic is that it always burns longer than expected (due to it throttling down at the end of the burn) and it always underburns by a slight amount. The upside is that the gradual throttle down eliminates overburn caused by the spool down time of the engines.

Performing Burns (4/4)

As another alternative, I came up with what I consider a rather ugly hack. It does work, and you could consider it to be a simple form of “closed loop guidance,” but… I still think of it as a hack. To implement, modify the initial algorithm as follows:

• During the setup for the burn (before you ignite the engines), add the following steps:
  1. Convert the PCI delta-v vector into a “Prograde, Normal, Radial” vector. The required unit vectors to perform this conversion are:
     
     

     Assuming you are following along with this guide, this will be incorrect for Juno (or any other system that uses left-handed cross products). The corrected equations are:
     

  2. Take the absolute value of each component and compare it to the others two components. Once you have identified the largest component, map it to one of the attributes (the “key attribute”) of the orbit as follows:
     • If the prograde component is largest, use “a” (the semi-major axis, which is the average of apoapsis and periapsis).
       
     • If the radial component is largest, use “e” (the eccentricity vector).
       
     • If the normal component is largest, use “h” (the angular momentum vector).
     For the key attribute, calculate a formula that compares the current value (as calculated from the current velocity and position vectors) and the desired value in the targeted orbit (as calculated by adding the desired delta-v to the velocity of the current orbit at the selected true anomaly) and store the result in a variable.
     

     This assumes that your programming environment supports evaluating a string as program code, which Vizzy does. If your programming environment does not, you will need to use another “if” statement in your monitor loop to select the appropriate (hard coded) formula.

     If the key attribute is a this is just a matter of subtracting the absolute value of the current “a” minus desired “a.”
     For “e” and “h,” find the angle between the current and desired vectors. This will always be a positive number, so no absolute value is required.
     

     In Juno, this will return a result in degrees. There is no need to convert it to radians, as we will see later.

• Replace step #7 with the following:
  1. The craft should be pointed in the direction of the original delta-v vector rather than the residual delta-v vector (you should still update the attitude with each iteration – but now you are setting it to the same value each iteration). This has the plus of preventing spinning altogether, but it also means that the you will never actually achieve the targeted delta-v. This is fine, because we no longer care about matching the input delta-v value.
     

     It may still be useful to calculate this value – you can use it to throttle down when the delta-v required gets low, for example, and it is useful for troubleshooting and display purposes. It no longer determines when the burn is completed, however.

  2. Calculate the current orbit (from the current velocity and position vectors) and calculate the difference between the current value of the key attribute and the desired value (same as before, but the current orbit has obviously changed).
     
  3. Subtract the value that you calculated in the previous step with the value calculated during the previous iteration. If:
     • The value is positive, then the error is shrinking and we should continue to burn.
       
     • The value is negative, then the error is increasing and we should stop burning.
       
     • The value is zero (very, very rare) then we have exactly matched this attribute and we should stop burning.
  4. Update the value that we calculated originally with the values from the current iteration.
• The exit condition for the monitor burn loop is the key value remaining the same or increasing between two iterations.

I refer to this as the “Complex” burn logic.

To make a long story short, rather than stopping the burn when we add a pre-determined amount of delta-v, we stop the burn when the key attribute of the current craft’s orbit gets as close as possible to the desired value while burning on a constant PCI vector.

The above method may produce better results if and only if the burn is dominated by a single purpose – that is, it is either aligning planes, rotating apse lines, or changing periapsis or apoapsis. If you are doing two or more of these things in a single burn, the logic above will optimize for one of the three, at the cost of the other two. It is possible to detect when this is happening (as two or more components of the PNR vector will have “significant” magnitudes) but I do not know what the proper standard for “significant” is.

Two examples of multi-attribute burns:

1. Later in this guide I cover a procedure to rotate apse lines and change apoapsis and / or periapsis at the same time. Obviously, this will result in an orbit that differs from the original in terms of both e and AP/PE (and therefore “a”), which may lead to problems. Whether or not it causes problem will depend on the magnitude of the change in “a” versus the amount of apse line rotation requested (as a fraction of a 360°). If a large rotation is requested and only a very small change in “a” then the burn will be almost exclusively radial and all will be well. On the flip side, if the change in “a” is very large, then even a large rotation of the apse lines will be “folded into” the prograde / retrograde maneuver and things are good as well. However, when a relatively small change in “a” is combined with a large rotation of the apse lines, you will end up with a burn that is, say, 20% radial and 80% prograde. In such a scenario, the above logic will optimize for “a” but introduce large errors in “e.”
   
2. It is possible to transfer from one orbit to another orbit that whose plane is not aligned (direction of h is different), apse lines are not aligned (direction of e is different), and whose AP / PE are not the same (“a” is different). All you need to do is use is the procedure described later (in “Transferring to a moon (or ‘More direct rendezvous’)”) but force the rendezvous to occur along the line of intersection between the two orbital planes (departing at the ascending node and arriving at the descending node, for example). If you did, the burn at the end (to change velocities to match the target orbit) could result in a significant values in all three axis of the PNR reference frame.

Which method of performing burns is best? (1/5)

Necessarily, this section will be far more Juno specific than any of the other sections. I decided to include it here, rather than in my other guide, for two reasons:

1. Without the specific details of how the burn termination logic varies between implementations the results are not very useful.
   
2. The advantages and disadvantages of different forms of burn termination logic should apply on any platform -- while the example is Juno specific, the lessons learned should be more generally applicable.

I do try to provide all the necessary background information so that you can follow along even if you do not own Juno.

For testing purposes, I attempted to setup a rendezvous between the active craft whose orbital parameters looked like this:

Be aware that the number in parenthesis (next to AP and PE) are produced by subtracting the radius of the planet then converting into friendly units. So "100 km" = "100 km above the surface." The raw PE value does not subtract the radius of the planet, so these two numbers will disagree with one another.

In addition, the components of e and h are the normalized components of these vectors. This makes it much easier to compare these values between two orbits.

This initial orbit will vary minutely from trial to trial, but it is very close.

The target will be the moon “T.T.” in Juno: New Origins, whose orbital parameters are:

These orbital parameters will be exactly the same in each trial, due to the way in which planets are implemented in Juno.

This is a very challenging target given the initial conditions, because:

1. T.T. is moving in a retrograde direction relative to the active craft. This means that the velocities of the active craft and the moon will be moving towards one another at rendezvous, instead of in the same direction.
   
2. Even setting aside #1, there is a significant difference in the inclination. This means that the rendezvous needs to occur exactly at the line of intersection between the two orbital planes to avoid a large error in the k direction (out of plane).
   
3. The difference between the periapsis of the initial orbit and the periapsis of the target orbit is very, very large – this means that a large change in orbital velocity will be required, and large changes are much less accurate than small changes.
   
4. The sphere of influence for T.T. is very small, so if there is any significant error in the position of the active craft or the moon at rendezvous then a miss is highly likely.

In order to get the required accuracy, I performed phasing over 6 orbits, and after each orbit I recalculated the required phasing and performed the calculated burn. Put another way, I determined how long I needed to wait to rendezvous in 6 orbits, performed the burn to enter that orbit, then calculated how long I need to wait to rendezvous in 5 orbits, and so forth. Thus, each burn had an opportunity to cancel out any error left over from previous burns.

The last burn places the craft on an orbit starts at the line of intersection between the orbital planes and ends on the other point of intersection between the orbital planes. This burn is "blind" -- it simply calculates a burn whose periapsis is a point on the current orbit and whose apoapsis matches the altitude of the target after traveling ½ the period of the transfer orbit. No effort is made to accommodate any remaining phasing error.

A large portion of the reason for the large number of phasing burns is the highly unfavorable position of T.T. and the active craft at the start of the rendezvous process. 4 orbits would be sufficient for accuracy purposes, but trying to achieve the required delay in only 4 orbits resulted in a hyperbolic orbit. Ooops. 😊

The numbers shown in the screenshot are only representative of the expected results. In addition to the errors that I was trying to isolate, there are numerous other sources of variability between separate trials. These include errors introduced by using numeric integration, changes in exactly when physics ticks occur, errors in the launch to orbit code, and likely dozens of other things.

When I repeated the first trial several times, I ended up with results between a collision (several times) and a fly-by at ~200 km. If you try to repeat these tests, do not expect to reproduce my results exactly.

Which method of performing burns is best? (2/5)

First up is the trivial burn termination logic – all this does is set the throttle and wait until the burn duration has finished.

Trivial burn #1

FYI: The implementation used to calculate Prograde / Normal / Radial Delta-V was incorrect when these screenshots were taken. I have corrected the unit vector definitions (above). I decided against re-doing all of these screenshots to fix the problem, as the incorrect definitions used are close to being correct.

All burns are carried out using the PCI Delta-V vector -- the PNR vector is only used in the "Complex" burn termination logic, and even then only to identify the "most important" component, for which the incorrect unit vectors are "good enough."

The key parameter to monitor here is “T” – this is the period of the orbit, and we want the error to be as close to zero as possible. We also need the error in e to be very low, as this indicates how closely we are aligned with the line of intersection between the orbital planes. The results are not horrible, but are not all that good either.

The comparison here is between "The orbit that the burn actually produced" and "The orbit the burn was intended to produce." The orbit of the target (T.T.) is irrelevant to this comparison.

Trivial burn #2



The error in T is very small now (3 seconds) and the error in e has dropped the resolution that I have chosen to display (4 decimal points).

Trivial burn #3




Trivial burn #4




Trivial burn #5




Trivial burn #6




Trivial burn #7




Result of trivial burns
The orbital prediction for the active craft immediately after completing the 7th burn. The brown line indicates that a sphere of influence change is expected.


My calculation of the orbital parameters after the sphere of influence change.


The game's projection of the craft's orbit immediately after the sphere of influence change.


And yes, it does end up colliding with the moon. It is likely that it would have missed had I not dropped out of time warp 60 seconds before the SOI change, as this introduced significant error in the position of the craft at SOI change. Realistically, only two phasing burns were required – the remaining burns were chasing error introduced by the numerical integrator, rather than real errors in the burn.

Which method of performing burns is best? (3/5)

Next up is the “Moderate” burn termination logic. This logic calculates how much delta-v has been added and ends the burn when all but 0.001 m/s of delta-v has been added. In addition, it throttles down once the amount of delta-v remaining drops below 10 and points in the direction of the residual delta-v (rather than the initial delta-v vector, as the trivial logic does).

Moderate Burn #1



This is slightly more accurate than the trivial burn (look at the error in T – 1380 vs 1367).

Moderate Burn #2



This is significantly worse – the error in T for the trivial burn after 2 phasing burns was ~3 seconds, but the error here is ~10 seconds.

Moderate Burn #3




Moderate Burn #4




Moderate Burn #5




Moderate Burn #6




Moderate Burn #7




This is much, much worse than the trivial burn logic. An error of T of 7260 is 121 minutes! The trivial burn logic produced an error of 183 seconds (~3 minutes). This error is catastrophic:

The error in T results in the active craft failing to enter the sphere of influence of the target. Note that the closest approach point occurs well before apoapsis due to the high error in T.
Why are the results so bad? I think the root causes are twofold:

1. The logic used always results in an under-burn of 0.01 m/s. This is necessary to avoid spinning, but it means that a small amount of error is guaranteed with this algorithm.
   
2. Throttling down before the end of the burn eliminates the error caused by spool-down time, but it also means that the calculated start time is less accurate.

Combined, these two sources of error produce a large error at the end in the worst case, which is a burn that occurs at periapsis in a highly elliptical orbit.

Which method of performing burns is best? (4/5)

Next up, we look at the “Complex” burn termination logic – here, we convert the burn into prograde / normal / radial components and continue to burn (along a constant vector) until the attribute associated with the dominate PNR component starts changing in the wrong direction. We also throttle down as the amount of delta-v required drops (using the same formula as the “Moderate” logic), but once the expected delta-v drops below 0.01 m/s we lock the throttle at a low setting (0.00001) for the remainder of the burn.

The goal here is to reduce the error introduced by non-instantaneous burns.

This is a somewhat unfair test, because it is known that this logic will introduce errors when two or more components of the PNR vector are of similar magnitudes. However, I did not try to adjust for this, in large part because it would default back to the “Moderate” burn termination logic for many of the burns, which would defeat the point of making all these screenshots. 😊

Complex Burn #1



The error in T is bit less than either of the two other cases that we looked at: 1362 for this burn, 1380 for the trivial logic, and 1367 for the moderate logic. However, it is not that much better, and this is pretty much the best case for this logic.

Complex Burn #2



This burn is worse than either of the two cases presented – it is close to the accuracy of the “moderate” logic (9.9 seconds vs. 9.76 seconds) but far inferior to the “trivial” logic (9.9 seconds vs. 3.1 seconds)

Complex Burn #3




Complex Burn #4




Complex Burn #5




Complex Burn #6




Complex Burn #7




And… That is catastrophically bad again. The T error is slightly less than with the moderate burn logic, but 7194 vs. 7260 is not anything to write home about. It still results in missing the target:

Which method of performing burns is best? (4/5)

So… Conclusions? Any algorithm that involves throttling down prior to the calculated end of the burn is likely to be less accurate than just burning at full throttle for the entire duration. The largest source of error when transitioning between two orbits appears to be caused by the non-instantaneous nature of the burn, so anything you can do to make the burn shorter is likely to reduce error.

The exception is when the burn length is very short (< 10 seconds is what I settled on) – in this case, extending the burn length reduces error, as it reduces the error introduced by spool up and spool down time.

So…

Yeah. That was easier than I had thought it would be.

Thoughts on further increasing burn accuracy

The major source of error in engine burns is caused by adding delta-v over time rather than in a single impulse. Doing so creates an unavoidable source of error – for example, if you thrust purely in a prograde direction starting 10 seconds before periapsis and ending 10 seconds after, you will adjust your periapsis by a small amount. Not very much, obviously, but there is no way to perform a 20 second burn such that it only increases your apoapsis while holding periapsis constant. And since some of the energy went towards increasing periapsis, your apoapsis value will also be incorrect.

In the “Real World”™ this is handled by using the Kelper’s equations as a starting point, but the actual burn length / duration / consequences is determined by an iterative simulation. You simulate performing the burn, observe the resulting orbit, then try altering the direction of the burn (which could change over the course of the burn) and magnitude of the burn (which could also change over the course of the burn) repeatedly until you find an optimal burn. Such a solution is required if you need to account for the gravitational effects of multiple bodies at the same time (the “N-Body Problem”), which is an important consideration with many real-world burns (the moon gravitation has a significant effect on the orbits of crafts in the orbit of the Earth, for example).

However, implementing such a solver would be very complicated and is far beyond the scope of this guide.

As a result of this, the longer a burn takes the less accurate the result of the burn will be. How long is too long? I honestly have no idea, as I have not performed any experiments. As a starting point, the following things are probably worthy of looking at:

• The lower the orbital altitude at the time of the burn the more error introduced by burning longer. This is especially true if the orbit is highly elliptical (or is parabolic or hyperbolic) and you are near the periapsis point. The faster the orbital velocity vector is changing due to orbital effects, the worse the effect of a longer burn will be.
  
• The lower the mass of the central body the more problematic longer burns will be. This combines with the previous point – a smaller body also likely supports a lower periapsis.
  
• If the actual burn duration is long relative to the period of the initial orbit, you will almost certainly run into problems. If the period of your orbit is 30 minutes, a 5 minute burn is likely to be highly inaccurate, and a 10 minute burn won’t do what you want at all.

Normal (and radial, if eccentricity is low) burns are much more sensitive to the above, because these directions can shift very quickly under the right circumstances and even in the best case will change quicker than the prograde direction. This is why you cannot setup a maneuver node to burn in the “Normal” direction and produce a retrograde orbit – attempting this will produce a hyperbolic orbit or an orbit that collides with the planet.

The math (in the section “Match planes between two orbiting objects”) will produce the optimal burn to shift inclination by 180°. The burn produced will not be in the “Normal” direction – it will be retrograde (relative to the initial orbit), and the magnitude will be twice the current orbital velocity at the point of the orbit. Unless special precautions are taken (performing the burn at a very high apoapsis), such a burn will most certainly result in the craft crashing, most likely after running out of fuel. But it would work if all the delta-v was delivered in a single impulse.

However, if you create a maneuver node and starting adding velocity in the normal direction, you the result will never be a retrograde orbit. The issue is that the normal direction that the maneuver node gizmo is adding delta-v in is fixed based on the position of the craft at the point where the maneuver will take place. However, as the burn proceeds, the normal direction will change, yet the gizmo does not take this into account.

Another way to view this effect in practice is to lock the autopilot to the normal direction, then start thrusting. This sort of maneuver will eventually put you in a retrograde orbit (although it will be less efficient than if you had just started thrusting retrograde) because the autopilot will adjust to the shifting normal direction.

If you do run into problems with burn duration the “best” solution is to increase the TWR of the stage that performs the burn with. Assuming that is not possible, the next best solution is to break the burn into several increments. So, for example, rather than increasing your apoapsis from 10 km to 5000 km in one burn, first increase apoapsis to 2500 km and once that completes further increase it to the desired 5000 km.

While out of scope for this guide, a third option would be to recalculate the burn after each “physics tick.” Let us say that you want to increase you apoapsis from 50 km to 5000 km. You start the same way as normal (calculating the true anomaly, delta-v, and start time as described above), but after the first physics tick completes, you recalculate the burn from scratch – that is, you ask the question “Given my orbit now, after one physics tick of acceleration, when and how much do I need to burn increase my apoapsis to 5000 km.” Then use that true anomaly / delta-v to calculate a new start time and repeat the process until you have completed the maneuver.

It is not quite as bad as it seems – if you have calculated the burn length correctly, the start time for the second “physics tick” should match the current time, so you can indeed burn continuously. However, if the start time ever falls into the past, you cannot complete the desired burn on this orbit and will need to wait another orbit to try again (presumably starting the burn earlier for the next try).

Why not use this as the standard method for performing burns? There are several disadvantages:

1. It would be very, very hard to implement without combining the "Perform Burn" and the "Calculate Maneuver" logic into one mega-instruction. In something like C# you could use delegates to maintain code separation, so it is more reasonable in a modding context.
   
2. There is lots of math going on while the burn is in progress. This will make it very hard (likely impossible) to start the burn on time. If you cannot stop the burn on time, then most, if not all the advantage disappears.
   
3. Floating point errors put an absolute limit on the accuracy of your calculations.
   
4. Finally, the numerical integrator adds errors to the burns, and that cannot be avoided either.

All in all, implementing a numerically integrated (vs analytically integrated) maneuver system is not worth it unless you are stuck with very low TWR – like, 0.05 or so. The Dawn spacecraft (and, in general, any spacecraft that uses ion engines for thrust) uses these sorts of techniques. If you look at the interplanetary trajectory chart in the linked page, you will see that Dawn spends most of its time thrusting, especially for the Vesta / Ceres transition.

If your TWR is, say, 0.25 (or higher), then I would recommend splitting your burn into two (or more) pieces rather than trying to implement numerically integrated maneuvers.

Calculating point of impact

Calculating the true anomaly where the craft impacts the planet

Evaluate this formula:


With r set to the radius of the planet. If this returns an error, the orbit does not intersect with the planet -- otherwise, it returns one of the two true anomaly values where the orbit intersects with the planet.

In order to interpret the results of this formula properly, it is very useful to look at the value of the operand for the arccosine function as a separate value. If

1. If , then there are no solutions because the orbit always remains above the surface of the planet.
   
2. If , then there is one solution (when true anomaly is 0°), and the periapsis of the orbit equals the radius of the planet.
   
3. If , then the orbit crosses the altitude that you are checking for at two locations (the second of which is given by 2π-x).
   
   Geometrically, the two solutions can be produced by mirroring the first across the x axis in the perifocal reference frame. If you are checking for impact with the planet, the true anomaly of the craft will always lie between these two values, so which one you are interested in depends on the direction the craft is traveling. Due to the way in which true anomaly is defined, you will always be interested in the solutions whose true anomaly value is both greater than the current true anomaly and closest to the current true anomaly value. In the vast majority of cases, this will be "2π-x."
   
4. If , then there is only one solution (when true anomaly is 180°), and the entire orbit is below the target altitude except apoapsis. This would be the expected result when the craft is resting on the surface of a featureless planet.
   
5. If , there are no solutions because the entire orbit is contained within the planet. This, obviously should never happen. 😊

Determine the location where the orbit intersects the ground on a non-rotating, featureless planet

As above to get θ, then plug the value of θ into the equation

to get the location in the perifocal reference frame. In almost all contexts, you will want to convert this to a PCI vector.

Determine the time to impact

As above to get θ, then use the process described above under “Converting from true anomaly (θ) to time since periapsis (t).” Then calculate the true anomaly of the craft's current position, convert that to a time since periapsis and subtract the two values to get "time from current position to impact."

This is a quite common operation -- you may want to implement a function that takes two true anomaly values as input and returns the time from the first to the second.

Determine the location where the orbit intersects the ground on a rotating featureless planet

As above, but you need to convert the target location into a latitude / longitude position, then offset the latitude component by an amount based on the rotation of the planet at the target’s longitude component and the time to impact.

Given that we know how long the planet takes to complete a single revolution (T) in seconds, we can calculate the number of degrees of rotation / second with this formula: , where Φ is the latitude of the point of interest. We can then multiply this value with the time of impact (calculated in the previous step) and add it to the longitude component of our calculated impact location.

Determine the location where the orbit intersects the ground on rotating planet with irregular terrain

As above, but once we have our adjusted impact location, determine the height of the terrain at the impact point. If this number, plus the radius of the planet, is not the same as the value we used for r in the first step, we need to repeat the whole process, this time setting r = terrain height + radius of the planet.

To get accurate results at low altitudes, you will want to cap the terrain height so that it is never larger than the altitude above sea level of the active craft. The issue is that the initial impact point (ignoring terrain) may put the impact point inside, say, a mountain. Adding that value to the radius will produce a result that indicates the craft is inside the planet, which will prevent you from calculating a new impact point. With the height capped at the craft’s altitude above sea level, however, it will return a result, and that result will have a true anomaly closer to 180°, and you can continue to refine the impact point.

How many iterations this will take depends on the terrain near the impact zone, but it should always converge to a single latitude / longitude position.

This still is not perfect – a fully robust algorithm would need to, after the initial impact point was selected, reduce θ by a small amount (say, 0.002 radians, or ~ 0.1°), find the perifocal position at the new true anomaly value, convert perifocal to PCI, convert PCI to latitude / longitude, adjust for the rotation of the planet, then get the terrain height at the latitude / longitude location. If the terrain height is greater than the altitude derived from the position vector then there is a mountain in the way and the process should be restarted with an updated r value for the new location. If no collision is detected, then continue to decrement θ until the measured altitude is large enough that you are no longer concerned about terrain collisions (say, 5000 meters). This is, of course, very expensive to implement, especially at low altitudes.

Adjust the Apoapsis / Periapsis to a specified value

If we want to adjust apoapsis, then the burn will occur at periapsis, and therefore at a true anomaly (θ) of 0°. If you are attempting to adjust periapsis with a burn at apoapsis, the example is the same, just set θ = 180°.

The initial perifocal velocity at apoapsis / periapsis is given by this formula:
With θ set to the proper value. Note that all the velocity will be in the q direction, or tangent to the orbit. This is expected at periapsis or apoapsis, as that is the point where the distance from the center of the planet (r) reaches its minimum or maximum value.

The new orbit will include the point where we made the burn, but the velocity at that point will be different and therefore the angular momentum, eccentricity, and semi-major axis will all be different.

• First, calculate the semi-major axis (a) of the new orbit using the formula:
  
• Then we use the formula for eccentricity (e) to calculate that:
  
• Then, we rearrange terms in the first definition of the semi-major axis to define angular momentum h in terms of a and e: and calculate that. While this formula returns two results (differing in sign), we are only interested in the positive result, as h is also defined as the magnitude of a vector (“h”), and the magnitude of a vector is always positive.
  
• Next we need to deal with the direction component of e. If the new apoapsis is greater than the new periapsis, then we can just use the normalized e from the current orbit. If not, we need to invert (multiply times -1) the normalized e vector in addition to flipping the apoapsis and periapsis values.
  
• Finally, we know the direction of h will match the direction of h for the current orbit as we are not attempting to change the inclination of the orbit.

With the above information, we have the new apoapsis, periapsis, magnitude of h, magnitude of e, direction of h (which is just the normalized h of the current orbit), and direction of e (which is the normalized e of the current orbit, possibly inverted). This is sufficient to fully define an orbit, and we should do that.

It is strongly recommended that you create a callable procedure that calculates an orbit in these terms that you can call as required. This will come up very frequently.

Next, using the new (“Transfer”) orbit, we calculate the velocity at the matching point in the new orbit using the same equation as previously. Again, the velocity will all be in the q direction (as defined by the transfer orbit’s perifocal reference frame).

If we flipped the e vector when we defined the "Transfer" orbit, you will need to add π (180°) to the true anomaly value we used earlier to calculate velocity.

Convert both velocity vectors into PCI velocity vectors (using the matching conversion factors) and subtract the velocity of the transfer orbit from the velocity of the initial orbit. The magnitude of the resulting vector will be the amount of delta-v required, and the direction of the vector will be the direction that we need to point the craft to achieve the required delta-v in the most efficient possible way.

In a very real sense, the guide could stop here. All possible changes to an orbit can be expressed as “Calculate the target orbit, find a point at which the current orbit and target orbit intersect, then subtract the PCI velocities at the point of intersection to determine a burn that needs to be performed at the true anomaly value of the intersection.” Sometimes it may not be easy to find the points of intersection between the two orbits, and sometimes it may not be easy to work out what the transfer orbit should be, but philosophically that is all that needs to be done.

Match planes between two orbiting objects

The first task is to find the two points of in which the active crafts plane intersects with the plane of the target vessel. This turns out not to be so bad: This formula will return a unit vector that points in the correct direction (using the PCI reference frame):


From this, we can define the two points on the orbit of our active craft that fall on the line described by L:

1. Convert L and -L into the perifocal reference frame of the active craft.
   
2. Convert the perifocal values from the previous step into true anomaly values for the active craft’s current orbit. This can be most easily done via the formula , although there are other ways to do it as well.
   
3. Select one of the true anomaly values from the previous step to perform the burn at. The best way to do this is to calculate the r value at each true anomaly candidate with this formula:
   
   And select the one with the larger r value. Plane change maneuvers are always cheaper the further away you are from the planet.

There are two ways to proceed from here.

A somewhat limited way (see below) way is to construct a transfer orbit with the following characteristics:

• The same apoapsis and periapsis as the current orbit. This implies that the semi-major axis (a), magnitude of eccentricity (e), and magnitude of angular momentum (h) will also be the same.
  
• The direction of h should be the same as the orbit of our target orbit.
  
• The direction of e should be the same direction as our current orbit projected onto the plane which is normal to h of the target orbit.
  

  For instructions on how to project a vector onto a plane defined by its normal vector, see the appendixes (at the end of the guide).

• The magnitude of h should be the same as the current orbit (and it can be calculated from the apoapsis and periapsis values of the current orbit).
  
• The magnitude of e should be the same as the current orbit (and it can be calculated from the apoapsis and periapsis values of the current orbit).

The calculated orbit is the orbit that the active craft is currently in – but it is now the planes are properly aligned.

Due to the way in which the perifocal reference frame is defined, the active orbit and the transfer orbit will appear identical in their respective perifocal reference frames -- that is, for any given true anomaly value the distance from the center of the planet (r) and magnitude of velocity (v) will be the same. This is expected, because one of the useful qualities of the perifocal reference frame is to eliminate differences that result from inclination.

In this case, we can take advantage of this by noting that the true anomaly values where L intersects with the new orbit will be the exact same values that we calculated previously. Due to the way in which we defined L these are also the only two points where the orbits intersect in 3 dimensional space.

Thus, we can finish up by calculating the PCI velocity for both the inital and transfer orbits at one of the two true anomaly values we calculated, then subtracting inital from transfer. The result will be the amount of delta-v we need to gain and the burn should occur at whichever true anomaly value we selected.

This is very slightly limited: when the relative inclination of the two orbits is 90° and the line of intersection is the same as the line of apse, then the result of projecting e onto the desired orbital plane will be the zero vector, which has no direction. An alternative solution bypasses this problem.

This process is much more similar to the process that you will see described for calculating orbital maneuvers online. The only change is that the online examples will perform all the operations in the perifocal reference frame, while the below procedure uses the PCI reference frame.

This is a case of being homework / exam friendly -- converting perifocal to PCI requires lots of math, none of which is relevant to the question "Do you understand what is going on" so it makes sense to stick to the perifocal reference frame.

On the other hand, if you are actually trying to calculate burns, you need to use the PCI reference frame and (in this case) there is no downside to sticking with the PCI reference frame throughout.

1. Find the PCI velocity vector (v) for the current orbit at the selected true anomaly value θ as described above.
   
2. Calculate the vector that is tangent to the current orbit at the selected θ. This can be done with the cross product:
   
   

   h is the PCI angular momentum vector of the current orbit.

   It does not matter:

   • Whether the cross product is defined using the left-handed rule or right handed rule,
   • Whether you fixed the h vector to point in the correct direction or not or,
   • What order you perform the cross product.

   The only important thing is that all of the above are the same in step #4.

3. Calculate the angle (in radians) between the tangent vector and the velocity vector at θ.
   
4. Apply the equation in step #2 again, this time using the h vector from the passive craft.
   
5. Normalize the tangent vector and rotate it the distance you calculated in step #3, using h of the target orbit to define the plane to perform the rotation in.
   

   For instructions on how to rotate a vector in a plane, see the appendixes.

6. Multiply the tangent vector by the magnitude of the velocity vector calculated in step #1.

The velocity vector from step #1 is the initial velocity, and the velocity vector from step #6 is the final velocity (both in the PCI reference frame). Subtract final velocity from initial velocity to produce delta-v, and the burn true anomaly will be the θ value you used in step #1.

If desired, you can calculate the parameters of the transfer orbit from a state vector composed of the PCI position of the current orbit at θ and the final velocity calculated in step #6. This is not, however, required.

Aligning apse lines between two co-planar orbits (1/2)

The “apse line” is the line that connects the periapsis and apoapsis. The apse lines between two orbits are “aligned” when the angle between them is zero. The e vector always points towards periapsis, so this vector suffices to define the apse lines. Due to the way in which it is defined, the p unit vector points in the same direction.

This section assumes that both the target and current orbit are in the same orbital plane. If they are not, you can still follow this procedure if you project the e vector of one orbit into the orbital plane of the other. This process will fail if the angle between the two planes is 90°, as the projection of the e vector will return the zero vector.

A pure “Align Apse” burn will always be either in the radial up or radial down direction. The delta-v cost of such a burn depends primarily on the eccentricity of the initial orbit – the closer it is to 0 (the smaller the difference between initial apoapsis and initial periapsis) the cheaper it will be. The cost increases very, very rapidly as e increases, and if e is, say, 0.1 the cost can exceed 1 km / second in delta-v.

You can dramatically reduce the cost of aligning apse by changing apoapsis and / or periapsis at the same time. If you have ever done a manual rendezvous, you have probably learned the trick of “Burn prograde at the point on your orbit that aligns most closely with the periapsis of the target,” which is a highly effective way of aligning apse lines. This same trick can be used here, but… The math involved is excessively complex to include in this section. Once you have implemented the code in this section, see “Aligning apse lines and adjusting apoapsis” for a way to align apse lines and adjust apoapsis and periapsis at the same time.

This is amazingly hard, and it is not something that a great deal has been written about on the Internet. The following is based on this[space.stackexchange.com].

It seems that this should be straightforward. Since we are assuming that the two orbits are co-planar, all we must do is find the points of intersection between the two orbits and then subtract the velocity of the initial orbit from the transfer orbit. Easy-peasy.

Not so fast! To find the points of intersection between two orbits, we need to be able to express the orbit equation for both orbits using a shared perifocal reference frame. If we do not, the true anomaly values for the two orbits will not be the same, and we will end up with a single equation with two unknowns, which cannot be solved analytically. So…

The first step is to create a transfer orbit that has the following characteristics:

• The apoapsis and periapsis is the same as the active orbit. By implication, the semi-major axis (a), angular momentum (h), and eccentricity (e) will also be the same.
  

  A later section will relax this constraint.

• The direction of the angular momentum vector (h) will also match the active orbit.
  
• The direction of the eccentricity vector (e) will match the target orbit.

For the rest of this section, the active orbit will be referred to as “initial” (or just “i”) and the transfer orbit will be referred to as “final” (or just “f”). Yes, this is somewhat confusing, but I wanted to stick with the nomenclature used by the example I was following.

We start by creating a formula that relates the true anomaly of the initial orbit with the final orbit. This formula will produce the required angle (∅):


There are several alternate methods to get this angle.

• Assuming your environment supports it, will return the correct value. Do not forget to convert degrees to radians, if required.
  

  Despite what the web page[wnp78.github.io] says (at the time this guide was written) the FUNK “angle” operator returns results in degrees, not radians. If you want to use this operator in FUNK, you will need to convert the result to radians (multiply by pi over 180).

• Calculate the PCI position vector at periapsis for the initial orbit. Then convert that PCI vector into a perifocal vector relative to the final orbit, and calculate a true anomaly from that perifocal position vector. The true anomaly returned will be the angle between the two orbits.
  

  This works because the formula for converting a perifocal position vector into a true anomaly only depends on the direction of the vector. Thus, it does not matter if the perifocal position vector points at a point on the orbit or not, only that it points in the correct direction.

  This is the best solution, because it returns an angle in the range [0°…360°], while the other two options only return angles between [0°…180°]. It can return the full range of values an orbit has a “direction” associated with it (counter-clockwise) and that allows the full range of values to be returned.

• The p unit vector is simply the normalized e vector, which is the PCI position of the craft at periapsis. Since we are only interested in the angle, this works just as well:
  

With this value in hand, we can rewrite the orbit equation into a more general form:

The important change is to add ∅, which allows us to rotate the apse line. Up to this point we could make the simplifying assumption that the apse line aligned with the x-axis, but now this will not be true for one of the two orbits.

The two orbits intersect when the r values are the same, so we have:


Now we use to express both true anomalies in terms of one, then use the angle sum trig identity[en.wikipedia.org] to further simplify:


To make the math easier in the following steps, we define some additional variables as follows:


Substituting back produces a much simpler equation (remember, we have not actually simplified anything. All we have done is mask the complexity of the original equation):


Further define:


Using another trig identity (here[en.wikipedia.org]), we can rewrite this into terms of α and finally solve for θ:

Aligning apse lines between two co-planar orbits (2/2)

This is the easiest form of the equation to implement in code, for two reasons:

1. We eliminate the need to calculate several intermediate values multiple times, with both costs CPU cycles (not a big deal) and provides numerous opportunities for bugs to develop in the code (which is a big deal).
   
2. The intermediate values (a, b, and c) are useful in and of themselves:
   1. If a is zero, the two orbits do not intersect at all (one is entirely contained within the other).
      
   2. If b is zero, the two orbits share the same apse line (∅ must be zero).
      
   3. If b and c are both zero, the two orbits are identical (they intersect at all points).
3. You can easily solve for the two points of intersection by implementing the final equation twice (once with a positive sign and once with a negative sign).

If you prefer, we can substitute a, b, c and α back to produce this equation:


From here, the task is simple – we need to use the process described in the next section to find the “Points of Interception Between Two Co-Planar Orbits,” find the PCI velocity of both the initial and final orbit at matching points, then calculate the difference between the final PCI velocity and initial PCI velocity to generate the required delta-v. The true anomaly the burn should occur at will be the true anomaly of the initial orbit at the intersection point.

Points of Interception Between Two Co-Planar Orbits

Note that this process will only work if the orbits are very close to being co-planar.

At first glance, this seems like a simple task – find the angle between the two e vectors (∅), then find the angle to one of the two true anomaly values. Assuming the angle between the e vectors is not 0° or 180°, then if we substitute these four values into the orbit equation the r values returned match in pairs. So, it seems like the equation is stating “At true anomaly value θ the final orbit passes through the same point in space when the true anomaly of the final orbit is “∅ ± θ.“ In practice, however, it does not seem to be that easy. ☹

All the sources that I could locate on this topic believe that this is true, and act as if it is true, and (I assume) produce correct results assuming that it is true, but… Experimentation reveals that only one of the ± values returns points that actually match up when you convert all the way back to PCI coordinates. At least I could not get it to work consistently via a simple, analytical method.

If you generate all 4 potential values returned by the above equation (2 for the ±, then you can invert the sign of each result to produce an equivalent result), evaluate “∅ + θ” for each of those results you end up with 4 true anomaly values in terms of the initial orbit, and 4 true anomaly values in terms of the final orbit. If you then calculate the PCI position for each of these points then calculate the distance between each of the possible combinations you do find “matched pairs” (values where the two true anomaly values produce the same PCI position).

A couple of examples of this in action:

Note: 1 & 2 are the results of the equation for both “active” and “transfer.” 3 & 4 are generated by negated 1 & 2 (respectively). If the process worked “as expected,” then there would always be a match between 1 & 1 and 2 & 2.


In the first screenshot, note that 2 & 4 and 1 & 3 match up.


In the second screenshot 2 & 4 and 1 & 3 match up.

I suspect that there is a clever way to get directly match up values – but I cannot figure it out if there is. Thus, I must recommend simply calculating all 4 values and comparing the distance between all possible pairs.

The above screenshots were taken before I switched to using rotations for perifocal -> PCI conversion, and the relatively large errors is caused by errors associated with using unit vectors for conversion. Once I switched the error at intersection dropped to < 0.0001 (but still not 0).

However, much later I discovered that, under certain circumstances, even using the rotation method for vector conversion could produce errors > 0.0001. I eventually resorted to storing all 8 distances values in a list, sorting the list, finding the "2nd smallest distance," then making another pass through the all the pairs to extract the two pairs with the lowest PCI distance at intersection.

Ugh, very inefficient, but it does work.

Adjust the phasing between two orbits

To produce useful results the two orbits must have aligned apse lines. If they do not, both craft will indeed arrive at their respective periapsis points at the same time – but those two points might be on opposite sides of the planet.

It is possible to align apse planes and adjust phasing in a single burn. However, “Aligning apse lines and adjusting apoapsis” is an advanced topic and will be covered later.

This section assumes both orbits are elliptical. In some cases it may be possible to adjust phasing between a hyperbolic orbit and an elliptical orbit to result in both craft arriving at periapsis at the same time. However, there is no real practical reason to do so and it would add unnecessary complexity to this section.

The “phase” between two orbiting objects is the difference in true anomaly of the two objects when they are both at periapsis. Note that each object has a unique definition of true anomaly (as this is measured relative to periapsis of that object).

1. We decide to perform the burn at periapsis (apoapsis will also work). We need to figure out how much time will pass until the active craft reaches periapsis, which is easy enough:
   
   
   Where T is the period of the current orbit, given by this equation:
   
   
   And the time from periapsis is given by the process described above under “Converting from true anomaly (θ) to time (t)”, with the true anomaly value calculated using the current position of the craft.
   
2. We need to increment this value by one full period – this gives us the time to the periapsis after next, which is when we want the periapsis times to align. This will be used in the next step to determine the earliest possible rendezvous time.
   
3. We need to determine where the target vehicle will be in terms of its true anomaly, so we follow the procedure in step #1 for the target object. If this value is less than the value calculated in the previous step, increment by the period of the target object as many times as required to place it in the future. This is the actual time that the “rendezvous” will occur.
   
4. Next, subtract the time value from step 2 from the time value from step 3 – this is how much longer we need to make the period of our next orbit.
   
5. Next, we need to solve the equation given in step #1 to determine the new semi-major axis (a) that will increase the period by the required amount. This formula will do the trick:
   
   

   To raise something to the “2/3” power, square it, then take the cube root.

6. Finally, we need to calculate what we need to set our new apoapsis to achieve the required value of a. This is just a matter of solving the alternate form of the equation for the semi-major axis for the apoapsis value. This formula will do the job:
   

Now, we just use the process described in the previous section “Adjust the Apoapsis / Periapsis to a specified value” to change apoapsis.

The two craft will reach their next periapsis after the burn at the same time. Unless you take further action (namely, performing another burn to adjust the apoapsis value of the active craft to match the target craft) at that point, however, they will not reach periapsis at the same time for any future orbits.

Orbital rendezvous (slow)

At this point you have all the tools to perform an orbital rendezvous:

1. Use the procedure outlined in “Adjust the Apoapsis / Periapsis to a specified value” to adjust the periapsis to match between the two orbits.
   
2. Use the procedure described in “Match planes between two orbiting objects” to align the planes.
   
3. Use the procedure described in “Aligning the apse lines between two co-planar orbits” to align the apse lines.
   
4. Use the procedure described in “Adjust phasing between two orbiting objects” to ensure both the target and the active craft arrive at periapsis at the same time.
   
5. Use the procedure outlined in “Adjust the Apoapsis / Periapsis to a specified value” to adjust the apoapsis values to be the same.

In theory, you should adjust the periapsis value you set in step #1 to something slightly different from the target object's periapsis -- otherwise, a collision is possible. In practice, I cannot make the burns accurate enough to make this an issue. Generally, separation at rendezvous falls in the range 50 km - 100 km. This can be dramatically improved by increasing the number of phasing orbits and performing burns to refine the period of the orbit with each iteration. With 2 phasing orbits I can reliably produce closest approaches < 5 km, with almost all the error resulting from an incorrect periapsis altitude.

This particular sort of rendezvous works poorly for getting an encounter with a moon -- it will work, but the requirement to increase periapsis in step #1 will make it a very, very slow process.

Adjusting apoapsis while rotating apse lines

This section includes far more math than the previous sections. That is because, as far as I can tell, this is original to me. While I am sure that someone has done this math before (it is not that obscure), as far as I know nobody has published it to the Internet, so if someone is looking for information on how to do this, this is the guide that they will find.
As such, it is more important to show how to reproduce my math than it is in other sections.

As mentioned in the “Aligning the apse lines between two co-planar orbits” section we could specify a new apoapsis and / or periapsis values when we construct a new orbit. Trivial, actually. Well, sometimes you can do that and everything still works properly – and sometimes you cannot.

A diagram will (hopefully!) help:


In this diagram, both the blue orbit and green orbit have the same apoapsis and periapsis values (2000 and 1000), and the green orbit is rotated 45° relative to the blue orbit. The origin of the plot is one of the two foci of the ellipse (and is also the location where the center of the planet would lie), and that is where the “45° rotation” is measured from.

For clarity, the rotation is counter-clockwise in these plots, so the green orbits periapsis is in the upper left-hand quadrant and the apoapsis is in the lower right quadrant.

It should be clear that, as long as we hold apoapsis and periapsis constant, we will always find at least one point of intersection no matter how we rotate the second (green) orbit, and as long as the orbits intersect we can perform a burn to move from one to the other.

So, great, we have proven that if you hold both apoapsis and periapsis constant there will always be two points of intersection.

What if we want to change the apoapsis to 2500 and the periapsis to 1500?

0° of rotation -- this will not work.

45° of rotation, will not work.

90° of rotation, will work.

135° of rotation, will work.

180° of rotation, will work.

We can clearly see that sometimes the two orbits intersect (and thus, you can move from one to the other with a single impulse) and sometimes they do not. It would be helpful if we could determine (analytically) when the two orbits will intersect in at least one place, given a certain set of starting conditions and partially defined ending conditions.

We can, and the following sections will document the math to do so, but a warning upfront:

If you only change apoapsis (holding periapsis constant), you are always guaranteed to have at least one point of intersection, no matter what rotation you perform, as long as the new apoapsis value is not less than the original periapsis value. Similarly, you can change periapsis while holding apoapsis and still guarantee at least one point of intersection as long as the new periapsis value is not larger than the original apoapsis value. This makes all the very fancy math in the following sections pretty much irrelevant, since none of the processes described in this guide require varying both apoapsis and periapsis and rotation at the same time.

Since the next few sections are rather optional, here is the algorithm to adjust apoapsis / periapsis / rotation at the same time:

1. Calculate the desired direction of e using the procedure in “Aligning the apse lines between two co-planar orbits.”
   
2. Calculate the transfer orbit with these orbital parameters:
   • The direction of h will be the same as the same as the current orbit.
     
   • The direction of e should be the value you calculated in step #1.
     
   • The new apoapsis must be larger than the current apoapsis and the new periapsis should be the same as the current periapsis.
     
   • From the new apoapsis and periapsis, calculate e
     
   • From the new apoapsis and periapsis and a
     
   • Calculate h from e and a

With the above information in hand, you should be able to calculate the new orbit, and from here, you can follow the procedure outlined in “Aligning the apse lines between two co-planar orbits,” confident that there will always be at least one point of intersection between the transfer orbit and the current orbit.

Calculating the limit values associated with apse rotation (1/3)

We start with the equation that we generated in the earlier section (“Aligning the apse lines between two co-planar orbits”):

Where:


The limit values (the values beyond which no real solution exists) of the above equation will occur when the operand of the arcosine function is either zero or π, so we can start with this:

And

We want to define the final orbit in terms of its apoapsis and periapsis (rather than e and p), so we start with the four equations that we defined earlier:

a refers to the semi-major axis in the above equations as well in the equations that follow!

a and e are already in the desired form, we start with h:

Then rewrite p into terms of apoapsis and periapsis:

Now we have the tools to rewrite “a”, “b”, “c”, and “α” that we defined earlier in this section into apoapsis and periapsis. Starting with “a”:

"b":

"c":

And finally, "α":

The minimum periapsis value (for a given combination of initial apoapsis / periapsis and desired apoapsis) occurs when one the following two equations is true:

So, we start by converting "c/a" into the desired from:

Then we substitute "c/a" and "α" back into the equation we are interested in, producing this:

I tried feeding this into Wolfram Alpha, and got back "Standard computation time exceeded," so...

But wait, Wolfram Alpha can simplify this equation:

And, shockingly, this simplifies into an expression with no trigonometric functions at all:

We can then define x as:

And substitute for x to produce:

This allows us to “simplify” the original equation to:

With a goal of converting the equation into something that can be easily fed into Wolfram Alpha, we define many substitutions:

Producing:

And Wolfram Alpha still times out. ☹

So, we proceed by hand:

Note that in order to get to this point, we squared both sides (in the third step). Due to this, the difference between the +1 equation that we are solving and the -1 equation that we were saving for later has disappeared. This will cause problems in the future, but for now it is a good thing.

Equation 1
Wolfram Alpha will solve this equation for whatever we want, so we will be coming back to it in a bit. For now, we will continue with manually solving:

The equation above is, believe it or not, a quadratic equation in terms of x (that is why I grouped the terms the way that I did). Therefore, we can apply the quadratic formula to solve for x:

Where:

Substitute back y and z (this both eliminates a variable and allows us to use trigonometric identities to simplify, which is important):
Note:

Group and simplify:





However, by looking at the final answer from Wolfram Alpha (as produced from equation 1) the definition of "b" above is incorrect. b must be:

I have no idea where the error is. For the sake of completing the work, I will proceed with the correct definition of "b."
The quadratic formula is much simpler when c=0:

We split the ± into two equations, one with a + and one with a -.

Or

Substitute:

Substitute back:

Q.E.D.

Going back to equation 1 which Wolfram Alpha will accept, first the "Math Input" that Wolfram Alpha requires for this equation:

solve Square[\(40)\(40)w-x\(41)*v*Sin[θ]\(41)]+Square[\(40)2*w*x*u-\(40)w-x\(41)*v*Cos[θ]\(41)]=Square[\(40)\(40)w+x\(41)*v-2*w*x\(41)] for x”

Using Wolfram Alpha we can solve for a number of potentially useful values starting with "If we vary the final periapsis, what values for final apoapsis produce a single point of intersection?"


And "If we know both the initial and final apoapsis and periapsis values, what apse rotation produces a single point of intersection between the two orbits?”


Now we move to figuring out how to interpret the results of the formula.

Calculating the limit values associated with apse rotation (2/3)

To validate the formula for periapsis, I created an Excel spreadsheet that produced minimum periapsis values for a given initial orbit and final apoapsis. The periapsis values returned are shown in this graph, while the constants are displayed in the title. For the sake of simplicity, all the values are much smaller than they would be the game, but that doesn't change the math.

We start with the most common case -- where the new apoapsis value is greater than the old apoapsis value.






This is straightforward to interpret -- as the rotation increases (up to 180°) the maximum periapsis value increases, reaching the original apoapsis value at 180°. Nice and straightforward. When the increase in apoapsis is small the maximum value of new periapsis increases to old apoapsis very rapidly and becomes much less sensitive to the apse line rotation.

What happens if we decrease apoapsis (but leave it above the initial periapsis)?






This is much more difficult to interpret. Comments:

• The Excel spreadsheet (in addition to calculating the periapsis values) also tries substituting the initial / final apoapsis / periapsis values in the original () equation. These graphs are valid.
  
• The distance between the two “spikes” decreases as the value of final apoapsis approaches the initial periapsis and vice versa.
  To be more precise, the spikes occur due to division by zero, and so occur at the following apse rotation:
  
  
• In the “saddle” region (between the two peaks) the minimum value is the initial apoapsis. The maximum value approaches infinity, which occurs at the spikes.
  
• Outside the saddle region, the maximum value is the initial periapsis and the minimum value is negative infinity.

Some plots of the orbits are required to decode what is happening when apoapsis is decreasing.

In all the following plots, the parameters are as follows:
• The blue orbit has an apoapsis of 2000, a periapsis of 1000 and no rotation.
• The green orbit has an apoapsis of 1500, the periapsis is the result of the equation that we came up with, and the rotation angle is vary from example to example.

A rotation of 70° returns a minimum periapsis of 38.0978, which looks like this:

Note that both orbits share a common focus, and that is what apoapsis and periapsis is measured from for both orbits. Clearly, there is a single point of intersection, which is the expected result.

If we change the green orbits rotation to 71°, we get a minimum periapsis of -35.7767, which is meaningless. If we force the negative value to a very small positive value (0.00000000001) then the plot looks like this:

If we continue to increase the rotation angle, at 89° we get a minimum periapsis of -27149.34425. Changing the negative number to 0.00000000001 produces this:

We are clearly overshooting a “minimum orbit” by a noticeable amount now (which is why the formula returns a negative periapsis), but all is still well – the “orbit” (really, a line) still intersects with the original orbit in at least one place.

At a rotation angle of 90°, the formula is undefined (due to division by zero). If we go ahead and assume the smallest possible periapsis and plot that, we get this:

Ehhh… Special case code could handle this, if need be.

The plot gets… Odd when the rotation angle reaches 91° -- the periapsis is now 30149.34!

If we move along to a rotation angle of 120°, we get a minimum periapsis of 2500, which is also unexpected, but easier to plot:

Ahaha – there is exactly one point of intersection, so the formula is working.

We note, however, that the apoapsis of the green orbit is 1500 and the periapsis is 2500!? That is not right. Let's try flipping apoapsis and periapsis values and plotting that:

Nope, that's clearly wrong. Since we flipped apoapsis and periapsis values, we also need to flip the e vector. We can do this in this context by adding 180° to the rotation. So, we set the apse rotation to 120° + 180° = 300° we produce this plot:

Which is what we generated in the first place, but now the apoapsis is greater than periapsis values are valid.

So whats going on? It all comes back to the original pair of equations:


And the fact that we ended up squaring both sides of the equation while we were solving it. This caused a single equation to cover both equations, and thus it returns to very different values. Sometimes it is returning the maximum periapsis value and sometimes it is returning a minimum periapsis value. And, no matter what, periapsis of 0 is a valid solution.

Geometrically, the root cause of the issue is, for certain rotations, the periapsis and apoapsis values flip (that is, the final periapsis value is greater than the final apoapsis value). When this occurs, the question changes from “If we hold apoapsis constant, what range of periapsis values produce a single point of intersection” to “If we hold periapsis constant, what apoapsis value produces a single point of intersection.”

If we always want to get back a minimum periapsis value (lower bound) we can do this:

1. If final apoapsis is greater than initial apoapsis, return 0.
   
2. Otherwise
   • if , return 0.
     
   • Otherwise, if the result of the formula is negative, return 0.
     
   • Otherwise, if the result of the formula is greater than the initial apoapsis, return 0.
     
   • Otherwise, return the result of the formula.

To determine the maximum periapsis value (upper bound), the following works:

• If final apoapsis is greater than initial apoapsis, return the result of the formula.
  
• Otherwise, return final apoapsis

The above procedure will always produce a range of periapsis values that, when combined with the specified initial orbit, apse rotation angle, and final apoapsis will produce at least one point of intersection.

Calculating the limit values associated with apse rotation (3/3)

While the above formula does produce the correct answer, the limits of floating-point math make it all but certain that if you plug the result back into the original equation you will get back a value minutely less than -1 or greater than +1, and thus get back a result of “No points of intersection found.”

The easiest solution is to “bump” the value by a small amount to the value returned -- “0.1” works. The direction of the “bump” will vary depending on the meaning of the result, so if you are looking for a lower bound you need to increase the value returned by the formula by a small amount and if you are looking for an upper bound you need to decrease the returned value by a small amount.

To sum up, this formula says:

• If you are increasing apoapsis you can perform any apse rotation that you wish as long as you hold periapsis constant. You may be able to increase periapsis as well if the apse rotation is sufficiently large or if the increase in apoapsis is relatively small. In this scenario you never be able to decrease periapsis.
  
• If you are decreasing apoapsis (but not below the initial periapsis) then you can increase periapsis (up to the original apoapsis or the new apoapsis, which ever is less) and can also decrease periapsis to some degree (a large decrease is acceptable if the apse rotation is large).

Orbital rendezvous (fast)

There is another means to achieve rendezvous between two objects:

1. Align planes between the target and the active craft.
   
2. Iteratively:
   1. Find the true anomaly and position of the passive craft at the current time plus a delta (which starts at zero). This will be slow, because you will have to use the procedure described in "To convert time since periapsis (t) to true anomaly (θ)" to advance the true anomaly by a time increment.
      
   2. Using the position of the passive craft from the previous step, find the true anomaly of the active craft when it is exactly 180° away from the target’s crafts position.
      
      Put another way, if you took the PCI position vector for the target craft, multiplied it by -1, what point on the active craft’s orbit would fall in the same direction.
      
      Put yet another way, this is the point on the active craft orbit where you can draw a line that passes through the passive craft’s current position, the center of the planet, and the desired active craft’s position.
      
   3. Find out how much time it will take the active craft to reach true anomaly value from the previous step. This value should always be positive -- if you get back a negative value, add the period of the active craft's orbit to make it positive.
      
   4. Calculate an orbit (this is an example of "Adjusting apoapsis while rotating apse lines."):
      • Whose periapsis value equals the r value of the active craft at the active craft’s true anomaly value from step #2,
        
      • Whose apoapsis value equals the target craft’s r value at its true anomaly value from step #1.
        
      • e and h should be calculated from the above apoapsis and periapsis.
        
      • The direction of h should be the same as the current orbit.
        
      • The direction of e should be normalized PCI vector you calculated in step 2.
      This is the transfer orbit, which you will enter at the true anomaly found in step 2.
      

      If the calculated apoapsis and periapsis values in this step are wrong (periapsis is greater than apoapsis), reverse the new apoapsis and periapsis values and do not multiply the vector from step 2 by -1 when calculating the direction of e.

   5. Calculate ½ the period of the orbit calculated in the previous step (this is how long it will take the active craft to move from periapsis to apoapsis along the transfer orbit).
      
   6. Add the time calculated in the previous step to the time calculated in step 3. This is how long it will take the active craft to reach the rendezvous point.
      
   7. Check to see if the difference between the time calculated in step 6 and the “delta time” value used in step 1 is close enough to zero (I use 0.1). If it is not, use the value from step 6 as the new time delta and try again.
3. Perform the burn to place the active craft on the transfer orbit you found.
   
4. When you reach apoapsis on the transfer orbit, the active craft will be there. At this point you can zero out relative velocity.

Be warned: This method will generally cost significantly more delta-v than the method presented earlier. In my testing, rendezvous with the same craft after highly similar launches that occurred at the same time, the slow method consumed ~600 m/s less delta-v than the fast method. The penalty can be eliminated by ensuring that the apse lines are aligned before the burn and the transfer orbit joins the initial periapsis and ends at the target’s apoapsis (when the period of the target’s orbit is larger than the active craft’s orbit – otherwise, you want to join initial apoapsis to the target’s periapsis). Of course, setting up such an optimal situation defeats the whole purpose of a fast rendezvous method… ☹

The big advantage of this method is that it only requires three burns (one to align planes, one to inject into the transfer orbit, and the last to zero out the relative velocity between the two craft at closest approach).

Additionally, this method is much more suitable to transferring to a moon. As written above (assuming your burns are sufficiently accurate) will result in a trajectory that collides with the moon. If this is not desired, you will want to adjust the apoapsis (or, rarely, periapsis) value in step d to be something less than (or greater than) the actual radial distance to the target’s orbit at the selected true anomaly. You should decrease the apoapsis by some fraction of the size of the “sphere of influence” of the target body – if you subtract the full value, the resulting transfer orbit will just barely touch the SOI of the target, which probably is not very useful. So, say, ½ the SOI radius works well.

The biggest downside of this method is that it is not guaranteed to converge if the initial orbits of the active and passive intersect with one another (prior to the transfer burn). In this case, the apoapsis / periapsis values calculated in step d will sometimes be “correct” (apoapsis greater than periapsis) and sometimes “incorrect.” When you flip the apoapsis / periapsis values you also must flip the apse line which changes where the burn will occur (by 180°) which has a dramatic impact on the total time to reach the rendezvous point (to be exact, it increases or decreases it by ½ the period of the active craft’s initial orbit). This prevents the delta time values from converging, which results in an infinite loop.

There are a couple of ways to address this problem:

• First, you could just advance the initial position of the active / passive craft by a certain amount of time, then retry. Assuming the area of overlap is small, delaying like this will eventually ensure that the apse line never intersects the portion of the orbit where the periapsis / apoapsis need to be flipped / do not need to be flipped. This is what I did.
  
• Second, you can restrict using this method to only orbits that do not overlap. The easiest way to do this is to check to see if there are any points of intersection between the two orbits before any burn is attempted at all. If there are, then you can perform a burn to lower / raise periapsis to eliminate the points of intersection before applying this method.

Finally, this method can be adapted for interplanetary transfers as well, which I am not going to cover. Sorry. ☹

Patched conics (1/2)

It is helpful to be able to calculate patched conics so that you can predict what your orbit will be after you pass into a new sphere of influence. There are two cases to consider:

• Entering in the sphere of influence of a child of the current central body.
  
• Exiting the sphere of influence of the current central body into the parent body's sphere of influence.

The basic idea is the same for both, however:

1. Determine when the true anomaly, time since periapsis, and PCI position and velocity, all in terms of the current orbit (and current parent body)
   
2. Using the information from #1, translate the PCI position and PCI velocity to the new central body.
   
3. Construct a new orbit using the translated PCI position and PCI velocity with the new parent body.

Before you proceed, you should check to see how your chosen platform handles the PCI reference frame for different central bodies. There are two options:

• The only thing that changes when moving from one sphere of influence to another is the origin of the PCI reference frame. The directions of the unit vectors remain constant.
  
• Both the origin of the PCI reference frame and the direction of the PCI unit vectors changes from one central body to another.

In the vast majority of cases, this will come down to one question: Does your system support planets with an axial tilt? If it does, you most likely fall into the #2 camp. Otherwise, you are almost certainly in the #1 camp. Since Juno does not support axial tilt I do not have a frame of reference to advise you on how to handle the changing PCI unit vector directions. Hopefully, your platform provides a mechanism to make these conversion easy.

Exiting the current sphere of influence to the parent’s sphere of influence

This is fairly easy -- if the apoapsis of the current orbit exceeds the radius of the current central bodies sphere of influence or the eccentricity of the current orbit is greater than 1, then the active craft will be exiting the current sphere of influence. The true anomaly at the sphere of influence change will be equal to the true anomaly at the sphere of influence radius, and the PCI velocity and PCI position can be easily calculated from this true anomaly.

With that information in hand, you then need to:

1. Calculate the orbit of the current parent body as seen from its parent (for example, the orbit of the Earth relative to the Sun).
   
2. Calculate how long it will take the active craft to reach the true anomaly at exit.
   
3. Calculate the position of the parent body after the amount of time calculated in the previous step has passed.
   
4. Calculate the PCI position and velocity of the parent body at the true anomaly calculated in the previous step.
   
5. Add the PCI position of the active craft at SOI change to the PCI position of the parent body (as measured relative to its parent) at SOI change to generate the new PCI position of the active craft at the SOI change.
   
6. Add the PCI velocity of the active craft at SOI change to the PCI velocity of the parent body (as measured relative to its parent) at SOI change to generated the new PCI velocity of the active craft at the SOI change.
   
7. Finally, use the vectors from step #5 and #6 as a state vector to calculate a new orbit around the parent’s parent body.

Patched Conics (2/2)

Entering the sphere of influence of a moon (or other child with a sphere of influence

This is much more difficult, and the slowest calculation described in this guide. The problem, as always, is that there is no way to relate the position of a moon in its orbit to elapsed time with algebra alone – it is necessary to find a root (or evaluate an infinite series) to make this relationship. This means that this is no simple formula that can be used to determine if the craft passes through the sphere of influence of a moon.

This is the process that I came up with, which works (but is slow):

1. For each moon (or other object who is a child of the current central body and has a sphere of influence):
   1. Calculate the orbit of the moon.
      
   2. Check to see if the active craft’s apoapsis is greater than the periapsis of the moon minus the sphere of influence radius of that moon or the active craft’s periapsis is less than the apoapsis of the moon plus the sphere of influence radius of that moon. If neither condition is met, the active craft’s orbit is either contained completely within the moon’s orbit or the moon’s orbit completely contains the active craft’s orbit and no further processing is required for this moon.
      
   3. Calculate the true anomaly of the active craft at a radial distance of the moon’s periapsis minus the moon’s sphere of influence. If this is undefined, then return 0.
      
   4. Calculate the true anomaly of the active craft at a radial distance of the moon’s apoapsis plus moon’s sphere of influence. If this is undefined, then return π.
      
   5. For each value “x” between the value in step 3 and 4, incrementing by some fraction (I used 1/10th) of the difference between 3 and 4.
      1. Calculate the time from periapsis for the active craft to reach the true anomaly value of “x.”
         
      2. Calculate the true anomaly of the moon after the amount of time elapsed in the previous step has passed from its initial position.
         
      3. Calculate the distance between the calculated position of the moon and the calculated position of the active craft. If this number is less than the sphere of influence of the moon, then skip immediately to step 9 and the active craft may pass into the sphere of influence of this moon (see below).
   6. Subtract the true anomaly value calculated in step 3 from 2π (for example, “2π – θ”).
      
   7. Subtract the true anomaly value calculated in step 4 from 2π.
      
   8. Repeat step 5 with the new lower and upper bound values.
      
   9. If an intersection with a sphere of influence was detected, see if it is closer to the initial true anomaly position (of the active craft) than any other sphere of influence change that has been detected so far. If it is, then save the calculated true anomaly at SoI change for later processing – otherwise, discard it since the theoretical SoI change occurs after an earlier SoI change, which makes it irrelevant.
2. Assuming that an SoI change was detected (the “back-up” step – the goal here is to find a value that just before the active craft enters the SoI of the moon):
   1. Recalculate (or restore) the increment used in step 1.5.
      
   2. Reduce the increment value by 50%.
      
   3. While the distance between the calculated active craft position and the position of the moon when the craft reaches that position is less than the SoI radius of the moon:
      1. Subtract the increment from the true anomaly of the active craft at SoI change.
         
      2. Recalculate the position of the active craft at the new true anomaly value.
         
      3. Recalculate how long it will take the active craft to reach the new true anomaly value.
         
      4. Recalculate the true anomaly of the moon from its initial position after the amount of time calculated in the previous step.
         
      5. Recalculate the distance between the active craft at the new true anomaly value and the moon at the new true anomaly value.
3. Refine the true anomaly value at SoI change (the “binary search” step – the goal here is to find an exact value, within a specified tolerance, of the true anomaly where the SoI change occurs). The exit condition for this loop is when the error value (calculated in step “b”, below) is below some specified value – say, 1 meter.
   1. Calculate the distance between the moon and the active craft at their respective true anomaly values.
      
   2. Subtract the value in the previous step from the SoI radius of the moon. This will produce a value that is positive if the active craft is outside of the SoI, negative if the craft is within the SoI, and zero if the craft is positioned exactly on the sphere of influence.
      
   3. If the sign of the value calculated in step 2 has changed from the value calculated on the previous iteration, then reduce the increment by 50%.
      
   4. If the value calculated in step 2 is positive then add the increment to the true anomaly of the active craft, otherwise subtract it.
4. Find the position and velocity of both the active craft and the moon at their respective true anomaly values as determined by step #3.
   
5. Subtract the active craft’s position from the moon’s position to find the active craft’s position in PCI coordinates centered on the moon.
   
6. Subtract the active craft’s velocity from the moon’s velocity to find the active craft’s velocity in PCI coordinates centered on the moon.
   
7. Using the vectors calculated in steps #5 and 6, find the orbit of the active craft in the new sphere of influence.

Whew, that is a mess and a half. Some comments:

1. The initial increment (as specified in step 1.5) determines how reliable this process will be. If this increment is very small, then you are almost certain to detect all SoI intersections – but the calculation will take longer. If this value is large, then you will miss some SoI changes, but the process will be a little faster.
   
2. The backup step is required – there is no guarantee that the value returned from step #1 will be closer to the “entrance” (which we want) value than the exit (which is invalid, so we certainly do not want that) value, so if you immediately start the binary search you might find the wrong zero.
   
3. By far the hardest part both in implementation and computationally is calculating the new true anomaly values of the moon. Each and every time you recalculate the true anomaly of the active craft you must determine how long it will take the active craft to reach that true anomaly from its current position (which is reasonably fast) and then offset the position of the moon by that amount of time.

Rendezvous without first aligning planes

We know that two non-coplanar orbits intersect at two points along a line given by the following formula:

If we could ensure that the rendezvous occurred along this line then we could rendezvous two craft without first aligning planes. This is indeed doable, although it is likely to take a large amount of game time to complete. The process looks like this:

1. Calculate the line of intersection, as described above.
   
2. Calculate the true anomaly of the passive craft at L and -L.
   
3. Determine how long it will take the passive craft to reach the two true anomaly values you calculated in the previous step. Find the value of L associated with the true anomaly whose "time to reach" value is the smaller positive value. Leave L unchanged or set it to -L to match the true anomaly value you selected. The true anomaly value you selected is "Passive craft's true anomaly at rendezvous."
   
4. Calculate the magnitude of r for the passive craft at the true anomaly value you found in the previous step.
   
5. Calculate the true anomaly of the active craft at -L.
   
6. Find the r value for the active craft at the true anomaly value from the previous step.
   
7. Create a transfer orbit whose apse line is rotated to match -L – that is, it has the following characteristics:
   • The direction of e is the same as the direction of -L (or L, see below)
     
   • The direction of h is the same as the active craft
     
   • Set apoapsis to the value found in step #4 (if it is larger than the current periapsis) or periapsis (otherwise) while leaving the other value alone. If this change results in a AP / PE swap, define e with L (instead of -L).
     
   • Calculate e and a and h from the apoapsis and periapsis values.
   This is the "transfer" orbit.
   
8. Find the points of intersection between the active craft's current orbit and the transfer orbit.
   
9. Determine how long it will take the active craft will take from its current position to reach one of the points of intersection of the transfer orbit. It does not, generally, matter which one, so you should select the point of intersection which will be reached first.
   
10. Calculate ½ the period of the transfer orbit.
    
11. Add the values in step #9 and #10 together – this is how long it will take the active craft, from its current position, to reach apoapsis (or periapsis) of the transfer orbit from its current position.
    
12. Find the difference between the value in step #3 and step #11. This is the delta time at rendezvous. If this value is zero (very, very unlikely), then calculate and perform the burn to change from the active craft’s current orbit and the calculated transfer orbit.
    
13. Otherwise, if the delta time is less than the period of the current orbit, restart the whole process but this time force the selection of the opposite value in step #3 -- that is, if step #3 selected L, retry and this time select -L. If even this does not suffice to produce a positive value, add the period of the passive orbit and toggle the sign on L again. Continue this process until you get a delta time that is greater than the period of the active craft's orbit.
    

    You will need to persist both the number of full passive orbits to stall and the selected value of L between iterations if you plan on calculating multiple phasing burns, as suggested below.
    
    The issue with recalculating this value is that you may (due to burn inaccuracy) increase the period of the active craft to be slightly larger than you intended. When this occurs, the delta time value will be slightly negative. You do not want to restart the phasing process from scratch when this occurs, so you need to "lock in" the L value.

14. Once you have a positive delta time, then create a phasing orbit with the following characteristics:
    • The direction of e is L (or -L, as you determined above, in step #7)
      
    • The direction of h is the same as the active craft
      
    • The length of the semi-major axis (“a”) is such that the period of the orbit (“T”) is exactly the delta time from step #13. However, there are a few things to consider here:
      • First, you strongly prefer to adjust apoapsis (if the initial transfer orbit’s apoapsis is increasing) or periapsis (if the initial transfer orbit’s periapsis is decreasing) to achieve the required change in “a.” Obviously, if you are increasing apoapsis, you ensure that you do not exit the current sphere of influence and if you are reducing periapsis you need to ensure you do not crash.
        
      • Secondly, you should prefer not to increase apoapsis beyond the passive craft’s apoapsis or decrease periapsis beyond the passive craft’s periapsis. If you do, you will be wasting delta-v.
        
      • In order to achieve both goals, you may need to set the desired period to a fraction (say, ½, or 1/3rd) of the desired period, then wait several number of orbits (2 if ½, 3 if 1/3rd, and so forth).
15. Perform a burn to achieve the phasing orbit.
    
16. Wait the appropriate number of orbits (at least one).
    

    Assuming that several phasing orbits are required, you can (and likely should) recalculate the required T value for phasing given the current orbit the active craft is in. This will allow you to compensate for inaccuracy in your initial phasing burn.

17. Once you have completed all of the required phasing orbits, recalculate and perform a burn to achieve the transfer orbit.

This works quite well at achieving flyby with moons without first aligning the planes -- the screenshots used in the "Which method of performing burns is best?" document this process in action. It is rather slow, in terms of game time, but this is rarely an issue.

It also works for setting up rendezvous between two craft that start in fairly similar orbits (say, one in a 100 x 102 km orbit and the other in a 250 x 200 km orbit) with one enormous caveat: There may be an absurd amount of relative velocity at rendezvous that needs to be canceled out.

The problem is that when you zero out relative velocity between two craft this is the same as saying "Change all the orbital parameters of the active craft to match the orbital parameters of the passive craft." This is, of course, possible -- the active craft and passive craft's orbit intersect at the rendezvous point, so there must be a burn to switch from one orbit to the other.

The problem is that the resulting burn will very likely have a large change in apoapsis, periapsis, e and h. The specific issue is the large change in e. Apse rotations cheaper when you change apoapsis or periapsis as part of the rotation, if you change both apoapsis, periapsis, and perform an apse rotation all at the same time much of this advantage disappears. This is exactly the situation the active craft will be in at the end of the non-coplanar rendezvous. In addition, it very, very likely that the transfer orbit will have a relatively high e, since it has to connect two orbits that are likely fairly circular.

As a consequence, the odds strongly favor ending up with a burn to zero out relative velocity that is heavily dominated by a radial component and is very expensive. How expensive? Several km of delta-v is not out of the question. The sky is quite literally the limit here.

This does not indicate a problem with the algorithm being used or the underlying math -- if you look at the relative velocity reported by your game when the two craft approach rendezvous, you will indeed see a relative velocity of (say) 1.2 km / second.

Thus, this method is not recommended for rendezvous between two craft in similar orbits.

Adjusting the periapsis of a hyperbolic orbit

This is another topic that the internet seems to be silent on, so this is all original work.

Changing apoapsis of a hyperbolic orbit is trivial – in fact, if you have already written code to adjust apoapsis of an elliptical orbit then you can just re-use that code without modification. By implication, you use the same code to transition from an elliptical orbit to a hyperbolic one, although you need to ensure that both apoapsis is negative and that the semi-major axis (“a”) is negative for this to work. The actual burn will occur at a true anomaly value of 0, which work perfectly well for a hyperbolic orbit.

Changing periapsis, on the other hand, requires more heroic methods – the apoapsis is a negative value, and the craft will never reach it in any case, so we need to do something else. The solution starts with constructing a tentative hyperbolic transfer orbit with the following characteristics:

• Periapsis is the desired value.
  
• Apoapsis is the value required to ensure the new orbit has the same value for “a” (the semi-major axis).
  
• e should be recalculated based on the new apoapsis and periapsis values.
  
• h should be recalculated based on the calculated e and a values.
  
• The direction of h will be the same as the original orbit.
  
• For the moment, the direction of e should be the same as the original orbit.

This is not the final transfer orbit -- there will be zero points of intersection between this orbit and the current orbit.

To fix this we need to perform an apse rotation of the tentative transfer orbit, which we can do as follows:

First, we need to select a true anomaly value for the active orbit where we want to perform the burn. This should be a point in the future (burn true anomaly should be greater than the current true anomaly) while remaining as far from periapsis as possible. I ended up selecting a point 60 seconds in the future.

Next, we need the answer to the question "Given two orbits, what are the rotations that produce at least one point of intersection between the two orbits?"

We start with this equation (from “Aligning the apse lines between two co-planar orbits”):


This is not the final form of the equation, but this works better for our current purposes. What we want to do is solve this equation for ∅, which is the rotation of the two orbits with a given true anomaly value (θ). Being rather lazy, we convert the formula into something that Wolfram Alpha will accept:

And ask Wolfram Alpha to solve for x, which produces:


After substituting back the terms and rearranging a bit, we get this:

This returns two values which have the following meanings:

1. One results in the craft still approaching periapsis (time since periapsis is a negative number).
   
2. The other results in the craft having already passed periapsis (time since periapsis is a positive number).

One way to think of it is that the second solution is the retrograde version of the first solution.

Obviously, we want the first solution. It seems like there should be a clever way to determine which one is which based on the rotation angle used, but if there is I could not find it. Instead, I recommend calculating the delta-v required to transfer to both orbits then choosing the option with the lower delta-v.

To calculate the required delta-v, we need to find the PCI velocity of the active craft at the selected true anomaly value, then find the PCI velocity on the rotated transfer orbit at "true anomaly of the active orbit" ± "apse rotation." The sign used to generate the transfer true anomaly should be the same as the sign used to rotate the transfer orbit. As mentioned above, I recommend calculating both, then selecting the option with the lowest delta-v requirements.

This is a mostly pure apse rotation operation, so the burn will be mostly in the radial direction, and could be quite expensive. It will be cheaper the lower the velocity of the active craft is, which is why it is important to place the burn true anomaly as close to the active position of the craft as possible.

Thoughts on escaping from a moon to a specified periapsis

This is really, really, really hard, far harder than you might expect. I believe that solving this problem is mostly equivalent to solving interplanetary transfers. After due consideration, I have decided not to try to solve this, at least for the initial release of the guide. However, I will present my thoughts on this topic if you want to take a crack at it.

This seems like it should be fairly straightforward -- calculate the desired orbit in the parent's sphere of influence using the moon as a proxy for the "active craft." Once you have done this, determine where the final orbit intersects with the moon's sphere of influence, then calculate the orbit that corresponds with that orbit within the moon's sphere of influence. This will be the ejection orbit, and all we need to do is enter that orbit.

Except...

The ejection orbit that is produced by this process almost certainly does not intersect with the active craft's orbit, and thus we cannot enter it. Even if it does, by chance, intersect with the current orbit, it is unlikely to intersect with it in the "right place" -- we want the periapsis of the new orbit to correspond with a position on the current active orbit that is roughly 180° away from the location where the craft leaves the sphere of influence.

Thus the problem is impossible to solve in a single impulse -- except MechJeb 2 does do this, so there must be something more here.

While there is only one ejection orbit that exactly corresponds with our desired final orbit, we are not all that picky about the parameters of the final orbit. All we need is for the periapsis to be a specified value which means that apoapsis can be changed to alter the ejection orbit. Further, inclination of the final orbit is not specified, so we can vary that as well.

So, we could try varying the parameters of the final orbit until we get an ejection orbit that we are satisfied with. Alternatively, we could start by defining the periapsis and e vector of the ejection orbit, then vary the velocity at periapsis until we get a final orbit with the desired periapsis. Either way, we are going to be doing a binary search across a sphere of influence change and it is going to be absurdly slow.

I feel that there should be a pure algebraic solution to this problem, but I cannot see how to get there and I feel that the binary search across the SOI change is going to be so slow as to be unusable.

Thoughts on interplanatary transfers

My thoughts on this topic are very similar to previous section.

First, modify the inclination of the calculated orbit of the target planet to match the inclination of current planet. While this will introduce error at rendezvous, it will be "good enough" thanks to the very large sphere's of influences associated with planets. Alternatively, you could use this simplifying assumption to generate an initial solution, then refine it with the actual target planet's orbit.

With the modified target planet's orbit, calculate a "Fast Rendezvous" between the current planet and the target planet. We want to enter this transfer orbit, so we calculate where it intersects with the sphere of influence of the current planet and produce an ejection orbit.

From here the problem is the same as the previous section -- except that there is much less flexibility in the parameters of the transfer orbit.

It is, obviously, a very hard problem.

Appendix: Projecting a vector onto a plane defined by its normal

This is a standard vector operation[www.maplesoft.com], and is defined as follows:

Where a is the vector to be projected and n is a vector that is normal to the plane onto which a should be projected.

I am unsure if this produces meaningful results if n is not normalized -- in any case you should normalize this vector to produce the expected results.

Appendix: Rotating a vector in a plane

Also known as "Rotating a vector around an arbitrary axis."

This can be done with Rodrigues' rotation formula[en.wikipedia.org]:


Where:

• v is the vector to be rotated
  
• k is a normalized vector that the v should be rotated around.
  
• ∅ is the amount of rotation desired (in radians)

The direction of a "positive rotation" is defined by the handiness of the cross product -- if it is right handed, a positive rotation will be counter-clockwise, otherwise a positive rotation will be clockwise.

Appendix: Calculating the North and East unit vectors at an arbitrary location

This section is fairly specific to Juno -- in this game, the built-in autopilot can accept either a PCI position vector or a "heading" and "pitch" values (both measured in degrees) as inputs. The heading and pitch values are measured relative to a non-inertial reference frame that I call the "North-East-Up" reference frame.

It would be helpful to be able to calculate "heading" and "pitch" values that the active craft in a particular orbit would occupy in the future. In order to do this, we must first calculate the NEU unit vectors for some point in the future. This is indeed doable, and the unit vectors should be defined as follows:



Where:

and "proj_on_plane" refers to the function defined in the appendix above.

With these unit vectors in hand, the remainder of the process (converting to spherical coordinates and then to heading and pitch) is covered in my other guide and will not be repeated here:
https://steamcommunity.com/sharedfiles/filedetails/?id=2944674093

[lucasgslarson]

I already have the impact point calculator, it’s very accurate for my craft. The reason I ask is for a boostback burn for a falcon 9, how much delta-V is needed to change the crafts impact point from down range all the way back to the Launch Site. I ask for Delta-V because I can use it to get an approximate “burn time”, which I then can use to set a heading and pitch so that when the impact point’s longitude gets to the landing site’s longitude, so does the latitude. I think this might be one of the more complicated ways of doing a boostback, but I want to learn, eh.

[mreed2]

The same process can be used on a planet with an atmosphere -- but you need to move the "Point of Impact" aimpoint downrange by an indeterminate amount, which depends upon many variables, including:
1) Angle of the velocity vector relative to the atmosphere when re-entry starts.
2) The drag of the reentering craft, including taking into account the attitude of the craft throughout re-entry.
3) The atmospheric model used.

It is possible to come up with a workable solution for a specific craft reentering on a highly similar trajectory, especially if you have active guidance during re-entry -- this is how the various Falcon-9 and Starship recreations work. But a general solution isn't possible.

[mreed2]

That's the "Calculating the point of impact" section. As stated in that section, this only works on an airless body. Calculating a point of impact on a planet with atmosphere does not have a good general solution, as the point of impact will depend on how much drag the re-entering vehicle has, which is both very vehicle specific and very implementation specific (e.g. Juno, KSP1, and KSP2 may use very different atmospheric models).

If you are trying to land at a specific spot and calculate a deorbit burn, the best solution i have to offer right now is an iterative process -- reduce the orbital velocity vectors magnitude until a point of impact exists, then rotate the velocity vector until the point of impact lines up with the orbital plane. Repeat until the point of impact is correct.

[lucasgslarson]

Is it possible to calculate how much delta-V is required to perform a burn on a suborbital trajectory, or changing the impact point of the craft to another specified spot, with this code? If so, how? Also if not with this code, do you have any idea how you could?

[mreed2]

It's hard to estimate, because alot of time was spent learning the material for myself rather than writing the guide.

All up, it was likely a month of calendar time, but i wasn't working on it full time. Maybe 40 man hours?

[KA-BOOOOM]

This is really cool ngl, how long did it take tho 💀. Big props tho its not everyday someone writes an orbital mechanics book in steam lol. Doing god's work here lol.

[mreed2]

Thanks!

[Crazy Th1ngs]

absolute legend

[mreed2]

Thanks. No issue with linking your guide to mine.

[Neobab]

Great job ! Thank you. Would you mind me adding a link to your guide in mine ?