surge

Advanced Orbital Mechanics using kOS

Recommended Posts

After a lot of research on various orbital manuevres, I've managed to create a fairly simple bit of code that allows one to set both perigee and apogee in one burn.

Most of the research I did showed that some people think this is impossible without complicated loops, and it's barely mentioned, if at all, in the textbooks I have access to, so I thought I'd share it here.

Obviously, there are some conditions:

1) The requested orbit and the current orbit must cross. It will crash if not - it would be mathematically and physically impossible to do in 1 burn.

2) If the burn is not done at an apse, it will cause the apse line to rotate.

Obviously there are more fuel efficient ways to do this, however it can be useful for minor correction burns or when time is of the essence (perhaps you're about to crash into a planet!). Probably the most important advantage is that it can be done anywhere in the orbit (depending on how much fuel you can spare), so no more endlessly circularising and faffing about with those stupid "hohmann" burns. Yay!

Anyway, on with the show:

Spoiler

 

function oneorbit
{
// Sets up a single burn in time t to effect an orbit with
// perigee p and apogee a.
// Necessarily will rotate the apse line, unless done at an apsis.
    parameter t, p, a.

    print "oneorbit(): t = " +round(t)
        +"s, p = " +round((p-BODY:RADIUS)/1000,2)
        +"km, a = " +round((a-BODY:RADIUS)/1000, 2) +"km  ".

    // initial orbit characteristics at burn point
    local Vr to positionat(SHIP, TIME +t) -BODY:POSITION.
    local Vi to velocityat(SHIP, TIME +t):ORBIT.
    local Vn to VCRS(Vr, Vi):normalized. // normal
    local Vt to VCRS(Vn, Vr):normalized. // tangent direction
    local frame to -SHIP:PROGRADE
        *rotatefromto(Vi, OBT:VELOCITY:ORBIT).

    // desired final orbit characteristics
    local smaf to (p +a) /2.
    local Sf to sqrt(BODY:MU *(2/Vr:mag -1/smaf)).
    local Sp to sqrt((2*BODY:MU *a)/(p*(a+p))).

    // calculate flight path angles...
    local zaf to arcsin((Sp *p)/(Vr:mag *Sf)).
    local fpaf to 90 -zaf.

    local Vf to (angleaxis(fpaf, Vn) *Vt):normalized *Sf.
    local dVf to frame *(Vf -Vi).

    local nd to node(TIME:SECONDS +t, dVf:y, dVf:x, dVf:z).
    add nd.
    
    return dVf:mag.
}

 

 

 

Share this post


Link to post
Share on other sites

Here's a simplified version that just circularises.

 

Spoiler

 

function circ {
    parameter lead.

    set Vr to positionat(SHIP, TIME +lead) -BODY:POSITION.
    set Vb to velocityat(SHIP, TIME +lead):ORBIT.
    // set a rotation frame to fix up wacky KSP coordinates...
    set frame to -SHIP:PROGRADE
        *rotatefromto(Vb, OBT:VELOCITY:ORBIT).

    set Sc to sqrt(BODY:MU /Vr:mag). // circular speed
    set onorm to VCRS(Vb, Vr).
    set Vc to VCRS(Vr, onorm):normalized *Sc.
    set dV to frame *(Vc -Vb).
    if dV:mag > 0.1 {
        print "circ(): circularizing in "
            +round(lead,1) +"s   ".
        set cn to node(TIME:SECONDS +lead, dV:y, dV:x, dV:z).
        add cn.
        return cn.
    }
    return False.
}


 

 

Share this post


Link to post
Share on other sites

Are we there yet, dad?

These will tell you how long it takes to get to another point in your orbit.

Spoiler

function eccan
{
// returns eccentric anomaly for true anomaly t
    parameter t.
    local ecc to OBT:ECCENTRICITY.
    local i1 to sqrt((1 -ecc) /(1 +ecc)).
    return arctan(i1 *tan(t/2)) *2.
}

function tofang
{
// calculates time to rotate angle a around orbit from now.
    parameter a.
    local ecc to OBT:ECCENTRICITY.

    local Eo to eccan(OBT:TRUEANOMALY).
    local Et to eccan(OBT:TRUEANOMALY +a).

    // not sure why we need *RadtoDeg, ecc is a ratio
    // and everything else is in degrees. it works though.
    local r2d to CONSTANT:RadtoDeg.
    local Mo to Eo -ecc *sin(Eo) *r2d.
    local Mt to Et -ecc *sin(Et) *r2d.

    // Account for angles > 360. Dont do backward angles if
    // we cross perigee
    set dM to mod(Mt -Mo +360, 360).
    return (dM *(OBT:PERIOD/360)
        +(OBT:PERIOD *floor(a/360)))
        *(a/abs(a)).
}

It could be modified to determine time of flight between random points in the orbit. I'll leave that as an exercise, or maybe present it later if the personal need arises.

Tip: The hard part about that, is that you can't get true anomaly from the other angles/"anomalys" easily, so you have to calculate the first point another way.

Edited by surge

Share this post


Link to post
Share on other sites
On 6.9.2016 at 3:16 AM, surge said:

After a lot of research on various orbital manuevres, I've managed to create a fairly simple bit of code that allows one to set both perigee and apogee in one burn.

Most of the research I did showed that some people think this is impossible without complicated loops, and it's barely mentioned, if at all, in the textbooks I have access to, so I thought I'd share it here.

Obviously, there are some conditions:

1) The requested orbit and the current orbit must cross. It will crash if not - it would be mathematically and physically impossible to do in 1 burn.

