/*
 * bounce.c
 * author: Ari Kardasis
 * calculate next impact point and draw parabolic arc to it
 * outputs to a ps file.
 */


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


float highQuadraticRoot(float a, float b, float c){
	return (-b - sqrt(b*b - 4*a*c))/(2.0*a);
}

int main(){
	FILE *file;
	file = fopen("bounce.ps","w");

	int NSTEPS=100;
	int pathSteps = 50;
	float damp = 0.9;
	float r = 5.0;
	float border = 10.0;
	float g = -9.8;

	float xmax = 572;
	float ymax = 820;

	//initial position and velocity
	float px, py, vx, vy;
	px = 200.0;
	py = 400.0;
	vx = 100.0;
	vy = 250.0;


	//draw page

	fprintf(file, "newpath\n");
	fprintf(file, "%f %f moveto\n", border, border);
	fprintf(file, "%f %f lineto\n", border, ymax + border);
	fprintf(file, "%f %f lineto\n", xmax + border, ymax + border);
	fprintf(file, "%f %f lineto\n", xmax + border, border);
	fprintf(file, "closepath\n");
	fprintf(file, "stroke\n");

	int i;
	for (i=0; i < NSTEPS ; i++){
		//draw current pos
		fprintf(file, "0.0 setgray\n");
		fprintf(file, "newpath\n");
		fprintf(file, "%f %f %f 0 360 arc\n", px+border, py+border, r);
		fprintf(file, "closepath\n");
		fprintf(file, "stroke\n");

		fprintf(file, "newpath\n");
		fprintf(file, "%f %f %f 0 360 arc\n", px+border, py+border, r/3);
		fprintf(file, "closepath\n");
		fprintf(file, "stroke\n");


		//calculate next impact point
		float t, xtimpact, ytimpact, newpx, newpy, newvx, newvy;
		if (vx > 0){
			xtimpact = (xmax-r-px)/vx;
		}else{
			xtimpact = -(px-r)/vx;
		}
		ytimpact = highQuadraticRoot(g/2, vy, py-r);

		if(xtimpact < ytimpact){
			t = xtimpact;
			newpx = px + t * vx;
			newpy = py + vy*t + g*t*t/2;
			newvx = - damp * vx;
			newvy = vy + g*t;
		}else{
			t = ytimpact;
			newpx = px + t * vx;
			newpy = py + vy*t + g*t*t/2;
			newvx = vx;
			newvy = -damp*(vy + g*t);
		}

		//draw trajectory
		fprintf(file, "0.5 setgray\n");
		fprintf(file, "newpath\n");
		fprintf(file, "%f %f moveto\n", px+border, py+border);
		int j;
		for (j=1; j<=pathSteps; j++){
			float patht = j*t/pathSteps;
			float pathx = px + patht*vx;
			float pathy = py + vy*patht + g*patht*patht/2.0;
			fprintf(file, "%f %f lineto\n", pathx+border, pathy+border);
		}
		fprintf(file, "stroke\n");
		//set new values
		px = newpx;
		py = newpy;
		vx = newvx;
		vy = newvy;
	}

	//close things up
	fclose(file);
	return 0;
}