/* a simulation of a cannonball flight, with subsonic drag */

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

float d2r(float);

int main(int argc, char *argv[]){

	float x, y, maxy = 0.0, xv, yv, deltaxv, deltayv, t, deltat, angle, velocity;
	float mass = 2.7, cd = 0.44, area = 0.0102070, density = 1.204, force_gravity_proj = 26.46;
	float force_x, drag_force_y, force_y;

	/* check for proper usage */
	if (argc != 4) {
		printf("Usage: cannon_drag <launch angle in degrees> <launch velocity in meters/second> <timestep in seconds>\n");
		exit(1);
	}

	/* parse variables */
	angle = atof(argv[1]);
	velocity = atof(argv[2]);
	deltat = atof(argv[3]);

	/* print out input */
	printf("You entered:\n");
	printf("Launch angle: %f\n", angle);
	printf("Launch velocity: %f\n", velocity);
	printf("Timestep: %f\n\n", deltat);

	/* calculate vectors */
	xv = cos( d2r (angle)) * velocity;
	yv = sin( d2r (angle)) * velocity;
	printf("Vectors:\n");
	printf("X velocity: %f\n", xv);
	printf("Y velocity: %f\n\n", yv);
	
	/* set initial conditions */
	x = 0.0;
	y = 0.00001; /* make sure it's not zero initially */
	t = 0.0;

	/* print headder */
	printf("Time, X position, X velocity, Y position, Y velocity\n");

	/* main loop */
	while (y > 0.0) {
		/* print the current state */
		printf("%f, %f, %f, %f, %f\n", t, x, xv, y, yv);

		/* update the forces, accelerations, and speeds */
		force_x = 0.5 * density * xv * xv * area * cd;
		xv = xv - (force_x / mass) * deltat;

		/* drag force in the Y dimension changes direction! */
		drag_force_y = 0.5 * density * yv * yv * area * cd;
		if (yv > 0.0) {
			drag_force_y = drag_force_y * -1.0;
		} 

		force_y = (-1.0 * force_gravity_proj) + drag_force_y;
		deltayv = (force_y / mass) * deltat;

		yv = yv + deltayv;

		/* update the other state variables based on the new speeds */
		t = t + deltat;
		x = x + xv * deltat;
		y = y + yv * deltat;
		if (y > maxy)
			maxy = y;
	}

	/* print some final numbers */
	printf("\n\n");
	printf("Maximum height: %f\n", maxy);
	printf("Maximum distance: %f\n", x);
	printf("Time to impact: %f\n", t);
}

float d2r(float degrees) {

	float radians;

	/* convert degrees to radians for trig calcs */
	radians = degrees * M_PI / 180.0;
	return(radians);
}