Calculating optimal ballistic trajectories

[jofwu]

Hey folks, I could use some help/input from anyone who likes a healthy dose of math and physics!

I'm building a base on the Mun with a lab and mining/ISRU capabilities. The idea is to use a shuttle to hop around and collect science, then bring it back to the lab and fill up the tanks on the shuttle. So a common problem I have is figuring an efficient path from base to destination (and back). Right now I basically just point in whatever feels like the right direction, give myself some altitude, and burn prograde until my trajectory I comes down in the right area. I got curious about the optimal ways to do this, so I posted on Reddit to see if anyone knew of an existing tool/calculator that would be of use. Not learning of one, I thought it might be a fun exercise to give the math a shot myself.

We've begun with a simplified problem: (1) no drag/atmosphere, (2) perfectly spherical planet, (3) non-rotating planet. I'd like to try tackling 2 and 3 as well (1 is more than I want to bite off), but we'll see where it goes.

Anyways, the optimal case seems to be an elliptical orbit with initial and final locations located on the latus rectum of the second focus. This geometry minimizes the semi-major axis and thus orbital energy. I've put everything so far in a spreadsheet: https://docs.google.com/spreadsheets/d/1JIUHIAujZuQdvPRIo-FTho82QcEqVLeuRe-lLiM_b-I/edit#gid=0.

The next step I'm trying to solve is how to determine alpha, the pitch angle for the initial burn. We've got an equation that solves it for the optimal case, but I'm struggling to work out a general equation for any sma/eccentricity. With that problem solved you'll be able to input a minimum apoapsis in order to clear terrain/obstacles.

If I can get that done, I'll look next differing elevations and then a rotating body.

Just a note on terms... I've been using θ to represent the angle between initial/final location (measured from center of planet), and α to represent the pitch (flight path) angle for the burn.

If anone wants to catch up on the discussion we've had so far, here it is: https://www.reddit.com/r/KerbalAcademy/61tm0f. I wanted to move the discussion here, partly because a forum format is more useful in this case and partly because I'm hoping to get some new voices involved.

If you can contribute I'd love the help! Bonus points if you take the end results and turn them into a program or mod.

 

Edited March 28, 2017 by jofwu

[wumpus]

Well step one looks like solving the travelling salesman problem, so you will need to substitute "optimal" for "really good" after a certain number of landings/biomes (actually I think it is solvable with the number of biomes on any non-kerbin mass.

For rotating plantets, I'd check just how long any flight would be on any planet and compare that to the planetary rotation.  I doubt it would be significant on the Mun, but it might through your initial plots off (which aren't quite exact anyway thanks to being NP-complete).  For non-smooth planets, you are stuck increasing your launch angle to miss mountain ranges.  For everything else, it looks like you have the equations already.

[Armisael]

18 minutes ago, wumpus said:

Well step one looks like solving the travelling salesman problem, so you will need to substitute "optimal" for "really good" after a certain number of landings/biomes (actually I think it is solvable with the number of biomes on any non-kerbin mass.

For rotating plantets, I'd check just how long any flight would be on any planet and compare that to the planetary rotation.  I doubt it would be significant on the Mun, but it might through your initial plots off (which aren't quite exact anyway thanks to being NP-complete).  For non-smooth planets, you are stuck increasing your launch angle to miss mountain ranges.  For everything else, it looks like you have the equations already.

The assumption is just two sites - a launch site and a landing site. We've been looking at figuring out the most dv-efficient suborbital path between those two points, not pathfinding algorithms.

The problem with simply figuring out how far the landing site moves mid-flight is that the duration of the flight is a function of the position of the landing site. It seems like it ought to be possible to figure out the motion precisely, though.

You might not be going very far on a non-smooth planet - maybe you want to travel 3km and land on top of that mountain.

[wumpus]

5 minutes ago, Armisael said:

The assumption is just two sites - a launch site and a landing site. We've been looking at figuring out the most dv-efficient suborbital path between those two points, not pathfinding algorithms.

