/*
* 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
#include
#include
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;
}