Compute first order, orbit averaged rate of change in orbital elements due to j2 perturbation

ProjectFiles.cmake

@@ -6,8 +6,9 @@
 set(TEST_SRC
   "${TEST_SRC_PATH}/testAstro.cpp"
   "${TEST_SRC_PATH}/testConstants.cpp"
+  "${TEST_SRC_PATH}/testJ2Gravity.cpp"
   "${TEST_SRC_PATH}/testOrbitalElementConversions.cpp"
   "${TEST_SRC_PATH}/testStateVectorIndices.cpp"
   "${TEST_SRC_PATH}/testSolarRadiationPressureAccelerationModel.cpp"
   "${TEST_SRC_PATH}/testTwoBodyMethods.cpp"
-)
+)

include/Astro/astro.hpp

@@ -8,8 +8,10 @@
 #define ASTRO_HPP
 
 #include "Astro/constants.hpp"
+#include "Astro/j2Gravity.hpp"
 #include "Astro/orbitalElementConversions.hpp"
 #include "Astro/twoBodyMethods.hpp"
 #include "Astro/solarRadiationPressureAccelerationModel.hpp"
+#include "Astro/stateVectorIndices.hpp"
 
 #endif // ASTRO_HPP

include/Astro/j2Gravity.hpp

@@ -0,0 +1,71 @@
+/*
+ * Copyright (c) 2014-2016 Kartik Kumar, Dinamica Srl (me@kartikkumar.com)
+ * Copyright (c) 2014-2016 Abhishek Agrawal, Delft University of Technology
+ *                                          (abhishek.agrawal@protonmail.com)
+ * Distributed under the MIT License.
+ * See accompanying file LICENSE.md or copy at http://opensource.org/licenses/MIT
+ */
+
+#ifndef J2_GRAVITY_HPP
+#define J2_GRAVITY_HPP
+
+#include <cmath>
+
+#include <Astro/stateVectorIndices.hpp>
+
+namespace astro
+{
+
+//! Compute orbit-averaged rate of change in the orbital elements due to J2 perturbation.
+/*!
+ * Computes the first order, orbit-averaged rate of change in longitude of ascending node and
+ * argument of periapsis due to the J2 perturbation.
+ *
+ * @tparam      Real                Real type
+ * @tparam      Vector6             Vector of 6 elements
+ *
+ * @param[in]   keplerianElements   Vector containing Keplerian elements
+ *                                  N.B.: Order of elements and units is important!
+ *                                  keplerianElements( 0 ) = semiMajorAxis                [km]  <br>
+ *                                  keplerianElements( 1 ) = eccentricity                 [-]   <br>
+ *                                  keplerianElements( 2 ) = inclination                  [rad] <br>
+ *                                  keplerianElements( 3 ) = argument of periapsis        [rad] <br>
+ *                                  keplerianElements( 4 ) = longitude of ascending node  [rad] <br>
+ *                                  keplerianElements( 5 ) = true anomaly                 [rad]
+ * @param[in]  meanMotion           meanMotion of the orbiting object [deg/day]
+ * @param[in]  earthEquatorialRadius     Earth equatorial radius [km]
+ *
+ * @param[out] longitudeAscendingNodeDot Stores rate of change in longitude of ascending node
+ *                                       due to J2 perturbation [deg/day]
+ * @param[out] argumentOfPeriapsisDot    Stores rate of change in argument of periapsis
+ *                                       due to J2 perturbation [deg/day]
+ */
+ template< typename Real, typename Vector6 >
+ void computeFirstOrderAveragedEffectJ2Perturbation(
+    const Vector6& keplerianElements,
+    const Real     meanMotion,
+    const Real     earthEquatorialRadius,
+    Real&          longitudeAscendingNodeDot,
+    Real&          argumentOfPeriapsisDot )
+ {
+    const Real J2 = 0.00108263;
+    const Real semiMajorAxis = keplerianElements[ astro::semiMajorAxisIndex ];
+    const Real inclination   = keplerianElements[ astro::inclinationIndex ];
+    const Real eccentricity  = keplerianElements[ astro::eccentricityIndex ];
+
+    longitudeAscendingNodeDot
+        = -1.5 * meanMotion * J2
+          * ( earthEquatorialRadius / semiMajorAxis ) * ( earthEquatorialRadius / semiMajorAxis )
+          * std::cos( inclination )
+          / ( ( 1 - eccentricity * eccentricity ) * ( 1 - eccentricity * eccentricity ) );
+
+    argumentOfPeriapsisDot
+        = 0.75 * meanMotion * J2
+          * ( earthEquatorialRadius / semiMajorAxis ) * ( earthEquatorialRadius / semiMajorAxis )
+          * ( 4 - 5 * std::sin( inclination ) * std::sin( inclination ) )
+          / ( ( 1 - eccentricity * eccentricity ) * ( 1 - eccentricity * eccentricity ) );
+ }
+
+} // namespace_astro
+
+#endif // J2_GRAVITY_HPP