The problem with simply figuring out how far the landing site moves mid-flight is that the duration of the flight is a function of the position of the landing site. It seems like it ought to be possible to figure out the motion precisely, though.

You might not be going very far on a non-smooth planet - maybe you want to travel 3km and land on top of that mountain.

If you are finding the values with a computer, it might be a lot easier simply to iterate: ignore rotation for your first attempt, then adjust your target by how much you missed.  You should have an exact value after a small number of iterations (I suspect the rotation speed needed for this algorithm to fail are sufficient to throw you off the moon while landed).

[jofwu]

53 minutes ago, wumpus said:

travelling salesman problem

@wumpus Yeah, I'm not trying to figure out the minimum delta-v to visit each biome. This is just looking at a single trip. The travelling salesman problem of visiting biomes would be an interesting one, but not important here because I'll be returning to base and refueling each time. Honestly, looking at early results, efficiency isn't even that significant. You can add unnecessary additional height for minimal extra DV cost. But I think this is useful anyways as a "manual maneuver node" of sorts. Rather than guessing at where to point and how long to burn, you just enter two sets of coordinates, point, and shoot. Then look at the map and make necessary adjustments. Just saves from having to do lots of in-the-moment guesswork for the hop.

I was thinking the same way for when the time comes to try introducing rotation. I'd calculate the orbit first assuming no rotation to get an estimate for T. Use that T to calculate the final position of the landing site. Re-calculate to this new site and get a new T. Repeat until the difference between two T's is within some margin of error. I think that would work, and seems very doable. The hard part will be considering initial/final surface velocities and working out how this plays into determining the optimal orbit. If nothing else it makes delta-v for the two burns (launch and landing) asymmetrical.

But anyways, for now I think I want to feel good about the simple case, allow for different apoapsis, and then look at different starting/ending elevations.

@Armisael I had a thought this morning...

We have coordinates for the focus of the ellipse (planet center, set to the origin) and the two sites make up two points on the ellipse. There are an infinite number of ellipses to fit these points, but can we use that to calculate eccentricity as a function of semi-major axis? You could then take the derivative with respect to eccentricity to get the minimum sma (and corresponding e). This would help a lot for the non-optimal case and for sites at different elevation.

