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