From 9f91afdacd4739ab271ff4761b5f77e5e6d2e9c6 Mon Sep 17 00:00:00 2001 From: Rene107 Date: Mon, 21 Mar 2016 16:07:21 +0100 Subject: [PATCH 1/9] CMakeList comment changed Find Cmake added folders --- Tudat/CMakeLists.txt | 2 +- Tudat/External/CMake/FindNRLMSISE00.cmake | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/Tudat/CMakeLists.txt b/Tudat/CMakeLists.txt index c8db9e97e7..bf3752323d 100755 --- a/Tudat/CMakeLists.txt +++ b/Tudat/CMakeLists.txt @@ -197,7 +197,7 @@ endif() # # NRLMSISE-00 # -# Set whether to use the SPICE library integration with Tudat or not. If it not supplied by the +# Set whether to use the NRLMSISE-00 library integration with Tudat or not. If it not supplied by the # user (either directly as an argument or through the "UserSettings.txt" file, the default setting # is "OFF". option(USE_NRLMSISE00 "build Tudat with NRLMSISE-00 enabled" OFF) diff --git a/Tudat/External/CMake/FindNRLMSISE00.cmake b/Tudat/External/CMake/FindNRLMSISE00.cmake index f54095d5b5..da1883c0b8 100644 --- a/Tudat/External/CMake/FindNRLMSISE00.cmake +++ b/Tudat/External/CMake/FindNRLMSISE00.cmake @@ -43,6 +43,8 @@ if(NOT NRLMSISE00_BASE_PATH) ${PROJECT_SOURCE_DIR}/../.. ${PROJECT_SOURCE_DIR}/../../.. ${PROJECT_SOURCE_DIR}/../../../.. + ${PROJECT_SOURCE_DIR}/../../../../.. + ${PROJECT_SOURCE_DIR}/../../../../../.. PATH_SUFFIXES nrlmsise-00 ) endif(NOT NRLMSISE00_BASE_PATH) From a219c39ae411e065ea3365dbcaba8cc3382cb6f2 Mon Sep 17 00:00:00 2001 From: Reneh107 Date: Wed, 4 May 2016 21:31:25 +0200 Subject: [PATCH 2/9] Added virtual function resetHashKey, which is used in NRLMSISE Add molar gas constant Added mean free path + speed of sound + molar mass + Gas properties datastructure Changed altitude input to meters and radians + Added mean molar mass + speed of sound + mean free path etc. Added stream output for nrlmsise input data structure --- .../unitTestNRLMSISE00Atmosphere.cpp | 24 +- .../Aerodynamics/atmosphereModel.h | 176 ++++++++++++ .../Aerodynamics/nrlmsise00Atmosphere.cpp | 117 ++++++-- .../Aerodynamics/nrlmsise00Atmosphere.h | 258 ++++++++++++++---- .../BasicAstrodynamics/physicalConstants.h | 6 + 5 files changed, 489 insertions(+), 92 deletions(-) diff --git a/Tudat/Astrodynamics/Aerodynamics/UnitTests/unitTestNRLMSISE00Atmosphere.cpp b/Tudat/Astrodynamics/Aerodynamics/UnitTests/unitTestNRLMSISE00Atmosphere.cpp index faab04d726..ad49582d8a 100644 --- a/Tudat/Astrodynamics/Aerodynamics/UnitTests/unitTestNRLMSISE00Atmosphere.cpp +++ b/Tudat/Astrodynamics/Aerodynamics/UnitTests/unitTestNRLMSISE00Atmosphere.cpp @@ -57,6 +57,8 @@ #include "Tudat/Astrodynamics/Aerodynamics/tabulatedAtmosphere.h" #include "Tudat/InputOutput/basicInputOutput.h" +#include "tudat/Mathematics/BasicMathematics/mathematicalConstants.h" + namespace tudat { @@ -70,12 +72,14 @@ BOOST_AUTO_TEST_SUITE( test_nrlmsise00_atmosphere ) using tudat::aerodynamics::NRLMSISE00Input; using tudat::aerodynamics::NRLMSISE00Atmosphere; +using tudat::mathematical_constants::PI; + // Global variable to be changed by tests and function. NRLMSISE00Input data; // Define input data generic (or almost completely) for all tests. NRLMSISE00Input gen_data(0, 172, 29000.0, 16.0, 150.0, 150.0, 4.0); -std::vector< double > gen_input = boost::assign::list_of(400.0)(-70.0)(60.0)(0.0); +std::vector< double > gen_input = boost::assign::list_of(400.0E3)(-70.0*PI/180.0)(60.0*PI/180.0)(0.0); //! Test function to update the NRLMSISE00 iunput data as a function of time and geometry. NRLMSISE00Input function( double altitude, double longitude, @@ -301,7 +305,7 @@ BOOST_AUTO_TEST_CASE( testNRLMSISE00AtmosphereTest3 ) std::vector< double > input = gen_input; // Define variations from the standard case data.secondOfTheDay = 75000.0; - input[ 0 ] = 1000.0; + input[0] = 1000.0E3; // We change parameters other than alt, lat, long, and/or time, // so it's best to reset the hash. @@ -346,7 +350,7 @@ BOOST_AUTO_TEST_CASE( testNRLMSISE00AtmosphereTest4 ) std::vector< double > input = gen_input; // Define variations from the standard case - input[ 0 ] = 100.0; + input[0] = 100.0E3; // We change parameters other than alt, lat, long, and/or time, // so it's best to reset the hash. @@ -706,7 +710,7 @@ BOOST_AUTO_TEST_CASE( testNRLMSISE00AtmosphereTest12 ) std::vector< double > input = gen_input; // Define variations from the standard case - input[ 0 ] = 10.0; + input[0] = 10.0E3; // We change parameters other than alt, lat, long, and/or time, // so it's best to reset the hash. @@ -751,7 +755,7 @@ BOOST_AUTO_TEST_CASE( testNRLMSISE00AtmosphereTest13 ) std::vector< double > input = gen_input; // Define variations from the standard case - input[ 0 ] = 30.0; + input[0] = 30.0E3; // We change parameters other than alt, lat, long, and/or time, // so it's best to reset the hash. @@ -796,7 +800,7 @@ BOOST_AUTO_TEST_CASE( testNRLMSISE00AtmosphereTest14 ) std::vector< double > input = gen_input; // Define variations from the standard case - input[ 0 ] = 50.0; + input[0] = 50.0E3; // We change parameters other than alt, lat, long, and/or time, // so it's best to reset the hash. @@ -841,7 +845,7 @@ BOOST_AUTO_TEST_CASE( testNRLMSISE00AtmosphereTest15 ) std::vector< double > input = gen_input; // Define variations from the standard case - input[ 0 ] = 70.0; + input[0] = 70.0E3; // We change parameters other than alt, lat, long, and/or time, // so it's best to reset the hash. @@ -932,9 +936,9 @@ BOOST_AUTO_TEST_CASE( testNRLMSISE00AtmosphereTest17 ) std::vector< double > input = gen_input; // Define variations from the standard case - input[ 0 ] = 100.0; - data.apVector = std::vector< double >( 7, 100.0 ); - data.switches[ 9 ] = -1; + input[0] = 100.0E3; + data.apVector = std::vector< double >(7, 100.0); + data.switches[9] = -1; // We change parameters other than alt, lat, long, and/or time, // so it's best to reset the hash. diff --git a/Tudat/Astrodynamics/Aerodynamics/atmosphereModel.h b/Tudat/Astrodynamics/Aerodynamics/atmosphereModel.h index b7a305a121..5cec49818b 100644 --- a/Tudat/Astrodynamics/Aerodynamics/atmosphereModel.h +++ b/Tudat/Astrodynamics/Aerodynamics/atmosphereModel.h @@ -42,12 +42,82 @@ #define TUDAT_ATMOSPHERE_MODEL_H #include +#include namespace tudat { namespace aerodynamics { +//! Gas component properties data structure +/*! + * This data structure contains the molar mass and collision diameter of + * the most frequently observed gasses in the atmosphere. + */ +struct GasComponentProperties +{ + //! default constructor + /*! + * Constructs a GasComponentProperties data structure with values obtained from the following sources: + * Collision diameters from: Tables of Physical & Chemical Constants Kaye & Laby Online, Kaye & Laby Online, 2016 + * (http://www.kayelaby.npl.co.uk/general_physics/2_2/2_2_4.html) + * Molar mass from: NIST,2016 (http://www.nist.gov/pml/data/images/illo_for_2014_PT_1.PNG) + * + * Data to be verified + * + */ + GasComponentProperties(): + diameterArgon(340E-12), diameterAtomicHydrogen(260E-12), diameterHelium(256E-12), + diameterNitrogen(370E-12), diameterOxygen(358E-12), diameterAtomicNitrogen(290E-12), + diameterAtomicOxygen(280E-12), molarMassArgon(39.948E-3), molarMassAtomicHydrogen(1.008E-3), + molarMassHelium(4.002602E-3), molarMassNitrogen(2.0*14.007E-3), molarMassOxygen(2.0*15.999E-3), + molarMassAtomicNitrogen(14.007E-3), molarMassAtomicOxygen(15.999E-3) + { } + + //! Molecular colision diameter of Argon in m + double diameterArgon; + + //! Molecular colision diameter of Atomic Hydrogen in m + double diameterAtomicHydrogen; + + //! Molecular colision diameter of Helium in m + double diameterHelium; + + //! Molecular colision diameter of Nitrogen in m + double diameterNitrogen; + + //! Molecular colision diameter of Oxygen in m + double diameterOxygen; + + //! Molecular colision diameter of Atomic Nitrogen in m + double diameterAtomicNitrogen; + + //! Molecular colision diameter of Atomic Oxygen in m + double diameterAtomicOxygen; + + + //! molar mass of Argon in kg/mole + double molarMassArgon; + + //! Molar mass of Atomic Hydrogen in kg/mole + double molarMassAtomicHydrogen; + + //! Molar mass of Helium in kg/mole + double molarMassHelium; + + //! Molar mass of Nitrogen in kg/mole + double molarMassNitrogen; + + //! Molar mass of Oxygen in kg/mole + double molarMassOxygen; + + //! Molar mass of Atomic Nitrogen in kg/mole + double molarMassAtomicNitrogen; + + //! Molar mass of Atomic Oxygen in kg/mole + double molarMassAtomicOxygen; +}; + //! Atmosphere model class. /*! * Base class for all atmosphere models. @@ -64,6 +134,14 @@ class AtmosphereModel */ virtual ~AtmosphereModel( ) { } + //! Set gas component properties. + /*! + * Sets the gas component properties, which contains the molar mass and collision diameter of the molecules. + * \param GasComponentProperties . + * \return void. + */ + virtual void setGasComponentProperties( GasComponentProperties GasComponentProperties ){} + //! Get local density. /*! * Returns the local density parameter of the atmosphere in kg per meter^3. @@ -100,6 +178,104 @@ class AtmosphereModel virtual double getTemperature( const double altitude, const double longitude, const double latitude, const double time ) = 0; + //! Get Speed of Sound + /*! + * This function returns the speed of sound. + * The function is not implemented by default and returns NaN if the atmosphere model doesn't implement this function. + * \param altitude Altitude. + * \param longitude Longitude. + * \param latitude Latitude. + * \param time Time. + * \return double SpeedOfSound + */ + virtual double getSpeedOfSound(const double altitude, const double longitude, + const double latitude, const double time){ + return std::numeric_limits::quiet_NaN(); + } + + //! Get mean free path. + /*! + * This function returns the mean free path, which can be calculated using the number density and the collision diameter. + * The function is not implemented by default and returns NaN if the atmosphere model doesn't implement this function. + * \param altitude Altitude. + * \param longitude Longitude. + * \param latitude Latitude. + * \param time Time. + * \return double mean free path. + */ + virtual double getMeanFreePath(const double altitude, const double longitude, + const double latitude, const double time){ + return std::numeric_limits::quiet_NaN(); + } + + //! Get number densities. + /*! + * This function returns a vector with number densities for all the components of the gas(air). + * The function is not implemented by default and returns null vector if the atmosphere model doesn't implement this function. + * \param altitude Altitude. + * \param longitude Longitude. + * \param latitude Latitude. + * \param time Time. + * \return std::vector number densities. + */ + virtual std::vector getNumberDensities(const double altitude, const double longitude, + const double latitude, const double time){ + std::vector null(0); + return null; + } + + //! Get mean molar mass. + /*! + * This function returns the mean molar mass at the current location and time. + * The function is not implemented by default and returns nan if the atmosphere model doesn't implement this function. + * \param altitude Altitude. + * \param longitude Longitude. + * \param latitude Latitude. + * \param time Time. + * \return double meanMolarMass. + */ + virtual double getMeanMolarMass(const double altitude, const double longitude, + const double latitude, const double time){ + return std::numeric_limits::quiet_NaN(); + } + + //! Get average number density + /*! + * This function returns the (unweighted) average number density at the current location and time. + * The function is not implemented by default and returns nan if the atmosphere model doesn't implement this function. + * \param altitude Altitude. + * \param longitude Longitude. + * \param latitude Latitude. + * \param time Time. + * \return double average number density. + */ + virtual double getAverageNumberDensity(const double altitude, const double longitude, + const double latitude, const double time){ + return std::numeric_limits::quiet_NaN(); + } + + //! Get weighted average collision diameter + /*! + * This function returns the weighted average collision diameter of the gas components at the current location and time. + * The function is not implemented by default and returns nan if the atmosphere model doesn't implement this function. + * \param altitude Altitude. + * \param longitude Longitude. + * \param latitude Latitude. + * \param time Time. + * \return double weighted average collision diameter. + */ + virtual double getWeightedAverageCollisionDiameter(const double altitude, const double longitude, + const double latitude, const double time){ + return std::numeric_limits::quiet_NaN(); + } + + //! Reset Hash key + /*! + * This function resets the hashKey, which is used to check if recalculation of the atmospheric properties is required or not. + * \return void + */ + virtual void resetHashKey(){} + //! Get local speed of sound. /*! * Returns the local speed of sound of the atmosphere in m/s. diff --git a/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.cpp b/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.cpp index 9a9acdefcc..b351f0d675 100644 --- a/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.cpp +++ b/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.cpp @@ -31,6 +31,7 @@ */ #include #include "Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.h" +#include "tudat/Mathematics/BasicMathematics/mathematicalConstants.h" //! Tudat library namespace. @@ -39,7 +40,8 @@ namespace tudat namespace aerodynamics { -//! Compute the local atmospheric properties. +using tudat::mathematical_constants::PI; + void NRLMSISE00Atmosphere::computeProperties(double altitude, double longitude, double latitude, double time) { // Compute the hash key @@ -58,28 +60,97 @@ void NRLMSISE00Atmosphere::computeProperties(double altitude, double longitude, std::copy( inputData.apVector.begin( ), inputData.apVector.end( ), aph_.a ); std::copy( inputData.switches.begin( ), inputData.switches.end(), flags_.switches); - // Set environment properties. - input_.g_lat = latitude; - input_.g_long = longitude; - input_.alt = altitude; - input_.year = inputData.year; - input_.doy = inputData.dayOfTheYear; - input_.sec = inputData.secondOfTheDay; - input_.lst = inputData.localSolarTime; - input_.f107 = inputData.f107; - input_.f107A = inputData.f107a; - input_.ap = inputData.apDaily; - input_.ap_a = &aph_; - - // Call Neutral Atmosphere Empircial Model from the surface to lower exosphere from NRLMSISE00. - gtd7( &input_, &flags_, &output_ ); - - // Compute current density and temperature. - density_ = output_.d[ 5 ] * 1000.0; - temperature_ = output_.t[ 1 ]; - if( useIdealGasLaw_ ) - { - pressure_ = density_ * physical_constants::SPECIFIC_GAS_CONSTANT_AIR * temperature_; + std::copy(inputData.apVector.begin(), + inputData.apVector.end(), aph_.a); + std::copy(inputData.switches.begin(), + inputData.switches.end(), flags_.switches); + + input_.g_lat = latitude*180.0/PI; // rad to deg + input_.g_long = longitude*180.0/PI; // rad to deg + input_.alt = altitude*1.0E-3; // m to km + input_.year = inputData.year; + input_.doy = inputData.dayOfTheYear; + input_.sec = inputData.secondOfTheDay; + input_.lst = inputData.localSolarTime; + input_.f107 = inputData.f107; + input_.f107A = inputData.f107a; + input_.ap = inputData.apDaily; + input_.ap_a = &aph_; + + gtd7(&input_, &flags_, &output_); + + density_ = output_.d[5]*1000.0; // GM/CM3 to kg/M3 + temperature_ = output_.t[1]; + pressure_ = TUDAT_NAN; + + // Get number densities + numberDensities_.resize(7); + numberDensities_[0] = output_.d[0] * 1.0E6 ; // HE NUMBER DENSITY (M-3) + numberDensities_[1] = output_.d[1] * 1.0E6 ; // O NUMBER DENSITY (M-3) + numberDensities_[2] = output_.d[2] * 1.0E6 ; // N2 NUMBER DENSITY (M-3) + numberDensities_[3] = output_.d[3] * 1.0E6 ; // O2 NUMBER DENSITY (M-3) + numberDensities_[4] = output_.d[4] * 1.0E6 ; // AR NUMBER DENSITY (M-3) + numberDensities_[5] = output_.d[6] * 1.0E6 ; // H NUMBER DENSITY (M-3) + numberDensities_[6] = output_.d[7] * 1.0E6 ; // N NUMBER DENSITY (M-3) + numberDensities_[7] = output_.d[8] * 1.0E6 ; // Anomalous oxygen NUMBER DENSITY (M-3) + + // Get average number density + double sumOfNumberDensity = 0.0 ; + for(unsigned int i = 0 ; i < numberDensities_.size() ; i++){ + sumOfNumberDensity += numberDensities_[i] ; + } + averageNumberDensity_ = sumOfNumberDensity / double( numberDensities_.size() ) ; + + // Mean molar mass (SOURCE) + meanMolarMass_ = numberDensities_[0] * gasComponentProperties_.molarMassHelium; + meanMolarMass_ += numberDensities_[1] * gasComponentProperties_.molarMassAtomicOxygen; + meanMolarMass_ += numberDensities_[2] * gasComponentProperties_.molarMassNitrogen; + meanMolarMass_ += numberDensities_[3] * gasComponentProperties_.molarMassOxygen; + meanMolarMass_ += numberDensities_[4] * gasComponentProperties_.molarMassArgon; + meanMolarMass_ += numberDensities_[5] * gasComponentProperties_.molarMassAtomicHydrogen; + meanMolarMass_ += numberDensities_[6] * gasComponentProperties_.molarMassAtomicNitrogen; + meanMolarMass_ += numberDensities_[7] * gasComponentProperties_.molarMassOxygen; + meanMolarMass_ = meanMolarMass_ / sumOfNumberDensity ; + + // Speed of sound (Anderson, Fundamentals of Aerodynamics, 2006) + speedOfSound_ = std::sqrt( specificHeatRatio_ * molarGasConstant_ * temperature_ / meanMolarMass_ ); + + // Collision diameter (SOURCE) + weightedAverageCollisionDiameter_ = numberDensities_[0]* gasComponentProperties_.diameterHelium ; + weightedAverageCollisionDiameter_ += numberDensities_[1]* gasComponentProperties_.diameterAtomicOxygen ; + weightedAverageCollisionDiameter_ += numberDensities_[2]* gasComponentProperties_.diameterNitrogen ; + weightedAverageCollisionDiameter_ += numberDensities_[3]* gasComponentProperties_.diameterOxygen ; + weightedAverageCollisionDiameter_ += numberDensities_[4]* gasComponentProperties_.diameterArgon ; + weightedAverageCollisionDiameter_ += numberDensities_[5]* gasComponentProperties_.diameterAtomicHydrogen ; + weightedAverageCollisionDiameter_ += numberDensities_[6]* gasComponentProperties_.diameterAtomicNitrogen ; + weightedAverageCollisionDiameter_ += numberDensities_[7]* gasComponentProperties_.diameterAtomicOxygen ; + weightedAverageCollisionDiameter_ = weightedAverageCollisionDiameter_ / sumOfNumberDensity; + + // Mean free path (Chapman, S. & Cowling, T. The mathematical theory of nonuniform gases Cambridge University Press, 1970) + meanFreePath_ = (std::pow( std::sqrt(2.0) * tudat::mathematical_constants::PI * std::pow(weightedAverageCollisionDiameter_,2.0) * + averageNumberDensity_ , -1.0 )) ; +} + +//! Overloaded ostream to print class information. +std::ostream& operator<<( std::ostream& stream, + NRLMSISE00Input& nrlmsiseInput ){ + stream << "This is a NRLMSISE Input data object." << std::endl; + stream << "The input data is stored as: " << std::endl; + + stream << "Year = " << nrlmsiseInput.year << std::endl; + stream << "Day of the year = " << nrlmsiseInput.dayOfTheYear << std::endl; + stream << "Second of the day = " << nrlmsiseInput.secondOfTheDay << std::endl; + stream << "Local solar time = " << nrlmsiseInput.localSolarTime << std::endl; + stream << "f107 = " << nrlmsiseInput.f107 << std::endl; + stream << "f107a = " << nrlmsiseInput.f107a << std::endl; + stream << "apDaily = " << nrlmsiseInput.apDaily << std::endl; + + for(unsigned int i = 0 ; i < nrlmsiseInput.apVector.size() ; i++){ + stream << "apVector[ " << i << " ] = " << nrlmsiseInput.apVector[i] << std::endl; + } + + for(unsigned int i = 0 ; i < nrlmsiseInput.switches.size() ; i++){ + stream << "switches[ " << i << " ] = " << nrlmsiseInput.switches[i] << std::endl; } else { diff --git a/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.h b/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.h index 6dc8ed6161..c185c4d2d1 100644 --- a/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.h +++ b/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.h @@ -137,24 +137,51 @@ class NRLMSISE00Atmosphere : public AtmosphereModel { * (altitude, longitude, latitude, time). * \param useIdealGasLaw Variable denoting whether to use the ideal gas law for computation of pressure. */ - NRLMSISE00Atmosphere( const NRLMSISE00InputFunction nrlmsise00InputFunction, - const bool useIdealGasLaw = false ) - :nrlmsise00InputFunction_( nrlmsise00InputFunction ), useIdealGasLaw_( useIdealGasLaw ) - { - if( useIdealGasLaw ) - { - std::cerr<<"Warning, using ideal gas law with R at standard atmospheric conditions in pressure computation of NRLMSISE00InputFunction."< getNumberDensities(const double altitude, const double longitude, + const double latitude, const double time) { + computeProperties(altitude, longitude, latitude, time); + return numberDensities_; + } + + //! Get local average number density. + /*! + * Returns the local average number density in m^-3 + * \param altitude Altitude [m]. + * \param longitude Longitude [rad]. + * \param latitude Latitude [rad]. + * \param time time. + * \return average number density. + */ + double getAverageNumberDensity(const double altitude, const double longitude, + const double latitude, const double time) { + computeProperties(altitude, longitude, latitude, time); + return averageNumberDensity_; + } + + //! Get local weighted average collision diameter. + /*! + * Returns the local weighted average collision diameter using the number densities as weights. + * \param altitude Altitude [m]. + * \param longitude Longitude [rad]. + * \param latitude Latitude [rad]. + * \param time time. + * \return weighted average collision diameter. + */ + double getWeightedAverageCollisionDiameter(const double altitude, const double longitude, + const double latitude, const double time) { + computeProperties(altitude, longitude, latitude, time); + return weightedAverageCollisionDiameter_; } //! Get the full model output @@ -232,10 +326,10 @@ class NRLMSISE00Atmosphere : public AtmosphereModel { * Gets the output directly from the model. This will return a * pair of double vectors containing density and temperature * values. - * \param altitude Altitude. - * \param longitude Longitude. - * \param latitude Latitude. - * \param time Time. + * \param altitude Altitude [m]. + * \param longitude Longitude [rad]. + * \param latitude Latitude [rad]. + * \param time time. * \return Full density and temperature values */ std::pair< std::vector< double >, std::vector< double > > getFullOutput( @@ -260,19 +354,65 @@ class NRLMSISE00Atmosphere : public AtmosphereModel { //! Current key hash size_t hashKey_; - //! Current local density + /*! + * Current local density (kg/m3) + */ double density_; - //! Current local temperature + /*! + * Current local temperature (K) + */ double temperature_; - //! Current local pressure + /*! + * Current local pressure (Not implemented!) + */ double pressure_; - //! Variable denoting whether to use the ideal gas law for computation of pressure. - bool useIdealGasLaw_; + /*! + * Current speed of sound (m/s) + */ + double speedOfSound_; + + /*! + * Current mean free path (m) + */ + double meanFreePath_; + + /*! + * Current number densities of gas components + * numberDensities_[0] - HE NUMBER DENSITY (M-3) + * numberDensities_[1] - O NUMBER DENSITY (M-3) + * numberDensities_[2] - N2 NUMBER DENSITY (M-3) + * numberDensities_[3] - O2 NUMBER DENSITY (M-3) + * numberDensities_[4] - AR NUMBER DENSITY (M-3) + * numberDensities_[5] - H NUMBER DENSITY (M-3) + * numberDensities_[6] - N NUMBER DENSITY (M-3) + * numberDensities_[7] - Anomalous oxygen NUMBER DENSITY (M-3) + */ + std::vector numberDensities_; + + //! Current average number density (M-3) + double averageNumberDensity_; - //! The file name of the solar activity data. + //! Current weighted average of the collision diameter using the number density as weights in (M) + double weightedAverageCollisionDiameter_; + + //! mean molar mass (kg/mole) + double meanMolarMass_; + + //! Data structure that contains the colision diameter + GasComponentProperties gasComponentProperties_; + + //! Specific heat ratio + double specificHeatRatio_; + + //! Molar gas constant (J/mol K) + double molarGasConstant_; + + /*! + * The file name of the solar activity data. + */ struct nrlmsise_flags flags_; //! Magnetic index structure of 7d array. @@ -297,8 +437,8 @@ class NRLMSISE00Atmosphere : public AtmosphereModel { * doy - day of the year (int) * sec - seconds in the day (UT) (double) * alt - altitude in kilometers (double) - * g_lat - geodetic latitude (double) - * g_long - geodetic longitude (double) + * g_lat - geodetic latitude in deg (double) + * g_long - geodetic longitude in deg (double) * lst - local apparent solar time (hours) (double) * f107A - 81 day average of F10.7 flux (centered on doy) (double) * f107 - daily F10.7 flux for previous day (double) @@ -346,10 +486,10 @@ class NRLMSISE00Atmosphere : public AtmosphereModel { //! Get Hash Key. /*! * Returns hash key value based on a vector of keys - * \param altitude Altitude. - * \param longitude Longitude. - * \param latitude Latitude. - * \param time Time. + * \param altitude Altitude [m]. + * \param longitude Longitude [rad]. + * \param latitude Latitude [rad]. + * \param time time. * \return hash key value. */ size_t hashFunc( const double altitude, const double longitude, @@ -366,10 +506,10 @@ class NRLMSISE00Atmosphere : public AtmosphereModel { //! Compute the local atmospheric properties. /*! * Computes the local atmospheric density, pressure and temperature. - * \param altitude Altitude. - * \param longitude Longitude. - * \param latitude Latitude. - * \param time Time. + * \param altitude Altitude [m]. + * \param longitude Longitude [rad]. + * \param latitude Latitude [rad]. + * \param time time. */ void computeProperties( const double altitude, const double longitude, const double latitude, const double time ); diff --git a/Tudat/Astrodynamics/BasicAstrodynamics/physicalConstants.h b/Tudat/Astrodynamics/BasicAstrodynamics/physicalConstants.h index 59e8baea6d..c326da9c35 100644 --- a/Tudat/Astrodynamics/BasicAstrodynamics/physicalConstants.h +++ b/Tudat/Astrodynamics/BasicAstrodynamics/physicalConstants.h @@ -129,6 +129,12 @@ const static double ASTRONOMICAL_UNIT = 1.49597870691e11; */ const static double SPECIFIC_GAS_CONSTANT_AIR = 2.87e2; +//! Molar gas constant. +/*! + * The molar gas constant in J per mole Kelvin (J/(mol K)) (NIST: http://physics.nist.gov/cgi-bin/cuu/Value?r, 2016). + */ +const static double MOLAR_GAS_CONSTANT = 8.3144598; + //! Planck constant. /*! * Planck's constant in m^{2} kg/s, (NIST, 2013). From 82826785df943acf88f442622f0644397962d564 Mon Sep 17 00:00:00 2001 From: Reneh107 Date: Wed, 4 May 2016 22:28:15 +0200 Subject: [PATCH 3/9] Minor changes + Computation atmospheric pressure --- .../Aerodynamics/atmosphereModel.h | 25 ++-------- .../Aerodynamics/nrlmsise00Atmosphere.cpp | 16 ++++--- .../Aerodynamics/nrlmsise00Atmosphere.h | 48 ++++++++++++------- 3 files changed, 47 insertions(+), 42 deletions(-) diff --git a/Tudat/Astrodynamics/Aerodynamics/atmosphereModel.h b/Tudat/Astrodynamics/Aerodynamics/atmosphereModel.h index 5cec49818b..204fa2ed26 100644 --- a/Tudat/Astrodynamics/Aerodynamics/atmosphereModel.h +++ b/Tudat/Astrodynamics/Aerodynamics/atmosphereModel.h @@ -178,20 +178,17 @@ class AtmosphereModel virtual double getTemperature( const double altitude, const double longitude, const double latitude, const double time ) = 0; - //! Get Speed of Sound + //! Get local speed of sound. /*! - * This function returns the speed of sound. - * The function is not implemented by default and returns NaN if the atmosphere model doesn't implement this function. + * Returns the local speed of sound of the atmosphere in m/s. * \param altitude Altitude. * \param longitude Longitude. * \param latitude Latitude. * \param time Time. - * \return double SpeedOfSound + * \return Atmospheric speed of sound. */ - virtual double getSpeedOfSound(const double altitude, const double longitude, - const double latitude, const double time){ - return std::numeric_limits::quiet_NaN(); - } + virtual double getSpeedOfSound( const double altitude, const double longitude, + const double latitude, const double time ) = 0; //! Get mean free path. /*! @@ -276,18 +273,6 @@ class AtmosphereModel */ virtual void resetHashKey(){} - //! Get local speed of sound. - /*! - * Returns the local speed of sound of the atmosphere in m/s. - * \param altitude Altitude. - * \param longitude Longitude. - * \param latitude Latitude. - * \param time Time. - * \return Atmospheric speed of sound. - */ - virtual double getSpeedOfSound( const double altitude, const double longitude, - const double latitude, const double time ) = 0; - protected: private: diff --git a/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.cpp b/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.cpp index b351f0d675..a75e08ff43 100644 --- a/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.cpp +++ b/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.cpp @@ -81,7 +81,6 @@ void NRLMSISE00Atmosphere::computeProperties(double altitude, double longitude, density_ = output_.d[5]*1000.0; // GM/CM3 to kg/M3 temperature_ = output_.t[1]; - pressure_ = TUDAT_NAN; // Get number densities numberDensities_.resize(7); @@ -101,7 +100,7 @@ void NRLMSISE00Atmosphere::computeProperties(double altitude, double longitude, } averageNumberDensity_ = sumOfNumberDensity / double( numberDensities_.size() ) ; - // Mean molar mass (SOURCE) + // Mean molar mass (Thermodynamics an Engineering Approach, Michael A. Boles) meanMolarMass_ = numberDensities_[0] * gasComponentProperties_.molarMassHelium; meanMolarMass_ += numberDensities_[1] * gasComponentProperties_.molarMassAtomicOxygen; meanMolarMass_ += numberDensities_[2] * gasComponentProperties_.molarMassNitrogen; @@ -129,6 +128,14 @@ void NRLMSISE00Atmosphere::computeProperties(double altitude, double longitude, // Mean free path (Chapman, S. & Cowling, T. The mathematical theory of nonuniform gases Cambridge University Press, 1970) meanFreePath_ = (std::pow( std::sqrt(2.0) * tudat::mathematical_constants::PI * std::pow(weightedAverageCollisionDiameter_,2.0) * averageNumberDensity_ , -1.0 )) ; + + // Calculate pressure using ideal gas law (Thermodynamics an Engineering Approach, Michael A. Boles) + if( useIdealGasLaw_ ){ + pressure_ = density_ * molarGasConstant_ * temperature_ / meanMolarMass_ ; + } + else{ + pressure_ = TUDAT_NAN; + } } //! Overloaded ostream to print class information. @@ -152,11 +159,8 @@ std::ostream& operator<<( std::ostream& stream, for(unsigned int i = 0 ; i < nrlmsiseInput.switches.size() ; i++){ stream << "switches[ " << i << " ] = " << nrlmsiseInput.switches[i] << std::endl; } - else - { - pressure_ = TUDAT_NAN; - } + return stream; } //! Get the full model output diff --git a/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.h b/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.h index c185c4d2d1..6233b43be7 100644 --- a/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.h +++ b/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.h @@ -117,15 +117,15 @@ struct NRLMSISE00Input * NRLMSISE-00 atmosphere model class. This class uses the NRLMSISE00 atmosphere model to calculate atmospheric * density and temperature. The GTD7 function is used: Neutral Atmosphere Empircial Model from the surface to lower * exosphere. - * Compositional data is currently not used (computation of pressure; speed of sound, etc.). - * + * Currently the ideal gas law is used to compute the speed of sound. + * The specific heat ratio is assumed to be constant and equal to 1.4. */ class NRLMSISE00Atmosphere : public AtmosphereModel { public: - //! Solar activity function + //! NRLMSISEInput function /*! - * Boost function that accepts (altitude, longitude, latitude, time) and returns solar activity data. + * Boost function that accepts (altitude, longitude, latitude, time) and returns NRLMSISEInput data. */ typedef boost::function< NRLMSISE00Input( double, double, double, double ) > NRLMSISE00InputFunction; @@ -137,13 +137,16 @@ class NRLMSISE00Atmosphere : public AtmosphereModel { * (altitude, longitude, latitude, time). * \param useIdealGasLaw Variable denoting whether to use the ideal gas law for computation of pressure. */ - NRLMSISE00Atmosphere(const NRLMSISE00InputFunction nrlmsise00InputFunction) - : nrlmsise00InputFunction_(nrlmsise00InputFunction) { + NRLMSISE00Atmosphere(const NRLMSISE00InputFunction nrlmsise00InputFunction, + bool useIdealGasLaw = false) + : nrlmsise00InputFunction_(nrlmsise00InputFunction) + { resetHashKey(); molarGasConstant_ = tudat::physical_constants::MOLAR_GAS_CONSTANT; specificHeatRatio_ = 1.4 ; GasComponentProperties gasProperties; gasComponentProperties_ = gasProperties; // Default gas properties + useIdealGasLaw_ = useIdealGasLaw; } //! Constructor @@ -156,12 +159,15 @@ class NRLMSISE00Atmosphere : public AtmosphereModel { */ NRLMSISE00Atmosphere(const NRLMSISE00InputFunction nrlmsise00InputFunction, double specificHeatRatio, - GasComponentProperties gasProperties) - : nrlmsise00InputFunction_(nrlmsise00InputFunction) { + GasComponentProperties gasProperties, + bool useIdealGasLaw = false) + : nrlmsise00InputFunction_(nrlmsise00InputFunction) + { resetHashKey(); molarGasConstant_ = tudat::physical_constants::MOLAR_GAS_CONSTANT; specificHeatRatio_ = specificHeatRatio ; - gasComponentProperties_ = gasProperties; // Default gas properties + gasComponentProperties_ = gasProperties; + useIdealGasLaw_ = useIdealGasLaw ; } //! Set gas component properties. @@ -171,7 +177,8 @@ class NRLMSISE00Atmosphere : public AtmosphereModel { * \param GasComponentProperties . * \return void. */ - void setGasComponentProperties(GasComponentProperties gasComponentProperties){ + void setGasComponentProperties(GasComponentProperties gasComponentProperties) + { gasComponentProperties_ = gasComponentProperties; } @@ -241,7 +248,8 @@ class NRLMSISE00Atmosphere : public AtmosphereModel { * \return speed of sound. */ double getSpeedOfSound(const double altitude, const double longitude, - const double latitude, const double time) { + const double latitude, const double time) + { computeProperties(altitude, longitude, latitude, time); return speedOfSound_; } @@ -256,7 +264,8 @@ class NRLMSISE00Atmosphere : public AtmosphereModel { * \return mean free path. */ double getMeanFreePath(const double altitude, const double longitude, - const double latitude, const double time) { + const double latitude, const double time) + { computeProperties(altitude, longitude, latitude, time); return meanFreePath_; } @@ -271,7 +280,8 @@ class NRLMSISE00Atmosphere : public AtmosphereModel { * \return mean molar mass. */ double getMeanMolarMass(const double altitude, const double longitude, - const double latitude, const double time) { + const double latitude, const double time) + { computeProperties(altitude, longitude, latitude, time); return meanMolarMass_; } @@ -286,7 +296,8 @@ class NRLMSISE00Atmosphere : public AtmosphereModel { * \return number densities of gas components */ std::vector getNumberDensities(const double altitude, const double longitude, - const double latitude, const double time) { + const double latitude, const double time) + { computeProperties(altitude, longitude, latitude, time); return numberDensities_; } @@ -301,7 +312,8 @@ class NRLMSISE00Atmosphere : public AtmosphereModel { * \return average number density. */ double getAverageNumberDensity(const double altitude, const double longitude, - const double latitude, const double time) { + const double latitude, const double time) + { computeProperties(altitude, longitude, latitude, time); return averageNumberDensity_; } @@ -316,7 +328,8 @@ class NRLMSISE00Atmosphere : public AtmosphereModel { * \return weighted average collision diameter. */ double getWeightedAverageCollisionDiameter(const double altitude, const double longitude, - const double latitude, const double time) { + const double latitude, const double time) + { computeProperties(altitude, longitude, latitude, time); return weightedAverageCollisionDiameter_; } @@ -351,6 +364,9 @@ class NRLMSISE00Atmosphere : public AtmosphereModel { //! Shared pointer to solar activity function NRLMSISE00InputFunction nrlmsise00InputFunction_; + //! Use the ideal gas law for the computation of the pressure. + bool useIdealGasLaw_; + //! Current key hash size_t hashKey_; From 1822f057c9c3fea395a44dcd5cb2d0c201f425d1 Mon Sep 17 00:00:00 2001 From: Reneh107 Date: Thu, 5 May 2016 09:29:00 +0200 Subject: [PATCH 4/9] Minor adjustments + Add tests: speed of sound , mean molar mass, mean free path --- .../Astrodynamics/Aerodynamics/CMakeLists.txt | 2 +- .../unitTestNRLMSISE00Atmosphere.cpp | 89 +++++++++++++++++++ .../Aerodynamics/atmosphereModel.h | 10 +-- .../Aerodynamics/nrlmsise00Atmosphere.h | 10 +-- 4 files changed, 100 insertions(+), 11 deletions(-) diff --git a/Tudat/Astrodynamics/Aerodynamics/CMakeLists.txt b/Tudat/Astrodynamics/Aerodynamics/CMakeLists.txt index 12b548d155..2f3a5322c9 100644 --- a/Tudat/Astrodynamics/Aerodynamics/CMakeLists.txt +++ b/Tudat/Astrodynamics/Aerodynamics/CMakeLists.txt @@ -98,5 +98,5 @@ target_link_libraries(test_TabulatedAtmosphere tudat_aerodynamics tudat_interpol if(USE_NRLMSISE00) add_executable(test_NRLMSISE00Atmosphere "${SRCROOT}${AERODYNAMICSDIR}/UnitTests/unitTestNRLMSISE00Atmosphere.cpp") setup_custom_test_program(test_NRLMSISE00Atmosphere "${SRCROOT}${AERODYNAMICSDIR}") - target_link_libraries(test_NRLMSISE00Atmosphere tudat_aerodynamics tudat_basic_mathematics nrlmsise00 ${Boost_LIBRARIES}) + target_link_libraries(test_NRLMSISE00Atmosphere tudat_aerodynamics tudat_interpolators tudat_basic_mathematics nrlmsise00 tudat_input_output ${Boost_LIBRARIES}) endif( ) diff --git a/Tudat/Astrodynamics/Aerodynamics/UnitTests/unitTestNRLMSISE00Atmosphere.cpp b/Tudat/Astrodynamics/Aerodynamics/UnitTests/unitTestNRLMSISE00Atmosphere.cpp index ad49582d8a..9899bfd28c 100644 --- a/Tudat/Astrodynamics/Aerodynamics/UnitTests/unitTestNRLMSISE00Atmosphere.cpp +++ b/Tudat/Astrodynamics/Aerodynamics/UnitTests/unitTestNRLMSISE00Atmosphere.cpp @@ -962,6 +962,95 @@ BOOST_AUTO_TEST_CASE( testNRLMSISE00AtmosphereTest17 ) } } +BOOST_AUTO_TEST_CASE( testSpeedOfSound ) +{ + // Construct model with default properties + NRLMSISE00Atmosphere model( boost::bind( &function, _1, _2, _3, _4, false, false ) ); + + double altitude = 0.0 ; + double longitude = 0.0 ; + double latitude = 0.0 ; + double time = 0.0 ; + + // Check speed of sound : (https://en.wikipedia.org/wiki/Speed_of_sound (at 20 deg celcius) ) + // 340.9 at sealevel (anderson, fundamentals of aerodynamics) + BOOST_CHECK_CLOSE_FRACTION( model.getSpeedOfSound(altitude,longitude,latitude,time) , 343.21 , 2E-2 ); // < 2% + BOOST_CHECK_CLOSE_FRACTION( model.getSpeedOfSound(altitude,longitude,latitude,time) , 340.9 , 1.8E-2 ); // < 1.8% + + // US Standard Atmosphere Model 1976 + std::string atmosphereTableFile = tudat::input_output::getTudatRootPath( ) + "/External/AtmosphereTables/" + + "USSA1976Until100kmPer100mUntil1000kmPer1000m.dat" ; + tudat::aerodynamics::TabulatedAtmosphere US76model(atmosphereTableFile); + + // Test wist US76 table (NRLMSISE should provide approximately the same results) + BOOST_CHECK_CLOSE_FRACTION( model.getSpeedOfSound(altitude,longitude,latitude,time) , US76model.getSpeedOfSound(altitude) , 2E-2 ); // < 2% + + altitude = 10.0E3 ; + BOOST_CHECK_CLOSE_FRACTION( model.getSpeedOfSound(altitude,longitude,latitude,time) , US76model.getSpeedOfSound(altitude) , 3.5E-2 ); // < 3.5% + + altitude = 40.0E3 ; + BOOST_CHECK_CLOSE_FRACTION( model.getSpeedOfSound(altitude,longitude,latitude,time) , US76model.getSpeedOfSound(altitude) , 6.5E-2 ); // < 6.5% + + altitude = 80.0E3 ; + BOOST_CHECK_CLOSE_FRACTION( model.getSpeedOfSound(altitude,longitude,latitude,time) , US76model.getSpeedOfSound(altitude) , 5E-2 ); // < 5% + + // Verify calculation + double speedOfSound = std::sqrt( 1.4 * 8.3144598 * model.getTemperature(altitude,longitude,latitude,time) / model.getMeanMolarMass(altitude,longitude,latitude,time) ) ; + BOOST_CHECK_CLOSE_FRACTION( speedOfSound , model.getSpeedOfSound(altitude,longitude,latitude,time), 1E-15 ); +} + +BOOST_AUTO_TEST_CASE( testMolarMass ) +{ + // Construct model with default properties + NRLMSISE00Atmosphere model( boost::bind( &function, _1, _2, _3, _4, false, false ) ); + + double altitude = 0.0 ; + double longitude = 0.0 ; + double latitude = 0.0 ; + double time = 0.0 ; + + // Check using textbook example: Dynamics of atmospheric re-entry, Regan , 1993 + BOOST_CHECK_CLOSE_FRACTION( model.getMeanMolarMass(altitude,longitude,latitude,time), 28.96643E-3 , 3E-4 ); // < 0.03% + + altitude = 50.0E3; + BOOST_CHECK_CLOSE_FRACTION( model.getMeanMolarMass(altitude,longitude,latitude,time), 28.96643E-3 , 3E-4 ); // < 0.03% + + altitude = 200.0E3; + BOOST_CHECK_CLOSE_FRACTION( model.getMeanMolarMass(altitude,longitude,latitude,time), 21.0E-3 , 3E-2 ); // < 3% + + altitude = 300.0E3; + BOOST_CHECK_CLOSE_FRACTION( model.getMeanMolarMass(altitude,longitude,latitude,time), 18.0E-3 , 3E-2 ); // < 3% +} + +BOOST_AUTO_TEST_CASE( testMeanFreePath ) +{ + // Construct model with default properties + NRLMSISE00Atmosphere model( boost::bind( &function, _1, _2, _3, _4, false, false ) ); + + double altitude = 20.0E3 ; + double longitude = 0.5 ; + double latitude = 0.2 ; + double time = 1.0E5 ; + + // Verify calculation + double meanFreePath = std::sqrt(2.0) * PI * std::pow( model.getWeightedAverageCollisionDiameter(altitude,longitude,latitude,time) , 2.0) + * model.getAverageNumberDensity(altitude,longitude,latitude,time) ; + meanFreePath = std::pow( meanFreePath , (-1.0) ); + + BOOST_CHECK_CLOSE_FRACTION( model.getMeanFreePath(altitude,longitude,latitude,time) , meanFreePath , 1.0E-15 ); + + // Verify using data - Logarithmic plot: physics of the Earth's space environment, Gerd W. Prolls (page 29) + // Test using approximate values obtained from figure + altitude = 0.0E3 ; + BOOST_CHECK_CLOSE_FRACTION( model.getMeanFreePath(altitude,longitude,latitude,time) , 1.0E-7 , 4.0 ); + + altitude = 60.0E3 ; + BOOST_CHECK_CLOSE_FRACTION( model.getMeanFreePath(altitude,longitude,latitude,time) , 1.0E-3 , 0.8 ); + + altitude = 300.0E3 ; + BOOST_CHECK_CLOSE_FRACTION( model.getMeanFreePath(altitude,longitude,latitude,time) , 1.0E4 , 0.6 ); +} + BOOST_AUTO_TEST_SUITE_END( ) } // namespace unit_tests diff --git a/Tudat/Astrodynamics/Aerodynamics/atmosphereModel.h b/Tudat/Astrodynamics/Aerodynamics/atmosphereModel.h index 204fa2ed26..91d5d0ef13 100644 --- a/Tudat/Astrodynamics/Aerodynamics/atmosphereModel.h +++ b/Tudat/Astrodynamics/Aerodynamics/atmosphereModel.h @@ -42,7 +42,7 @@ #define TUDAT_ATMOSPHERE_MODEL_H #include -#include +#include "tudat/Mathematics/BasicMathematics/mathematicalConstants.h" namespace tudat { @@ -202,7 +202,7 @@ class AtmosphereModel */ virtual double getMeanFreePath(const double altitude, const double longitude, const double latitude, const double time){ - return std::numeric_limits::quiet_NaN(); + return TUDAT_NAN; } //! Get number densities. @@ -233,7 +233,7 @@ class AtmosphereModel */ virtual double getMeanMolarMass(const double altitude, const double longitude, const double latitude, const double time){ - return std::numeric_limits::quiet_NaN(); + return TUDAT_NAN; } //! Get average number density @@ -248,7 +248,7 @@ class AtmosphereModel */ virtual double getAverageNumberDensity(const double altitude, const double longitude, const double latitude, const double time){ - return std::numeric_limits::quiet_NaN(); + return TUDAT_NAN; } //! Get weighted average collision diameter @@ -263,7 +263,7 @@ class AtmosphereModel */ virtual double getWeightedAverageCollisionDiameter(const double altitude, const double longitude, const double latitude, const double time){ - return std::numeric_limits::quiet_NaN(); + return TUDAT_NAN; } //! Reset Hash key diff --git a/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.h b/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.h index 6233b43be7..299e2587da 100644 --- a/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.h +++ b/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.h @@ -138,7 +138,7 @@ class NRLMSISE00Atmosphere : public AtmosphereModel { * \param useIdealGasLaw Variable denoting whether to use the ideal gas law for computation of pressure. */ NRLMSISE00Atmosphere(const NRLMSISE00InputFunction nrlmsise00InputFunction, - bool useIdealGasLaw = false) + const bool useIdealGasLaw = true) : nrlmsise00InputFunction_(nrlmsise00InputFunction) { resetHashKey(); @@ -158,9 +158,9 @@ class NRLMSISE00Atmosphere : public AtmosphereModel { * the molecule collision diameters and the molar mass. */ NRLMSISE00Atmosphere(const NRLMSISE00InputFunction nrlmsise00InputFunction, - double specificHeatRatio, - GasComponentProperties gasProperties, - bool useIdealGasLaw = false) + const double specificHeatRatio, + const GasComponentProperties gasProperties, + const bool useIdealGasLaw = true) : nrlmsise00InputFunction_(nrlmsise00InputFunction) { resetHashKey(); @@ -177,7 +177,7 @@ class NRLMSISE00Atmosphere : public AtmosphereModel { * \param GasComponentProperties . * \return void. */ - void setGasComponentProperties(GasComponentProperties gasComponentProperties) + void setGasComponentProperties(const GasComponentProperties gasComponentProperties) { gasComponentProperties_ = gasComponentProperties; } From f5d142a8199dc43c2449bdc4e3ee945462604876 Mon Sep 17 00:00:00 2001 From: Reneh107 Date: Mon, 8 Aug 2016 12:19:00 +0200 Subject: [PATCH 5/9] Add nrlmsise input functions --- .../Astrodynamics/Aerodynamics/CMakeLists.txt | 6 +- .../Aerodynamics/nrlmsise00InputFunctions.cpp | 96 +++++++++++++++++++ .../Aerodynamics/nrlmsise00InputFunctions.h | 63 ++++++++++++ 3 files changed, 163 insertions(+), 2 deletions(-) create mode 100644 Tudat/Astrodynamics/Aerodynamics/nrlmsise00InputFunctions.cpp create mode 100644 Tudat/Astrodynamics/Aerodynamics/nrlmsise00InputFunctions.h diff --git a/Tudat/Astrodynamics/Aerodynamics/CMakeLists.txt b/Tudat/Astrodynamics/Aerodynamics/CMakeLists.txt index 2f3a5322c9..715bc753a6 100644 --- a/Tudat/Astrodynamics/Aerodynamics/CMakeLists.txt +++ b/Tudat/Astrodynamics/Aerodynamics/CMakeLists.txt @@ -65,9 +65,11 @@ set(AERODYNAMICS_HEADERS if(USE_NRLMSISE00) set(AERODYNAMICS_SOURCES "${AERODYNAMICS_SOURCES}" - "${SRCROOT}${AERODYNAMICSDIR}/nrlmsise00Atmosphere.cpp") + "${SRCROOT}${AERODYNAMICSDIR}/nrlmsise00Atmosphere.cpp" + "${SRCROOT}${AERODYNAMICSDIR}/nrlmsise00InputFunctions.cpp") set(AERODYNAMICS_HEADERS "${AERODYNAMICS_HEADERS}" - "${SRCROOT}${AERODYNAMICSDIR}/nrlmsise00Atmosphere.h") + "${SRCROOT}${AERODYNAMICSDIR}/nrlmsise00Atmosphere.h" + "${SRCROOT}${AERODYNAMICSDIR}/nrlmsise00InputFunctions.h") endif( ) # Add static libraries. diff --git a/Tudat/Astrodynamics/Aerodynamics/nrlmsise00InputFunctions.cpp b/Tudat/Astrodynamics/Aerodynamics/nrlmsise00InputFunctions.cpp new file mode 100644 index 0000000000..eee1f7013e --- /dev/null +++ b/Tudat/Astrodynamics/Aerodynamics/nrlmsise00InputFunctions.cpp @@ -0,0 +1,96 @@ +/* Copyright (c) 2010-2015, Delft University of Technology + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without modification, are + * permitted provided that the following conditions are met: + * - Redistributions of source code must retain the above copyright notice, this list of + * conditions and the following disclaimer. + * - Redistributions in binary form must reproduce the above copyright notice, this list of + * conditions and the following disclaimer in the documentation and/or other materials + * provided with the distribution. + * - Neither the name of the Delft University of Technology nor the names of its contributors + * may be used to endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS + * OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE + * COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE + * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED + * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED + * OF THE POSSIBILITY OF SUCH DAMAGE. + * + * Changelog + * + * References + * + * Notes + * + */ +#include + +#include "Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.h" +#include "Tudat/Astrodynamics/Aerodynamics/nrlmsise00InputFunctions.h" +#include "Tudat/Mathematics/BasicMathematics/mathematicalConstants.h" +#include "Tudat/Astrodynamics/BasicAstrodynamics/timeConversions.h" + + +//! Tudat library namespace. +namespace tudat +{ +namespace aerodynamics +{ + +//! NRLMSISE00Input function +NRLMSISE00Input nrlmsiseInputFunction(double altitude, double longitude, + double latitude, double time, + tudat::input_output::solar_activity::SolarActivityDataMap& solarActivityMap) { + // Declare input data class member + NRLMSISE00Input nrlmsiseInputData ; + + // Julian dates + double julianDate = tudat::basic_astrodynamics::convertSecondsSinceEpochToJulianDay( time ); + double julianDay = std::floor(julianDate-0.5)+0.5 ; + + // Move not always find, only if day changes + // Only use datamap with this day + some extra days.. + // Find solar activity data for current date + SolarActivityDataPtr solarActivity = solarActivityMap[ julianDay ] ; + if(!solarActivity){ + std::cerr << "Solar activity data could not be found for this date.." << std::endl; + } + + // Compute julian date at the first of januari + double julianDate1Jan = tudat::basic_astrodynamics::convertCalendarDateToJulianDay( + solarActivity->year, 1, 1, 0, 0, 0.0) ; + + nrlmsiseInputData.year = solarActivity->year ; // int + nrlmsiseInputData.dayOfTheYear = julianDay - julianDate1Jan + 1 ; + nrlmsiseInputData.secondOfTheDay = time - + tudat::basic_astrodynamics::convertJulianDayToSecondsSinceEpoch( julianDay, + tudat::basic_astrodynamics::JULIAN_DAY_ON_J2000 ) ; + + if( solarActivity->fluxQualifier == 1 ){ // requires adjustment + nrlmsiseInputData.f107 = solarActivity->solarRadioFlux107Adjusted ; + nrlmsiseInputData.f107a = solarActivity->centered81DaySolarRadioFlux107Adjusted ; + } + else{ // no adjustment required + nrlmsiseInputData.f107 = solarActivity->solarRadioFlux107Observed ; + nrlmsiseInputData.f107a = solarActivity->centered81DaySolarRadioFlux107Observed ; + } + nrlmsiseInputData.apDaily = solarActivity->planetaryEquivalentAmplitudeAverage ; + nrlmsiseInputData.apVector = EigenToStdVector( solarActivity->planetaryEquivalentAmplitudeVector ) ; // std vector + + // Compute local solar time + // Hrs since begin of the day at longitude 0 (GMT) + Hrs passed at current longitude +// nrlmsiseInputData.localSolarTime = 16 ; // use this for unit test + nrlmsiseInputData.localSolarTime = nrlmsiseInputData.secondOfTheDay/3600.0 + + longitude/(tudat::mathematical_constants::PI/12.0); + + return nrlmsiseInputData; +} + +} // namespace aerodynamics +} // namespace tudat diff --git a/Tudat/Astrodynamics/Aerodynamics/nrlmsise00InputFunctions.h b/Tudat/Astrodynamics/Aerodynamics/nrlmsise00InputFunctions.h new file mode 100644 index 0000000000..21a7c93255 --- /dev/null +++ b/Tudat/Astrodynamics/Aerodynamics/nrlmsise00InputFunctions.h @@ -0,0 +1,63 @@ +/*! Copyright (c) 2010-2012 Delft University of Technology. + * + * This software is protected by national and international copyright. + * Any unauthorized use, reproduction or modification is unlawful and + * will be prosecuted. Commercial and non-private application of the + * software in any form is strictly prohibited unless otherwise granted + * by the authors. + * + * The code is provided without any warranty; without even the implied + * warranty of merchantibility or fitness for a particular purpose. + * + * Changelog + * YYMMDD Author Comment + * 110224 F.M. Engelen File created. + * 110324 J. Melman Added overloaded get functions. + * 110427 F.M. Engelen Changed input parameter to altitude, longitude and latitude. + * 110629 F.M. Engelen Added predefined function. + * 110705 F.M. Engelen Changed to passing by reference. + * 110810 J. Leloux Corrected doxygen documentation. + * 151104 J. Geul Complete rewrite + * + * References + * + */ + +#ifndef TUDAT_NRLMSISE00_INPUT_FUNCTIONS_H +#define TUDAT_NRLMSISE00_INPUT_FUNCTIONS_H + +#include +#include + +#include + +#include "Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.h" +#include "Tudat/InputOutput/solarActivityData.h" + + +namespace tudat +{ + +namespace aerodynamics +{ + + +//! NRLMSISE00 Input function +/*! + * This function is used to define the input for the NRLMSISE model. + * This function reads solar activity data and defines the input using this data. + * \param altitude in m + * \param longitude in RAD + * \param latitude in RAD + * \param time in seconds since J2000 epoch + * \param solarActivityData SolarActivityData structure + * \return NRLMSISE00Input nrlmsiseInputFunction + */ +NRLMSISE00Input nrlmsiseInputFunction(double altitude, double longitude, + double latitude, double time, + tudat::input_output::solar_activity::SolarActivityDataMap& solarActivityMap); + +} // namespace aerodynamics +} // namespace tudat + +#endif // TUDAT_NRLMSISE00_ATMOSPHERE_H_ From 6553307ce07774feb0983dea54c3ddf353f5c5a1 Mon Sep 17 00:00:00 2001 From: Reneh107 Date: Mon, 8 Aug 2016 12:20:16 +0200 Subject: [PATCH 6/9] Update vector size (fix data leak) --- Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.cpp b/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.cpp index a75e08ff43..281be44c8f 100644 --- a/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.cpp +++ b/Tudat/Astrodynamics/Aerodynamics/nrlmsise00Atmosphere.cpp @@ -83,7 +83,7 @@ void NRLMSISE00Atmosphere::computeProperties(double altitude, double longitude, temperature_ = output_.t[1]; // Get number densities - numberDensities_.resize(7); + numberDensities_.resize(8); numberDensities_[0] = output_.d[0] * 1.0E6 ; // HE NUMBER DENSITY (M-3) numberDensities_[1] = output_.d[1] * 1.0E6 ; // O NUMBER DENSITY (M-3) numberDensities_[2] = output_.d[2] * 1.0E6 ; // N2 NUMBER DENSITY (M-3) From 40083e3201f0f6516d9c918965317219b94b9339 Mon Sep 17 00:00:00 2001 From: Reneh107 Date: Mon, 8 Aug 2016 12:24:56 +0200 Subject: [PATCH 7/9] base class functions update --- .../Aerodynamics/atmosphereModel.h | 64 +------------------ 1 file changed, 2 insertions(+), 62 deletions(-) diff --git a/Tudat/Astrodynamics/Aerodynamics/atmosphereModel.h b/Tudat/Astrodynamics/Aerodynamics/atmosphereModel.h index 91d5d0ef13..abf7ae3b2c 100644 --- a/Tudat/Astrodynamics/Aerodynamics/atmosphereModel.h +++ b/Tudat/Astrodynamics/Aerodynamics/atmosphereModel.h @@ -201,68 +201,8 @@ class AtmosphereModel * \return double mean free path. */ virtual double getMeanFreePath(const double altitude, const double longitude, - const double latitude, const double time){ - return TUDAT_NAN; - } - - //! Get number densities. - /*! - * This function returns a vector with number densities for all the components of the gas(air). - * The function is not implemented by default and returns null vector if the atmosphere model doesn't implement this function. - * \param altitude Altitude. - * \param longitude Longitude. - * \param latitude Latitude. - * \param time Time. - * \return std::vector number densities. - */ - virtual std::vector getNumberDensities(const double altitude, const double longitude, - const double latitude, const double time){ - std::vector null(0); - return null; - } - - //! Get mean molar mass. - /*! - * This function returns the mean molar mass at the current location and time. - * The function is not implemented by default and returns nan if the atmosphere model doesn't implement this function. - * \param altitude Altitude. - * \param longitude Longitude. - * \param latitude Latitude. - * \param time Time. - * \return double meanMolarMass. - */ - virtual double getMeanMolarMass(const double altitude, const double longitude, - const double latitude, const double time){ - return TUDAT_NAN; - } - - //! Get average number density - /*! - * This function returns the (unweighted) average number density at the current location and time. - * The function is not implemented by default and returns nan if the atmosphere model doesn't implement this function. - * \param altitude Altitude. - * \param longitude Longitude. - * \param latitude Latitude. - * \param time Time. - * \return double average number density. - */ - virtual double getAverageNumberDensity(const double altitude, const double longitude, - const double latitude, const double time){ - return TUDAT_NAN; - } - - //! Get weighted average collision diameter - /*! - * This function returns the weighted average collision diameter of the gas components at the current location and time. - * The function is not implemented by default and returns nan if the atmosphere model doesn't implement this function. - * \param altitude Altitude. - * \param longitude Longitude. - * \param latitude Latitude. - * \param time Time. - * \return double weighted average collision diameter. - */ - virtual double getWeightedAverageCollisionDiameter(const double altitude, const double longitude, - const double latitude, const double time){ + const double latitude, const double time) + { return TUDAT_NAN; } From 10be9329e3f1c9f00e188c8bd71cb8dc601403e0 Mon Sep 17 00:00:00 2001 From: Reneh107 Date: Mon, 8 Aug 2016 13:00:45 +0200 Subject: [PATCH 8/9] Added unit tests for nrlmsise input functions --- .../Astrodynamics/Aerodynamics/CMakeLists.txt | 2 +- .../UnitTests/swAtmosTestNoAdjust.txt | 49 +++++++++ .../UnitTests/swAtmosTestWithAdjust.txt | 48 +++++++++ .../unitTestNRLMSISE00Atmosphere.cpp | 99 +++++++++++++++++++ .../Aerodynamics/nrlmsise00InputFunctions.cpp | 32 +++++- .../Aerodynamics/nrlmsise00InputFunctions.h | 3 +- 6 files changed, 226 insertions(+), 7 deletions(-) create mode 100644 Tudat/Astrodynamics/Aerodynamics/UnitTests/swAtmosTestNoAdjust.txt create mode 100644 Tudat/Astrodynamics/Aerodynamics/UnitTests/swAtmosTestWithAdjust.txt diff --git a/Tudat/Astrodynamics/Aerodynamics/CMakeLists.txt b/Tudat/Astrodynamics/Aerodynamics/CMakeLists.txt index 715bc753a6..1e469387e8 100644 --- a/Tudat/Astrodynamics/Aerodynamics/CMakeLists.txt +++ b/Tudat/Astrodynamics/Aerodynamics/CMakeLists.txt @@ -100,5 +100,5 @@ target_link_libraries(test_TabulatedAtmosphere tudat_aerodynamics tudat_interpol if(USE_NRLMSISE00) add_executable(test_NRLMSISE00Atmosphere "${SRCROOT}${AERODYNAMICSDIR}/UnitTests/unitTestNRLMSISE00Atmosphere.cpp") setup_custom_test_program(test_NRLMSISE00Atmosphere "${SRCROOT}${AERODYNAMICSDIR}") - target_link_libraries(test_NRLMSISE00Atmosphere tudat_aerodynamics tudat_interpolators tudat_basic_mathematics nrlmsise00 tudat_input_output ${Boost_LIBRARIES}) + target_link_libraries(test_NRLMSISE00Atmosphere tudat_aerodynamics tudat_interpolators tudat_basic_mathematics nrlmsise00 tudat_input_output tudat_basic_astrodynamics ${Boost_LIBRARIES}) endif( ) diff --git a/Tudat/Astrodynamics/Aerodynamics/UnitTests/swAtmosTestNoAdjust.txt b/Tudat/Astrodynamics/Aerodynamics/UnitTests/swAtmosTestNoAdjust.txt new file mode 100644 index 0000000000..fd5a7555ee --- /dev/null +++ b/Tudat/Astrodynamics/Aerodynamics/UnitTests/swAtmosTestNoAdjust.txt @@ -0,0 +1,49 @@ +DATATYPE CssiSpaceWeather +VERSION 1.2 +UPDATED 2016 Mar 24 18:35:11 UTC +# -------------------------------------------------------------------------------------------------------------------------------- +# ANALYTICAL GRAPHICS, INC. +# CENTER FOR SPACE STANDARDS & INNOVATION +# SPACE WEATHER DATA +# -------------------------------------------------------------------------------------------------------------------------------- +# +# See http://celestrak.com/SpaceData/SpaceWx-format.asp for format details. +# +# FORMAT(I4,I3,I3,I5,I3,8I3,I4,8I4,I4,F4.1,I2,I4,F6.1,I2,5F6.1) +# -------------------------------------------------------------------------------------------------------------------------------- +# Adj Adj Adj Obs Obs Obs +# yy mm dd BSRN ND Kp Kp Kp Kp Kp Kp Kp Kp Sum Ap Ap Ap Ap Ap Ap Ap Ap Avg Cp C9 ISN F10.7 Q Ctr81 Lst81 F10.7 Ctr81 Lst81 +# -------------------------------------------------------------------------------------------------------------------------------- +# +NUM_OBSERVED_POINTS 21360 +BEGIN OBSERVED +1957 10 01 1700 19 43 40 30 20 37 23 43 37 273 32 27 15 7 22 9 32 22 21 1.1 5 236 268.0 0 265.2 230.6 269.3 266.6 230.9 +1957 10 02 1700 20 37 37 17 17 27 23 17 30 203 22 22 6 6 12 9 6 15 12 0.7 3 234 252.0 0 266.0 231.4 253.3 267.4 231.7 +1957 10 03 1700 21 27 20 13 33 37 47 43 30 250 12 7 5 18 22 39 32 15 19 1.0 5 242 265.0 0 266.7 232.3 266.3 268.1 232.7 +1957 10 04 1700 22 30 30 23 27 23 27 30 27 217 15 15 9 12 9 12 15 12 12 0.7 3 217 237.0 0 267.4 232.9 238.2 268.8 233.3 +1957 10 05 1700 23 30 30 17 23 20 27 27 20 193 15 15 6 9 7 12 12 7 10 0.6 3 219 245.0 0 267.8 233.5 246.2 269.3 233.9 +1957 10 06 1700 24 17 3 10 7 0 0 3 3 43 6 2 4 3 0 0 2 2 2 0.0 0 227 250.0 0 268.1 233.9 251.2 269.6 234.3 +1957 10 07 1700 25 10 20 17 17 17 7 20 13 120 4 7 6 6 6 3 7 5 6 0.2 1 234 253.0 0 268.4 234.2 254.3 269.9 234.6 +1957 10 08 1700 26 17 0 10 7 3 3 3 3 47 6 0 4 3 2 2 2 2 3 0.0 0 244 261.0 0 268.2 234.7 262.6 269.8 235.2 +1957 10 09 1700 27 3 17 23 17 23 23 20 23 150 2 6 9 6 9 9 7 9 7 0.4 2 267 275.0 0 268.3 235.2 276.6 270.0 235.7 +1957 10 10 1701 1 30 40 40 30 23 30 37 27 257 15 27 27 15 9 15 22 12 18 1.0 5 264 275.0 0 268.2 235.6 276.6 269.8 236.1 +1957 10 11 1701 2 33 43 33 43 37 30 23 23 267 18 32 18 32 22 15 9 9 19 1.0 5 232 278.0 0 267.8 236.1 279.7 269.5 236.7 +1957 10 12 1701 3 23 30 27 37 27 27 17 33 220 9 15 12 22 12 12 6 18 13 0.8 4 236 284.0 0 267.8 236.6 285.7 269.5 237.1 +1957 10 13 1701 4 50 53 30 23 27 23 43 40 290 48 56 15 9 12 9 32 27 26 1.2 6 244 281.0 0 268.1 236.9 282.7 269.9 237.5 +1957 10 14 1701 5 43 63 57 50 43 53 40 50 400 32 94 67 48 32 56 27 48 50 1.6 7 232 289.0 0 268.2 237.4 290.7 270.0 238.0 +1957 10 15 1701 6 33 20 17 23 30 27 33 20 203 18 7 6 9 15 12 18 7 12 0.7 3 264 277.0 0 268.7 238.0 278.9 270.5 238.6 +1957 10 16 1701 7 13 3 3 7 13 13 17 3 73 5 2 2 3 5 5 6 2 4 0.1 0 268 289.0 0 269.1 238.7 291.0 271.0 239.3 +1957 10 17 1701 8 13 13 13 30 27 20 10 13 140 5 5 5 15 12 7 4 5 7 0.4 2 251 272.0 0 269.5 239.4 273.9 271.5 240.1 +1957 10 18 1701 9 17 13 17 13 7 17 27 13 123 6 5 6 5 3 6 12 5 6 0.3 1 222 294.0 0 270.0 240.6 296.1 271.9 241.3 +1957 10 19 1701 10 7 13 27 33 17 23 27 17 163 3 5 12 18 6 9 12 6 9 0.5 2 217 291.0 0 270.3 241.8 293.0 272.3 242.5 +1957 10 20 1701 11 20 23 13 20 20 30 27 40 193 7 9 5 7 7 15 12 27 11 0.6 3 230 303.0 0 270.5 243.2 305.1 272.5 243.9 +1957 10 21 1701 12 33 27 27 17 27 40 37 67 273 18 12 12 6 12 27 22 111 28 1.2 6 237 282.0 0 270.9 244.1 284.0 273.0 244.9 +1957 10 22 1701 13 43 30 37 37 23 33 27 37 267 32 15 22 22 9 18 12 22 19 1.0 5 241 274.0 0 271.2 245.1 276.2 273.4 245.9 +2030 06 21 3892 1 1 1 1 1 1 1 1 1 8 0 0 0 0 0 0 0 0 4 220.0 0 220.0 210.0 150.0 150.0 180.0 +2030 06 22 3892 1 1 1 1 1 1 1 1 1 8 0 0 0 0 0 0 0 0 4 100.0 0 100.0 100.0 100.0 100.0 100.0 +1957 10 24 1701 15 20 17 20 23 27 17 37 10 170 7 6 7 9 12 6 22 4 9 0.5 2 276 280.0 0 271.3 247.3 282.2 273.6 248.2 +1957 10 25 1701 16 20 17 20 30 20 13 20 23 163 7 6 7 15 7 5 7 9 8 0.4 2 240 259.0 0 271.6 248.2 261.1 273.8 249.1 +1957 10 26 1701 17 30 20 13 23 20 13 17 20 157 15 7 5 9 7 5 6 7 8 0.4 2 293 266.0 0 271.8 249.2 268.1 274.1 250.1 +1957 10 27 1701 18 20 20 20 27 40 33 20 23 203 7 7 7 12 27 18 7 9 12 0.7 3 280 296.0 0 271.9 250.6 298.4 274.2 251.5 +END OBSERVED + diff --git a/Tudat/Astrodynamics/Aerodynamics/UnitTests/swAtmosTestWithAdjust.txt b/Tudat/Astrodynamics/Aerodynamics/UnitTests/swAtmosTestWithAdjust.txt new file mode 100644 index 0000000000..b37938b8e9 --- /dev/null +++ b/Tudat/Astrodynamics/Aerodynamics/UnitTests/swAtmosTestWithAdjust.txt @@ -0,0 +1,48 @@ +DATATYPE CssiSpaceWeather +VERSION 1.2 +UPDATED 2016 Mar 24 18:35:11 UTC +# -------------------------------------------------------------------------------------------------------------------------------- +# ANALYTICAL GRAPHICS, INC. +# CENTER FOR SPACE STANDARDS & INNOVATION +# SPACE WEATHER DATA +# -------------------------------------------------------------------------------------------------------------------------------- +# +# See http://celestrak.com/SpaceData/SpaceWx-format.asp for format details. +# +# FORMAT(I4,I3,I3,I5,I3,8I3,I4,8I4,I4,F4.1,I2,I4,F6.1,I2,5F6.1) +# -------------------------------------------------------------------------------------------------------------------------------- +# Adj Adj Adj Obs Obs Obs +# yy mm dd BSRN ND Kp Kp Kp Kp Kp Kp Kp Kp Sum Ap Ap Ap Ap Ap Ap Ap Ap Avg Cp C9 ISN F10.7 Q Ctr81 Lst81 F10.7 Ctr81 Lst81 +# -------------------------------------------------------------------------------------------------------------------------------- +# +NUM_OBSERVED_POINTS 21360 +BEGIN OBSERVED +1957 10 01 1700 19 43 40 30 20 37 23 43 37 273 32 27 15 7 22 9 32 22 21 1.1 5 236 268.0 0 265.2 230.6 269.3 266.6 230.9 +1957 10 02 1700 20 37 37 17 17 27 23 17 30 203 22 22 6 6 12 9 6 15 12 0.7 3 234 252.0 0 266.0 231.4 253.3 267.4 231.7 +1957 10 03 1700 21 27 20 13 33 37 47 43 30 250 12 7 5 18 22 39 32 15 19 1.0 5 242 265.0 0 266.7 232.3 266.3 268.1 232.7 +1957 10 04 1700 22 30 30 23 27 23 27 30 27 217 15 15 9 12 9 12 15 12 12 0.7 3 217 237.0 0 267.4 232.9 238.2 268.8 233.3 +1957 10 05 1700 23 30 30 17 23 20 27 27 20 193 15 15 6 9 7 12 12 7 10 0.6 3 219 245.0 0 267.8 233.5 246.2 269.3 233.9 +1957 10 06 1700 24 17 3 10 7 0 0 3 3 43 6 2 4 3 0 0 2 2 2 0.0 0 227 250.0 0 268.1 233.9 251.2 269.6 234.3 +1957 10 07 1700 25 10 20 17 17 17 7 20 13 120 4 7 6 6 6 3 7 5 6 0.2 1 234 253.0 0 268.4 234.2 254.3 269.9 234.6 +1957 10 08 1700 26 17 0 10 7 3 3 3 3 47 6 0 4 3 2 2 2 2 3 0.0 0 244 261.0 0 268.2 234.7 262.6 269.8 235.2 +1957 10 09 1700 27 3 17 23 17 23 23 20 23 150 2 6 9 6 9 9 7 9 7 0.4 2 267 275.0 0 268.3 235.2 276.6 270.0 235.7 +1957 10 10 1701 1 30 40 40 30 23 30 37 27 257 15 27 27 15 9 15 22 12 18 1.0 5 264 275.0 0 268.2 235.6 276.6 269.8 236.1 +1957 10 11 1701 2 33 43 33 43 37 30 23 23 267 18 32 18 32 22 15 9 9 19 1.0 5 232 278.0 0 267.8 236.1 279.7 269.5 236.7 +1957 10 12 1701 3 23 30 27 37 27 27 17 33 220 9 15 12 22 12 12 6 18 13 0.8 4 236 284.0 0 267.8 236.6 285.7 269.5 237.1 +1957 10 13 1701 4 50 53 30 23 27 23 43 40 290 48 56 15 9 12 9 32 27 26 1.2 6 244 281.0 0 268.1 236.9 282.7 269.9 237.5 +1957 10 14 1701 5 43 63 57 50 43 53 40 50 400 32 94 67 48 32 56 27 48 50 1.6 7 232 289.0 0 268.2 237.4 290.7 270.0 238.0 +1957 10 15 1701 6 33 20 17 23 30 27 33 20 203 18 7 6 9 15 12 18 7 12 0.7 3 264 277.0 0 268.7 238.0 278.9 270.5 238.6 +1957 10 16 1701 7 13 3 3 7 13 13 17 3 73 5 2 2 3 5 5 6 2 4 0.1 0 268 289.0 0 269.1 238.7 291.0 271.0 239.3 +1957 10 17 1701 8 13 13 13 30 27 20 10 13 140 5 5 5 15 12 7 4 5 7 0.4 2 251 272.0 0 269.5 239.4 273.9 271.5 240.1 +1957 10 18 1701 9 17 13 17 13 7 17 27 13 123 6 5 6 5 3 6 12 5 6 0.3 1 222 294.0 0 270.0 240.6 296.1 271.9 241.3 +1957 10 19 1701 10 7 13 27 33 17 23 27 17 163 3 5 12 18 6 9 12 6 9 0.5 2 217 291.0 0 270.3 241.8 293.0 272.3 242.5 +1957 10 20 1701 11 20 23 13 20 20 30 27 40 193 7 9 5 7 7 15 12 27 11 0.6 3 230 303.0 0 270.5 243.2 305.1 272.5 243.9 +1957 10 21 1701 12 33 27 27 17 27 40 37 67 273 18 12 12 6 12 27 22 111 28 1.2 6 237 282.0 0 270.9 244.1 284.0 273.0 244.9 +1957 10 22 1701 13 43 30 37 37 23 33 27 37 267 32 15 22 22 9 18 12 22 19 1.0 5 241 274.0 0 271.2 245.1 276.2 273.4 245.9 +2030 06 21 3892 1 1 1 1 1 1 1 1 1 8 0 0 0 0 0 0 0 0 4 150.0 1 150.0 210.0 250.0 250.0 180.0 +1957 10 24 1701 15 20 17 20 23 27 17 37 10 170 7 6 7 9 12 6 22 4 9 0.5 2 276 280.0 0 271.3 247.3 282.2 273.6 248.2 +1957 10 25 1701 16 20 17 20 30 20 13 20 23 163 7 6 7 15 7 5 7 9 8 0.4 2 240 259.0 0 271.6 248.2 261.1 273.8 249.1 +1957 10 26 1701 17 30 20 13 23 20 13 17 20 157 15 7 5 9 7 5 6 7 8 0.4 2 293 266.0 0 271.8 249.2 268.1 274.1 250.1 +1957 10 27 1701 18 20 20 20 27 40 33 20 23 203 7 7 7 12 27 18 7 9 12 0.7 3 280 296.0 0 271.9 250.6 298.4 274.2 251.5 +END OBSERVED + diff --git a/Tudat/Astrodynamics/Aerodynamics/UnitTests/unitTestNRLMSISE00Atmosphere.cpp b/Tudat/Astrodynamics/Aerodynamics/UnitTests/unitTestNRLMSISE00Atmosphere.cpp index 9899bfd28c..2a1b71de5a 100644 --- a/Tudat/Astrodynamics/Aerodynamics/UnitTests/unitTestNRLMSISE00Atmosphere.cpp +++ b/Tudat/Astrodynamics/Aerodynamics/UnitTests/unitTestNRLMSISE00Atmosphere.cpp @@ -57,6 +57,8 @@ #include "Tudat/Astrodynamics/Aerodynamics/tabulatedAtmosphere.h" #include "Tudat/InputOutput/basicInputOutput.h" +#include "Tudat/Astrodynamics/Aerodynamics/nrlmsise00InputFunctions.h" + #include "tudat/Mathematics/BasicMathematics/mathematicalConstants.h" @@ -962,6 +964,8 @@ BOOST_AUTO_TEST_CASE( testNRLMSISE00AtmosphereTest17 ) } } +//! Perform NRLMSISE-00 test 18 +// Check computation of speed of sound BOOST_AUTO_TEST_CASE( testSpeedOfSound ) { // Construct model with default properties @@ -999,6 +1003,8 @@ BOOST_AUTO_TEST_CASE( testSpeedOfSound ) BOOST_CHECK_CLOSE_FRACTION( speedOfSound , model.getSpeedOfSound(altitude,longitude,latitude,time), 1E-15 ); } +//! Perform NRLMSISE-00 test 19 +// Check computation of molar mass BOOST_AUTO_TEST_CASE( testMolarMass ) { // Construct model with default properties @@ -1022,6 +1028,8 @@ BOOST_AUTO_TEST_CASE( testMolarMass ) BOOST_CHECK_CLOSE_FRACTION( model.getMeanMolarMass(altitude,longitude,latitude,time), 18.0E-3 , 3E-2 ); // < 3% } +//! Perform NRLMSISE-00 test 20 +// Check computation mean free path BOOST_AUTO_TEST_CASE( testMeanFreePath ) { // Construct model with default properties @@ -1051,6 +1059,97 @@ BOOST_AUTO_TEST_CASE( testMeanFreePath ) BOOST_CHECK_CLOSE_FRACTION( model.getMeanFreePath(altitude,longitude,latitude,time) , 1.0E4 , 0.6 ); } + +//! Perform NRLMSISE-00 test 21 - Test input function +// Corresponds to NRLMSISE test 1 but space weather file is used to set solar activity data. +// No adjustment values in space weather file are used. +BOOST_AUTO_TEST_CASE( test_nrlmise_InputFunction_no_adjustment ) +{ + +// NRLMSISE00Input gen_data(2030, 172, 29000.0, 16.0, 150.0, 150.0, 4.0); + // DD-MM-YY = 21-06-2030 , Hrs-Min-Sec = 8-3-20 + double julianDate = tudat::basic_astrodynamics::convertCalendarDateToJulianDay(2030,6,21,8,3,20) ; + double time = tudat::basic_astrodynamics::convertJulianDayToSecondsSinceEpoch( + julianDate , tudat::basic_astrodynamics::JULIAN_DAY_ON_J2000) ; + + // std::vector< double > gen_input = boost::assign::list_of(400.0)(-70.0)(60.0)(0.0); + double altitude = 400.0E3 ; // km + double longitude = -70.0 * PI / 180.0 ; + double latitude = 60.0 * PI / 180.0 ; + + // find space weather file + std::string cppPath( __FILE__ ); + std::string folder = cppPath.substr( 0, cppPath.find_last_of("/\\")+1); + std::string spaceWeatherFilePath = folder + "swAtmosTestNoAdjust.txt"; + + tudat::input_output::solar_activity::SolarActivityDataMap solarActivityData = + tudat::input_output::solar_activity::readSolarActivityData(spaceWeatherFilePath) ; + + // Create atmosphere model using NRLMISE00 input function + // Local solar time is set to 16.0 to correspond with testing data! + boost::function< tudat::aerodynamics::NRLMSISE00Input (double,double,double,double) > inputFunction = + boost::bind(&tudat::aerodynamics::nrlmsiseInputFunction,_1,_2,_3,_4, solarActivityData , true , 16.0 ); + + // Create Pointer to NRLMSISE model + tudat::aerodynamics::NRLMSISE00Atmosphere atmosphereModel( inputFunction ); + + // Verification data + std::vector verificationData = boost::assign::list_of + (6.665176904952E+05)(1.138805559752E+08)(1.998210925573E+07) + (4.022763585713E+05)(3.557464994516E+03)(4.074713532757E-15) + (3.475312399717E+04)(4.095913268293E+06)(2.667273209336E+04) + (1.250539943561E+03)(1.241416130019E+03); + + double computedDensity = atmosphereModel.getDensity( altitude , longitude, latitude , time ) ; + + BOOST_CHECK_CLOSE_FRACTION(verificationData[5]*1000 , computedDensity , 1E-11); +} + +//! Perform NRLMSISE-00 test 22 - Test input function +// Corresponds to NRLMSISE test 1 but space weather file is used to set solar activity data. +// Adjustment values (flux qualifier is 1) in space weather file are used. +BOOST_AUTO_TEST_CASE( test_nrlmise_InputFunction_with_adjustment ) +{ + +// NRLMSISE00Input gen_data(2030, 172, 29000.0, 16.0, 150.0, 150.0, 4.0); + // DD-MM-YY = 21-06-2030 , Hrs-Min-Sec = 8-3-20 + double julianDate = tudat::basic_astrodynamics::convertCalendarDateToJulianDay(2030,6,21,8,3,20) ; + double time = tudat::basic_astrodynamics::convertJulianDayToSecondsSinceEpoch( + julianDate , tudat::basic_astrodynamics::JULIAN_DAY_ON_J2000) ; + + // std::vector< double > gen_input = boost::assign::list_of(400.0)(-70.0)(60.0)(0.0); + double altitude = 400.0E3 ; // km + double longitude = -70.0 * PI / 180.0 ; + double latitude = 60.0 * PI / 180.0 ; + + // find space weather file + std::string cppPath( __FILE__ ); + std::string folder = cppPath.substr( 0, cppPath.find_last_of("/\\")+1); + std::string spaceWeatherFilePath = folder + "swAtmosTestWithAdjust.txt"; + + tudat::input_output::solar_activity::SolarActivityDataMap solarActivityData = + tudat::input_output::solar_activity::readSolarActivityData(spaceWeatherFilePath) ; + + // Create atmosphere model using NRLMISE00 input function + // Local solar time is set to 16.0 to correspond with testing data! + boost::function< tudat::aerodynamics::NRLMSISE00Input (double,double,double,double) > inputFunction = + boost::bind(&tudat::aerodynamics::nrlmsiseInputFunction,_1,_2,_3,_4, solarActivityData , true , 16.0 ); + + // Create Pointer to NRLMSISE model + tudat::aerodynamics::NRLMSISE00Atmosphere atmosphereModel( inputFunction ); + + // Verification data + std::vector verificationData = boost::assign::list_of + (6.665176904952E+05)(1.138805559752E+08)(1.998210925573E+07) + (4.022763585713E+05)(3.557464994516E+03)(4.074713532757E-15) + (3.475312399717E+04)(4.095913268293E+06)(2.667273209336E+04) + (1.250539943561E+03)(1.241416130019E+03); + + double computedDensity = atmosphereModel.getDensity( altitude , longitude, latitude , time ) ; + + BOOST_CHECK_CLOSE_FRACTION(verificationData[5]*1000 , computedDensity , 1E-11); +} + BOOST_AUTO_TEST_SUITE_END( ) } // namespace unit_tests diff --git a/Tudat/Astrodynamics/Aerodynamics/nrlmsise00InputFunctions.cpp b/Tudat/Astrodynamics/Aerodynamics/nrlmsise00InputFunctions.cpp index eee1f7013e..41594a22d1 100644 --- a/Tudat/Astrodynamics/Aerodynamics/nrlmsise00InputFunctions.cpp +++ b/Tudat/Astrodynamics/Aerodynamics/nrlmsise00InputFunctions.cpp @@ -37,16 +37,31 @@ #include "Tudat/Astrodynamics/BasicAstrodynamics/timeConversions.h" -//! Tudat library namespace. +// Tudat library namespace. namespace tudat { namespace aerodynamics { +//! Function to convert Eigen::VectorXd to std::vector +std::vector EigenToStdVector(Eigen::VectorXd vector) +{ + std::vector stdVector(vector.rows()) ; + for(int i = 0 ; i < vector.rows() ; i++) + { + stdVector[i] = vector(i) ; + } + return stdVector; +} + //! NRLMSISE00Input function NRLMSISE00Input nrlmsiseInputFunction(double altitude, double longitude, double latitude, double time, - tudat::input_output::solar_activity::SolarActivityDataMap& solarActivityMap) { + tudat::input_output::solar_activity::SolarActivityDataMap& solarActivityMap, + bool adjustSolarTime, + double localSolarTime) { + using namespace tudat::input_output::solar_activity; + // Declare input data class member NRLMSISE00Input nrlmsiseInputData ; @@ -57,6 +72,7 @@ NRLMSISE00Input nrlmsiseInputFunction(double altitude, double longitude, // Move not always find, only if day changes // Only use datamap with this day + some extra days.. // Find solar activity data for current date + SolarActivityDataPtr solarActivity = solarActivityMap[ julianDay ] ; if(!solarActivity){ std::cerr << "Solar activity data could not be found for this date.." << std::endl; @@ -85,9 +101,15 @@ NRLMSISE00Input nrlmsiseInputFunction(double altitude, double longitude, // Compute local solar time // Hrs since begin of the day at longitude 0 (GMT) + Hrs passed at current longitude -// nrlmsiseInputData.localSolarTime = 16 ; // use this for unit test - nrlmsiseInputData.localSolarTime = nrlmsiseInputData.secondOfTheDay/3600.0 - + longitude/(tudat::mathematical_constants::PI/12.0); + if( adjustSolarTime ) + { + nrlmsiseInputData.localSolarTime = localSolarTime; + } + else + { + nrlmsiseInputData.localSolarTime = nrlmsiseInputData.secondOfTheDay/3600.0 + + longitude/(tudat::mathematical_constants::PI/12.0); + } return nrlmsiseInputData; } diff --git a/Tudat/Astrodynamics/Aerodynamics/nrlmsise00InputFunctions.h b/Tudat/Astrodynamics/Aerodynamics/nrlmsise00InputFunctions.h index 21a7c93255..575daac265 100644 --- a/Tudat/Astrodynamics/Aerodynamics/nrlmsise00InputFunctions.h +++ b/Tudat/Astrodynamics/Aerodynamics/nrlmsise00InputFunctions.h @@ -55,7 +55,8 @@ namespace aerodynamics */ NRLMSISE00Input nrlmsiseInputFunction(double altitude, double longitude, double latitude, double time, - tudat::input_output::solar_activity::SolarActivityDataMap& solarActivityMap); + tudat::input_output::solar_activity::SolarActivityDataMap& solarActivityMap, + bool adjustSolarTime = false, double localSolarTime = 0.0); } // namespace aerodynamics } // namespace tudat From a36a76762b3ea29b140990e21b4569d14e77d771 Mon Sep 17 00:00:00 2001 From: Reneh107 Date: Mon, 8 Aug 2016 13:32:41 +0200 Subject: [PATCH 9/9] update unit test mean free path --- .../Aerodynamics/UnitTests/unitTestNRLMSISE00Atmosphere.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Tudat/Astrodynamics/Aerodynamics/UnitTests/unitTestNRLMSISE00Atmosphere.cpp b/Tudat/Astrodynamics/Aerodynamics/UnitTests/unitTestNRLMSISE00Atmosphere.cpp index 2a1b71de5a..560490b4fd 100644 --- a/Tudat/Astrodynamics/Aerodynamics/UnitTests/unitTestNRLMSISE00Atmosphere.cpp +++ b/Tudat/Astrodynamics/Aerodynamics/UnitTests/unitTestNRLMSISE00Atmosphere.cpp @@ -1050,13 +1050,13 @@ BOOST_AUTO_TEST_CASE( testMeanFreePath ) // Verify using data - Logarithmic plot: physics of the Earth's space environment, Gerd W. Prolls (page 29) // Test using approximate values obtained from figure altitude = 0.0E3 ; - BOOST_CHECK_CLOSE_FRACTION( model.getMeanFreePath(altitude,longitude,latitude,time) , 1.0E-7 , 4.0 ); + BOOST_CHECK_CLOSE_FRACTION( model.getMeanFreePath(altitude,longitude,latitude,time) , 1.0E-7 , 10.0 ); altitude = 60.0E3 ; - BOOST_CHECK_CLOSE_FRACTION( model.getMeanFreePath(altitude,longitude,latitude,time) , 1.0E-3 , 0.8 ); + BOOST_CHECK_CLOSE_FRACTION( model.getMeanFreePath(altitude,longitude,latitude,time) , 1.0E-3 , 1.2 ); altitude = 300.0E3 ; - BOOST_CHECK_CLOSE_FRACTION( model.getMeanFreePath(altitude,longitude,latitude,time) , 1.0E4 , 0.6 ); + BOOST_CHECK_CLOSE_FRACTION( model.getMeanFreePath(altitude,longitude,latitude,time) , 1.0E4 , 0.9 ); }