Escape velocity calculated vs measured

[surge]

I'm trying to calculate an escape from a moon with kOS. At the end of http://www.braeunig.us/space/orbmech.htm there is an equation for escape velocity. As I understand it, doing such a burn (minus the current orbital speed) should result in a parabolic orbit with eccentricity of exactly 1. When I try it though, it's way too much. Where rb is the radius at the burn point, vb is the speed at the burn point, and _eta is time to the burn point, I have:

    set vesc to sqrt(2*BODY:MU/rb).
    print "Total escape velocity required: " +round(vesc) +"m/s ".
    // lastly we subtract orbital velocity at that point to find dv
    set dv to vesc -vb.

    set n to node(TIME:SECONDS +_eta,0,0,dv).
    add n.

This gives me a manuevre node for which:

1) the ejection angle is not 90degrees from the burn point (which it should be, according to cos(n) = -1/e on the same page).

2) I can manually take off nearly 1/2 of the dv from the node and still get an ejection, (at which point the ejection point is 90degrees) so its clearly not a minimum.

Is there some kind of fudge factor in KSP that I don't know about? Ive seen the equation on other websites, so it's probably not wrong. I can't think of any other reason it's not working other than I'm misunderstanding something about the maths...

 

Edited July 2, 2016 by surge

[pincushionman]

What do you mean by "is way too much?" Is your measured e way off, or are you simply way past so-called escape?

Keep in mind, in KSP "escape" isn't a paranolic orbit where e=1, it's an orbit such that Ap is outside the (arbitrary) sphere of influence of the parent body, and as such is often still elliptical.

The game may not even have a special case for a parabolic escape orbit, since you have to go far out of your way to hit e=1 exactly. Rounding errors all but guarantee you're either in a really big ellipse or a really tight hyperbola. Both of which would be indistinguishable from a parabola over the short section rendered.

Edited July 2, 2016 by pincushionman

[PB666]

Its not true, escape velocity from the surface is theoretical. 

Case 1: no atmosphere, v = a * t, peak efficiency is when v is at the surface, rotation velocity at surface is pretty much trivial. since a is not inifinite you have gravity loses while accelerating. If it takes 5 minutes to raech MECO, then you are fighting and loosing energy due to 'hovering' over the lauchpad while trying to reach a velocity at that sltitude that allows escape. 

Case 2. Atmosphere, same acceleration constraints, but you ha a maximum dynamic pessure and at that you are losing x amount of energy to drag. 

So lets say Vescape is 11500 to escape the total energy you have is 1/2 11500^2. You can add to that energy loses due tongravity which is going to be something like gh for the period of time under acceleration, then the energy loses due to drag. You add that up multiply by two and take the square root and you get a delta v. 

 

If you are in orbit escape energy is the circukar orbital energy times two. So that if you are in a circular orbit, a single ideal burn represents a dV 0.414 x stable orbital velocity. So lets say you are orbiting the mun at 550 m/s and you want to leave the mun you need about 230 dV to leave muns orbit. Now lets say you want a deep orbit, you don't burn all of that and you burn a few dV to circularize to remain. Of course this is fantasy because you have la grange points and other troughs in the gravity wells you can slide through. So the n-body problem shifts this. The game it a bunch of serialized 2 body problems.

[Bill Phil]

https://en.wikipedia.org/wiki/Escape_velocity#Calculating_an_escape_velocity

You could use the circular orbit velocity calculation and just multiply by sqrt(2). So, your equation should be fine, I think, but KSP is not that realistic.

[surge]

Hmm, so I actually did the calculated burn just for the hell of it, and lo and behold, it did end up with eccentricity 0.999. Close enough to 1 for my purposes. Clearly, the problem is that I'm not understanding what they're talking about in the section about ejection angles and impact parameters at:

http://www.braeunig.us/space/orbmech.htm

What then, is cos(n) = -1/e meant to be representing, as I clearly have it wrong?

