#include
#include
#define NDIM (2)
typedef struct {
double pos[NDIM];
double vel[NDIM];
} body_t;
const double G = 6.67428e-8; /* cm^3/g/s^2 */
void initargs(int argc, char *argv[]);
double acc(double m, double oneor);
int main(int argc, char *argv[])
{
FILE *ofp;
double m, r, oneor, v, a, t, dt = 1.0;
unsigned int i;
body_t body;
initargs(argc, argv);
ofp = fopen(argv[1], "w");
m = strtod(argv[2], NULL);
r = strtod(argv[3], NULL);
v = strtod(argv[4], NULL);
if (v == 0.0) v = sqrt(G*m/r);
body.pos[0] = r;
body.pos[1] = 0.0;
body.vel[0] = 0.0;
body.vel[1] = v;
for (t = 0.0, i = 0; t < 1.844e6; ++i, t += dt) {
oneor = 1.0/sqrt(body.pos[0] * body.pos[0] +
body.pos[1] * body.pos[1]);
a = acc(m, oneor);
body.vel[0] += dt * a * body.pos[0] * oneor;
body.vel[1] += dt * a * body.pos[1] * oneor;
body.pos[0] += dt * body.vel[0];
body.pos[1] += dt * body.vel[1];
if (i % 10000 == 0) {
fprintf(ofp, "%g\t%g\n", (t / 60) - 2050, body.pos[0]);
}
}
fclose(ofp);
return 0;
}
void initargs(int argc, char *argv[])
{
if (argc != 5) {
fprintf(stderr, "Usage: %s outfile M r v\n", argv[0]);
exit(1);
}
}
double acc(double m, double oneor)
{
return -G*m*oneor*oneor;
}
the orbit code we ran at earthwatch!
