C++ Distance Calculator

[B3anss]

I’ve been working on this c++ program for a little while and I’ve finally got it to run. Unfortunately, it does not work and the values seem to be not changing. Does anyone see what’s wrong? Thanks

#include <iostream>
#include <cmath>
#include <iomanip>
#include <vector>
#include <math.h>

using namespace std;

///Use structs possibly to not have to declare all these variables globally

vector<double> acceleration(double, double, double, double, double);

vector<double> distance_final(double, double, double, double, double, double, double);

//initial condition
double mass = 5.125; //Oz
double circumference = 9.125; //inches
double x0 = 0; //ft
double y0 = 2.000; //ft
double z0 = 3.0; //ft
double exit_speed = 100.0; //mph
double launch_angle = 30.0; //degree
double direction = 0; //degree
double sign = 1.0;
double backspin = -763 + 120*launch_angle + 21*direction * sign;//rpm
double sidespin = (-1*sign)*849-94*direction; //rpm
double wg = 0;//rpm
double tau = 25.0; //sec
double dt = 0.01;
double pi = 3.1415;
double e = 2.7182818;


double Temp_F = 78.0; //F
double elev_ft = 0; //ft
double vwind = 0; //mph
double phiwind = 0; //degree
double hwind = 0; //mph
double relative_humidity = 50.0; //%
double barometric_pressure = 29.92; //Hg

double Temp_C = (5.0/9.0)*(Temp_F-32.0);
double elev_m = elev_ft / 3.2808;

double vxw = vwind*1.467*sin(phiwind*pi/180.0); //ft/s wind speed in x direction
double vyw = vwind*1.467*cos(phiwind*pi/180.0); //ft/s wind speed in y direction

double beta = 0.0001217; //1/m
double SVP = 4.5841*(pow(e,((18.687-Temp_C /234.5)*Temp_C/(257.14+Temp_C )))); //barometric calculations incorporated into distance and involving temperature
double barometric_pressure_mmHg = barometric_pressure*1000.0/39.37;
double rho_kg_m3 = 1.2929*(273/(Temp_C+273.0)*(barometric_pressure_mmHg*(pow(e,(-beta*elev_ft)))-0.3783*relative_humidity*SVP/100.0)/760.0); //kg/m^3
double rho_lb_ft3 = rho_kg_m3*0.06261; //lb/ft^3
// constant value, non editable variables
double c0 = 0.07182*rho_lb_ft3*(5.125/mass)*(pow(2, (circumference/9.125)));

double v0 = 1.467*exit_speed;
double v0x = v0*cos(launch_angle*pi/180.0)*sin(direction*pi/180.0);
double v0y = v0*cos(launch_angle*pi/180.0)*cos(direction*pi/180.0);
double v0z = v0*sin(launch_angle*pi/180.0);
double wx = (backspin*cos(direction*pi/180.0)-sidespin*sin(launch_angle*pi/180.0)*sin(direction*pi/180.0)+wg*v0x/v0)*pi/30.0;
double wy = (-backspin*sin(direction*pi/180.0)-sidespin*sin(launch_angle*pi/180.0)*cos(direction*pi/180.0)+wg*v0y/v0)*pi/30.0;
double wz = (sidespin*cos(launch_angle*pi/180.0)+wg*v0z/v0)*pi/30.0;
double omega = sqrt(wx*wx+wy*wy+wz*wz);
double romega = (circumference/2.0/pi)*omega/12.0;




///All the constants/variables declared above...