For e = 1, it's 180 degrees, which I thought meant the 'xxx escape' circle in the maneuvre is 90 degrees around the orbit from the burn point; essentially, it exits the moon going the opposite direction it would have came from if it was never in proper orbit. Why does it look more like 45 degrees? Is it some optical illusion? Who knows \o/

 

Edited July 2, 2016 by surge
Cant 'resolve' it, wrong thread. oops.

[PB666]

for any given circular orbital velocity like he siad escape velocity is 1.4143 x that velocity but the deltaV needed is that velocity minus the circular orbital velocity. therefore the dV needed is 0.4143 x orbital velocity in that circular orbit. 

[YNM]

Of course it's flawed - escape velocity as defined by sqrt(2) times circular orbital velocity assumes the end zero velocity to be located at infinity - we have finite SOIs in-game and even that's not really applicable IRL...

Another possible error is due to imperfection in continuity of the orbit. KSP, outside of warp and celestial bodies, slightly (hmm, not so sure the exact method) uses the euclidean integration, which over time results in errors in the velocity from perfect. It's also because of the slight phantom forces between attached parts. One way to tell this is to rotate a fairly one-side-heavy rocket around in space - you'll notice the Ap and Pe, as well as SMA, to change.

Yeah, it's not a perfect game. But the world isn't either - so just assume, like IRL, you gotta do the mid-course correction burns. (this part is why I refrain from using kOS. I tried but I'm never coder. And then I can't imagine how guys IRL manages to do it !)

Edited July 3, 2016 by YNM

[surge]

I see. Then does anyone know a way of calculating the (true|mean)anomaly at which to  burn so that I end up going a particular direction (relative to the parent body) at SOI crossing? Obviously I have the speed nailed, but the weird direction that it crosses the boundary is making the point at which to start the burn hard. Since I know the SOI boundary radius, it must be doable, but I'm waaay too dumb to figure out the equation for it. I suck at maths, and I suck even harder at polar geometry. In my head, I just turn all circles into squares  - 0, 90, 180, 270, 360

In other words, I'm trying to eject from the Muen and get into a specific orbit around Kerbin. Wildly gesticulating and yelling "THAT WAY!!!" just wont cut it.

[PB666]

4 hours ago, surge said:

I see. Then does anyone know a way of calculating the (true|mean)anomaly at which to  burn so that I end up going a particular direction (relative to the parent body) at SOI crossing? Obviously I have the speed nailed, but the weird direction that it crosses the boundary is making the point at which to start the burn hard. Since I know the SOI boundary radius, it must be doable, but I'm waaay too dumb to figure out the equation for it. I suck at maths, and I suck even harder at polar geometry. In my head, I just turn all circles into squares  - 0, 90, 180, 270, 360

In other words, I'm trying to eject from the Muen and get into a specific orbit around Kerbin. Wildly gesticulating and yelling "THAT WAY!!!" just wont cut it.

If you want to go prograde burn at the sunrise termination, angle to prograde for Kerbin is 180' ( in mechJeb) for most of the outer planets, but you are only accurate to the SOI, if you want to create a holmann transfer from LKO then you need to center your burn slightly after. The furher out in the system you want you apogee to be the later your burn. If you have some inclination corrections to make LKO is a good place to 'start'. If you really plan ahead you can launch with those inclinations and just expand upon them, in general however you will only be able to do a few degrees before the reward no longer pays off, this is because the sin aspect increases rapidly with only small changes in the cosine aspect, and also because there is a limit on how efficient the plane matching will be without an alignment of a transfer window start point time and an inclination node. 

For the inner planets or objects its a burn at angle to prograde of zero, in general, because the energy required to get to moho is so great and because its orbital period is so small you want to wait for a closer overlap of the inclination node, this then requires a inclination shift at launch and a considerable burn at sunset termination, AtP > 0, in this has you want to launch and not establish orbit but instead keep the verticle ascent positive until you have achieve your hohmann intercept preferably most below 75k alt. Once established you then correct inclination to optimal and make final inclination corrections in the transit. Moho is hard, not because of the time, but because of the dV required to transfer (you can use kerbins gravity well to improve) and moho has no atmosphere so you have to retro into moho at just about the lowest altitude posible, there is no way to soften this expense. So planning on the inclination side is the way to go, That 6 degree inclination change close to moho can be extremely costly. 

