Fix several bugs

source/calc.hpp

@@ -32,8 +32,7 @@ inline DegFT CalculateTemperatureAtAltitudeMcCoy(FeetT altitude,
 // https://wikipedia.org/wiki/Barometric_formula
 inline InHgT BarometricFormula(FeetT altitude, InHgT pressure,
                                DegFT temperature) {
-  const double kGasConstant = 1716.49;     // ft-lb / slug^{-1}R^{-1}
-  const double kMolarMassOfAir = 28.9644;  // lb/lb-mol
+  const double kGasConstant = 1716.49;  // ft-lb / slug^{-1}R^{-1}
   const FeetT kHeight = std::min(altitude, FeetT(kIsaTropopauseAltitudeFt));
 
   const double kExponent =
@@ -46,12 +45,9 @@ inline InHgT BarometricFormula(FeetT altitude, InHgT pressure,
 
   if (altitude > FeetT(kIsaTropopauseAltitudeFt)) {
     const double kNumberator = -1.0 * kStandardGravityFtPerSecSq *
-                               kMolarMassOfAir *
                                (altitude - kIsaTropopauseAltitudeFt).Value();
-
     const double kDemoninator =
         kGasConstant * DegRT(DegFT(kIsaMinimumTempDegF)).Value();
-
     output *= std::exp(kNumberator / kDemoninator);
   }
 

source/lob_builder.cpp

@@ -68,15 +68,8 @@ Builder::~Builder() { pimpl_->~Impl(); }
 Builder::Builder(const Builder& other)
     : pimpl_(new(buffer_.data()) Impl(*other.pimpl_)) {}
 
-Builder::Builder(Builder&& other) noexcept {
-  if (this != &other) {
-    if (pimpl_ != nullptr) {
-      pimpl_->~Impl();
-    }
-    pimpl_ = new (buffer_.data()) Impl();
-    std::swap(pimpl_, other.pimpl_);
-  }
-}
+Builder::Builder(Builder&& other) noexcept
+    : pimpl_(new(buffer_.data()) Impl(std::move(*other.pimpl_))) {}
 
 Builder& Builder::operator=(const Builder& rhs) {
   if (this != &rhs) {
@@ -93,8 +86,7 @@ Builder& Builder::operator=(Builder&& rhs) noexcept {
     if (pimpl_ != nullptr) {
       pimpl_->~Impl();
     }
-    pimpl_ = new (buffer_.data()) Impl();
-    std::swap(pimpl_, rhs.pimpl_);
+    pimpl_ = new (buffer_.data()) Impl(std::move(*rhs.pimpl_));
   }
   return *this;
 }
@@ -141,10 +133,6 @@ Builder& Builder::BCAtmosphere(AtmosphereReferenceT type) {
 
 Builder& Builder::BCDragFunction(DragFunctionT type) {
   switch (type) {
-    case DragFunctionT::kG1: {
-      pimpl_->pdrag_lut = &kG1Drags;
-      break;
-    }
     case DragFunctionT::kG2: {
       pimpl_->pdrag_lut = &kG2Drags;
       break;
@@ -165,7 +153,9 @@ Builder& Builder::BCDragFunction(DragFunctionT type) {
       pimpl_->pdrag_lut = &kG8Drags;
       break;
     }
+    case DragFunctionT::kG1:
     default: {
+      pimpl_->pdrag_lut = &kG1Drags;
       break;
     }
   }
@@ -549,6 +539,40 @@ void BuildStability(Impl* pimpl) {
                  FpsT(pimpl->build.velocity));
 }
 
+void BuildCoriolis(Impl* pimpl) {
+  assert(pimpl != nullptr);
+
+  if (!std::isnan(pimpl->azimuth_rad) && !std::isnan(pimpl->latitude_rad)) {
+    const DegreesT kAzimuthLimit(kDegreesPerTurn);
+    if (pimpl->azimuth_rad > kAzimuthLimit ||
+        pimpl->azimuth_rad < kAzimuthLimit * -1) {
+      pimpl->build.error = ErrorT::kAzimuthOOR;
+      return;
+    }
+    const DegreesT kLatitudeLimit(90);
+    if (pimpl->latitude_rad > kLatitudeLimit ||
+        pimpl->latitude_rad < kLatitudeLimit * -1) {
+      pimpl->build.error = ErrorT::kLatitudeOOR;
+      return;
+    }
+    // Coriolis Effect Page 178 of Modern Exterior Ballistics - McCoy
+    const double kCosL = std::cos(pimpl->latitude_rad).Value();
+    const double kSinA = std::sin(pimpl->azimuth_rad).Value();
+    const double kSinL = std::sin(pimpl->latitude_rad).Value();
+    const double kCosA = std::cos(pimpl->azimuth_rad).Value();
+
+    pimpl->build.corilolis.cos_l_sin_a =
+        2 * kAngularVelocityOfEarthRadPerSec * kCosL * kSinA;
+    pimpl->build.corilolis.sin_l = 2 * kAngularVelocityOfEarthRadPerSec * kSinL;
+    pimpl->build.corilolis.cos_l_cos_a =
+        2 * kAngularVelocityOfEarthRadPerSec * kCosL * kCosA;
+  } else {
+    pimpl->build.corilolis.cos_l_sin_a = 0;
+    pimpl->build.corilolis.sin_l = 0;
+    pimpl->build.corilolis.cos_l_cos_a = 0;
+  }
+}
+
 void BuildBoatright(Impl* pimpl) {
   assert(pimpl != nullptr);
   assert(pimpl->pdrag_lut != nullptr);
@@ -643,7 +667,6 @@ void BuildBoatright(Impl* pimpl) {
   SecT t(0.0);
 
   static const FpsT kTransonicBarrier(MachT(1.2), kSos);
-
   while (s.V().X() > kTransonicBarrier) {
     assert(t < SecT(60) && "This is taking too long");
     SolveStep(&s, &t, pimpl->build);
@@ -701,40 +724,6 @@ void BuildLitzAerodynamicJump(Impl* pimpl) {
   }
 }
 
-void BuildCoriolis(Impl* pimpl) {
-  assert(pimpl != nullptr);
-
-  if (!std::isnan(pimpl->azimuth_rad) && !std::isnan(pimpl->latitude_rad)) {
-    const DegreesT kAzimuthLimit(kDegreesPerTurn);
-    if (pimpl->azimuth_rad > kAzimuthLimit ||
-        pimpl->azimuth_rad < kAzimuthLimit * -1) {
-      pimpl->build.error = ErrorT::kAzimuthOOR;
-      return;
-    }
-    const DegreesT kLatitudeLimit(90);
-    if (pimpl->latitude_rad > kLatitudeLimit ||
-        pimpl->latitude_rad < kLatitudeLimit * -1) {
-      pimpl->build.error = ErrorT::kLatitudeOOR;
-      return;
-    }
-    // Coriolis Effect Page 178 of Modern Exterior Ballistics - McCoy
-    const double kCosL = std::cos(pimpl->latitude_rad).Value();
-    const double kSinA = std::sin(pimpl->azimuth_rad).Value();
-    const double kSinL = std::sin(pimpl->latitude_rad).Value();
-    const double kCosA = std::cos(pimpl->azimuth_rad).Value();
-
-    pimpl->build.corilolis.cos_l_sin_a =
-        2 * kAngularVelocityOfEarthRadPerSec * kCosL * kSinA;
-    pimpl->build.corilolis.sin_l = 2 * kAngularVelocityOfEarthRadPerSec * kSinL;
-    pimpl->build.corilolis.cos_l_cos_a =
-        2 * kAngularVelocityOfEarthRadPerSec * kCosL * kCosA;
-  } else {
-    pimpl->build.corilolis.cos_l_sin_a = 0;
-    pimpl->build.corilolis.sin_l = 0;
-    pimpl->build.corilolis.cos_l_cos_a = 0;
-  }
-}
-
 void BuildZeroAngle(Impl* pimpl) {
   assert(pimpl != nullptr);
 
@@ -839,15 +828,15 @@ Input Builder::Build() {
   if (pimpl_->build.error != ErrorT::kNotFormed) {
     return pimpl_->build;
   }
-  BuildBoatright(pimpl_);
+  BuildCoriolis(pimpl_);
   if (pimpl_->build.error != ErrorT::kNotFormed) {
     return pimpl_->build;
   }
-  BuildLitzAerodynamicJump(pimpl_);
-  BuildCoriolis(pimpl_);
+  BuildBoatright(pimpl_);
   if (pimpl_->build.error != ErrorT::kNotFormed) {
     return pimpl_->build;
   }
+  BuildLitzAerodynamicJump(pimpl_);
   BuildZeroAngle(pimpl_);
   if (pimpl_->build.error != ErrorT::kNotFormed) {
     return pimpl_->build;

source/lob_solve.cpp

@@ -44,15 +44,15 @@ Output LerpOutput(const TrajectoryStateT& s_now, const SecT t_now,
 void ApplyGyroscopicSpinDrift(const Input& in, Output* pouts, size_t size) {
   assert(pouts != nullptr);
   // If we can apply Boatright-Ruiz spin drift, prefer it
-  if (in.spindrift_factor > 0) {
+  if (!std::isnan(in.spindrift_factor)) {
     for (size_t i = 0; i < size; i++) {
       pouts[i].deflection +=
           in.spindrift_factor * std::fabs(pouts[i].elevation);
     }
     return;
   }
   // If we can apply Litz spin drift, use that
-  if (std::fabs(in.stability_factor) > 0) {
+  if (std::fabs(in.stability_factor) > 0.0) {
     for (size_t i = 0; i < size; i++) {
       pouts[i].deflection +=
           litz::CalculateGyroscopicSpinDrift(in.stability_factor,

test/source/calc_test.cpp

@@ -17,7 +17,7 @@ namespace tests {
 
 TEST(CalcTests, CalculateTemperatureAtAltitude) {
   // Test data from page 167 of Modern Exterior Ballistics - McCoy
-  const std::vector<uint16_t> kAltitudesFt = {
+  const std::vector<uint32_t> kAltitudesFt = {
       0,    500,  1000, 1500,  2000,  3000,  4000,  5000,  6000,
       7000, 8000, 9000, 10000, 15000, 20000, 25000, 30000, 35000};
   const std::vector<double> kExpectedResultsDegF = {
@@ -37,7 +37,7 @@ TEST(CalcTests, CalculateTemperatureAtAltitude) {
 
 TEST(CalcTests, CalculateTemperatureAtAltitudeMcCoy) {
   // Test data from page 167 of Modern Exterior Ballistics - McCoy
-  const std::vector<uint16_t> kAltitudesFt = {
+  const std::vector<uint32_t> kAltitudesFt = {
       0,    500,  1000, 1500,  2000,  3000,  4000,  5000,  6000,
       7000, 8000, 9000, 10000, 15000, 20000, 25000, 30000, 35000};
   const std::vector<double> kExpectedResultsDegF = {
@@ -59,7 +59,7 @@ TEST(CalcTests, CalculateTemperatureAtAltitudeMcCoy) {
 
 TEST(CalcTests, BarometricFormula) {
   // Test data from page 167 of Modern Exterior Ballistics - McCoy
-  const std::vector<uint16_t> kAltitudesFt = {
+  const std::vector<uint32_t> kAltitudesFt = {
       0,    500,  1000, 1500,  2000,  3000,  4000,  5000,  6000,
       7000, 8000, 9000, 10000, 15000, 20000, 25000, 30000, 35000};
   const std::vector<double> kExpectedResultsInHg = {
@@ -77,6 +77,27 @@ TEST(CalcTests, BarometricFormula) {
   }
 }
 
+TEST(CalcTests, BarometricFormulaStratosphere) {
+  // test data from
+  // https://www.sensorsone.com/altitude-pressure-units-conversion/
+  const std::vector<uint32_t> kAltitudesFt = {35000, 40000, 45000, 50000, 55000,
+                                              60000, 65000, 70000, 75000, 80000,
+                                              85000, 90000, 95000, 100000};
+  const std::vector<double> kExpectedResultsInHg = {
+      7.04, 5.54, 4.36, 3.42, 2.69, 2.12, 1.67,
+      1.31, 1.03, 0.82, 0.64, 0.51, 0.41, 0.32};
+  const double kError = 0.025;
+  for (uint32_t i = 0; i < kAltitudesFt.size(); i++) {
+    EXPECT_NEAR(
+        kExpectedResultsInHg.at(i),
+        lob::BarometricFormula(lob::FeetT(kAltitudesFt.at(i)),
+                               lob::InHgT(lob::kIsaSeaLevelPressureInHg),
+                               lob::DegFT(lob::kIsaSeaLevelDegF))
+            .Value(),
+        kError);
+  }
+}
+
 TEST(CalcTests, BarometricFormulaNegative) {
   constexpr int16_t kAltitude = -1000;
   constexpr double kExpectedResult = 31.02;
@@ -90,7 +111,7 @@ TEST(CalcTests, BarometricFormulaNegative) {
 
 TEST(CalcTests, CalculateAirDensityAtAltitude) {
   // Test data from page 167 of Modern Exterior Ballistics - McCoy
-  const std::vector<uint16_t> kAltitudesFt = {
+  const std::vector<uint32_t> kAltitudesFt = {
       0,    500,  1000, 1500,  2000,  3000,  4000,  5000,  6000,
       7000, 8000, 9000, 10000, 15000, 20000, 25000, 30000, 35000};
   const double kP0 = lob::kIsaSeaLevelAirDensityLbsPerCuFt;

test/source/lob_api_test.cpp

@@ -191,6 +191,12 @@ TEST(LobAPITest, InchToDeg) {
   EXPECT_DOUBLE_EQ(kA.Value(), lob::IphyT(kB).Value());
 }
 
+TEST(LobAPITest, InchToDegZero) {
+  const auto kA = lob::IphyT(5);
+  const auto kB = lob::DegreesT(lob::InchToDeg(kA.Value(), 0.0));
+  EXPECT_DOUBLE_EQ(0.0, lob::IphyT(kB).Value());
+}
+
 TEST(LobAPITest, JToFtLbs) {
   const auto kA = lob::JouleT(100);
   const auto kB = lob::FtLbsT(lob::JToFtLbs(kA.Value()));

test/source/lob_builder_test.cpp

@@ -75,7 +75,6 @@ TEST_F(BuilderTestFixture, MoveConstructor) {
   EXPECT_EQ(kVal.velocity, kTestMuzzleVelocity);
 }
 
-// Test for copy assignment operator
 TEST_F(BuilderTestFixture, CopyAssignmentOperator) {
   const double kTestBC = 0.425;
   const double kTestDiameter = 0.308;
@@ -95,7 +94,6 @@ TEST_F(BuilderTestFixture, CopyAssignmentOperator) {
   EXPECT_DOUBLE_EQ(kVal2.velocity, kTestMuzzleVelocity);
 }
 
-// Test for move assignment operator
 TEST_F(BuilderTestFixture, MoveAssignmentOperator) {
   const double kTestBC = 0.425;
   const double kTestDiameter = 0.308;
@@ -156,6 +154,21 @@ TEST_F(BuilderTestFixture, BuildMissingZeroInput) {
   EXPECT_EQ(kResult.error, lob::ErrorT::kZeroDataRequired);
 }
 
+TEST_F(BuilderTestFixture, InvalidDragFunctionIsG1) {
+  const double kTestBC = 1.0;
+  const uint16_t kTestMuzzleVelocity = 2500U;
+  const double kTestZeroAngle = 5.59;
+  const auto kInvalidDragFunction = static_cast<lob::DragFunctionT>(0xFF);
+  const lob::Input kResult = puut->BallisticCoefficientPsi(kTestBC)
+                                 .BCDragFunction(kInvalidDragFunction)
+                                 .InitialVelocityFps(kTestMuzzleVelocity)
+                                 .ZeroAngleMOA(kTestZeroAngle)
+                                 .Build();
+  for (size_t i = 0; i < lob::kG1Drags.size(); i++) {
+    EXPECT_EQ(kResult.drags.at(i), lob::kG1Drags.at(i));
+  }
+}
+
 TEST_F(BuilderTestFixture, BuildG1UsingCustomTable) {
   const double kTestBC = 1.0;
   const uint16_t kTestMuzzleVelocity = 2500U;