Simulate ascents in KSP

ksp_ascent_sim.js

var THRUST = 5;
var PLANET_RADIUS = 200000;
var PLANET_MASS = 9.7600236e20;
var PLANET_SIDEREAL_VELOCITY = 9.0416;
var G = 6.67384e-11;
var ISP = 350;
var FUEL_RATE = THRUST / (ISP * 9.81);
var TARGET_ALT = 40000;
var DRY_MASS = 0.39875;
var FUEL_MASS = 0.191025;
var START_TURN = 1000;
var TURN_RATE = -0.174532925; // 10 degrees per second

function gravity(state) {
    var r = PLANET_RADIUS + state.h;
    var F = G * ((PLANET_MASS * state.mass) / (r * r));
    return F / state.mass;
}

function calculateAngle(state, dt) {
    // Vector of our current position
    var ax = 0;
    var ay = PLANET_RADIUS + state.h;

    // Vector of our change in position
    var bmag = state.v * dt;
    var btheta = state.theta;
    var bx = bmag * Math.cos(btheta);
    var by = bmag * Math.sin(btheta);

    // Vector of our new position
    var rx = ax + bx;
    var ry = ay + by;
    var rtheta = Math.atan2(ry, rx);

    return Math.PI / 2 - rtheta;
}

/**
 * Return the change in velocity away from the planet
 */
function acceleration(state) {
    var a = THRUST / state.mass;

    // Acceleration components
    var ax = a * Math.cos(state.pitch);
    var ay = a * Math.sin(state.pitch);
    ay -= gravity(state); // Gravity is a bitch

    // Velocity components
    var vx = state.v * Math.cos(state.theta);
    var vy = state.v * Math.sin(state.theta);

    // Calculate the new velocity vector
    var rx = ax + vx;
    var ry = ay + vy;

    // Convert to polar form
    var magnitude = Math.sqrt(rx * rx + ry * ry);
    var theta = Math.atan2(ry, rx);

    return {
        "dv": magnitude - state.v,
        "dtheta": theta - state.theta
    };
}

function evaluate(initialState, t, dt, d) {
    if (evaluate.arguments.length === 2) {
        var a = acceleration(initialState);
        var output = {};
        output.dh = initialState.v;
        output.dv = a.dv;
        output.dtheta = a.dtheta;
        return output;
    } else {
        var state = {};
        state.pitch = initialState.pitch;
        state.h = initialState.h + d.dh * dt;
        state.v = initialState.v + d.dv * dt;
        state.theta = initialState.theta + d.dtheta * dt;
        state.mass = initialState.mass - FUEL_RATE * dt;
        state.fuelMass = initialState.fuelMass - FUEL_RATE * dt;
        var a = acceleration(state);
        var output = {};
        output.dh = state.v;
        output.dv = a.dv;
        output.dtheta = a.dtheta;
        return output;
    }
}

function integrate(state, t, dt) {
    var a = evaluate(state, t);
    var b = evaluate(state, t, dt * 0.5, a);
    var c = evaluate(state, t, dt * 0.5, b);
    var d = evaluate(state, t, dt, c);

    var dhdt = 1.0/6.0 * (a.dh + 2.0 * (b.dh + c.dh) + d.dh);
    var dvdt = 1.0/6.0 * (a.dv + 2.0 * (b.dv + c.dv) + d.dv);
    var dtheta_dt = 1.0/6.0 * (a.dtheta + 2.0 * (b.dtheta + c.dtheta) + d.dtheta);

    state.h += dhdt * dt;
    state.v += dvdt * dt;
    state.theta += dtheta_dt * dt;
    state.mass -= FUEL_RATE * dt;
    state.fuelMass -= FUEL_RATE * dt;
    state.angle += calculateAngle(state, dt);

    // Perform turn
    if (state.h >= START_TURN) {
        if (state.pitch >= 0.0) {
            state.pitch = state.pitch + TURN_RATE;
            // Clamp to horizon
            if (state.pitch < 0.0) {
                state.pitch = 0.0;
            }
        }
    }

    var ap = ap_sim(state);
    state.eta = ap.t;
    state.ap = ap.h;
    state.totalAngle = ap.angle;

    return state;
}

