//http://www.timestocome.com

//linda Cobb

//adapted from 'numerical recipes in c'

//winter 98/99





//falling ball w/ drag calc. trajectory



#include <stdio.h>

#define NRANSI

#include "nr.h"

#include "nrutil.h"



#define N 2



static double g = 9.8;

static double c = 0.01;

static double yi = 0.0;

static double vi = +20.0;





//dydt[N] = f';

void derivs(double t, double y[], double dydt[])

{



  dydt[1] = y[2];

  dydt[2] = -g - c * y[2] * abs( y[2] );

}



int main(void)

{

	int    i, k, nsteps;

	double t, dt, *y, *dydx, *yout, j;

	double data[3][100];



	y = dvector(1,N);          //data before rk4 call

	dydx = dvector(1,N);       //derivative info

	yout = dvector(1,N);       //data after rk4 call



	k = 1;                       //loop control used when storing data

	nsteps = 0;                  //max time steps to take

	t = 0.0;                     //initial time

	dt = 0.01;                   //time step size

	j = 0.0;                     //data storage loop control

	

	y[1] = yi;                   //initial position

	y[2] = vi;                   //initial velocity

	derivs(t, y, dydx);          //initial value of derivatives



	while(nsteps < 5000){



	  t += dt;

	  rk4( y, dydx, N, t, dt, yout, derivs);

	  y[1] = yout[1];

	  y[2] = yout[2];

	

	  nsteps++;

	 

	  if(j >= 1){

	    data[0][k] = k;

	    data[1][k] = yout[2];

	    data[2][k] = yout[1];

	    j = 0.0;

	     if(((data[1][k-1]-data[1][k])<0.01)&&(nsteps > 500))break;

	   k++;

	  }

	  

	  j+= dt;



	}

	for(i=0; i<k; i++){

	    printf( "\n time(sec) %.1lf, vel. %.3lf, dist. %.1lf", data[0][i], data[1][i], data[2][i]);

	  }







	free_dvector(yout,1,N);

	free_dvector(dydx,1,N);

	free_dvector(y,1,N);

	

	printf("\n");



	return 0;

}

#undef NRANSI