# Suicide Burn Code

## Recommended Posts

Kobymaru    784

Hi guys!

I wanted to update Kerbal Engineer Redux to handle "real" suicide burns, namely suicide burns that start out with a horizontal velocity component. I tried looking around the interwebs for examples or math for this, but all I ever found was "vertical-only"-suicide burns.

Then I stumbled upon this gem in MechJeb2:

```public static double SuicideBurnCountdown(Orbit orbit, VesselState vesselState, Vessel vessel)
{
if (vesselState.mainBody == null) return 0;
if (orbit.PeA > 0) return Double.PositiveInfinity;

double angleFromHorizontal = 90 - Vector3d.Angle(-vessel.srf_velocity, vesselState.up);
angleFromHorizontal = MuUtils.Clamp(angleFromHorizontal, 0, 90);
double sine = Math.Sin(angleFromHorizontal * UtilMath.Deg2Rad);
double g = vesselState.localg;
double T = vesselState.limitedMaxThrustAccel;

double effectiveDecel = 0.5 * (-2 * g * sine + Math.Sqrt((2 * g * sine) * (2 * g * sine) + 4 * (T * T - g * g)));
double decelTime = vesselState.speedSurface / effectiveDecel;

Vector3d estimatedLandingSite = vesselState.CoM + 0.5 * decelTime * vessel.srf_velocity;
double impactTime = 0;
try
{
}
catch (ArgumentException)
{
return 0;
}
return impactTime - decelTime / 2 - vesselState.time;
}```

So this piece of code here works like a charm. Except for minor errors, this provides a rock-solid suicide burn countdown timer. And obviously there is some math behind it (I mean look at all the sine's and square roots!). I don't really understand the math here though, and I haven't seen *anything* that looks remotely similar.

So instead of copy-pasting it as a black box and hoping for the best, I thought I would ask you guys, if someone could help me understanding this fine piece of math.

I looked at the git blame, and it seems that a guy named @The_Duck or @Meumeu is the original author. Are you still around?

And obviously, @sarbian is the current MechJeb maintainer, so he might know what's going on?

Any help would be appreciated.

Cheers

Edited by Kobymaru

##### Share on other sites

effectiveDecel makes no sense even after drawing a scetch.

Edited by Boris-Barboris
grammar

##### Share on other sites
Kobymaru    784
1 hour ago, Boris-Barboris said:

effectiveDecel makes no sence even after drawing a scetch.

Yeah. But it works well enough, so there's gotta be at least some sense.

##### Share on other sites
Diazo    933

Hmmm, not sure I have an answer, but see my comments in line:

```public static double SuicideBurnCountdown(Orbit orbit, VesselState vesselState, Vessel vessel)
{
//make sure we have a vessel reference
if (vesselState.mainBody == null) return 0;
//make sure our orbit actually will hit the ground. Note this check assumes the planet is a perfect sphere with ground always at zero altitdue.
if (orbit.PeA > 0) return Double.PositiveInfinity;

//how far off horizontal are we for our thrust vector?
double angleFromHorizontal = 90 - Vector3d.Angle(-vessel.srf_velocity, vesselState.up);
//sanity check that doesn't really work, see my comments below
angleFromHorizontal = MuUtils.Clamp(angleFromHorizontal, 0, 90);
//Untiy does math in radians for some reason so convert
double sine = Math.Sin(angleFromHorizontal * UtilMath.Deg2Rad);
//current force of gravity
double g = vesselState.localg;
//current available max thrust. This is calculated elsewhere.
double T = vesselState.limitedMaxThrustAccel;

//now this is a mess, see my comments below
double effectiveDecel = 0.5 * (-2 * g * sine + Math.Sqrt((2 * g * sine) * (2 * g * sine) + 4 * (T * T - g * g)));
//we know how fast we can decel, find out how long it takes us to stop
double decelTime = vesselState.speedSurface / effectiveDecel;

//take a guess at our landing site, this won't be too accurate but close enough
Vector3d estimatedLandingSite = vesselState.CoM + 0.5 * decelTime * vessel.srf_velocity;
//get height of terrain at our estimated impact site, note no check for ocean
double impactTime = 0;
try
{
//set impactTime to equal the time our *current* orbit intersects that altitude, so without a burn.
}
catch (ArgumentException)
{
return 0;
}
//actually return how long our suicide burn will take
return impactTime - decelTime / 2 - vesselState.time;
}```

First, this code makes a lot of assumptions. For just displaying a "Suicide Burn Time" number it gets away with them, but I'm not sure I'd tie this directly into a control circuit. You say it's for KER so no issue there.

Notably this code assumes we are already in a mostly vertical descent on final approach. Some of the numbers returned are "close enough" in that case, but go out of wack otherwise.

if (orbit.PeA > 0) return Double.PositiveInfinity; -> No check for mountains. Are you in a low orbit around the Mun and have the bad luck of intersecting a mountain range several KM high? Sorry, this code thinks you will miss them as it checks current orbit against sea level to see if you are going to hit the ground.

double angleFromHorizontal = 90 - Vector3d.Angle(-vessel.srf_velocity, vesselState.up); -> Assumes our engines are always point in a downwards direction even if they are pointing up. So if your engines are point 45° above the horizon, this line will think your engines are 45° below the horizon.

Vector3d estimatedLandingSite = vesselState.CoM + 0.5 * decelTime * vessel.srf_velocity; -> Guess our landing site based on our *currently available* horizontal thrust based on how far off vertical we are. Are we almost vertical on final approach? Close enough. Are we higher up with a significant sideways speed remaining with vessel off vertical? Will be way out of wack. Player doesn't see this so probably not an issue.

Now for the kicker, and I can't fully parse this line:

double effectiveDecel = 0.5 * (-2 * g * sine + Math.Sqrt((2 * g * sine) * (2 * g * sine) + 4 * (T * T - g * g)));

Breaking it down, but I'm not sure of my logic here.

First, get rid of all the numbers, they may be needed for accuracy in the math but for figuring the logic out we can ignore them, so:

double effectiveDecel =  -1 * g * sine + Math.Sqrt((g * sine) * (g * sine) + (T * T - g * g)));