int main()
{


exit_speed = 100.0; //mph
launch_angle = 30.0; //degree
direction = 0.0; //degree
sign = 1.0;
Temp_F = 78.0; //F
elev_ft = 0; //ft
vwind = 0; //mph
phiwind = 0; //degree
hwind = 0; //mph
relative_humidity = 50.0; //%


vector<int> x_path_RK4; ///Making the path in the x, y, and z coordinates.
vector<int> y_path_RK4;
vector<int> z_path_RK4;

double x = 0;
double y = y0;
double z = z0;

double vx = v0x;
double vy = v0y;
double vz = v0z;
double bruhArray[6] = {};

for (int t = 1; t < 1000; t++)3
{

x_path_RK4.push_back(x);
y_path_RK4.push_back(y);
z_path_RK4.push_back(z);



bruhArray[0] = distance_final(x,y,z,vx,vy,vz,t)[0];
bruhArray[1] = distance_final(x,y,z,vx,vy,vz,t)[1];
bruhArray[2] = distance_final(x,y,z,vx,vy,vz,t)[2];
bruhArray[3] = distance_final(x,y,z,vx,vy,vz,t)[3];
bruhArray[4] = distance_final(x,y,z,vx,vy,vz,t)[4];
bruhArray[5] = distance_final(x,y,z,vx,vy,vz,t)[5];


cout << " & " << bruhArray[2] << " & " << endl;
if(z <= 0)
{
t = dt * t;
cout << "Hang time: " << t << endl;
cout << "Distance Traveled: " << y << endl;
break;
}
}

///Sfml code down here...
///Use vector math using horizontal angle and magnitude(distance) of ball vector find the coordinate points in the x and y direction


}

vector<double> acceleration(double t, double vx, double vy, double vz, double z)
{
double v = sqrt(vx*vx+vy*vy+vz*vz);
double vw, vyw_sim;
if(z > hwind)
{
vw = sqrt((pow(2, (vx-vxw))) + (pow(2, (vy-vyw))) + (pow(2, vz)));
vyw_sim = vyw;
}
else
{
vw = v;
vyw_sim = 0;
}

double S = (romega/vw)*(pow(2, (-t/(tau*146.7/v))));
double Cd = 0.4105*(1.0+0.2017*S*S);
double Cl = 1.0/(2.32+0.4/S);
double vxw_sim = 0;
double adragx = -c0*Cd*vw*(vx-vxw);
double adragy = -c0*Cd*vw*(vy-vyw_sim);
double adragz = -c0*Cd*vw*vz;
double w = omega*(pow(e, (-t/tau)))*30.0/pi;
double aMagx = c0*(Cl/omega)*vw*(wy*vz-wz*(vy-vyw));
double aMagy = c0*(Cl/omega)*vw*(wz*(vx-vw)-wx*vz);
cout << " { " << aMagy << " } " << endl;
double aMagz = c0*(Cl/omega)*vw*(wx*(vy-vyw)-wy*(vx-vxw));
double ax = adragx+aMagx;
double ay = adragy+aMagy;
double az = adragz+aMagz-32.174;


vector<double> myArray;
myArray.push_back(ax);
myArray.push_back(ay);
myArray.push_back(az);

return myArray;

}


vector<double>distance_final(double x, double y, double z, double vx, double vy, double vz, double t)
{
///Final distances are calculated here in this function using all previous variable
vector<double> result = acceleration(t,vx,vy,vz,z); ///Result of the acceleration function, gets three values from the acceleration function and puts them into this vector
double a1x = result[0];
double a1y = result[1];
double a1z = result[2];
double k1x = dt*a1x;
double k1z = dt*a1z;
double k1y = dt*a1y;



vector<double> result2 = acceleration(t+dt/2.0,vx+k1x/2.0,vy+k1y/2.0,vz+k1z/2.0,z);
double a2x = result2[0];
double a2y = result2[1];
double a2z = result2[2];
double k2x = dt*a2x;
double k2y = dt*a2y;
double k2z = dt*a2z;


vector<double> result3 = acceleration(t+dt/2.0,vx+k2x/2.0,vy+k2y/2.0,vz+k2z/2.0,z);
double a3x = result3[0];
double a3y = result3[1];
double a3z = result3[2];
double k3x = dt*a3x;
double k3y = dt*a3y;
double k3z = dt*a3z;

vector<double> result4 = acceleration(t+dt,vx+k3x,vy+k3y,vz+k3z,z);
double a4x = result4[0];
double a4y = result4[1];
double a4z = result4[2];
double k4x = dt*a4x;
double k4y = dt*a4y;
double k4z = dt*a4z;



vx += (1.0/6.0)*(k1x+2.0*k2x+2.0*k3x+k4x);
vy += (1.0/6.0)*(k1y+2.0*k2y+2.0*k3y+k4y);
vz += (1.0/6.0)*(k1z+2.0*k2z+2.0*k3z+k4z);

x += vx * dt;
y += vy * dt;///total distance calculated
z += vz * dt;

vector<double> bruhVector;
bruhVector.push_back(x);
bruhVector.push_back(y);
bruhVector.push_back(z);
bruhVector.push_back(vx);
bruhVector.push_back(vy);
bruhVector.push_back(vz);
return bruhVector;

}

