Aleixpinardell tabulated aerodynamic coefficients

Tudat/Astrodynamics/Aerodynamics/CMakeLists.txt

@@ -104,6 +104,10 @@ add_executable(test_TabulatedAtmosphere "${SRCROOT}${AERODYNAMICSDIR}/UnitTests/
 setup_custom_test_program(test_TabulatedAtmosphere "${SRCROOT}${AERODYNAMICSDIR}")
 target_link_libraries(test_TabulatedAtmosphere tudat_aerodynamics tudat_interpolators tudat_basic_mathematics tudat_input_output ${Boost_LIBRARIES})
 
+add_executable(test_TabulatedAerodynamicCoefficients "${SRCROOT}${AERODYNAMICSDIR}/UnitTests/unitTestTabulatedAerodynamicCoefficients.cpp")
+setup_custom_test_program(test_TabulatedAerodynamicCoefficients "${SRCROOT}${AERODYNAMICSDIR}")
+target_link_libraries(test_TabulatedAerodynamicCoefficients ${TUDAT_PROPAGATION_LIBRARIES} ${Boost_LIBRARIES})
+
 if(USE_NRLMSISE00)
     add_executable(test_NRLMSISE00Atmosphere "${SRCROOT}${AERODYNAMICSDIR}/UnitTests/unitTestNRLMSISE00Atmosphere.cpp")
     setup_custom_test_program(test_NRLMSISE00Atmosphere "${SRCROOT}${AERODYNAMICSDIR}")

Tudat/Astrodynamics/Aerodynamics/UnitTests/applicationOutput.h

@@ -0,0 +1,21 @@
+#ifndef TUDAT_APPLICATIONOUTPUT_H
+#define TUDAT_APPLICATIONOUTPUT_H
+
+namespace tudat_applications
+{
+
+//! Get path for output directory.
+static inline std::string getOutputPath( )
+{
+    // Declare file path string assigned to filePath.
+    // __FILE__ only gives the absolute path of the header file!
+    std::string filePath_( __FILE__ );
+
+    // Strip filename from temporary string and return root-path string.
+    return filePath_.substr( 0, filePath_.length( ) -
+                                std::string( "applicationOutput.h" ).length( ) );
+}
+
+}
+
+#endif // TUDAT_APPLICATIONOUTPUT_H

Tudat/Astrodynamics/Aerodynamics/UnitTests/tabulatedDragCoefficient.txt

@@ -0,0 +1,63 @@
+% Altitude [km]	Drag coefficient
+121.73913	2.0711009
+126.6087	2.0733945
+136.34783	2.0779817
+146.08696	2.0802752
+155.82609	2.0848624
+165.56522	2.0894495
+177.73913	2.0917431
+197.21739	2.0940367
+204.52174	2.0986239
+221.56522	2.0986239
+236.17391	2.103211
+250.78261	2.1055046
+270.26087	2.1055046
+292.17391	2.1077982
+318.95652	2.1100917
+343.30435	2.1100917
+372.52174	2.1123853
+394.43478	2.1169725
+423.65217	2.1215596
+448	2.1307339
+474.78261	2.1422018
+494.26087	2.1513761
+511.30435	2.1651376
+530.78261	2.1788991
+547.82609	2.190367
+562.43478	2.2041284
+579.47826	2.2224771
+598.95652	2.2408257
+611.13043	2.2545872
+620.86957	2.2683486
+633.04348	2.2889908
+645.21739	2.3050459
+654.95652	2.3233945
+667.13043	2.3348624
+676.86957	2.353211
+689.04348	2.3761468
+703.65217	2.4036697
+715.82609	2.4311927
+742.6087	2.4816514
+776.69565	2.5504587
+803.47826	2.603211
+832.69565	2.6605505
+861.91304	2.7087156
+888.69565	2.75
+908.17391	2.7821101
+927.65217	2.8004587
+947.13043	2.8211009
+969.04348	2.8440367
+993.3913	2.8600917
+1022.6087	2.8784404
+1039.6522	2.8899083
+1061.5652	2.9013761
+1090.7826	2.912844
+1122.4348	2.9243119
+1158.9565	2.9334862
+1188.1739	2.940367
+1224.6957	2.9449541
+1258.7826	2.9518349
+1295.3043	2.956422
+1326.9565	2.9587156
+1358.6087	2.9633028
+1385.3913	2.9633028

Tudat/Astrodynamics/Aerodynamics/UnitTests/unitTestTabulatedAerodynamicCoefficients.cpp

@@ -0,0 +1,259 @@
+/*    Copyright (c) 2010-2016, Delft University of Technology
+ *    All rigths reserved
+ *
+ *    This file is part of the Tudat. Redistribution and use in source and
+ *    binary forms, with or without modification, are permitted exclusively
+ *    under the terms of the Modified BSD license. You should have received
+ *    semiMajorAixscopy of the license with this file. If not, please or visit:
+ *    http://tudat.tudelft.nl/LICENSE.
+ */
+
+#define BOOST_TEST_MAIN
+
+#include <boost/test/unit_test.hpp>
+
+#include <Tudat/SimulationSetup/tudatSimulationHeader.h>
+
+#include "applicationOutput.h"
+#include "Tudat/InputOutput/basicInputOutput.h"
+
+
+namespace tudat
+{
+
+namespace unit_tests
+{
+
+BOOST_AUTO_TEST_SUITE( test_aerodynamic_acceleration_force_moment_models )
+
+BOOST_AUTO_TEST_CASE( testTabulatedDragCoefficient )
+{
+    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    ///////////////////////            USING STATEMENTS              //////////////////////////////////////////////////////
+    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+    using namespace tudat;
+    using namespace simulation_setup;
+    using namespace propagators;
+    using namespace numerical_integrators;
+    using namespace orbital_element_conversions;
+    using namespace basic_mathematics;
+    using namespace basic_astrodynamics;
+    using namespace gravitation;
+    using namespace numerical_integrators;
+    using namespace interpolators;
+    using namespace input_output;
+
+
+    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    ///////////////////////     CREATE ENVIRONMENT AND VEHICLE       //////////////////////////////////////////////////////
+    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+    // Load Spice kernels.
+    spice_interface::loadSpiceKernelInTudat( input_output::getSpiceKernelPath( ) + "pck00009.tpc" );
+    spice_interface::loadSpiceKernelInTudat( input_output::getSpiceKernelPath( ) + "de-403-masses.tpc" );
+    spice_interface::loadSpiceKernelInTudat( input_output::getSpiceKernelPath( ) + "de421.bsp" );
+
+    // Set simulation time settings.
+    const double simulationStartEpoch = 0.0;
+    const double simulationEndEpoch = tudat::physical_constants::JULIAN_DAY;
+
+    // Define body settings for simulation.
+    std::vector< std::string > bodiesToCreate;
+    bodiesToCreate.push_back( "Sun" );
+    bodiesToCreate.push_back( "Earth" );
+    bodiesToCreate.push_back( "Moon" );
+
+    // Create body objects.
+    std::map< std::string, boost::shared_ptr< BodySettings > > bodySettings =
+            getDefaultBodySettings( bodiesToCreate, simulationStartEpoch - 300.0, simulationEndEpoch + 300.0 );
+
+    for( unsigned int i = 0; i < bodiesToCreate.size( ); i++ ) {
+        bodySettings[ bodiesToCreate.at( i ) ]->ephemerisSettings->resetFrameOrientation( "J2000" );
+        bodySettings[ bodiesToCreate.at( i ) ]->rotationModelSettings->resetOriginalFrame( "J2000" );
+    }
+
+    // EARTH
+    bodySettings[ "Earth" ]->gravityFieldSettings = boost::make_shared< GravityFieldSettings >( central_spice );
+    bodySettings[ "Earth" ]->atmosphereSettings = boost::make_shared< AtmosphereSettings >( nrlmsise00 );
+
+    // MOON
+    bodySettings[ "Moon" ]->gravityFieldSettings = boost::make_shared< GravityFieldSettings >( central_spice );
+
+    NamedBodyMap bodyMap = createBodies( bodySettings );
+
+    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    ///////////////////////             CREATE VEHICLE            /////////////////////////////////////////////////////////
+    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+    // Create spacecraft object.
+    bodyMap[ "Vehicle" ] = boost::make_shared< simulation_setup::Body >( );
+    bodyMap[ "Vehicle" ]->setConstantBodyMass( 400.0 );
+    double referenceArea = 10.0;
+
+    boost::shared_ptr< AerodynamicCoefficientSettings > aerodynamicCoefficientSettings;
+
+    // Read CD as semiMajorAixsfunction of altitude
+    Eigen::MatrixXd aerodynamicsDataFromFile = readMatrixFromFile(
+                tudat_applications::getOutputPath( ) + "tabulatedDragCoefficient.txt" );
+    Eigen::VectorXd altitudeInKm = aerodynamicsDataFromFile.col(0);
+    std::vector< double > altitudes;
+    for ( unsigned int i = 0; i < altitudeInKm.size( ); i++ )
+    {
+        altitudes.push_back( altitudeInKm( i ) * 1.0E3 );
+    }
+    Eigen::VectorXd dragCoefficientsFromFile = aerodynamicsDataFromFile.col( 1 );
+    std::vector< double > dragCoefficients(
+                dragCoefficientsFromFile.data( ), dragCoefficientsFromFile.data( ) + dragCoefficientsFromFile.size( ) );
+    std::vector< Eigen::Vector3d > aerodynamicCoefficients;
+    for ( unsigned int i = 0; i < dragCoefficients.size( ); i++ )
+    {
+        aerodynamicCoefficients.push_back( dragCoefficients[ i ] * Eigen::Vector3d::UnitX( ) );
+    }
+
+    // Create interpolator
+    boost::shared_ptr< InterpolatorSettings > interpolatorSettings =
+            boost::make_shared< InterpolatorSettings >( OneDimensionalInterpolatorTypes::linear_interpolator );
+
+    // Tabulated aerodynamic settings
+    aerodynamicCoefficientSettings = boost::make_shared< TabulatedAerodynamicCoefficientSettings< 1 > >(
+                altitudes, aerodynamicCoefficients, referenceArea,
+                aerodynamics::altitude_dependent, interpolatorSettings, 1, 1 );
+
+    // Aerodynamics interface
+    bodyMap[ "Vehicle" ]->setAerodynamicCoefficientInterface(
+                createAerodynamicCoefficientInterface( aerodynamicCoefficientSettings, "Vehicle" ) );
+
+    // Finalize body creation.
+    setGlobalFrameBodyEphemerides( bodyMap, "SSB", "J2000" );
+
+    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    ///////////////////////            CREATE ACCELERATIONS          //////////////////////////////////////////////////////
+    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+    // Define propagator settings variables.
+    SelectedAccelerationMap accelerationMap;
+    std::vector< std::string > bodiesToPropagate;
+    std::vector< std::string > centralBodies;
+
+    // Define propagation settings.
+    std::map< std::string, std::vector< boost::shared_ptr< AccelerationSettings > > > accelerationsOfVehicle;
+
+    accelerationsOfVehicle[ "Earth" ].push_back( boost::make_shared< AccelerationSettings >(
+                                                     basic_astrodynamics::central_gravity ) );
+    accelerationsOfVehicle[ "Sun" ].push_back( boost::make_shared< AccelerationSettings >(
+                                                   basic_astrodynamics::central_gravity ) );
+    accelerationsOfVehicle[ "Moon" ].push_back( boost::make_shared< AccelerationSettings >(
+                                                    basic_astrodynamics::central_gravity ) );
+    accelerationsOfVehicle[ "Earth" ].push_back( boost::make_shared< AccelerationSettings >(
+                                                     basic_astrodynamics::aerodynamic ) );
+
+    accelerationMap[ "Vehicle" ] = accelerationsOfVehicle;
+
+    bodiesToPropagate.push_back( "Vehicle" );
+    centralBodies.push_back( "Earth" );
+
+    basic_astrodynamics::AccelerationMap accelerationModelMap =
+            createAccelerationModelsMap( bodyMap, accelerationMap, bodiesToPropagate, centralBodies );
+
+
+    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    ///////////////////////             CREATE PROPAGATION SETTINGS            ////////////////////////////////////////////
+    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+    double meanEarthRadius = 6371.0E3;
+    double perigeeAltitude = 150.0E3;
+    double apogeeAltitude = 35780.0E3;
+    double semiMajorAixs= (apogeeAltitude + perigeeAltitude)/2 + meanEarthRadius;
+    double eccentricity = (apogeeAltitude - perigeeAltitude)/(apogeeAltitude + perigeeAltitude + 2*meanEarthRadius);
+
+    // Set Keplerian elements for Vehicle.
+    Vector6d vehicleInitialStateInKeplerianElements;
+    vehicleInitialStateInKeplerianElements( semiMajorAxisIndex ) = semiMajorAixs;
+    vehicleInitialStateInKeplerianElements( eccentricityIndex ) = eccentricity;
+    vehicleInitialStateInKeplerianElements( inclinationIndex ) = mathematical_constants::PI / 180.0 * 23.4;
+
+    double earthGravitationalParameter = bodyMap.at( "Earth" )->getGravityFieldModel( )->getGravitationalParameter( );
+    const Vector6d vehicleInitialState = convertKeplerianToCartesianElements(
+                vehicleInitialStateInKeplerianElements, earthGravitationalParameter );
+
+    // Hybrid termination conditions
+    std::vector< boost::shared_ptr< propagators::PropagationTerminationSettings > > constituentSettings;
+
+    // Time limit
+    constituentSettings.push_back( boost::make_shared< propagators::PropagationTimeTerminationSettings >(
+                                       simulationEndEpoch ) );
+
+    // Altitude limit
+    boost::shared_ptr< PropagationTerminationSettings > altitudeTerminationSettings =
+            boost::make_shared< propagators::PropagationDependentVariableTerminationSettings >(
+                boost::make_shared< propagators::SingleDependentVariableSaveSettings >(
+                    propagators::altitude_dependent_variable, "Vehicle" ), 100.0E3, 1 );
+    constituentSettings.push_back( altitudeTerminationSettings );
+
+    // Stop if ANY of the two is met
+    boost::shared_ptr< PropagationTerminationSettings > terminationSettings =
+            boost::make_shared< propagators::PropagationHybridTerminationSettings >( constituentSettings, 1 );
+
+    // Save dependent variables
+    std::vector< boost::shared_ptr< SingleDependentVariableSaveSettings > > dependentVariablesToSave;
+    dependentVariablesToSave.push_back( boost::make_shared< SingleDependentVariableSaveSettings >(
+                                            altitude_dependent_variable, "Vehicle" ) );
+    dependentVariablesToSave.push_back( boost::make_shared< SingleDependentVariableSaveSettings >(
+                                            aerodynamic_moment_coefficients_dependent_variable, "Vehicle" ) );
+
+    boost::shared_ptr< DependentVariableSaveSettings > dependentVariableSaveSettings =
+            boost::make_shared< DependentVariableSaveSettings >( dependentVariablesToSave, 0 ) ;
+
+    // Translational propagator settings
+    boost::shared_ptr< TranslationalStatePropagatorSettings< double > > translationalPropagatorSettings =
+            boost::make_shared< TranslationalStatePropagatorSettings< double > > (
+                centralBodies, accelerationModelMap, bodiesToPropagate, vehicleInitialState, terminationSettings );
+
+
+    boost::shared_ptr< IntegratorSettings< > > integratorSettings = boost::make_shared< IntegratorSettings< > > (
+                rungeKutta4, simulationStartEpoch, 300.0 );
+
+
+    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    ///////////////////////             PROPAGATE ORBIT            ////////////////////////////////////////////////////////
+    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+    // Create simulation object and propagate dynamics.
+    SingleArcDynamicsSimulator< > dynamicsSimulator(
+                bodyMap, integratorSettings, translationalPropagatorSettings, true, false, false );
+    std::map< double, Eigen::Matrix< double, Eigen::Dynamic, 1 > > dependentVariableOutput =
+            dynamicsSimulator.getDependentVariableHistory( );
+
+
+    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    ///////////////////////       CHECK DRAG COEFFICIENTS         ////////////////////////////////////////////////////////
+    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+    // Retrieve altitudes and drag coefficients from propagation results
+
+    Eigen::Vector3d currentCoefficients;
+    double currentAltitude;
+    double computedDragCoefficient;
+    LinearInterpolatorDoublePointer dragCoefficientInterpolator =
+            boost::make_shared<LinearInterpolatorDouble>( LinearInterpolatorDouble( altitudes, dragCoefficients ) );
+    for ( std::map< double, Eigen::Matrix< double, Eigen::Dynamic, 1 > >::iterator outputIterator =
+          dependentVariableOutput.begin( ); outputIterator != dependentVariableOutput.end( ); ++outputIterator )
+    {
+        currentAltitude = outputIterator->second( 0 );
+        computedDragCoefficient = dragCoefficientInterpolator->interpolate( currentAltitude );
+        currentCoefficients = outputIterator->second.segment( 1, 3 );
+        BOOST_CHECK_CLOSE_FRACTION( computedDragCoefficient, currentCoefficients( 0 ),
+                                    4.0 * std::numeric_limits< double >::epsilon( ) );
+        BOOST_CHECK_SMALL( std::fabs( currentCoefficients( 1 ) ), std::numeric_limits< double >::epsilon( ) );
+        BOOST_CHECK_SMALL( std::fabs( currentCoefficients( 2 ) ), std::numeric_limits< double >::epsilon( ) );
+
+    }
+
+}
+
+BOOST_AUTO_TEST_SUITE_END( )
+
+}
+
+}

Tudat/Astrodynamics/Aerodynamics/aerodynamicCoefficientInterface.h

@@ -66,7 +66,8 @@ enum AerodynamicCoefficientsIndependentVariables
     mach_number_dependent = 0,
     angle_of_attack_dependent = 1,
     angle_of_sideslip_dependent = 2,
-    undefined_independent_variable = 3
+    altitude_dependent = 3,
+    undefined_independent_variable = 4
 };
 
 //! Base class to hold an aerodynamic coefficient interface.

