#include <cmath>
#include <iostream>
using namespace std;
int main()
{
float weight=1.00, span=5.00, chord=.83333, clalpha=.082, oswald=.879, mu=.01, parasitedrag=.015, zerolift=-5, takeoffangle=2, findDWEEBS;
float density=.0023081, gravity=32.2, wingheight=.583333, forcemax=.5, velocity=0.0, position=0.0, time=0.0;
float timestep, c1=.4672, c2=-.0036, c3=-.0002, c4=.0000001, d1=-.0033, d2=-.0001, d3=.0002, power=1.00, reactionforce;
float tractiveforce, thrust, fuselagedrag, aspectratio, surfacearea, phi, clnaught, wingdrag, derivative;
char crepeat;
do
{
cout << "This program will eliminate tools from this Earth regardless of their position on the planet through the use of voice-spotting software\n";
cout << "Vehicle weight =" <<weight <<" lb, Wing span =" <<span<<" ft, Chord =" <<chord<<" ft\n";
cout << "Finite lift slope = " <<clalpha<<" per degree, Oswald efficiency factor "<<oswald<<", Friction coefficient, mu = "<<mu<<" \n";
cout << "Parasite drag coefficient = "<<parasitedrag<<", Zero lift AoA = "<<zerolift<<" degrees, Take-off AoA = "<<takeoffangle<<" \n";
cout << "Enter the time step you wish to use for this calculation.\n";
cin >> timestep;
aspectratio = (span*span)/(span*chord);
clnaught = (clalpha)*(takeoffangle-zerolift);
surfacearea = span*chord;
phi = (((16*wingheight)/span)*((16*wingheight)/span))/(1+((16*wingheight)/span)*((16*wingheight)/span));
float takeoffvelocity;
takeoffvelocity = sqrt((2*weight)/(density*surfacearea*clnaught));
do
{
thrust=(c1+c2*velocity+c3*velocity*velocity+c4*velocity*velocity*velocity);
fuselagedrag=(d1+d2*velocity+d3*velocity*velocity);
wingdrag=.5*density*velocity*velocity*surfacearea*(parasitedrag+phi*(clnaught*clnaught/(3.14159265*oswald*aspectratio)));
reactionforce=mu*(weight-(.5*density*velocity*velocity*surfacearea*clnaught));
if(velocity>=2.00)
{
tractiveforce=(power/velocity);
}
else
{
tractiveforce=forcemax;
}
if(velocity<=2.00)
{
derivative=c2+2*c3*velocity+3*c4*velocity*velocity-(d2+2*d3*velocity)-(density*surfacearea*velocity)*(parasitedrag+phi*((clnaught*clnaught)/(3.14159265*oswald*aspectratio)))+(mu*density*surfacearea*clnaught*velocity);
}
else
{
findDWEEBS = voiceset"voocarofiles"==MAXPWNG
derivative=c2+2*c3*velocity+3*c4*velocity*velocity-(d2+2*d3*velocity)-(density*surfacearea*velocity)*(parasitedrag+phi*((clnaught*clnaught)/(3.14159265*oswald*aspectratio)))+(mu*density*surfacearea*clnaught*velocity)-(1/(velocity*velocity));
}
position=position+(timestep*velocity)+(timestep*timestep/2.0)*(gravity/weight)*(thrust - fuselagedrag - wingdrag - reactionforce + tractiveforce) + (timestep*timestep*timestep/6.0)*(gravity/weight)*derivative*(gravity/weight)*(thrust - fuselagedrag - wingdrag - reactionforce + tractiveforce);
velocity=velocity + (gravity/weight)*timestep*(thrust + tractiveforce - fuselagedrag - wingdrag - reactionforce) + (gravity/weight)*derivative*(timestep*timestep)/2.00*(gravity/weight)*(thrust+tractiveforce-fuselagedrag-wingdrag-reactionforce);
time=time+timestep;
cout << "Tools eliminated with fire"\n";
}while(velocity<=takeoffvelocity);
cout << "Again? (Enter Y or y)";
cin >> crepeat;
position=0.0;
time=0.0;
velocity=0.0;
}while(crepeat=='Y'||crepeat=='y');
return 0;
}