Kerbal Space Program satellite orbit phasing using kOS and MechJeb

phase.ks

// Change the phase of your orbit, that is, advance yourself a certain
// amount of time along the same orbit path.
//
// See https://en.wikipedia.org/wiki/Orbit_phasing for the math and lots
// of explanations.
//
// AUTHOR: Eric Kidd.
//
// REQUIREMENTS: To run this, you'll need Kerbel Space Program, the kOS mod,
// and a computer on your craft.
//
// HOW IT WORKS: When run, this will set up two maneuver nodes, but not actually
// pilot the ship. If you want an autopilot, considering installing MechJeb
// and using it to execute the planned maneuver nodes.
//
// ASSUMPTIONS: This script assumes that your orbit already has the right period
// and inclination, but it has the wrong phase. It's mostly useful for spacing
// out communications relay satellites evenly.
//
// This script uses a single phasing orbit, starting at apoapsis and going
// _down_, because I use it for setting up relay satellites in high orbits
// around moons after releasing them from my "satellite bus". So I'm more
// concerned by completely exiting the moon's weak sphere of influence than
// I am by crashing into it. But if you're setting up satellites in very low
// orbits, then you'll want to edit the code to burn at periapsis.
//
// CODING STYLE: Whenever possible, all functions and variables are marked with
// the metric units of the values they contain. This helped eliminate a lot
// of surprises.

// PARAMETERS

// Set this to the number of degrees you want to phase forward. Large
// positive numbers may cause you to crash into the planet if you're low.
// Negative numbers may cause you to escape if you're near escape velocity.
set phase_forward_deg to 150.

// IMPLEMENTATION

print "Calculating nodes to phase orbit " + phase_forward_deg + " degrees forward.".
add_phasing_nodes(phase_forward_deg).

// Calculate the period of a phasing orbit that allows us to make the entire
// rephase in one orbit. This may not be optimal, or it may launch you
// into the outer darkness, or it may set you up for surprise lithobraking.
// Please carefully examine all maneuver nodes before executing. I mean, you
// don't _really_ want to send a bunch of new satellites to the
function phasing_orbit_period_s {
    parameter vessel.
    parameter phase_forward_deg.

    local original_orbit_period_s is vessel:orbit:period.
    local time_shift_s is (phase_forward_deg/360) * original_orbit_period_s.
    return original_orbit_period_s - time_shift_s.
}

// Calculate the "standard gravitational parameter" for a celestial body.
// See https://en.wikipedia.org/wiki/Standard_gravitational_parameter
function standard_gravitational_parameter_m3_per_s2 {
    parameter body.

    local body_mass_kg is body:mass.
    print "Celestial body mass is " + body_mass_kg + "kg".
    local gravitational_constant_m3_per_kg_s2 is constant:G.
    return body_mass_kg * gravitational_constant_m3_per_kg_s2.
}

// Calculate the semimajor axis of a circular orbit with the specified period
// around a cellestial body.
function circular_orbit_semimajor_axis_m {
    parameter orbit_period_s.
    parameter body.

    local mu_m3_per_s2 is standard_gravitational_parameter_m3_per_s2(body).
    local m1_5 is sqrt(mu_m3_per_s2) * orbit_period_s.
    return (m1_5 / (2*constant:PI)) ^ (2/3).
}

// Get the _other_ apsis of an orbit.
function other_apsis_m {
    parameter apsis_m.
    parameter circular_semimajor_axis_m.
    return 2*circular_semimajor_axis_m - apsis_m.
}

// Vessel mass is supposedly in metric tons, not kg. But body masses are
// supposedly in kg. I'm filled with fear and doubt, and fully expect to wind
// up lithobraking into Mun again.
//
// However, after messing with this code, it turns out I get the right answers
// if I assume vessel masses are in kg.
function vessel_kg {
    parameter vessel.
    return vessel:mass.
}

// Get the actual apoapsis, including the sea level of the body.
function real_apoapsis_m {
    parameter orbit.
    return orbit:apoapsis + orbit:body:radius.
}

// Get the actual periapsis, including the sea level of the body.
function real_periapsis_m {
    parameter orbit.
    return orbit:periapsis + orbit:body:radius.
}