Tudat/Astrodynamics/Aerodynamics/flightConditions.cpp

@@ -60,7 +60,8 @@ void FlightConditions::setAerodynamicCoefficientsIndependentVariableFunction(
 {
     if( ( independentVariable == mach_number_dependent ) ||
             ( independentVariable == angle_of_attack_dependent ) ||
-            ( independentVariable == angle_of_sideslip_dependent ) )
+            ( independentVariable == angle_of_sideslip_dependent )||
+            ( independentVariable == altitude_dependent ) )
     {
         throw std::runtime_error(
                     std::string( "Error when setting aerodynamic coefficient function dependency, value of parameter " ) +
@@ -162,6 +163,13 @@ void FlightConditions::updateAerodynamicCoefficientInput( )
                         aerodynamicAngleCalculator_->getAerodynamicAngle(
                             reference_frames::angle_of_sideslip ) );
             break;
+        case altitude_dependent:
+            if( aerodynamicAngleCalculator_== NULL )
+            {
+                throw std::runtime_error( "Error, aerodynamic angle calculator is null, but require angle of sideslip" );
+            }
+            aerodynamicCoefficientIndependentVariables_.push_back( currentAltitude_ );
+            break;
         default:
             if( customCoefficientDependencies_.count(
                         aerodynamicCoefficientInterface_->getIndependentVariableName( i ) ) == 0 )