Jump to content

Suicide Burn Code


Recommended Posts

Just now, Kobymaru said:

Incidentally, it's me who got stock with maintaining that mode so I know at least a little bit. 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?

Yes, I saw that after I posted.

Trajectories calculates free-fall trajectories, with that info I can do the suicide burn calculations.

Where is your github repo, there are some additions I'd like to add to the mod

 

Link to comment
Share on other sites

1 hour ago, Boris-Barboris said:

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.
Goes without saying to spread the load on multiple frames.

Vector3D x = vessel.GetWorldPos3D();

Double FlightGlobals.ActiveVessel.srfSpeed;

What are the names with (t), (x), etc?

Link to comment
Share on other sites

4 minutes ago, Boris-Barboris said:

explicit declaration that it's a function of at least time or position.

So, this:

acc(t) = gravity(x) + aero(x, t, v_surf) / mass(t)  

would mean (in english):  Acceleration at time "t" = gravity(vessel.3dposition) + aero(vessel.3dposition, "t", vessel.surfaceVelocity) / mass("t")

where mass(t) = the mass of the vessel at that point in time

Correct?

Also, what is "d" as in d(t)

Edited by linuxgurugamer
Link to comment
Share on other sites

3 minutes ago, linuxgurugamer said:

would mean (in english):  Acceleration at time "t" = gravity(vessel.3dposition) + aero(vessel.3dposition, "t", vessel.surfaceVelocity) / mass("t")

Acceleration at time "t" = gravity(vessel.3dposition, ... unimportant known stuff like planet position and gravity constants) + SumOfAllAeroForces(vessel.3dposition that is needed to get temperature, pressure etc., at the moment of time "t", vessel.surfaceVelocity to calculate forces, ... lots of stock aero-related things) / mass(at the moment of time "t")
 

6 minutes ago, linuxgurugamer said:

Also, what is "d" as in d(t)

differential.

Don't look too much into the symbols, it's not like they are supposed to be the exact formulas. Just a sketch of the solution.

Link to comment
Share on other sites

To try to solve this problem in a kOS script, I eventually gave up on an analytical answer and just went with a brute force numeric approximation by running a timeslice simulation loop, like so:

// (Pseudo code here !! Not real - just using code snippet for indent formatting.)

Keep a set of state variables that will change during the loop:
  XYZ position of ship,
  ship velocity (velocity in surface reference frame, NOT orbital).
  ship mass.

Variables that will stay constant:
  set dT to time delta per loop iteration (I went with 0.5 s.  This
    is a compromise because too coarse-grain and the answer is inaccurate,
    but too fine-grain and the answer takes too long to obtain and you
    need to know the answer fast in real-time to make the decision when
    to burn.)

Initialize state variables to current situation