/**
 * Return the change in velocity from freefall
 */
function freefall(state) {
    // Acceleration components
    var ax = 0;
    var ay = -gravity(state); // Gravity is a bitch

    // Velocity components
    var vx = state.v * Math.cos(state.theta);
    var vy = state.v * Math.sin(state.theta);

    // Calculate the new velocity vector
    var rx = ax + vx;
    var ry = ay + vy;

    // Convert to polar form
    var magnitude = Math.sqrt(rx * rx + ry * ry);
    var theta = Math.atan2(ry, rx);

    return {
        "dv": magnitude - state.v,
        "dtheta": theta - state.theta
    };
}

function ap_evaluate(initialState, t, dt, d) {
    if (ap_evaluate.arguments.length === 2) {
        var a = freefall(initialState);
        var output = {};
        output.dh = initialState.v;
        output.dv = a.dv;
        output.dtheta = a.dtheta;
        return output;
    } else {
        var state = {};
        state.h = initialState.h + d.dh * dt;
        state.v = initialState.v + d.dv * dt;
        state.theta = initialState.theta + d.dtheta * dt;
        state.mass = initialState.mass;
        var a = freefall(state);
        var output = {};
        output.dh = state.v;
        output.dv = a.dv;
        output.dtheta = a.dtheta;
        return output;
    }
}

function ap_integrate(state, t, dt) {
    var a = ap_evaluate(state, t);
    var b = ap_evaluate(state, t, dt * 0.5, a);
    var c = ap_evaluate(state, t, dt * 0.5, b);
    var d = ap_evaluate(state, t, dt, c);

    var dhdt = 1.0/6.0 * (a.dh + 2.0 * (b.dh + c.dh) + d.dh);
    var dvdt = 1.0/6.0 * (a.dv + 2.0 * (b.dv + c.dv) + d.dv);
    var dtheta_dt = 1.0/6.0 * (a.dtheta + 2.0 * (b.dtheta + c.dtheta) + d.dtheta);

    state.h += dhdt * dt;
    state.v += dvdt * dt;
    state.theta += dtheta_dt * dt;
    state.angle += calculateAngle(state, dt);

    return state;
}

/* Simulate to AP of orbit */
function ap_sim(initialState) {
    var state = {};
    state.mass = initialState.mass;
    state.h = initialState.h;
    state.v = initialState.v;
    state.theta = initialState.theta;
    state.angle = initialState.angle;
    var dt = 1.0;
    var t = 0;

    while (state.v * Math.sin(state.theta) >= 0.0) {
        state = ap_integrate(state, t, dt);
        //console.log(JSON.parse(JSON.stringify(state)));
        t += dt;
    }

    return { "t": t, "h": state.h, "angle": state.angle };
}

function sim() {
    var totalMass = DRY_MASS + FUEL_MASS;
    var state = {};
    state.mass = totalMass;
    state.fuelMass = FUEL_MASS;
    state.h = 0;
    state.pitch = Math.PI / 2; // Straight up, 90 degrees from horizon
    state.theta = 0; // Travelling along horizon while on the surface
    state.v = PLANET_SIDEREAL_VELOCITY; // at the speed that the planet rotates
    state.ap = 0;
    state.angle = 0;
    var t = 0;
    var dt = 0.1;

    while (state.ap <= TARGET_ALT && state.fuelMass >= 0.0) {
        //console.log(JSON.parse(JSON.stringify(state)));
        state = integrate(state, t, dt);
        t += dt;
    }

    state.t = t;
    state.deltaV = Math.log(totalMass / state.mass) * ISP * 9.81;
    return state;
}