[surge]

PB666: launching from a surface is not relevant here, nor are hohmann transfers. If you're still having trouble:

Spoiler

 

parameter talt, tinc, poang.
// Target Altitude, Target Inclination, Pitchover Angle

clearscreen.
print "== Launch Ascent Guidance System (LAGS) ==".

set talt to talt *1000. // target orbital perigee
set safetyspeed to 50. // speed at which safety turn ends.
set safestart to 20. // Altitude at which safety turn starts.
set pitchhold to 10. // halt zero-lift turn at this pitch
set debugmode to true.

set critq to 15. // maxQ Not implemented yet.
// fudge factor - how 'pointy' rocket is
set dragf to 0.4. // Not implemented yet

// ===END CONFIG OPTIONS===

// release controls early in case we get interrupted.
SET SHIP:CONTROL:PILOTMAINTHROTTLE TO 0.

lock LAc to 90 -vang(up:vector, ship:velocity:surface).
lock dynp to SHIP:DYNAMICPRESSURE *constant:AtmToKPa.

function printscreen {
// were going to use a 15 line screen total.
// Debug messages start on line 20
print " " at(0,1).
print "Target altitude: " +round(talt/1000, 1) +"km   " at(0,2).
print "Target inclination: " +round(tinc, 1) +"o   " at(25,2).
// Vector start, Vector end
// Launch mass, twr
print "Apogee: " +round(SHIP:OBT:APOAPSIS/1000,1) +"km     "at(0,6).
print "Perigee: " +round(SHIP:OBT:PERIAPSIS/1000,1) +"km    "at(25,6).
print "Time to apogee: " +round(ETA:APOAPSIS) +"s   "at(0,7).
print "Ecc: " +round(SHIP:OBT:ECCENTRICITY,2) +"   " at(0,8).
print "Inc: " +round(SHIP:OBT:INCLINATION,2)  +"   " at(25,8).
print "Dynamic pressure: " +round(dynp, 2) +"/" +round(critq, 2)
    +"kPa     " at(0,9).
print "Flight path angle: " +round(LAc, 1) +"   " at(0,10).
print " " at(0,11).
print " " at(0,12).
print " " at(0,13).
print " " at(0,14).
}

wait 0.5.
printscreen().
print "STATUS: calculating...                    " at(0,20).

function clamped {
// returns stage number of clamps any
    local stagecmp to STAGE:number.
    local highclamp to -1.
    
    if SHIP:STATUS <> "PRELAUNCH" and SHIP:STATUS <> "LANDED"
        return -1.

    local clamps to ship:partsnamed("launchClamp1").
    for c in clamps {
        if (c:stage > highclamp) set highclamp to c:stage.
    }
    return highclamp.
}

set coremotors to list().
set boostmotors to list().
function detect_motors {
// parse and sort motors in current stage into
// globals coremotors and boostmotors
    
    local el to list().
    list engines in el.
    
    set coremotors to list().
    set boostmotors to list().
    
    local stagehigh to STAGE:number.
    local stagelow to clamped().
    if (SHIP:STATUS = "PRELAUNCH") and (stagelow = -1) {
        set stagelow to stagehigh -1.
    } else if (stagelow = -1)
        set stagelow to stagehigh.

    for e in el {
        // we have to detect if stage above has launch
        // clamps AND if we are in launch ready state
        if (e:stage >= stagelow) and (e:stage <= stagehigh)
        {
            // Check if the engine is
            // on centre stack.
            if (e:tag = "boostMotor") {

                boostmotors:add(e).
            } else {
                // everything else should be cores.
                coremotors:add(e).
            }
        }
    }
    return coremotors:LENGTH +boostmotors:LENGTH.
}