Loop:
  set a_grav = (scalar) How much gravity accel do we have at this position( i.e. get radius
     to body center by vector subtraction of position and body position, and
     use that with the body's MU to calculate the exact 'g' here.)
  set F_engine = (scalar) How much thrust do we get at max throttle?
  set a_engine = (scalar) How much accelleration does F_engine give? (divide by current mass)
  set down_unit_vec = ((vector)body position) - ((vector)ship position).normalized
  set pointing_unit_vec = (-1) * ((vector)velocity).Normalized // assumes flying butt-first.

  // total accelleration this timestep is vector addition of gravity effect plus engine effect:
  set a_total = (vector) (a_grav*down_unit_vec + a_engine*pointing_unit_vec)

  set velocity = old velocity + a_total/dT.
  set position = old position + velocity/dT.

  // See if this is the point where it would stop and go upward if you kept thrusting,
  // using dot products to decide if the old velocity and the new velocity are opposing directions:
  if vector_dot(velocity, (old velocity)) < 0
    break from loop.
    
  set old velocity = velocity.
  set old position = position.
end loop.

set suicide_margin = position's altitude - terrain's altitude directly under position

If the above loop calculates suicide marge as positive, then you don't have to burn yet.  If it calculates it as negative, you should have already started burning and you're dead.  If it calculates it as exactly zero, that's exactly when to start burning.

Obviously you need a bit of wiggle room because the loop is running in realtime while the game is progressing so by the time it gives you the answer of zero it's too late :-), so you can fudge it a bit by learning how many milliseconds the loop iterations seem to be taking, and calculate ahead for what the position *would be* if you didn't start the burn until that many milliseconds from now, making the first iteration of the loop have a predicted engine thrust of zero and not start it until the second iteration...

 

Link to comment
Share on other sites

10 hours ago, linuxgurugamer said:

Ok, so I'll fork the main and move from there

Ok. What did you want to change anyway?

BTW, I would prefer to continue the discussion about Trajectories in the trajectories thread: 

 

Edited by Kobymaru
Link to comment
Share on other sites

On ‎26‎.‎01‎.‎2017 at 1:32 PM, Kobymaru said:

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.

 

Any help would be appreciated.

Cheers

1. The way from a vertical to a diagonal suecide burn is relative trivial, as you just take the Integration of rocket equation and then you have the distance nedded to come to a stop. (doesn't matter which direction). the only real difference are cosine losses of your engine against the gravity, which are 0 ich you are horizontal, and 1 when vertical.

2. You really need to stop the surface velocity and not the orbital velocity.

3. The coriolis force is a poodle and will ruin your day. (and spaceship)

Edited by Ger_space
spelling
Link to comment
Share on other sites

35 minutes ago, Ger_space said:

The way from a vertiacal to a diagonal suecide burn is relative trivil,

The authors of about 6 scientific papers that I've read on this topic tend disagree with you

 

Quote

as you just take the Integration of rocket equation and then you have the distance nedded to come to a stop.

You can't "just take the Integration", because sometimes the Integration doesn't even have a closed form (one that can be written as a formula).

 

Quote

(doesn't matter which direction).

It matters in which direction, because ...

 

Quote

the only real difference are cosine losses of your engine against the gravity, which are 0 ich you are horizontal, and 1 when vertical.

That's cool and all, but what about the cosine losses in between horizontal and vertical?

The optimal path is not a) a straight line or b) a quarter-circle or a cosine or any simple form. The optimal path is a strange curve that looks approximately like SQRT but isn't. That's what it looks like:

7xQdroR.png

Along this path, the cosine loss depends on your radial velocity around the body (because of the centrifugal force) and your current velocity angle. Your current velocity and angle depends on how much you have slowed down on the path before, and on the previous cosine loss. Your previous cosine loss depends on the previous radial velocity around the body and the previous velocity angle. And those depend on the cosine loss before that...

This kind of recursive relationship can be expressed mathematically as a system of coupled differential equations. In its simplest form, these are:

dv/dt == g * (TWR - cos(beta))
v * d beta/dt == g* sin(beta)

This already is hard to solve, but possible (only for angle and velocity - not for altitude/downrange distance!), but mind you: this assumes 

  • Flat Earth
  • Constant g
  • Constant TWR
  • No Rotation

Which obviously breaks in multiple ways in KSP. Now any attempt to expand the model a bit to account for a spherical body, changing g, changing TWR, and rotation, you end up with shenanigans like this one:

MPYplRv.png

Containing mathematical functions that I haven't even heard of. Mind you, even these monstrosities come with a huge list of assumptions like "but only if the velocity is much lower than orbital velocity and rotation is low and it's a monday evening and it's full moon".

Believe me guys, I *wish* it were as simple as "just add the second dimension". And I have read a lot of comments and a lot of threads that say "oh that should be trivial". But if you look at the physics, and if you try it out in practise you realize that it's not at all simple.

edit: Click here for more fun stuff

Spoiler

5omhaUH.png

 

Edited by Kobymaru
More fun stuff
Link to comment
Share on other sites

39 minutes ago, Ger_space said:

2. You really need to stop the surface velocity and not the orbital velocity.

3. The coriolis force is a poodle and will ruin your day. (and spaceship)

The coriolis force doesn't exist. Surface velocity doesn't exist. They're all just artifacts of a rotating reference frame. Depending on the problem, sometimes it's easier to pretend that they do exist. I think this is not one of those cases. Or maybe it is. Who knows, I'm not a physicist.

Link to comment
Share on other sites

20 hours ago, Kobymaru said:

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.

On second thought, you are right, it must be possible.

So, here's the kOS lander code: https://drive.google.com/open?id=0Bw2BzA4x17Y2Tl9OcW5GNkVkcDA. It does landing by trying to maintain constant horizontal and vertical acceleration.

So, on the problem in OP.

As far as I understand, the number needn't be very exact unless it's close to zero. So, you may integrate retrograde burn until all velocity is killed, then assume distance traveled through burn proportional to vertical speed and estimate when to do suicide burn by that condition.

The integration may be done by, say, Verlet scheme because it's computationally cheap and accurate enough, to my experience. Also, timestep of about 2 seconds should be enough, except for final moments of braking.

Some concerns:

  1. Equations of motion must either be done in inertial frame or centrifugal and Coriolis forces must be accounted for.
  2. Unlike vertical burn, the ground height at landing site is generally not the same as right beneath the ship.
  3. In the case of vertical burn, if you missed the right moment - you are certainly going to crash. When you have initial horizontal speed, there are two cases - one is when you can't kill all vertical speed at touching the ground even burning straight up all the time, the other when you don't have enough time to kill horizontal velocity before touching the ground. In the former case, you are doomed. In the latter case, the situation can still be fixed by burning a bit upwards from retrograde. Thus, the time shown may be quite a bit misleading.
Link to comment
Share on other sites

2 hours ago, Kobymaru said:

The coriolis force doesn't exist. Surface velocity doesn't exist. They're all just artifacts of a rotating reference frame. Depending on the problem, sometimes it's easier to pretend that they do exist. I think this is not one of those cases. Or maybe it is. Who knows, I'm not a physicist.

you pick up horizontal speed, when you drop vertical on a planet, which will flip you over when you dont compensate vor it..

this are the formulas from my kOS scripts.

	

	ch_rate = mass burned per second
	v_e  = Engine exhaus velocity
	M = SHIP:MASS.
	v_burn = SHIP:VELOCITY:SURFACE:MAG.

	  t = burn time needed to stop: M*(1 - e^(-dV/v_e))/ch_rate.

stop_distance_horizontal = v_burn*t - (v_e*(t - M/ch_rate) *ln(M/(M-t*ch_rate)) + v_e*t) .



here you need an other t, which cannot be simple computed, so you need an interatice solution. (or you have a solver for a lambert_W0 function.)

	until dv_change < 0.001 {
 
		set dv_total to dv_total + dv_change.
		set t to SHIP:MASS*(1 - e^(-dv_total/v_e))/ch_rate.
		set dv_new to g_mean * t.
		set dv_change to abs(dv_new - dv_old).
		set dv_old to dv_new.
	
	}


stopping_distance_versical = (grav_mean_b*t^2)/2  + v_burn*t  - (v_e*(t - M/ch_rate) *ln(M/(M-t*ch_rate)) + v_e*t) .

(grav_mean_b*t^2)/2 is the the gravitation loss. 

this only works for vacuum. If you have an athomsphere... good luck, because the engine will have different stats depending on your air pressure and the ship will be deaccelerated by the air.

Link to comment
Share on other sites

On ‎1‎/‎26‎/‎2017 at 4:32 AM, Kobymaru said:

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 terrainRadius = vesselState.mainBody.Radius + vesselState.mainBody.TerrainAltitude(estimatedLandingSite);
	double impactTime = 0;
	try
	{
		impactTime = orbit.NextTimeOfRadius(vesselState.time, terrainRadius);
	}
	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

"Real" suicide burns?

Link to comment
Share on other sites

9 minutes ago, The Moose In Your House said:

"Real" suicide burns?

Suicide burns with a horizontal component.

When you want to land on a body, you usually come from orbit with a high horizontal velocity. The most efficient way to land from there is a real suicide burn, like a reverse gravity turn. "Stopping in mid air" and then dropping to the ground is significantly less efficient and isn't even a proper suicide burn. That's why I consider the current implementation of "suicide burn" not very useful.

Edited by Kobymaru
Link to comment
Share on other sites

Just now, Kobymaru said:

Suicide burns with a horizontal component.

When you want to land on a body, you usually come from orbit with a high horizontal velocity. The most efficient way to land from there is a suicide burn. "Stopping in mid air" and then dropping to the ground is significantly less efficient. That's why I consider the current implementation of "suicide burn" not very useful.

Can you do a suicide burn on any planet/moon?

Link to comment
Share on other sites

45 minutes ago, The Moose In Your House said:

Can you do a suicide burn on any planet/moon?

I'm not sure I understand your question, but yes. Provided you have a sufficient TWR (which you need for landing anyway), you can do a suicide burn on any planet/moon.

On the ones with atmosphere, the trajectory will look a little different, but the priciple is still pretty much the same.

Link to comment
Share on other sites

  • 3 weeks later...

@Kobymaru

Got a new idea whilst writing in Lgg's topic: what if, instead of trying to guess a perfect retrograde burn, just estimate requirements for a "gamma-burn", i.e. first purely horizontal, after killing all h-speed purely vertical. It's less efficient than a perfect suicide burn, but we may call it a built-in safety margin.

Won't really work in thick atmospheres, but flies rockets sideways there anyway?

Link to comment
Share on other sites

For information to all who may be interested, I now have a fully developed and working math model to compute the ideal suicide burn on any body without an atmosphere. "Ideal" meaning it is the most efficient possible. While the routine simulating the suicide burn can solve with any starting condition, the model computes how to arrive at the perfect starting altitude for the suicide burn from any orbit around the body, and how to synchronize the vessel to arrive exactly at a designed landing location. All of the above using the least possible DV.

The model is NOT meant to work on bodies with an atmosphere, also because the ideal altitude for the suicide burn would be that low a vessel would almost certainly burn before getting there.

Link to comment
Share on other sites

1 minute ago, diomedea said:

For information to all who may be interested, I now have a fully developed and working math model to compute the ideal suicide burn

Very interested! Where can I find it?

Link to comment
Share on other sites

This thread is quite old. Please consider starting a new thread rather than reviving this one.

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.
Note: Your post will require moderator approval before it will be visible.

Guest
Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

×
×
  • Create New...