// N-Body simulation benchmark - floating point compute
// Simplified 2-body problem

fn simulate(steps: Int): Float = {
    // Two bodies: sun and planet
    let sunMass = 1000.0
    let planetMass = 1.0

    // Initial positions and velocities
    simulateLoop(steps, 0.0, 100.0, 0.0, 6.28, 0.01)
}

fn simulateLoop(steps: Int, x: Float, y: Float, vx: Float, vy: Float, dt: Float): Float = {
    if steps == 0 then sqrt(x * x + y * y)
    else {
        // Calculate distance and force
        let r2 = x * x + y * y
        let r = sqrt(r2)
        let f = 1000.0 / r2  // G * M / r^2

        // Acceleration
        let ax = 0.0 - f * x / r
        let ay = 0.0 - f * y / r

        // Update velocity
        let nvx = vx + ax * dt
        let nvy = vy + ay * dt

        // Update position
        let nx = x + nvx * dt
        let ny = y + nvy * dt

        simulateLoop(steps - 1, nx, ny, nvx, nvy, dt)
    }
}

fn sqrt(x: Float): Float = {
    sqrtIter(x, x / 2.0, 0)
}

fn sqrtIter(x: Float, guess: Float, iter: Int): Float = {
    if iter > 20 then guess
    else {
        let newGuess = (guess + x / guess) / 2.0
        sqrtIter(x, newGuess, iter + 1)
    }
}

fn main(): Unit = {
    let finalDistance = simulate(100000)
    Console.print("Final orbital distance: " + toString(finalDistance))
}