function setcorethrottle {
    parameter lvl.
    for e in coremotors {
        set e:THRUSTLIMIT to lvl *100.
    }
}

function getcorethrottle {
    if (coremotors:length = 0) return 0.
    return coremotors[0]:THRUSTLIMIT /100.
}

function getctwr {
    // returns current TWR.
    // *sigh* :MAXTHRUST doesnt work unless the engine is
    // staged. This SHOULD return the TWR upon launch, so
    // we only have to calculate it once.
    local gf to (SHIP:BODY:MU
        /(SHIP:BODY:RADIUS +SHIP:ALTITUDE)^2).
    local thrust to 0.

    list engines in el.
    for e in el
        if e:IGNITION set thrust to thrust +e:THRUST.

//    return SHIP:AVAILABLETHRUST /(SHIP:MASS*gf).
    return thrust /(SHIP:MASS*gf).
    
}

function inc2hdg {
// does the wierd maths when converting from an orbital inclination
// to a heading to fly for it.
    parameter inc.
    set rt to 180 -inc -90.
    if (rt < 0) set rt to 180 -rt.
    return rt.
}
set thdg to inc2hdg(tinc).

set east to up *R(0, -90, 0).
function getvhdg {
    local r1 to vcrs(up:vector,
        ship:velocity:orbit):normalized.
    local vhdg to vang(r1, east:vector).
    if (r1:y) > 0.0 set vhdg to 360 -vhdg.
    return vhdg.
}

set laststage to STAGE:NUMBER.
function stagectl {
    if not STAGE:READY return false.
    if THROTTLE <= 0 return false.

    // detect_motors() takes long and only needs to be
    // called once after staging has finished
    if (laststage <> STAGE:NUMBER) {
        detect_motors().
        set laststage to STAGE:NUMBER.
        return false.
    }
    
    // also check for no motors
    if (coremotors:length = 0) and (boostmotors:length = 0){
        lock throttle to 1.
        stage.
        return true.
    }

    for e in boostmotors {
        if (e:maxthrust <= 0) {
            lock throttle to 1.
            stage.
            return true.
        }
    }

    if (coremotors:length > 0) {
        // check that boost motors are gone.
        // but still need to do the fuel check thing.
        if (boostmotors:length = 0) setcorethrottle(1).

        local corethrust to 0.
        for e in coremotors {
            set corethrust to corethrust +e:maxthrust.
        }
        if (corethrust = 0) {
            stage.
            return true.
        }
    }
    return false.
}

function hdgerr {
    // Calculates the heading offset for the given orbital
    // inclination.
    parameter targethdg.
    
    // ffs. SHIP:OBT:INCLINATION doesnt do negatives?
    set shipinc to SHIP:OBT:INCLINATION.
    local shiphdg to getvhdg().
    if (shiphdg > 90 and shiphdg < 270) set shipinc to -shipinc.
    return targethdg -inc2hdg(shipinc).
}

detect_motors().
if boostmotors:LENGTH = 0 {
    print "WARNING NO BOOST MOTORS DETECTED!" at (0,14).
//    exit.
}

//
// ===GUIDANCE CALCULATION===
//

set poalt to SHIP:BODY:ATM:HEIGHT/7.

// FIXME: This needs to settle faster.
SET thpid to PIDLOOP(0.1, 0.8, 0.15, 0, 1).
set thpid:setpoint to 0.0.

SET pchpid to PIDLOOP(0.5, 0.05, 0.5, 0, 20).
set pchpid:setpoint to 15.

//set STEERINGMANAGER:PITCHPID:KD to 0.5.
//set STEERINGMANAGER:YAWPID:KD to 0.5.

FROM {local countdown is 5.}
UNTIL countdown = 0 STEP {SET countdown to countdown - 1.} DO {
    PRINT "STATUS: pre-launch... T-"
    + countdown +"s            " at(0,20).
    WAIT 1.
}

