Final polish to the Naval Gun

And now I want to study that interpolation function, so I have reference to return to in case I forget (again) how I did it. What happens?

double ref_Mach_table[MACH_CD_ELEMENTS] =
{-0.10, 0.00, 0.40, 0.60, 0.87, 1.02, 1.15, 2.40, 3.40};

double ref_Cd_table[MACH_CD_ELEMENTS] =
{0.000, 0.231, 0.229, 0.172, 0.139, 0.333, 0.341, 0.268, 0.255};

//............

double Cd_Calculator(double Mach)
{
int index_count;

for(index_count = 0; index_count < MACH_CD_ELEMENTS; index_count++)
{
if(ref_Mach_table[index_count] > Mach)
break;
}

double inter_factor =
(Mach - ref_Mach_table[index_count - 1]) / (ref_Mach_table[index_count] - ref_Mach_table[index_count - 1]);

double output =
((ref_Cd_table[index_count] - ref_Cd_table[index_count - 1]) * inter_factor) + ref_Cd_table[index_count - 1];

return output;
}

The function is called using the calculated Mach number as the parameter. The for loop counts up (increasing the index_count) until the value in the Mach table (the ref_Mach_table array element of the index_count) exceeds the value of the passed Mach number, then breaks out with the corresponding index_count.

The range of the table values, the upper and the lower, that enclose the passed Mach value is then calculated;

ref_Mach_table[index_count] - ref_Mach_table[index_count - 1]

For example, a value of Mach 1.9 is passed. That being between 2.40 and 1.15, the formula will calculate the Mach range of;

2.40 - 1.15 = 1.25

Then the difference of the passed Mach value from the lower table value is calculated;

1.9 - 1.15 = 0.75

Then, in coefficient (percentage, if you like) of the relation of the relation of the true difference to the available table range is calculated;

0.75 / 1.25 = 0.6

And that is the interpolation factor for the ref_Cd_table array, and what is happening in the line that calculates inter_factor. I am just doing a straight forward linear interpolation here, by the way. Armed with inter_factor and index count, I then move to the ref_Cd_table and find the range of the two Cd values that correspond to the index_count and index_count - 1.

ref_Cd_table[index_count] - ref_Cd_table[index_count - 1]

These values are 0.268 and 0.341. So we get;

0.268 - 0.341 = -0.073

A negative value, as the Cd is reducing with increase in Mach at these velocities. I multiply that by the inter_factor to get what the Cd change is between the two Cd table values at 0.6 (60%) between them.

0.6 x -0.073 = -0.0438

I now add that to the Cd value that corresponds to the [index_count - 1] element, which is 0.341;

0.341 + -0.0438 = 0.2972

And that is the new interpolated Cd returned by the function (result of the computations for the variable output). This is the Cd value that gets passed to the drag calculation formula.

Now, finally, a refinement of the drag calculations by obtaining the air density at the shell's altitude. As air density reduces with increasing altitude, the shell is going to be traveling through a less dense medium at the top of its trajectory. I believe it was this ship, that ended her days as a war trophy and was subjected to nuclear bomb tests at Kwajalein atoll, that had high elevation guns (37º, where most maximum elevations were in the region of 25º to 30º), which greatly enhanced the fighting capabilities of the already formidable 8" main weapon by sending the shell high into the lower density atmosphere.

I got the density ratios for ISA out of a book I have, Aerodynamics for Naval Aviators, though I am sure it can be found also with a good internet search.

const int REF_RHO_ELEMENTS = 12;

double ref_rhos[REF_RHO_ELEMENTS] =
{1.000, 0.8617, 0.7385, 0.6292, 0.5328, 0.4481, 0.3741, 0.3099, 0.2462, 0.1936, 0.1522, 0.1197};
// These are percentage values of the air density, relative to the SL density of 1.225 kg/m³
// The resolution of the density ratio is for every 1525 meters.

// Related declarations...
double rho_gradient = 1525; // Meters
double actual_rho;

double Rho_Calculator(double altitude, double gradient, double sl_rho)
{
double alt_ref = altitude / gradient;
// Eg; 1872 / 1525 = 1.227541

int index_count = int(alt_ref);
// Eg; int(1.227541) = 1. This is the index position in ref_rhos[]
// As this is a constant increment of altitude gradient,
// the interpolation factor is pretty direct, as in;

double inter_factor = alt_ref - index_count;

double rho_coeff =
((ref_rhos[index_count + 1] - ref_rhos[index_count]) * inter_factor) + ref_rhos[index_count];

return sl_rho * rho_coeff;
}

This is actually a simpler interpolation than the previous, as it does not have a variable primary interpolation range, as the first one did. The altitude increments are constant, and only the density ratio lapse relaxes with increasing altitude. It is fairly self explanatory. It is called, BEFORE the drag calculation, like this...

actual_rho = Rho_Calculator(ypos, rho_gradient, ISA_Rho);

And the drag calculation, to recap, looks like this...

double Velocity_Calculator(double vel, double mass, double rho, double d_coeff, double FA)
{
double energy = vel * mass;
double drag = 0.5 * rho * pow(vel, 2) * d_coeff * FA;
energy -= drag;
return energy / mass;
}

Note that I now have it as a separate function, called from inside the update loop, after all the Rho and Mach calculations are done, like this...

velocity =
Velocity_Calculator(velocity, mass_kg, actual_rho, Cd, frontal_area);

That would be all, for now.