// Transcendental equation solution
// C. Thaxton
//
// This is the solution to exercise 5.7.
//
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "comphys.h"
#include "comphys.c"

#define pi 3.14159265358979323846264338327950288
#define g 9.81
#define v0 600

double theta=60.0;
double x (double t);
double y (double t);
double T0();
double findT (double k);

double T,k,k1=0.0,k2=1.0;
double Rtarget=1500.0;
double R;
double Ttol=1e-12;
double Rtol=1e-12;
int i,j;

int main() 
{
   theta=theta*pi/180.0;
   printf("\n********************************************************");
   printf("\nThis program computes first the range and time of flight for a");
   printf("\nprojectile given a damping coefficient, k, of 0.005.");
   printf("\nIt then finds the k value required for a range of %lf [m].\n",Rtarget);
   printf("\nFor this problem...");
   printf("\n  v-initial=%lf [m/s]",v0);
   printf("\n  theta=%lf degrees\n",theta/pi*180.0);

// Part 1: set k=0.005 and find T

   k=0.005;
   T=findT(k);
   printf("\nRESULTS:\n");
   printf("\nThe time of flight for k=0.005 is: %lf [s],",T);
   printf("\nwith a Range of %lf [m].",x(T));
   printf("\nThis took %d iterations to complete.\n\n",i);

// Part 2: Find the k that allows R=1500

   R=x(0);j=0;
   while ( fabs(R-Rtarget) > Rtol )
   {
       k=(k1+k2)/2.0;j++;
       R=x(findT(k));
       if ((R-Rtarget) > 0.0) k1=k;
       else k2=k;
   }
   printf("For a range of %lf [m],\n",R);  
   printf("this requires a 'k' value of %lf.",k); 
   printf("\nThis took %d iterations to complete.\n\n",j);    
   return(0);
}

double x (double t)
{
   return(v0*cos(theta)*(1.0-exp(-k*t))/k);
}

double y (double t)
{
   return( -g*t/k + (k*v0*sin(theta)+g)*(1.0-exp(-k*t))/pow(k,2.0));
}

double T0 ()
{
   return(2*v0*sin(theta)/g);
}

double findT (double k)
{  
   double Tnew,Told=0.0;
   Tnew=T0();i=0;
   while (fabs(Tnew-Told) > Ttol) {
      Told=Tnew;i++;
      Tnew=(k*v0*sin(theta)+g)*(1-exp(-k*Told))/g/k;
   }   
   return(Tnew);
}