2) If the burn is not done at an apse, it will cause the apse line to rotate.

Obviously there are more fuel efficient ways to do this, however it can be useful for minor correction burns or when time is of the essence (perhaps you're about to crash into a planet!). Probably the most important advantage is that it can be done anywhere in the orbit (depending on how much fuel you can spare), so no more endlessly circularising and faffing about with those stupid "hohmann" burns. Yay!

Anyway, on with the show:

  Hide contents

 

function oneorbit
{
// Sets up a single burn in time t to effect an orbit with
// perigee p and apogee a.
// Necessarily will rotate the apse line, unless done at an apsis.
    parameter t, p, a.

    print "oneorbit(): t = " +round(t)
        +"s, p = " +round((p-BODY:RADIUS)/1000,2)
        +"km, a = " +round((a-BODY:RADIUS)/1000, 2) +"km  ".

    // initial orbit characteristics at burn point
    local Vr to positionat(SHIP, TIME +t) -BODY:POSITION.
    local Vi to velocityat(SHIP, TIME +t):ORBIT.
    local Vn to VCRS(Vr, Vi):normalized. // normal
    local Vt to VCRS(Vn, Vr):normalized. // tangent direction
    local frame to -SHIP:PROGRADE
        *rotatefromto(Vi, OBT:VELOCITY:ORBIT).

    // desired final orbit characteristics
    local smaf to (p +a) /2.
    local Sf to sqrt(BODY:MU *(2/Vr:mag -1/smaf)).
    local Sp to sqrt((2*BODY:MU *a)/(p*(a+p))).

    // calculate flight path angles...
    local zaf to arcsin((Sp *p)/(Vr:mag *Sf)).
    local fpaf to 90 -zaf.

    local Vf to (angleaxis(fpaf, Vn) *Vt):normalized *Sf.
    local dVf to frame *(Vf -Vi).

    local nd to node(TIME:SECONDS +t, dVf:y, dVf:x, dVf:z).
    add nd.
    
    return dVf:mag.
}

 

 

 

Nice! What do you mean by "rotating apse line"? I'm not familiar with that term. Longitude of the apse changes?

Edited by kurja

Share this post


Link to post
Share on other sites

It's fairly self-explanatory. According to textbooks, the 'apse line' is the line between periapsis and apoapsis. If you don't understand what's happening, I guess this maneuvre isn't for you!

Edited by surge

Share this post


Link to post
Share on other sites

I am curious as to your source for the conversions of the Keplerian elements, since they are a bit different (and seem much more concise) than the ones I have found or have derived using my own rusty algebra, for example your eccentric anomaly from true anomaly formula. If you have a reference, I would love to know which one you use.

Thanks in advance.

 

Share this post


Link to post
Share on other sites

The word 'Inflammable' was officially removed from real dictionaries at about the same time Pluto was declassified as a planet. Both of them had been in use for, lets average it out to about 1000 years.

Narcissistic rocket scientists keep making up new names for the same concept just because they can.

Does that help?

If it doesn't, I chose the least mathematically computative of them. I'm a computer scientist, not a rocket scientist, and I have very little idea or interest in what you're talking about.

If you can improve it computationally, go ahead, this is copyright free software!

Edited by surge

Share this post


Link to post
Share on other sites
Quote

After a lot of research on various orbital manuevres, I've managed to create a fairly simple bit of code...

This is what I was hoping you could help with, based on your research. If you found a good reference, I would be grateful if you could share it, since, as I said, your could seems better than what I have been able to find or derive on my own. I am not sure why you say you don't know what I mean about True Anomaly or Eccentric Anomaly, since you use those terms in your code.

Since I am not a rocket scientist nor a computer scientist, I doubt I can improve on it, which is why I asked the question. 

Share this post


Link to post
Share on other sites

Well, ok. I've read parts of books such as "Orbital Mechanics for Engineering Students", "Fundamentals of Astrodynamics", Bate Mueller, and "Fundamentals of Astrodynamics and Applications", D. Vallado.

But actually there is a video on youtube that explains all the "anomalies" (there is of course Mean Anomoly too), and that was in part, my real source. The guy that made such an eloquent video is called OrbitNerd

And of course braeunig's site which had all the fundamental mathematical equations in non-ridiculous terms that all those books above didn't. I essentially just condensed his work into the above code.

Happy reading.

 

Share this post


Link to post
Share on other sites

Reiterating over Lambert problems is taking forever? Tired of the sysadmins yelling at you for clogging up the cpus? Fear not, the laziest solution is relatively simple!

Spoiler

function lambertMinE {
//
// Calculate V1 at R1 to get to R2.
//
    parameter R1, R2.

    local R1R2 to R1:mag *R2:mag.
    local cosdA to VDOT(R1, R2) /R1R2.
    local sindA to VCRS(R1, R2):mag /R1R2.

    // some weird excrements to do with geometry.
    local c to sqrt(R1:mag^2 +R2:mag^2 -2 *R1R2 *cosdA).
    local s to (R1:mag +R2:mag +c)/2.

    // orbital parameters
    local amin to (R1:mag +R2:mag +c) /4.
    local pmin to R1R2 /c *(1 -cosdA).
    local emin to sqrt(1 -2 *pmin /s).

    local k to sqrt(BODY:MU *pmin) /(R1R2 *sindA).
    local V1 to k *(R2 -(1 -R2:mag/pmin * (1 -cosdA)) *R1).
    return V1.
    
    // FIXME: we could give TOF & V2 too, i think? (untested)
    local b to arcsin(sqrt((s -c) /s)) *2.
    local dT to sqrt(amin^3/BODY:MU) *(CONSTANT:PI -b +sin(b)).
//    return list(V1, V2, dT).
}

It will give you V1, which is how fast you need to be going (as in the dV and direction of your burn) at a radius of R1 to hit something at a radius of R2.

As you can see, you only need to supply the two "positions". It can potentially do the rest; V2 is the vector you will to be going when you get there, and dT is an *untested* calculation of the time it will take.

Useful if you want to land at an exact point on a planet, or rendevous with something else out there.

As always there are caveats:

1) It doesn't take into account movement/rotation of targets or planets; you must pre-calculate that.