So inside the square root we are left with two pythagoreas calcs. The first is our available thrust after accounting for gravity: t^2 = T^2 - g^2

Then the second is the g*sine, so x^2 = (g*sine)^2 + t^2

That is our effective thrust, so we then subtract gravity (-1 * g *sine) from it again?

That makes no sense.

It's also the best explanation I can come up with so I'm stumped.

D.

##### Share on other sites
sarbian    5000

That code predates me

##### Share on other sites
Kobymaru    784
Quote

First, this code makes a lot of assumptions. [...] You say it's for KER so no issue there.

It is for KER, but that doesn't mean we can't deconstruct the assumptions and make it more precise as we go.

Quote

double angleFromHorizontal = 90 - Vector3d.Angle(-vessel.srf_velocity, vesselState.up); -> Assumes our engines are always point in a downwards direction even if they are pointing up. So if your engines are point 45° above the horizon, this line will think your engines are 45° below the horizon.

I am not sure about that. I think vesselState.up is not the vessels "up" direction, it's the direction of "up" in general, away from the current main bodies center of gravity.

Quote

It's also the best explanation I can come up with so I'm stumped.

That's alright, apparently everyone is. Thanks for your insight!

1 hour ago, sarbian said:

That code predates me

Interesting. I guess the hunt for the truth begins!

Edited by Kobymaru

##### Share on other sites
Poodmund    785

The fact that your acceleration and mass are changing constantly, its very difficult to calculate this way.

Coming from a different angle, would there be any point in raycasting from the CoT along the thrust vector to see if it hits the currently orbited body and readout the distance to surface and angle of incidence on the fly? Is that even possible?

##### Share on other sites
flywlyx    104
8 hours ago, Diazo said:

That is our effective thrust, so we then subtract gravity (-1 * g *sine) from it again?

That is the gravity during your suicide burn, integral of cosine *g.

##### Share on other sites
diomedea    695

I'm not sure why the solution in MJ2 code, it spells to me about computing the root values for time in an equation derived from newtonian mechanics. However there are a number of issues with such approach, mainly due to values for gravity being a function of altitude, and vessel acceleration a function of vessel mass. To achieve a better solution for the powered descent problem (suicide burn being the limit with a powered descent when no margin exists), it is possible to run a simulation in time of the descent (e.g. as done here).

