simple-sim.zig

// via o1
const std = @import("std");

// Some physical constants (SI units).
// NOTE: Values are approximate for demonstration.
const G = 6.67430e-11;           // Universal Gravitational Constant (m^3 kg^-1 s^-2)
const M_EARTH = 5.97219e24;      // Mass of Earth (kg)
const R_EARTH = 6_371_000.0;     // Mean radius of Earth (meters)
const MU_EARTH = G * M_EARTH;    // Standard gravitational parameter of Earth (m^3 / s^2)
const M_MOON = 7.3477e22;        // Mass of the Moon (kg) (for future expansions)
const R_MOON_ORBIT = 384_400_000.0; // Average Earth-Moon distance (m)

// Default simulation step
const DEFAULT_TIME_STEP = 0.1;   // seconds

// A small epsilon for floating comparisons
const EPSILON = 1e-9;

// Simple struct to hold skyhook parameters
pub const SkyhookParams = struct {
    // Tether geometry
    tetherLength: f64,    // length of tether in meters
    altitudeCenter: f64,  // altitude of tether's rotation center above Earth's surface (meters)
    
    // Tether material & rotation
    tetherDensity: f64,   // density of tether material (kg/m^3)
    tensileStrength: f64, // tensile strength (Pa = N/m^2)
    rotationRate: f64,    // angular velocity (rad/s)

    // Payload properties
    payloadMass: f64,     // mass of the payload (kg)
    
    // Simulation
    timeStep: f64,        // time step for integrator (seconds)
};

// A struct holding the "state" of the tether tip and payload
pub const State = struct {
    time: f64,            // simulation time (seconds)
    tetherAngle: f64,     // current rotational angle of tether about center (rad)
    payloadPosX: f64,     // x-position of payload (assuming center at Earth center = 0, for 2D)
    payloadPosY: f64,     // y-position of payload
    payloadVelX: f64,     // x-velocity
    payloadVelY: f64,     // y-velocity
};

// Result of the simulation at each step
pub const SimResult = struct {
    states: std.ArrayList(State),
};

// Check input parameters for basic physical validity
fn validateParams(params: SkyhookParams) !void {
    // Basic checks on geometry
    if (params.tetherLength <= 0.0) {
        return error.InvalidTetherLength;
    }
    if (params.altitudeCenter <= 0.0) {
        return error.InvalidAltitude;
    }

    // Check GPa or Pa for the tether
    // For example, Kevlar ~ 3.6 GPa, Dyneema ~3.6 GPa, CNT theoretical ~ 100+ GPa, etc.
    // We'll treat the user-provided value as Pascals.
    if (params.tensileStrength <= 1e6) { // e.g. if less than 1 MPa, that’s invalid for a space tether
        return error.InvalidTensileStrength;
    }

    // Check rotation rate validity
    // The tether tip speed shouldn't exceed orbital velocity at that altitude (rough rule).
    // Orbital velocity at altitude h is sqrt(MU_EARTH / (R_EARTH + h)).
    // Tether tip speed ~ rotationRate * (tetherLength/2).
    const orbitalV = std.math.sqrt(MU_EARTH / (R_EARTH + params.altitudeCenter));
    const tetherTipSpeed = params.rotationRate * (params.tetherLength / 2.0);
    if (tetherTipSpeed > orbitalV * 2.0) {
        // If the tip speed is 2x the orbital speed, that’s suspiciously high.
        return error.TipSpeedTooHigh;
    }

    // Basic check on density
    if (params.tetherDensity <= 0.0) {
        return error.InvalidDensity;
    }

    // Check time step
    if (params.timeStep <= 0.0) {
        return error.InvalidTimeStep;
    }
}

// Estimate tension in the tether due to rotation and gravity.
// This is an extremely simplified approach that looks at the
// maximum tension at the tether root (center) if spinning in uniform circular motion.
// T = mass_of_tether_segment * (omega^2 * r) + mg(r).
// We do a rough integrative or approximate approach for demonstration.
fn estimateTetherMaxTension(params: SkyhookParams) f64 {
    // Approximate tether as small segments from center to tip.
    var steps = 100;
    var dr = params.tetherLength / @intToFloat(f64, steps);
    var maxTension = 0.0;

    var r = 0.0;
    while (r < params.tetherLength) : (r += dr) {
        // mass of small segment
        const segmentVolume = dr * 1.0 * 1.0; // assume cross-sectional area = 1 m^2 for simplicity
        // better approach: cross-sectional area = design based on stress constraints, etc.
        const segmentMass = segmentVolume * params.tetherDensity;

        // centripetal force on that segment
        const centripetalF = segmentMass * (params.rotationRate * params.rotationRate) * r;

        // gravitational force at that radius from Earth's center
        // distance from Earth center = R_EARTH + altitudeCenter ± r
        // we approximate that the difference is not huge for tension calculation (rough).
        const distance = R_EARTH + params.altitudeCenter + r;
        const gravityF = segmentMass * (MU_EARTH / (distance * distance));

        // total force
        const f = centripetalF - gravityF; 
        // note: direction is slightly more complicated, but we'll do a rough scalar approach.

        // Accumulate tension from tip down to center
        // For a real model, tension accumulates from the tip to the root. 
        maxTension += f * dr; // again, extremely simplified
    }

    return maxTension;
}

