// Project01: Physics 420/520 Spring 2007
/*
  Find the angle corresponding to the maximum horizontal range 
  of the projectile as a function of the initial height. 
  Calculate for this angle the maximum range, and time of flight. 
  Select the initial speed of the projectile as of 100m/s, 
  and the free fall acceleration as of 9.81m/s2.
*/
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
using namespace std;

int main()
{
    const double pi=3.1415926;
    const double g = 9.81;
    double y0, v0, angle, dangle, angler, x, y, t, dt;
    double range, theta, xlast, ylast, xrange;
    double ystep, ymax, tflight;

    ofstream out2disk;
    out2disk.open ("table02.dat");
    cout.setf(ios::fixed | ios::showpoint);
    cout.precision(2);    

    dt     = 0.1;           // step over time       (seconds)
    dangle = 0.1;           // step over the angle  (degrees)
    v0     = 100.0;         // initial speed        (m/s)
     
    y0     =   0.0;         // initial height       (m)
    ystep  = 100.0;         // step over height     (m)
    ymax   = 2000.0;        // max height in calculations (m)
    
    cout << setw(12) <<"y0" << setw(12) << "angle" 
         << setw(12) <<"time"<<setw(12) << "range"  << endl;
//* loop over y0
    while (y0 <= ymax)
    {
      range  = 0.0;         // initial maximum range is 0.0
//** loop over the initial angles from 0.0 to 90.0 degrees
      angle  = 0.0;
      while (angle < 90.0)
      {
        t = 0.0;            // initial time is 0.0
        x = 0.0;            // initial x position is 0.0
        y = y0;             // initial y position is y0
        angle = angle + dangle;     // angle in degrees
        angler = angle*pi/180.0;    // angle in radians
//*** loop over time till the projectile hits the ground        
        while (y >= 0.0)
        {
            t = t + dt;
            xlast = x;      // last value of x (for linear interpolation later)
            ylast = y;      // last value of y (for linear interpolation later)
            x =      v0*cos(angler)*t;
            y = y0 + v0*sin(angler)*t - g*t*t/2.0;
        }
//*** end loop over time        
        xrange = xlast + (x - xlast)*ylast/(ylast - y); //linear interpolation 
        if(range < xrange) 
            {range = xrange;
             theta = angle;
             tflight = range/(v0*cos(angler));}
       }
//** end loop over angles
//    out2disk
    cout
    << setw(12) << y0 << setw(12) << theta << setw(12) << tflight 
    << setw(12) << range << endl;
    y0 = y0 + ystep;
    }
//* end loop over y0 
    system ("pause");   
    return 0;
}