[lastchance]

@B3anss,
Please put your code in code tags so that we can see its structure, and also so that we can run it in CPP shell without having to download it.

At the moment, for each timestep ... you are calling the Runge-Kutta update (routine distance_final) SIX times! You should just call it once per timestep.

Moreover, you are not updating your generalised coordinates (x,y,z,vx,vy,vz) at all - you are just putting them in an unused array bruhVector[].

To be honest, it would be simpler just to pass x,y,z,vx,vy,vz as reference parameters; then they will be updated in the Runge-Kutta routine and returned to the calling procedure. You don't need to push_back and then return a 6-element vector.


For the long term, if you used a valarray to hold your dependent variables (positions and velocities), then your Runge-Kutta routine would be only one sixth of the size. Just saying.

[Niccolo]

Please edit your post and put your code between code tags.

Are you using a debugger, and can you trace through the code using one?

Any idea at all what is happening and where things start to go wrong?

edit:

Like lastchance, I see a great deal of confusion here.

I can't really deduce what the plan is (there are few comments and hardly more than a "it's broke, what's wrong" question here).

For you to see what's wrong, I re-iterate the question: Do you have a debugger, and have you traced through to see what is happening?

This may be your best option at the moment. You'll see things you don't expect happen (which we can't see from here since we can't read your mind).

In order to be of real help (despite the fact you may get lots of helpful observations, none will likely fix the problem(s) ), we need to better recognize what the plan really is here.

Why, for example, do you have x_path_RK4 (and another for y and then z), as opposed to one container of a 3d vector object?

[cpp2019]

Hi,

In your code above, what does the "3" stand for?

for (int t = 1; t < 1000; t++)3
{
...

I'm currently learning C++. And may post questions here in the future.

[ne555]

$ g++ foo.cpp -ggdb
$ gdb ./a.out
(gdb) break 178
Breakpoint 1 at 0x2edf: file foo.cpp, line 178.
(gdb) run
Breakpoint 1, acceleration (t=1, vx=0, vy=127.04705941016442, vz=73.348038114391002, z=3) at foo.cpp:178
178		return myArray;
(gdb) backtrace
#0  acceleration (t=1, vx=0, vy=127.04705941016442, vz=73.348038114391002, z=3) at foo.cpp:178
#1  0x000055555555626e in distance_final (x=0, y=2, z=3, vx=0, vy=127.04705941016442, vz=73.348038114391002, 
    t=1) at foo.cpp:185
#2  0x0000555555555dc3 in main () at foo.cpp:122
(gdb) info args
t = 1
vx = 0
vy = 127.04705941016442
vz = 73.348038114391002
z = 3
(gdb) info locals
v = 146.70000000000002
vw = 1.3258300914694328e+19
vyw_sim = 0
S = 2.7533124740282869e-18
Cd = 0.41049999999999998
Cl = 6.8832811850707175e-18
vxw_sim = 0
adragx = -0
adragy = -7.306416212882902e+18
adragz = -4.2182109318404562e+18
w = 2845.197497765796
aMagx = 40.557926828860509
aMagy = 3.1744464756010721e+18
aMagz = 117.37129505436201
ax = 40.557926828860509
ay = -4.1319697372818299e+18
az = -4.2182109318404562e+18
myArray = std::vector of length 3, capacity 4 = {40.557926828860509, -4.1319697372818299e+18, 
  -4.2182109318404562e+18}
(gdb) quit

notice that some of your variables are quite big, in the order of 1e18, check out why
others are almost 0, ¿is that good?


by the way vw = sqrt((pow(2, (vx - vxw))) + (pow(2, (vy - vyw))) + (pow(2, vz)));
pow(5, 2) == 5*5
pow(2, 5) == 2*2*2*2*2
¿which one do you want?

[lastchance]

Oh, and while we look at it, you do realise that
pow(e, (-t/tau))
is more obviously coded as
exp(-t/tau)

There are quite a few gems in that code when we delve into it!
double az = adragz+aMagz-32.174;
I have this horrible feeling that you are working in imperial units!

Might be best if you started by setting the drag coefficient and Magnus-effect coefficient to zero, so that you can compare the results for basic projectile motion without the effects of friction or spin.