// 
// Torbert, 8.20.2008, revised 8.14.2009
// 
// Projectile Motion with Air Resistance, Version 1.0
// 
//   Input: none, but wind's velocity and launch angle must be specified
// Process: uses turbulent flow model (high velocity) for air resistance
//          contribution to acceleration (thus, v^2 term) and explicit
//          Euler for time-marching (like left-hand rule for approximating
//          integrals) until projectile hits the ground (y<=0.0)
//  Output: position (x,y) of projectile for each time step, for plotting
// 

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

double sign(double arg) // example of a funtion (method)
{
	if(arg>0.0)
		return -1.0;
	else
		return 1.0;
}
double Fax(double vw,double vx,double rf) // turbulent flow
{
	return sign(vx-vw)*(vx-vw)*(vx-vw)*rf;
}
double Fay(double g,double vy,double rf)
{
	return -g+sign(vy)*vy*vy*rf;
}
int main(int argc, char* argv[])
{
	double	m      = 1.0;
	double	v0     = 134.11200;			// meters/sec = 300.0 mph
	double	rf     = 0.0034;				// air resistance factor
	double	pi     = 3.1415926;
	double	g      = 9.80665;
	double	dt     = 0.1;
	double	vw,theta,t;
	double	x,vx,ax;
	double	y,vy,ay;

	vw        = -10.0;                  // headwind
	theta     = pi/4.0;
	y         = 0.0;

	// vw        = 1.0;                 // freefall
	// theta     = 0.0;
	// v0        = 0.0;
	// y         = 2500.0;
	
	vx        = v0*cos(theta);
	vy        = v0*sin(theta);
	t = x     = 0.0;
	do{
		ax     = Fax(vw,vx,rf)/m;
		ay     = Fay(g,vy,rf)/m;
		vx    += ax*dt;                  // explicit Euler
		vy    += ay*dt;
		x     += vx*dt;
		y     += vy*dt;
		t     += dt;
		printf("%9.2f %9.2f\n",x,y);
		// printf("%9.2f %9.2f\n",t,vy); // freefall, terminal velocity
	}while(y>0.0);

	return 0;
}