test/testJ2Gravity.cpp

@@ -0,0 +1,156 @@
+/*
+ * Copyright (c) 2014-2016 Kartik Kumar, Dinamica Srl (me@kartikkumar.com)
+ * Copyright (c) 2014-2016 Abhishek Agrawal, Delft University of Technology
+ *                                           (abhishek.agrawal@protonmail.com)
+ * Distributed under the MIT License.
+ * See accompanying file LICENSE.md or copy at http://opensource.org/licenses/MIT
+ */
+
+#include <limits>
+#include <vector>
+
+#include <catch.hpp>
+
+#include "Astro/j2Gravity.hpp"
+#include "Astro/twoBodyMethods.hpp"
+
+namespace astro
+{
+namespace tests
+{
+
+typedef double Real;
+typedef std::vector< Real > Vector;
+
+TEST_CASE( "First order, orbit-averaged Keplerian element rate of change due to J2 perturbation",
+           "[first-order-orbit-averaged-j2]" )
+{
+    const Real pi = 3.14159265358979323846;
+    // Earth's equatorial radius in [km].
+    const Real earthEquatorialRadius = 6378.13649;
+    // Earth's gravitational parameter [km^3 s^-2].
+    const Real earthGravitationalParameter = 398600.4418;
+
+    SECTION( "Test for Shuttle (LEO) orbit" )
+    {
+        Vector keplerianElements( 6 );
+        keplerianElements[ 0 ] = 6700.0;
+        keplerianElements[ 1 ] = 0.0;
+        keplerianElements[ 2 ] = 28.0 * pi / 180.0 ;
+        keplerianElements[ 3 ] = 0.0;
+        keplerianElements[ 4 ] = 0.0;
+        keplerianElements[ 5 ] = 0.0;
+
+        // Mean motion or orbiting object in radians per second.
+        const Real semiMajorAxis = keplerianElements[ 0 ];
+        const Real massOrbitingBody = 0.0;
+        const Real meanMotion = astro::computeKeplerMeanMotion( semiMajorAxis,
+                                                                earthGravitationalParameter,
+                                                                massOrbitingBody );
+        // Mean motion in degrees per day.
+        const Real meanMotionDegreesPerDay = ( meanMotion * 86400.0 ) * ( 180.0 / pi );
+
+        Real longitudeAscendingNodeDot = 0.0;
+        Real argumentOfPeriapsisDot = 0.0;
+
+        const Real expectedLongitudeAscendingNodeDot = -7.35;
+        const Real expectedArgumentOfPeriapsisDot = 12.05;
+
+        computeFirstOrderAveragedEffectJ2Perturbation(
+            keplerianElements,
+            meanMotionDegreesPerDay,
+            earthEquatorialRadius,
+            longitudeAscendingNodeDot,
+            argumentOfPeriapsisDot );
+
+        REQUIRE( expectedLongitudeAscendingNodeDot
+                    == Approx( longitudeAscendingNodeDot ).epsilon( 1.0e-2 ) );
+        REQUIRE( expectedArgumentOfPeriapsisDot
+                    == Approx( argumentOfPeriapsisDot ).epsilon( 1.0e-2 ) );
+    }
+
+    SECTION( "Test for GPS (HEO) orbit" )
+    {
+        Vector keplerianElements( 6 );
+        keplerianElements[ 0 ] = 26600.0;
+        keplerianElements[ 1 ] = 0.0;
+        keplerianElements[ 2 ] = 60.0 * pi / 180.0 ;
+        keplerianElements[ 3 ] = 0.0;
+        keplerianElements[ 4 ] = 0.0;
+        keplerianElements[ 5 ] = 0.0;
+
+        // Mean motion or orbiting object in radians per second.
+        const Real semiMajorAxis = keplerianElements[ 0 ];
+        const Real massOrbitingBody = 0.0;
+        const Real meanMotion = astro::computeKeplerMeanMotion( semiMajorAxis,
+                                                                earthGravitationalParameter,
+                                                                massOrbitingBody );
+        // Mean motion in degrees per day.
+        const Real meanMotionDegreesPerDay = ( meanMotion * 86400.0 ) * ( 180.0 / pi );
+
+        Real longitudeAscendingNodeDot = 0.0;
+        Real argumentOfPeriapsisDot = 0.0;
+
+        const Real expectedLongitudeAscendingNodeDot = -0.033;
+        const Real expectedArgumentOfPeriapsisDot = 0.008;
+
+        computeFirstOrderAveragedEffectJ2Perturbation(
+            keplerianElements,
+            meanMotionDegreesPerDay,
+            earthEquatorialRadius,
+            longitudeAscendingNodeDot,
+            argumentOfPeriapsisDot );
+
+        REQUIRE( expectedLongitudeAscendingNodeDot
+                    == Approx( longitudeAscendingNodeDot ).epsilon( 1.0e-3 ) );
+        REQUIRE( expectedArgumentOfPeriapsisDot
+                    == Approx( argumentOfPeriapsisDot ).epsilon( 1.0e-3 ) );
+    }
+
+    SECTION( "Test for Geo-Stationary (GEO) orbit" )
+    {
+        Vector keplerianElements( 6 );
+        keplerianElements[ 0 ] = 42160.0;
+        keplerianElements[ 1 ] = 0.0;
+        keplerianElements[ 2 ] = 0.0;
+        keplerianElements[ 3 ] = 0.0;
+        keplerianElements[ 4 ] = 0.0;
+        keplerianElements[ 5 ] = 0.0;
+
+        // Mean motion or orbiting object in radians per second.
+        const Real semiMajorAxis = keplerianElements[ 0 ];
+        const Real massOrbitingBody = 0.0;
+        const Real meanMotion = astro::computeKeplerMeanMotion( semiMajorAxis,
+                                                                earthGravitationalParameter,
+                                                                massOrbitingBody );
+        // Mean motion in degrees per day.
+        const Real meanMotionDegreesPerDay = ( meanMotion * 86400.0 ) * ( 180.0 / pi );
+
+        Real longitudeAscendingNodeDot = 0.0;
+        Real argumentOfPeriapsisDot = 0.0;
+
+        const Real expectedLongitudeAscendingNodeDot = -0.013;
+        const Real expectedArgumentOfPeriapsisDot = 0.026;
+
+        computeFirstOrderAveragedEffectJ2Perturbation(
+            keplerianElements,
+            meanMotionDegreesPerDay,
+            earthEquatorialRadius,
+            longitudeAscendingNodeDot,
+            argumentOfPeriapsisDot );
+
+        REQUIRE( expectedLongitudeAscendingNodeDot
+                    == Approx( longitudeAscendingNodeDot ).epsilon( 1.0e-3 ) );
+        REQUIRE( expectedArgumentOfPeriapsisDot
+                    == Approx( argumentOfPeriapsisDot ).epsilon( 1.0e-3 ) );
+    }
+}
+
+} // namespace tests
+} // namespace astro
+
+/*!
+ * References
+ * Wertz, J.R., et al. Space mission analysis and design, Third Edition,
+ * ISBN 1-881883-10-8, Space Technology Library.
+ */