print "STATUS: Launch...                    " at(0,20).
lock throttle to 1.
lock steering to up.
// Since engine:MAXTHRUST is broken, we have to actually fire
// up the engines (possibly launch if no clamps) to get a value.
set launchtwr to 0.
set sa to SHIP:ALTITUDE -SHIP:GEOPOSITION:TERRAINHEIGHT.
until SHIP:ALTITUDE -SHIP:GEOPOSITION:TERRAINHEIGHT-sa > safestart {

    set launchtwr to getctwr().

    local dTWR to getctwr() -1.3.
    setcorethrottle(thpid:update(TIME:SECONDS*10, dTWR)).
    
    print "Launch TWR: " +round(launchtwr, 2) +"     " at(25,3).
    print "Launch Mass: " +round(SHIP:MASS, 2) +"    " at(0,3).
    stagectl().
    
    printscreen().
    wait 0.1.
}
unset sa.
//set poang to 90-(-100*(1/dragf*0.2)*(launchtwr -1.8)).
print "Zero-lift start: " +round(poalt) +"m   " at(0,4).
print "Vectoring angle: " +round(poang, 1) +"o    "at(25,4).
unset startalt.

print "STATUS: safety turn...                    " at(0,20).
lock steering to heading(thdg, 85).
set rcs to true.
until (SHIP:VELOCITY:SURFACE:MAG >= safetyspeed) {
    stagectl().

    printscreen().
    wait 0.5.
}

print "STATUS: vector guidance...                    " at(0,20).
lock pitcherr to SHIP:ALTITUDE *(90-poang)/poalt.
lock steering to heading(thdg, 90 -pitcherr).
until (SHIP:ALTITUDE > poalt)
{
    stagectl().

    printscreen().
    wait 0.5.
}

print "STATUS: zero-lift turn...                    " at(0,20).
lock steering to heading(thdg, LAc).
// FIXME: do throttle control
lock steering to heading(thdg +hdgerr(thdg), LAc).
until (LAc <= pitchhold) or (SHIP:APOAPSIS >= talt) {
    stagectl().

    printscreen().
    wait 0.5.
}

print "STATUS: pitch hold...                    " at(0,20).
lock steering to heading(thdg +hdgerr(thdg), pitchhold).
until (SHIP:ALTITUDE >= SHIP:BODY:ATM:HEIGHT) {
    stagectl().
    lock throttle to thpid:update(TIME:SECONDS,
        SHIP:APOAPSIS -talt).

    local pitchang to pchpid:update(TIME:SECONDS*10, LAc).
    lock steering to heading(thdg +hdgerr(thdg), pitchang).

    printscreen().
    wait 0.5.
}
lock throttle to 0.
unlock steering.
unlock LAc.
unlock dynp.
RCS off.

print "STATUS: Orbital Insertion...                 " at(0,20).
set apR to APOAPSIS +BODY:RADIUS.
set dV to sqrt(BODY:MU /apR)
    -sqrt(BODY:MU* (2/apR -1/OBT:SEMIMAJORAXIS)).
set circNode to node(TIME:SECONDS +ETA:APOAPSIS, 0, 0, dV).
add circNode.

// we could check the node is correct by comparing
// circNode:ORBIT:APOAPSIS and circNode:ORBIT:PERIAPSIS
// to make sure they are a) nearly equal, and b) close to talt

run xnode(1, 0).
print "STATUS: Program ended.                " at(0,20).

 

 

[surge]

For the googlers: since 'escaped' is mathematically just an elliptical orbit, to calculate KSP escape velocity is then:

v = sqrt(BODY:MU *(2/r -1/(r +BODY:SOIRADIUS)/2)) +0.01.

where r is BODY:RADIUS +SHIP:ALTITUDE (or whatever altitude you want to launch from)

The 0.01 is not mathematically correct, but just to 'nudge' it over the line.

Still working on figuring out the transit angle so I can figure out *when* to burn, but this is a start.

Edited July 21, 2016 by surge