/// Get the angular momentum of a vessel's orbit.
function orbit_angular_momentum_kg_m2_per_s {
    parameter vessel.

    local apoapsis_m is real_apoapsis_m(vessel:orbit).
    local periapsis_m is real_periapsis_m(vessel:orbit).
    return orbit_angular_momentum_for_apsides_kg_m2_per_s(periapsis_m, apoapsis_m, vessel).
}

// Calculate the angular momentum of an orbit, given its apoapsis and periapsis.
// Normally, angular momentum is kg*m^2/s, but the formula on Wikipedia gave
// m^s/s. I think they just cancelled it out in a later part of the equations.
function orbit_angular_momentum_for_apsides_kg_m2_per_s {
    parameter apoapsis_m.
    parameter periapsis_m.
    parameter vessel.

    local body is vessel:orbit:body.
    local mu_m3_per_s2 is standard_gravitational_parameter_m3_per_s2(body).
    print "Standard gravitational parameter is " + mu_m3_per_s2 + "m^3/s^2".
    local m2_per_s is sqrt(2*mu_m3_per_s2 * apoapsis_m * periapsis_m / (apoapsis_m + periapsis_m)).
    return vessel_kg(vessel) * m2_per_s.
}

// Get the delta V for one our two burns (which will be identical) in m/s.
function phasing_delta_v_m_per_s {
    parameter vessel.
    parameter phase_forward_deg.

    // We need to calculate how much angular momentum we'll want in our
    // phasing orbit.
    local new_period_s is phasing_orbit_period_s(vessel, phase_forward_deg).
    print "Changing period from " + vessel:orbit:period + "s to " + new_period_s.
    local new_circ_semimajor_m is circular_orbit_semimajor_axis_m(new_period_s, vessel:orbit:body).
    print "(Circular semimajor axis " + new_circ_semimajor_m + "m)".
    local both_apoapsis_m is real_apoapsis_m(vessel:orbit).
    local phasing_peripasis_m is other_apsis_m(both_apoapsis_m, new_circ_semimajor_m).
    print "Phasing periapsis will be " + phasing_peripasis_m + "m".
    local new_angular_momentum_kg_m2_per_s is
        orbit_angular_momentum_for_apsides_kg_m2_per_s(both_apoapsis_m, phasing_peripasis_m, vessel).
    print "Phasing orbit speed " + new_angular_momentum_kg_m2_per_s / (vessel_kg(vessel) * both_apoapsis_m) + "m/s".

    // And we can calculate how much angular moment we have now.
    local current_angular_momentum_kg_m2_per_s is orbit_angular_momentum_kg_m2_per_s(vessel).
    print "Current orbit speed " + current_angular_momentum_kg_m2_per_s / (vessel_kg(vessel) * both_apoapsis_m) + "m/s".

    // Now we need to convert our change in momentum into a change in velocity.
    local delta_angular_momentum_kg_m2_per_s is
        new_angular_momentum_kg_m2_per_s - current_angular_momentum_kg_m2_per_s.
    return delta_angular_momentum_kg_m2_per_s / (vessel_kg(vessel) * both_apoapsis_m).
}

// Add two nodes which will change our orbit phase. This only operates on the
// current ship, which we need to be officially piloting according to the
// tracking center, because that's the only one with maneuver nodes.
function add_phasing_nodes {
    parameter phase_forward_deg.

    // Remove all existing maneuver nodes.
    for node in allnodes {
        remove(node).
    }

    // The first burn shifts us into a faster phasing orbit.
    local start_first_burn_s is time:seconds + eta:apoapsis.
    local delta_v_m_per_s is phasing_delta_v_m_per_s(ship, phase_forward_deg).
    print "The delta V for each burn will be " + abs(delta_v_m_per_s) + "m/s".
    add(node(start_first_burn_s, 0, 0, delta_v_m_per_s)).

    // The second burn shifts us back into our original orbit, but with a
    // new phase.
    local phasing_period_s is phasing_orbit_period_s(ship, phase_forward_deg).
    local start_second_burn_s is start_first_burn_s + phasing_period_s.
    add(node(start_second_burn_s, 0, 0, -delta_v_m_per_s)).
}