Air drag model and test

CMakeLists.txt

@@ -44,6 +44,8 @@ CMAKE_DEPENDENT_OPTION(BUILD_COVERAGE_ANALYSIS "Build code coverage analysis"
 set(TEST_SRC
   "${TEST_SRC_PATH}/testAstro.cpp"
   "${TEST_SRC_PATH}/testConstants.cpp"
+  "${TEST_SRC_PATH}/testDragAccelerationModel.cpp"
+  "${TEST_SRC_PATH}/testEddyTorqueModel.cpp"
   "${TEST_SRC_PATH}/testOrbitalElementConversions.cpp"
   "${TEST_SRC_PATH}/testTwoBodyMethods.cpp"
 )

include/Astro/dragAccelerationModel.hpp

@@ -0,0 +1,69 @@
+/*
+ * Copyright (c) 2014-2015 Kartik Kumar (me@kartikkumar.com)
+ * Distributed under the MIT License.
+ * See accompanying file LICENSE.md or copy at http://opensource.org/licenses/MIT
+ */
+
+#ifndef ASTRO_DRAG_ACCELERATION_MODEL_HPP
+#define ASTRO_DRAG_ACCELERATION_MODEL_HPP
+
+#include <SML/linearAlgebra.hpp>
+
+namespace astro
+{
+
+//! Compute drag acceleration on a cannonball.
+/*!
+ * Computes drag acceleration using a cannonball model. The model for the acceleration is given by
+ *
+ * \f[
+ *      a_{drag} = \frac{1}{2} \frac{C_{d}}{m} \rho S V \vec{V}
+ * \f]
+ * 
+ * where \f$a_{drag}\f$ is the drag acceleration, \f$C_{d}\f$ is the drag coefficient, \f$m\f$ is
+ * the mass, \f$\rho\f$ is the atmospheric density, \f$\vec{V}\f$ is the relative velocity with 
+ * respect to the body-fixed frame and \f$S\f$ is the drag area, that is, the projected area
+ * of the object perpendicular to \f$\vec{V}\f$.
+ *
+ * @param[in]  dragCoefficient      Drag coefficient               [adim]
+ * @param[in]  atmosphericDensity   Atmospheric density            [kg m^-3]
+ * @param[in]  velocity             Velocity vector (3x1)          [m s^-1]
+ * @param[in]  dragArea             Drag area                      [m^2]
+ * @param[in]  mass                 Mass                           [kg]
+ * @return                          Drag acceleration vector (3x1) [m s^-2]                                                      
+ */
+template< typename Real, typename Vector3 >
+Vector3 computeDragAcceleration( const Real     dragCoefficient, 
+                                 const Real     atmosphericDensity, 
+                                 const Vector3& velocity,
+                                 const Real     dragArea,
+                                 const Real     mass ); 
+
+//! Compute drag acceleration on a cannonball.
+template< typename Real, typename Vector3 >
+Vector3 computeDragAcceleration( const Real     dragCoefficient, 
+                                 const Real     atmosphericDensity, 
+                                 const Vector3& velocity,
+                                 const Real     dragArea,
+                                 const Real     mass )
+{
+    Vector3 dragAcceleration = velocity;
+
+    // Compute the squared norm of the velocity.
+    const Real normVelocity = sml::norm< Real >( velocity );
+
+    // Compute a premultiplier so that it does not have to be written several times.
+    const Real preMultiplier = 0.5 * dragCoefficient * atmosphericDensity 
+                               * dragArea * normVelocity / mass;
+
+    // Compute the drag acceleration vector.
+    dragAcceleration[ 0 ] = preMultiplier * velocity[ 0 ];
+    dragAcceleration[ 1 ] = preMultiplier * velocity[ 1 ];
+    dragAcceleration[ 2 ] = preMultiplier * velocity[ 2 ];
+
+    return dragAcceleration;
+}
+
+} // namespace astro
+
+#endif // ASTRO_DRAG_ACCELERATION_MODEL_HPP

include/Astro/eddyCurrentModel.hpp

@@ -0,0 +1,45 @@
+/*
+ * Copyright (c) 
+ */
+
+#ifndef ASTRO_EDDY_CURRENT_HPP
+#define ASTRO_EDDY_CURRENT_HPP
+
+#include <SML/linearAlgebra.hpp>
+
+namespace astro
+{
+
+//! Compute eddy current torque.
+/*!
+ * Computes the eddy current torque on a certain object generated by an external source. The model
+ * for the torque is given by
+ *
+ * \f[
+ *      \vec{T}_{eddy} = \vec{m} \times \vec{B}
+ * \f]
+ * 
+ * where \f$\vec{m}\f$ is the magnetic moment of the object and \f$\vec{B}\f$ is the magnetic field
+ * generated by an external source.
+ *      
+ * @param  magneticMoment    Magnetic Moment         [A * m^2]
+ * @param  magneticField     Magnetic Field          [T]
+ * @return                   Eddy Current Torque     [N * m]                     
+ */
+ template< typename Vector3 >
+ Vector3 computeEddyTorque( const Vector3&  magneticMoment, 
+                            const Vector3&  magneticField ); 
+
+  //! Compute eddy current torque.
+ template< typename Vector3 >
+  Vector3 computeEddyTorque( const Vector3&  magneticMoment, 
+                             const Vector3&  magneticField )
+
+ {
+    // Compute the eddy current torque.
+    return sml::cross< Vector3 >( magneticMoment, magneticField);
+ }
+
+ } // namespace astro
+
+ #endif // ASTRO_EDDY_CURRENT_HPP

test/testDragAccelerationModel.cpp

@@ -0,0 +1,109 @@
+/*
+ * Copyright (c) 2015 Kartik Kumar (me@kartikkumar.com)
+ * Distributed under the MIT License.
+ * See accompanying file LICENSE.md or copy at http://opensource.org/licenses/MIT
+ */
+
+#include <cmath>
+#include <vector>
+
+#include <catch.hpp>
+
+#include "Astro/dragAccelerationModel.hpp"
+
+namespace astro
+{
+namespace tests
+{
+
+typedef double Real;
+typedef std::vector< Real > Vector;
+
+TEST_CASE( "Obtain drag acceleration: test 1", "[drag, acceleration, models]" )
+{
+    // Set expected drag acceleration vector [m/s].
+    Vector expectedDragAcceleration( 3 );
+    expectedDragAcceleration[ 0 ] = 0.107800109999944e-4;
+    expectedDragAcceleration[ 1 ] = 0.0;
+    expectedDragAcceleration[ 2 ] = 0.000154000157143e-4;
+
+    // Set epsilon = error between expected value and computed value.
+    const Real epsilon = 1.0e-10;
+
+    // Set drag coefficient.
+    const Real dragCoefficient = 2.2;
+
+    // Set atmospheric density [kg/m^3].
+    const Real atmosphericDensity = 2.0e-11;
+
+    // Velocity vector in [m/s].
+    Vector velocity( 3 );
+    velocity[ 0 ] = 7000.0;
+    velocity[ 1 ] = 0.0;
+    velocity[ 2 ] = 10.0;
+
+    // Set drag area [m^2].
+    const Real dragArea = 5.0;
+
+    // Set mass [kg].
+    const Real mass = 500.0;
+
+     //! Compute drag acceleration.
+    const Vector dragAcceleration = computeDragAcceleration( dragCoefficient, 
+                                                             atmosphericDensity, 
+                                                             velocity,
+                                                             dragArea,
+                                                             mass ); 
+
+    // Check if computed mean motion matches expected value.
+    REQUIRE( std::fabs( dragAcceleration[ 0 ] - expectedDragAcceleration[ 0 ] ) <= epsilon );
+    REQUIRE( std::fabs( dragAcceleration[ 1 ] - expectedDragAcceleration[ 1 ] ) <= epsilon );
+    REQUIRE( std::fabs( dragAcceleration[ 2 ] - expectedDragAcceleration[ 2 ] ) <= epsilon );
+}
+
+TEST_CASE( "Obtain drag acceleration: test 2", "[obtain-drag-acceleration-2]" )
+{
+    // Reference: http://tudat.tudelft.nl/
+
+    // Set expected drag acceleration vector [m/s].
+    Vector expectedDragAcceleration( 3 );
+    expectedDragAcceleration[ 0 ] = 0.0;
+    expectedDragAcceleration[ 1 ] = 0.0;
+    expectedDragAcceleration[ 2 ] =  267.4211815284975;
+
+    // Set epsilon = error between expected value and computed value.
+    const Real epsilon = 1.0e-10;
+
+    // Set drag coefficient.
+    const Real dragCoefficient = 1.1;
+
+    // Set atmospheric density [kg/m3].
+    const Real atmosphericDensity = 3.5e-5;
+
+    // Velocity vector in [m/s].
+    Vector velocity( 3 );
+    velocity[ 0 ] = 0.0;
+    velocity[ 1 ] = 0.0;
+    velocity[ 2 ] = 3491.0;
+
+    // Set drag area [m^2]
+    const Real dragArea = 2.2;
+
+    // Set mass [kg]
+    const Real mass = 1.93;
+
+     //! Compute drag acceleration.
+    const Vector dragAcceleration = computeDragAcceleration( dragCoefficient, 
+                                                             atmosphericDensity, 
+                                                             velocity,
+                                                             dragArea,
+                                                             mass ); 
+
+    // Check if computed mean motion matches expected value.
+    REQUIRE( std::fabs(dragAcceleration[ 0 ] - expectedDragAcceleration[ 0 ]) <= epsilon );
+    REQUIRE( std::fabs(dragAcceleration[ 1 ] - expectedDragAcceleration[ 1 ]) <= epsilon );
+    REQUIRE( std::fabs(dragAcceleration[ 2 ] - expectedDragAcceleration[ 2 ]) <= epsilon );
+}
+
+} // namespace tests
+} // namespace astro

test/testEddyTorqueModel.cpp

@@ -0,0 +1,91 @@
+/*
+ * Copyright (c) 
+ */
+#include <cmath>
+
+#include <vector>
+
+#include <catch.hpp>
+
+#include <Astro/eddyCurrentModel.hpp>
+
+using namespace std;
+
+namespace astro
+{
+
+namespace tests
+{
+
+typedef double Real;
+typedef std::vector< Real > Vector3;
+
+TEST_CASE( "Obtain eddy current torque: test 1", "[obtain-eddy-torque-1]" )
+{
+    // Set expected eddy current torque vector [N * m].
+    Vector3 expectedEddyTorque( 3 );
+    expectedEddyTorque[ 0 ] =  0.095000000000000;
+    expectedEddyTorque[ 1 ] =  0.065000000000000;
+    expectedEddyTorque[ 2 ] = -0.149000000000000;
+
+    // Set magnetic moment vector [A * m^2].
+    Vector3 magneticMoment( 3 );
+    magneticMoment[ 0 ] = 100.0;
+    magneticMoment[ 1 ] = 1000.0;
+    magneticMoment[ 2 ] = 500.0;
+
+    // Set magnetic field vector [T].
+    Vector3 magneticField( 3 );
+    magneticField[ 0 ] = 150e-6;
+    magneticField[ 1 ] = 10e-6;
+    magneticField[ 2 ] = 100e-6;
+
+    // Set epsilon = error between expected value and computed value.
+    const Real epsilon = 1.0e-10;
+
+     //! Compute eddy current torque.
+    const Vector3 eddyTorque = computeEddyTorque( magneticMoment, 
+                                                  magneticField ); 
+
+    // Check if computed torque matches expected value.
+    REQUIRE( std::fabs(eddyTorque[ 0 ] - expectedEddyTorque[ 0 ]) <= epsilon );
+    REQUIRE( std::fabs(eddyTorque[ 1 ] - expectedEddyTorque[ 1 ]) <= epsilon );
+    REQUIRE( std::fabs(eddyTorque[ 2 ] - expectedEddyTorque[ 2 ]) <= epsilon );
+}
+
+TEST_CASE( "Obtain eddy current torque: test 2", "[obtain-eddy-torque-2]" )
+{
+    // Set expected eddy current torque vector [N * m].
+    Vector3 expectedEddyTorque( 3 );
+    expectedEddyTorque[ 0 ] = 0.0;
+    expectedEddyTorque[ 1 ] = 0.0;
+    expectedEddyTorque[ 2 ] = 0.0;
+
+    // Set magnetic moment vector [A * m^2].
+    Vector3 magneticMoment( 3 );
+    magneticMoment[ 0 ] = 0.0;
+    magneticMoment[ 1 ] = 0.0;
+    magneticMoment[ 2 ] = 1150.0;
+
+    // Set magnetic field vector [T].
+    Vector3 magneticField( 3 );
+    magneticField[ 0 ] = 0.0;
+    magneticField[ 1 ] = 0.0;
+    magneticField[ 2 ] = 127e-6;
+
+    // Set epsilon = error between expected value and computed value.
+    const Real epsilon = 1.0e-10;
+
+     //! Compute eddy current torque.
+    const Vector3 eddyTorque = computeEddyTorque( magneticMoment, 
+                                                  magneticField ); 
+
+    // Check if computed torque matches expected value.
+    REQUIRE( std::fabs(eddyTorque[ 0 ] - expectedEddyTorque[ 0 ]) <= epsilon );
+    REQUIRE( std::fabs(eddyTorque[ 1 ] - expectedEddyTorque[ 1 ]) <= epsilon );
+    REQUIRE( std::fabs(eddyTorque[ 2 ] - expectedEddyTorque[ 2 ]) <= epsilon );
+}
+
+
+} // namespace tests
+} // namespace astro