#include <stdio.h>

double Mm = 0.02897; //kg/mol molar mass of air
double G  = 6.67300e-11; /* %in SI */
double R  = 8.31; /* %gas constant */
double rho_e = 5515; /* %kg/m^3, density of Earth */
double rho_h = 1400; /* %kg/m^3, density of human body */
double R0    = 6371e3; /* %m, mean Earth radius */

double Cd    = 1.3; /* %drag coeff. of human body */
double A     = 0.7; /* %m^2, incident suface area of human body */
double T     = 300; /* %K, temperature */
double p0    = 1e5; //Pa, 1 atmosphere
double m     = 80;  //kg, human weight

const double pi = 3.14159;

double f1(double x1, double x2)
{
  return x2;
}

double f2(double x1, double x2)
{
  return (-4/3*pi*G*rho_e*x1 + Cd*A/2/m*Mm/R/T*(p0 + 2/3*G*rho_e*pi*(R0*R0 - x1*x1))*x2*x2);
}

int main()
{
  double dt = 0.5; /* seconds */

  /* initial conditions */
  double x1 = R0;
  double x2 = 0;
  double t  = 0;
  while (x1 > 0)
    {
      double x1dot = f1(x1, x2);
      double x2dot = f2(x1, x2);

      double hstep1 = x1 + x1dot*dt/2;
      double hstep2 = x2 + x2dot*dt/2;

      printf("%f\t%f\t%f\n", t, x1, x2);

      x1 += f1(hstep1, hstep2)*dt;
      x2 += f2(hstep1, hstep2)*dt;
      t  += dt;
    }

  return 0;
}