/*
 *  pebble.cpp  -  Don Cross  -  http://cosinekitty.com
 *
 *  On Linux, build using:
 *  g++ -o pebble -Wall -Ofast pebble.cpp
 *
 *  Simulates the math of a pebble falling from outer space onto the Earth.
 *
 *  See:  http://cosinekitty.com/pebble
 *
 *  The goal of this program is to verify the math on that page is valid
 *  by comparing its prediction for the amount of time it takes the
 *  pebble to strike the Earth against the time determined by a numerical
 *  simulation of the pebble falling.
 */

#include <iostream>
#include <cmath>
#include <stdlib.h>

const double G = 6.67384e-11;   // universal gravitation constant [m^3 kg^(-1) s^(-2)]
const double M = 5.97219e+24;   // mass of the Earth [kg]
const double R = 6.3781e+6;     // radius of the Earth [m]


// Uses the formula at the bottom of http://cosinekitty.com/pebble
// to calculate the amount of time a pebble dropped from 
// a distance H (in meters) from the center of the Earth
// takes to strike the surface of the Earth.
// The return value is expressed in seconds.
double CalcDropTime(double H)
{
    return sqrt(H/(2*M*G)) * ((H/2) * acos((2*R - H) / H) + sqrt(R*(H - R)));
}


// Calculate the instantaneous speed of a particle that starts falling
// from distance H at the moment it passes through distance x.
// This formula derives directly from conserving kinetic energy + potential energy,
// with v=0 when x=H.
inline double Speed(double H, double x)
{
    return -sqrt(2*M*G * (1/x - 1/H));
}


// The same idea as CalcDropTime, only this is a numerical simulation.
double SimDropTime(double H)
{
    // We will keep this simulation as simple and direct as possible.
    // There is no attempt to do any fancy tricks here for optimization.
    // We just want to know if the formula in CalcDropTime is correct.

    // t = time after release of pebble from rest [s]
    // x = Distance from center of Earth [m]
    // v = speed [m/s]

    double t = 0;   // no time has elapsed (yet)
    double x = H;   // the pebble starts at distance H from center of Earth
    double v = 0;
    while (x > R) 
    {
        // We need to choose a delta-x value that is large enough
        // that this loop won't go forever, but small enough that
        // the answer will be reasonably accurate.
        double dx = 0.01 * v;    // approximate how far we would fall in 10 milliseconds
        if (-dx < 0.001)
        {
            dx = -0.001;        // but fall at least 1 mm each time, especially when speed is very low
        }

        // Calculate the instantaneous speed at the midpoint between x and x+dx.
        // We will use this as the average speed over the interval x ... x+dx.
        v = Speed(H, x + dx/2);

        // Estimate the change in time dt using dx/dt = v.
        // Add up each increment of time.
        t += dx / v;

        // Update position x for the next iteration.
        x += dx;
    }

    // In general, when we exit the loop, the pebble has gone below the surface: x < R.
    // We need to "back up" to correct the overshoot and find the exact moment of impact.
    double xerror = x - R;          // how far underground we went, in negative meters
    v = Speed(H, R + xerror/2);     // average speed for the underground movement (also negative)
    t -= xerror / v;                // subtract out the amount of time spent underground

    return t;
}


int main(int argc, const char *argv[])
{
    using namespace std;

    int rc = 1;

    if (argc == 2)
    {
        double height = atof(argv[1]);
        if (height > 0.0)
        {
            cout << endl;
            cout << "Height = " << height << " meters." << endl;
            cout << endl;

            // The variable H we need is actually the distance from the center of the Earth,
            // NOT the height above the ground.

            double H = height + R;   // add Earth radius to get distance from center of Earth.

            // For fun, show the impact speed, since we need that function anyway.
            double vi = Speed(H, R);
            cout << "Impact speed = " << vi << " m/s." << endl;
            cout << endl;
            cout << "Time to impact:" << endl;

            // Use my derived formula for calculating impact time.
            double t1 = CalcDropTime(H);
            cout << "    Calculated = " << t1 << " seconds = " << (t1/(60*60*24)) << " days." << endl;

            // Run a numerical simulation of the pebble falling until it hits the Earth.
            double t2 = SimDropTime(H);
            cout << "    Simulated  = " << t2 << " seconds." << endl;
            double error = t2 - t1;
            cout << "    Error      = " << error << " seconds" << endl;
            cout << endl;

            rc = 0;
        }
        else
        {
            cout << "The height must be a positive real number of meters." << endl;
        }
    }
    else
    {
        cout << "USAGE:  pebble height" << endl;
        cout << endl;
        cout << "Where 'height' is the distance above the Earth in meters." << endl;
        cout << endl;
    }

    return rc;
}