diff --git a/CMakeLists.txt b/CMakeLists.txt index d86bae1..8eb16d1 100644 --- a/CMakeLists.txt +++ b/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" ) diff --git a/include/Astro/dragAccelerationModel.hpp b/include/Astro/dragAccelerationModel.hpp new file mode 100755 index 0000000..4f783f8 --- /dev/null +++ b/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 + +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 diff --git a/include/Astro/eddyCurrentModel.hpp b/include/Astro/eddyCurrentModel.hpp new file mode 100755 index 0000000..e43445e --- /dev/null +++ b/include/Astro/eddyCurrentModel.hpp @@ -0,0 +1,45 @@ +/* + * Copyright (c) + */ + +#ifndef ASTRO_EDDY_CURRENT_HPP +#define ASTRO_EDDY_CURRENT_HPP + +#include + +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 \ No newline at end of file diff --git a/test/testDragAccelerationModel.cpp b/test/testDragAccelerationModel.cpp new file mode 100755 index 0000000..eba6019 --- /dev/null +++ b/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 +#include + +#include + +#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 diff --git a/test/testEddyTorqueModel.cpp b/test/testEddyTorqueModel.cpp new file mode 100755 index 0000000..a1f3954 --- /dev/null +++ b/test/testEddyTorqueModel.cpp @@ -0,0 +1,91 @@ +/* + * Copyright (c) + */ +#include + +#include + +#include + +#include + +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