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

double g=32, pi=3.141593, c=0.0025; // define simulation constants

inline double fnf(double p, double q) {
  return -c*p*sqrt(p*p+q*q);
}

inline double fng(double p, double q) {
  return -c*q*sqrt(p*p+q*q)-g;
}

int main(void)
{
  double t=0; // start time
  double x=0, y=0; // initial position
  double v=160; // initial velocity
  double a=30; // fireing angle
  double h=0.05; // step size
  double m=80; // number of steps
  double k=10; // steps between prints

  a=pi*a/180; // convert angle to rad
  double p=v*cos(a), q=v*sin(a); // calc vert & horis. speeds.

  printf("t\tx\ty\tv\ta\n");
  printf("%1.1f\t%3.2f\t%3.2f\t%3.2f\t%+2.0f\n", t, x, y, sqrt(p*p+q*q), 180*atan(q/p)/pi);

  for(int n=1; n<=m; n++){
    double f1=fnf(p,q);
    double g1=fng(p,q);
    double f2=fnf(p+h*f1/2, q+h*g1/2);
    double g2=fng(p+h*f1/2, q+h*g1/2);
    double f3=fnf(p+h*f2/2, q+h*g2/2);
    double g3=fng(p+h*f2/2, q+h*g2/2);
    double f4=fnf(p+h*f3, q+h*g3);
    double g4=fng(p+h*f3, q+h*g3);

    double dp=(f1+2*(f2+f3)+f4)/6;
    double dq=(g1+2*(g2+g3)+g4)/6;

    x+=p*h+0.5*dp*h*h; // three-term-Taylor approximation
    y+=q*h+0.5*dq*h*h;

    p+=dp*h; q+=dq*h; t+=h;

    if((int)(n/k)==n/k)
      printf("%1.1f\t%3.2f\t%3.2f\t%3.2f\t%+2.0f\n", t, x, y, sqrt(p*p+q*q), 180*atan(q/p)/pi);
  }
}