I found an answer to the question here (http://math.stackexchange.com/questions/547045/ellipses-given-focus-and-two-points), but I don't quite understand how to derive or apply the solution given.

Edited March 29, 2017 by jofwu

[diomedea]

The optimal hop is the one that requires less energy. With the hypothesis of perfectly spherical body and no rotation, all energy involved in the hop is the kinetic energy at launch (and conversely, at landing). Kinetic energy is function (given mass of the vessel a constant) of only speed², so the elliptic trajectory characterized by lower speed at the point of launch is your solution. As shown here, at any point in a elliptical trajectory  v² = [k/p](1+e²-2 cos(θ)) (k = MG, p = semi-latus rectum, e = eccentricity, θ = true anomaly)

[Aethon]

This may help you.

http://hopsblog-hop.blogspot.com/2014/06/travel-on-airless-worlds.html

http://hopsblog-hop.blogspot.com/2014/06/travel-on-airless-worlds-part-ii.html

[Armisael]

@jofwu Yeah, I think you're right. If you run through r = p/[1 + e cos(ν)] using both r_Ap (where ν = 180°) and r_Moon (where ν = 180° - θ/2) you get an eccentricity of [r_Ap - r_Moon]/[r_Ap - r_Moon cos(θ/2)]. From there it's pretty trivial to calculate the semi-latus rectum and SMA, and then you're basically done - we have easy expressions for any other interesting quantity in terms of those.

[jofwu]

1 hour ago, Armisael said:

@jofwu Yeah, I think you're right. If you run through r = p/[1 + e cos(ν)] using both r_Ap (where ν = 180°) and r_Moon (where ν = 180° - θ/2) you get an eccentricity of [r_Ap - r_Moon]/[r_Ap - r_Moon cos(θ/2)]. From there it's pretty trivial to calculate the semi-latus rectum and SMA, and then you're basically done - we have easy expressions for any other interesting quantity in terms of those.

The problem with this is that for different elevations at each site (r1 != r2) I don't think you can assume θ/2 is the correct angle. The points would not be mirrored across the SMA, so θ isn't perfectly bisected.

Edit: Found this: http://sketchexchange.keypress.com/forum/showthread.php?117-How-can-I-construct-an-ellipse-given-one-focus-and-two-points-on-the-ellipse&s=1257b4aa7ad910355cc8a4319036ffd7

Hard to say, but it might be useful...

@Aethon, didn't have a chance to see your links right away. This is great, thanks! The second isn't relevant, but the first explains very well how to come to the optimal solution for the basic case.

Edited March 29, 2017 by jofwu

[jofwu]

Awesome! I think I'm on the right track for getting the solution for different start/end elevations. The last link I shared confused me because there should be multiple solutions. I realized that it shows the solution where the second focus, F2, is located directly between the two points, A and B. And as I explored how you'd get other solutions everything clicked into place from what the link before was saying. The set of possible solutions for F2 forms a hyperbola where A and B are the foci.

So I worked out the math with an example, to wrap my mind around it

I drew an ugly looking moon with center of mass F1 and two points on the surface A and B. I conveniently drew them as 3-4-5/6-8-10 triangles so that I could easily locate a second focus, F2o. So the path either way from F1 to F2o is 15 giving a=7.5.

Next draw a line through A and B (which conveniently has a slope of -2). This is the x' axis for the hyperbola. The center is located directly between A and B at point M. Then draw a line perpendicular to the last at slope 1/2 through M. This is the y' axis.

Since we know the location of F2o we can calculate the hyperbola's semi-major axis as 1/2 the difference of the distances from there to A and B. You get a=2.5. And the distance from A to B is 2c, so you easily get c=5.59. c² = a² + b² gives you b=5. Now we have the equation for the hyperbola on the x'-y' axis: x'² / 2.5² - y'² / 5² = 1

If you want to vertex (on the B side of y') just use y'=0 and you get x'=2.5. I rotated and shifted back to x-y so I could plot the point, called F2v, accurately. I also eyeballed the mirror image of F2o accross the x' axis. And there's a hyperbola taking shape!

This hyperbola represents all possible solutions for the second focus of an ellipse with known focus at F1 and two points (A & B) on the ellipse. After thinking about it, I'm pretty sure F2v is the optimal solution. The distance from F1 to A and B is fixed for each. As you slide F2 either direction from the vertex the point is going to get further from both. And that means the semi-major axis of our ellipse is getting bigger. So just as with the Ra=Rb=R solution, the optimal orbit happens when the second focus is located between the two points! Pretty cool!

Here's what some of those orbits look like for my example:

Need to process the implications of this, but wanted to share...

Edit for some idle thoughts... It's neat to see how the orbits for F2o and F21 in that last image have the same energy cost. It's probably due to the concept that getting into orbit is best done with a Hohmann transfer from the surface. So there's this natural tendency for me to think that shallow is better than high. But those paths have the same energy and the steeper one sure would be easier to land.

Also realizing that if you want to have the option of trying a less than optimal SMA for these conditions, you have to deal with the issue that there's two solutions for each SMA. (one with a focus on either side of the hyperbola)

Lastly, I wonder if the other side of the parabola (the half on the other side of y') has any meaning...

Edited March 29, 2017 by jofwu

[jofwu]

I've updated my spreadsheet with what I worked out yesterday. There's a tab for the super simple case as we had it, for the record, but the main page now handles sites at different altitudes. I'm going to explain it here for the sake of anyone who wants to check it, or for myself if I ever take some time away and need to figure out what I did again.

Theta and bearing are easy equations to look up. I've converted the equation I found for bearing so that it gives degrees as 0 to 360 rather than -180 to 180, so that it's friendly with the KSP navball.

Next are some vectors inspired by this explanation. Rb is the vector from center of planet to the point at the higher altitude, Ra is the vector to the other point. Rc is the difference between these vectors and Rf is the vector from center of planet to the (optimal) second focus point, which lies along Rc. We use the geometry to work out the point that gives an equal distance (2a) from origin to each site and then to the other focus. Calculate a, c, and apoapsis for the resulting ellipse. Then it's just a matter of using the vis-viva equation to calculate the velocity at each point. Total delta-v requirement is the sum of these.

I used the dot products of the vectors above to work out the angles between site locations and the major axis. Phi (formerly alpha) is the pitch angle (more properly, flight path angle). The formula is from here: http://www.bogan.ca/orbits/kepler/orbteqtn.html. For some reason it's giving different results compared to what we had before for the simple case (see simple case sheet for comparison). But I've looked closely and I definitely think this equation is getting the right results.

Eccentric, mean anomaly, and time equations are from that last link as well. These are all just needed to calculate the travel time between the sites. I had to introduce a few conditions to keep the equations correct for some weird sign flip cases and such. I tried a few other numbers and didn't see any other problems with the results, but I could still have missed something.

I started to look at what happens when we allow an optionally higher apoapsis, but it's not proving simple. I can't find a sleek way to work out the ellipse given the apoapsis (and other info we have). I started looking at a brute force method, using the hyperbolic equation for foci, but it's looking like the equation I'm headed towards (that gives the location of the required focus) is a real monster. I'm afraid it won't have an easy solution.

Apparently Google Sheets now allows you to create functions using JavaScript. So my next attempt I think will be to write a function to guess and check until it finds the answer. I'll just pick a focus point, see what apoapsis it gives, and then make the next guess higher or lower until I find a solution within a meter of the desired apoapsis.

This should be good practice for the next step, of considering surface rotation.

My thought here would be to first calculate an orbit assuming the ground is static. I'd use calculated travel time to work out the final position of the target as it rotates. Then recalculate an orbit using the final target position. Compare travel times, adjust the final target position, and repeat. You're essentially just taking shots, measuring your travel time, and comparing to the target's movement. Eventually you'll be within some reasonable margin of error where your ship and target arrive at the same point at the same time.

The REAL trick will be accounting for surface velocity. It plays into your orbital velocity and ultimately changes the game of what's optimal. Consider a hop between two points on the equator. For a non-rotating body, the problem is symmetrical either way. But for a rotating body... If you launch in the direction of rotation, your rotational velocity decreases as you go up higher (you have to travel faster to keep up with the ground) so that the ground is passing beneath you slightly in the wrong direction. This suggests a lower orbit would be more optimal. And you can use that same problem to your advantage if you launch retrograde. You want the ground to pass beneath you and bring the target closer. So a higher hop is better. I'm not even sure what to think yet about north/south components. It's going to be tricky to see how it plays in.

[Armisael]

I don't think you're calculating the eccentricity correctly for the pitch angle in the sample case - you need to wrap a sqrt around the (1-sin(θ/2))/(1+sin(θ/2)).

The numbers you have to the right of that definitely aren't right, since the pitch angle should never go above 45 degrees (not in the simple case, at least).

[jofwu]

7 hours ago, Armisael said:

I don't think you're calculating the eccentricity correctly for the pitch angle in the sample case - you need to wrap a sqrt around the (1-sin(θ/2))/(1+sin(θ/2)).

The numbers you have to the right of that definitely aren't right, since the pitch angle should never go above 45 degrees (not in the simple case, at least).

Thanks for taking a look at that!

I  think I had a slightly different  equation for e, but regardless I've changed it now to the much simpler c/a. Tried theta for 0, 90, 180, and a few in between and it looks correct now if it wasn't before.

I was hoping that would fix my new pitch angle equation, but no. Turns out I should have been using ATAN rather than ATAN2. That change gives good answers. Undefined for 0 theta, limit of 45 as approaching 0 theta, 22.5 for 90 theta, 0 for 180 theta.  I need to look more closely at how that plays out in my new sheet though...

The old equation is giving different results between 0 and 180. I think the new equation is right on those now that I've made the tangent fix. And I like how the new one approaches 45 at 0 but is undefined at that exact angle.

Edited March 31, 2017 by jofwu

[Aethon]

On 3/30/2017 at 0:58 PM, jofwu said:

These are all just needed to calculate the travel time between the sites.

Travel time is more related to the density of the bodies.  

Mass = volume * density.

"Period of a circular orbit (T) is 2 π * sqrt(r³/(Gm)) where m is mass of planet.

Mass can be expressed as ρ * volume where ρ is density.

2 π * sqrt(r³/Gm) = 2 π * sqrt(r³/(G * 4/3 *π * r^{3 *} ρ).

The r³ cancels out and we're left with T = sqrt(π/(G * 4/3 * ρ).

Thus trip times for a given separation rely solely on density. Period scales with inverse square root of density."

[Armisael]

@Aethon That's only true if your trajectory is always exactly at sea level.

Edited April 1, 2017 by Armisael

[jofwu]

@Aethon

Not sure I understand what you're trying to say. You need those equations to calculate the time between two points along an elliptical orbit.

[jofwu]

Updates to spreadsheet.

You can now choose a minimum apoapsis (higher than optimal). The second focus falls on a hyperbola with launch/landing sites as the foci (or for same launch/landing altitudes it falls on a line). I wrote a user-defined function that takes the geometry and tries out different options for focus 2 until it narrows down one that gives the desired apoapsis. Seems to work correctly except for a bug that sometimes pops up when starting and final altitudes are the same.

Also something's wrong with the pitch/flight path angle at the landing site. Can't figure out why, because it looks perfect at launch and it's the same equation either way. But I don't see that number as useful anyways. When you want to land you're just going to burn retrograde.

Next step I'd like to take it to try out rotation, using the approach described above. But doing that essentially requires performing all of the operations on the spreadsheet repeatedly. That means I either have to pick a set number of iterations or I have to ditch the spreadsheet and program the entire thing. I have the skills to work out the guts of the program, but not to make an executable, build an interface, etc. So I think I'll cast out a net and see if anyone wants to help with that.

[surge]

function launchangMinE {
    parameter Alt, ra, Rt.

    // Optimal launch angle will vary from 45 (closest) - 0
    // (other side of planet).
    local Rr to (Alt +Re)/Rt.
    local EminA to arctan(sin(ra)/(Rr -cos(ra)))/2.
    return EminA.
}

...in kOS script is minimum energy launch angle, where:

Re: Radius of body (otherwise known as sealevel)

Alt = current ship altitude above Re

Rt = Re + targets altitude above Re

ra = range angle, the angle from you through the centre of the body to the target.

The only equation I have for range angle is related to geodesic coordinates, and probably useless to you, so you should find another way to calculate it. But just in case, I use this:

function rangeangle {
    parameter from, to. // geopositions

    local _ra to arccos(sin(from:lat) *sin(to:lat)
        +cos(from:lat) *cos(to:lat)
        *cos(to:lng -from:lng)).
    return _ra.
}

"from" and "to" should be obvious.

This is minimum energy launch angle (those yanks worked reeeeaaallly hard on this, so it's very accurate). It can never be more than 45°, never less than 0°.

Spreadsheets are a terrible way to do this, but good luck trying to do the above maths in them.

And most people just launch at 45° because you don't really lose that much energy in a no-atmosphere environment and it's safer. The only reason I bothered to learn this is because I was trying to emulate ICBMs on Kerbin where pummelling through the atmosphere at demented angles is a big energy loss.

 

Edited April 21, 2017 by surge