/*
  Randall J. Scalise - 26 October 1997
  Given the time increment dt, numerically integrate to find the distance
  a particle has fallen in gravity with aerodynamic drag -cv^2.
*/

#include<stdio.h>
#include<math.h>

#define m   70.0       /* mass of the particle in kg             */
#define c    0.3       /* drag constant (dimensionless)          */
#define g    9.8       /* acceleration due to gravity in m/s^2   */
#define T   10         /* time of fall in seconds                */
#define tol  0.00001   /* tolerance for comparing times          */

double theory(double t)
{
double vterm = sqrt(m*g/c);
return vterm * (t + vterm/g*log((1 + exp(-2*g*t/vterm))/2));
}


main(argc, argv)
int argc;
char **argv;
{
if(argc==1||argc>2) fprintf(stderr, "Usage: %s <dt>\n", argv[0]);
else
  {
  double dt = atof(argv[1]);
  int i;
  double t=0,x=0,v=0,a;

  printf("time increment dt = %f\n",dt);
  printf("t(sec)  x(meters)     theory\n");

  for(i=1; i<=T/dt; i++)
    {
    a = -c*v*v/m + g;

    t = t + dt;
    x = x + v*dt;
    v = v + a*dt;

    if(fabs(t-floor(t+tol))<=tol) 
      printf("%2.0f      %8.4f    %8.4f\n",t,x,theory(t));
    }
  }
}