xeger / kos-ramp

Relatively Adequate Mission Planner: a rudimentary, scriptable autopilot for kOS.
Other
55 stars 23 forks source link

Allow node_apo() to place node at a time other than periapsis #30

Closed mathuin closed 6 years ago

mathuin commented 6 years ago

The comments below apply to both node_apo() and node_peri(), obviously.

I recently wanted to change my apoapsis at a point in my orbit that was not my periapsis.

My script lifts heavily from node_apo.ks and looks like this:

local mu is body:mu.
local br is body:radius.
local vom is ship:obt:velocity:orbit:mag.
// local r is ship:body:altitudeof(positionat(ship, time:seconds)) + br.
local r is br + altitude.
local v1 is velocityat(ship, nodetime):orbit:mag.

// ugh v2 is hard.
local sma1 is obt:semimajoraxis.
local r2 is ship:body:altitudeof(positionat(ship, nodetime)) + br.

local sma2 is (tgtapo + br + r2)/2.
local v2 is sqrt(vom^2 + (mu * (2/r2 - 2/r + 1/sma1 - 1/sma2))).

local deltav is v2 - v1.
local nd is node(nodetime, 0, 0, deltav).
add nd.

I know that this code "works" in the sense that the node created will place the apoapsis at the right distance at the right time. Would it be worth merging this code into node_apo() and creating a second optional argument for node time (with default value TIME:SECONDS+ETA:PERIAPSIS)? I'm willing to submit a PR if it has a chance of being accepted.

xeger commented 6 years ago

Seems reasonable to me. Out of curiosity, what was your need to raise apoapsis at a specific point?

mathuin commented 6 years ago

I am automating the construction of a relay network. The network design includes two highly elliptical polar orbits, one over each pole, offset such that one satellite is far while the other is near and vice versa. These orbits require the apoapsis burn to occur as close to the opposite pole as possible. Here's a drawing of the network, taken from the KSP forum thread which inspired me:

29389037634_6d304685ac_o

My code to calculate the appropriate burn time is rather naive:

function timeToNode {
    parameter tgtlat. // could be 90 or -90!
    set step to 1. // 0.1 is overkill
    set epsilon to 0.01. // probably never reached
    set toobig to time:seconds+ship:orbit:period.
    set nodetime to time:seconds.
    // closest has to approach 0 when time is best.
    // lock closest to abs(vectorangle(positionat(ship, nodetime), -ship:body:angularvel)-90.0).
    lock closest to abs(ship:body:geopositionof(positionat(ship, nodetime)):lat-tgtlat).
    set bestsofar to abs(tgtlat).
    set besttime to toobig+step.
    until closest < epsilon or nodetime > toobig {
        if closest < bestsofar {
            set bestsofar to closest.
            set besttime to nodetime.
        }
        set nodetime to nodetime + step.
    }
    if nodetime > toobig {
        // print "using least bad value (closest " + bestsofar + ")".
        return besttime.
    }
    return nodetime.
}

This same algorithm will be handy when I want to place the periapsis of a resonant orbit over somewhere specific, say KSC, in order to deploy a network of synchronous satellites with one directly over the center. It's just so ... inelegant. There has got to be a better way to do this. I hope to read enough of the Orbital Mechanics text to find a better way to calculate this value. :laughing:

Anyway, with this change, my script will call run node_apo(body:soiradius*0.95, nodetime). instead of the code in the first comment.

mathuin commented 6 years ago

Okay, now I feel dumb for not seeing node_alt.ks, which does most of the heavy lifting that my script needs. I'm a little surprised that node_apo() and node_peri() don't (yet) call node_alt() as they're both special cases of that script once the 120 second limit becomes a parameter.

ghost commented 6 years ago

@mathuin Look at utilDtTrue in lib_util I have added recently, just give it angle from periapsis and it will return ETA (unless you are on escape trajectory, TBD, works for any eccentricity below 1).

And you may also look at burnTimeForDv in lib_staging in my latest pull request. It currently cannot take staging into account, but may in the future. It would be good to have all scripts use it that we can implement the logic in one place.

mathuin commented 6 years ago

I am super-excited about your lib_staging changes. Once they're stabilized, I'd like to adjust the node times for apoapsis and periapsis to account for the actual burn times.

I will try to work your utilDtTrue function into my timing code but I'm not exactly sure how. :-)

ghost commented 6 years ago

I am not forcing you to use utilDtTrue, but it really is the ultimate ETA function. utilDtTrue(0) = eta:periapsis utilDtTrue(180) = eta:apoapsis 180-orbit:argumentOfPeriapsis and 360-orbit:argumentOfPeriapsis are eq. AN/DN nodes. you may look into node_inc how I use it, orbit:trueAnomaly is your current angle from periapsis. https://en.wikipedia.org/wiki/Argument_of_periapsis

but if orbit:eccentricity = 0 then it is same as utilDtMean (true and mean become the same, no conversion needed) which is really simple, you move with same speed everywhere: https://en.wikipedia.org/wiki/Orbital_speed v = sqrt(body:mu/r) which you can see in circ.ks in my latest PR.

Nothing against iterative solution, but I would think about something like binary search - halving intervals, start with orbit:period, see if positionAt(ship,orbit:period/4) is closer than positionAt(ship,orbit:period*3/4) and you should know in which half you need to be. Repeat until you reach interval span of second or less. Just an idea. This could potentially shrink a loop of 1000 iterations to just 10 iterations, 1.000.000 iterations to just 20 (log2, 2^10 = 1024, 2^20 > 1.000.000).

mathuin commented 6 years ago

My first draft of the algorithm was binary search, but I had trouble getting it to settle down due to the way that the distance from orbit changes -- I wasn't handling the inflection change correctly. This worked, and didn't take too much time. :-)

Since I'm looking for polar orbits, I could definitely optimize by searching between AN and DN for North Pole and DN and AN for South Pole, respectively. That's clearly helpful. What I don't yet see is how to use the function as a magic bullet for all problems, but that will come with time and experience.

mathuin commented 6 years ago

Oh hey, the code to fix this has been committed. :-)