2) It calculates "minimum energy", but it's minimum energy from an absolute 0 starting velocity. So if you try to use this to descend from a higher orbit, it will often give a very inefficient path that "jumps up" first, then descends to the target radius. It really works best from a lower orbit, or even the ground, like an ICBM.

Now you too, can act like america and threaten everyone on Kerbin with this new technology!

 

Edited by surge
i wonder why the censorship pluralised 's_h_i_t', singular

Share this post


Link to post
Share on other sites

If you're left wondering how to turn these vectors into something tangable, this function converts body-relative vectors into the strange ship-relative triplets required by node():

FUNCTION nodeFromVector
// from a reddit post somewhere (I've forgotten)
{
  PARAMETER vec, n_time IS TIME:SECONDS.
  LOCAL s_pro IS VELOCITYAT(SHIP,n_time):ORBIT.
  // The following assumes you do not change sphere of influence
  // between now and n_time
  LOCAL s_pos IS POSITIONAT(SHIP,n_time) - BODY:POSITION.
  LOCAL s_nrm IS VCRS(s_pro,s_pos).
  LOCAL s_rad IS VCRS(s_nrm,s_pro).

  LOCAL pro IS VDOT(vec,s_pro:NORMALIZED).
  LOCAL nrm IS VDOT(vec,s_nrm:NORMALIZED).
  LOCAL rad IS VDOT(vec,s_rad:NORMALIZED).

  RETURN NODE(n_time, rad, nrm, pro).
}

As noted I don't take credit for this (nor do I really take credit for anything here, since it's all just regurgitated equations from other sources). I found it on reddit after the much worse version I invented didn't work properly (see the oneorbit() post). Anyway, you can now do things like:

set V1 to lambertMinE(R1,R2) -velocityat(SHIP, burntime).
set n to nodeFromVector(V1, TIME:seconds +burntime).
add n.

to add the calculated maneuvres into the map. Making them actually happen is a less advanced topic not for this thread.

Edited by surge

Share this post


Link to post
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.