// Perform a single step of the simulation, returning updated state
fn stepSimulation(state: State, params: SkyhookParams) State {
    // We model the center of rotation of the tether at altitudeCenter above Earth's surface.
    // We'll place Earth at (0,0), tether center at (0, R_EARTH + altitudeCenter).
    // We'll treat the tether tip rotating around that point.
    // Meanwhile, the payload is at some (x,y) with some velocity, 
    // with gravitational acceleration from Earth.

    const dt = params.timeStep;

    var newState = state;

    // 1) Update tether rotation
    newState.tetherAngle = state.tetherAngle + params.rotationRate * dt;

    // 2) Gravitational acceleration on payload
    const rPayload = std.math.sqrt(state.payloadPosX * state.payloadPosX +
                                   state.payloadPosY * state.payloadPosY);
    // Avoid divide-by-zero
    if (rPayload < EPSILON) {
        // If the payload is essentially at the center, skip or clamp
        return newState;
    }
    const gAccel = MU_EARTH / (rPayload * rPayload);

    // Direction toward Earth center
    const ax = -gAccel * (state.payloadPosX / rPayload);
    const ay = -gAccel * (state.payloadPosY / rPayload);

    // 3) Update velocity from gravitational acceleration
    newState.payloadVelX = state.payloadVelX + ax * dt;
    newState.payloadVelY = state.payloadVelY + ay * dt;

    // 4) Update position
    newState.payloadPosX = state.payloadPosX + newState.payloadVelX * dt;
    newState.payloadPosY = state.payloadPosY + newState.payloadVelY * dt;

    // Additional logic: if the payload is "close" to the tether tip,
    // one might simulate a "grappling" event to attach the payload to the tether,
    // or to impart velocity to the payload. This is highly simplified.
    // Let's define "tetherTipPosX, tetherTipPosY" for the rotating tip.
    const centerX = 0.0;
    const centerY = R_EARTH + params.altitudeCenter;
    const tipR = params.tetherLength / 2.0;
    const tipAngle = newState.tetherAngle;
    const tetherTipPosX = centerX + tipR * std.math.cos(tipAngle);
    const tetherTipPosY = centerY + tipR * std.math.sin(tipAngle);

    const dx = newState.payloadPosX - tetherTipPosX;
    const dy = newState.payloadPosY - tetherTipPosY;
    const distToTip = std.math.sqrt(dx * dx + dy * dy);

    // If within capture distance, "hook" the payload
    if (distToTip < 10.0) { 
        // Very simplistic: match velocity with tether tip
        // Tether tip linear velocity: v_tip = cross(omega, radius)
        // We can approximate as the rotational velocity
        const tipSpeed = params.rotationRate * tipR;
        // direction of tip velocity is tangential to the circle
        // We'll attempt to align with that direction
        // Tip tangential direction = rotate radial vector by 90 degrees
        const vx = -tipSpeed * std.math.sin(tipAngle);
        const vy =  tipSpeed * std.math.cos(tipAngle);

        // Impart that velocity to the payload
        newState.payloadVelX = vx;
        newState.payloadVelY = vy;
    }

    newState.time = state.time + dt;
    return newState;
}

// Main simulation function
fn runSimulation(params: SkyhookParams, maxTime: f64) !SimResult {
    // Validate parameters
    try validateParams(params);

    // Check tether tension vs. tensile strength
    const tension = estimateTetherMaxTension(params);
    // cross-sectional area is "1 m^2" in our simplistic code. So stress = T / A.
    // Real designs must ensure the tether cross-sectional area is adequate to keep stress < tensileStrength.
    if (tension > params.tensileStrength) {
        return error.TetherWouldBreak;
    }

    // Initialize state
    // We place the payload at 100km altitude for "pickup"
    // i.e., at altitude = 100,000m from Earth's surface => R_EARTH + 100,000
    // Start it stationary w.r.t Earth for demonstration (or we can give it some suborbital velocity).
    var state = State{
        .time = 0.0,
        .tetherAngle = 0.0,
        .payloadPosX = 0.0,
        .payloadPosY = R_EARTH + 100_000.0,
        .payloadVelX = 0.0,
        .payloadVelY = 0.0,
    };

    var result = SimResult{
        .states = std.ArrayList(State).init(std.heap.page_allocator),
    };

    // Reserve memory to avoid frequent reallocations
    try result.states.ensureCapacity(@as(usize, maxTime / params.timeStep) + 10);

    // Main loop
    while (state.time < maxTime) {
        // Record current state
        try result.states.append(state);

        // Step simulation
        state = stepSimulation(state, params);
    }

    return result;
}

// Example usage
pub fn main() !void {
    const params = SkyhookParams{
        .tetherLength = 100_000.0,        // 100 km tether
        .altitudeCenter = 300_000.0,      // Tether center 300 km above Earth's surface
        .tetherDensity = 1_450.0,         // ~ Kevlar density (kg/m^3)
        .tensileStrength = 3.6e9,         // ~ Kevlar (Pa)
        .rotationRate = 0.001,           // rad/s (small rotation for example)
        .payloadMass = 1000.0,           // 1 tonne payload
        .timeStep = 0.5,                 // 0.5 s step
    };

    // Run for, say, 12 hours to see what happens
    var simResult = try runSimulation(params, 12.0 * 3600.0);

    // Print a few lines of results
    for (simResult.states.items) |state, i| {
        if (i % 1000 == 0) {
            std.debug.print(
                "t={}s, tetherAngle={}, payloadPos=({},{})\n",
                .{ state.time, state.tetherAngle, state.payloadPosX, state.payloadPosY }
            );
        }
    }

    // In a more advanced code, you'd check final payload velocity
    // and see if it is approaching Earth escape or a lunar transfer orbit, etc.
    // That portion would require additional orbital mechanics steps 
    // such as computing orbital elements, verifying injection conditions, etc.
}