However, I believe a simpler solution exists based on orbital energy (though haven't yet found it applied to this problem anywhere). Let me consider just the vertical components for sake of simplicity. Given a specific initial altitude (zi) and initial vertical speed component (vi) of a craft with mass m, is possible to easily compute the total initial orbital energy Ei = m * (1/2 *vi^2 - Gk*M/(R+zi)), (Gk = universal gravitational constant = 6.67408E-11, M = mass of mainbody, R = radius of mainbody). Of course we also know the final energy at altitude = 0 and speed = 0 (apart from the horizontal speed component with a rotating body that we can dismiss for now), Ef = - Gk M m /R.

Now the difference in energy Ef - Ei is exactly what we need to provide to have the craft stop at altitude 0. Such energy difference is also known as Work and is easily applied when a constant Force (Thrust in our case) is applied along the direction the vessel moves: W = F * (zb - zf) (zb being the altitude where 100% thrust must be applied, or burn to begin). That way is immediate to find zb from Work; having zb and using newtonian mechanic the speed at burn start vb = RADQ(2/m*(Ei - potential energy at zb))  = RADQ(2/m*(Ei +GkMm/(R+zb)). With vb is then easy to compute time of burn tb = vb / (Thrust/m - GkM/(R+zb)^2). Of course, if Time to Impact is already known (simple newtonian mechanic again while no thrust is applied), subtracting the time of burn tb gives the time left before the suicide burn. If solving for horizontal as well, there is need to consider thrust applied retrograde, craft pitch changes while horizontal speed is matched to the body surface speed and thrust is applied by sine and cosine of pitch to the vertical and horizontal portions of the problem. Again a simple solution goes with energy and work, the total work being the composition of horizontal work required to match horizontal speed and vertical work as above.

##### Share on other sites
Kobymaru    784

Hi, a little update on my endevour:

I have implemented the "MechJeb2-version" into KerbalEngineer in this branch:

The timer works reasonably well, but I had to change the definition of Altitude and Distance

Altitude: "altitude when to start a suicide burn" -> "Altitude lost during burn until velocity reaches 0"
Distance: "distance to the point at which to start a suicide burn." -> "Distance above ground after a suicide burn"

Note that both definitions are equivalent at the time of the start of the burn - but only then, not before. Turns out it's pretty hard to calculate the Altitude lost.

I am also currently preparing a different approach: there is a method of calculating the velocity and time of a gravity turn as a function of the angle, for constant TWR:

There is a closed form for the burn duration from the initial conditions, but unfortunately not closed form for altitude lost. So I decided to implement that algorithm in KerbalEngineer, and do a little numerical integration. Maybe this will be a bit costly in terms of performance, but I really really want proper suicide burn info  I'll keep you guys posted on my progress, and when I'm done, I'd love to submit my code for "peer review"

@diomedeaI'm afraid I do not fully understand your approach, but I see a bit of an issue here:

On 27.1.2017 at 3:25 PM, diomedea said:

Let me consider just the vertical components for sake of simplicity.

On 27.1.2017 at 3:25 PM, diomedea said:

If solving for horizontal as well, there is need to consider thrust applied retrograde, craft pitch changes while horizontal speed is matched to the body surface speed and thrust is applied by sine and cosine of pitch to the vertical and horizontal portions of the problem.

On 27.1.2017 at 3:25 PM, diomedea said:

Again a simple solution goes with energy and work, the total work being the composition of horizontal work required to match horizontal speed and vertical work as above.

The vertical-only problem is already fully solved, and in fact is already implemented in KerbalEngineer (and keeps crashing Scott Manley because he's not aware that it's vertical-only ).

The horizontal component is where the money is at, beause that's what's required for entry from orbit, and for minimal dV usage. Now I believe that the complicated part lies in "consider thrust applied retrograde" and in "craft pitch changes while horizontal speed is matched to the body surface speed and thrust is applied by sine and cosine of pitch to the vertical and horizontal portions of the problem."

For infinite TWR, everything is very simple: the Energy we use is: Energy of vertical component, Energy of horizontal component and Potential Energy of altitude. So far so good. But for finite TWR, we have to consider gravity losses: while we are killing the horizontal velocity, we need to keep our head afloat to not crash into the ground. For TWR -> 1, the work that needs to be expended tends towards Infinity. Somewhere in between lies the actual value, but for sure, the Energy expenditure is not at all constant.

To put it in your terms:

"Work [...] is applied along the direction the vessel moves: W = F * (zb - zf)" - in the 1-D case, everything is clear. But in the 2-D case, the thrust is applied along a curve, and the shape and length of the curve is dependent on the intial speed, the angle to vertical and the TWR.

Please do correct me if I'm wrong! If you can, I would be glad if you could flesh out your idea to include horizontal+vertical. If you could package it in a nice, simple formula for altitude and burn duration as a function of velocity and TWR, you would make me a very happy man

##### Share on other sites
diomedea    695

@Kobymaru, let me start with some considerations. You have implemented a MechJeb2-like algorithm. Believe you know the equation it embeds (as you showed in the OP) is analytical, unaware of changes in gravity and TWR during the descent. Both me and you linked in the posts above two different approaches for a numerical computation, using integration in time to handle changing both gravity and TWR.

The reason I showed an approach based on energy instead was to still be able to solve this problem with an analytical approach. At least, it allows that in regard to gravity (as potential energy is computed considering the different gravity at altitude and landing); I have no solution for a change in vessel mass, therefore TWR (so, my approach works best if the mass of fuel burnt during descent is a small fraction of the total vessel mass).

Before getting into more equations, let me show one consideration about "gravity losses" that is often lost. While almost everybody seems to agree about the Oberth effect, which happens when thrust is applied to increment speed (so, energy changes the more the current speed is, going by the dEk/dt = F * v equation), very few seem to make that the change in energy required to bring a suborbital in a orbital trajectory is also described in the same guise. The loss in speed known as "gravity losses" isn't but a Oberth effect in reverse (less then enough speed to achieve circular orbit). Exactly as the Oberth effect, it is totally explained when dealing with energy, instead of with speed. Therefore, the reason you showed requiring to compute gravity losses because of a finite burntime, doesn't exist at all when computing energy instead: we only need to know the amount of energy we have at start and the amount we want at end of the maneuver, how that energy is changed (Thrust/m = acceleration) makes the time of burn change but not the energy equation.

Another important consideration. Energy is a scalar, which means by itself can't show the correct direction for a speed change. On the vertical direction (exposed in my previous post), the same total energy applies to a same vessel that, at one specific altitude, is climbing or is descending by the same absolute vertical speed. Of course, without atmospheric drag, a climbing vessel will sooner or later (barring exiting SOI) be descending at exactly the same vertical speed at the same altitude, so the equations don't really need consider the case (we compute burntime only with the descending part of this trajectory). When dealing with only horizontal components, we don't need to include potential but only kinetic energy; however to solve the problem requires to ensure the final horizontal speed is oriented as the radial component due to the body rotation (meaning, final horizontal speed has to exist only in the E/W component, positive with E, parallel to the equator, and be exactly = 2*π* radius* cos(latitude)/(rotation period)). The above "final horizontal speed" can clearly be used to provide the "final kinetic energy due to horizontal speed", which is the value we need to achieve (but, due to the scalar nature, we can't just subtract final from initial kinetic energy, we need to consider the whole change in energy, or "work", performed). It seems easy to consider that, however steep or mild the descent profile, all changes in energy in both N/S and E/W components need be completed before (not necessarily at) impact time. "Before" meaning the last part of the descent will be pretty vertical, while "after" is more alike what a plane does at landing.

Now, different coders have different approaches to solve this kind of problems. I like to handle calculations separate on each of 3D components; however the vectorial sum of required thrust when time is equal on all axes gives exactly the retrograde direction (which is neat). Components wise, total thrust would be made on each axis: for vertical, effective thrust = Thrust * sin(pitch); for E/W, effective thrust = Thrust *cos(pitch)*sin(heading); for N/S, effective thrust = Thrust *cos(pitch)*cos(heading). Of course effective thrust is what would be used in each of the components equations, e.g. for E/W: tbE/W = ΔvE/W / (effective thrust E/W / m). Given tb should be the same (at least with both horizontal components) means the heading is invariant during descent; if pitch was invariant too, tb would be the same in all three components (that is what would happen if the trajectory before the descent burn was already suborbital, crossing the body surface exactly at the landing site).

This last consideration shows (at least to me) why changing pitch during descent isn't really an issue: the end result still is to achieve the final energy state at the right time, the descent can mix a late deorbiting burn with the landing burn to perform a "reversed gravity turn" which is pretty efficient, but the most efficient of all trajectories would require the deorbiting burn be done well before so to have an almost flat trajectory during descent (very risky, however). Just, exactly, the reverse profile of a gravity turn launch trajectory from an airless body. Having pitch change during descent only makes our "gravity losses" more prominent, but when  we so choose, the effect is actually to have a longer burn time in horizontal than with vertical.

Now, with all the above shown, I'm in doubt if this method actually fits your needs. No TWR variance is considered, and method works best with flat descents instead of allowing pitching (but, should either pitch or thrust be actively changed during descent, all precomputed analytical solutions would be invalid and need be recomputed in real time). In the end, what makes the MechJeb2 approach work in practice is just this, even if not totally accurate it is kept updated in real time (and, given TWR increases due to burnt fuel, the resulting error is towards safety).

##### Share on other sites
linuxgurugamer    5567
On 3/24/2017 at 1:44 PM, Kobymaru said:

I have implemented the "MechJeb2-version" into KerbalEngineer in this branch:

The timer works reasonably well, but I had to change the definition of Altitude and Distance

Altitude: "altitude when to start a suicide burn" -> "Altitude lost during burn until velocity reaches 0"
Distance: "distance to the point at which to start a suicide burn." -> "Distance above ground after a suicide burn"

Note that both definitions are equivalent at the time of the start of the burn - but only then, not before. Turns out it's pretty hard to calculate the Altitude lost.

I am also currently preparing a different approach: there is a method of calculating the velocity and time of a gravity turn as a function of the angle, for constant TWR:

There is a closed form for the burn duration from the initial conditions, but unfortunately not closed form for altitude lost. So I decided to implement that algorithm in KerbalEngineer, and do a little numerical integration. Maybe this will be a bit costly in terms of performance, but I really really want proper suicide burn info  I'll keep you guys posted on my progress, and when I'm done, I'd love to submit my code for "peer review"

Wondering how this is going?  I'm working on a Suicide Burn mod, and would prefer to not have to reinvent the wheel.

My mod is working with BetterBurnTime, but due to a licensing issue, I won't be able to use that when I release.  So I was wondering if I could help you with this, to both of our benefits.

##### Share on other sites
Kobymaru    784
1 minute ago, linuxgurugamer said:

Wondering how this is going?

Meh

So I did implement the Gravity Turn algorithm in reverse. It "works" in the sense that the numeric values come out to be what is expected by solving the Differential Equations in [1] numerically.

It does not work particularly well in practice, because the Suicide Burn altitude is greatly overestimated. I believe it's partially because the "flat earth" assumption breaks the centrifugal force. However, I tried to account for the centrifugal force in another model [2] and still came up short (or rather too high  ) Right now I'm not quite sure where to go from here.

I recall that the timer itself being pretty much OK though.

You can check out my progress on the "suicide-numeric" branch of KerbalEngineer here:

Like I said, it's not quite there yet

BTW, I have Mathematica Notebooks with 2 models of suicide Burns, I can post these if there is interest. Papers "available" on sci-hub.io, or from me on request. What I also have is a collection of papers on this subject that may or may not be useful - I can't even tell yet.

1 minute ago, linuxgurugamer said:

I'm working on a Suicide Burn mod, and would prefer to not have to reinvent the wheel.

My mod is working with BetterBurnTime, but due to a licensing issue, I won't be able to use that when I release.  So I was wondering if I could help you with this, to both of our benefits.

Oh that would be wonderful. I'm not really doing this for KerbalEngineer, I just chose it because of the nice framework it provides.  What I reall want are functioning and precise Suicide Burn aids. If it's in a separate Mod, so be it.

[1] Culler, G. J., & Fried, B. D. (1957). Universal Gravity Turn Trajectories. Journal of Applied Physics, 28(6), 672–676. https://doi.org/10.1063/1.1722828
[2] McInnes, C. R. (2003). Gravity-Turn Descent from Low Circular Orbit Conditions. Journal of Guidance, Control, and Dynamics, 26(1), 183–185. https://doi.org/10.2514/2.5033

##### Share on other sites
Pand5461    91

I can guess why no one is even bothering about suicide burns that start with significant horizontal component. The difficulty is to keep the direction of the burn that won't be purely retrograde in such case.

As for SQRT in time to burn computations - that might come from the approximate solution of the Tsiolkovsky's equation with gravity, I have to compare it against the formulae I obtained while writing kOS landing autopilot.

##### Share on other sites
Kobymaru    784
1 minute ago, Pand5461 said:

I can guess why no one is even bothering about suicide burns that start with significant horizontal component. The difficulty is to keep the direction of the burn that won't be purely retrograde in such case.

Since we got really nice SAS functions including Retrograde hold in Stock KSP, I don't think this is such a big issue anymore. After digging around in ways to actually calculate the burn, I would say the difficulty lies in the implementation side

##### Share on other sites
linuxgurugamer    5567
11 minutes ago, Kobymaru said:

It does not work particularly well in practice, because the Suicide Burn altitude is greatly overestimated. I believe it's partially because the "flat earth" assumption breaks the centrifugal force. However, I tried to account for the centrifugal force in another model [2] and still came up short (or rather too high  ) Right now I'm not quite sure where to go from her

That may be because the calculation of the time to start the burn may not take into account the fact that the vessel is slowing down.  Ignoring orbital speed for the moment, look at it this way:

Time to impact is calculated based on the vessel accelerating due to gravity

Assumptions:

Vessel is moving at 100 m/s

Engine thrust gives an acceleration of 20m/s in space.

Gravity = 10m/s/s

Net deceleration at full thrust = 10m/s/s

So, at an altitude of 1500, you have a time to impact of 10 seconds.

Now, if you start burning at 1500m, you will stop in 10 seconds, at an altitude of about 1000m.

This is getting in calculcus, which I was never good and and have no memory of anymore.  But, with the proper calculations, you should be able to figure the correct altitude to start burning at, which is going to be somewhere below 1000m.

Any math geniuses around?

##### Share on other sites
Kobymaru    784
4 minutes ago, linuxgurugamer said:

That may be because the calculation of the time to start the burn may not take into account the fact that the vessel is slowing down.

Nope, that should all be accounted for.

##### Share on other sites
linuxgurugamer    5567

Thinking about it, add a calculation of how far the vessel will travel during the entire burn.  Then, compare that to the altitude, if not below a certain threshold, continue falling.

For example, in the example above, starting the burn at 1500, the vessel will travel 500m.  At an altitude of 1500, it's obviously too early to start burning.

So,continue falling.  If we have a value we are aiming for, say, 0 velocity at 10m in height, then keep falling until the altitude - decelDistance <=10

##### Share on other sites
linuxgurugamer    5567

So, to do this, I think I need the following, can you provide it:

double timeUntilImpact()

double distanceToImpact() // Take angular velocity into account, if vertical, then should equal altitude

I already have the calculations for the vessel speed and engine thrust, etc.

##### Share on other sites
Pand5461    91
37 minutes ago, Kobymaru said:

Since we got really nice SAS functions including Retrograde hold in Stock KSP, I don't think this is such a big issue anymore. After digging around in ways to actually calculate the burn, I would say the difficulty lies in the implementation side

Erm... Holding Retrograde at full thrust most probably won't work - you won't be able to bleed off horizontal and vertical speed at the same moment. You are going to need either fairly complicated steering or throttling. Stock SAS doesn't do either of that automatically.

I can post my continuous-burn landing script somewhere. It is not exactly this exact problem but might give some ideas.

##### Share on other sites
Kobymaru    784

24 minutes ago, linuxgurugamer said:

So, to do this, I think I need the following, can you provide it:

double timeUntilImpact()

double distanceToImpact() // Take angular velocity into account, if vertical, then should equal altitude

This is the hard part. At the moment, I can't provide it accurately. My current model is so inaccurate that it's almost useless.

17 minutes ago, Pand5461 said:

Erm... Holding Retrograde at full thrust most probably won't work - you won't be able to bleed off horizontal and vertical speed at the same moment. Stock SAS doesn't do either of that automatically.

When was the last time you tried? I do this on a regular basis. Set NavBall to surface mode, press retrograde button, press Z and wait until your vehicle has killed all velocity.

Quote

I can post my continuous-burn landing script somewhere. It is not exactly this exact problem but might give some ideas.

If it's not too much effort, that would be nice.

##### Share on other sites
linuxgurugamer    5567

I just realized, that the work has already been done, and it's under an MIT license.

Mod is Trajectories, github repo is here:  https://github.com/neuoy/KSPTrajectories

Should have everything I need to get those two values.

##### Share on other sites
Posted (edited)

It's not a planar problem. I'd suggest basic shooting method over burn start time (Tstart), RK4 or some other NI scheme, applied on the system:

// x - 3d world-space position vector. Note that transform.position in KSP is noninertial surface-relative space, so correct the forces if you want to work in Unity space.
// v_whatever - 3d velocities
// acc - 3d acceleration

acc(t) = gravity(x) + aero(x, t, v_surf) / mass(t)                                                                                // engines off, t < Tstart. If no aero, use vessel.Orbit to get x(Tstart) and v_orbital(Tstart) analitically.
acc(t) = - v_orbital(t).normalized * max_thrust(x) / mass(t) + gravity(x) + aero(x, t, v_surf) / mass(t)     // SAS in orbital retro mode, t >= Tstart
acc(t) = - v_surf(t, x).normalized * max_thrust(x) / mass(t) + gravity(x) + aero(x, t, v_surf) / mass(t)   // SAS in surface retro mode, idk exactly when the switch happens, ask devs or sekrit documents
d(mass(t)/dt) = magic_vessel_function(x, t)                                                                                  // vessel-specific
v_surf(t, x) = v_orbital(t) - cross(planet_ang_vel_vector, (x - planet_center))
d(v_orbital(t))/dt = acc(t)
d(x(t))/dt = v_orbital(t)
v_surf(Tend, x_end) = zero_vector   // Tend - time of full stop. Criteria: t = Tend if v_surf(t, x).length < Eps, where you can choose Eps like, for example: time_step * max_thrust / mass(t) * 0.7424242...
planet_ang_vel_vector = planet_ang_vel * north_vector
(x - planet_center).length < terrain_height(x, t)         // you crashed, decrease Tstart. Check this condition during NI, not necessarily every step.
(x(Tend) - planet_center).length > terrain_height(x, t)   // undershoot, increase Tstart
dot(v_surf(t, x), (x - planet_center)) > 0                      // you're literally flying away.

v_orbital - 3 variables
x - 3 variables
mass - 1 variable
you'll integrate 7 state variables simultaneously. All others are either constant, or can be calculated from them directly.

Substract a couple of seconds from Tstart to form a safety margin.

Edited by Boris-Barboris

##### Share on other sites
Kobymaru    784
Posted (edited)
2 hours ago, linuxgurugamer said:

I just realized, that the work has already been done, and it's under an MIT license.

Mod is Trajectories, github repo is here:  https://github.com/neuoy/KSPTrajectories

Should have everything I need to get those two values.

Incidentally, it's me who got stock with maintaining that mod so I know at least a little bit about it. Care to enlighten me how it calculates the required values? Because Trajectories doesn't calculate the trajectory with an engine burn, it calculates it under free fall or atmospheric. How can it help you for suicide burn?

1 hour ago, Boris-Barboris said:

It's not a planar problem.

Why not? I don't see the 3rd dimension. I see only vertical and longitudinal

Quote

I'd suggest basic shooting method over burn start time (Tstart), RK4 or some other NI scheme, applied on the system:

Thanks for the code! Full numeric solution was my last resort, I was hoping for "half-analytical".

Especially this looks like a nasty bit of coding that I don't know how to do nicely yet:

1 hour ago, Boris-Barboris said:

Edited by Kobymaru

##### Share on other sites
Just now, Kobymaru said:

Why not? I don't see the 3rd dimension. I see only vertical and longitudinal

Planet rotation. ~170m\s on equator for Kerbin. Now imagine a landing from polar orbit.