From abb561360954206e6c7acf57417e81bf1db4d93b Mon Sep 17 00:00:00 2001 From: Dominic Dirkx Date: Thu, 19 Oct 2023 12:13:27 +0200 Subject: [PATCH 01/25] Skeleton for Mars DTM atmosphere --- .../aerodynamics/marsDtmAtmosphereModel.h | 106 ++++++++++++++++++ src/astro/aerodynamics/CMakeLists.txt | 6 +- .../aerodynamics/marsDtmAtmosphereModel.cpp | 11 ++ tests/src/astro/aerodynamics/CMakeLists.txt | 5 + .../unitTestMarsDtmAtmosphere.cpp | 43 +++++++ 5 files changed, 169 insertions(+), 2 deletions(-) create mode 100644 include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h create mode 100644 src/astro/aerodynamics/marsDtmAtmosphereModel.cpp create mode 100644 tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp diff --git a/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h b/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h new file mode 100644 index 0000000000..dd64cfadda --- /dev/null +++ b/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h @@ -0,0 +1,106 @@ +/* Copyright (c) 2010-2019, 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 + * a copy of the license with this file. If not, please or visit: + * http://tudat.tudelft.nl/LICENSE. + * + */ + +#ifndef TUDAT_MARS_DTM_ATMOSPHERE_MODEL_H +#define TUDAT_MARS_DTM_ATMOSPHERE_MODEL_H + +#include "tudat/astro/aerodynamics/atmosphereModel.h" + +namespace tudat +{ + +namespace aerodynamics +{ + +class MarsDtmAtmosphereModel: public AtmosphereModel +{ +public: + + MarsDtmAtmosphereModel( ){ } + + //! Default destructor. + /*! + * Default destructor. + */ + virtual ~MarsDtmAtmosphereModel( ) { } + + //! Get local density. + /*! + * Returns the local density parameter of the atmosphere in kg per meter^3. + * \param altitude Altitude. + * \param longitude Longitude. + * \param latitude Latitude. + * \param time Time. + * \return Atmospheric density. + */ + virtual double getDensity( const double altitude, const double longitude, + const double latitude, const double time ) + { + + } + + //! Get local pressure. + /*! + * Returns the local pressure of the atmosphere parameter in Newton per meter^2. + * \param altitude Altitude. + * \param longitude Longitude. + * \param latitude Latitude. + * \param time Time. + * \return Atmospheric pressure. + */ + virtual double getPressure( const double altitude, const double longitude, + const double latitude, const double time ) + { + + } + + //! Get local temperature. + /*! + * Returns the local temperature of the atmosphere parameter in Kelvin. + * \param altitude Altitude. + * \param longitude Longitude. + * \param latitude Latitude. + * \param time Time. + * \return Atmospheric temperature. + */ + virtual double getTemperature( const double altitude, const double longitude, + const double latitude, const double time ) + { + + } + + //! 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 ) + { + + } + +protected: + +private: + +}; + + +} // namespace aerodynamics + +} // namespace tudat + +#endif // TUDAT_MARS_DTM_ATMOSPHERE_MODEL_H diff --git a/src/astro/aerodynamics/CMakeLists.txt b/src/astro/aerodynamics/CMakeLists.txt index 00fd5a0401..543a806378 100644 --- a/src/astro/aerodynamics/CMakeLists.txt +++ b/src/astro/aerodynamics/CMakeLists.txt @@ -23,7 +23,8 @@ set(aerodynamics_SOURCES "flightConditions.cpp" "trimOrientation.cpp" "equilibriumWallTemperature.cpp" - ) + "marsDtmAtmosphereModel.cpp" +) # Set the header files. set(aerodynamics_HEADERS @@ -46,7 +47,8 @@ set(aerodynamics_HEADERS "aerodynamicGuidance.h" "equilibriumWallTemperature.h" "windModel.h" - ) + "marsDtmAtmosphereModel.h" +) if (TUDAT_BUILD_WITH_NRLMSISE00) set(aerodynamics_SOURCES diff --git a/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp b/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp new file mode 100644 index 0000000000..6a22f18d1a --- /dev/null +++ b/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp @@ -0,0 +1,11 @@ + +#include "tudat/astro/aerodynamics/marsDtmAtmosphereModel.h" + +namespace tudat +{ +namespace aerodynamics +{ + +} // namespace aerodynamics + +} // namespace tudat diff --git a/tests/src/astro/aerodynamics/CMakeLists.txt b/tests/src/astro/aerodynamics/CMakeLists.txt index 7980e45977..5eef35b1d9 100644 --- a/tests/src/astro/aerodynamics/CMakeLists.txt +++ b/tests/src/astro/aerodynamics/CMakeLists.txt @@ -72,6 +72,11 @@ TUDAT_ADD_TEST_CASE(WindModel ${Tudat_PROPAGATION_LIBRARIES} ) +TUDAT_ADD_TEST_CASE(MarsDtmAtmosphere + PRIVATE_LINKS + ${Tudat_PROPAGATION_LIBRARIES} +) + if (TUDAT_BUILD_WITH_NRLMSISE00) TUDAT_ADD_TEST_CASE(NRLMSISE00Atmosphere PRIVATE_LINKS diff --git a/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp b/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp new file mode 100644 index 0000000000..35f5311c0c --- /dev/null +++ b/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp @@ -0,0 +1,43 @@ +/* Copyright (c) 2010-2019, 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 + * a copy of the license with this file. If not, please or visit: + * http://tudat.tudelft.nl/LICENSE. + * + * References + * + */ + +//#define BOOST_TEST_DYN_LINK +//#define BOOST_TEST_MAIN + +#include + +#include + +#include "tudat/astro/aerodynamics/marsDtmAtmosphereModel.h" + +int main( ) +{ + +} +// +//namespace tudat +//{ +//namespace unit_tests +//{ +// +//BOOST_AUTO_TEST_SUITE( test_mars_dtm_atmosphere ) +// +//BOOST_AUTO_TEST_CASE( testMarsDtmAtmosphere ) +//{ +// +//} +// +//BOOST_AUTO_TEST_SUITE_END( ) +// +//} // namespace unit_tests +//} // namespace tudat From 7187ee49777e96bcf92ccb592f1f600aa86f9638 Mon Sep 17 00:00:00 2001 From: RivaAlkahal Date: Tue, 5 Dec 2023 21:43:31 +0100 Subject: [PATCH 02/25] First version of Mars DTM working --- .../aerodynamics/marsDtmAtmosphereModel.h | 144 ++++++- .../aerodynamics/marsDtmAtmosphereModel.cpp | 394 ++++++++++++++++++ src/math/basic/legendrePolynomials.cpp | 32 ++ .../unitTestMarsDtmAtmosphere.cpp | 41 ++ 4 files changed, 608 insertions(+), 3 deletions(-) diff --git a/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h b/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h index dd64cfadda..70e3ea47d8 100644 --- a/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h +++ b/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h @@ -12,19 +12,98 @@ #ifndef TUDAT_MARS_DTM_ATMOSPHERE_MODEL_H #define TUDAT_MARS_DTM_ATMOSPHERE_MODEL_H +#include +#include #include "tudat/astro/aerodynamics/atmosphereModel.h" +#include "tudat/astro/basic_astro/timeConversions.h" +#include "tudat/astro/basic_astro/dateTime.h" namespace tudat { namespace aerodynamics { + class marsDate { + public: + int year; + int month; + int day; + double hours; + double minutes; + double seconds; + + //constructor + //structure of the marsDate to be defined later + marsDate(int y, int m, int d, double hh, double mm, double ss) : year(y), month(m), day(d), hours(hh), minutes(mm), seconds(ss) {} + + //distructor + virtual ~marsDate( ) { } + + // takes the absolute difference between the two dates after converting them to Julian date. + double dateDifference(const marsDate& date1, const marsDate& date2) { + return std::abs( basic_astrodynamics::convertCalendarDateToJulianDaysSinceEpoch( + date1.year, date1.month, date1.day, date1.hours, date1.minutes, date1.seconds, basic_astrodynamics::JULIAN_DAY_ON_J2000) - + basic_astrodynamics::convertCalendarDateToJulianDaysSinceEpoch( + date2.year, date2.month, date2.day, date2.hours, date2.minutes, date2.seconds, basic_astrodynamics::JULIAN_DAY_ON_J2000)); + } + + // function to find the nearest date in the list of dates provided in the next variable + marsDate findNearestDate(const marsDate &targetDate); + + double marsDayofYear(const marsDate &targetDate); + }; + + //from: https://www.planetary.org/articles/mars-calendar + static std::vector marsYears = { + marsDate(1955, 4, 11, 0, 0, 0), + marsDate(1957, 2, 26, 0, 0, 0), + marsDate(1959, 1, 14, 0, 0, 0), + marsDate(1960, 12, 1, 0, 0, 0), + marsDate(1962, 10, 19, 0, 0, 0), + marsDate(1964, 9, 5, 0, 0, 0), + marsDate(1966, 7, 24, 0, 0, 0), + marsDate(1968, 6, 10, 0, 0, 0), + marsDate(1970, 4, 28, 0, 0, 0), + marsDate(1972, 3, 15, 0, 0, 0), + marsDate(1974, 1, 31, 0, 0, 0), + marsDate(1975, 12, 19, 0, 0, 0), + marsDate(1977, 11, 5 ,0, 0, 0), + marsDate(1979, 9, 23, 0, 0, 0), + marsDate(1981, 8, 10, 0, 0, 0), + marsDate(1983, 6, 28, 0, 0, 0), + marsDate(1985, 5, 15, 0, 0, 0), + marsDate(1987, 4, 1, 0, 0, 0), + marsDate(1989, 2, 16, 0, 0, 0), + marsDate(1991, 1, 4, 0, 0, 0), + marsDate(1992, 11, 21 , 0, 0, 0), + marsDate(1994, 10, 9 , 0, 0, 0), + marsDate(1996, 8, 26 , 0, 0, 0), + marsDate(1998, 7, 14 , 0, 0, 0), + marsDate(2000, 5, 31 , 0, 0, 0), + marsDate(2002, 4, 18 , 0, 0, 0), + marsDate(2004, 3, 5 , 0, 0, 0), + marsDate(2006, 1, 21 , 0, 0, 0), + marsDate(2007, 12, 9 , 0, 0, 0), + marsDate(2009, 10, 26 , 0, 0, 0), + marsDate(2011, 9, 13 , 0, 0, 0), + marsDate(2013, 7, 31 , 0, 0, 0), + marsDate(2015, 6, 18 , 0, 0, 0), + marsDate(2017, 5, 5 , 0, 0, 0), + marsDate(2019, 3, 23 , 0, 0, 0), + marsDate(2021, 2, 7 , 0, 0, 0), + marsDate(2022, 12, 26 , 0, 0, 0), + marsDate(2024, 11, 12 , 0, 0, 0), + marsDate(2026, 9, 30 , 0, 0, 0), + marsDate(2028, 8, 17 , 0, 0, 0) + }; + class MarsDtmAtmosphereModel: public AtmosphereModel { public: - MarsDtmAtmosphereModel( ){ } + + MarsDtmAtmosphereModel(const double polarRadius, const std::string &filename); //! Default destructor. /*! @@ -32,6 +111,7 @@ class MarsDtmAtmosphereModel: public AtmosphereModel */ virtual ~MarsDtmAtmosphereModel( ) { } + //! Get local density. /*! * Returns the local density parameter of the atmosphere in kg per meter^3. @@ -44,7 +124,8 @@ class MarsDtmAtmosphereModel: public AtmosphereModel virtual double getDensity( const double altitude, const double longitude, const double latitude, const double time ) { - + computeProperties( altitude, longitude, latitude, time ); + return currentDensity_; } //! Get local pressure. @@ -74,7 +155,8 @@ class MarsDtmAtmosphereModel: public AtmosphereModel virtual double getTemperature( const double altitude, const double longitude, const double latitude, const double time ) { - + computeProperties( altitude, longitude, latitude, time ); + return currentTemperature_z; } //! Get local speed of sound. @@ -92,10 +174,66 @@ class MarsDtmAtmosphereModel: public AtmosphereModel } + + + const double Omega = 2.0*mathematical_constants::PI/686.98;//(686.98*24.63);//rad/h//2*M_PI/(659355072); //rad/s + const double omega = 2.0*mathematical_constants::PI/88668;//24.63;//2*M_PI/(88668); // rad/s + + double computeLocalSolarTime(const double longitude, const int day, const int month, const int year, const double hours, + const double minutes); + void updateLegendrePolynomials( const double latitude ); + double computeGl(const double latitude, const double longitude, const double minutes, const double hours_, const int day_, + const int month_, const int year_, const int indexg); + double computeGeopotentialAltitude( const double altitude ); + + double computeCurrentTemperature( const double latitude, const double longitude, const double minutes_, const double hours_, const int day_ ,const int month_,const int year_, const int indexg); + + double computeGamma(const double latitude, const double longitude, const double minutes_, const double hours_, const int day_ ,const int month_,const int year_, const int indexm); + + double heightDistributionFunction(const double altitude, const double latitude, + const double longitude, const double minutes_, const double hours_, const int day_ ,const int month_,const int year_,const int indexm); + + double getTotalDensity( const double altitude, const double latitude, + const double longitude, const double minutes_, const double hours_, const int day_ , const int month_, const int year_); + + double getCurrentLegendrePolynomial( const int degree ) + { + return currentLegendrePolynomials_.at( degree ); + } + + double computeGl_Subr(const double latitude, const double longitude, const double minutes_, const double hours_, + const int day_, const int month_, const int year_, const int indexg); + + void defineDustPars(const int rows); + + double computeDustStorm(const double Ls); + protected: private: + void computeProperties( const double altitude, const double longitude, + const double latitude, const double time ); + + double polarRadius_; + std::basic_string, std::allocator> filename_; + + std::vector< double > alpha_; + std::vector< double > mmass_; + + std::vector> coefficients_; + std::vector taus; + + double Ls; + std::vector< double > currentLegendrePolynomials_; + double currentGeopotentialAltitude_; + double currentTemperature_138; + double currentTemperature_inf; + double currentdTemperature_; + double currentDensity_; + double sigma; + double currentTemperature_z; + }; diff --git a/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp b/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp index 6a22f18d1a..21c61bff8a 100644 --- a/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp +++ b/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp @@ -1,11 +1,405 @@ +#include +#include +#include #include "tudat/astro/aerodynamics/marsDtmAtmosphereModel.h" +#include "tudat/astro/basic_astro/timeConversions.h" +#include "tudat/math/basic/legendrePolynomials.h" +#include "tudat/astro/basic_astro/unitConversions.h" +#include "tudat/astro/basic_astro/celestialBodyConstants.h" + namespace tudat { namespace aerodynamics { + void MarsDtmAtmosphereModel::computeProperties( + const double altitude, const double longitude, + const double latitude, const double time ) + { + // Compute the hash key?? + basic_astrodynamics::DateTime currentDateTime_ = basic_astrodynamics::getCalendarDateFromTime( time ); + currentDensity_ = getTotalDensity( + altitude, longitude, latitude, + currentDateTime_.getMinute(), currentDateTime_.getHour(), currentDateTime_.getDay(), currentDateTime_.getMonth(), currentDateTime_.getYear() ); + currentGeopotentialAltitude_ = computeGeopotentialAltitude( altitude ); + currentTemperature_138 = computeCurrentTemperature( + latitude, longitude, currentDateTime_.getMinute(), currentDateTime_.getHour(), currentDateTime_.getDay(), currentDateTime_.getMonth(), currentDateTime_.getYear(), 3 ); + currentTemperature_inf = computeCurrentTemperature( + latitude, longitude, currentDateTime_.getMinute(), currentDateTime_.getHour(), currentDateTime_.getDay(), currentDateTime_.getMonth(), currentDateTime_.getYear(), 1 ); + currentdTemperature_ = computeCurrentTemperature( + latitude, longitude, currentDateTime_.getMinute(), currentDateTime_.getHour(), currentDateTime_.getDay(), currentDateTime_.getMonth(), currentDateTime_.getYear(), 5 ); + sigma = currentdTemperature_/(currentTemperature_inf-currentTemperature_138); + currentTemperature_z = currentTemperature_inf - (currentTemperature_inf - currentTemperature_138)* exp(-sigma*currentGeopotentialAltitude_); + } + // function to load the coefficients from the file + std::vector> loadCoefficients(const std::string& filename){ + std::ifstream file(filename); + // define the number of rows and columns in the file + const int rows = 70; + //const int rows = 74; + const int cols = 25; + std::vector> matrix(rows, std::vector(cols)); + + for (int i = 0; i < rows; ++i) { + for (int j = 0; j < cols; ++j) { + if (!(file >> matrix[i][j])) { + // if the file is not read correctly, print an error message + std::cerr << "Error reading data from the file." << std::endl; + // return an error code + } + } + } + //std::cout << "Loaded Matrix values:" << std::endl; + //std::cout << matrix[0][3] << "\t"; + file.close(); + return matrix; + } + + // function to provide the nearest date in a list of dates to a target date (i + marsDate tudat::aerodynamics::marsDate::findNearestDate(const marsDate& targetDate){ + if (targetDate.year < marsYears[0].year){ + std::cerr << "The target date is before the first date of the Martian years: 1955, 4, 11. Day of year can't be calculated. PROGRAM TERMINATED!!" << std::endl; + exit(EXIT_FAILURE); + } + marsDate nearestDate = marsYears[0]; + double minDifference = marsDate::dateDifference(targetDate, nearestDate); + for (const marsDate& date : marsYears){ + if (date.year < targetDate.year || (date.year == targetDate.year && date.month < targetDate.month) || (date.year == targetDate.year && date.month == targetDate.month && date.day < targetDate.day)){ + double difference = marsDate::dateDifference(targetDate, date); + if (difference < minDifference){ + minDifference = difference; + nearestDate = date; + } + } + + } + return nearestDate; + + } + + // function to compute the day of martian year + double tudat::aerodynamics::marsDate::marsDayofYear(const marsDate& targetDate){ + marsDate nearestDate = findNearestDate(targetDate); + double difference = marsDate::dateDifference(targetDate, nearestDate); + // convert the difference in days to seconds, and devide them by the number of seconds in a martian year to get the day of year + double dayOfYear = difference*24*3600 / 88668.0; // 88668 seconds is the length of a martian day (24.63 h * 3600 s/h) + return dayOfYear; + } + + MarsDtmAtmosphereModel::MarsDtmAtmosphereModel(const double polarRadius, const std::string &filename) : + polarRadius_( polarRadius ), // polar radius of Mars + filename_ ( filename ), // file name of the coefficients + alpha_( {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.38, -0.40, 0.0}), // thermal diffusion coefficients + mmass_( {44.01, 16.00, 28.0, 40.0, 28.0, 32.0, 4.0, 1.0, 2.0} )// molar mass of the species + { + coefficients_ = loadCoefficients( filename ); + currentLegendrePolynomials_.resize( 6 ); + currentLegendrePolynomials_[ 0 ] = 0.0; + std::cout << "Loaded Matrix values:" << std::endl; + std::cout << coefficients_[0][3] << "\t"; + } + + // function to compute the local true solar time + // equations from: https://www.giss.nasa.gov/tools/mars24/help/algorithm.html + double MarsDtmAtmosphereModel::computeLocalSolarTime(const double longitude, const int day ,const int month,const int year, const double hours, const double minutes) + { + //Day since J2000 epoch in Terrestrial Time (or TAI) + double timeElaspsedJD =basic_astrodynamics::convertCalendarDateToJulianDaysSinceEpoch( + year, month, day, hours, minutes, 0.0, basic_astrodynamics::JULIAN_DAY_ON_J2000 ); + //std::cout << "timeElaspsedJD: " << timeElaspsedJD << std::endl; + //Mars orbital parameters: + //Mars Mean Anomaly + double meanAnomaly = 19.3870 + 0.52402073 * timeElaspsedJD; //degrees + //std::cout << "MeanAnomaly: " << MeanAnomaly << std::endl; + //Angle of Fiction Mean Sun + double omegaFMS = 270.3863 + 0.52403840 * timeElaspsedJD; + //std::cout << "OmegaFMS: " << OmegaFMS << std::endl; + //Perturbers + double PBS = 0.0071 * cos( unit_conversions::convertDegreesToRadians(0.985626*timeElaspsedJD/2.2353 + 49.409)) + + 0.0057 * cos( unit_conversions::convertDegreesToRadians(0.985626*timeElaspsedJD/2.7543 + 168.173)) + + 0.0039 * cos( unit_conversions::convertDegreesToRadians(0.985626*timeElaspsedJD/1.1177 + 191.837)) + + 0.0037 * cos( unit_conversions::convertDegreesToRadians(0.985626*timeElaspsedJD/15.7866 + 21.736)) + + 0.0021 * cos( unit_conversions::convertDegreesToRadians(0.985626*timeElaspsedJD/2.1354 + 15.704)) + + 0.0020 * cos( unit_conversions::convertDegreesToRadians(0.985626*timeElaspsedJD/2.4694 + 95.528)) + + 0.0018 * cos( unit_conversions::convertDegreesToRadians(0.985626*timeElaspsedJD/32.8493 + 49.095)); + //std::cout << "PBS: " << PBS << std::endl; + //Equation of Center + double EoC = (10.691 + 3.0E-7 * timeElaspsedJD) * sin( unit_conversions::convertDegreesToRadians(meanAnomaly)) + + 0.623 * sin( unit_conversions::convertDegreesToRadians(2*meanAnomaly)) + + 0.050 * sin( unit_conversions::convertDegreesToRadians(3*meanAnomaly)) + + 0.005 * sin( unit_conversions::convertDegreesToRadians(4*meanAnomaly)) + + 0.0005 * sin( unit_conversions::convertDegreesToRadians(5*meanAnomaly)) + PBS;//Equation of Center + //std::cout << "EoC: " << EoC << std::endl; + //Areocentric Solar Longitude + double Ls1 = omegaFMS + EoC; + Ls = basic_mathematics::computeModulo( Ls1,360.0); + std::cout << "Ls: " << Ls << std::endl; + //Equation of Time + double EoT = 2.861 * sin( unit_conversions::convertDegreesToRadians(2*Ls)) + - 0.071 * sin( unit_conversions::convertDegreesToRadians(4*Ls)) + + 0.002 * sin( unit_conversions::convertDegreesToRadians(6*Ls)) - EoC; // Equation of Time + //std::cout << "EoT: " << EoT << std::endl; + //Mean Solar Time at Mars's prime merdian, i.e., Airy Mean Time + double MST = 24*((basic_astrodynamics::convertCalendarDateToJulianDay(year,month,day,hours,minutes,0.0) - 2451549.5)/1.0274912517 + 44796.0 - 0.0009626); //double MST + //std::cout << "MST1: " << MST << std::endl; + MST = basic_mathematics::computeModulo( MST,24.0);// Modulo MST + std::cout << "MST: " << MST << std::endl; + //Local Mean Solar Time + double LMST = basic_mathematics::computeModulo( MST-longitude*(24.0/360.0),24.0); + //Local True Solar Time + double LTST = basic_mathematics::computeModulo(LMST + EoT*(24.0/360.0),24.0); + std::cout << "LTST: " << LTST << std::endl; + //Subsolar Longitude + double subSolarLongitude = MST*(360.0/24.0) + EoT + 180.0; + // return LTST in seconds. + return LTST*3600.0; + } + // Function to define the dust storm parameters + void MarsDtmAtmosphereModel::defineDustPars(const int rows ) + { + double tau = 0.0; // dust opacity + double lsdeb = 0.0; // dust storm onset + double lsfin = 0.0; // dust storm end + double taubg = 0.0; // background dust opacity + + if (rows >= 74){ + tau = coefficients_[1][rows-3]; + lsdeb = coefficients_[1][rows-2]; + lsfin = coefficients_[1][rows-1]; + taubg = coefficients_[1][rows]; + + if (taubg<0.15){ + taubg = 0.0; + } + + } + + taus = {tau, lsdeb, lsfin, taubg}; + + + } + + // function to compute the dust storm + double MarsDtmAtmosphereModel::computeDustStorm(const double Ls){ + MarsDtmAtmosphereModel::defineDustPars(70); + + double dells = 15.0; //delta Ls of Stewart + double xlsmaxd = taus[1] + 6.0; //tau max reached after 6 degree, and xlsdev is storm onset + double dectau = 12.0; // decay of tau = max to normal in 12 degree of Ls + double xlsmaxf = taus[2] - dectau; // start os storm decay (tau) + double deuxse = 2.0/2.78128; + double correc = pow(deuxse,dectau); + double amr = 0.0; + double da41 = 0.0; + //Dust storm + if (Ls >= taus[1] && Ls >= 180){ + if (Ls < xlsmaxd){ + double q = exp(-(Ls - taus[1])/dells); + double da41 = q*(1.0-pow(q,4))*taus[0]; + } + if (Ls >= xlsmaxd && Ls < xlsmaxf){ + double q = exp(-(xlsmaxd - Ls)/dells); + double da41 = q*(1.0-pow(q,4))*taus[0]; + } + if (Ls >= xlsmaxf && Ls < taus[2]){ + double qmax = exp(-(xlsmaxd - Ls)/dells); + double taumax = qmax*(1.0-pow(qmax,4))*taus[0]; + double xmult = taumax/deuxse; + double ipow = (Ls - xlsmaxf) + 1.0; + double denom = pow(deuxse,ipow); + double da41 = xmult*(denom-correc); + } + + } + return da41; + + } + + void MarsDtmAtmosphereModel::updateLegendrePolynomials( const double latitude ) + { + for( int i = 1; i <= 6; i++ ) + { + currentLegendrePolynomials_[ i ] = basic_mathematics::computeLegendrePolynomialExplicit( i, 0, sin( latitude ) ); + } + } + + // Function to compute the spherical harmonic expansions defined as G(l) as in the paper (Bruinsma and Lemoine, 2002) + double MarsDtmAtmosphereModel::computeGl(const double latitude, const double longitude, const double minutes_, const double hours_, const int day_ ,const int month_,const int year_, const int indexg) + { + updateLegendrePolynomials( latitude ); + using basic_mathematics::computeLegendrePolynomialExplicit; + marsDate date2 = marsDate(year_, month_, day_, hours_ , minutes_, 0.0); + double doy = date2.marsDayofYear(date2); //day of year + std::cout << "doy: " << doy << std::endl; + double t = computeLocalSolarTime(longitude, day_, month_, year_, hours_, minutes_); //seconds + int F = 0.0; // flux parameter is set to 0 for now. + std::cout<<"coefficients_[14][indexg]" << coefficients_[14][indexg] << std::endl; + // Non-periodic term + double NP = (coefficients_[1][indexg] * F) + coefficients_[30][indexg] * currentLegendrePolynomials_[1] + + coefficients_[31][indexg] * currentLegendrePolynomials_[ 2 ] + + coefficients_[32][indexg] * currentLegendrePolynomials_[ 3 ] + + coefficients_[33][indexg] * currentLegendrePolynomials_[ 4 ] + + coefficients_[34][indexg] * currentLegendrePolynomials_[ 5 ]; + // Annual-Periodic terms + double PA = (coefficients_[2][indexg]+ coefficients_[4][indexg] * currentLegendrePolynomials_[2] + coefficients_[6][indexg]*F)* cos(Omega * doy) + + (coefficients_[3][indexg]+ coefficients_[5][indexg] * currentLegendrePolynomials_[2] + coefficients_[7][indexg]*F)* sin(Omega * doy) + + (coefficients_[8][indexg] * currentLegendrePolynomials_[1] + coefficients_[10][indexg] * currentLegendrePolynomials_[3] + coefficients_[12][indexg] *currentLegendrePolynomials_[5] ) * cos(Omega * doy) + + (coefficients_[9][indexg] * currentLegendrePolynomials_[1] + coefficients_[11][indexg] * currentLegendrePolynomials_[3] + coefficients_[13][indexg] * currentLegendrePolynomials_[5]) * sin(Omega * doy); + //Semi-Annual-Periodic terms + double PSA = (coefficients_[14][indexg] + coefficients_[16][indexg] * currentLegendrePolynomials_[2] + coefficients_[18][indexg]*F) * cos(2.0*Omega * doy) + + (coefficients_[15][indexg]+ coefficients_[17][indexg] * currentLegendrePolynomials_[2] + coefficients_[19][indexg]*F)* sin(2.0*Omega * doy) + + (coefficients_[20][indexg] * currentLegendrePolynomials_[1] + coefficients_[22][indexg] * currentLegendrePolynomials_[3]) * cos(2.0*Omega * doy) + + (coefficients_[21][indexg] * currentLegendrePolynomials_[1] + coefficients_[23][indexg] * currentLegendrePolynomials_[3]) * sin(2.0*Omega * doy); + //Diurnal terms + double PD = (coefficients_[24][indexg]*computeLegendrePolynomialExplicit(1,1,sin(latitude)) + coefficients_[26][indexg]*computeLegendrePolynomialExplicit(2,1,sin(latitude)) + coefficients_[28][indexg]*computeLegendrePolynomialExplicit(3,1,sin(latitude)))*cos(omega*t) + + (coefficients_[25][indexg]*computeLegendrePolynomialExplicit(1,1,sin(latitude)) + coefficients_[27][indexg]*computeLegendrePolynomialExplicit(2,1,sin(latitude)) + coefficients_[29][indexg]*computeLegendrePolynomialExplicit(3,1,sin(latitude)))*sin(omega*t); + double Gl = NP + PA + PSA + PD; + return Gl; + } + + //Function to compute the spherical harmonic expansions defined as G(l) as in the subroutine provided by Bruinsma + double MarsDtmAtmosphereModel::computeGl_Subr(const double latitude, const double longitude, const double minutes_, const double hours_, const int day_ ,const int month_,const int year_, const int indexg) + { + //Get the zonal Legendre polynomials for the current latitude + updateLegendrePolynomials( latitude ); + using basic_mathematics::computeLegendrePolynomialExplicit; + //Get the dust storm parameters + MarsDtmAtmosphereModel::defineDustPars(70); + //Get the day of year + marsDate date2 = marsDate(year_, month_, day_, hours_ , minutes_, 0.0); + double doy = date2.marsDayofYear(date2); + //Get the local solar time + double t = computeLocalSolarTime(longitude, day_, month_, year_, hours_, minutes_); //hours + //Get the local solar time in radians + pi. + double hl0 = omega*t+ mathematical_constants::PI; + //std::cout << "hl: " << omega*t << std::endl; + double cos2h = cos(hl0)* cos(hl0) - sin(hl0)*sin(hl0); + double sin2h = 2.0*sin(hl0)*cos(hl0); + //flux terms: + double ff0 = 0.0; + double F107 = 65.0; + double Fbar = 65.0; + double F = Fbar - F107; + double F2 = F*F; + double f0 = coefficients_[4][indexg]*F + coefficients_[39][indexg] * F2; + double f1f = 1.0 + f0*ff0; // coupling terms + //flux + latitude terms + f0 = f0+ coefficients_[1][indexg]*currentLegendrePolynomials_[2] + coefficients_[2][indexg]*currentLegendrePolynomials_[3] + + coefficients_[3][indexg]*currentLegendrePolynomials_[4] + coefficients_[58][indexg]*currentLegendrePolynomials_[6] + + coefficients_[59][indexg]*currentLegendrePolynomials_[1] + coefficients_[60][indexg]*currentLegendrePolynomials_[5]; + // symmetrical and seasonal annual terms + double PA = coefficients_[5][indexg]* cos(Omega*doy) + coefficients_[6][indexg]* sin(Omega*doy) + + coefficients_[7][indexg]* cos(Omega*doy) * F + coefficients_[8][indexg]* sin(Omega*doy) * F + + coefficients_[9][indexg]* cos(Omega*doy) * currentLegendrePolynomials_[2] + coefficients_[10][indexg]* sin(Omega*doy) * currentLegendrePolynomials_[2] + + coefficients_[11][indexg] * cos(Omega*doy) * currentLegendrePolynomials_[1] + coefficients_[12][indexg] * sin(Omega*doy) * currentLegendrePolynomials_[1] + + coefficients_[17][indexg] * cos(Omega*doy) * currentLegendrePolynomials_[3] + coefficients_[18][indexg] * sin(Omega*doy) * currentLegendrePolynomials_[3] + + coefficients_[61][indexg] * cos(Omega*doy) * currentLegendrePolynomials_[1] *F + coefficients_[62][indexg] * sin(Omega*doy) * currentLegendrePolynomials_[1] *F + + coefficients_[63][indexg] * cos(Omega*doy) * currentLegendrePolynomials_[5] + coefficients_[64][indexg] * sin(Omega*doy) * currentLegendrePolynomials_[5]; + //symmetrical and seasonal semi-annual terms + double PSA = coefficients_[13][indexg] * cos(2.0*Omega*doy) + coefficients_[14][indexg] * sin(2.0*Omega*doy) + + coefficients_[15][indexg] * cos(2.0*Omega*doy) * F + coefficients_[16][indexg] * sin(2.0*Omega*doy) * F + + coefficients_[19][indexg] * cos(2.0*Omega*doy) * currentLegendrePolynomials_[1] + coefficients_[20][indexg] * sin(2.0*Omega*doy) * currentLegendrePolynomials_[1] + + coefficients_[65][indexg] * cos(2.0*Omega*doy) * currentLegendrePolynomials_[2] + coefficients_[66][indexg] * sin(2.0*Omega*doy) * currentLegendrePolynomials_[2] + + coefficients_[67][indexg] * cos(2.0*Omega*doy) * currentLegendrePolynomials_[3] + coefficients_[68][indexg] * sin(2.0*Omega*doy) * currentLegendrePolynomials_[3]; + // diurnal terms + double PD = coefficients_[21][indexg] * computeLegendrePolynomialExplicit(1,1,sin(latitude)) * cos(hl0) + coefficients_[22][indexg] * computeLegendrePolynomialExplicit(1,1,sin(latitude)) * sin(hl0) + + coefficients_[23][indexg] * computeLegendrePolynomialExplicit(2,1,sin(latitude)) * cos(hl0) + coefficients_[24][indexg] * computeLegendrePolynomialExplicit(2,1,sin(latitude)) * sin(hl0) + + coefficients_[25][indexg] * computeLegendrePolynomialExplicit(3,1,sin(latitude)) * cos(hl0) + coefficients_[26][indexg] * computeLegendrePolynomialExplicit(3,1,sin(latitude)) * sin(hl0) + + coefficients_[27][indexg] * computeLegendrePolynomialExplicit(1,1,sin(latitude)) * cos(hl0) * F + coefficients_[28][indexg] * computeLegendrePolynomialExplicit(1,1,sin(latitude)) * sin(hl0) * F + + coefficients_[30][indexg] * computeLegendrePolynomialExplicit(2,2,sin(latitude)) * (cos2h) + coefficients_[31][indexg] * computeLegendrePolynomialExplicit(2,2,sin(latitude)) * sin2h + + coefficients_[32][indexg] * computeLegendrePolynomialExplicit(3,2,sin(latitude)) * (cos2h) + coefficients_[33][indexg] * computeLegendrePolynomialExplicit(3,2,sin(latitude)) * sin2h + + coefficients_[34][indexg] * computeLegendrePolynomialExplicit(4,2,sin(latitude)) * (cos2h) + coefficients_[35][indexg] * computeLegendrePolynomialExplicit(4,2,sin(latitude)) * sin2h + + coefficients_[36][indexg] * computeLegendrePolynomialExplicit(2,2,sin(latitude)) * (cos2h) * F + coefficients_[37][indexg] * computeLegendrePolynomialExplicit(2,2,sin(latitude)) * sin2h * F; + //dust terms + double da41 = computeDustStorm( Ls ); + double fpds = coefficients_[40][indexg]*da41 + coefficients_[69][indexg]*taus[3]; + double Gl = fpds + f0 + PA + PSA + PD; + std::cout << "f0: " << f0 << std::endl; + std::cout << "PA: " << PA << std::endl; + std::cout << "PSA: " << PSA << std::endl; + std::cout << "PD: " << PD << std::endl; + return Gl; + } + + // Function to compute the geopotential altitude + double MarsDtmAtmosphereModel::computeGeopotentialAltitude( const double altitude ) + { + return ((altitude - 138.0E3)*(polarRadius_ + 138.0E3 )/(polarRadius_+altitude))/1000.0; //km //removed the 80 from Bruinsma's formula, because in Bates, 1959, it is not there, results seem to be more realistic withoout the 80 + } + + // Function to compute the current temperature (exospheric, 138 km, and partial temperatures) + double MarsDtmAtmosphereModel::computeCurrentTemperature( const double latitude, const double longitude, const double minutes_, const double hours_, const int day_ ,const int month_,const int year_, const int indexg) + { + double T0 = coefficients_[0][indexg]; + double Ti; + double Gl = computeGl_Subr(latitude, longitude, minutes_, hours_, day_ , month_, year_, indexg); + Ti = T0*(1.0 + Gl); + return Ti; + } + // Function to compute gamma parameter + double MarsDtmAtmosphereModel::computeGamma( const double latitude, const double longitude, const double minutes_, const double hours_, const int day_ ,const int month_,const int year_, const int indexm) + { + double universalGasConstant = 8.314;//physical_constants::MOLAR_GAS_CONSTANT; //J/mol/K //kg m^2 / s^2 / K / mol + //double g138 =celestial_body_constants::MARS_GRAVITATIONAL_PARAMETER/((138.0E3+3376.78E3)*(138.0E3+3376.78E3)); //m/s2 + //Get gravity at 138 km + double g138 = 3.727/pow((1+(138.0/3376.78)),2); //m/s2 + currentTemperature_138 = computeCurrentTemperature( latitude, longitude, minutes_, hours_, day_ , month_, year_, 3); //K + currentTemperature_inf = computeCurrentTemperature( latitude, longitude, minutes_, hours_, day_ , month_, year_, 1); //K + currentdTemperature_ = computeCurrentTemperature( latitude, longitude, minutes_, hours_, day_ , month_, year_, 5); //K + sigma = currentdTemperature_/(currentTemperature_inf-currentTemperature_138); + //std::cout << "sigma: " << sigma << std::endl; + double gamma = mmass_[indexm] * g138 / ( universalGasConstant * sigma * currentTemperature_inf); //km^-1 + /* + std::cout << "gamma: " << gamma << std::endl; + + std::cout << "currentTemperature_138: " << currentTemperature_138 << std::endl; + std::cout << "currentTemperature_inf: " << currentTemperature_inf << std::endl; + std::cout << "currentdTemperature_: " << currentdTemperature_ << std::endl; + */ + return gamma; + } + + // Function to compute the height distribution function + double MarsDtmAtmosphereModel::heightDistributionFunction(const double altitude, const double latitude, + const double longitude, const double minutes_, const double hours_, const int day_ ,const int month_,const int year_,const int indexm) + { + currentGeopotentialAltitude_ = computeGeopotentialAltitude( altitude ); //km + double gamma= computeGamma( latitude, longitude, minutes_, hours_, day_ , month_, year_, indexm); //km^-1 + currentTemperature_z = currentTemperature_inf - (currentTemperature_inf - currentTemperature_138)* exp(-sigma*currentGeopotentialAltitude_); + //std::cout << "Tz: " << Tz << std::endl; + double fi = pow((currentTemperature_138/currentTemperature_z),(1+alpha_[indexm]+gamma))* exp(-sigma*gamma*currentGeopotentialAltitude_); + //std::cout << "fi: " << fi << std::endl; + return fi; + } + + // Function to compute the total density + double MarsDtmAtmosphereModel::getTotalDensity( const double altitude, const double latitude, + const double longitude, const double minutes_, const double hours_, const int day_ ,const int month_,const int year_) + { + double avogadroConstant_ = 6.022E23;//physical_constants::AVOGADRO_CONSTANT; + double rho0; + double rho = 0; + int indexm = 0; + for (int col = 7; col < 25; col += 2) { + rho0 = coefficients_[0][col]*mmass_[indexm]/avogadroConstant_; //g/cm3 + //std::cout << "coefficients_[0][col]: " << coefficients_[0][col] << std::endl; + double fi = heightDistributionFunction(altitude, latitude, longitude, minutes_, hours_, day_ , month_, year_, indexm); + if (col == 7) + std::cout << "fi CO2: " << fi << std::endl; + + double gl = computeGl_Subr(latitude, longitude, minutes_, hours_, day_ , month_, year_, col); + std::cout << "gl: " << gl << std::endl; + std::cout << "col" << col << std::endl; + rho += rho0*fi*exp(gl); + indexm +=1; + std::cout << "rho: " << indexm << " " << rho << std::endl; + } + return rho; + } + + } // namespace aerodynamics } // namespace tudat diff --git a/src/math/basic/legendrePolynomials.cpp b/src/math/basic/legendrePolynomials.cpp index df23f9a63b..8b44716108 100644 --- a/src/math/basic/legendrePolynomials.cpp +++ b/src/math/basic/legendrePolynomials.cpp @@ -695,6 +695,38 @@ double computeLegendrePolynomialExplicit( const int degree, } } break; + case 5: + switch( order ) + { + case 0: + return ( 63.0 * polynomialParameter * polynomialParameter + * polynomialParameter * polynomialParameter * polynomialParameter + - 70.0 * polynomialParameter * polynomialParameter * polynomialParameter + + 15.0 * polynomialParameter ) / 8.0; + default: + { + std::string errorMessage = "Error, explicit legendre polynomial not possible for " + + std::to_string( degree ) + ", " + + std::to_string( order ); + throw std::runtime_error( errorMessage ); + } + } + case 6: + switch( order ) + { + case 0: + return ( 231.0 * polynomialParameter * polynomialParameter + * polynomialParameter * polynomialParameter * polynomialParameter * polynomialParameter + - 315.0 * polynomialParameter * polynomialParameter * polynomialParameter * polynomialParameter + + 105.0 * polynomialParameter * polynomialParameter - 5.0 ) / 16.0; + default: + { + std::string errorMessage = "Error, explicit legendre polynomial not possible for " + + std::to_string( degree ) + ", " + + std::to_string( order ); + throw std::runtime_error( errorMessage ); + } + } default: { std::string errorMessage = "Error, explicit legendre polynomial not possible for " + diff --git a/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp b/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp index 35f5311c0c..54dcd2ae62 100644 --- a/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp +++ b/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp @@ -15,14 +15,55 @@ //#define BOOST_TEST_MAIN #include +#include "fstream" +#include "iostream" #include +#include "tudat/astro/basic_astro/timeConversions.h" #include "tudat/astro/aerodynamics/marsDtmAtmosphereModel.h" +using namespace tudat::aerodynamics; int main( ) { + //std::string filename = "/Users/ralkahal/Documents/PhD/atmodensitypds/dtm_mars"; + std::string filename = "/Users/ralkahal/Documents/PhD/Programs/atmodensitydtm/dtm_mars"; + MarsDtmAtmosphereModel atmosphereModel = MarsDtmAtmosphereModel(3376.78E3, filename); + //atmosphereModel.computelocalsolartime(0.0,16 ,12,2000, 0.0, 1.0697333); + std::cout<(altitude); + //double rho = atmosphereModel.getTotalDensity( alt_km, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000); + + // Write altitude and corresponding density to the file + // outputFile << alt_km << " " << rho << "\n"; + //} + + // Close the file when you are done + //outputFile.close(); + //std::cout << "Density computation and output written to file successfully." << std::endl; + std::cout << atmosphereModel.getTotalDensity( 400.0E3, 0.785398, 0.0, 1.0697333, 0.0, 16 ,12, 2000) << std::endl; + +// return 0; + //std::cout << atmosphereModel.getTotalDensity( 138.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000) < Date: Fri, 8 Dec 2023 13:37:09 +0100 Subject: [PATCH 03/25] Intermediate commit --- src/astro/aerodynamics/marsDtmAtmosphereModel.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp b/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp index 21c61bff8a..754ae1ce35 100644 --- a/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp +++ b/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp @@ -334,7 +334,7 @@ namespace aerodynamics { double T0 = coefficients_[0][indexg]; double Ti; - double Gl = computeGl_Subr(latitude, longitude, minutes_, hours_, day_ , month_, year_, indexg); + double Gl = computeGl(latitude, longitude, minutes_, hours_, day_ , month_, year_, indexg); Ti = T0*(1.0 + Gl); return Ti; } @@ -389,7 +389,7 @@ namespace aerodynamics if (col == 7) std::cout << "fi CO2: " << fi << std::endl; - double gl = computeGl_Subr(latitude, longitude, minutes_, hours_, day_ , month_, year_, col); + double gl = computeGl(latitude, longitude, minutes_, hours_, day_ , month_, year_, col); std::cout << "gl: " << gl << std::endl; std::cout << "col" << col << std::endl; rho += rho0*fi*exp(gl); From 9bc756118b52b66d8c9d16a979fd8f87b757c3c6 Mon Sep 17 00:00:00 2001 From: RivaAlkahal Date: Fri, 8 Dec 2023 13:37:34 +0100 Subject: [PATCH 04/25] Intermediate commit --- tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp b/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp index 54dcd2ae62..42286697cd 100644 --- a/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp +++ b/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp @@ -60,7 +60,7 @@ int main( ) // Close the file when you are done //outputFile.close(); //std::cout << "Density computation and output written to file successfully." << std::endl; - std::cout << atmosphereModel.getTotalDensity( 400.0E3, 0.785398, 0.0, 1.0697333, 0.0, 16 ,12, 2000) << std::endl; + std::cout << atmosphereModel.getTotalDensity( 1000.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000) << std::endl; // return 0; //std::cout << atmosphereModel.getTotalDensity( 138.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000) < Date: Mon, 11 Dec 2023 11:31:23 +0100 Subject: [PATCH 05/25] Added link of Mars DTM to environment settings --- .../environment_setup/createAtmosphereModel.h | 28 +++++++++++++++++++ .../createAtmosphereModel.cpp | 16 +++++++++++ 2 files changed, 44 insertions(+) diff --git a/include/tudat/simulation/environment_setup/createAtmosphereModel.h b/include/tudat/simulation/environment_setup/createAtmosphereModel.h index 5777e9f0fe..bcdaca076d 100644 --- a/include/tudat/simulation/environment_setup/createAtmosphereModel.h +++ b/include/tudat/simulation/environment_setup/createAtmosphereModel.h @@ -172,6 +172,7 @@ enum AtmosphereTypes custom_constant_temperature_atmosphere, tabulated_atmosphere, nrlmsise00, + mars_dtm_atmosphere, scaled_atmosphere }; @@ -528,6 +529,26 @@ class NRLMSISE00AtmosphereSettings: public AtmosphereSettings }; +class MarsDtmAtmosphereSettings: public AtmosphereSettings +{ +public: + + MarsDtmAtmosphereSettings( const std::string& marsDtmFile, const double polarRadius ): + AtmosphereSettings( mars_dtm_atmosphere ), marsDtmFile_( marsDtmFile ), polarRadius_( polarRadius ){ } + + + std::string getMarsDtmFile( ){ return marsDtmFile_; } + + double getPolarRadius( ){ return polarRadius_; } + +private: + + std::string marsDtmFile_; + + double polarRadius_; +}; + + // AtmosphereSettings for defining an atmosphere with tabulated data from file. // //! @get_docstring(TabulatedAtmosphereSettings.__docstring__) class TabulatedAtmosphereSettings: public AtmosphereSettings @@ -936,6 +957,13 @@ inline std::shared_ptr< AtmosphereSettings > nrlmsise00AtmosphereSettings( return std::make_shared< NRLMSISE00AtmosphereSettings >( dataFile ); } +inline std::shared_ptr< AtmosphereSettings > marsDtmAtmosphereSettings( + const std::string& marsDtmFile, const double polarRadius ) +{ + return std::make_shared< MarsDtmAtmosphereSettings >( marsDtmFile, polarRadius ); +} + + typedef std::function< double( const double, const double, const double, const double ) > DensityFunction; //! @get_docstring(customConstantTemperatureAtmosphereSettings,0) inline std::shared_ptr< AtmosphereSettings > customConstantTemperatureAtmosphereSettings( diff --git a/src/simulation/environment_setup/createAtmosphereModel.cpp b/src/simulation/environment_setup/createAtmosphereModel.cpp index 84783fa885..cac6ed344d 100644 --- a/src/simulation/environment_setup/createAtmosphereModel.cpp +++ b/src/simulation/environment_setup/createAtmosphereModel.cpp @@ -13,6 +13,7 @@ #include "tudat/astro/aerodynamics/exponentialAtmosphere.h" #include "tudat/astro/aerodynamics/tabulatedAtmosphere.h" +#include "tudat/astro/aerodynamics/marsDtmAtmosphereModel.h" #if TUDAT_BUILD_WITH_NRLMSISE #include "tudat/astro/aerodynamics/nrlmsise00Atmosphere.h" #include "tudat/astro/aerodynamics/nrlmsise00InputFunctions.h" @@ -180,6 +181,21 @@ std::shared_ptr< aerodynamics::AtmosphereModel > createAtmosphereModel( } break; } + case mars_dtm_atmosphere: + { + std::shared_ptr< MarsDtmAtmosphereSettings > marsDtmAtmosphereSettings = + std::dynamic_pointer_cast< MarsDtmAtmosphereSettings >( atmosphereSettings ); + + if( marsDtmAtmosphereSettings == nullptr ) + { + throw std::runtime_error( "Error when creating Mars DTM atmosphere model, model settings are incompatible." ); + } + + // Create atmosphere model using NRLMISE00 input function + atmosphereModel = std::make_shared< aerodynamics::MarsDtmAtmosphereModel >( + marsDtmAtmosphereSettings->getPolarRadius( ), marsDtmAtmosphereSettings->getMarsDtmFile( ) ); + break; + } #if TUDAT_BUILD_WITH_NRLMSISE case nrlmsise00: { From 5a836e98365be6d136cd314877b6ba8ab76be1c7 Mon Sep 17 00:00:00 2001 From: Dominic Dirkx Date: Mon, 11 Dec 2023 12:27:49 +0100 Subject: [PATCH 06/25] Adding exceptions to unimplemented functions --- include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h b/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h index 70e3ea47d8..6526a198c7 100644 --- a/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h +++ b/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h @@ -140,7 +140,8 @@ class MarsDtmAtmosphereModel: public AtmosphereModel virtual double getPressure( const double altitude, const double longitude, const double latitude, const double time ) { - + throw std::runtime_error( "Error wen getting Mars DTM pressure; model not yet implemented." ); + return TUDAT_NAN; } //! Get local temperature. @@ -171,7 +172,8 @@ class MarsDtmAtmosphereModel: public AtmosphereModel virtual double getSpeedOfSound( const double altitude, const double longitude, const double latitude, const double time ) { - + throw std::runtime_error( "Error wen getting Mars DTM speed of sounds; model not yet implemented." ); + return TUDAT_NAN; } From 998a4669d1c9209a5eeba54207c2f3844a135b76 Mon Sep 17 00:00:00 2001 From: Dominic Dirkx Date: Fri, 15 Dec 2023 15:42:09 +0100 Subject: [PATCH 07/25] Quick first try at unit test --- .../unitTestMarsDtmAtmosphere.cpp | 137 +++++++++++++----- 1 file changed, 97 insertions(+), 40 deletions(-) diff --git a/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp b/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp index 42286697cd..cf9f5ef552 100644 --- a/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp +++ b/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp @@ -19,51 +19,108 @@ #include "iostream" #include -#include "tudat/astro/basic_astro/timeConversions.h" -#include "tudat/astro/aerodynamics/marsDtmAtmosphereModel.h" +#include "tudat/basics/testMacros.h" +#include "tudat/astro/aerodynamics/customAerodynamicCoefficientInterface.h" +#include "tudat/astro/aerodynamics/aerodynamicAcceleration.h" +#include "tudat/astro/reference_frames/aerodynamicAngleCalculator.h" +#include "tudat/simulation/propagation_setup/dynamicsSimulator.h" +#include "tudat/interface/spice/spiceEphemeris.h" +#include "tudat/interface/spice/spiceRotationalEphemeris.h" +#include "tudat/io/basicInputOutput.h" +#include "tudat/simulation/environment_setup/body.h" +#include "tudat/simulation/estimation_setup/createNumericalSimulator.h" +#include "tudat/simulation/environment_setup/defaultBodies.h" + using namespace tudat::aerodynamics; int main( ) { - //std::string filename = "/Users/ralkahal/Documents/PhD/atmodensitypds/dtm_mars"; - std::string filename = "/Users/ralkahal/Documents/PhD/Programs/atmodensitydtm/dtm_mars"; - MarsDtmAtmosphereModel atmosphereModel = MarsDtmAtmosphereModel(3376.78E3, filename); - //atmosphereModel.computelocalsolartime(0.0,16 ,12,2000, 0.0, 1.0697333); - std::cout<(altitude); - //double rho = atmosphereModel.getTotalDensity( alt_km, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000); - - // Write altitude and corresponding density to the file - // outputFile << alt_km << " " << rho << "\n"; - //} - - // Close the file when you are done - //outputFile.close(); - //std::cout << "Density computation and output written to file successfully." << std::endl; - std::cout << atmosphereModel.getTotalDensity( 1000.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000) << std::endl; - -// return 0; - //std::cout << atmosphereModel.getTotalDensity( 138.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000) <atmosphereSettings = marsDtmAtmosphereSettings( TODO ); + SystemOfBodies bodies = createSystemOfBodies( defaultBodySettings ); + + // Create vehicle object. + double vehicleMass = 5.0E3; + bodies.createEmptyBody( "Vehicle" ); + bodies.at( "Vehicle" )->setConstantBodyMass( vehicleMass ); + + // Set aerodynamic coefficients. + std::shared_ptr aerodynamicCoefficientSettings = + std::make_shared( + 2.0, 4.0, Eigen::Vector3d::Zero( ), Eigen::Vector3d::UnitX( ), Eigen::Vector3d::Zero( ), + negative_aerodynamic_frame_coefficients, negative_aerodynamic_frame_coefficients ); + bodies.at( "Vehicle" )->setAerodynamicCoefficientInterface( + createAerodynamicCoefficientInterface( aerodynamicCoefficientSettings, "Vehicle", bodies )); + + // Define acceleration model settings. + SelectedAccelerationMap accelerationMap; + std::vector bodiesToPropagate; + std::vector centralBodies; + std::map > > accelerationsOfVehicle; + accelerationsOfVehicle[ "Mars" ].push_back( std::make_shared( point_mass_gravity )); + accelerationsOfVehicle[ "Mars" ].push_back( std::make_shared( aerodynamic )); + accelerationMap[ "Vehicle" ] = accelerationsOfVehicle; + bodiesToPropagate.push_back( "Vehicle" ); + centralBodies.push_back( "Mars" ); + + // Set initial state + Eigen::Vector6d systemInitialState = TODO (add Cartesian initial state w.r.t. Mars ) + + // Create acceleration models and propagation settings. + basic_astrodynamics::AccelerationMap accelerationModelMap = createAccelerationModelsMap( + bodies, accelerationMap, bodiesToPropagate, centralBodies ); + + // Set variables to save + std::vector > dependentVariables; + dependentVariables.push_back( + std::make_shared( + altitude_dependent_variable, "Vehicle", "Earth" )); + dependentVariables.push_back( + std::make_shared( + "Vehicle", reference_frames::longitude_angle )); + dependentVariables.push_back( + std::make_shared( + "Vehicle", reference_frames::latitude_angle )); + dependentVariables.push_back( + std::make_shared( + local_density_dependent_variable, "Vehicle", "Earth" )); + dependentVariables.push_back( + std::make_shared( + aerodynamic, "Vehicle", "Earth", 0 )); + + // Set propagation/integration settings + std::shared_ptr terminationSettings = + std::make_shared( 1000.0 ); + std::shared_ptr > integratorSettings = + std::make_shared > + ( rungeKutta4, 0.0, 10.0 ); + std::shared_ptr > + translationalPropagatorSettings = + std::make_shared > + ( centralBodies, accelerationModelMap, bodiesToPropagate, systemInitialState, 0.0, + integratorSettings, terminationSettings, + cowell, dependentVariables ); + + // Create simulation object and propagate dynamics. + SingleArcDynamicsSimulator<> dynamicsSimulator( + bodies, translationalPropagatorSettings ); + } // //namespace tudat From 3e27b0d84526dba7baf28a8c819c9b7ec5bb8854 Mon Sep 17 00:00:00 2001 From: RivaAlkahal Date: Tue, 19 Dec 2023 11:35:06 +0100 Subject: [PATCH 08/25] test_round_ --- .../aerodynamics/marsDtmAtmosphereModel.cpp | 31 ++++++++----------- .../unitTestMarsDtmAtmosphere.cpp | 12 ++++--- 2 files changed, 20 insertions(+), 23 deletions(-) diff --git a/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp b/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp index 754ae1ce35..e5f3af6ff3 100644 --- a/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp +++ b/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp @@ -97,8 +97,8 @@ namespace aerodynamics coefficients_ = loadCoefficients( filename ); currentLegendrePolynomials_.resize( 6 ); currentLegendrePolynomials_[ 0 ] = 0.0; - std::cout << "Loaded Matrix values:" << std::endl; - std::cout << coefficients_[0][3] << "\t"; + //std::cout << "Loaded Matrix values:" << std::endl; + //std::cout << coefficients_[0][3] << "\t"; } // function to compute the local true solar time @@ -135,7 +135,7 @@ namespace aerodynamics //Areocentric Solar Longitude double Ls1 = omegaFMS + EoC; Ls = basic_mathematics::computeModulo( Ls1,360.0); - std::cout << "Ls: " << Ls << std::endl; + //std::cout << "Ls: " << Ls << std::endl; //Equation of Time double EoT = 2.861 * sin( unit_conversions::convertDegreesToRadians(2*Ls)) - 0.071 * sin( unit_conversions::convertDegreesToRadians(4*Ls)) @@ -145,12 +145,12 @@ namespace aerodynamics double MST = 24*((basic_astrodynamics::convertCalendarDateToJulianDay(year,month,day,hours,minutes,0.0) - 2451549.5)/1.0274912517 + 44796.0 - 0.0009626); //double MST //std::cout << "MST1: " << MST << std::endl; MST = basic_mathematics::computeModulo( MST,24.0);// Modulo MST - std::cout << "MST: " << MST << std::endl; + //std::cout << "MST: " << MST << std::endl; //Local Mean Solar Time double LMST = basic_mathematics::computeModulo( MST-longitude*(24.0/360.0),24.0); //Local True Solar Time double LTST = basic_mathematics::computeModulo(LMST + EoT*(24.0/360.0),24.0); - std::cout << "LTST: " << LTST << std::endl; + //std::cout << "LTST: " << LTST << std::endl; //Subsolar Longitude double subSolarLongitude = MST*(360.0/24.0) + EoT + 180.0; // return LTST in seconds. @@ -232,10 +232,10 @@ namespace aerodynamics using basic_mathematics::computeLegendrePolynomialExplicit; marsDate date2 = marsDate(year_, month_, day_, hours_ , minutes_, 0.0); double doy = date2.marsDayofYear(date2); //day of year - std::cout << "doy: " << doy << std::endl; + //std::cout << "doy: " << doy << std::endl; double t = computeLocalSolarTime(longitude, day_, month_, year_, hours_, minutes_); //seconds int F = 0.0; // flux parameter is set to 0 for now. - std::cout<<"coefficients_[14][indexg]" << coefficients_[14][indexg] << std::endl; + //std::cout<<"coefficients_[14][indexg]" << coefficients_[14][indexg] << std::endl; // Non-periodic term double NP = (coefficients_[1][indexg] * F) + coefficients_[30][indexg] * currentLegendrePolynomials_[1] + coefficients_[31][indexg] * currentLegendrePolynomials_[ 2 ] @@ -316,10 +316,6 @@ namespace aerodynamics double da41 = computeDustStorm( Ls ); double fpds = coefficients_[40][indexg]*da41 + coefficients_[69][indexg]*taus[3]; double Gl = fpds + f0 + PA + PSA + PD; - std::cout << "f0: " << f0 << std::endl; - std::cout << "PA: " << PA << std::endl; - std::cout << "PSA: " << PSA << std::endl; - std::cout << "PD: " << PD << std::endl; return Gl; } @@ -334,7 +330,7 @@ namespace aerodynamics { double T0 = coefficients_[0][indexg]; double Ti; - double Gl = computeGl(latitude, longitude, minutes_, hours_, day_ , month_, year_, indexg); + double Gl = computeGl_Subr(latitude, longitude, minutes_, hours_, day_ , month_, year_, indexg); Ti = T0*(1.0 + Gl); return Ti; } @@ -386,15 +382,14 @@ namespace aerodynamics rho0 = coefficients_[0][col]*mmass_[indexm]/avogadroConstant_; //g/cm3 //std::cout << "coefficients_[0][col]: " << coefficients_[0][col] << std::endl; double fi = heightDistributionFunction(altitude, latitude, longitude, minutes_, hours_, day_ , month_, year_, indexm); - if (col == 7) - std::cout << "fi CO2: " << fi << std::endl; + //if (col == 7) + // std::cout << "fi CO2: " << fi << std::endl; + + double gl = computeGl_Subr(latitude, longitude, minutes_, hours_, day_ , month_, year_, col); - double gl = computeGl(latitude, longitude, minutes_, hours_, day_ , month_, year_, col); - std::cout << "gl: " << gl << std::endl; - std::cout << "col" << col << std::endl; rho += rho0*fi*exp(gl); indexm +=1; - std::cout << "rho: " << indexm << " " << rho << std::endl; + //std::cout << "rho: " << indexm << " " << rho << std::endl; } return rho; } diff --git a/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp b/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp index cf9f5ef552..7ee0e9cc86 100644 --- a/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp +++ b/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp @@ -52,7 +52,8 @@ int main( ) // Create Earth object BodyListSettings defaultBodySettings = getDefaultBodySettings( { "Mars" } ); - defaultBodySettings.at( "Mars" )->atmosphereSettings = marsDtmAtmosphereSettings( TODO ); + std::string filename = "/Users/ralkahal/Documents/PhD/Programs/atmodensitydtm/dtm_mars"; + defaultBodySettings.at( "Mars" )->atmosphereSettings = marsDtmAtmosphereSettings( filename, 3378.0E3); SystemOfBodies bodies = createSystemOfBodies( defaultBodySettings ); // Create vehicle object. @@ -80,7 +81,8 @@ int main( ) centralBodies.push_back( "Mars" ); // Set initial state - Eigen::Vector6d systemInitialState = TODO (add Cartesian initial state w.r.t. Mars ) + Eigen::Vector6d systemInitialState; + systemInitialState // Create acceleration models and propagation settings. basic_astrodynamics::AccelerationMap accelerationModelMap = createAccelerationModelsMap( @@ -90,7 +92,7 @@ int main( ) std::vector > dependentVariables; dependentVariables.push_back( std::make_shared( - altitude_dependent_variable, "Vehicle", "Earth" )); + altitude_dependent_variable, "Vehicle", "Mars" )); dependentVariables.push_back( std::make_shared( "Vehicle", reference_frames::longitude_angle )); @@ -99,10 +101,10 @@ int main( ) "Vehicle", reference_frames::latitude_angle )); dependentVariables.push_back( std::make_shared( - local_density_dependent_variable, "Vehicle", "Earth" )); + local_density_dependent_variable, "Vehicle", "Mars" )); dependentVariables.push_back( std::make_shared( - aerodynamic, "Vehicle", "Earth", 0 )); + aerodynamic, "Vehicle", "Mars", 0 )); // Set propagation/integration settings std::shared_ptr terminationSettings = From dd182ddf1f78307a45f1c495fbfd9d7e6dd04603 Mon Sep 17 00:00:00 2001 From: RivaAlkahal Date: Tue, 19 Dec 2023 11:40:19 +0100 Subject: [PATCH 09/25] test_round_ --- tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp b/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp index 7ee0e9cc86..84a21a4091 100644 --- a/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp +++ b/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp @@ -82,7 +82,7 @@ int main( ) // Set initial state Eigen::Vector6d systemInitialState; - systemInitialState + systemInitialState << 3378.0E3, 3000E3, 4200E3, 0.0, 0.0, 0.0; // Create acceleration models and propagation settings. basic_astrodynamics::AccelerationMap accelerationModelMap = createAccelerationModelsMap( From 0fb3689c057391f5e7b111c9385215a0f8831124 Mon Sep 17 00:00:00 2001 From: RivaAlkahal Date: Tue, 19 Dec 2023 22:49:56 +0100 Subject: [PATCH 10/25] test_round_3 --- src/astro/aerodynamics/marsDtmAtmosphereModel.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp b/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp index e5f3af6ff3..9769b9c326 100644 --- a/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp +++ b/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp @@ -95,7 +95,7 @@ namespace aerodynamics mmass_( {44.01, 16.00, 28.0, 40.0, 28.0, 32.0, 4.0, 1.0, 2.0} )// molar mass of the species { coefficients_ = loadCoefficients( filename ); - currentLegendrePolynomials_.resize( 6 ); + currentLegendrePolynomials_.resize( 7 ); currentLegendrePolynomials_[ 0 ] = 0.0; //std::cout << "Loaded Matrix values:" << std::endl; //std::cout << coefficients_[0][3] << "\t"; From fb8cece2b80c0cbcfb42800514e58022a09b0f1b Mon Sep 17 00:00:00 2001 From: RivaAlkahal Date: Tue, 19 Dec 2023 23:02:01 +0100 Subject: [PATCH 11/25] test_round_3 --- src/astro/aerodynamics/marsDtmAtmosphereModel.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp b/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp index 9769b9c326..5c7922700f 100644 --- a/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp +++ b/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp @@ -391,7 +391,7 @@ namespace aerodynamics indexm +=1; //std::cout << "rho: " << indexm << " " << rho << std::endl; } - return rho; + return rho*1000; } From 64c85d8137a2b706e7eac2754ea1d5dccfdbb2d5 Mon Sep 17 00:00:00 2001 From: Dominic Dirkx Date: Mon, 15 Jan 2024 14:50:15 +0100 Subject: [PATCH 12/25] Added f107 interface to mars DTM model --- .../aerodynamics/marsDtmAtmosphereModel.h | 9 ++++++- .../environment_setup/createAtmosphereModel.h | 14 +++++++++-- .../aerodynamics/marsDtmAtmosphereModel.cpp | 11 +++++---- .../createAtmosphereModel.cpp | 24 ++++++++++++++++++- 4 files changed, 49 insertions(+), 9 deletions(-) diff --git a/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h b/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h index 6526a198c7..3133b96c0f 100644 --- a/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h +++ b/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h @@ -103,7 +103,10 @@ class MarsDtmAtmosphereModel: public AtmosphereModel public: - MarsDtmAtmosphereModel(const double polarRadius, const std::string &filename); + MarsDtmAtmosphereModel( + const double polarRadius, + const std::string &filename, + const std::function< double( const double ) > f107Function = [](const double){return 65.0;} ); //! Default destructor. /*! @@ -218,8 +221,11 @@ class MarsDtmAtmosphereModel: public AtmosphereModel const double latitude, const double time ); double polarRadius_; + std::basic_string, std::allocator> filename_; + std::function< double( const double ) > f107Function_; + std::vector< double > alpha_; std::vector< double > mmass_; @@ -227,6 +233,7 @@ class MarsDtmAtmosphereModel: public AtmosphereModel std::vector taus; double Ls; + double currentF107_; std::vector< double > currentLegendrePolynomials_; double currentGeopotentialAltitude_; double currentTemperature_138; diff --git a/include/tudat/simulation/environment_setup/createAtmosphereModel.h b/include/tudat/simulation/environment_setup/createAtmosphereModel.h index bcdaca076d..293fbbd48b 100644 --- a/include/tudat/simulation/environment_setup/createAtmosphereModel.h +++ b/include/tudat/simulation/environment_setup/createAtmosphereModel.h @@ -533,19 +533,29 @@ class MarsDtmAtmosphereSettings: public AtmosphereSettings { public: - MarsDtmAtmosphereSettings( const std::string& marsDtmFile, const double polarRadius ): - AtmosphereSettings( mars_dtm_atmosphere ), marsDtmFile_( marsDtmFile ), polarRadius_( polarRadius ){ } + MarsDtmAtmosphereSettings( const std::string& marsDtmFile, const double polarRadius, const std::string& spaceWeatherFile = "" ): + AtmosphereSettings( mars_dtm_atmosphere ), marsDtmFile_( marsDtmFile ), polarRadius_( polarRadius ), spaceWeatherFile_( spaceWeatherFile ){ } std::string getMarsDtmFile( ){ return marsDtmFile_; } double getPolarRadius( ){ return polarRadius_; } + std::string getSpaceWeatherFile( ){ return spaceWeatherFile_; } + private: std::string marsDtmFile_; double polarRadius_; + + // File containing space weather data. + /* + * File containing space weather data, as in https://celestrak.com/SpaceData/sw19571001.txt + */ + std::string spaceWeatherFile_; + + }; diff --git a/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp b/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp index 5c7922700f..52fc4c2507 100644 --- a/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp +++ b/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp @@ -20,6 +20,7 @@ namespace aerodynamics { // Compute the hash key?? basic_astrodynamics::DateTime currentDateTime_ = basic_astrodynamics::getCalendarDateFromTime( time ); + currentF107_ = f107Function_( time ); currentDensity_ = getTotalDensity( altitude, longitude, latitude, currentDateTime_.getMinute(), currentDateTime_.getHour(), currentDateTime_.getDay(), currentDateTime_.getMonth(), currentDateTime_.getYear() ); @@ -88,9 +89,11 @@ namespace aerodynamics return dayOfYear; } - MarsDtmAtmosphereModel::MarsDtmAtmosphereModel(const double polarRadius, const std::string &filename) : + MarsDtmAtmosphereModel::MarsDtmAtmosphereModel(const double polarRadius, const std::string &filename, + const std::function< double( const double ) > f107Function ) : polarRadius_( polarRadius ), // polar radius of Mars filename_ ( filename ), // file name of the coefficients + f107Function_( f107Function ), alpha_( {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.38, -0.40, 0.0}), // thermal diffusion coefficients mmass_( {44.01, 16.00, 28.0, 40.0, 28.0, 32.0, 4.0, 1.0, 2.0} )// molar mass of the species { @@ -234,7 +237,7 @@ namespace aerodynamics double doy = date2.marsDayofYear(date2); //day of year //std::cout << "doy: " << doy << std::endl; double t = computeLocalSolarTime(longitude, day_, month_, year_, hours_, minutes_); //seconds - int F = 0.0; // flux parameter is set to 0 for now. + double F = currentF107_ - 65.0; //std::cout<<"coefficients_[14][indexg]" << coefficients_[14][indexg] << std::endl; // Non-periodic term double NP = (coefficients_[1][indexg] * F) + coefficients_[30][indexg] * currentLegendrePolynomials_[1] @@ -279,9 +282,7 @@ namespace aerodynamics double sin2h = 2.0*sin(hl0)*cos(hl0); //flux terms: double ff0 = 0.0; - double F107 = 65.0; - double Fbar = 65.0; - double F = Fbar - F107; + double F = currentF107_ - 65.0; double F2 = F*F; double f0 = coefficients_[4][indexg]*F + coefficients_[39][indexg] * F2; double f1f = 1.0 + f0*ff0; // coupling terms diff --git a/src/simulation/environment_setup/createAtmosphereModel.cpp b/src/simulation/environment_setup/createAtmosphereModel.cpp index cac6ed344d..0a4f2ed52c 100644 --- a/src/simulation/environment_setup/createAtmosphereModel.cpp +++ b/src/simulation/environment_setup/createAtmosphereModel.cpp @@ -191,9 +191,31 @@ std::shared_ptr< aerodynamics::AtmosphereModel > createAtmosphereModel( throw std::runtime_error( "Error when creating Mars DTM atmosphere model, model settings are incompatible." ); } + std::string spaceWeatherFilePath; + + if( marsDtmAtmosphereSettings->getSpaceWeatherFile( ) == "" ) + { + // Use default space weather file stored in tudatBundle. + spaceWeatherFilePath = paths::getSpaceWeatherDataPath( ) + "/sw19571001.txt"; + } + else + { + // Use space weather file specified by user. + spaceWeatherFilePath = marsDtmAtmosphereSettings->getSpaceWeatherFile( ); + } + + tudat::input_output::solar_activity::SolarActivityDataMap solarActivityData = + tudat::input_output::solar_activity::readSolarActivityData( spaceWeatherFilePath ); + std::shared_ptr< input_output::solar_activity::SolarActivityContainer > solarActivityContainer = + std::make_shared< input_output::solar_activity::SolarActivityContainer >( solarActivityData ); + std::function< double( const double ) > f107Function = [=](const double time) + { + return solarActivityContainer->getSolarActivityData( time )->solarRadioFlux107Observed / 2.25; + }; + // Create atmosphere model using NRLMISE00 input function atmosphereModel = std::make_shared< aerodynamics::MarsDtmAtmosphereModel >( - marsDtmAtmosphereSettings->getPolarRadius( ), marsDtmAtmosphereSettings->getMarsDtmFile( ) ); + marsDtmAtmosphereSettings->getPolarRadius( ), marsDtmAtmosphereSettings->getMarsDtmFile( ), f107Function ); break; } #if TUDAT_BUILD_WITH_NRLMSISE From 6f428fbbf72e522d32f939ea843133a4e4e3a14e Mon Sep 17 00:00:00 2001 From: Riva Date: Thu, 25 Jan 2024 16:44:34 +0100 Subject: [PATCH 13/25] test_mars_dtm added --- CMakeLists.txt | 2 +- .../aerodynamics/marsDtmAtmosphereModel.h | 7 + .../aerodynamics/marsDtmAtmosphereModel.cpp | 1 + .../unitTestMarsDtmAtmosphere.cpp | 2 +- tests_riva/CMakeLists.txt | 4 + tests_riva/TestMarsDtm.cpp | 120 ++++++++++++++++++ 6 files changed, 134 insertions(+), 2 deletions(-) create mode 100644 tests_riva/CMakeLists.txt create mode 100644 tests_riva/TestMarsDtm.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index ba7e2f837e..28f0d24b55 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -347,7 +347,7 @@ add_subdirectory(src/simulation) add_subdirectory(src/io) add_subdirectory(src/utils) #add_subdirectory(examples) -#add_subdirectory(applications) +add_subdirectory(tests_riva) if (TUDAT_BUILD_TESTS) diff --git a/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h b/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h index 3133b96c0f..de7eb7cf4a 100644 --- a/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h +++ b/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h @@ -213,6 +213,13 @@ class MarsDtmAtmosphereModel: public AtmosphereModel double computeDustStorm(const double Ls); + double getSolarLongitude() const { + return Ls; + } + + double getSolarFluxIndex() const{ + return currentF107_; + } protected: private: diff --git a/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp b/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp index 52fc4c2507..f1af6f5d6b 100644 --- a/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp +++ b/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp @@ -284,6 +284,7 @@ namespace aerodynamics double ff0 = 0.0; double F = currentF107_ - 65.0; double F2 = F*F; + std::cout << "F: " << F << std::endl; double f0 = coefficients_[4][indexg]*F + coefficients_[39][indexg] * F2; double f1f = 1.0 + f0*ff0; // coupling terms //flux + latitude terms diff --git a/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp b/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp index 84a21a4091..1812b6910b 100644 --- a/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp +++ b/tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere.cpp @@ -52,7 +52,7 @@ int main( ) // Create Earth object BodyListSettings defaultBodySettings = getDefaultBodySettings( { "Mars" } ); - std::string filename = "/Users/ralkahal/Documents/PhD/Programs/atmodensitydtm/dtm_mars"; + std::string filename ="/Users/ralkahal/OneDrive - Delft University of Technology/PhD/Programs/atmodensitydtm/dtm_mars";; defaultBodySettings.at( "Mars" )->atmosphereSettings = marsDtmAtmosphereSettings( filename, 3378.0E3); SystemOfBodies bodies = createSystemOfBodies( defaultBodySettings ); diff --git a/tests_riva/CMakeLists.txt b/tests_riva/CMakeLists.txt new file mode 100644 index 0000000000..ec23312573 --- /dev/null +++ b/tests_riva/CMakeLists.txt @@ -0,0 +1,4 @@ +TUDAT_ADD_EXECUTABLE(application_test_mars_dtm + "TestMarsDtm.cpp" + ${Tudat_ESTIMATION_LIBRARIES} +) diff --git a/tests_riva/TestMarsDtm.cpp b/tests_riva/TestMarsDtm.cpp new file mode 100644 index 0000000000..2743c0a81c --- /dev/null +++ b/tests_riva/TestMarsDtm.cpp @@ -0,0 +1,120 @@ +/* Copyright (c) 2010-2019, 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 + * a copy of the license with this file. If not, please or visit: + * http://tudat.tudelft.nl/LICENSE. + * + * References + * + */ + +//#define BOOST_TEST_DYN_LINK +//#define BOOST_TEST_MAIN + +#include +#include "fstream" +#include "iostream" + +#include + +#include "tudat/basics/testMacros.h" +#include "tudat/astro/aerodynamics/customAerodynamicCoefficientInterface.h" +#include "tudat/astro/aerodynamics/aerodynamicAcceleration.h" +#include "tudat/astro/reference_frames/aerodynamicAngleCalculator.h" +#include "tudat/simulation/propagation_setup/dynamicsSimulator.h" +#include "tudat/interface/spice/spiceEphemeris.h" +#include "tudat/interface/spice/spiceRotationalEphemeris.h" +#include "tudat/io/basicInputOutput.h" +#include "tudat/simulation/environment_setup/body.h" +#include "tudat/simulation/estimation_setup/createNumericalSimulator.h" +#include "tudat/simulation/environment_setup/defaultBodies.h" + +using namespace tudat::aerodynamics; + + +#include "tudat/astro/basic_astro/timeConversions.h" + +#include "tudat/astro/aerodynamics/marsDtmAtmosphereModel.h" + +int main( ) { + //std::string filename = "/Users/ralkahal/Documents/PhD/atmodensitypds/dtm_mars"; + std::string filename = "/Users/ralkahal/OneDrive - Delft University of Technology/PhD/Programs/atmodensitydtm/dtm_mars"; + MarsDtmAtmosphereModel atmosphereModel = MarsDtmAtmosphereModel(3376.78E3, filename); + //atmosphereModel.computelocalsolartime(0.0,16 ,12,2000, 0.0, 1.0697333); + //std::cout<(altitude); + + double rho = atmosphereModel.getTotalDensity(alt_km, 0.0, 0.0, 1.0697333, 0.0, 16, 12, 2000); + + // Write altitude and corresponding density to the file + outputFile << alt_km << " " << rho << "\n"; + } + */ + //for (int altitude = 138E3; altitude <= 1000E3; altitude += 10E3) { + // double alt_km = static_cast(altitude); + int alt_km = 138E3; + for (int time = 0; time <= 5 * 24 * 60 * 60; time += 60) { + + int minutes = time % (60 * 60); + int hours = (time / (60 * 60)) % 24; + int days = ((time / (60 * 60 * 24)) % 31) + 1; // Adding 1 to ensure days are between 1 and 31 + int months = ((time / (60 * 60 * 24 * 31)) % 12) + 1; // Adding 1 to ensure months are between 1 and 12 + int years = (time / (60 * 60 * 24 * 31 * 12)) + 2000; + + double rho = atmosphereModel.getTotalDensity(alt_km, 0.0, 0.0, minutes, hours, days, months, years); + double Ls = atmosphereModel.getSolarLongitude(); + double currentF107 = atmosphereModel.getSolarFluxIndex(); + // Write altitude and corresponding density to the file + std::cout << alt_km << "altitude" << rho << "density" << Ls << "Ls" << currentF107 << "F107" << std::endl; + outputFile << alt_km << " " << rho << "\n"; + + //} + }// Close the file when you are done + outputFile.close(); + std::cout << "Density computation and output written to file successfully." << std::endl; + //std::cout << atmosphereModel.getTotalDensity( 1000.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000) << std::endl; + +// return 0; + //std::cout << atmosphereModel.getTotalDensity( 138.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000) < Date: Mon, 5 Feb 2024 13:55:47 +0100 Subject: [PATCH 14/25] negative sigma bug fixed --- .../aerodynamics/marsDtmAtmosphereModel.h | 28 +-- .../propagation_setup/.idea/.gitignore | 8 + .../propagation_setup/.idea/modules.xml | 8 + .../.idea/propagation_setup.iml | 8 + .../propagation_setup/.idea/vcs.xml | 6 + .../propagation_setup/propagationOutput.h | 7 + .../propagationOutputSettings.h | 3 +- .../aerodynamics/marsDtmAtmosphereModel.cpp | 117 ++++++---- .../createEnvironmentUpdater.cpp | 5 +- .../propagation_setup/propagationOutput.cpp | 3 + .../propagationOutputSettings.cpp | 3 + tests_riva/CMakeLists.txt | 4 + tests_riva/TestMarsDtm.cpp | 127 +++++++---- tests_riva/TestMarsDtm_.cpp | 204 ++++++++++++++++++ 14 files changed, 432 insertions(+), 99 deletions(-) create mode 100644 include/tudat/simulation/propagation_setup/.idea/.gitignore create mode 100644 include/tudat/simulation/propagation_setup/.idea/modules.xml create mode 100644 include/tudat/simulation/propagation_setup/.idea/propagation_setup.iml create mode 100644 include/tudat/simulation/propagation_setup/.idea/vcs.xml create mode 100644 tests_riva/TestMarsDtm_.cpp diff --git a/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h b/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h index de7eb7cf4a..b57401ff35 100644 --- a/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h +++ b/include/tudat/astro/aerodynamics/marsDtmAtmosphereModel.h @@ -97,6 +97,8 @@ namespace aerodynamics marsDate(2028, 8, 17 , 0, 0, 0) }; +// Function to compute the solar longitude of Mars +std::tuple computeSolarLongitude(const double longitude, const int day ,const int month,const int year, const int hours, const int minutes); class MarsDtmAtmosphereModel: public AtmosphereModel { @@ -184,42 +186,42 @@ class MarsDtmAtmosphereModel: public AtmosphereModel const double Omega = 2.0*mathematical_constants::PI/686.98;//(686.98*24.63);//rad/h//2*M_PI/(659355072); //rad/s const double omega = 2.0*mathematical_constants::PI/88668;//24.63;//2*M_PI/(88668); // rad/s - double computeLocalSolarTime(const double longitude, const int day, const int month, const int year, const double hours, - const double minutes); + double computeLocalSolarTime(const double longitude, const int day, const int month, const int year, const int hours, + const int minutes); void updateLegendrePolynomials( const double latitude ); - double computeGl(const double latitude, const double longitude, const double minutes, const double hours_, const int day_, + double computeGl(const double latitude, const double longitude, const int minutes, const int hours_, const int day_, const int month_, const int year_, const int indexg); double computeGeopotentialAltitude( const double altitude ); - double computeCurrentTemperature( const double latitude, const double longitude, const double minutes_, const double hours_, const int day_ ,const int month_,const int year_, const int indexg); + double computeCurrentTemperature( const double latitude, const double longitude, const int minutes_, const int hours_, const int day_ ,const int month_,const int year_, const int indexg); - double computeGamma(const double latitude, const double longitude, const double minutes_, const double hours_, const int day_ ,const int month_,const int year_, const int indexm); + double computeGamma(const double latitude, const double longitude, const int minutes_, const int hours_, const int day_ ,const int month_,const int year_, const int indexm); double heightDistributionFunction(const double altitude, const double latitude, - const double longitude, const double minutes_, const double hours_, const int day_ ,const int month_,const int year_,const int indexm); + const double longitude, const int minutes_, const int hours_, const int day_ ,const int month_,const int year_,const int indexm); double getTotalDensity( const double altitude, const double latitude, - const double longitude, const double minutes_, const double hours_, const int day_ , const int month_, const int year_); + const double longitude, const int minutes_, const int hours_, const int day_ , const int month_, const int year_); double getCurrentLegendrePolynomial( const int degree ) { return currentLegendrePolynomials_.at( degree ); } - double computeGl_Subr(const double latitude, const double longitude, const double minutes_, const double hours_, + double computeGl_Subr(const double latitude, const double longitude, const int minutes_, const int hours_, const int day_, const int month_, const int year_, const int indexg); void defineDustPars(const int rows); double computeDustStorm(const double Ls); - double getSolarLongitude() const { - return Ls; - } - double getSolarFluxIndex() const{ return currentF107_; } + double getSolarLongitude() const{ + std::cout << "Solar Longitude: " << Ls_ << std::endl; + return Ls_; + } protected: private: @@ -239,7 +241,7 @@ class MarsDtmAtmosphereModel: public AtmosphereModel std::vector> coefficients_; std::vector taus; - double Ls; + double Ls_; double currentF107_; std::vector< double > currentLegendrePolynomials_; double currentGeopotentialAltitude_; diff --git a/include/tudat/simulation/propagation_setup/.idea/.gitignore b/include/tudat/simulation/propagation_setup/.idea/.gitignore new file mode 100644 index 0000000000..13566b81b0 --- /dev/null +++ b/include/tudat/simulation/propagation_setup/.idea/.gitignore @@ -0,0 +1,8 @@ +# Default ignored files +/shelf/ +/workspace.xml +# Editor-based HTTP Client requests +/httpRequests/ +# Datasource local storage ignored files +/dataSources/ +/dataSources.local.xml diff --git a/include/tudat/simulation/propagation_setup/.idea/modules.xml b/include/tudat/simulation/propagation_setup/.idea/modules.xml new file mode 100644 index 0000000000..a78a4f7fd7 --- /dev/null +++ b/include/tudat/simulation/propagation_setup/.idea/modules.xml @@ -0,0 +1,8 @@ + + + + + + + + \ No newline at end of file diff --git a/include/tudat/simulation/propagation_setup/.idea/propagation_setup.iml b/include/tudat/simulation/propagation_setup/.idea/propagation_setup.iml new file mode 100644 index 0000000000..bc2cd87409 --- /dev/null +++ b/include/tudat/simulation/propagation_setup/.idea/propagation_setup.iml @@ -0,0 +1,8 @@ + + + + + + + + \ No newline at end of file diff --git a/include/tudat/simulation/propagation_setup/.idea/vcs.xml b/include/tudat/simulation/propagation_setup/.idea/vcs.xml new file mode 100644 index 0000000000..4fce1d86b4 --- /dev/null +++ b/include/tudat/simulation/propagation_setup/.idea/vcs.xml @@ -0,0 +1,6 @@ + + + + + + \ No newline at end of file diff --git a/include/tudat/simulation/propagation_setup/propagationOutput.h b/include/tudat/simulation/propagation_setup/propagationOutput.h index accd25da1d..cc077b4dbf 100644 --- a/include/tudat/simulation/propagation_setup/propagationOutput.h +++ b/include/tudat/simulation/propagation_setup/propagationOutput.h @@ -25,6 +25,7 @@ #include "tudat/simulation/propagation_setup/propagationSettings.h" #include "tudat/simulation/environment_setup/createFlightConditions.h" #include "tudat/math/basic/rotationRepresentations.h" +#include "tudat/astro/aerodynamics/marsDtmAtmosphereModel.h" namespace tudat { @@ -2514,6 +2515,12 @@ std::function< double( ) > getDoubleDependentVariableFunction( } break; } + case solar_longitude: + { + variableFunction = std::bind( &::tudat::aerodynamics::MarsDtmAtmosphereModel::getSolarLongitude, + std::dynamic_pointer_cast< aerodynamics::MarsDtmAtmosphereModel >( bodies.at( bodyWithProperty )->getAtmosphereModel( ) ) ); + break; + } default: std::string errorMessage = "Error, did not recognize double dependent variable type when making variable function: " + diff --git a/include/tudat/simulation/propagation_setup/propagationOutputSettings.h b/include/tudat/simulation/propagation_setup/propagationOutputSettings.h index bd98a67475..fa3f33b703 100644 --- a/include/tudat/simulation/propagation_setup/propagationOutputSettings.h +++ b/include/tudat/simulation/propagation_setup/propagationOutputSettings.h @@ -136,7 +136,8 @@ enum PropagationDependentVariables received_irradiance = 63, received_fraction = 64, visible_and_emitting_source_panel_count = 65, - visible_source_area = 66 + visible_source_area = 66, + solar_longitude = 67 }; // Functional base class for defining settings for dependent variables that are to be saved during propagation diff --git a/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp b/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp index f1af6f5d6b..b58984ac2e 100644 --- a/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp +++ b/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp @@ -7,7 +7,7 @@ #include "tudat/math/basic/legendrePolynomials.h" #include "tudat/astro/basic_astro/unitConversions.h" #include "tudat/astro/basic_astro/celestialBodyConstants.h" - +#include namespace tudat { @@ -89,6 +89,45 @@ namespace aerodynamics return dayOfYear; } + std::tuple computeSolarLongitude(const double longitude, const int day ,const int month,const int year, const int hours, const int minutes){ + + //Day since J2000 epoch in Terrestrial Time (or TAI) + double timeElaspsedJD = basic_astrodynamics::convertCalendarDateToJulianDaysSinceEpoch( + year, month, day, hours, minutes, 0.0, basic_astrodynamics::JULIAN_DAY_ON_J2000); + //std::cout << "timeElaspsedJD: " << timeElaspsedJD << std::endl; + //Mars orbital parameters: + //Mars Mean Anomaly + double meanAnomaly = 19.3870 + 0.52402073 * timeElaspsedJD; //degrees + //std::cout << "MeanAnomaly: " << MeanAnomaly << std::endl; + //Angle of Fiction Mean Sun + double omegaFMS = 270.3863 + 0.52403840 * timeElaspsedJD; + //std::cout << "OmegaFMS: " << OmegaFMS << std::endl; + //Perturbers + double PBS = + 0.0071 * cos(unit_conversions::convertDegreesToRadians(0.985626 * timeElaspsedJD / 2.2353 + 49.409)) + + 0.0057 * cos(unit_conversions::convertDegreesToRadians(0.985626 * timeElaspsedJD / 2.7543 + 168.173)) + + 0.0039 * cos(unit_conversions::convertDegreesToRadians(0.985626 * timeElaspsedJD / 1.1177 + 191.837)) + + 0.0037 * cos(unit_conversions::convertDegreesToRadians(0.985626 * timeElaspsedJD / 15.7866 + 21.736)) + + 0.0021 * cos(unit_conversions::convertDegreesToRadians(0.985626 * timeElaspsedJD / 2.1354 + 15.704)) + + 0.0020 * cos(unit_conversions::convertDegreesToRadians(0.985626 * timeElaspsedJD / 2.4694 + 95.528)) + + 0.0018 * cos(unit_conversions::convertDegreesToRadians(0.985626 * timeElaspsedJD / 32.8493 + 49.095)); + //std::cout << "PBS: " << PBS << std::endl; + //Equation of Center + double EoC = (10.691 + 3.0E-7 * timeElaspsedJD) * sin(unit_conversions::convertDegreesToRadians(meanAnomaly)) + + 0.623 * sin(unit_conversions::convertDegreesToRadians(2 * meanAnomaly)) + + 0.050 * sin(unit_conversions::convertDegreesToRadians(3 * meanAnomaly)) + + 0.005 * sin(unit_conversions::convertDegreesToRadians(4 * meanAnomaly)) + + 0.0005 * sin(unit_conversions::convertDegreesToRadians(5 * meanAnomaly)) + + PBS;//Equation of Center + //std::cout << "EoC: " << EoC << std::endl; + //Areocentric Solar Longitude + double Ls1 = omegaFMS + EoC; + double Ls = basic_mathematics::computeModulo(Ls1, 360.0); + return std::make_tuple(Ls, EoC); + } + + + MarsDtmAtmosphereModel::MarsDtmAtmosphereModel(const double polarRadius, const std::string &filename, const std::function< double( const double ) > f107Function ) : polarRadius_( polarRadius ), // polar radius of Mars @@ -100,49 +139,24 @@ namespace aerodynamics coefficients_ = loadCoefficients( filename ); currentLegendrePolynomials_.resize( 7 ); currentLegendrePolynomials_[ 0 ] = 0.0; + Ls_ = 0.0; //std::cout << "Loaded Matrix values:" << std::endl; //std::cout << coefficients_[0][3] << "\t"; } // function to compute the local true solar time // equations from: https://www.giss.nasa.gov/tools/mars24/help/algorithm.html - double MarsDtmAtmosphereModel::computeLocalSolarTime(const double longitude, const int day ,const int month,const int year, const double hours, const double minutes) + double MarsDtmAtmosphereModel::computeLocalSolarTime(const double longitude, const int day ,const int month,const int year, const int hours, const int minutes) { - //Day since J2000 epoch in Terrestrial Time (or TAI) - double timeElaspsedJD =basic_astrodynamics::convertCalendarDateToJulianDaysSinceEpoch( - year, month, day, hours, minutes, 0.0, basic_astrodynamics::JULIAN_DAY_ON_J2000 ); - //std::cout << "timeElaspsedJD: " << timeElaspsedJD << std::endl; - //Mars orbital parameters: - //Mars Mean Anomaly - double meanAnomaly = 19.3870 + 0.52402073 * timeElaspsedJD; //degrees - //std::cout << "MeanAnomaly: " << MeanAnomaly << std::endl; - //Angle of Fiction Mean Sun - double omegaFMS = 270.3863 + 0.52403840 * timeElaspsedJD; - //std::cout << "OmegaFMS: " << OmegaFMS << std::endl; - //Perturbers - double PBS = 0.0071 * cos( unit_conversions::convertDegreesToRadians(0.985626*timeElaspsedJD/2.2353 + 49.409)) - + 0.0057 * cos( unit_conversions::convertDegreesToRadians(0.985626*timeElaspsedJD/2.7543 + 168.173)) - + 0.0039 * cos( unit_conversions::convertDegreesToRadians(0.985626*timeElaspsedJD/1.1177 + 191.837)) - + 0.0037 * cos( unit_conversions::convertDegreesToRadians(0.985626*timeElaspsedJD/15.7866 + 21.736)) - + 0.0021 * cos( unit_conversions::convertDegreesToRadians(0.985626*timeElaspsedJD/2.1354 + 15.704)) - + 0.0020 * cos( unit_conversions::convertDegreesToRadians(0.985626*timeElaspsedJD/2.4694 + 95.528)) - + 0.0018 * cos( unit_conversions::convertDegreesToRadians(0.985626*timeElaspsedJD/32.8493 + 49.095)); - //std::cout << "PBS: " << PBS << std::endl; - //Equation of Center - double EoC = (10.691 + 3.0E-7 * timeElaspsedJD) * sin( unit_conversions::convertDegreesToRadians(meanAnomaly)) - + 0.623 * sin( unit_conversions::convertDegreesToRadians(2*meanAnomaly)) - + 0.050 * sin( unit_conversions::convertDegreesToRadians(3*meanAnomaly)) - + 0.005 * sin( unit_conversions::convertDegreesToRadians(4*meanAnomaly)) - + 0.0005 * sin( unit_conversions::convertDegreesToRadians(5*meanAnomaly)) + PBS;//Equation of Center - //std::cout << "EoC: " << EoC << std::endl; - //Areocentric Solar Longitude - double Ls1 = omegaFMS + EoC; - Ls = basic_mathematics::computeModulo( Ls1,360.0); + + auto outResults = computeSolarLongitude(longitude, day, month, year, hours, minutes); + double EoC = std::get<1>(outResults); + Ls_ = std::get<0>(outResults); //std::cout << "Ls: " << Ls << std::endl; //Equation of Time - double EoT = 2.861 * sin( unit_conversions::convertDegreesToRadians(2*Ls)) - - 0.071 * sin( unit_conversions::convertDegreesToRadians(4*Ls)) - + 0.002 * sin( unit_conversions::convertDegreesToRadians(6*Ls)) - EoC; // Equation of Time + double EoT = 2.861 * sin( unit_conversions::convertDegreesToRadians(2*Ls_)) + - 0.071 * sin( unit_conversions::convertDegreesToRadians(4*Ls_)) + + 0.002 * sin( unit_conversions::convertDegreesToRadians(6*Ls_)) - EoC; // Equation of Time //std::cout << "EoT: " << EoT << std::endl; //Mean Solar Time at Mars's prime merdian, i.e., Airy Mean Time double MST = 24*((basic_astrodynamics::convertCalendarDateToJulianDay(year,month,day,hours,minutes,0.0) - 2451549.5)/1.0274912517 + 44796.0 - 0.0009626); //double MST @@ -229,7 +243,7 @@ namespace aerodynamics } // Function to compute the spherical harmonic expansions defined as G(l) as in the paper (Bruinsma and Lemoine, 2002) - double MarsDtmAtmosphereModel::computeGl(const double latitude, const double longitude, const double minutes_, const double hours_, const int day_ ,const int month_,const int year_, const int indexg) + double MarsDtmAtmosphereModel::computeGl(const double latitude, const double longitude, const int minutes_, const int hours_, const int day_ ,const int month_,const int year_, const int indexg) { updateLegendrePolynomials( latitude ); using basic_mathematics::computeLegendrePolynomialExplicit; @@ -263,7 +277,7 @@ namespace aerodynamics } //Function to compute the spherical harmonic expansions defined as G(l) as in the subroutine provided by Bruinsma - double MarsDtmAtmosphereModel::computeGl_Subr(const double latitude, const double longitude, const double minutes_, const double hours_, const int day_ ,const int month_,const int year_, const int indexg) + double MarsDtmAtmosphereModel::computeGl_Subr(const double latitude, const double longitude, const int minutes_, const int hours_, const int day_ ,const int month_,const int year_, const int indexg) { //Get the zonal Legendre polynomials for the current latitude updateLegendrePolynomials( latitude ); @@ -273,6 +287,7 @@ namespace aerodynamics //Get the day of year marsDate date2 = marsDate(year_, month_, day_, hours_ , minutes_, 0.0); double doy = date2.marsDayofYear(date2); + //std::cout << "doy: " << doy << std::endl; //Get the local solar time double t = computeLocalSolarTime(longitude, day_, month_, year_, hours_, minutes_); //hours //Get the local solar time in radians + pi. @@ -284,7 +299,7 @@ namespace aerodynamics double ff0 = 0.0; double F = currentF107_ - 65.0; double F2 = F*F; - std::cout << "F: " << F << std::endl; + //std::cout << "F: " << F << std::endl; double f0 = coefficients_[4][indexg]*F + coefficients_[39][indexg] * F2; double f1f = 1.0 + f0*ff0; // coupling terms //flux + latitude terms @@ -315,7 +330,7 @@ namespace aerodynamics + coefficients_[34][indexg] * computeLegendrePolynomialExplicit(4,2,sin(latitude)) * (cos2h) + coefficients_[35][indexg] * computeLegendrePolynomialExplicit(4,2,sin(latitude)) * sin2h + coefficients_[36][indexg] * computeLegendrePolynomialExplicit(2,2,sin(latitude)) * (cos2h) * F + coefficients_[37][indexg] * computeLegendrePolynomialExplicit(2,2,sin(latitude)) * sin2h * F; //dust terms - double da41 = computeDustStorm( Ls ); + double da41 = computeDustStorm( Ls_ ); double fpds = coefficients_[40][indexg]*da41 + coefficients_[69][indexg]*taus[3]; double Gl = fpds + f0 + PA + PSA + PD; return Gl; @@ -328,16 +343,18 @@ namespace aerodynamics } // Function to compute the current temperature (exospheric, 138 km, and partial temperatures) - double MarsDtmAtmosphereModel::computeCurrentTemperature( const double latitude, const double longitude, const double minutes_, const double hours_, const int day_ ,const int month_,const int year_, const int indexg) + double MarsDtmAtmosphereModel::computeCurrentTemperature( const double latitude, const double longitude, const int minutes_, const int hours_, const int day_ ,const int month_,const int year_, const int indexg) { double T0 = coefficients_[0][indexg]; double Ti; double Gl = computeGl_Subr(latitude, longitude, minutes_, hours_, day_ , month_, year_, indexg); + //std::cout << "Gl: " << Gl << std::endl; + //std::cout << "T0: " << T0 << std::endl; Ti = T0*(1.0 + Gl); return Ti; } // Function to compute gamma parameter - double MarsDtmAtmosphereModel::computeGamma( const double latitude, const double longitude, const double minutes_, const double hours_, const int day_ ,const int month_,const int year_, const int indexm) + double MarsDtmAtmosphereModel::computeGamma( const double latitude, const double longitude, const int minutes_, const int hours_, const int day_ ,const int month_,const int year_, const int indexm) { double universalGasConstant = 8.314;//physical_constants::MOLAR_GAS_CONSTANT; //J/mol/K //kg m^2 / s^2 / K / mol //double g138 =celestial_body_constants::MARS_GRAVITATIONAL_PARAMETER/((138.0E3+3376.78E3)*(138.0E3+3376.78E3)); //m/s2 @@ -346,6 +363,10 @@ namespace aerodynamics currentTemperature_138 = computeCurrentTemperature( latitude, longitude, minutes_, hours_, day_ , month_, year_, 3); //K currentTemperature_inf = computeCurrentTemperature( latitude, longitude, minutes_, hours_, day_ , month_, year_, 1); //K currentdTemperature_ = computeCurrentTemperature( latitude, longitude, minutes_, hours_, day_ , month_, year_, 5); //K + double tmt = currentTemperature_inf - currentTemperature_138; + if (tmt <= 0.0) { + currentTemperature_138 = 0.99 * currentTemperature_inf; + } sigma = currentdTemperature_/(currentTemperature_inf-currentTemperature_138); //std::cout << "sigma: " << sigma << std::endl; double gamma = mmass_[indexm] * g138 / ( universalGasConstant * sigma * currentTemperature_inf); //km^-1 @@ -361,12 +382,17 @@ namespace aerodynamics // Function to compute the height distribution function double MarsDtmAtmosphereModel::heightDistributionFunction(const double altitude, const double latitude, - const double longitude, const double minutes_, const double hours_, const int day_ ,const int month_,const int year_,const int indexm) + const double longitude, const int minutes_, const int hours_, const int day_ ,const int month_,const int year_,const int indexm) { currentGeopotentialAltitude_ = computeGeopotentialAltitude( altitude ); //km + //std::cout << "currentGeopotentialAltitude_: " << currentGeopotentialAltitude_ << std::endl; double gamma= computeGamma( latitude, longitude, minutes_, hours_, day_ , month_, year_, indexm); //km^-1 currentTemperature_z = currentTemperature_inf - (currentTemperature_inf - currentTemperature_138)* exp(-sigma*currentGeopotentialAltitude_); - //std::cout << "Tz: " << Tz << std::endl; + // std::cout << "T138: " << currentTemperature_138 << std::endl; + // std::cout << "Tinf: " << currentTemperature_inf << std::endl; + // std::cout << "exp(-sigma*currentGeopotentialAltitude_): " << exp(-sigma*currentGeopotentialAltitude_) << std::endl; + // std::cout << "Tz: " << currentTemperature_z << std::endl; + double fi = pow((currentTemperature_138/currentTemperature_z),(1+alpha_[indexm]+gamma))* exp(-sigma*gamma*currentGeopotentialAltitude_); //std::cout << "fi: " << fi << std::endl; return fi; @@ -374,7 +400,7 @@ namespace aerodynamics // Function to compute the total density double MarsDtmAtmosphereModel::getTotalDensity( const double altitude, const double latitude, - const double longitude, const double minutes_, const double hours_, const int day_ ,const int month_,const int year_) + const double longitude, const int minutes_, const int hours_, const int day_ ,const int month_,const int year_) { double avogadroConstant_ = 6.022E23;//physical_constants::AVOGADRO_CONSTANT; double rho0; @@ -391,7 +417,8 @@ namespace aerodynamics rho += rho0*fi*exp(gl); indexm +=1; - //std::cout << "rho: " << indexm << " " << rho << std::endl; + // std::cout << "fi: " << indexm << " " << fi << std::endl; + // std::cout << "rho: " << indexm << " " << rho << std::endl; } return rho*1000; } diff --git a/src/simulation/propagation_setup/createEnvironmentUpdater.cpp b/src/simulation/propagation_setup/createEnvironmentUpdater.cpp index 60c68ac385..31e95380ef 100644 --- a/src/simulation/propagation_setup/createEnvironmentUpdater.cpp +++ b/src/simulation/propagation_setup/createEnvironmentUpdater.cpp @@ -1098,7 +1098,10 @@ std::vector< std::string > > createEnvironmentUpdaterSettingsForDependentVariabl case body_inertia_tensor: variablesToUpdate[ body_mass_distribution_update ].push_back( dependentVariableSaveSettings->associatedBody_ ); break; - default: + case solar_longitude: + variablesToUpdate[ body_mass_distribution_update ].push_back( dependentVariableSaveSettings->associatedBody_ ); + break; + default: throw std::runtime_error( "Error when getting environment updates for dependent variables, parameter " + std::to_string( dependentVariableSaveSettings->dependentVariableType_ ) + " not found." ); } diff --git a/src/simulation/propagation_setup/propagationOutput.cpp b/src/simulation/propagation_setup/propagationOutput.cpp index 0013751a74..fc72341937 100644 --- a/src/simulation/propagation_setup/propagationOutput.cpp +++ b/src/simulation/propagation_setup/propagationOutput.cpp @@ -486,6 +486,9 @@ int getDependentVariableSize( case visible_source_area: variableSize = 1; break; + case solar_longitude: + variableSize = 1; + break; default: std::string errorMessage = "Error, did not recognize dependent variable size of type: " + std::to_string( dependentVariableSettings->dependentVariableType_ ); diff --git a/src/simulation/propagation_setup/propagationOutputSettings.cpp b/src/simulation/propagation_setup/propagationOutputSettings.cpp index 752f188743..0456af6e19 100644 --- a/src/simulation/propagation_setup/propagationOutputSettings.cpp +++ b/src/simulation/propagation_setup/propagationOutputSettings.cpp @@ -307,6 +307,9 @@ std::string getDependentVariableName( case visible_source_area: variableName = "Visible area"; break; + case solar_longitude: + variableName = "Solar longitude"; + break; default: std::string errorMessage = "Error, dependent variable " + std::to_string( propagationDependentVariables ) + diff --git a/tests_riva/CMakeLists.txt b/tests_riva/CMakeLists.txt index ec23312573..191c0cbb53 100644 --- a/tests_riva/CMakeLists.txt +++ b/tests_riva/CMakeLists.txt @@ -1,4 +1,8 @@ TUDAT_ADD_EXECUTABLE(application_test_mars_dtm + "TestMarsDtm_.cpp" + ${Tudat_ESTIMATION_LIBRARIES} +) +TUDAT_ADD_EXECUTABLE(application_test_mars_dtm2 "TestMarsDtm.cpp" ${Tudat_ESTIMATION_LIBRARIES} ) diff --git a/tests_riva/TestMarsDtm.cpp b/tests_riva/TestMarsDtm.cpp index 2743c0a81c..38f3bfe04c 100644 --- a/tests_riva/TestMarsDtm.cpp +++ b/tests_riva/TestMarsDtm.cpp @@ -31,74 +31,122 @@ #include "tudat/simulation/environment_setup/body.h" #include "tudat/simulation/estimation_setup/createNumericalSimulator.h" #include "tudat/simulation/environment_setup/defaultBodies.h" - -using namespace tudat::aerodynamics; - - +#include "tudat/simulation/environment_setup/body.h" #include "tudat/astro/basic_astro/timeConversions.h" - #include "tudat/astro/aerodynamics/marsDtmAtmosphereModel.h" +#include "tudat/io/solarActivityData.h" int main( ) { + using namespace tudat; + using namespace aerodynamics; + using namespace simulation_setup; + using namespace numerical_integrators; + using namespace simulation_setup; + using namespace basic_astrodynamics; + using namespace propagators; + using namespace basic_mathematics; + using namespace basic_astrodynamics; + + //std::string filename = "/Users/ralkahal/Documents/PhD/atmodensitypds/dtm_mars"; std::string filename = "/Users/ralkahal/OneDrive - Delft University of Technology/PhD/Programs/atmodensitydtm/dtm_mars"; - MarsDtmAtmosphereModel atmosphereModel = MarsDtmAtmosphereModel(3376.78E3, filename); - //atmosphereModel.computelocalsolartime(0.0,16 ,12,2000, 0.0, 1.0697333); - //std::cout<( mars_dtm_atmosphere ); + std::shared_ptr< AtmosphereSettings > marsDtmAtmosphereSettings; + //std::shared_ptr< AtmosphereModel > atmosphereModel = + // std::dynamic_pointer_cast< MarsDtmAtmosphereModel >( + // createAtmosphereModel( marsDtmAtmosphereSettings, "Mars" ) ); + std::string spaceWeatherFilePath = paths::getSpaceWeatherDataPath( ) + "/sw19571001.txt"; + tudat::input_output::solar_activity::SolarActivityDataMap solarActivityData = + tudat::input_output::solar_activity::readSolarActivityData( spaceWeatherFilePath ); + std::shared_ptr< input_output::solar_activity::SolarActivityContainer > solarActivityContainer = + std::make_shared< input_output::solar_activity::SolarActivityContainer >( solarActivityData ); + std::function< double( const double ) > f107Function = [=](const double time) + { + return solarActivityContainer->getSolarActivityData( time )->solarRadioFlux107Observed / 2.25; + }; + + marsDtmAtmosphereSettings = std::make_shared< MarsDtmAtmosphereSettings >( + filename, 3378.0E3); + + std::shared_ptr< MarsDtmAtmosphereModel > atmosphereModel = + std::dynamic_pointer_cast< MarsDtmAtmosphereModel >( + createAtmosphereModel( marsDtmAtmosphereSettings, "Mars" ) ); + + /*if( atmosphereModel == nullptr ) + { + std::cerr<< "Atmosphere model is null" << std::endl; + } + */ + //MarsDtmAtmosphereModel(3376.78E3, filename); + //atmosphereModel->computelocalsolartime(0.0,16 ,12,2000, 0.0, 1.0697333); + //std::cout<computeGl(0.0,0.0,1.0697333,0.0,16,12,2000,1)<computeGeopotentialAltitude( 255.0E3 )<computeCurrentTemperature( 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000, 1)<computeGamma( 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000, 1)<heightDistributionFunction(255.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000, 1)<(altitude); - double rho = atmosphereModel.getTotalDensity(alt_km, 0.0, 0.0, 1.0697333, 0.0, 16, 12, 2000); + double rho = atmosphereModel->getTotalDensity(alt_km, 0.0, 0.0, 1.0697333, 0.0, 16, 12, 2000); // Write altitude and corresponding density to the file outputFile << alt_km << " " << rho << "\n"; } */ //for (int altitude = 138E3; altitude <= 1000E3; altitude += 10E3) { - // double alt_km = static_cast(altitude); - int alt_km = 138E3; - for (int time = 0; time <= 5 * 24 * 60 * 60; time += 60) { - - int minutes = time % (60 * 60); - int hours = (time / (60 * 60)) % 24; - int days = ((time / (60 * 60 * 24)) % 31) + 1; // Adding 1 to ensure days are between 1 and 31 - int months = ((time / (60 * 60 * 24 * 31)) % 12) + 1; // Adding 1 to ensure months are between 1 and 12 - int years = (time / (60 * 60 * 24 * 31 * 12)) + 2000; - - double rho = atmosphereModel.getTotalDensity(alt_km, 0.0, 0.0, minutes, hours, days, months, years); - double Ls = atmosphereModel.getSolarLongitude(); - double currentF107 = atmosphereModel.getSolarFluxIndex(); - // Write altitude and corresponding density to the file - std::cout << alt_km << "altitude" << rho << "density" << Ls << "Ls" << currentF107 << "F107" << std::endl; - outputFile << alt_km << " " << rho << "\n"; - - //} - }// Close the file when you are done - outputFile.close(); - std::cout << "Density computation and output written to file successfully." << std::endl; - //std::cout << atmosphereModel.getTotalDensity( 1000.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000) << std::endl; + // double alt_km = static_cast(altitude); + int alt_km = 400E3; + //for (int time = 0; time <= 1 * 24 * 60 * 60; time += 60) { +/* + int minutes = (time / 60) % 60; + int hours = (time / (60 * 60)) % 24; + int days = ((time / (60 * 60 * 24)) % 31) + 1; // Adding 1 to ensure days are between 1 and 31 + int months = ((time / (60 * 60 * 24 * 31)) % 12) + 1; // Adding 1 to ensure months are between 1 and 12 + int years = (time / (60 * 60 * 24 * 31 * 12)) + 2000; +*/ + // std::cout << "Time: " << time << " minutes: " << minutes << " hours: " << hours << " days: " << days << " months: " << months << " years: " << years << std::endl; + + int minutes = 0; + int hours =0; + int days = 1; + int months = 1; + int years = 2001; + + double rho = atmosphereModel->getTotalDensity(alt_km, 0.0, 0.0, minutes, hours, days, months, years); + auto results = computeSolarLongitude( 0.0, days, months, years, hours, minutes); + double Ls = std::get<0>(results); + + double currentF107 = atmosphereModel->getSolarFluxIndex(); + // Write altitude and corresponding density to the file + std::cout << rho << " " << Ls << " " << currentF107 << std::endl; + //outputFile << rho << " " << Ls << "\n"; + + // } + // Close the file when you are done + //outputFile.close(); + //std::cout << "Density computation and output written to file successfully." << std::endl; + //std::cout << atmosphereModel->getTotalDensity( 1000.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000) << std::endl; // return 0; - //std::cout << atmosphereModel.getTotalDensity( 138.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000) <getTotalDensity( 138.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000) < +#include "fstream" +#include "iostream" + +#include + +#include "tudat/basics/testMacros.h" +#include "tudat/astro/aerodynamics/customAerodynamicCoefficientInterface.h" +#include "tudat/astro/aerodynamics/aerodynamicAcceleration.h" +#include "tudat/astro/reference_frames/aerodynamicAngleCalculator.h" +#include "tudat/simulation/propagation_setup/dynamicsSimulator.h" +#include "tudat/interface/spice/spiceEphemeris.h" +#include "tudat/interface/spice/spiceRotationalEphemeris.h" +#include "tudat/io/basicInputOutput.h" +#include "tudat/simulation/environment_setup/body.h" +#include "tudat/simulation/estimation_setup/createNumericalSimulator.h" +#include "tudat/simulation/environment_setup/defaultBodies.h" +#include "tudat/astro/basic_astro/celestialBodyConstants.h" + +using namespace tudat::aerodynamics; + +int main( ) +{ + using namespace tudat; + using namespace aerodynamics; + using namespace simulation_setup; + using namespace numerical_integrators; + using namespace simulation_setup; + using namespace basic_astrodynamics; + using namespace propagators; + using namespace basic_mathematics; + using namespace basic_astrodynamics; + using namespace orbital_element_conversions; + using namespace unit_conversions; + + //Load spice kernels. + spice_interface::loadStandardSpiceKernels( ); + + // Create Earth object + BodyListSettings defaultBodySettings = + getDefaultBodySettings( { "Mars" } ); + std::string filename ="/Users/ralkahal/OneDrive - Delft University of Technology/PhD/Programs/atmodensitydtm/dtm_mars";; + defaultBodySettings.at( "Mars" )->atmosphereSettings = marsDtmAtmosphereSettings( filename, 3378.0E3); + SystemOfBodies bodies = createSystemOfBodies( defaultBodySettings ); + + // Create vehicle object. + double vehicleMass = 400; + bodies.createEmptyBody( "Vehicle" ); + bodies.at( "Vehicle" )->setConstantBodyMass( vehicleMass ); + + // Set aerodynamic coefficients. + std::shared_ptr aerodynamicCoefficientSettings = + std::make_shared( + 2.0, 4.0, Eigen::Vector3d::Zero( ), Eigen::Vector3d::UnitX( ), Eigen::Vector3d::Zero( ), + negative_aerodynamic_frame_coefficients, negative_aerodynamic_frame_coefficients ); + bodies.at( "Vehicle" )->setAerodynamicCoefficientInterface( + createAerodynamicCoefficientInterface( aerodynamicCoefficientSettings, "Vehicle", bodies )); + + // Define acceleration model settings. + SelectedAccelerationMap accelerationMap; + std::vector bodiesToPropagate; + std::vector centralBodies; + std::map > > accelerationsOfVehicle; + accelerationsOfVehicle[ "Mars" ].push_back( std::make_shared( point_mass_gravity )); + accelerationsOfVehicle[ "Mars" ].push_back( std::make_shared( aerodynamic )); + accelerationMap[ "Vehicle" ] = accelerationsOfVehicle; + bodiesToPropagate.push_back( "Vehicle" ); + centralBodies.push_back( "Mars" ); + + // Set initial state + //Eigen::Vector6d systemInitialState; + //systemInitialState << 3378.0E3, 3000E3, 4200E3, 0.0, 0.0, 0.0; + + Eigen::Vector6d vehicleInitialStateInKeplerianElements; + vehicleInitialStateInKeplerianElements( semiMajorAxisIndex ) = 3389.0E3+400.0E3; + vehicleInitialStateInKeplerianElements( eccentricityIndex ) = 0; + vehicleInitialStateInKeplerianElements( inclinationIndex ) = convertDegreesToRadians( 85.3 ); + vehicleInitialStateInKeplerianElements( argumentOfPeriapsisIndex ) = + convertDegreesToRadians( 235.7 ); + vehicleInitialStateInKeplerianElements( longitudeOfAscendingNodeIndex ) = + convertDegreesToRadians( 23.4 ); + vehicleInitialStateInKeplerianElements( trueAnomalyIndex ) = convertDegreesToRadians( 139.87 ); + + double marsGravitationalParameter = celestial_body_constants::MARS_GRAVITATIONAL_PARAMETER; + + Eigen::VectorXd systemInitialState = convertKeplerianToCartesianElements( + vehicleInitialStateInKeplerianElements, + marsGravitationalParameter ); + + // Create acceleration models and propagation settings. + basic_astrodynamics::AccelerationMap accelerationModelMap = createAccelerationModelsMap( + bodies, accelerationMap, bodiesToPropagate, centralBodies ); + + // Set variables to save + std::vector > dependentVariables; + dependentVariables.push_back( + std::make_shared( + altitude_dependent_variable, "Vehicle", "Mars" )); + dependentVariables.push_back( + std::make_shared( + "Vehicle", reference_frames::longitude_angle )); + dependentVariables.push_back( + std::make_shared( + "Vehicle", reference_frames::latitude_angle )); + dependentVariables.push_back( + std::make_shared( + local_density_dependent_variable, "Vehicle", "Mars" )); + dependentVariables.push_back( + std::make_shared( + solar_longitude,"Mars")); + dependentVariables.push_back( + std::make_shared( + aerodynamic, "Vehicle", "Mars", 0 )); + + + // Set propagation/integration settings + std::shared_ptr terminationSettings = + std::make_shared( 365*24*3600.0 ); + std::shared_ptr > integratorSettings = + std::make_shared > + ( rungeKutta4, 0.0, 60.0 ); + std::shared_ptr > + translationalPropagatorSettings = + std::make_shared > + ( centralBodies, accelerationModelMap, bodiesToPropagate, systemInitialState, 0.0, + integratorSettings, terminationSettings, + cowell, dependentVariables ); + + std::cout << "Propagation started" << std::endl; + // Create simulation object and propagate dynamics. + SingleArcDynamicsSimulator<> dynamicsSimulator( + bodies, translationalPropagatorSettings ); + + std::map< double, Eigen::VectorXd > dependentVariableHistory = dynamicsSimulator.getDependentVariableHistory( ); + //std::map< double, Eigen::VectorXd > stateHistory = dynamicsSimulator.getEquationsOfMotionNumericalSolution( ); + std::map< double, Eigen::Matrix< double, Eigen::Dynamic, 1 > > dependentVariableOutput = + dynamicsSimulator.getDependentVariableHistory( ); + + std::vector times; + std::vector densityValues; + std::vector solarLongValues; + std::cout << "Saving densities started" << std::endl; + // Iterate over the data + std::for_each(dependentVariableOutput.begin(), dependentVariableOutput.end(), + [&](const auto& it) { + // Extract the density value + double density = it.second(3); // Assuming density is at index 4 + double solarLong = it.second(4); // Assuming density is at index 4 + + // Store time and density value + //convert time to days + times.push_back(it.first/86400.0); + + densityValues.push_back(density); + solarLongValues.push_back(solarLong); + + }); + std::ofstream outputFile( + "/Users/ralkahal/OneDrive - Delft University of Technology/PhD/Programs/atmodensitydtm/density_output_test.txt"); + if (!outputFile.is_open()) { + std::cerr << "Unable to open the file!" << std::endl; + return 1; // return an error code + } + for (int i = 0; i < times.size(); ++i) { + outputFile << times[i] << " " << densityValues[i] << " " << solarLongValues[i] << std::endl; + } + outputFile.close(); + std::cout << "Density computation and output written to file successfully." << std::endl; + + +} +// +//namespace tudat +//{ +//namespace unit_tests +//{ +// +//BOOST_AUTO_TEST_SUITE( test_mars_dtm_atmosphere ) +// +//BOOST_AUTO_TEST_CASE( testMarsDtmAtmosphere ) +//{ +// +//} +// +//BOOST_AUTO_TEST_SUITE_END( ) +// +//} // namespace unit_tests +//} // namespace tudat From b19b245a0e7f208af22f3b24cb658e694b2832e1 Mon Sep 17 00:00:00 2001 From: RivaAlkahal Date: Wed, 14 Feb 2024 16:07:15 +0100 Subject: [PATCH 15/25] another test done --- .../aerodynamics/marsDtmAtmosphereModel.cpp | 2 +- tests_riva/CMakeLists.txt | 4 +- tests_riva/TestMarsDtm.cpp | 84 +++++++++++-------- 3 files changed, 52 insertions(+), 38 deletions(-) diff --git a/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp b/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp index b58984ac2e..bb4a7a829c 100644 --- a/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp +++ b/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp @@ -297,7 +297,7 @@ namespace aerodynamics double sin2h = 2.0*sin(hl0)*cos(hl0); //flux terms: double ff0 = 0.0; - double F = currentF107_ - 65.0; + double F = 65.0-65.0;//currentF107_ - 65.0; double F2 = F*F; //std::cout << "F: " << F << std::endl; double f0 = coefficients_[4][indexg]*F + coefficients_[39][indexg] * F2; diff --git a/tests_riva/CMakeLists.txt b/tests_riva/CMakeLists.txt index 191c0cbb53..f85989c5c5 100644 --- a/tests_riva/CMakeLists.txt +++ b/tests_riva/CMakeLists.txt @@ -1,8 +1,8 @@ TUDAT_ADD_EXECUTABLE(application_test_mars_dtm - "TestMarsDtm_.cpp" + "TestMarsDtm.cpp" ${Tudat_ESTIMATION_LIBRARIES} ) TUDAT_ADD_EXECUTABLE(application_test_mars_dtm2 - "TestMarsDtm.cpp" + "TestMarsDtm_.cpp" ${Tudat_ESTIMATION_LIBRARIES} ) diff --git a/tests_riva/TestMarsDtm.cpp b/tests_riva/TestMarsDtm.cpp index 38f3bfe04c..34620a4872 100644 --- a/tests_riva/TestMarsDtm.cpp +++ b/tests_riva/TestMarsDtm.cpp @@ -35,6 +35,8 @@ #include "tudat/astro/basic_astro/timeConversions.h" #include "tudat/astro/aerodynamics/marsDtmAtmosphereModel.h" #include "tudat/io/solarActivityData.h" +#include "tudat/astro/basic_astro/unitConversions.h" + int main( ) { using namespace tudat; @@ -52,7 +54,7 @@ int main( ) { std::string filename = "/Users/ralkahal/OneDrive - Delft University of Technology/PhD/Programs/atmodensitydtm/dtm_mars"; //marsDtmAtmosphereSettings = - std::make_shared< AtmosphereSettings >( mars_dtm_atmosphere ); + //std::shared_ptr< AtmosphereSettings > atmosphereSettings; std::shared_ptr< AtmosphereSettings > marsDtmAtmosphereSettings; //std::shared_ptr< AtmosphereModel > atmosphereModel = // std::dynamic_pointer_cast< MarsDtmAtmosphereModel >( @@ -92,7 +94,7 @@ int main( ) { //std::cout<computeCurrentTemperature( 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000, 1)<computeGamma( 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000, 1)<heightDistributionFunction(255.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000, 1)<(altitude); int alt_km = 400E3; - //for (int time = 0; time <= 1 * 24 * 60 * 60; time += 60) { -/* - int minutes = (time / 60) % 60; - int hours = (time / (60 * 60)) % 24; - int days = ((time / (60 * 60 * 24)) % 31) + 1; // Adding 1 to ensure days are between 1 and 31 - int months = ((time / (60 * 60 * 24 * 31)) % 12) + 1; // Adding 1 to ensure months are between 1 and 12 - int years = (time / (60 * 60 * 24 * 31 * 12)) + 2000; -*/ - // std::cout << "Time: " << time << " minutes: " << minutes << " hours: " << hours << " days: " << days << " months: " << months << " years: " << years << std::endl; - - int minutes = 0; - int hours =0; - int days = 1; - int months = 1; - int years = 2001; - - double rho = atmosphereModel->getTotalDensity(alt_km, 0.0, 0.0, minutes, hours, days, months, years); - auto results = computeSolarLongitude( 0.0, days, months, years, hours, minutes); - double Ls = std::get<0>(results); - - double currentF107 = atmosphereModel->getSolarFluxIndex(); - // Write altitude and corresponding density to the file - std::cout << rho << " " << Ls << " " << currentF107 << std::endl; - //outputFile << rho << " " << Ls << "\n"; - - // } - // Close the file when you are done - //outputFile.close(); - //std::cout << "Density computation and output written to file successfully." << std::endl; - //std::cout << atmosphereModel->getTotalDensity( 1000.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000) << std::endl; - -// return 0; + //668 * 24 * 60 * 60 + for (int time = 0.0; time <= 345600 ; time += 86400) { + //for (int time = 0; time <= 365 * 24 * 60 * 60; time += 60) { + + basic_astrodynamics::DateTime currentDateTime_ = basic_astrodynamics::getCalendarDateFromTime(time); + int minutes = currentDateTime_.getMinute(); + int hours = currentDateTime_.getHour(); + int days = currentDateTime_.getDay(); + int months = currentDateTime_.getMonth(); + int years = currentDateTime_.getYear(); + + std::string filename = + "/Users/ralkahal/OneDrive - Delft University of Technology/PhD/atmodensitypds/densityFiles/" + + std::to_string(time) + ".txt"; + std::ofstream outfile(filename); + + //std::cout << "Time: " << time << " minutes: " << minutes << " hours: " << hours << " days: " << days << " months: " << months << " years: " << years << std::endl; + //double latitude = 0.0; + //double longitude = 0.0; + if (outfile.is_open()) { + for (double latitude = -90.0; latitude <= 90.0; latitude += 1.0) { + for (double longitude = 0.0; longitude <= 360.0; longitude += 1.0) { + double rho = atmosphereModel->getTotalDensity(alt_km, unit_conversions::convertDegreesToRadians(latitude), unit_conversions::convertDegreesToRadians(longitude), minutes, hours, days, + months, years); + auto results = computeSolarLongitude(longitude, days, months, years, hours, minutes); + double Ls = std::get<0>(results); + double currentF107 = 0.0; // Assuming this is defined somewhere + outfile << latitude << " " << longitude << " " << rho << " " << Ls << "\n"; + } + } + + } else { + std::cerr << "Error creating file " << filename << std::endl; + } + outfile.close(); + + + + + //std::cout << "Density computation and output written to file successfully." << std::endl; + //std::cout << atmosphereModel->getTotalDensity( 1000.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000) << std::endl; + } + return 0; //std::cout << atmosphereModel->getTotalDensity( 138.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000) < Date: Thu, 4 Jul 2024 12:55:25 +0200 Subject: [PATCH 16/25] Correcting exception in gravity field file reading --- .../environment_setup/createGravityField.cpp | 64 +++++++++++-------- 1 file changed, 37 insertions(+), 27 deletions(-) diff --git a/src/simulation/environment_setup/createGravityField.cpp b/src/simulation/environment_setup/createGravityField.cpp index be0398e94c..e83ba7d512 100644 --- a/src/simulation/environment_setup/createGravityField.cpp +++ b/src/simulation/environment_setup/createGravityField.cpp @@ -220,47 +220,57 @@ std::pair< double, double > readGravityFieldFile( Eigen::MatrixXd sineCoefficients = Eigen::MatrixXd( maximumDegree + 1, maximumOrder + 1 ); sineCoefficients.setZero( ); + int counter = 0; + std::string previousLine = ""; // Read coefficients up to required maximum degree and order. while ( !stream.fail( ) && !stream.eof( ) && ( currentDegree <= maximumDegree && currentOrder < maximumOrder ) ) { // Read current line std::getline( stream, line ); - - // Trim input string (removes all leading and trailing whitespaces). - boost::algorithm::trim( line ); - - // Split string into multiple strings, each containing one element from a line from the - // data file. - boost::algorithm::split( vectorOfIndividualStrings, - line, - boost::algorithm::is_any_of( ", \t" ), - boost::algorithm::token_compress_on ); - - // Check current line for consistency - if( vectorOfIndividualStrings.size( ) != 0 ) + if( !line.empty( ) ) { - if( vectorOfIndividualStrings.size( ) < 4 ) + // Trim input string (removes all leading and trailing whitespaces). + boost::algorithm::trim( line ); + + // Split string into multiple strings, each containing one element from a line from the + // data file. + boost::algorithm::split( vectorOfIndividualStrings, + line, + boost::algorithm::is_any_of( ", \t" ), + boost::algorithm::token_compress_on ); + + // Check current line for consistency + if ( vectorOfIndividualStrings.size( ) != 0 ) { - std::string errorMessage = "Error when reading pds gravity field file, number of fields is " + - std::to_string( vectorOfIndividualStrings.size( ) ); - throw std::runtime_error( errorMessage ); - } - else - { - // Read current degree and orde from line. - currentDegree = static_cast< int >( std::round( std::stod( vectorOfIndividualStrings[ 0 ] ) ) ); - currentOrder = static_cast< int >( std::round( std::stod( vectorOfIndividualStrings[ 1 ] ) ) ); - // Set cosine and sine coefficients for current degree and order. - if( currentDegree <= maximumDegree && currentOrder <= maximumOrder ) + if ( vectorOfIndividualStrings.size( ) < 4 && !( + vectorOfIndividualStrings.size( ) == 1 && vectorOfIndividualStrings.at( 0 ).size( ) == 0 )) { - cosineCoefficients( currentDegree, currentOrder ) = + std::string errorMessage = fileName + "; " + line + "; " + previousLine + "; " + + std::to_string( counter ) + + "Error when reading pds gravity field file, number of fields is " + + std::to_string( vectorOfIndividualStrings.size( )) + "; full line is " + + line + "; previouse line was " + previousLine; + throw std::runtime_error( errorMessage ); + } + else + { + // Read current degree and orde from line. + currentDegree = static_cast< int >( std::round( std::stod( vectorOfIndividualStrings[ 0 ] ))); + currentOrder = static_cast< int >( std::round( std::stod( vectorOfIndividualStrings[ 1 ] ))); + // Set cosine and sine coefficients for current degree and order. + if ( currentDegree <= maximumDegree && currentOrder <= maximumOrder ) + { + cosineCoefficients( currentDegree, currentOrder ) = std::stod( vectorOfIndividualStrings[ 2 ] ); - sineCoefficients( currentDegree, currentOrder ) = + sineCoefficients( currentDegree, currentOrder ) = std::stod( vectorOfIndividualStrings[ 3 ] ); + } } } } + previousLine = line; + counter++; } // Set cosine coefficient at (0,0) to 1. From a91e8c6234c25f0d29cba47cfd90969ab75337d0 Mon Sep 17 00:00:00 2001 From: Dominic Dirkx Date: Fri, 5 Jul 2024 12:15:22 +0200 Subject: [PATCH 17/25] Polynomial gravity field variation parameter mostly working --- .../polynomialGravityFieldVariations.h | 34 ++- .../estimatableParameter.h | 3 +- .../gravityFieldVariationParameters.h | 195 ++++++++++++++++++ .../createEstimatableParameters.h | 48 +++++ .../polynomialGravityFieldVariations.cpp | 21 ++ 5 files changed, 298 insertions(+), 3 deletions(-) create mode 100644 include/tudat/astro/orbit_determination/estimatable_parameters/gravityFieldVariationParameters.h diff --git a/include/tudat/astro/gravitation/polynomialGravityFieldVariations.h b/include/tudat/astro/gravitation/polynomialGravityFieldVariations.h index 5dc749e108..b0533b2fc0 100644 --- a/include/tudat/astro/gravitation/polynomialGravityFieldVariations.h +++ b/include/tudat/astro/gravitation/polynomialGravityFieldVariations.h @@ -36,12 +36,42 @@ class PolynomialGravityFieldVariations: public GravityFieldVariations std::pair< Eigen::MatrixXd, Eigen::MatrixXd > calculateSphericalHarmonicsCorrections( const double time ); + std::map< int, Eigen::MatrixXd > getCosineAmplitudes( ) + { + return cosineAmplitudes_; + } + + std::map< int, Eigen::MatrixXd > getSineAmplitudes( ) + { + return sineAmplitudes_; + } + + void resetCosineAmplitudes( const std::map< int, Eigen::MatrixXd > cosineAmplitudes ) + { + cosineAmplitudes_ = cosineAmplitudes; + } + + void resetSineAmplitudes( const std::map< int, Eigen::MatrixXd > sineAmplitudes ) + { + sineAmplitudes_ = sineAmplitudes; + } + + + std::map< int, Eigen::MatrixXd >& getCosineAmplitudesReference( ) + { + return cosineAmplitudes_; + } + + std::map< int, Eigen::MatrixXd >& getSineAmplitudesReference( ) + { + return sineAmplitudes_; + } protected: - const std::map< int, Eigen::MatrixXd > cosineAmplitudes_; + std::map< int, Eigen::MatrixXd > cosineAmplitudes_; - const std::map< int, Eigen::MatrixXd > sineAmplitudes_; + std::map< int, Eigen::MatrixXd > sineAmplitudes_; const double referenceEpoch_; diff --git a/include/tudat/astro/orbit_determination/estimatable_parameters/estimatableParameter.h b/include/tudat/astro/orbit_determination/estimatable_parameters/estimatableParameter.h index 659bc92800..cc99083853 100644 --- a/include/tudat/astro/orbit_determination/estimatable_parameters/estimatableParameter.h +++ b/include/tudat/astro/orbit_determination/estimatable_parameters/estimatableParameter.h @@ -77,7 +77,8 @@ enum EstimatebleParametersEnum inverse_tidal_quality_factor, yarkovsky_parameter, custom_estimated_parameter, - reference_point_position + reference_point_position, + poynomial_gravity_field_variation_amplitudes }; std::string getParameterTypeString( const EstimatebleParametersEnum parameterType ); diff --git a/include/tudat/astro/orbit_determination/estimatable_parameters/gravityFieldVariationParameters.h b/include/tudat/astro/orbit_determination/estimatable_parameters/gravityFieldVariationParameters.h new file mode 100644 index 0000000000..6f82398011 --- /dev/null +++ b/include/tudat/astro/orbit_determination/estimatable_parameters/gravityFieldVariationParameters.h @@ -0,0 +1,195 @@ +/* Copyright (c) 2010-2019, 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 + * a copy of the license with this file. If not, please or visit: + * http://tudat.tudelft.nl/LICENSE. + */ + +#ifndef TUDAT_GRAVITYFIELDVARIATIONPARAMETERS_H +#define TUDAT_GRAVITYFIELDVARIATIONPARAMETERS_H + +#include +#include +#include +#include +#include + + +#include +#include +#include + +#include "tudat/astro/gravitation/polynomialGravityFieldVariations.h" +#include "tudat/astro/orbit_determination/estimatable_parameters/estimatableParameter.h" + +namespace tudat +{ + +namespace estimatable_parameters +{ + +class PolynomialGravityFieldVariationsParameters: public EstimatableParameter< Eigen::VectorXd > +{ + +public: + //! Constructor. + /*! + * Constructor taking parameter name and associated body. All parameters are identified by a these two variables. + * Any additional information that may be required for uniquely defining a parameter is to be defined in the derived class. + * \param parameterName Enum value defining the type of the parameter. + * \param associatedBody Name of body associated with patameters + * \param pointOnBodyId Reference point on body associated with parameter (empty by default). + */ + PolynomialGravityFieldVariationsParameters( + const std::shared_ptr< gravitation::PolynomialGravityFieldVariations > polynomialVariationModel, + const std::map< int, std::vector< std::pair< int, int > > >& cosineBlockIndices, + const std::map< int, std::vector< std::pair< int, int > > >& sineBlockIndices, + const std::string& bodyName ): + EstimatableParameter< Eigen::VectorXd >( poynomial_gravity_field_variation_amplitudes, bodyName ), + polynomialVariationModel_( polynomialVariationModel ) + { + std::map< int, Eigen::MatrixXd > cosineVariations = polynomialVariationModel->getCosineAmplitudes( ); + for( auto it : cosineBlockIndices ) + { + if( cosineVariations.count( it.first ) == 0 ) + { + throw std::runtime_error( "Error when estimationg gravity field polynomial corrections of body " + bodyName + + ", not polynomial term of order " + std::to_string( it.first ) + " found for cosine coefficients" ); + } + + for( unsigned int i = 0; i < it.second.size( ); i++ ) + { + if( it.second.at( i ).first > polynomialVariationModel_->getMaximumDegree( ) || + it.second.at( i ).first < polynomialVariationModel_->getMinimumDegree( ) ) + { + throw std::runtime_error( "Error when estimating gravity field polynomial corrections of body " + bodyName + + " of polynomial term of order " + std::to_string( it.first ) + + ", no cosine coefficient variation of degree " + std::to_string( it.second.at( i ).first ) + " found." ); + } + + if( it.second.at( i ).second > polynomialVariationModel_->getMaximumOrder( ) || + it.second.at( i ).second < polynomialVariationModel_->getMaximumOrder( ) ) + { + throw std::runtime_error( "Error when estimating gravity field polynomial corrections of body " + bodyName + + " of polynomial term of order " + std::to_string( it.first ) + + ", no cosine coefficient variation of order " + std::to_string( it.second.at( i ).second ) + " found." ); + } + cosineCorrectionIndices_.push_back( std::make_tuple( it.first, it.second.at( i ).first, it.second.at( i ).second ) ); + } + } + + std::map< int, Eigen::MatrixXd > sineVariations = polynomialVariationModel->getSineAmplitudes( ); + for( auto it : sineBlockIndices ) + { + if( sineVariations.count( it.first ) == 0 ) + { + throw std::runtime_error( "Error when estimationg gravity field polynomial corrections of body " + bodyName + + ", not polynomial term of order " + std::to_string( it.first ) + " found for sine coefficients" ); + } + + for( unsigned int i = 0; i < it.second.size( ); i++ ) + { + if( it.second.at( i ).first > polynomialVariationModel_->getMaximumDegree( ) || + it.second.at( i ).first < polynomialVariationModel_->getMinimumDegree( ) ) + { + throw std::runtime_error( "Error when estimating gravity field polynomial corrections of body " + bodyName + + " of polynomial term of order " + std::to_string( it.first ) + + ", no sine coefficient variation of degree " + std::to_string( it.second.at( i ).first ) + " found." ); + } + + if( it.second.at( i ).second > polynomialVariationModel_->getMaximumOrder( ) || + it.second.at( i ).second < polynomialVariationModel_->getMaximumOrder( ) ) + { + throw std::runtime_error( "Error when estimating gravity field polynomial corrections of body " + bodyName + + " of polynomial term of order " + std::to_string( it.first ) + + ", no sine coefficient variation of order " + std::to_string( it.second.at( i ).second ) + " found." ); + } + } + sineCorrectionIndices_.push_back( std::make_tuple( it.first, it.second.at( i ).first, it.second.at( i ).second ) ); + } + } + + //! Virtual destructor. + ~PolynomialGravityFieldVariationsParameters( ) { } + + //! Pure virtual function to retrieve the value of the parameter + /*! + * Pure virtual function to retrieve the value of the parameter + * \return Current value of parameter. + */ + Eigen::VectorXd getParameterValue( ) + { + Eigen::VectorXd polynomialCorrections = Eigen::VectorXd::Zero( getParameterSize( ) ); + for( unsigned int i = 0; i < cosineCorrectionIndices_.size( ); i++ ) + { + polynomialCorrections( i ) = polynomialVariationModel_->getCosineAmplitudesReference( ).at( std::get< 0 >( cosineCorrectionIndices_.at( i ) ) ) + ( std::get< 1 >( cosineCorrectionIndices_.at( i ) ), std::get< 2 >( cosineCorrectionIndices_.at( i ) ) ); + } + + for( unsigned int i = cosineCorrectionIndices_.size( ); i < polynomialCorrections.rows( ); i++ ) + { + polynomialCorrections( i ) = polynomialVariationModel_->getSineAmplitudesReference( ).at( std::get< 0 >( sineCorrectionIndices_.at( i ) ) ) + ( std::get< 1 >( sineCorrectionIndices_.at( i ) ), std::get< 2 >( cosineCorrectionIndices_.at( i ) ) ); + } + } + + //! Pure virtual function to (re)set the value of the parameter. + /*! + * Pure virtual function to (re)set the value of the parameter. + * \param parameterValue to which the parameter is to be set. + */ + void setParameterValue( const Eigen::VectorXd parameterValue ) + { + if( parameterValue.rows( ) != getParameterSize( ) ) + { + throw std::runtime_error( "Error when resetting polynomial gravity field variation parameter; sizes are incompatible" ); + } + std::map< int, Eigen::MatrixXd > cosineVariations = polynomialVariationModel_->getCosineAmplitudes( ); + std::map< int, Eigen::MatrixXd > sineVariations = polynomialVariationModel_->getSineAmplitudes( ); + + for( unsigned int i = 0; i < cosineCorrectionIndices_.size( ); i++ ) + { + cosineVariations[ std::get< 0 >( cosineCorrectionIndices_.at( i ) ) ]( std::get< 1 >( cosineCorrectionIndices_.at( i ) ), std::get< 2 >( cosineCorrectionIndices_.at( i ) ) ) = + parameterValue( i ); + } + + for( unsigned int i = cosineCorrectionIndices_.size( ); i < parameterValue.rows( ); i++ ) + { + sineVariations[ std::get< 0 >( sineCorrectionIndices_.at( i ) ) ]( std::get< 1 >( sineCorrectionIndices_.at( i ) ), std::get< 2 >( sineCorrectionIndices_.at( i ) ) ) = + parameterValue( i ); + } + polynomialVariationModel_->resetCosineAmplitudes( cosineVariations ); + polynomialVariationModel_->resetSineAmplitudes( sineVariations ); + } + + //! Function to retrieve the size of the parameter + /*! + * Pure virtual function to retrieve the size of the parameter (i.e. 1 for double parameters) + * \return Size of parameter value. + */ + int getParameterSize( ) + { + return cosineCorrectionIndices_.size( ) + sineCorrectionIndices_.size( ); + } + + + +protected: + + std::shared_ptr< gravitation::PolynomialGravityFieldVariations > polynomialVariationModel_; + + std::vector< std::tuple< int, int, int > > cosineCorrectionIndices_; + + std::vector< std::tuple< int, int, int > > sineCorrectionIndices_; + +}; + + +} // namespace estimatable_parameters + +} // namespace tudat + +#endif // TUDAT_GRAVITYFIELDVARIATIONPARAMETERS_H diff --git a/include/tudat/simulation/estimation_setup/createEstimatableParameters.h b/include/tudat/simulation/estimation_setup/createEstimatableParameters.h index b7d77b7907..a070d05c7c 100644 --- a/include/tudat/simulation/estimation_setup/createEstimatableParameters.h +++ b/include/tudat/simulation/estimation_setup/createEstimatableParameters.h @@ -43,6 +43,7 @@ #include "tudat/astro/orbit_determination/estimatable_parameters/constantThrust.h" #include "tudat/astro/orbit_determination/estimatable_parameters/yarkovskyParameter.h" #include "tudat/astro/orbit_determination/estimatable_parameters/referencePointPosition.h" +#include "tudat/astro/orbit_determination/estimatable_parameters/gravityFieldVariationParameters.h" #include "tudat/astro/relativity/metric.h" #include "tudat/astro/basic_astro/accelerationModelTypes.h" #include "tudat/simulation/estimation_setup/estimatableParameterSettings.h" @@ -1981,6 +1982,53 @@ std::shared_ptr< estimatable_parameters::EstimatableParameter< Eigen::VectorXd > } break; } + case poynomial_gravity_field_variation_amplitudes: + { + // Check consistency of body gravity field + std::shared_ptr< GravityFieldModel > gravityField = currentBody->getGravityFieldModel( ); + std::shared_ptr< TimeDependentSphericalHarmonicsGravityField > timeDepGravityField = + std::dynamic_pointer_cast< TimeDependentSphericalHarmonicsGravityField >( gravityField ); + if( timeDepGravityField == nullptr ) + { + throw std::runtime_error( + "Error, requested polynomial gravity field variation parameter of " + + vectorParameterName->parameterType_.second.first + + ", but body does not have a time dependent spherical harmonic gravity field." ); + } + else if( currentBody->getGravityFieldVariationSet( ) == nullptr ) + { + throw std::runtime_error( "Error, requested polynomial gravity field variation parameter of " + + vectorParameterName->parameterType_.second.first + + ", but body does not have gravity field variations" ); + } + else + { + + // Get associated gravity field variation + std::pair< bool, std::shared_ptr< gravitation::GravityFieldVariations > > polynomialVariaton = + currentBody->getGravityFieldVariationSet( )->getGravityFieldVariation( polynomial_variation ); + if( polynomialVariaton.first == 0 ) + { + throw std::runtime_error( "Error when creating polynomial gravity field variation parameter; associated gravity field model not found." ); + } + std::shared_ptr< gravitation::PolynomialGravityFieldVariations > gravityFieldVariation = + std::dynamic_pointer_cast< gravitation::PolynomialGravityFieldVariations >( + polynomialVariaton.second ); + + // Create parameter object + if( gravityFieldVariation != nullptr ) + { + vectorParameterToEstimate = std::make_shared< PolynomialGravityFieldVariationsParameters >( + polynomialVariaton, cosineVariationsToEstimate, sineVariationsToEstimate, currentBodyName ); + } + else + { + throw std::runtime_error( + "Error, expected PolynomialGravityFieldVariations when creating polynomial gravity field variation parameter" ); + } + } + break; + } case custom_estimated_parameter: { std::shared_ptr< CustomEstimatableParameterSettings > customParameterSettings = diff --git a/src/astro/gravitation/polynomialGravityFieldVariations.cpp b/src/astro/gravitation/polynomialGravityFieldVariations.cpp index 1998dfe6e0..938da11a68 100644 --- a/src/astro/gravitation/polynomialGravityFieldVariations.cpp +++ b/src/astro/gravitation/polynomialGravityFieldVariations.cpp @@ -85,6 +85,27 @@ PolynomialGravityFieldVariations::PolynomialGravityFieldVariations( sineAmplitudes_( sineAmplitudes ), referenceEpoch_( referenceEpoch ) { + for( auto it : cosineAmplitudes_ ) + { + if( cosineAmplitudes_.count( it.first ) != 0 ) + { + if( ( it.second.rows( ) != sineAmplitudes_.at( it.first ).rows( ) ) || ( it.second.cols( ) != sineAmplitudes_.at( it.first ).cols( ) ) ) + { + throw std::runtime_error( "Error when creating polynomial gravity field variation, sine and cosine sizes are inconsistent" ); + } + } + } + + for( auto it : sineAmplitudes_ ) + { + if( sineAmplitudes_.count( it.first ) != 0 ) + { + if( ( it.second.rows( ) != cosineAmplitudes_.at( it.first ).rows( ) ) || ( it.second.cols( ) != cosineAmplitudes_.at( it.first ).cols( ) ) ) + { + throw std::runtime_error( "Error when creating polynomial gravity field variation, sine and cosine sizes are inconsistent" ); + } + } + } } From 82440c5dca5032c591ba0f1a22a1efd60f60ff60 Mon Sep 17 00:00:00 2001 From: Dominic Dirkx Date: Fri, 5 Jul 2024 21:21:11 +0200 Subject: [PATCH 18/25] Added estimation needs for polynomial gravity field variation amplitudes (completely untested; only compiling) --- .../polynomialGravityFieldVariations.h | 5 ++ .../sphericalHarmonicAccelerationPartial.h | 8 ++ .../estimatableParameter.h | 5 +- .../gravityFieldVariationParameters.h | 41 ++++++++-- .../createEstimatableParameters.h | 24 ++++-- .../estimatableParameterSettings.h | 25 +++++++ .../sphericalHarmonicAccelerationPartial.cpp | 75 +++++++++++++++++++ .../estimatableParameter.cpp | 16 ++++ 8 files changed, 185 insertions(+), 14 deletions(-) diff --git a/include/tudat/astro/gravitation/polynomialGravityFieldVariations.h b/include/tudat/astro/gravitation/polynomialGravityFieldVariations.h index b0533b2fc0..745e561997 100644 --- a/include/tudat/astro/gravitation/polynomialGravityFieldVariations.h +++ b/include/tudat/astro/gravitation/polynomialGravityFieldVariations.h @@ -67,6 +67,11 @@ class PolynomialGravityFieldVariations: public GravityFieldVariations return sineAmplitudes_; } + double getReferenceEpoch( ) + { + return referenceEpoch_; + } + protected: std::map< int, Eigen::MatrixXd > cosineAmplitudes_; diff --git a/include/tudat/astro/orbit_determination/acceleration_partials/sphericalHarmonicAccelerationPartial.h b/include/tudat/astro/orbit_determination/acceleration_partials/sphericalHarmonicAccelerationPartial.h index 9fc1957899..d7191a4948 100644 --- a/include/tudat/astro/orbit_determination/acceleration_partials/sphericalHarmonicAccelerationPartial.h +++ b/include/tudat/astro/orbit_determination/acceleration_partials/sphericalHarmonicAccelerationPartial.h @@ -332,6 +332,14 @@ class SphericalHarmonicsGravityPartial: public AccelerationPartial } } + void wrtPolynomialGravityFieldVariations( + const std::vector< std::pair< int, int > >& cosineBlockIndices, + const std::vector< std::pair< int, int > >& sineBlockIndices, + const std::vector< std::vector< std::pair< int, int > > > powersPerCosineBlockIndex, + const std::vector< std::vector< std::pair< int, int > > > powersPerSineBlockIndex, + const double referenceEpoch, + Eigen::MatrixXd& partialDerivatives ); + //! Function to calculate the partial of the acceleration wrt a set of cosine coefficients. /*! * Function to calculate the partial of the acceleration wrt a set of cosine coefficients. diff --git a/include/tudat/astro/orbit_determination/estimatable_parameters/estimatableParameter.h b/include/tudat/astro/orbit_determination/estimatable_parameters/estimatableParameter.h index cc99083853..827465232c 100644 --- a/include/tudat/astro/orbit_determination/estimatable_parameters/estimatableParameter.h +++ b/include/tudat/astro/orbit_determination/estimatable_parameters/estimatableParameter.h @@ -78,7 +78,7 @@ enum EstimatebleParametersEnum yarkovsky_parameter, custom_estimated_parameter, reference_point_position, - poynomial_gravity_field_variation_amplitudes + polynomial_gravity_field_variation_amplitudes }; std::string getParameterTypeString( const EstimatebleParametersEnum parameterType ); @@ -132,6 +132,9 @@ bool isParameterObservationLinkTimeProperty( const EstimatebleParametersEnum par */ bool isParameterTidalProperty( const EstimatebleParametersEnum parameterType ); +//! Function to determine whether the given parameter influences a body's non-tidal gravity field variations. +bool isParameterNonTidalGravityFieldVariationProperty( const EstimatebleParametersEnum parameterType ); + //! Function to determine whether the given parameter represents an arc-wise initial dynamical state. /*! * Function to determine whether the given parameter represents an arc-wise initial dynamical state. diff --git a/include/tudat/astro/orbit_determination/estimatable_parameters/gravityFieldVariationParameters.h b/include/tudat/astro/orbit_determination/estimatable_parameters/gravityFieldVariationParameters.h index 6f82398011..6e8cbd1b8e 100644 --- a/include/tudat/astro/orbit_determination/estimatable_parameters/gravityFieldVariationParameters.h +++ b/include/tudat/astro/orbit_determination/estimatable_parameters/gravityFieldVariationParameters.h @@ -45,14 +45,15 @@ class PolynomialGravityFieldVariationsParameters: public EstimatableParameter< E */ PolynomialGravityFieldVariationsParameters( const std::shared_ptr< gravitation::PolynomialGravityFieldVariations > polynomialVariationModel, - const std::map< int, std::vector< std::pair< int, int > > >& cosineBlockIndices, - const std::map< int, std::vector< std::pair< int, int > > >& sineBlockIndices, + const std::map< int, std::vector< std::pair< int, int > > >& cosineBlockIndicesPerPower, + const std::map< int, std::vector< std::pair< int, int > > >& sineBlockIndicesPerPower, const std::string& bodyName ): - EstimatableParameter< Eigen::VectorXd >( poynomial_gravity_field_variation_amplitudes, bodyName ), + EstimatableParameter< Eigen::VectorXd >( polynomial_gravity_field_variation_amplitudes, bodyName ), polynomialVariationModel_( polynomialVariationModel ) { std::map< int, Eigen::MatrixXd > cosineVariations = polynomialVariationModel->getCosineAmplitudes( ); - for( auto it : cosineBlockIndices ) + int cosineIndexCounter = 0; + for( auto it : cosineBlockIndicesPerPower ) { if( cosineVariations.count( it.first ) == 0 ) { @@ -78,11 +79,15 @@ class PolynomialGravityFieldVariationsParameters: public EstimatableParameter< E ", no cosine coefficient variation of order " + std::to_string( it.second.at( i ).second ) + " found." ); } cosineCorrectionIndices_.push_back( std::make_tuple( it.first, it.second.at( i ).first, it.second.at( i ).second ) ); + indexAndPowerPerCosineBlockIndex_[ std::make_pair( it.second.at( i ).first, it.second.at( i ).second ) ].push_back( + std::make_pair( cosineIndexCounter, it.first ) ); + cosineIndexCounter++; } } std::map< int, Eigen::MatrixXd > sineVariations = polynomialVariationModel->getSineAmplitudes( ); - for( auto it : sineBlockIndices ) + int sineIndexCounter = 0; + for( auto it : sineBlockIndicesPerPower ) { if( sineVariations.count( it.first ) == 0 ) { @@ -107,8 +112,13 @@ class PolynomialGravityFieldVariationsParameters: public EstimatableParameter< E " of polynomial term of order " + std::to_string( it.first ) + ", no sine coefficient variation of order " + std::to_string( it.second.at( i ).second ) + " found." ); } + + sineCorrectionIndices_.push_back( std::make_tuple( it.first, it.second.at( i ).first, it.second.at( i ).second ) ); + indexAndPowerPerSineBlockIndex_[ std::make_pair( it.second.at( i ).first, it.second.at( i ).second ) ].push_back( + std::make_pair( sineIndexCounter, it.first ) ); + sineIndexCounter++; } - sineCorrectionIndices_.push_back( std::make_tuple( it.first, it.second.at( i ).first, it.second.at( i ).second ) ); + } } @@ -134,6 +144,8 @@ class PolynomialGravityFieldVariationsParameters: public EstimatableParameter< E polynomialCorrections( i ) = polynomialVariationModel_->getSineAmplitudesReference( ).at( std::get< 0 >( sineCorrectionIndices_.at( i ) ) ) ( std::get< 1 >( sineCorrectionIndices_.at( i ) ), std::get< 2 >( cosineCorrectionIndices_.at( i ) ) ); } + + return polynomialCorrections; } //! Pure virtual function to (re)set the value of the parameter. @@ -175,6 +187,19 @@ class PolynomialGravityFieldVariationsParameters: public EstimatableParameter< E return cosineCorrectionIndices_.size( ) + sineCorrectionIndices_.size( ); } + std::map< std::pair< int, int >, std::vector< std::pair< int, int > > > getIndexAndPowerPerCosineBlockIndex( ) + { + return indexAndPowerPerCosineBlockIndex_; + } + + std::map< std::pair< int, int >, std::vector< std::pair< int, int > > > getIndexAndPowerPerSineBlockIndex( ) + { + return indexAndPowerPerSineBlockIndex_; + } + std::shared_ptr< gravitation::PolynomialGravityFieldVariations > getPolynomialVariationModel( ) + { + return polynomialVariationModel_; + } protected: @@ -185,6 +210,10 @@ class PolynomialGravityFieldVariationsParameters: public EstimatableParameter< E std::vector< std::tuple< int, int, int > > sineCorrectionIndices_; + std::map< std::pair< int, int >, std::vector< std::pair< int, int > > > indexAndPowerPerCosineBlockIndex_; + + std::map< std::pair< int, int >, std::vector< std::pair< int, int > > > indexAndPowerPerSineBlockIndex_; + }; diff --git a/include/tudat/simulation/estimation_setup/createEstimatableParameters.h b/include/tudat/simulation/estimation_setup/createEstimatableParameters.h index a070d05c7c..c8052b443c 100644 --- a/include/tudat/simulation/estimation_setup/createEstimatableParameters.h +++ b/include/tudat/simulation/estimation_setup/createEstimatableParameters.h @@ -1982,8 +1982,15 @@ std::shared_ptr< estimatable_parameters::EstimatableParameter< Eigen::VectorXd > } break; } - case poynomial_gravity_field_variation_amplitudes: + case polynomial_gravity_field_variation_amplitudes: { + std::shared_ptr< PolynomialGravityFieldVariationEstimatableParameterSettings > gravityFieldVariationSettings = + std::dynamic_pointer_cast< PolynomialGravityFieldVariationEstimatableParameterSettings >( vectorParameterName ); + if( gravityFieldVariationSettings == nullptr ) + { + throw std::runtime_error( "Error, expected polynomial gravity field variation parameter settings " ); + } + // Check consistency of body gravity field std::shared_ptr< GravityFieldModel > gravityField = currentBody->getGravityFieldModel( ); std::shared_ptr< TimeDependentSphericalHarmonicsGravityField > timeDepGravityField = @@ -2005,21 +2012,24 @@ std::shared_ptr< estimatable_parameters::EstimatableParameter< Eigen::VectorXd > { // Get associated gravity field variation - std::pair< bool, std::shared_ptr< gravitation::GravityFieldVariations > > polynomialVariaton = + std::pair< bool, std::shared_ptr< gravitation::GravityFieldVariations > > gravityFieldVariation = currentBody->getGravityFieldVariationSet( )->getGravityFieldVariation( polynomial_variation ); - if( polynomialVariaton.first == 0 ) + if( gravityFieldVariation.first == 0 ) { throw std::runtime_error( "Error when creating polynomial gravity field variation parameter; associated gravity field model not found." ); } - std::shared_ptr< gravitation::PolynomialGravityFieldVariations > gravityFieldVariation = + std::shared_ptr< gravitation::PolynomialGravityFieldVariations > polynomialVariaton = std::dynamic_pointer_cast< gravitation::PolynomialGravityFieldVariations >( - polynomialVariaton.second ); + gravityFieldVariation.second ); // Create parameter object - if( gravityFieldVariation != nullptr ) + if( polynomialVariaton != nullptr ) { vectorParameterToEstimate = std::make_shared< PolynomialGravityFieldVariationsParameters >( - polynomialVariaton, cosineVariationsToEstimate, sineVariationsToEstimate, currentBodyName ); + polynomialVariaton, + gravityFieldVariationSettings->cosineBlockIndicesPerPower_, + gravityFieldVariationSettings->sineBlockIndicesPerPower_, + currentBodyName ); } else { diff --git a/include/tudat/simulation/estimation_setup/estimatableParameterSettings.h b/include/tudat/simulation/estimation_setup/estimatableParameterSettings.h index a31052511c..e23ba0a70d 100644 --- a/include/tudat/simulation/estimation_setup/estimatableParameterSettings.h +++ b/include/tudat/simulation/estimation_setup/estimatableParameterSettings.h @@ -988,6 +988,31 @@ class InverseTidalQualityFactorEstimatableParameterSettings: public EstimatableP }; + +class PolynomialGravityFieldVariationEstimatableParameterSettings: public EstimatableParameterSettings +{ +public: + + //! Constructor + /*! + * Constructor + * \param associatedBody Body being deformed + * \param deformingBody Body causing deformed + */ + PolynomialGravityFieldVariationEstimatableParameterSettings( + const std::string &associatedBody, + const std::map > > &cosineBlockIndicesPerPower, + const std::map > > &sineBlockIndicesPerPower ) : + EstimatableParameterSettings( associatedBody, polynomial_gravity_field_variation_amplitudes ), + cosineBlockIndicesPerPower_( cosineBlockIndicesPerPower ), + sineBlockIndicesPerPower_( sineBlockIndicesPerPower ){ } + + std::map< int, std::vector< std::pair< int, int > > > cosineBlockIndicesPerPower_; + std::map< int, std::vector< std::pair< int, int > > > sineBlockIndicesPerPower_; + +}; + + class CustomEstimatableParameterSettings: public EstimatableParameterSettings { public: diff --git a/src/astro/orbit_determination/acceleration_partials/sphericalHarmonicAccelerationPartial.cpp b/src/astro/orbit_determination/acceleration_partials/sphericalHarmonicAccelerationPartial.cpp index 1548b50b94..1f53012ef3 100644 --- a/src/astro/orbit_determination/acceleration_partials/sphericalHarmonicAccelerationPartial.cpp +++ b/src/astro/orbit_determination/acceleration_partials/sphericalHarmonicAccelerationPartial.cpp @@ -14,6 +14,7 @@ #include "tudat/astro/orbit_determination/acceleration_partials/centralGravityAccelerationPartial.h" #include "tudat/astro/orbit_determination/estimatable_parameters/sphericalHarmonicCosineCoefficients.h" #include "tudat/astro/orbit_determination/estimatable_parameters/sphericalHarmonicSineCoefficients.h" +#include "tudat/astro/orbit_determination/estimatable_parameters/gravityFieldVariationParameters.h" namespace tudat { @@ -238,6 +239,31 @@ std::pair< std::function< void( Eigen::MatrixXd& ) >, int > SphericalHarmonicsGr } } } + else if( estimatable_parameters::isParameterNonTidalGravityFieldVariationProperty( parameter->getParameterName( ).first ) ) + { + switch( parameter->getParameterName( ).first ) + { + case polynomial_gravity_field_variation_amplitudes: + { + std::shared_ptr< PolynomialGravityFieldVariationsParameters > polynomialVariationParameter = + std::dynamic_pointer_cast< PolynomialGravityFieldVariationsParameters >( parameter ); + std::map< std::pair< int, int >, std::vector< std::pair< int, int > > > indexAndPowerPerCosineBlockIndex = + polynomialVariationParameter->getIndexAndPowerPerCosineBlockIndex( ); + std::map< std::pair< int, int >, std::vector< std::pair< int, int > > > indexAndPowerPerSineBlockIndex = + polynomialVariationParameter->getIndexAndPowerPerSineBlockIndex( ); + + partialFunction = std::bind( &SphericalHarmonicsGravityPartial::wrtPolynomialGravityFieldVariations, this, + utilities::createVectorFromMapKeys( indexAndPowerPerCosineBlockIndex ), + utilities::createVectorFromMapKeys( indexAndPowerPerSineBlockIndex ), + utilities::createVectorFromMapValues( indexAndPowerPerCosineBlockIndex ), + utilities::createVectorFromMapValues( indexAndPowerPerSineBlockIndex ), + polynomialVariationParameter->getPolynomialVariationModel( )->getReferenceEpoch( ), + std::placeholders::_1 ); + } + default: + break; + } + } // Check non-rotational parameters. else { @@ -400,6 +426,54 @@ void SphericalHarmonicsGravityPartial::update( const double currentTime ) tidalLoveNumberPartialInterfaces_.at( i )->update( currentTime ); } } + +} + +void SphericalHarmonicsGravityPartial::wrtPolynomialGravityFieldVariations( + const std::vector< std::pair< int, int > >& cosineBlockIndices, + const std::vector< std::pair< int, int > >& sineBlockIndices, + const std::vector< std::vector< std::pair< int, int > > > powersPerCosineBlockIndex, + const std::vector< std::vector< std::pair< int, int > > > powersPerSineBlockIndex, + const double referenceEpoch, + Eigen::MatrixXd& partialDerivatives ) +{ + Eigen::MatrixXd staticCosinePartialsMatrix = Eigen::MatrixXd::Zero( cosineBlockIndices.size( ), 1 ); + Eigen::MatrixXd staticSinePartialsMatrix = Eigen::MatrixXd::Zero( sineBlockIndices.size( ), 1 ); + + calculateSphericalHarmonicGravityWrtCCoefficients( + bodyFixedSphericalPosition_, bodyReferenceRadius_( ), gravitationalParameterFunction_( ), + sphericalHarmonicCache_, + cosineBlockIndices, coordinate_conversions::getSphericalToCartesianGradientMatrix( + bodyFixedPosition_ ), fromBodyFixedToIntegrationFrameRotation_( ), staticCosinePartialsMatrix, + maximumDegree_, maximumOrder_ ); + + calculateSphericalHarmonicGravityWrtSCoefficients( + bodyFixedSphericalPosition_, bodyReferenceRadius_( ), gravitationalParameterFunction_( ), + sphericalHarmonicCache_, + sineBlockIndices, coordinate_conversions::getSphericalToCartesianGradientMatrix( + bodyFixedPosition_ ), fromBodyFixedToIntegrationFrameRotation_( ), staticSinePartialsMatrix, + maximumDegree_, maximumOrder_ ); + + partialDerivatives.setZero( ); + + for( unsigned int i = 0; i < powersPerCosineBlockIndex.size( ); i++ ) + { + for( unsigned int j = 0; j < powersPerCosineBlockIndex.at( i ).size( ); j++ ) + { + partialDerivatives.block( 0, powersPerCosineBlockIndex.at( i ).at( j ).first, 3, 1 ) += + staticCosinePartialsMatrix.block( 0, i, 3, 1 ) * std::pow( ( currentTime_ - referenceEpoch), powersPerCosineBlockIndex.at( i ).at( j ).second ); + } + } + + for( unsigned int i = 0; i < powersPerSineBlockIndex.size( ); i++ ) + { + for( unsigned int j = 0; j < powersPerSineBlockIndex.at( i ).size( ); j++ ) + { + partialDerivatives.block( 0, powersPerSineBlockIndex.at( i ).at( j ).first, 3, 1 ) += + staticSinePartialsMatrix.block( 0, i, 3, 1 ) * std::pow( ( currentTime_ - referenceEpoch), powersPerSineBlockIndex.at( i ).at( j ).second ); + } + } + } //! Function to calculate the partial of the acceleration wrt a set of cosine coefficients. @@ -415,6 +489,7 @@ void SphericalHarmonicsGravityPartial::wrtCosineCoefficientBlock( maximumDegree_, maximumOrder_ ); } + //! Function to calculate the partial of the acceleration wrt a set of sine coefficients. void SphericalHarmonicsGravityPartial::wrtSineCoefficientBlock( const std::vector< std::pair< int, int > >& blockIndices, diff --git a/src/astro/orbit_determination/estimatable_parameters/estimatableParameter.cpp b/src/astro/orbit_determination/estimatable_parameters/estimatableParameter.cpp index 039002a292..4b1ccb541c 100644 --- a/src/astro/orbit_determination/estimatable_parameters/estimatableParameter.cpp +++ b/src/astro/orbit_determination/estimatable_parameters/estimatableParameter.cpp @@ -425,6 +425,22 @@ bool isParameterTidalProperty( const EstimatebleParametersEnum parameterType ) return flag; } +//! Function to determine whether the given parameter influences a body's tidal gravity field variations. +bool isParameterNonTidalGravityFieldVariationProperty( const EstimatebleParametersEnum parameterType ) +{ + bool flag; + switch( parameterType ) + { + case polynomial_gravity_field_variation_amplitudes: + flag = true; + break; + default: + flag = false; + break; + } + return flag; +} + //! Function to determine whether the given parameter represents an arc-wise initial dynamical state. bool isParameterArcWiseInitialStateProperty( const EstimatebleParametersEnum parameterType ) { From e64441f6abca935a320d858c140fa174b9ea0e9d Mon Sep 17 00:00:00 2001 From: Dominic Dirkx Date: Sat, 6 Jul 2024 19:46:31 +0200 Subject: [PATCH 19/25] Polynomial gravity field variation partials working --- .../gravityFieldVariationParameters.h | 23 ++++++------ .../polynomialGravityFieldVariations.cpp | 4 +- .../sphericalHarmonicAccelerationPartial.cpp | 10 +++-- .../estimatableParameter.cpp | 6 +++ .../unitTestSphericalHarmonicPartials.cpp | 37 +++++++++++++++++++ 5 files changed, 64 insertions(+), 16 deletions(-) diff --git a/include/tudat/astro/orbit_determination/estimatable_parameters/gravityFieldVariationParameters.h b/include/tudat/astro/orbit_determination/estimatable_parameters/gravityFieldVariationParameters.h index 6e8cbd1b8e..9932c39e83 100644 --- a/include/tudat/astro/orbit_determination/estimatable_parameters/gravityFieldVariationParameters.h +++ b/include/tudat/astro/orbit_determination/estimatable_parameters/gravityFieldVariationParameters.h @@ -72,7 +72,7 @@ class PolynomialGravityFieldVariationsParameters: public EstimatableParameter< E } if( it.second.at( i ).second > polynomialVariationModel_->getMaximumOrder( ) || - it.second.at( i ).second < polynomialVariationModel_->getMaximumOrder( ) ) + it.second.at( i ).second < polynomialVariationModel_->getMinimumOrder( ) ) { throw std::runtime_error( "Error when estimating gravity field polynomial corrections of body " + bodyName + " of polynomial term of order " + std::to_string( it.first ) + @@ -106,7 +106,7 @@ class PolynomialGravityFieldVariationsParameters: public EstimatableParameter< E } if( it.second.at( i ).second > polynomialVariationModel_->getMaximumOrder( ) || - it.second.at( i ).second < polynomialVariationModel_->getMaximumOrder( ) ) + it.second.at( i ).second < polynomialVariationModel_->getMinimumOrder( ) ) { throw std::runtime_error( "Error when estimating gravity field polynomial corrections of body " + bodyName + " of polynomial term of order " + std::to_string( it.first ) + @@ -136,15 +136,14 @@ class PolynomialGravityFieldVariationsParameters: public EstimatableParameter< E for( unsigned int i = 0; i < cosineCorrectionIndices_.size( ); i++ ) { polynomialCorrections( i ) = polynomialVariationModel_->getCosineAmplitudesReference( ).at( std::get< 0 >( cosineCorrectionIndices_.at( i ) ) ) - ( std::get< 1 >( cosineCorrectionIndices_.at( i ) ), std::get< 2 >( cosineCorrectionIndices_.at( i ) ) ); + ( std::get< 1 >( cosineCorrectionIndices_.at( i ) ) - polynomialVariationModel_->getMinimumDegree( ), std::get< 2 >( cosineCorrectionIndices_.at( i ) ) - polynomialVariationModel_->getMinimumOrder( ) ); } - for( unsigned int i = cosineCorrectionIndices_.size( ); i < polynomialCorrections.rows( ); i++ ) + for( unsigned int i = 0; i < sineCorrectionIndices_.size( ); i++ ) { - polynomialCorrections( i ) = polynomialVariationModel_->getSineAmplitudesReference( ).at( std::get< 0 >( sineCorrectionIndices_.at( i ) ) ) - ( std::get< 1 >( sineCorrectionIndices_.at( i ) ), std::get< 2 >( cosineCorrectionIndices_.at( i ) ) ); + polynomialCorrections( i + cosineCorrectionIndices_.size( ) ) = polynomialVariationModel_->getSineAmplitudesReference( ).at( std::get< 0 >( sineCorrectionIndices_.at( i ) ) ) + ( std::get< 1 >( sineCorrectionIndices_.at( i ) ) - polynomialVariationModel_->getMinimumDegree( ), std::get< 2 >( cosineCorrectionIndices_.at( i ) ) - polynomialVariationModel_->getMinimumOrder( ) ); } - return polynomialCorrections; } @@ -164,14 +163,16 @@ class PolynomialGravityFieldVariationsParameters: public EstimatableParameter< E for( unsigned int i = 0; i < cosineCorrectionIndices_.size( ); i++ ) { - cosineVariations[ std::get< 0 >( cosineCorrectionIndices_.at( i ) ) ]( std::get< 1 >( cosineCorrectionIndices_.at( i ) ), std::get< 2 >( cosineCorrectionIndices_.at( i ) ) ) = + cosineVariations[ std::get< 0 >( cosineCorrectionIndices_.at( i ) ) ]( std::get< 1 >( cosineCorrectionIndices_.at( i ) ) - polynomialVariationModel_->getMinimumDegree( ), + std::get< 2 >( cosineCorrectionIndices_.at( i ) ) - polynomialVariationModel_->getMinimumOrder( ) ) = parameterValue( i ); } - for( unsigned int i = cosineCorrectionIndices_.size( ); i < parameterValue.rows( ); i++ ) + for( unsigned int i = 0; i < sineCorrectionIndices_.size( ); i++ ) { - sineVariations[ std::get< 0 >( sineCorrectionIndices_.at( i ) ) ]( std::get< 1 >( sineCorrectionIndices_.at( i ) ), std::get< 2 >( sineCorrectionIndices_.at( i ) ) ) = - parameterValue( i ); + sineVariations[ std::get< 0 >( sineCorrectionIndices_.at( i ) ) ]( std::get< 1 >( sineCorrectionIndices_.at( i ) ) - polynomialVariationModel_->getMinimumDegree( ), + std::get< 2 >( sineCorrectionIndices_.at( i ) ) - polynomialVariationModel_->getMinimumOrder( ) ) = + parameterValue( i + cosineCorrectionIndices_.size( ) ); } polynomialVariationModel_->resetCosineAmplitudes( cosineVariations ); polynomialVariationModel_->resetSineAmplitudes( sineVariations ); diff --git a/src/astro/gravitation/polynomialGravityFieldVariations.cpp b/src/astro/gravitation/polynomialGravityFieldVariations.cpp index 938da11a68..4846762b1c 100644 --- a/src/astro/gravitation/polynomialGravityFieldVariations.cpp +++ b/src/astro/gravitation/polynomialGravityFieldVariations.cpp @@ -87,7 +87,7 @@ PolynomialGravityFieldVariations::PolynomialGravityFieldVariations( { for( auto it : cosineAmplitudes_ ) { - if( cosineAmplitudes_.count( it.first ) != 0 ) + if( sineAmplitudes_.count( it.first ) != 0 ) { if( ( it.second.rows( ) != sineAmplitudes_.at( it.first ).rows( ) ) || ( it.second.cols( ) != sineAmplitudes_.at( it.first ).cols( ) ) ) { @@ -98,7 +98,7 @@ PolynomialGravityFieldVariations::PolynomialGravityFieldVariations( for( auto it : sineAmplitudes_ ) { - if( sineAmplitudes_.count( it.first ) != 0 ) + if( cosineAmplitudes_.count( it.first ) != 0 ) { if( ( it.second.rows( ) != cosineAmplitudes_.at( it.first ).rows( ) ) || ( it.second.cols( ) != cosineAmplitudes_.at( it.first ).cols( ) ) ) { diff --git a/src/astro/orbit_determination/acceleration_partials/sphericalHarmonicAccelerationPartial.cpp b/src/astro/orbit_determination/acceleration_partials/sphericalHarmonicAccelerationPartial.cpp index 1f53012ef3..9b5d60b07c 100644 --- a/src/astro/orbit_determination/acceleration_partials/sphericalHarmonicAccelerationPartial.cpp +++ b/src/astro/orbit_determination/acceleration_partials/sphericalHarmonicAccelerationPartial.cpp @@ -259,6 +259,8 @@ std::pair< std::function< void( Eigen::MatrixXd& ) >, int > SphericalHarmonicsGr utilities::createVectorFromMapValues( indexAndPowerPerSineBlockIndex ), polynomialVariationParameter->getPolynomialVariationModel( )->getReferenceEpoch( ), std::placeholders::_1 ); + + numberOfRows = parameter->getParameterSize( ); } default: break; @@ -437,8 +439,8 @@ void SphericalHarmonicsGravityPartial::wrtPolynomialGravityFieldVariations( const double referenceEpoch, Eigen::MatrixXd& partialDerivatives ) { - Eigen::MatrixXd staticCosinePartialsMatrix = Eigen::MatrixXd::Zero( cosineBlockIndices.size( ), 1 ); - Eigen::MatrixXd staticSinePartialsMatrix = Eigen::MatrixXd::Zero( sineBlockIndices.size( ), 1 ); + Eigen::MatrixXd staticCosinePartialsMatrix = Eigen::MatrixXd::Zero( 3, cosineBlockIndices.size( ) ); + Eigen::MatrixXd staticSinePartialsMatrix = Eigen::MatrixXd::Zero( 3, sineBlockIndices.size( ) ); calculateSphericalHarmonicGravityWrtCCoefficients( bodyFixedSphericalPosition_, bodyReferenceRadius_( ), gravitationalParameterFunction_( ), @@ -456,12 +458,14 @@ void SphericalHarmonicsGravityPartial::wrtPolynomialGravityFieldVariations( partialDerivatives.setZero( ); + int counter = 0; for( unsigned int i = 0; i < powersPerCosineBlockIndex.size( ); i++ ) { for( unsigned int j = 0; j < powersPerCosineBlockIndex.at( i ).size( ); j++ ) { partialDerivatives.block( 0, powersPerCosineBlockIndex.at( i ).at( j ).first, 3, 1 ) += staticCosinePartialsMatrix.block( 0, i, 3, 1 ) * std::pow( ( currentTime_ - referenceEpoch), powersPerCosineBlockIndex.at( i ).at( j ).second ); + counter++; } } @@ -469,7 +473,7 @@ void SphericalHarmonicsGravityPartial::wrtPolynomialGravityFieldVariations( { for( unsigned int j = 0; j < powersPerSineBlockIndex.at( i ).size( ); j++ ) { - partialDerivatives.block( 0, powersPerSineBlockIndex.at( i ).at( j ).first, 3, 1 ) += + partialDerivatives.block( 0, powersPerSineBlockIndex.at( i ).at( j ).first + counter, 3, 1 ) += staticSinePartialsMatrix.block( 0, i, 3, 1 ) * std::pow( ( currentTime_ - referenceEpoch), powersPerSineBlockIndex.at( i ).at( j ).second ); } } diff --git a/src/astro/orbit_determination/estimatable_parameters/estimatableParameter.cpp b/src/astro/orbit_determination/estimatable_parameters/estimatableParameter.cpp index 4b1ccb541c..25c4ba6035 100644 --- a/src/astro/orbit_determination/estimatable_parameters/estimatableParameter.cpp +++ b/src/astro/orbit_determination/estimatable_parameters/estimatableParameter.cpp @@ -150,6 +150,9 @@ std::string getParameterTypeString( const EstimatebleParametersEnum parameterTyp case custom_estimated_parameter: parameterDescription = " Custom parameter "; break; + case polynomial_gravity_field_variation_amplitudes: + parameterDescription = " Polynomial gravity field variations "; + break; default: std::string errorMessage = "Error when getting parameter string, did not recognize parameter " + std::to_string( parameterType ); @@ -306,6 +309,9 @@ bool isDoubleParameter( const EstimatebleParametersEnum parameterType ) case custom_estimated_parameter: isDoubleParameter = false; break; + case polynomial_gravity_field_variation_amplitudes: + isDoubleParameter = false; + break; default: throw std::runtime_error( "Error, parameter type " + std::to_string( parameterType ) + " not found when getting parameter type" ); diff --git a/tests/src/astro/orbit_determination/acceleration_partials/unitTestSphericalHarmonicPartials.cpp b/tests/src/astro/orbit_determination/acceleration_partials/unitTestSphericalHarmonicPartials.cpp index 04916f7748..ba5ac46c2e 100644 --- a/tests/src/astro/orbit_determination/acceleration_partials/unitTestSphericalHarmonicPartials.cpp +++ b/tests/src/astro/orbit_determination/acceleration_partials/unitTestSphericalHarmonicPartials.cpp @@ -427,6 +427,16 @@ std::vector< std::shared_ptr< GravityFieldVariationSettings > > getEarthGravityF std::shared_ptr< GravityFieldVariationSettings > singleGravityFieldVariation = std::make_shared< BasicSolidBodyGravityFieldVariationSettings >( deformingBodies, loveNumbers ); gravityFieldVariations.push_back( singleGravityFieldVariation ); + + std::map cosineAmplitudes; + cosineAmplitudes[ 1 ] = Eigen::Matrix< double, 3, 5 >::Zero( ); + cosineAmplitudes[ 2 ] = Eigen::Matrix< double, 3, 5 >::Zero( ); + std::map sineAmplitudes; + sineAmplitudes[ 1 ] = Eigen::Matrix< double, 3, 5 >::Zero( ); + std::shared_ptr< GravityFieldVariationSettings > polynomialGravityFieldVariations = + std::make_shared< PolynomialGravityFieldVariationsSettings >( + cosineAmplitudes, sineAmplitudes, -1.0E4, 2, 0 ); + gravityFieldVariations.push_back( polynomialGravityFieldVariations ); return gravityFieldVariations; } @@ -578,6 +588,22 @@ BOOST_AUTO_TEST_CASE( testSphericalHarmonicAccelerationPartial ) "Earth", 3, "", false ) ); parameterNames.push_back( std::make_shared< SingleDegreeVariableTidalLoveNumberEstimatableParameterSettings >( "Earth", 3, std::vector< int >{ 0, 3 }, "", true ) ); + parameterNames.push_back( std::make_shared< SingleDegreeVariableTidalLoveNumberEstimatableParameterSettings >( + "Earth", 3, std::vector< int >{ 0, 3 }, "", true ) ); + + std::map > > cosineBlockIndicesPerPower; + cosineBlockIndicesPerPower[ 1 ].push_back( std::make_pair( 2, 0 ) ); + cosineBlockIndicesPerPower[ 1 ].push_back( std::make_pair( 3, 2 ) ); + cosineBlockIndicesPerPower[ 1 ].push_back( std::make_pair( 4, 4 ) ); + cosineBlockIndicesPerPower[ 2 ].push_back( std::make_pair( 2, 0 ) ); + cosineBlockIndicesPerPower[ 2 ].push_back( std::make_pair( 2, 2 ) ); + + std::map > > sineBlockIndicesPerPower; + sineBlockIndicesPerPower[ 1 ].push_back( std::make_pair( 3, 2 ) ); + parameterNames.push_back( std::make_shared< PolynomialGravityFieldVariationEstimatableParameterSettings >( + "Earth", cosineBlockIndicesPerPower, sineBlockIndicesPerPower ) ); + + std::shared_ptr< estimatable_parameters::EstimatableParameterSet< double > > parameterSet = createParametersToEstimate( parameterNames, bodies ); @@ -779,11 +805,20 @@ BOOST_AUTO_TEST_CASE( testSphericalHarmonicAccelerationPartial ) vectorParametersIterator->second->setParameterValue( nominalTidalParameter - Eigen::VectorXd::Constant( nominalTidalParameter.rows( ), 1.0 ) ); earthGravityField->update( testTime ); + Eigen::MatrixXd downperturbedCosineCoefficients = earthGravityField->getCosineCoefficients( ).block( 0, 0, 3, 3 ); Eigen::MatrixXd downperturbedSineCoefficients = earthGravityField->getSineCoefficients( ).block( 0, 0, 3, 3 ); + vectorParametersIterator++; + vectorParametersIterator++; + + Eigen::MatrixXd partialWrtPolynomialVariations = accelerationPartial->wrtParameter( + vectorParametersIterator->second ); + Eigen::MatrixXd testPartialWrtPolynomialVariations = calculateAccelerationWrtParameterPartials( + vectorParametersIterator->second, gravitationalAcceleration, Eigen::VectorXd::Constant( 6, 1.0E-6 ), sphericalHarmonicFieldUpdate ); + TUDAT_CHECK_MATRIX_CLOSE_FRACTION( testPartialWrtVehiclePosition, partialWrtVehiclePosition, 1.0E-6 ); TUDAT_CHECK_MATRIX_CLOSE_FRACTION( testPartialWrtVehicleVelocity, partialWrtVehicleVelocity, 1.0E-6 ); TUDAT_CHECK_MATRIX_CLOSE_FRACTION( testPartialWrtEarthPosition, partialWrtEarthPosition, 1.0E-6 ); @@ -820,6 +855,8 @@ BOOST_AUTO_TEST_CASE( testSphericalHarmonicAccelerationPartial ) TUDAT_CHECK_MATRIX_CLOSE_FRACTION( partialWrtDegreeThreeLoveNumber, testPartialWrtDegreeThreeLoveNumber, 1.0E-6 ); TUDAT_CHECK_MATRIX_CLOSE_FRACTION( partialWrtComplexDegreeThreeLoveNumberAtSeparateOrder, testPartialWrtComplexDegreeThreeLoveNumberAtSeparateOrder, 1.0E-6 ); + TUDAT_CHECK_MATRIX_CLOSE_FRACTION( partialWrtPolynomialVariations, testPartialWrtPolynomialVariations, 1.0E-12 ); + } //! Unit test to check working onf spherical harmonic state partial for synchronously rotating body (and rotation depending on state) From 0d7f1b68412e1750fabb9fda4467e901b226af1b Mon Sep 17 00:00:00 2001 From: Dominic Dirkx Date: Sat, 6 Jul 2024 20:20:45 +0200 Subject: [PATCH 20/25] Added additional parameter factory functions --- .../estimatableParameterSettings.h | 44 ++++++++++++++++--- 1 file changed, 39 insertions(+), 5 deletions(-) diff --git a/include/tudat/simulation/estimation_setup/estimatableParameterSettings.h b/include/tudat/simulation/estimation_setup/estimatableParameterSettings.h index e23ba0a70d..4012063ddb 100644 --- a/include/tudat/simulation/estimation_setup/estimatableParameterSettings.h +++ b/include/tudat/simulation/estimation_setup/estimatableParameterSettings.h @@ -1001,8 +1001,8 @@ class PolynomialGravityFieldVariationEstimatableParameterSettings: public Estima */ PolynomialGravityFieldVariationEstimatableParameterSettings( const std::string &associatedBody, - const std::map > > &cosineBlockIndicesPerPower, - const std::map > > &sineBlockIndicesPerPower ) : + const std::map > >& cosineBlockIndicesPerPower, + const std::map > >& sineBlockIndicesPerPower ) : EstimatableParameterSettings( associatedBody, polynomial_gravity_field_variation_amplitudes ), cosineBlockIndicesPerPower_( cosineBlockIndicesPerPower ), sineBlockIndicesPerPower_( sineBlockIndicesPerPower ){ } @@ -1464,6 +1464,43 @@ inline std::shared_ptr< EstimatableParameterSettings > yarkovskyParameter( const return std::make_shared< EstimatableParameterSettings >( bodyName, yarkovsky_parameter, centralBodyName ); } +inline std::shared_ptr< EstimatableParameterSettings > polynomialGravityFieldVariationParameter( + const std::string bodyName, + const std::map > >& cosineBlockIndicesPerPower, + const std::map > >& sineBlockIndicesPerPower ) +{ + return std::make_shared< PolynomialGravityFieldVariationEstimatableParameterSettings >( bodyName, cosineBlockIndicesPerPower, sineBlockIndicesPerPower ); +} + +inline std::shared_ptr< EstimatableParameterSettings > polynomialSinglePowerGravityFieldVariationParameter( + const std::string bodyName, + const int power, + const std::vector >& cosineBlockIndices, + const std::vector >& sineBlockIndices ) +{ + return std::make_shared< PolynomialGravityFieldVariationEstimatableParameterSettings >( + bodyName, + std::map > >( { { power, cosineBlockIndices } } ), + std::map > >( { { power, sineBlockIndices } } ) ); +} + +inline std::shared_ptr< EstimatableParameterSettings > polynomialSinglePowerFullBlockGravityFieldVariationParameter( + const std::string bodyName, + const int power, + const int minimumDegree, + const int minimumOrder, + const int maximumDegree, + const int maximumOrder ) +{ + + std::vector< std::pair< int, int > > blockIndices = getSphericalHarmonicBlockIndices( minimumDegree, minimumOrder, maximumDegree, maximumOrder ); + + return std::make_shared< PolynomialGravityFieldVariationEstimatableParameterSettings >( + bodyName, + std::map > >( { { power, blockIndices } } ), + std::map > >( { { power, blockIndices } } ) ); +} + inline std::shared_ptr< EstimatableParameterSettings > customParameterSettings( const std::string& customId, const int parameterSize, @@ -1474,9 +1511,6 @@ inline std::shared_ptr< EstimatableParameterSettings > customParameterSettings( customId, parameterSize, getParameterFunction, setParameterFunction ); } - - - } // namespace estimatable_parameters } // namespace tudat From 9aaa28d80bf749a3283734f4b5ba7a5dc1c51f32 Mon Sep 17 00:00:00 2001 From: Dominic Dirkx Date: Sat, 6 Jul 2024 21:50:34 +0200 Subject: [PATCH 21/25] Estimation of periodic gravity field varitions now added --- .../periodicGravityFieldVariations.h | 67 +++++- .../sphericalHarmonicAccelerationPartial.h | 9 + .../estimatableParameter.h | 3 +- .../gravityFieldVariationParameters.h | 209 +++++++++++++++++- .../createGravityFieldVariations.h | 80 ++++--- .../createEstimatableParameters.h | 57 +++++ .../estimatableParameterSettings.h | 23 ++ .../periodicGravityFieldVariations.cpp | 43 ++-- .../sphericalHarmonicAccelerationPartial.cpp | 78 +++++++ .../estimatableParameter.cpp | 9 + .../createGravityFieldVariations.cpp | 26 ++- .../unitTestGravityFieldVariations.cpp | 86 ++++--- .../unitTestSphericalHarmonicPartials.cpp | 73 ++++++ 13 files changed, 660 insertions(+), 103 deletions(-) diff --git a/include/tudat/astro/gravitation/periodicGravityFieldVariations.h b/include/tudat/astro/gravitation/periodicGravityFieldVariations.h index ac524ad3df..1034a46d93 100644 --- a/include/tudat/astro/gravitation/periodicGravityFieldVariations.h +++ b/include/tudat/astro/gravitation/periodicGravityFieldVariations.h @@ -25,10 +25,11 @@ class PeriodicGravityFieldVariations: public GravityFieldVariations { public: PeriodicGravityFieldVariations( - const std::vector< Eigen::MatrixXd >& cosineAmplitudes, - const std::vector< Eigen::MatrixXd >& sineAmplitudes, + const std::vector< Eigen::MatrixXd >& cosineShAmplitudesCosineTime, + const std::vector< Eigen::MatrixXd >& cosineShAmplitudesSineTime, + const std::vector< Eigen::MatrixXd >& sineShAmplitudesCosineTime, + const std::vector< Eigen::MatrixXd >& sineShAmplitudesSineTime, const std::vector< double >& frequencies, - const std::vector< double >& phases, const double referenceEpoch, const int minimumDegree = 2, const int minimumOrder = 0 ); @@ -39,15 +40,67 @@ class PeriodicGravityFieldVariations: public GravityFieldVariations const double time ); + std::vector< Eigen::MatrixXd > getCosineShAmplitudesCosineTime( ) + { + return cosineShAmplitudesCosineTime_; + } + + std::vector< Eigen::MatrixXd > getCosineShAmplitudesSineTime( ) + { + return cosineShAmplitudesSineTime_; + } + + std::vector< Eigen::MatrixXd > getSineShAmplitudesCosineTime( ) + { + return sineShAmplitudesCosineTime_; + } + + std::vector< Eigen::MatrixXd > getSineShAmplitudesSineTime( ) + { + return sineShAmplitudesSineTime_; + } + + void resetCosineShAmplitudesCosineTime( const std::vector< Eigen::MatrixXd >& cosineShAmplitudesCosineTime ) + { + cosineShAmplitudesCosineTime_ = cosineShAmplitudesCosineTime; + } + + void resetCosineShAmplitudesSineTime( const std::vector< Eigen::MatrixXd >& cosineShAmplitudesSineTime ) + { + cosineShAmplitudesSineTime_ = cosineShAmplitudesSineTime; + } + + void resetSineShAmplitudesCosineTime( const std::vector< Eigen::MatrixXd >& sineShAmplitudesCosineTime ) + { + sineShAmplitudesCosineTime_ = sineShAmplitudesCosineTime; + } + + void resetSineShAmplitudesSineTime( const std::vector< Eigen::MatrixXd >& sineShAmplitudesSineTime ) + { + sineShAmplitudesSineTime_ = sineShAmplitudesSineTime; + } + + std::vector< double > getFrequencies( ) + { + return frequencies_; + } + + double getReferenceEpoch( ) + { + return referenceEpoch_; + } + protected: - const std::vector< Eigen::MatrixXd > cosineAmplitudes_; + std::vector< Eigen::MatrixXd > cosineShAmplitudesCosineTime_; - const std::vector< Eigen::MatrixXd > sineAmplitudes_; + std::vector< Eigen::MatrixXd > cosineShAmplitudesSineTime_; - const std::vector< double > frequencies_; + std::vector< Eigen::MatrixXd > sineShAmplitudesCosineTime_; - const std::vector< double > phases_; + std::vector< Eigen::MatrixXd > sineShAmplitudesSineTime_; + + const std::vector< double > frequencies_; const double referenceEpoch_; diff --git a/include/tudat/astro/orbit_determination/acceleration_partials/sphericalHarmonicAccelerationPartial.h b/include/tudat/astro/orbit_determination/acceleration_partials/sphericalHarmonicAccelerationPartial.h index d7191a4948..6404fb6ed2 100644 --- a/include/tudat/astro/orbit_determination/acceleration_partials/sphericalHarmonicAccelerationPartial.h +++ b/include/tudat/astro/orbit_determination/acceleration_partials/sphericalHarmonicAccelerationPartial.h @@ -340,6 +340,15 @@ class SphericalHarmonicsGravityPartial: public AccelerationPartial const double referenceEpoch, Eigen::MatrixXd& partialDerivatives ); + void wrtPeriodicGravityFieldVariations( + const std::vector< std::pair< int, int > >& cosineBlockIndices, + const std::vector< std::pair< int, int > >& sineBlockIndices, + const std::vector< std::vector< std::pair< int, int > > > powersPerCosineBlockIndex, + const std::vector< std::vector< std::pair< int, int > > > powersPerSineBlockIndex, + const std::vector< double >& frequencies, + const double referenceEpoch, + Eigen::MatrixXd& partialDerivatives ); + //! Function to calculate the partial of the acceleration wrt a set of cosine coefficients. /*! * Function to calculate the partial of the acceleration wrt a set of cosine coefficients. diff --git a/include/tudat/astro/orbit_determination/estimatable_parameters/estimatableParameter.h b/include/tudat/astro/orbit_determination/estimatable_parameters/estimatableParameter.h index 827465232c..0c4e86e642 100644 --- a/include/tudat/astro/orbit_determination/estimatable_parameters/estimatableParameter.h +++ b/include/tudat/astro/orbit_determination/estimatable_parameters/estimatableParameter.h @@ -78,7 +78,8 @@ enum EstimatebleParametersEnum yarkovsky_parameter, custom_estimated_parameter, reference_point_position, - polynomial_gravity_field_variation_amplitudes + polynomial_gravity_field_variation_amplitudes, + periodic_gravity_field_variation_amplitudes }; std::string getParameterTypeString( const EstimatebleParametersEnum parameterType ); diff --git a/include/tudat/astro/orbit_determination/estimatable_parameters/gravityFieldVariationParameters.h b/include/tudat/astro/orbit_determination/estimatable_parameters/gravityFieldVariationParameters.h index 9932c39e83..227d26ea34 100644 --- a/include/tudat/astro/orbit_determination/estimatable_parameters/gravityFieldVariationParameters.h +++ b/include/tudat/astro/orbit_determination/estimatable_parameters/gravityFieldVariationParameters.h @@ -24,6 +24,7 @@ #include "tudat/astro/gravitation/polynomialGravityFieldVariations.h" #include "tudat/astro/orbit_determination/estimatable_parameters/estimatableParameter.h" +#include "tudat/astro/gravitation/periodicGravityFieldVariations.h" namespace tudat { @@ -57,7 +58,7 @@ class PolynomialGravityFieldVariationsParameters: public EstimatableParameter< E { if( cosineVariations.count( it.first ) == 0 ) { - throw std::runtime_error( "Error when estimationg gravity field polynomial corrections of body " + bodyName + + throw std::runtime_error( "Error when estimating gravity field polynomial corrections of body " + bodyName + ", not polynomial term of order " + std::to_string( it.first ) + " found for cosine coefficients" ); } @@ -91,7 +92,7 @@ class PolynomialGravityFieldVariationsParameters: public EstimatableParameter< E { if( sineVariations.count( it.first ) == 0 ) { - throw std::runtime_error( "Error when estimationg gravity field polynomial corrections of body " + bodyName + + throw std::runtime_error( "Error when estimating gravity field polynomial corrections of body " + bodyName + ", not polynomial term of order " + std::to_string( it.first ) + " found for sine coefficients" ); } @@ -218,6 +219,210 @@ class PolynomialGravityFieldVariationsParameters: public EstimatableParameter< E }; +class PeriodicGravityFieldVariationsParameters: public EstimatableParameter< Eigen::VectorXd > +{ + +public: + //! Constructor. + /*! + * Constructor taking parameter name and associated body. All parameters are identified by a these two variables. + * Any additional information that may be required for uniquely defining a parameter is to be defined in the derived class. + * \param parameterName Enum value defining the type of the parameter. + * \param associatedBody Name of body associated with patameters + * \param pointOnBodyId Reference point on body associated with parameter (empty by default). + */ + PeriodicGravityFieldVariationsParameters( + const std::shared_ptr periodicVariationModel, + const std::map > > &cosineBlockIndicesPerPeriod, + const std::map > > &sineBlockIndicesPerPeriod, + const std::string &bodyName ) : + EstimatableParameter( periodic_gravity_field_variation_amplitudes, bodyName ), + periodicVariationModel_( periodicVariationModel ) + { + std::vector< Eigen::MatrixXd > cosineVariations = periodicVariationModel->getCosineShAmplitudesCosineTime( ); + int cosineIndexCounter = 0; + for( auto it : cosineBlockIndicesPerPeriod ) + { + if( static_cast< int >( cosineVariations.size( ) ) < it.first ) + { + throw std::runtime_error( "Error when estimating gravity field periodic corrections of body " + bodyName + + ", not periodic term of order " + std::to_string( it.first ) + " found for cosine coefficients" ); + } + + for( unsigned int i = 0; i < it.second.size( ); i++ ) + { + if( it.second.at( i ).first > periodicVariationModel_->getMaximumDegree( ) || + it.second.at( i ).first < periodicVariationModel_->getMinimumDegree( ) ) + { + throw std::runtime_error( "Error when estimating gravity field periodic corrections of body " + bodyName + + " of periodic term of period index " + std::to_string( it.first ) + + ", no cosine coefficient variation of degree " + std::to_string( it.second.at( i ).first ) + " found." ); + } + + if( it.second.at( i ).second > periodicVariationModel_->getMaximumOrder( ) || + it.second.at( i ).second < periodicVariationModel_->getMinimumOrder( ) ) + { + throw std::runtime_error( "Error when estimating gravity field periodic corrections of body " + bodyName + + " of periodic term of period index " + std::to_string( it.first ) + + ", no cosine coefficient variation of order " + std::to_string( it.second.at( i ).second ) + " found." ); + } + cosineCorrectionIndices_.push_back( std::make_tuple( it.first, it.second.at( i ).first, it.second.at( i ).second ) ); + indexAndPowerPerCosineBlockIndex_[ std::make_pair( it.second.at( i ).first, it.second.at( i ).second ) ].push_back( + std::make_pair( cosineIndexCounter, it.first ) ); + cosineIndexCounter++; + } + } + + std::vector< Eigen::MatrixXd > sineVariations = periodicVariationModel->getSineShAmplitudesCosineTime( ); + int sineIndexCounter = 0; + for( auto it : sineBlockIndicesPerPeriod ) + { + if( static_cast< int >( sineVariations.size( ) ) < it.first ) + { + throw std::runtime_error( "Error when estimating gravity field periodic corrections of body " + bodyName + + ", not periodic term of order " + std::to_string( it.first ) + " found for sine coefficients" ); + } + + for( unsigned int i = 0; i < it.second.size( ); i++ ) + { + if( it.second.at( i ).first > periodicVariationModel_->getMaximumDegree( ) || + it.second.at( i ).first < periodicVariationModel_->getMinimumDegree( ) ) + { + throw std::runtime_error( "Error when estimating gravity field periodic corrections of body " + bodyName + + " of periodic term of period index " + std::to_string( it.first ) + + ", no sine coefficient variation of degree " + std::to_string( it.second.at( i ).first ) + " found." ); + } + + if( it.second.at( i ).second > periodicVariationModel_->getMaximumOrder( ) || + it.second.at( i ).second < periodicVariationModel_->getMinimumOrder( ) ) + { + throw std::runtime_error( "Error when estimating gravity field periodic corrections of body " + bodyName + + " of periodic term of period index " + std::to_string( it.first ) + + ", no sine coefficient variation of order " + std::to_string( it.second.at( i ).second ) + " found." ); + } + sineCorrectionIndices_.push_back( std::make_tuple( it.first, it.second.at( i ).first, it.second.at( i ).second ) ); + indexAndPowerPerSineBlockIndex_[ std::make_pair( it.second.at( i ).first, it.second.at( i ).second ) ].push_back( + std::make_pair( sineIndexCounter, it.first ) ); + sineIndexCounter++; + } + } + } + + //! Virtual destructor. + ~PeriodicGravityFieldVariationsParameters( ) { } + + //! Pure virtual function to retrieve the value of the parameter + /*! + * Pure virtual function to retrieve the value of the parameter + * \return Current value of parameter. + */ + Eigen::VectorXd getParameterValue( ) + { + Eigen::VectorXd periodicCorrections = Eigen::VectorXd::Zero( getParameterSize( ) ); + for( unsigned int i = 0; i < cosineCorrectionIndices_.size( ); i++ ) + { + int periodIndex = std::get< 0 >( cosineCorrectionIndices_.at( i ) ); + int degree = std::get< 1 >( cosineCorrectionIndices_.at( i ) ) - periodicVariationModel_->getMinimumDegree( ); + int order = std::get< 2 >( cosineCorrectionIndices_.at( i ) ) - periodicVariationModel_->getMinimumOrder( ); + + periodicCorrections( 2 * i ) = periodicVariationModel_->getCosineShAmplitudesCosineTime( ).at( periodIndex )( degree, order ); + periodicCorrections( 2 * i + 1 ) = periodicVariationModel_->getCosineShAmplitudesSineTime( ).at( periodIndex )( degree, order ); + } + + for( unsigned int i = 0; i < sineCorrectionIndices_.size( ); i++ ) + { + int periodIndex = std::get< 0 >( sineCorrectionIndices_.at( i ) ); + int degree = std::get< 1 >( sineCorrectionIndices_.at( i ) ) - periodicVariationModel_->getMinimumDegree( ); + int order = std::get< 2 >( sineCorrectionIndices_.at( i ) ) - periodicVariationModel_->getMinimumOrder( ); + + periodicCorrections( 2 * ( i + cosineCorrectionIndices_.size( ) ) ) = periodicVariationModel_->getSineShAmplitudesCosineTime( ).at( periodIndex )( degree, order ); + periodicCorrections( 2 * ( i + cosineCorrectionIndices_.size( ) ) + 1 ) = periodicVariationModel_->getSineShAmplitudesSineTime( ).at( periodIndex )( degree, order ); + } + return periodicCorrections; + } + + //! Pure virtual function to (re)set the value of the parameter. + /*! + * Pure virtual function to (re)set the value of the parameter. + * \param parameterValue to which the parameter is to be set. + */ + void setParameterValue( const Eigen::VectorXd parameterValue ) + { + if( parameterValue.rows( ) != getParameterSize( ) ) + { + throw std::runtime_error( "Error when resetting periodic gravity field variation parameter; sizes are incompatible" ); + } + std::vector< Eigen::MatrixXd > cosineVariationsCosineTime = periodicVariationModel_->getCosineShAmplitudesCosineTime( ); + std::vector< Eigen::MatrixXd > cosineVariationsSineTime = periodicVariationModel_->getCosineShAmplitudesSineTime( ); + std::vector< Eigen::MatrixXd > sineVariationsCosineTime = periodicVariationModel_->getSineShAmplitudesCosineTime( ); + std::vector< Eigen::MatrixXd > sineVariationsSineTime = periodicVariationModel_->getSineShAmplitudesSineTime( ); + + + for( unsigned int i = 0; i < cosineCorrectionIndices_.size( ); i++ ) + { + int periodIndex = std::get< 0 >( cosineCorrectionIndices_.at( i ) ); + int degree = std::get< 1 >( cosineCorrectionIndices_.at( i ) ) - periodicVariationModel_->getMinimumDegree( ); + int order = std::get< 2 >( cosineCorrectionIndices_.at( i ) ) - periodicVariationModel_->getMinimumOrder( ); + + cosineVariationsCosineTime[ periodIndex ]( degree, order ) = parameterValue( 2 * i ); + cosineVariationsSineTime[ periodIndex ]( degree, order ) = parameterValue( 2 * i + 1 ); + } + + for( unsigned int i = 0; i < sineCorrectionIndices_.size( ); i++ ) + { + int periodIndex = std::get< 0 >( sineCorrectionIndices_.at( i ) ); + int degree = std::get< 1 >( sineCorrectionIndices_.at( i ) ) - periodicVariationModel_->getMinimumDegree( ); + int order = std::get< 2 >( sineCorrectionIndices_.at( i ) ) - periodicVariationModel_->getMinimumOrder( ); + + sineVariationsCosineTime[ periodIndex ]( degree, order ) = parameterValue( 2 * ( i + cosineCorrectionIndices_.size( ) ) ); + sineVariationsSineTime[ periodIndex ]( degree, order ) = parameterValue( 2 * ( i + cosineCorrectionIndices_.size( ) ) + 1 ); + } + periodicVariationModel_->resetCosineShAmplitudesCosineTime( cosineVariationsCosineTime ); + periodicVariationModel_->resetCosineShAmplitudesSineTime( cosineVariationsSineTime ); + periodicVariationModel_->resetSineShAmplitudesCosineTime( sineVariationsCosineTime ); + periodicVariationModel_->resetSineShAmplitudesSineTime( sineVariationsSineTime ); + + } + + //! Function to retrieve the size of the parameter + /*! + * Pure virtual function to retrieve the size of the parameter (i.e. 1 for double parameters) + * \return Size of parameter value. + */ + int getParameterSize( ) + { + return 2 * ( cosineCorrectionIndices_.size( ) + sineCorrectionIndices_.size( ) ); + } + + std::map< std::pair< int, int >, std::vector< std::pair< int, int > > > getIndexAndPowerPerCosineBlockIndex( ) + { + return indexAndPowerPerCosineBlockIndex_; + } + + std::map< std::pair< int, int >, std::vector< std::pair< int, int > > > getIndexAndPowerPerSineBlockIndex( ) + { + return indexAndPowerPerSineBlockIndex_; + } + std::shared_ptr< gravitation::PeriodicGravityFieldVariations > getPeriodicVariationModel( ) + { + return periodicVariationModel_; + } + + +protected: + + std::shared_ptr< gravitation::PeriodicGravityFieldVariations > periodicVariationModel_; + + std::vector< std::tuple< int, int, int > > cosineCorrectionIndices_; + + std::vector< std::tuple< int, int, int > > sineCorrectionIndices_; + + std::map< std::pair< int, int >, std::vector< std::pair< int, int > > > indexAndPowerPerCosineBlockIndex_; + + std::map< std::pair< int, int >, std::vector< std::pair< int, int > > > indexAndPowerPerSineBlockIndex_; + +}; + } // namespace estimatable_parameters } // namespace tudat diff --git a/include/tudat/simulation/environment_setup/createGravityFieldVariations.h b/include/tudat/simulation/environment_setup/createGravityFieldVariations.h index f416ed256b..be70f42d94 100644 --- a/include/tudat/simulation/environment_setup/createGravityFieldVariations.h +++ b/include/tudat/simulation/environment_setup/createGravityFieldVariations.h @@ -270,40 +270,47 @@ class PeriodicGravityFieldVariationsSettings: public GravityFieldVariationSettin public: PeriodicGravityFieldVariationsSettings( - const std::vector< Eigen::MatrixXd >& cosineAmplitudes, - const std::vector< Eigen::MatrixXd >& sineAmplitudes, + const std::vector< Eigen::MatrixXd >& cosineShAmplitudesCosineTime, + const std::vector< Eigen::MatrixXd >& cosineShAmplitudesSineTime, + const std::vector< Eigen::MatrixXd >& sineShAmplitudesCosineTime, + const std::vector< Eigen::MatrixXd >& sineShAmplitudesSineTime, const std::vector< double >& frequencies, - const std::vector< double >& phases, const double referenceEpoch, const int minimumDegree = 2, const int minimumOrder = 0 ): GravityFieldVariationSettings( gravitation::periodic_variation ), - cosineAmplitudes_( cosineAmplitudes ), - sineAmplitudes_( sineAmplitudes ), - frequencies_( frequencies ), - phases_( phases ), - referenceEpoch_( referenceEpoch ), - minimumDegree_( minimumDegree ), - minimumOrder_( minimumOrder ){ } - - std::vector< Eigen::MatrixXd > getCosineAmplitudes( ) + cosineShAmplitudesCosineTime_( cosineShAmplitudesCosineTime ), + cosineShAmplitudesSineTime_( cosineShAmplitudesSineTime ), + sineShAmplitudesCosineTime_( sineShAmplitudesCosineTime ), + sineShAmplitudesSineTime_( sineShAmplitudesSineTime ), + frequencies_( frequencies ), + referenceEpoch_( referenceEpoch ), + minimumDegree_( minimumDegree ), + minimumOrder_( minimumOrder ){ } + + std::vector< Eigen::MatrixXd > getCosineShAmplitudesCosineTime( ) { - return cosineAmplitudes_; + return cosineShAmplitudesCosineTime_; } - std::vector< Eigen::MatrixXd > getSineAmplitudes( ) + std::vector< Eigen::MatrixXd > getCosineShAmplitudesSineTime( ) { - return sineAmplitudes_; + return cosineShAmplitudesSineTime_; } - std::vector< double > getFrequencies( ) + std::vector< Eigen::MatrixXd > getSineShAmplitudesCosineTime( ) { - return frequencies_; + return sineShAmplitudesCosineTime_; } - std::vector< double > getPhases( ) + std::vector< Eigen::MatrixXd > getSineShAmplitudesSineTime( ) { - return phases_; + return sineShAmplitudesSineTime_; + } + + std::vector< double > getFrequencies( ) + { + return frequencies_; } double getReferenceEpoch( ) @@ -323,13 +330,15 @@ class PeriodicGravityFieldVariationsSettings: public GravityFieldVariationSettin private: - std::vector< Eigen::MatrixXd > cosineAmplitudes_; + std::vector< Eigen::MatrixXd > cosineShAmplitudesCosineTime_; - std::vector< Eigen::MatrixXd > sineAmplitudes_; + std::vector< Eigen::MatrixXd > cosineShAmplitudesSineTime_; - std::vector< double > frequencies_; + std::vector< Eigen::MatrixXd > sineShAmplitudesCosineTime_; + + std::vector< Eigen::MatrixXd > sineShAmplitudesSineTime_; - std::vector< double > phases_; + std::vector< double > frequencies_; double referenceEpoch_; @@ -544,32 +553,35 @@ inline std::shared_ptr< GravityFieldVariationSettings > tabulatedGravityFieldVar } inline std::shared_ptr< GravityFieldVariationSettings > periodicGravityFieldVariationsSettings( - const std::vector< Eigen::MatrixXd >& cosineAmplitudes, - const std::vector< Eigen::MatrixXd >& sineAmplitudes, + const std::vector< Eigen::MatrixXd >& cosineShAmplitudesCosineTime, + const std::vector< Eigen::MatrixXd >& cosineShAmplitudesSineTime, + const std::vector< Eigen::MatrixXd >& sineShAmplitudesCosineTime, + const std::vector< Eigen::MatrixXd >& sineShAmplitudesSineTime, const std::vector< double >& frequencies, - const std::vector< double >& phases, const double referenceEpoch, const int minimumDegree = 2, const int minimumOrder = 0 ) { return std::make_shared< PeriodicGravityFieldVariationsSettings >( - cosineAmplitudes, sineAmplitudes, frequencies, phases, referenceEpoch, minimumDegree, minimumOrder ); + cosineShAmplitudesCosineTime, cosineShAmplitudesSineTime, sineShAmplitudesCosineTime, sineShAmplitudesSineTime, frequencies, referenceEpoch, minimumDegree, minimumOrder ); } -inline std::shared_ptr< GravityFieldVariationSettings > periodicGravityFieldVariationsSettings( - const Eigen::MatrixXd& cosineAmplitudes, - const Eigen::MatrixXd& sineAmplitudes, +inline std::shared_ptr< GravityFieldVariationSettings > periodicGravityFieldVariationsSettingsSingleFrequency( + const Eigen::MatrixXd& cosineShAmplitudesCosineTime, + const Eigen::MatrixXd& cosineShAmplitudesSineTime, + const Eigen::MatrixXd& sineShAmplitudesCosineTime, + const Eigen::MatrixXd& sineShAmplitudesSineTime, const double frequency, - const double phase, const double referenceEpoch, const int minimumDegree = 2, const int minimumOrder = 0 ) { return std::make_shared< PeriodicGravityFieldVariationsSettings >( - std::vector< Eigen::MatrixXd >( { cosineAmplitudes } ), - std::vector< Eigen::MatrixXd >( { sineAmplitudes } ), + std::vector< Eigen::MatrixXd >( { cosineShAmplitudesCosineTime } ), + std::vector< Eigen::MatrixXd >( { cosineShAmplitudesSineTime } ), + std::vector< Eigen::MatrixXd >( { sineShAmplitudesCosineTime } ), + std::vector< Eigen::MatrixXd >( { sineShAmplitudesSineTime } ), std::vector< double >( { frequency } ), - std::vector< double >( { phase } ), referenceEpoch, minimumDegree, minimumOrder ); } diff --git a/include/tudat/simulation/estimation_setup/createEstimatableParameters.h b/include/tudat/simulation/estimation_setup/createEstimatableParameters.h index c8052b443c..c1503202d6 100644 --- a/include/tudat/simulation/estimation_setup/createEstimatableParameters.h +++ b/include/tudat/simulation/estimation_setup/createEstimatableParameters.h @@ -2039,6 +2039,63 @@ std::shared_ptr< estimatable_parameters::EstimatableParameter< Eigen::VectorXd > } break; } + case periodic_gravity_field_variation_amplitudes: + { + std::shared_ptr< PeriodicGravityFieldVariationEstimatableParameterSettings > gravityFieldVariationSettings = + std::dynamic_pointer_cast< PeriodicGravityFieldVariationEstimatableParameterSettings >( vectorParameterName ); + if( gravityFieldVariationSettings == nullptr ) + { + throw std::runtime_error( "Error, expected periodic gravity field variation parameter settings " ); + } + + // Check consistency of body gravity field + std::shared_ptr< GravityFieldModel > gravityField = currentBody->getGravityFieldModel( ); + std::shared_ptr< TimeDependentSphericalHarmonicsGravityField > timeDepGravityField = + std::dynamic_pointer_cast< TimeDependentSphericalHarmonicsGravityField >( gravityField ); + if( timeDepGravityField == nullptr ) + { + throw std::runtime_error( + "Error, requested periodic gravity field variation parameter of " + + vectorParameterName->parameterType_.second.first + + ", but body does not have a time dependent spherical harmonic gravity field." ); + } + else if( currentBody->getGravityFieldVariationSet( ) == nullptr ) + { + throw std::runtime_error( "Error, requested periodic gravity field variation parameter of " + + vectorParameterName->parameterType_.second.first + + ", but body does not have gravity field variations" ); + } + else + { + + // Get associated gravity field variation + std::pair< bool, std::shared_ptr< gravitation::GravityFieldVariations > > gravityFieldVariation = + currentBody->getGravityFieldVariationSet( )->getGravityFieldVariation( periodic_variation ); + if( gravityFieldVariation.first == 0 ) + { + throw std::runtime_error( "Error when creating periodic gravity field variation parameter; associated gravity field model not found." ); + } + std::shared_ptr< gravitation::PeriodicGravityFieldVariations > periodicVariaton = + std::dynamic_pointer_cast< gravitation::PeriodicGravityFieldVariations >( + gravityFieldVariation.second ); + + // Create parameter object + if( periodicVariaton != nullptr ) + { + vectorParameterToEstimate = std::make_shared< PeriodicGravityFieldVariationsParameters >( + periodicVariaton, + gravityFieldVariationSettings->cosineBlockIndicesPerPower_, + gravityFieldVariationSettings->sineBlockIndicesPerPower_, + currentBodyName ); + } + else + { + throw std::runtime_error( + "Error, expected PeriodicGravityFieldVariations when creating periodic gravity field variation parameter" ); + } + } + break; + } case custom_estimated_parameter: { std::shared_ptr< CustomEstimatableParameterSettings > customParameterSettings = diff --git a/include/tudat/simulation/estimation_setup/estimatableParameterSettings.h b/include/tudat/simulation/estimation_setup/estimatableParameterSettings.h index 4012063ddb..e80676d503 100644 --- a/include/tudat/simulation/estimation_setup/estimatableParameterSettings.h +++ b/include/tudat/simulation/estimation_setup/estimatableParameterSettings.h @@ -1012,6 +1012,29 @@ class PolynomialGravityFieldVariationEstimatableParameterSettings: public Estima }; +class PeriodicGravityFieldVariationEstimatableParameterSettings: public EstimatableParameterSettings +{ +public: + + //! Constructor + /*! + * Constructor + * \param associatedBody Body being deformed + * \param deformingBody Body causing deformed + */ + PeriodicGravityFieldVariationEstimatableParameterSettings( + const std::string &associatedBody, + const std::map > >& cosineBlockIndicesPerPower, + const std::map > >& sineBlockIndicesPerPower ) : + EstimatableParameterSettings( associatedBody, periodic_gravity_field_variation_amplitudes ), + cosineBlockIndicesPerPower_( cosineBlockIndicesPerPower ), + sineBlockIndicesPerPower_( sineBlockIndicesPerPower ){ } + + std::map< int, std::vector< std::pair< int, int > > > cosineBlockIndicesPerPower_; + std::map< int, std::vector< std::pair< int, int > > > sineBlockIndicesPerPower_; + +}; + class CustomEstimatableParameterSettings: public EstimatableParameterSettings { diff --git a/src/astro/gravitation/periodicGravityFieldVariations.cpp b/src/astro/gravitation/periodicGravityFieldVariations.cpp index 24654fb8e4..02b8dfb854 100644 --- a/src/astro/gravitation/periodicGravityFieldVariations.cpp +++ b/src/astro/gravitation/periodicGravityFieldVariations.cpp @@ -18,34 +18,41 @@ namespace gravitation //! Class constructor PeriodicGravityFieldVariations::PeriodicGravityFieldVariations( - const std::vector< Eigen::MatrixXd >& cosineAmplitudes, - const std::vector< Eigen::MatrixXd >& sineAmplitudes, + const std::vector< Eigen::MatrixXd >& cosineShAmplitudesCosineTime, + const std::vector< Eigen::MatrixXd >& cosineShAmplitudesSineTime, + const std::vector< Eigen::MatrixXd >& sineShAmplitudesCosineTime, + const std::vector< Eigen::MatrixXd >& sineShAmplitudesSineTime, const std::vector< double >& frequencies, - const std::vector< double >& phases, const double referenceEpoch, const int minimumDegree, const int minimumOrder ): - GravityFieldVariations( minimumDegree, minimumOrder, minimumDegree + cosineAmplitudes.at( 0 ).rows( ) - 1, - minimumOrder + cosineAmplitudes.at( 0 ).cols( ) - 1 ), - cosineAmplitudes_( cosineAmplitudes ), - sineAmplitudes_( sineAmplitudes ), + GravityFieldVariations( minimumDegree, minimumOrder, minimumDegree + cosineShAmplitudesCosineTime.at( 0 ).rows( ) - 1, + minimumOrder + cosineShAmplitudesCosineTime.at( 0 ).cols( ) - 1 ), + cosineShAmplitudesCosineTime_( cosineShAmplitudesCosineTime ), + cosineShAmplitudesSineTime_( cosineShAmplitudesSineTime ), + sineShAmplitudesCosineTime_( sineShAmplitudesCosineTime ), + sineShAmplitudesSineTime_( sineShAmplitudesSineTime ), frequencies_( frequencies ), - phases_( phases ), referenceEpoch_( referenceEpoch ) { - if( cosineAmplitudes.size( ) != sineAmplitudes.size( ) ) + if( cosineShAmplitudesCosineTime_.size( ) != frequencies_.size( ) ) { - throw std::runtime_error( "Error when making periodic gravity field variations, amplitude input sizes inconsistent" ); + throw std::runtime_error( "Error when making periodic gravity field variations, frequency input sizes (C_lm * cos time) inconsistent" ); } - if( cosineAmplitudes.size( ) != frequencies.size( ) ) + if( cosineShAmplitudesSineTime_.size( ) != frequencies_.size( ) ) { - throw std::runtime_error( "Error when making periodic gravity field variations, frequency input size inconsistent" ); + throw std::runtime_error( "Error when making periodic gravity field variations, frequency input size (C_lm * sin time) inconsistent" ); } - if( cosineAmplitudes.size( ) != phases.size( ) ) + if( sineShAmplitudesCosineTime_.size( ) != frequencies_.size( ) ) { - throw std::runtime_error( "Error when making periodic gravity field variations, phase input size inconsistent" ); + throw std::runtime_error( "Error when making periodic gravity field variations, frequency input size (S_lm * cos time) inconsistent" ); + } + + if( sineShAmplitudesSineTime_.size( ) != frequencies_.size( ) ) + { + throw std::runtime_error( "Error when making periodic gravity field variations, frequency input size (S_lm * sin time) inconsistent" ); } } @@ -58,10 +65,10 @@ std::pair< Eigen::MatrixXd, Eigen::MatrixXd > PeriodicGravityFieldVariations::ca for( unsigned int i = 0; i < frequencies_.size( ); i++ ) { - cosineCorrections += cosineAmplitudes_.at( i ) * std::cos( - frequencies_.at( i ) * ( time - referenceEpoch_ ) + phases_.at( i ) ); - sineCorrections += sineAmplitudes_.at( i ) * std::cos( - frequencies_.at( i ) * ( time - referenceEpoch_ ) + phases_.at( i ) ); + cosineCorrections += cosineShAmplitudesCosineTime_.at( i ) * std::cos( frequencies_.at( i ) * ( time - referenceEpoch_ ) ) + + cosineShAmplitudesSineTime_.at( i ) * std::sin( frequencies_.at( i ) * ( time - referenceEpoch_ ) ) ; + sineCorrections += sineShAmplitudesCosineTime_.at( i ) * std::cos( frequencies_.at( i ) * ( time - referenceEpoch_ ) ) + + sineShAmplitudesSineTime_.at( i ) * std::sin( frequencies_.at( i ) * ( time - referenceEpoch_ ) ) ; } return std::make_pair( cosineCorrections, sineCorrections ); diff --git a/src/astro/orbit_determination/acceleration_partials/sphericalHarmonicAccelerationPartial.cpp b/src/astro/orbit_determination/acceleration_partials/sphericalHarmonicAccelerationPartial.cpp index 9b5d60b07c..6d410e342c 100644 --- a/src/astro/orbit_determination/acceleration_partials/sphericalHarmonicAccelerationPartial.cpp +++ b/src/astro/orbit_determination/acceleration_partials/sphericalHarmonicAccelerationPartial.cpp @@ -261,6 +261,28 @@ std::pair< std::function< void( Eigen::MatrixXd& ) >, int > SphericalHarmonicsGr std::placeholders::_1 ); numberOfRows = parameter->getParameterSize( ); + break; + } + case periodic_gravity_field_variation_amplitudes: + { + std::shared_ptr< PeriodicGravityFieldVariationsParameters > periodicVariationParameter = + std::dynamic_pointer_cast< PeriodicGravityFieldVariationsParameters >( parameter ); + std::map< std::pair< int, int >, std::vector< std::pair< int, int > > > indexAndPowerPerCosineBlockIndex = + periodicVariationParameter->getIndexAndPowerPerCosineBlockIndex( ); + std::map< std::pair< int, int >, std::vector< std::pair< int, int > > > indexAndPowerPerSineBlockIndex = + periodicVariationParameter->getIndexAndPowerPerSineBlockIndex( ); + + partialFunction = std::bind( &SphericalHarmonicsGravityPartial::wrtPeriodicGravityFieldVariations, this, + utilities::createVectorFromMapKeys( indexAndPowerPerCosineBlockIndex ), + utilities::createVectorFromMapKeys( indexAndPowerPerSineBlockIndex ), + utilities::createVectorFromMapValues( indexAndPowerPerCosineBlockIndex ), + utilities::createVectorFromMapValues( indexAndPowerPerSineBlockIndex ), + periodicVariationParameter->getPeriodicVariationModel( )->getFrequencies( ), + periodicVariationParameter->getPeriodicVariationModel( )->getReferenceEpoch( ), + std::placeholders::_1 ); + + numberOfRows = parameter->getParameterSize( ); + break; } default: break; @@ -477,9 +499,65 @@ void SphericalHarmonicsGravityPartial::wrtPolynomialGravityFieldVariations( staticSinePartialsMatrix.block( 0, i, 3, 1 ) * std::pow( ( currentTime_ - referenceEpoch), powersPerSineBlockIndex.at( i ).at( j ).second ); } } +} + + +void SphericalHarmonicsGravityPartial::wrtPeriodicGravityFieldVariations( + const std::vector< std::pair< int, int > >& cosineBlockIndices, + const std::vector< std::pair< int, int > >& sineBlockIndices, + const std::vector< std::vector< std::pair< int, int > > > powersPerCosineBlockIndex, + const std::vector< std::vector< std::pair< int, int > > > powersPerSineBlockIndex, + const std::vector< double >& frequencies, + const double referenceEpoch, + Eigen::MatrixXd& partialDerivatives ) +{ + Eigen::MatrixXd staticCosinePartialsMatrix = Eigen::MatrixXd::Zero( 3, cosineBlockIndices.size( ) ); + Eigen::MatrixXd staticSinePartialsMatrix = Eigen::MatrixXd::Zero( 3, sineBlockIndices.size( ) ); + + calculateSphericalHarmonicGravityWrtCCoefficients( + bodyFixedSphericalPosition_, bodyReferenceRadius_( ), gravitationalParameterFunction_( ), + sphericalHarmonicCache_, + cosineBlockIndices, coordinate_conversions::getSphericalToCartesianGradientMatrix( + bodyFixedPosition_ ), fromBodyFixedToIntegrationFrameRotation_( ), staticCosinePartialsMatrix, + maximumDegree_, maximumOrder_ ); + + calculateSphericalHarmonicGravityWrtSCoefficients( + bodyFixedSphericalPosition_, bodyReferenceRadius_( ), gravitationalParameterFunction_( ), + sphericalHarmonicCache_, + sineBlockIndices, coordinate_conversions::getSphericalToCartesianGradientMatrix( + bodyFixedPosition_ ), fromBodyFixedToIntegrationFrameRotation_( ), staticSinePartialsMatrix, + maximumDegree_, maximumOrder_ ); + + partialDerivatives.setZero( ); + + int counter = 0; + for( unsigned int i = 0; i < powersPerCosineBlockIndex.size( ); i++ ) + { + for( unsigned int j = 0; j < powersPerCosineBlockIndex.at( i ).size( ); j++ ) + { + double frequency = frequencies.at( powersPerCosineBlockIndex.at( i ).at( j ).second ); + partialDerivatives.block( 0, 2 * powersPerCosineBlockIndex.at( i ).at( j ).first, 3, 1 ) += + staticCosinePartialsMatrix.block( 0, i, 3, 1 ) * std::cos( frequency * ( currentTime_ - referenceEpoch ) ); + partialDerivatives.block( 0, 2 * powersPerCosineBlockIndex.at( i ).at( j ).first + 1, 3, 1 ) += + staticCosinePartialsMatrix.block( 0, i, 3, 1 ) * std::sin( frequency * ( currentTime_ - referenceEpoch ) ); + counter++; + } + } + for( unsigned int i = 0; i < powersPerSineBlockIndex.size( ); i++ ) + { + for( unsigned int j = 0; j < powersPerSineBlockIndex.at( i ).size( ); j++ ) + { + double frequency = frequencies.at( powersPerSineBlockIndex.at( i ).at( j ).second ); + partialDerivatives.block( 0, 2 * ( powersPerSineBlockIndex.at( i ).at( j ).first + counter ), 3, 1 ) += + staticSinePartialsMatrix.block( 0, i, 3, 1 ) * std::cos( frequency * ( currentTime_ - referenceEpoch ) ); + partialDerivatives.block( 0, 2 * ( powersPerSineBlockIndex.at( i ).at( j ).first + counter ) + 1, 3, 1 ) += + staticSinePartialsMatrix.block( 0, i, 3, 1 ) * std::sin( frequency * ( currentTime_ - referenceEpoch ) ); + } + } } + //! Function to calculate the partial of the acceleration wrt a set of cosine coefficients. void SphericalHarmonicsGravityPartial::wrtCosineCoefficientBlock( const std::vector< std::pair< int, int > >& blockIndices, diff --git a/src/astro/orbit_determination/estimatable_parameters/estimatableParameter.cpp b/src/astro/orbit_determination/estimatable_parameters/estimatableParameter.cpp index 25c4ba6035..a9b5c23d99 100644 --- a/src/astro/orbit_determination/estimatable_parameters/estimatableParameter.cpp +++ b/src/astro/orbit_determination/estimatable_parameters/estimatableParameter.cpp @@ -153,6 +153,9 @@ std::string getParameterTypeString( const EstimatebleParametersEnum parameterTyp case polynomial_gravity_field_variation_amplitudes: parameterDescription = " Polynomial gravity field variations "; break; + case periodic_gravity_field_variation_amplitudes: + parameterDescription = " Periodic gravity field variations "; + break; default: std::string errorMessage = "Error when getting parameter string, did not recognize parameter " + std::to_string( parameterType ); @@ -312,6 +315,9 @@ bool isDoubleParameter( const EstimatebleParametersEnum parameterType ) case polynomial_gravity_field_variation_amplitudes: isDoubleParameter = false; break; + case periodic_gravity_field_variation_amplitudes: + isDoubleParameter = false; + break; default: throw std::runtime_error( "Error, parameter type " + std::to_string( parameterType ) + " not found when getting parameter type" ); @@ -440,6 +446,9 @@ bool isParameterNonTidalGravityFieldVariationProperty( const EstimatebleParamete case polynomial_gravity_field_variation_amplitudes: flag = true; break; + case periodic_gravity_field_variation_amplitudes: + flag = true; + break; default: flag = false; break; diff --git a/src/simulation/environment_setup/createGravityFieldVariations.cpp b/src/simulation/environment_setup/createGravityFieldVariations.cpp index 9149a4d750..784d3eee52 100644 --- a/src/simulation/environment_setup/createGravityFieldVariations.cpp +++ b/src/simulation/environment_setup/createGravityFieldVariations.cpp @@ -268,21 +268,25 @@ std::shared_ptr< gravitation::GravityFieldVariations > createGravityFieldVariati } else { - if( periodicGravityFieldVariationSettings->getCosineAmplitudes( ).size( ) > 0 ) + if( periodicGravityFieldVariationSettings->getCosineShAmplitudesCosineTime( ).size( ) > 0 || + periodicGravityFieldVariationSettings->getCosineShAmplitudesSineTime( ).size( ) > 0 || + periodicGravityFieldVariationSettings->getSineShAmplitudesCosineTime( ).size( ) > 0 || + periodicGravityFieldVariationSettings->getSineShAmplitudesSineTime( ).size( ) > 0) { - // Create variation. - gravityFieldVariationModel = std::make_shared< PeriodicGravityFieldVariations > - ( periodicGravityFieldVariationSettings->getCosineAmplitudes( ), - periodicGravityFieldVariationSettings->getSineAmplitudes( ), - periodicGravityFieldVariationSettings->getFrequencies( ), - periodicGravityFieldVariationSettings->getPhases( ), - periodicGravityFieldVariationSettings->getReferenceEpoch( ), - periodicGravityFieldVariationSettings->getMinimumDegree( ), - periodicGravityFieldVariationSettings->getMinimumOrder( ) ); + // Create variation. + gravityFieldVariationModel = std::make_shared< PeriodicGravityFieldVariations > + ( periodicGravityFieldVariationSettings->getCosineShAmplitudesCosineTime(), + periodicGravityFieldVariationSettings->getCosineShAmplitudesSineTime( ), + periodicGravityFieldVariationSettings->getSineShAmplitudesCosineTime( ), + periodicGravityFieldVariationSettings->getSineShAmplitudesSineTime( ), + periodicGravityFieldVariationSettings->getFrequencies( ), + periodicGravityFieldVariationSettings->getReferenceEpoch( ), + periodicGravityFieldVariationSettings->getMinimumDegree( ), + periodicGravityFieldVariationSettings->getMinimumOrder( ) ); } else { - throw std::runtime_error( "Error when creating periodic gravity field variations; no cosine amplitudes found" ); + throw std::runtime_error( "Error when creating periodic gravity field variations; no amplitudes found" ); } } diff --git a/tests/src/astro/gravitation/unitTestGravityFieldVariations.cpp b/tests/src/astro/gravitation/unitTestGravityFieldVariations.cpp index ad471f819a..e0f887a61e 100644 --- a/tests/src/astro/gravitation/unitTestGravityFieldVariations.cpp +++ b/tests/src/astro/gravitation/unitTestGravityFieldVariations.cpp @@ -317,45 +317,66 @@ BOOST_AUTO_TEST_CASE( testGravityFieldVariations ) void getPeriodicGravityFieldVariationSettings( - std::vector< Eigen::MatrixXd >& cosineAmplitudes, - std::vector< Eigen::MatrixXd >& sineAmplitudes, + std::vector< Eigen::MatrixXd >& cosineShAmplitudesCosineTime, + std::vector< Eigen::MatrixXd >& cosineShAmplitudesSineTime, + std::vector< Eigen::MatrixXd >& sineShAmplitudesCosineTime, + std::vector< Eigen::MatrixXd >& sineShAmplitudesSineTime, std::vector< double >& frequencies, - std::vector< double >& phases, double& referenceEpoch, int& minimumDegree, int& minimumOrder, const int type = 0 ) { - cosineAmplitudes.clear(); - sineAmplitudes.clear(); + cosineShAmplitudesCosineTime.clear(); + cosineShAmplitudesSineTime.clear(); + sineShAmplitudesCosineTime.clear(); + sineShAmplitudesSineTime.clear(); frequencies.clear(); - phases.clear(); - cosineAmplitudes.push_back( + cosineShAmplitudesCosineTime.push_back( ( Eigen::MatrixXd( 2, 3 )<<2.0E-8, -2.0E-8, 4.2E-9, 4.0E-8, -3.7E-9, 2.9E-7 ).finished( ) ); - cosineAmplitudes.push_back( + cosineShAmplitudesCosineTime.push_back( ( Eigen::MatrixXd( 2, 3 )<<5.0E-8, -6.0E-8, 2.2E-9, 4.0E-7, -3.4E-7, 8.9E-7 ).finished( ) ); - cosineAmplitudes.push_back( + cosineShAmplitudesCosineTime.push_back( ( Eigen::MatrixXd( 2, 3 )<<8.0E-8, 2.0E-7, 7.2E-8, 7.5E-7, -1.9E-7, 2.5E-7 ).finished( ) ); - sineAmplitudes.push_back( + cosineShAmplitudesSineTime.push_back( + ( Eigen::MatrixXd( 2, 3 )<<3.0E-8, -6.0E-8, 2.2E-9, + 9.0E-8, -6.7E-9, 3.9E-7 ).finished( ) ); + cosineShAmplitudesSineTime.push_back( + ( Eigen::MatrixXd( 2, 3 )<<8.0E-8, -9.0E-8, 4.2E-9, + 2.0E-7, -6.4E-7, 4.9E-7 ).finished( ) ); + cosineShAmplitudesSineTime.push_back( + ( Eigen::MatrixXd( 2, 3 )<<3.0E-8, 6.0E-7, 8.2E-8, + 1.5E-7, 5.9E-7, 9.5E-7 ).finished( ) ); + + + sineShAmplitudesCosineTime.push_back( ( Eigen::MatrixXd( 2, 3 )<<0.0, 6.0E-8, 5.2E-9, 0.0, 7.4E-8, 5.6E-7 ).finished( ) ); - sineAmplitudes.push_back( + sineShAmplitudesCosineTime.push_back( ( Eigen::MatrixXd( 2, 3 )<<0.0, 5.5E-8, 9.2E-9, 0.0, 4.4E-8, 1.1E-7 ).finished( ) ); - sineAmplitudes.push_back( + sineShAmplitudesCosineTime.push_back( ( Eigen::MatrixXd( 2, 3 )<<0.0, 7.8E-7, 4.2E-8, 0.0, 8.5E-7, -2.5E-8 ).finished( ) ); + sineShAmplitudesSineTime.push_back( + ( Eigen::MatrixXd( 2, 3 )<<0.0, 5.0E-8, -4.2E-9, + 0.0, 4.4E-8, -1.6E-7 ).finished( ) ); + sineShAmplitudesSineTime.push_back( + ( Eigen::MatrixXd( 2, 3 )<<0.0, -4.5E-8, 6.2E-9, + 0.0, 6.4E-8, 2.1E-7 ).finished( ) ); + sineShAmplitudesSineTime.push_back( + ( Eigen::MatrixXd( 2, 3 )<<0.0, 5.8E-7, 6.2E-8, + 0.0, 2.5E-7, 5.5E-8 ).finished( ) ); + frequencies.resize( 3 ); - phases.resize( 3 ); frequencies = { 2.0E-5, 7.0E-5, 4.0E-6 }; - phases = { 1.0, 2.4, 5.0 }; referenceEpoch = 5.0E8; if( type == 0 ) @@ -392,16 +413,17 @@ BOOST_AUTO_TEST_CASE( testPeriodicGravityFieldVariations ) double testTime = 1.0E7; for( int k = 0; k < 4; k++ ) { - std::vector< Eigen::MatrixXd > cosineAmplitudes; - std::vector< Eigen::MatrixXd > sineAmplitudes; + std::vector< Eigen::MatrixXd > cosineShAmplitudesCosineTime; + std::vector< Eigen::MatrixXd > cosineShAmplitudesSineTime; + std::vector< Eigen::MatrixXd > sineShAmplitudesCosineTime; + std::vector< Eigen::MatrixXd > sineShAmplitudesSineTime; std::vector< double > frequencies; - std::vector< double > phases; double referenceEpoch; int minimumDegree; int minimumOrder; int testIndex = ( k < 3 ) ? k : 0; getPeriodicGravityFieldVariationSettings( - cosineAmplitudes, sineAmplitudes, frequencies, phases, + cosineShAmplitudesCosineTime, cosineShAmplitudesSineTime, sineShAmplitudesCosineTime, sineShAmplitudesSineTime, frequencies, referenceEpoch, minimumDegree, minimumOrder, testIndex ); std::shared_ptr< PeriodicGravityFieldVariationsSettings > variationSettings; @@ -409,16 +431,18 @@ BOOST_AUTO_TEST_CASE( testPeriodicGravityFieldVariations ) { variationSettings = std::dynamic_pointer_cast< PeriodicGravityFieldVariationsSettings >( periodicGravityFieldVariationsSettings( - cosineAmplitudes, sineAmplitudes, - frequencies, phases, referenceEpoch, + cosineShAmplitudesCosineTime, cosineShAmplitudesSineTime, + sineShAmplitudesCosineTime, sineShAmplitudesSineTime, + frequencies, referenceEpoch, minimumDegree, minimumOrder ) ); } else { variationSettings = std::dynamic_pointer_cast< PeriodicGravityFieldVariationsSettings >( - periodicGravityFieldVariationsSettings( - cosineAmplitudes.at( 0 ), sineAmplitudes.at( 0 ), - frequencies.at( 0 ), 0.0, referenceEpoch, + periodicGravityFieldVariationsSettingsSingleFrequency( + cosineShAmplitudesCosineTime.at( 0 ), cosineShAmplitudesSineTime.at( 0 ), + sineShAmplitudesCosineTime.at( 0 ), sineShAmplitudesSineTime.at( 0 ), + frequencies.at( 0 ), referenceEpoch, minimumDegree, minimumOrder ) ); } @@ -458,12 +482,14 @@ BOOST_AUTO_TEST_CASE( testPeriodicGravityFieldVariations ) Eigen::MatrixXd expectedCosineCoefficientsCorrections = Eigen::MatrixXd::Zero( 5, 5 ); Eigen::MatrixXd expectedSineCoefficientsCorrections = Eigen::MatrixXd::Zero( 5, 5 ); - for( unsigned int j = 0; j < cosineAmplitudes.size( ); j++ ) + for( unsigned int j = 0; j < cosineShAmplitudesCosineTime.size( ); j++ ) { - expectedCosineCoefficientsCorrections.block( minimumDegree, minimumOrder, 2, 3 ) += cosineAmplitudes.at( j ) * std::cos( - frequencies.at( j ) * ( testTime - referenceEpoch ) + phases.at( j ) ); - expectedSineCoefficientsCorrections.block( minimumDegree, minimumOrder, 2, 3 ) += sineAmplitudes.at( j ) * std::cos( - frequencies.at( j ) * ( testTime - referenceEpoch ) + phases.at( j ) ); + expectedCosineCoefficientsCorrections.block( minimumDegree, minimumOrder, 2, 3 ) += + cosineShAmplitudesCosineTime.at( j ) * std::cos( frequencies.at( j ) * ( testTime - referenceEpoch ) ) + + cosineShAmplitudesSineTime.at( j ) * std::sin( frequencies.at( j ) * ( testTime - referenceEpoch ) ); + expectedSineCoefficientsCorrections.block( minimumDegree, minimumOrder, 2, 3 ) += + sineShAmplitudesCosineTime.at( j ) * std::cos( frequencies.at( j ) * ( testTime - referenceEpoch ) ) + + sineShAmplitudesSineTime.at( j ) * std::sin( frequencies.at( j ) * ( testTime - referenceEpoch ) ); } // Compare corrections against expected variations. @@ -497,9 +523,9 @@ BOOST_AUTO_TEST_CASE( testPeriodicGravityFieldVariations ) Eigen::MatrixXd expectedSineCoefficientsCorrections = Eigen::MatrixXd::Zero( 5, 5 ); expectedCosineCoefficientsCorrections.block( minimumDegree, minimumOrder, 2, 3 ) = - cosineAmplitudes.at( 0 ); + cosineShAmplitudesCosineTime.at( 0 ); expectedSineCoefficientsCorrections.block( minimumDegree, minimumOrder, 2, 3 ) = - sineAmplitudes.at( 0 ); + sineShAmplitudesCosineTime.at( 0 ); // Compare corrections against expected variations. for( unsigned int i = 0; i < 5; i++ ) diff --git a/tests/src/astro/orbit_determination/acceleration_partials/unitTestSphericalHarmonicPartials.cpp b/tests/src/astro/orbit_determination/acceleration_partials/unitTestSphericalHarmonicPartials.cpp index ba5ac46c2e..98c506652d 100644 --- a/tests/src/astro/orbit_determination/acceleration_partials/unitTestSphericalHarmonicPartials.cpp +++ b/tests/src/astro/orbit_determination/acceleration_partials/unitTestSphericalHarmonicPartials.cpp @@ -437,6 +437,64 @@ std::vector< std::shared_ptr< GravityFieldVariationSettings > > getEarthGravityF std::make_shared< PolynomialGravityFieldVariationsSettings >( cosineAmplitudes, sineAmplitudes, -1.0E4, 2, 0 ); gravityFieldVariations.push_back( polynomialGravityFieldVariations ); + + + std::vector< Eigen::MatrixXd > cosineShAmplitudesCosineTime; + std::vector< Eigen::MatrixXd > cosineShAmplitudesSineTime; + std::vector< Eigen::MatrixXd > sineShAmplitudesCosineTime; + std::vector< Eigen::MatrixXd > sineShAmplitudesSineTime; + std::vector< double > frequencies; + + cosineShAmplitudesCosineTime.push_back( + ( Eigen::MatrixXd( 2, 3 )<<2.0E-8, -2.0E-8, 4.2E-9, + 4.0E-8, -3.7E-9, 2.9E-7 ).finished( ) ); + cosineShAmplitudesCosineTime.push_back( + ( Eigen::MatrixXd( 2, 3 )<<5.0E-8, -6.0E-8, 2.2E-9, + 4.0E-7, -3.4E-7, 8.9E-7 ).finished( ) ); + cosineShAmplitudesCosineTime.push_back( + ( Eigen::MatrixXd( 2, 3 )<<8.0E-8, 2.0E-7, 7.2E-8, + 7.5E-7, -1.9E-7, 2.5E-7 ).finished( ) ); + + cosineShAmplitudesSineTime.push_back( + ( Eigen::MatrixXd( 2, 3 )<<3.0E-8, -6.0E-8, 2.2E-9, + 9.0E-8, -6.7E-9, 3.9E-7 ).finished( ) ); + cosineShAmplitudesSineTime.push_back( + ( Eigen::MatrixXd( 2, 3 )<<8.0E-8, -9.0E-8, 4.2E-9, + 2.0E-7, -6.4E-7, 4.9E-7 ).finished( ) ); + cosineShAmplitudesSineTime.push_back( + ( Eigen::MatrixXd( 2, 3 )<<3.0E-8, 6.0E-7, 8.2E-8, + 1.5E-7, 5.9E-7, 9.5E-7 ).finished( ) ); + + sineShAmplitudesCosineTime.push_back( + ( Eigen::MatrixXd( 2, 3 )<<0.0, 6.0E-8, 5.2E-9, + 0.0, 7.4E-8, 5.6E-7 ).finished( ) ); + sineShAmplitudesCosineTime.push_back( + ( Eigen::MatrixXd( 2, 3 )<<0.0, 5.5E-8, 9.2E-9, + 0.0, 4.4E-8, 1.1E-7 ).finished( ) ); + sineShAmplitudesCosineTime.push_back( + ( Eigen::MatrixXd( 2, 3 )<<0.0, 7.8E-7, 4.2E-8, + 0.0, 8.5E-7, -2.5E-8 ).finished( ) ); + + sineShAmplitudesSineTime.push_back( + ( Eigen::MatrixXd( 2, 3 )<<0.0, 5.0E-8, -4.2E-9, + 0.0, 4.4E-8, -1.6E-7 ).finished( ) ); + sineShAmplitudesSineTime.push_back( + ( Eigen::MatrixXd( 2, 3 )<<0.0, -4.5E-8, 6.2E-9, + 0.0, 6.4E-8, 2.1E-7 ).finished( ) ); + sineShAmplitudesSineTime.push_back( + ( Eigen::MatrixXd( 2, 3 )<<0.0, 5.8E-7, 6.2E-8, + 0.0, 2.5E-7, 5.5E-8 ).finished( ) ); + + frequencies.resize( 3 ); + frequencies = { 2.0E-5, 7.0E-5, 4.0E-6 }; + + std::shared_ptr< GravityFieldVariationSettings > periodicGravityFieldVariations = + std::make_shared< PeriodicGravityFieldVariationsSettings >( + cosineShAmplitudesCosineTime, cosineShAmplitudesSineTime, sineShAmplitudesCosineTime, sineShAmplitudesSineTime, + frequencies, 1.0E3, 2, 0 ); + + gravityFieldVariations.push_back( periodicGravityFieldVariations ); + return gravityFieldVariations; } @@ -603,7 +661,15 @@ BOOST_AUTO_TEST_CASE( testSphericalHarmonicAccelerationPartial ) parameterNames.push_back( std::make_shared< PolynomialGravityFieldVariationEstimatableParameterSettings >( "Earth", cosineBlockIndicesPerPower, sineBlockIndicesPerPower ) ); + std::map > > cosineBlockIndicesPerPeriod; + cosineBlockIndicesPerPeriod[ 0 ].push_back( std::make_pair( 2, 1 ) ); + cosineBlockIndicesPerPeriod[ 2 ].push_back( std::make_pair( 2, 0 ) ); + cosineBlockIndicesPerPeriod[ 1 ].push_back( std::make_pair( 3, 2 ) ); + std::map > > sineBlockIndicesPerPeriod; + sineBlockIndicesPerPeriod[ 1 ].push_back( std::make_pair( 3, 1 ) ); + parameterNames.push_back( std::make_shared< PeriodicGravityFieldVariationEstimatableParameterSettings >( + "Earth", cosineBlockIndicesPerPeriod, sineBlockIndicesPerPeriod ) ); std::shared_ptr< estimatable_parameters::EstimatableParameterSet< double > > parameterSet = createParametersToEstimate( parameterNames, bodies ); @@ -819,6 +885,12 @@ BOOST_AUTO_TEST_CASE( testSphericalHarmonicAccelerationPartial ) Eigen::MatrixXd testPartialWrtPolynomialVariations = calculateAccelerationWrtParameterPartials( vectorParametersIterator->second, gravitationalAcceleration, Eigen::VectorXd::Constant( 6, 1.0E-6 ), sphericalHarmonicFieldUpdate ); + vectorParametersIterator++; + Eigen::MatrixXd partialWrtPeriodicVariations = accelerationPartial->wrtParameter( + vectorParametersIterator->second ); + Eigen::MatrixXd testPartialWrtPeriodicVariations = calculateAccelerationWrtParameterPartials( + vectorParametersIterator->second, gravitationalAcceleration, Eigen::VectorXd::Constant( vectorParametersIterator->second->getParameterSize( ), 1.0E-2 ), sphericalHarmonicFieldUpdate ); + TUDAT_CHECK_MATRIX_CLOSE_FRACTION( testPartialWrtVehiclePosition, partialWrtVehiclePosition, 1.0E-6 ); TUDAT_CHECK_MATRIX_CLOSE_FRACTION( testPartialWrtVehicleVelocity, partialWrtVehicleVelocity, 1.0E-6 ); TUDAT_CHECK_MATRIX_CLOSE_FRACTION( testPartialWrtEarthPosition, partialWrtEarthPosition, 1.0E-6 ); @@ -856,6 +928,7 @@ BOOST_AUTO_TEST_CASE( testSphericalHarmonicAccelerationPartial ) TUDAT_CHECK_MATRIX_CLOSE_FRACTION( partialWrtComplexDegreeThreeLoveNumberAtSeparateOrder, testPartialWrtComplexDegreeThreeLoveNumberAtSeparateOrder, 1.0E-6 ); TUDAT_CHECK_MATRIX_CLOSE_FRACTION( partialWrtPolynomialVariations, testPartialWrtPolynomialVariations, 1.0E-12 ); + TUDAT_CHECK_MATRIX_CLOSE_FRACTION( partialWrtPeriodicVariations, testPartialWrtPeriodicVariations, 1.0E-12 ); } From 54564fab704cc5f7a458c43ab722b6e45254452e Mon Sep 17 00:00:00 2001 From: Dominic Dirkx Date: Mon, 8 Jul 2024 09:32:04 +0200 Subject: [PATCH 22/25] Minor update to test --- .../unitTestSphericalHarmonicPartials.cpp | 25 ------------------- 1 file changed, 25 deletions(-) diff --git a/tests/src/astro/orbit_determination/acceleration_partials/unitTestSphericalHarmonicPartials.cpp b/tests/src/astro/orbit_determination/acceleration_partials/unitTestSphericalHarmonicPartials.cpp index 98c506652d..a6dd32a1af 100644 --- a/tests/src/astro/orbit_determination/acceleration_partials/unitTestSphericalHarmonicPartials.cpp +++ b/tests/src/astro/orbit_determination/acceleration_partials/unitTestSphericalHarmonicPartials.cpp @@ -646,8 +646,6 @@ BOOST_AUTO_TEST_CASE( testSphericalHarmonicAccelerationPartial ) "Earth", 3, "", false ) ); parameterNames.push_back( std::make_shared< SingleDegreeVariableTidalLoveNumberEstimatableParameterSettings >( "Earth", 3, std::vector< int >{ 0, 3 }, "", true ) ); - parameterNames.push_back( std::make_shared< SingleDegreeVariableTidalLoveNumberEstimatableParameterSettings >( - "Earth", 3, std::vector< int >{ 0, 3 }, "", true ) ); std::map > > cosineBlockIndicesPerPower; cosineBlockIndicesPerPower[ 1 ].push_back( std::make_pair( 2, 0 ) ); @@ -857,27 +855,6 @@ BOOST_AUTO_TEST_CASE( testSphericalHarmonicAccelerationPartial ) Eigen::MatrixXd testPartialWrtComplexDegreeThreeLoveNumberAtSeparateOrder = calculateAccelerationWrtParameterPartials( vectorParametersIterator->second, gravitationalAcceleration, Eigen::VectorXd::Constant( 4, 10.0 ), sphericalHarmonicFieldUpdate ); - - - - Eigen::VectorXd nominalTidalParameter = vectorParametersIterator->second->getParameterValue( ); - - vectorParametersIterator->second->setParameterValue( nominalTidalParameter + Eigen::VectorXd::Constant( nominalTidalParameter.rows( ), 1.0 ) ); - earthGravityField->update( testTime ); - Eigen::MatrixXd upperturbedCosineCoefficients = - earthGravityField->getCosineCoefficients( ).block( 0, 0, 3, 3 ); - Eigen::MatrixXd upperturbedSineCoefficients = - earthGravityField->getSineCoefficients( ).block( 0, 0, 3, 3 ); - - vectorParametersIterator->second->setParameterValue( nominalTidalParameter - Eigen::VectorXd::Constant( nominalTidalParameter.rows( ), 1.0 ) ); - earthGravityField->update( testTime ); - - Eigen::MatrixXd downperturbedCosineCoefficients = - earthGravityField->getCosineCoefficients( ).block( 0, 0, 3, 3 ); - Eigen::MatrixXd downperturbedSineCoefficients = - earthGravityField->getSineCoefficients( ).block( 0, 0, 3, 3 ); - - vectorParametersIterator++; vectorParametersIterator++; Eigen::MatrixXd partialWrtPolynomialVariations = accelerationPartial->wrtParameter( @@ -919,11 +896,9 @@ BOOST_AUTO_TEST_CASE( testSphericalHarmonicAccelerationPartial ) BOOST_CHECK_EQUAL( testPartialWrtCosineCoefficients.cols( ), 17 ); BOOST_CHECK_EQUAL( testPartialWrtSineCoefficients.cols( ), 13 ); - TUDAT_CHECK_MATRIX_CLOSE_FRACTION( partialWrtDegreeTwoLoveNumberAtSeparateOrders, testPartialWrtDegreeTwoOrderTwoLoveNumberAtSeparateOrders, 1.0E-6 ); TUDAT_CHECK_MATRIX_CLOSE_FRACTION( partialWrtComplexDegreeTwoLoveNumber, testPartialWrtComplexDegreeTwoLoveNumber, 1.0E-6 ); - TUDAT_CHECK_MATRIX_CLOSE_FRACTION( partialWrtDegreeThreeLoveNumber, testPartialWrtDegreeThreeLoveNumber, 1.0E-6 ); TUDAT_CHECK_MATRIX_CLOSE_FRACTION( partialWrtComplexDegreeThreeLoveNumberAtSeparateOrder, testPartialWrtComplexDegreeThreeLoveNumberAtSeparateOrder, 1.0E-6 ); From ca742d1d0095b78a94a8b4ceea90d423d5d8a159 Mon Sep 17 00:00:00 2001 From: RivaAlkahal Date: Tue, 23 Jul 2024 20:23:33 +0200 Subject: [PATCH 23/25] unitTestMars tested --- cmake_modules/TudatFindBoost.cmake | 5 +- .../aerodynamics/marsDtmAtmosphereModel.cpp | 20 +- .../unitTestMarsDtmAtmosphere_.cpp | 96 +++++++++ tests_riva/TestMarsDtm.cpp | 72 ++++--- tests_riva/unitTestMarsDtm.cpp | 195 ++++++++++++++++++ 5 files changed, 347 insertions(+), 41 deletions(-) create mode 100644 tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere_.cpp create mode 100644 tests_riva/unitTestMarsDtm.cpp diff --git a/cmake_modules/TudatFindBoost.cmake b/cmake_modules/TudatFindBoost.cmake index 4c83493a41..e44ef2d725 100644 --- a/cmake_modules/TudatFindBoost.cmake +++ b/cmake_modules/TudatFindBoost.cmake @@ -44,7 +44,10 @@ foreach(_TUDAT_BOOST_COMPONENT ${_TUDAT_REQUIRED_BOOST_LIBS}) if(Boost_USE_STATIC_LIBS) add_library(Boost::${_TUDAT_BOOST_COMPONENT} STATIC IMPORTED) else() - add_library(Boost::${_TUDAT_BOOST_COMPONENT} UNKNOWN IMPORTED ../tests/src/simulation/unitTestPolyhedron.cpp) + add_library(Boost::${_TUDAT_BOOST_COMPONENT} UNKNOWN IMPORTED ../tests/src/simulation/unitTestPolyhedron.cpp + ../tests_riva/unitTestMarsDtm.cpp + ../tests_riva/unitTestMarsDtm.cpp + ../tests_riva/unitTestMarsDtm.cpp) endif() set_target_properties(Boost::${_TUDAT_BOOST_COMPONENT} PROPERTIES INTERFACE_INCLUDE_DIRECTORIES "${Boost_INCLUDE_DIRS}") diff --git a/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp b/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp index bb4a7a829c..7a9cf30704 100644 --- a/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp +++ b/src/astro/aerodynamics/marsDtmAtmosphereModel.cpp @@ -21,9 +21,11 @@ namespace aerodynamics // Compute the hash key?? basic_astrodynamics::DateTime currentDateTime_ = basic_astrodynamics::getCalendarDateFromTime( time ); currentF107_ = f107Function_( time ); + //std::cout << "F107: " << currentF107_ << std::endl; currentDensity_ = getTotalDensity( - altitude, longitude, latitude, + altitude, latitude, longitude, currentDateTime_.getMinute(), currentDateTime_.getHour(), currentDateTime_.getDay(), currentDateTime_.getMonth(), currentDateTime_.getYear() ); + //std::cout << "current Density: " << currentDensity_ << std::endl; currentGeopotentialAltitude_ = computeGeopotentialAltitude( altitude ); currentTemperature_138 = computeCurrentTemperature( latitude, longitude, currentDateTime_.getMinute(), currentDateTime_.getHour(), currentDateTime_.getDay(), currentDateTime_.getMonth(), currentDateTime_.getYear(), 3 ); @@ -136,6 +138,7 @@ namespace aerodynamics alpha_( {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.38, -0.40, 0.0}), // thermal diffusion coefficients mmass_( {44.01, 16.00, 28.0, 40.0, 28.0, 32.0, 4.0, 1.0, 2.0} )// molar mass of the species { + //std::cout<<"MarsDtmAtmosphereModel constructor"<(outResults); Ls_ = std::get<0>(outResults); - //std::cout << "Ls: " << Ls << std::endl; + //std::cout << "Ls: " << Ls_ << std::endl; //Equation of Time double EoT = 2.861 * sin( unit_conversions::convertDegreesToRadians(2*Ls_)) - 0.071 * sin( unit_conversions::convertDegreesToRadians(4*Ls_)) @@ -251,8 +254,8 @@ namespace aerodynamics double doy = date2.marsDayofYear(date2); //day of year //std::cout << "doy: " << doy << std::endl; double t = computeLocalSolarTime(longitude, day_, month_, year_, hours_, minutes_); //seconds + //std::cout<<"F107: "< +#include "fstream" +#include "iostream" + +#include + +#include "tudat/basics/testMacros.h" +#include "tudat/astro/aerodynamics/customAerodynamicCoefficientInterface.h" +#include "tudat/astro/aerodynamics/aerodynamicAcceleration.h" +#include "tudat/astro/reference_frames/aerodynamicAngleCalculator.h" +#include "tudat/simulation/propagation_setup/dynamicsSimulator.h" +#include "tudat/interface/spice/spiceEphemeris.h" +#include "tudat/interface/spice/spiceRotationalEphemeris.h" +#include "tudat/io/basicInputOutput.h" +#include "tudat/simulation/environment_setup/body.h" +#include "tudat/simulation/estimation_setup/createNumericalSimulator.h" +#include "tudat/simulation/environment_setup/defaultBodies.h" + +using namespace tudat::aerodynamics; + + +namespace tudat { + namespace unit_tests { + using namespace tudat; + using namespace aerodynamics; + using namespace simulation_setup; + using namespace numerical_integrators; + using namespace simulation_setup; + using namespace basic_astrodynamics; + using namespace propagators; + using namespace basic_mathematics; + using namespace basic_astrodynamics; + + BOOST_AUTO_TEST_SUITE(test_mars_dtm_atmosphere) + + BOOST_AUTO_TEST_CASE(testMarsDtmAtmosphere) + { + // Define tolerance for equality + double tolerance = 1.0E-15; + std::string filename = "/Users/ralkahal/OneDrive - Delft University of Technology/PhD/Programs/atmodensitydtm/dtm_mars"; + std::shared_ptr marsDtmAtmosphereSettings; + marsDtmAtmosphereSettings = std::make_shared( + filename, 3378.0E3); + std::shared_ptr marsAtmosphereModel = createAtmosphereModel( + marsDtmAtmosphereSettings, "Mars"); + std::shared_ptr atmosphereModel = + std::dynamic_pointer_cast( + createAtmosphereModel(marsDtmAtmosphereSettings, "Mars")); + int alt_km = 400E3; + int time = 86400; + double latitude = unit_conversions::convertDegreesToRadians(15.0); + double longitude = unit_conversions::convertDegreesToRadians(10.0); + double rho = atmosphereModel->getDensity(alt_km, longitude, latitude, time); + + // reference density for this specific inputs + double rho_ref = 2.57862799E-14; + + // Check density + BOOST_CHECK_CLOSE_FRACTION(rho, rho_ref, tolerance); + } + + BOOST_AUTO_TEST_SUITE_END() + + } // namespace unit_tests +}//namespace tudat +//{ +//namespace unit_tests +//{ +// +//BOOST_AUTO_TEST_SUITE( test_mars_dtm_atmosphere ) +// +//BOOST_AUTO_TEST_CASE( testMarsDtmAtmosphere ) +//{ +// +//} +// +//BOOST_AUTO_TEST_SUITE_END( ) +// +//} // namespace unit_tests +//} // namespace tudat diff --git a/tests_riva/TestMarsDtm.cpp b/tests_riva/TestMarsDtm.cpp index 34620a4872..939a2be6f8 100644 --- a/tests_riva/TestMarsDtm.cpp +++ b/tests_riva/TestMarsDtm.cpp @@ -60,22 +60,28 @@ int main( ) { // std::dynamic_pointer_cast< MarsDtmAtmosphereModel >( // createAtmosphereModel( marsDtmAtmosphereSettings, "Mars" ) ); std::string spaceWeatherFilePath = paths::getSpaceWeatherDataPath( ) + "/sw19571001.txt"; - tudat::input_output::solar_activity::SolarActivityDataMap solarActivityData = - tudat::input_output::solar_activity::readSolarActivityData( spaceWeatherFilePath ); - std::shared_ptr< input_output::solar_activity::SolarActivityContainer > solarActivityContainer = - std::make_shared< input_output::solar_activity::SolarActivityContainer >( solarActivityData ); - std::function< double( const double ) > f107Function = [=](const double time) - { - return solarActivityContainer->getSolarActivityData( time )->solarRadioFlux107Observed / 2.25; - }; + //tudat::input_output::solar_activity::SolarActivityDataMap solarActivityData = + // tudat::input_output::solar_activity::readSolarActivityData( spaceWeatherFilePath ); + //std::shared_ptr< input_output::solar_activity::SolarActivityContainer > solarActivityContainer = + // std::make_shared< input_output::solar_activity::SolarActivityContainer >( solarActivityData ); + //std::function< double( const double ) > f107Function = [=](const double time) + //{ + // return solarActivityContainer->getSolarActivityData( time )->solarRadioFlux107Observed / 2.25; + //}; marsDtmAtmosphereSettings = std::make_shared< MarsDtmAtmosphereSettings >( filename, 3378.0E3); + //marsDtmAtmosphereSettings->setSpaceWeatherFile(spaceWeatherFilePath); + std::cout<<"Atmosphere settings created"< marsAtmosphereModel = createAtmosphereModel(marsDtmAtmosphereSettings, "Mars"); + std::cout<<"Mars atmosphere model created"< atmosphereModel = std::dynamic_pointer_cast< MarsDtmAtmosphereModel >( createAtmosphereModel( marsDtmAtmosphereSettings, "Mars" ) ); - + std::cout<<"Atmosphere model created"<(altitude); int alt_km = 400E3; //668 * 24 * 60 * 60 - for (int time = 0.0; time <= 345600 ; time += 86400) { + //for (int time = 0.0; time <= 345600 ; time += 86400) { //for (int time = 0; time <= 365 * 24 * 60 * 60; time += 60) { + //for (int time = 86400; time <= 86400 ; time += 86400) { + int time = 86400; + basic_astrodynamics::DateTime currentDateTime_ = basic_astrodynamics::getCalendarDateFromTime(time); - basic_astrodynamics::DateTime currentDateTime_ = basic_astrodynamics::getCalendarDateFromTime(time); - int minutes = currentDateTime_.getMinute(); - int hours = currentDateTime_.getHour(); - int days = currentDateTime_.getDay(); - int months = currentDateTime_.getMonth(); - int years = currentDateTime_.getYear(); + int minutes = currentDateTime_.getMinute(); + int hours = currentDateTime_.getHour(); + int days = currentDateTime_.getDay(); + int months = currentDateTime_.getMonth(); + int years = currentDateTime_.getYear(); - std::string filename = - "/Users/ralkahal/OneDrive - Delft University of Technology/PhD/atmodensitypds/densityFiles/" + - std::to_string(time) + ".txt"; - std::ofstream outfile(filename); + std::string filename1 = "/Users/ralkahal/OneDrive - Delft University of Technology/PhD/Programs/atmodensitydtm/" + std::to_string(time) + ".txt"; + std::ofstream outfile(filename1); //std::cout << "Time: " << time << " minutes: " << minutes << " hours: " << hours << " days: " << days << " months: " << months << " years: " << years << std::endl; //double latitude = 0.0; //double longitude = 0.0; if (outfile.is_open()) { - for (double latitude = -90.0; latitude <= 90.0; latitude += 1.0) { - for (double longitude = 0.0; longitude <= 360.0; longitude += 1.0) { - double rho = atmosphereModel->getTotalDensity(alt_km, unit_conversions::convertDegreesToRadians(latitude), unit_conversions::convertDegreesToRadians(longitude), minutes, hours, days, + std::cout<<"File opened successfully"<getDensity(alt_km, unit_conversions::convertDegreesToRadians(longitude), unit_conversions::convertDegreesToRadians(latitude), time); + std::cout<<"Density: "<getTotalDensity(alt_km, unit_conversions::convertDegreesToRadians(latitude), unit_conversions::convertDegreesToRadians(longitude), minutes, hours, days, months, years); - auto results = computeSolarLongitude(longitude, days, months, years, hours, minutes); - double Ls = std::get<0>(results); - double currentF107 = 0.0; // Assuming this is defined somewhere - outfile << latitude << " " << longitude << " " << rho << " " << Ls << "\n"; - } - } + std::cout<<"Total Density: "<(results); + //double currentF107 = 0.0; // Assuming this is defined somewhere + outfile << latitude << " " << longitude << " " << rho << " " << Ls << "\n"; + // } + //} } else { std::cerr << "Error creating file " << filename << std::endl; @@ -157,7 +169,7 @@ int main( ) { //std::cout << "Density computation and output written to file successfully." << std::endl; //std::cout << atmosphereModel->getTotalDensity( 1000.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000) << std::endl; - } + // } return 0; //std::cout << atmosphereModel->getTotalDensity( 138.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000) < +#include "fstream" +#include "iostream" + +#include + +#include "tudat/basics/testMacros.h" +#include "tudat/astro/aerodynamics/customAerodynamicCoefficientInterface.h" +#include "tudat/astro/aerodynamics/aerodynamicAcceleration.h" +#include "tudat/astro/reference_frames/aerodynamicAngleCalculator.h" +#include "tudat/simulation/propagation_setup/dynamicsSimulator.h" +#include "tudat/interface/spice/spiceEphemeris.h" +#include "tudat/interface/spice/spiceRotationalEphemeris.h" +#include "tudat/io/basicInputOutput.h" +#include "tudat/simulation/environment_setup/body.h" +#include "tudat/simulation/estimation_setup/createNumericalSimulator.h" +#include "tudat/simulation/environment_setup/defaultBodies.h" +#include "tudat/simulation/environment_setup/body.h" +#include "tudat/astro/basic_astro/timeConversions.h" +#include "tudat/astro/aerodynamics/marsDtmAtmosphereModel.h" +#include "tudat/io/solarActivityData.h" +#include "tudat/astro/basic_astro/unitConversions.h" + + +int main( ) { + using namespace tudat; + using namespace aerodynamics; + using namespace simulation_setup; + using namespace numerical_integrators; + using namespace simulation_setup; + using namespace basic_astrodynamics; + using namespace propagators; + using namespace basic_mathematics; + using namespace basic_astrodynamics; + + + //std::string filename = "/Users/ralkahal/Documents/PhD/atmodensitypds/dtm_mars"; + std::string filename = "/Users/ralkahal/OneDrive - Delft University of Technology/PhD/Programs/atmodensitydtm/dtm_mars"; + + //marsDtmAtmosphereSettings = + //std::shared_ptr< AtmosphereSettings > atmosphereSettings; + std::shared_ptr< AtmosphereSettings > marsDtmAtmosphereSettings; + //std::shared_ptr< AtmosphereModel > atmosphereModel = + // std::dynamic_pointer_cast< MarsDtmAtmosphereModel >( + // createAtmosphereModel( marsDtmAtmosphereSettings, "Mars" ) ); + std::string spaceWeatherFilePath = paths::getSpaceWeatherDataPath( ) + "/sw19571001.txt"; + //tudat::input_output::solar_activity::SolarActivityDataMap solarActivityData = + // tudat::input_output::solar_activity::readSolarActivityData( spaceWeatherFilePath ); + //std::shared_ptr< input_output::solar_activity::SolarActivityContainer > solarActivityContainer = + // std::make_shared< input_output::solar_activity::SolarActivityContainer >( solarActivityData ); + //std::function< double( const double ) > f107Function = [=](const double time) + //{ + // return solarActivityContainer->getSolarActivityData( time )->solarRadioFlux107Observed / 2.25; + //}; + + marsDtmAtmosphereSettings = std::make_shared< MarsDtmAtmosphereSettings >( + filename, 3378.0E3); + //marsDtmAtmosphereSettings->setSpaceWeatherFile(spaceWeatherFilePath); + std::cout<<"Atmosphere settings created"< marsAtmosphereModel = createAtmosphereModel(marsDtmAtmosphereSettings, "Mars"); + std::cout<<"Mars atmosphere model created"< atmosphereModel = + std::dynamic_pointer_cast< MarsDtmAtmosphereModel >( + createAtmosphereModel( marsDtmAtmosphereSettings, "Mars" ) ); + std::cout<<"Atmosphere model created"<computelocalsolartime(0.0,16 ,12,2000, 0.0, 1.0697333); + //std::cout<computeGl(0.0,0.0,1.0697333,0.0,16,12,2000,1)<computeGeopotentialAltitude( 255.0E3 )<computeCurrentTemperature( 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000, 1)<computeGamma( 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000, 1)<heightDistributionFunction(255.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000, 1)<(altitude); + + double rho = atmosphereModel->getTotalDensity(alt_km, 0.0, 0.0, 1.0697333, 0.0, 16, 12, 2000); + + // Write altitude and corresponding density to the file + outputFile << alt_km << " " << rho << "\n"; + } + */ + //for (int altitude = 138E3; altitude <= 1000E3; altitude += 10E3) { + // double alt_km = static_cast(altitude); + int alt_km = 400E3; + //668 * 24 * 60 * 60 + //for (int time = 0.0; time <= 345600 ; time += 86400) { + //for (int time = 0; time <= 365 * 24 * 60 * 60; time += 60) { + //for (int time = 86400; time <= 86400 ; time += 86400) { + int time = 86400; + basic_astrodynamics::DateTime currentDateTime_ = basic_astrodynamics::getCalendarDateFromTime(time); + + int minutes = currentDateTime_.getMinute(); + int hours = currentDateTime_.getHour(); + int days = currentDateTime_.getDay(); + int months = currentDateTime_.getMonth(); + int years = currentDateTime_.getYear(); + + std::string filename1 = "/Users/ralkahal/OneDrive - Delft University of Technology/PhD/Programs/atmodensitydtm/" + std::to_string(time) + ".txt"; + std::ofstream outfile(filename1); + + //std::cout << "Time: " << time << " minutes: " << minutes << " hours: " << hours << " days: " << days << " months: " << months << " years: " << years << std::endl; + //double latitude = 0.0; + //double longitude = 0.0; + if (outfile.is_open()) { + std::cout<<"File opened successfully"<getDensity(alt_km, unit_conversions::convertDegreesToRadians(longitude), unit_conversions::convertDegreesToRadians(latitude), time); + std::cout<<"Density: "<getTotalDensity(alt_km, unit_conversions::convertDegreesToRadians(latitude), unit_conversions::convertDegreesToRadians(longitude), minutes, hours, days, + months, years); + std::cout<<"Total Density: "<(results); + //double currentF107 = 0.0; // Assuming this is defined somewhere + outfile << latitude << " " << longitude << " " << rho << " " << Ls << "\n"; + // } + //} + + } else { + std::cerr << "Error creating file " << filename << std::endl; + } + outfile.close(); + + + + + //std::cout << "Density computation and output written to file successfully." << std::endl; + //std::cout << atmosphereModel->getTotalDensity( 1000.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000) << std::endl; + // } + return 0; + //std::cout << atmosphereModel->getTotalDensity( 138.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000) < Date: Wed, 24 Jul 2024 15:25:46 +0200 Subject: [PATCH 24/25] tests_riva removed --- cmake_modules/TudatFindBoost.cmake | 3 +- .../.cmake/api/v1/query/cache-v2 | 0 .../.cmake/api/v1/query/cmakeFiles-v1 | 0 .../.cmake/api/v1/query/codemodel-v2 | 0 .../.cmake/api/v1/query/toolchains-v1 | 0 tests_riva/CMakeLists.txt | 8 - tests_riva/TestMarsDtm.cpp | 195 ----------------- tests_riva/TestMarsDtm_.cpp | 204 ------------------ tests_riva/unitTestMarsDtm.cpp | 195 ----------------- 9 files changed, 2 insertions(+), 603 deletions(-) create mode 100644 tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/cache-v2 create mode 100644 tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/cmakeFiles-v1 create mode 100644 tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/codemodel-v2 create mode 100644 tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/toolchains-v1 delete mode 100644 tests_riva/CMakeLists.txt delete mode 100644 tests_riva/TestMarsDtm.cpp delete mode 100644 tests_riva/TestMarsDtm_.cpp delete mode 100644 tests_riva/unitTestMarsDtm.cpp diff --git a/cmake_modules/TudatFindBoost.cmake b/cmake_modules/TudatFindBoost.cmake index e44ef2d725..ae4ffd71fe 100644 --- a/cmake_modules/TudatFindBoost.cmake +++ b/cmake_modules/TudatFindBoost.cmake @@ -47,7 +47,8 @@ foreach(_TUDAT_BOOST_COMPONENT ${_TUDAT_REQUIRED_BOOST_LIBS}) add_library(Boost::${_TUDAT_BOOST_COMPONENT} UNKNOWN IMPORTED ../tests/src/simulation/unitTestPolyhedron.cpp ../tests_riva/unitTestMarsDtm.cpp ../tests_riva/unitTestMarsDtm.cpp - ../tests_riva/unitTestMarsDtm.cpp) + ../tests_riva/unitTestMarsDtm.cpp + ../tests/src/astro/aerodynamics/unitTestMarsDtmAtmosphere_.cpp) endif() set_target_properties(Boost::${_TUDAT_BOOST_COMPONENT} PROPERTIES INTERFACE_INCLUDE_DIRECTORIES "${Boost_INCLUDE_DIRS}") diff --git a/tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/cache-v2 b/tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/cache-v2 new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/cmakeFiles-v1 b/tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/cmakeFiles-v1 new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/codemodel-v2 b/tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/codemodel-v2 new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/toolchains-v1 b/tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/toolchains-v1 new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests_riva/CMakeLists.txt b/tests_riva/CMakeLists.txt deleted file mode 100644 index f85989c5c5..0000000000 --- a/tests_riva/CMakeLists.txt +++ /dev/null @@ -1,8 +0,0 @@ -TUDAT_ADD_EXECUTABLE(application_test_mars_dtm - "TestMarsDtm.cpp" - ${Tudat_ESTIMATION_LIBRARIES} -) -TUDAT_ADD_EXECUTABLE(application_test_mars_dtm2 - "TestMarsDtm_.cpp" - ${Tudat_ESTIMATION_LIBRARIES} -) diff --git a/tests_riva/TestMarsDtm.cpp b/tests_riva/TestMarsDtm.cpp deleted file mode 100644 index 939a2be6f8..0000000000 --- a/tests_riva/TestMarsDtm.cpp +++ /dev/null @@ -1,195 +0,0 @@ -/* Copyright (c) 2010-2019, 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 - * a copy of the license with this file. If not, please or visit: - * http://tudat.tudelft.nl/LICENSE. - * - * References - * - */ - -//#define BOOST_TEST_DYN_LINK -//#define BOOST_TEST_MAIN - -#include -#include "fstream" -#include "iostream" - -#include - -#include "tudat/basics/testMacros.h" -#include "tudat/astro/aerodynamics/customAerodynamicCoefficientInterface.h" -#include "tudat/astro/aerodynamics/aerodynamicAcceleration.h" -#include "tudat/astro/reference_frames/aerodynamicAngleCalculator.h" -#include "tudat/simulation/propagation_setup/dynamicsSimulator.h" -#include "tudat/interface/spice/spiceEphemeris.h" -#include "tudat/interface/spice/spiceRotationalEphemeris.h" -#include "tudat/io/basicInputOutput.h" -#include "tudat/simulation/environment_setup/body.h" -#include "tudat/simulation/estimation_setup/createNumericalSimulator.h" -#include "tudat/simulation/environment_setup/defaultBodies.h" -#include "tudat/simulation/environment_setup/body.h" -#include "tudat/astro/basic_astro/timeConversions.h" -#include "tudat/astro/aerodynamics/marsDtmAtmosphereModel.h" -#include "tudat/io/solarActivityData.h" -#include "tudat/astro/basic_astro/unitConversions.h" - - -int main( ) { - using namespace tudat; - using namespace aerodynamics; - using namespace simulation_setup; - using namespace numerical_integrators; - using namespace simulation_setup; - using namespace basic_astrodynamics; - using namespace propagators; - using namespace basic_mathematics; - using namespace basic_astrodynamics; - - - //std::string filename = "/Users/ralkahal/Documents/PhD/atmodensitypds/dtm_mars"; - std::string filename = "/Users/ralkahal/OneDrive - Delft University of Technology/PhD/Programs/atmodensitydtm/dtm_mars"; - - //marsDtmAtmosphereSettings = - //std::shared_ptr< AtmosphereSettings > atmosphereSettings; - std::shared_ptr< AtmosphereSettings > marsDtmAtmosphereSettings; - //std::shared_ptr< AtmosphereModel > atmosphereModel = - // std::dynamic_pointer_cast< MarsDtmAtmosphereModel >( - // createAtmosphereModel( marsDtmAtmosphereSettings, "Mars" ) ); - std::string spaceWeatherFilePath = paths::getSpaceWeatherDataPath( ) + "/sw19571001.txt"; - //tudat::input_output::solar_activity::SolarActivityDataMap solarActivityData = - // tudat::input_output::solar_activity::readSolarActivityData( spaceWeatherFilePath ); - //std::shared_ptr< input_output::solar_activity::SolarActivityContainer > solarActivityContainer = - // std::make_shared< input_output::solar_activity::SolarActivityContainer >( solarActivityData ); - //std::function< double( const double ) > f107Function = [=](const double time) - //{ - // return solarActivityContainer->getSolarActivityData( time )->solarRadioFlux107Observed / 2.25; - //}; - - marsDtmAtmosphereSettings = std::make_shared< MarsDtmAtmosphereSettings >( - filename, 3378.0E3); - //marsDtmAtmosphereSettings->setSpaceWeatherFile(spaceWeatherFilePath); - std::cout<<"Atmosphere settings created"< marsAtmosphereModel = createAtmosphereModel(marsDtmAtmosphereSettings, "Mars"); - std::cout<<"Mars atmosphere model created"< atmosphereModel = - std::dynamic_pointer_cast< MarsDtmAtmosphereModel >( - createAtmosphereModel( marsDtmAtmosphereSettings, "Mars" ) ); - std::cout<<"Atmosphere model created"<computelocalsolartime(0.0,16 ,12,2000, 0.0, 1.0697333); - //std::cout<computeGl(0.0,0.0,1.0697333,0.0,16,12,2000,1)<computeGeopotentialAltitude( 255.0E3 )<computeCurrentTemperature( 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000, 1)<computeGamma( 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000, 1)<heightDistributionFunction(255.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000, 1)<(altitude); - - double rho = atmosphereModel->getTotalDensity(alt_km, 0.0, 0.0, 1.0697333, 0.0, 16, 12, 2000); - - // Write altitude and corresponding density to the file - outputFile << alt_km << " " << rho << "\n"; - } - */ - //for (int altitude = 138E3; altitude <= 1000E3; altitude += 10E3) { - // double alt_km = static_cast(altitude); - int alt_km = 400E3; - //668 * 24 * 60 * 60 - //for (int time = 0.0; time <= 345600 ; time += 86400) { - //for (int time = 0; time <= 365 * 24 * 60 * 60; time += 60) { - //for (int time = 86400; time <= 86400 ; time += 86400) { - int time = 86400; - basic_astrodynamics::DateTime currentDateTime_ = basic_astrodynamics::getCalendarDateFromTime(time); - - int minutes = currentDateTime_.getMinute(); - int hours = currentDateTime_.getHour(); - int days = currentDateTime_.getDay(); - int months = currentDateTime_.getMonth(); - int years = currentDateTime_.getYear(); - - std::string filename1 = "/Users/ralkahal/OneDrive - Delft University of Technology/PhD/Programs/atmodensitydtm/" + std::to_string(time) + ".txt"; - std::ofstream outfile(filename1); - - //std::cout << "Time: " << time << " minutes: " << minutes << " hours: " << hours << " days: " << days << " months: " << months << " years: " << years << std::endl; - //double latitude = 0.0; - //double longitude = 0.0; - if (outfile.is_open()) { - std::cout<<"File opened successfully"<getDensity(alt_km, unit_conversions::convertDegreesToRadians(longitude), unit_conversions::convertDegreesToRadians(latitude), time); - std::cout<<"Density: "<getTotalDensity(alt_km, unit_conversions::convertDegreesToRadians(latitude), unit_conversions::convertDegreesToRadians(longitude), minutes, hours, days, - months, years); - std::cout<<"Total Density: "<(results); - //double currentF107 = 0.0; // Assuming this is defined somewhere - outfile << latitude << " " << longitude << " " << rho << " " << Ls << "\n"; - // } - //} - - } else { - std::cerr << "Error creating file " << filename << std::endl; - } - outfile.close(); - - - - - //std::cout << "Density computation and output written to file successfully." << std::endl; - //std::cout << atmosphereModel->getTotalDensity( 1000.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000) << std::endl; - // } - return 0; - //std::cout << atmosphereModel->getTotalDensity( 138.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000) < -#include "fstream" -#include "iostream" - -#include - -#include "tudat/basics/testMacros.h" -#include "tudat/astro/aerodynamics/customAerodynamicCoefficientInterface.h" -#include "tudat/astro/aerodynamics/aerodynamicAcceleration.h" -#include "tudat/astro/reference_frames/aerodynamicAngleCalculator.h" -#include "tudat/simulation/propagation_setup/dynamicsSimulator.h" -#include "tudat/interface/spice/spiceEphemeris.h" -#include "tudat/interface/spice/spiceRotationalEphemeris.h" -#include "tudat/io/basicInputOutput.h" -#include "tudat/simulation/environment_setup/body.h" -#include "tudat/simulation/estimation_setup/createNumericalSimulator.h" -#include "tudat/simulation/environment_setup/defaultBodies.h" -#include "tudat/astro/basic_astro/celestialBodyConstants.h" - -using namespace tudat::aerodynamics; - -int main( ) -{ - using namespace tudat; - using namespace aerodynamics; - using namespace simulation_setup; - using namespace numerical_integrators; - using namespace simulation_setup; - using namespace basic_astrodynamics; - using namespace propagators; - using namespace basic_mathematics; - using namespace basic_astrodynamics; - using namespace orbital_element_conversions; - using namespace unit_conversions; - - //Load spice kernels. - spice_interface::loadStandardSpiceKernels( ); - - // Create Earth object - BodyListSettings defaultBodySettings = - getDefaultBodySettings( { "Mars" } ); - std::string filename ="/Users/ralkahal/OneDrive - Delft University of Technology/PhD/Programs/atmodensitydtm/dtm_mars";; - defaultBodySettings.at( "Mars" )->atmosphereSettings = marsDtmAtmosphereSettings( filename, 3378.0E3); - SystemOfBodies bodies = createSystemOfBodies( defaultBodySettings ); - - // Create vehicle object. - double vehicleMass = 400; - bodies.createEmptyBody( "Vehicle" ); - bodies.at( "Vehicle" )->setConstantBodyMass( vehicleMass ); - - // Set aerodynamic coefficients. - std::shared_ptr aerodynamicCoefficientSettings = - std::make_shared( - 2.0, 4.0, Eigen::Vector3d::Zero( ), Eigen::Vector3d::UnitX( ), Eigen::Vector3d::Zero( ), - negative_aerodynamic_frame_coefficients, negative_aerodynamic_frame_coefficients ); - bodies.at( "Vehicle" )->setAerodynamicCoefficientInterface( - createAerodynamicCoefficientInterface( aerodynamicCoefficientSettings, "Vehicle", bodies )); - - // Define acceleration model settings. - SelectedAccelerationMap accelerationMap; - std::vector bodiesToPropagate; - std::vector centralBodies; - std::map > > accelerationsOfVehicle; - accelerationsOfVehicle[ "Mars" ].push_back( std::make_shared( point_mass_gravity )); - accelerationsOfVehicle[ "Mars" ].push_back( std::make_shared( aerodynamic )); - accelerationMap[ "Vehicle" ] = accelerationsOfVehicle; - bodiesToPropagate.push_back( "Vehicle" ); - centralBodies.push_back( "Mars" ); - - // Set initial state - //Eigen::Vector6d systemInitialState; - //systemInitialState << 3378.0E3, 3000E3, 4200E3, 0.0, 0.0, 0.0; - - Eigen::Vector6d vehicleInitialStateInKeplerianElements; - vehicleInitialStateInKeplerianElements( semiMajorAxisIndex ) = 3389.0E3+400.0E3; - vehicleInitialStateInKeplerianElements( eccentricityIndex ) = 0; - vehicleInitialStateInKeplerianElements( inclinationIndex ) = convertDegreesToRadians( 85.3 ); - vehicleInitialStateInKeplerianElements( argumentOfPeriapsisIndex ) = - convertDegreesToRadians( 235.7 ); - vehicleInitialStateInKeplerianElements( longitudeOfAscendingNodeIndex ) = - convertDegreesToRadians( 23.4 ); - vehicleInitialStateInKeplerianElements( trueAnomalyIndex ) = convertDegreesToRadians( 139.87 ); - - double marsGravitationalParameter = celestial_body_constants::MARS_GRAVITATIONAL_PARAMETER; - - Eigen::VectorXd systemInitialState = convertKeplerianToCartesianElements( - vehicleInitialStateInKeplerianElements, - marsGravitationalParameter ); - - // Create acceleration models and propagation settings. - basic_astrodynamics::AccelerationMap accelerationModelMap = createAccelerationModelsMap( - bodies, accelerationMap, bodiesToPropagate, centralBodies ); - - // Set variables to save - std::vector > dependentVariables; - dependentVariables.push_back( - std::make_shared( - altitude_dependent_variable, "Vehicle", "Mars" )); - dependentVariables.push_back( - std::make_shared( - "Vehicle", reference_frames::longitude_angle )); - dependentVariables.push_back( - std::make_shared( - "Vehicle", reference_frames::latitude_angle )); - dependentVariables.push_back( - std::make_shared( - local_density_dependent_variable, "Vehicle", "Mars" )); - dependentVariables.push_back( - std::make_shared( - solar_longitude,"Mars")); - dependentVariables.push_back( - std::make_shared( - aerodynamic, "Vehicle", "Mars", 0 )); - - - // Set propagation/integration settings - std::shared_ptr terminationSettings = - std::make_shared( 365*24*3600.0 ); - std::shared_ptr > integratorSettings = - std::make_shared > - ( rungeKutta4, 0.0, 60.0 ); - std::shared_ptr > - translationalPropagatorSettings = - std::make_shared > - ( centralBodies, accelerationModelMap, bodiesToPropagate, systemInitialState, 0.0, - integratorSettings, terminationSettings, - cowell, dependentVariables ); - - std::cout << "Propagation started" << std::endl; - // Create simulation object and propagate dynamics. - SingleArcDynamicsSimulator<> dynamicsSimulator( - bodies, translationalPropagatorSettings ); - - std::map< double, Eigen::VectorXd > dependentVariableHistory = dynamicsSimulator.getDependentVariableHistory( ); - //std::map< double, Eigen::VectorXd > stateHistory = dynamicsSimulator.getEquationsOfMotionNumericalSolution( ); - std::map< double, Eigen::Matrix< double, Eigen::Dynamic, 1 > > dependentVariableOutput = - dynamicsSimulator.getDependentVariableHistory( ); - - std::vector times; - std::vector densityValues; - std::vector solarLongValues; - std::cout << "Saving densities started" << std::endl; - // Iterate over the data - std::for_each(dependentVariableOutput.begin(), dependentVariableOutput.end(), - [&](const auto& it) { - // Extract the density value - double density = it.second(3); // Assuming density is at index 4 - double solarLong = it.second(4); // Assuming density is at index 4 - - // Store time and density value - //convert time to days - times.push_back(it.first/86400.0); - - densityValues.push_back(density); - solarLongValues.push_back(solarLong); - - }); - std::ofstream outputFile( - "/Users/ralkahal/OneDrive - Delft University of Technology/PhD/Programs/atmodensitydtm/density_output_test.txt"); - if (!outputFile.is_open()) { - std::cerr << "Unable to open the file!" << std::endl; - return 1; // return an error code - } - for (int i = 0; i < times.size(); ++i) { - outputFile << times[i] << " " << densityValues[i] << " " << solarLongValues[i] << std::endl; - } - outputFile.close(); - std::cout << "Density computation and output written to file successfully." << std::endl; - - -} -// -//namespace tudat -//{ -//namespace unit_tests -//{ -// -//BOOST_AUTO_TEST_SUITE( test_mars_dtm_atmosphere ) -// -//BOOST_AUTO_TEST_CASE( testMarsDtmAtmosphere ) -//{ -// -//} -// -//BOOST_AUTO_TEST_SUITE_END( ) -// -//} // namespace unit_tests -//} // namespace tudat diff --git a/tests_riva/unitTestMarsDtm.cpp b/tests_riva/unitTestMarsDtm.cpp deleted file mode 100644 index 939a2be6f8..0000000000 --- a/tests_riva/unitTestMarsDtm.cpp +++ /dev/null @@ -1,195 +0,0 @@ -/* Copyright (c) 2010-2019, 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 - * a copy of the license with this file. If not, please or visit: - * http://tudat.tudelft.nl/LICENSE. - * - * References - * - */ - -//#define BOOST_TEST_DYN_LINK -//#define BOOST_TEST_MAIN - -#include -#include "fstream" -#include "iostream" - -#include - -#include "tudat/basics/testMacros.h" -#include "tudat/astro/aerodynamics/customAerodynamicCoefficientInterface.h" -#include "tudat/astro/aerodynamics/aerodynamicAcceleration.h" -#include "tudat/astro/reference_frames/aerodynamicAngleCalculator.h" -#include "tudat/simulation/propagation_setup/dynamicsSimulator.h" -#include "tudat/interface/spice/spiceEphemeris.h" -#include "tudat/interface/spice/spiceRotationalEphemeris.h" -#include "tudat/io/basicInputOutput.h" -#include "tudat/simulation/environment_setup/body.h" -#include "tudat/simulation/estimation_setup/createNumericalSimulator.h" -#include "tudat/simulation/environment_setup/defaultBodies.h" -#include "tudat/simulation/environment_setup/body.h" -#include "tudat/astro/basic_astro/timeConversions.h" -#include "tudat/astro/aerodynamics/marsDtmAtmosphereModel.h" -#include "tudat/io/solarActivityData.h" -#include "tudat/astro/basic_astro/unitConversions.h" - - -int main( ) { - using namespace tudat; - using namespace aerodynamics; - using namespace simulation_setup; - using namespace numerical_integrators; - using namespace simulation_setup; - using namespace basic_astrodynamics; - using namespace propagators; - using namespace basic_mathematics; - using namespace basic_astrodynamics; - - - //std::string filename = "/Users/ralkahal/Documents/PhD/atmodensitypds/dtm_mars"; - std::string filename = "/Users/ralkahal/OneDrive - Delft University of Technology/PhD/Programs/atmodensitydtm/dtm_mars"; - - //marsDtmAtmosphereSettings = - //std::shared_ptr< AtmosphereSettings > atmosphereSettings; - std::shared_ptr< AtmosphereSettings > marsDtmAtmosphereSettings; - //std::shared_ptr< AtmosphereModel > atmosphereModel = - // std::dynamic_pointer_cast< MarsDtmAtmosphereModel >( - // createAtmosphereModel( marsDtmAtmosphereSettings, "Mars" ) ); - std::string spaceWeatherFilePath = paths::getSpaceWeatherDataPath( ) + "/sw19571001.txt"; - //tudat::input_output::solar_activity::SolarActivityDataMap solarActivityData = - // tudat::input_output::solar_activity::readSolarActivityData( spaceWeatherFilePath ); - //std::shared_ptr< input_output::solar_activity::SolarActivityContainer > solarActivityContainer = - // std::make_shared< input_output::solar_activity::SolarActivityContainer >( solarActivityData ); - //std::function< double( const double ) > f107Function = [=](const double time) - //{ - // return solarActivityContainer->getSolarActivityData( time )->solarRadioFlux107Observed / 2.25; - //}; - - marsDtmAtmosphereSettings = std::make_shared< MarsDtmAtmosphereSettings >( - filename, 3378.0E3); - //marsDtmAtmosphereSettings->setSpaceWeatherFile(spaceWeatherFilePath); - std::cout<<"Atmosphere settings created"< marsAtmosphereModel = createAtmosphereModel(marsDtmAtmosphereSettings, "Mars"); - std::cout<<"Mars atmosphere model created"< atmosphereModel = - std::dynamic_pointer_cast< MarsDtmAtmosphereModel >( - createAtmosphereModel( marsDtmAtmosphereSettings, "Mars" ) ); - std::cout<<"Atmosphere model created"<computelocalsolartime(0.0,16 ,12,2000, 0.0, 1.0697333); - //std::cout<computeGl(0.0,0.0,1.0697333,0.0,16,12,2000,1)<computeGeopotentialAltitude( 255.0E3 )<computeCurrentTemperature( 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000, 1)<computeGamma( 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000, 1)<heightDistributionFunction(255.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000, 1)<(altitude); - - double rho = atmosphereModel->getTotalDensity(alt_km, 0.0, 0.0, 1.0697333, 0.0, 16, 12, 2000); - - // Write altitude and corresponding density to the file - outputFile << alt_km << " " << rho << "\n"; - } - */ - //for (int altitude = 138E3; altitude <= 1000E3; altitude += 10E3) { - // double alt_km = static_cast(altitude); - int alt_km = 400E3; - //668 * 24 * 60 * 60 - //for (int time = 0.0; time <= 345600 ; time += 86400) { - //for (int time = 0; time <= 365 * 24 * 60 * 60; time += 60) { - //for (int time = 86400; time <= 86400 ; time += 86400) { - int time = 86400; - basic_astrodynamics::DateTime currentDateTime_ = basic_astrodynamics::getCalendarDateFromTime(time); - - int minutes = currentDateTime_.getMinute(); - int hours = currentDateTime_.getHour(); - int days = currentDateTime_.getDay(); - int months = currentDateTime_.getMonth(); - int years = currentDateTime_.getYear(); - - std::string filename1 = "/Users/ralkahal/OneDrive - Delft University of Technology/PhD/Programs/atmodensitydtm/" + std::to_string(time) + ".txt"; - std::ofstream outfile(filename1); - - //std::cout << "Time: " << time << " minutes: " << minutes << " hours: " << hours << " days: " << days << " months: " << months << " years: " << years << std::endl; - //double latitude = 0.0; - //double longitude = 0.0; - if (outfile.is_open()) { - std::cout<<"File opened successfully"<getDensity(alt_km, unit_conversions::convertDegreesToRadians(longitude), unit_conversions::convertDegreesToRadians(latitude), time); - std::cout<<"Density: "<getTotalDensity(alt_km, unit_conversions::convertDegreesToRadians(latitude), unit_conversions::convertDegreesToRadians(longitude), minutes, hours, days, - months, years); - std::cout<<"Total Density: "<(results); - //double currentF107 = 0.0; // Assuming this is defined somewhere - outfile << latitude << " " << longitude << " " << rho << " " << Ls << "\n"; - // } - //} - - } else { - std::cerr << "Error creating file " << filename << std::endl; - } - outfile.close(); - - - - - //std::cout << "Density computation and output written to file successfully." << std::endl; - //std::cout << atmosphereModel->getTotalDensity( 1000.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000) << std::endl; - // } - return 0; - //std::cout << atmosphereModel->getTotalDensity( 138.0E3, 0.0, 0.0, 1.0697333, 0.0, 16 ,12, 2000) < Date: Wed, 24 Jul 2024 15:35:18 +0200 Subject: [PATCH 25/25] deleted additional cmake folder --- .../cmake-build-relwithdebinfo/.cmake/api/v1/query/cache-v2 | 0 .../cmake-build-relwithdebinfo/.cmake/api/v1/query/cmakeFiles-v1 | 0 .../cmake-build-relwithdebinfo/.cmake/api/v1/query/codemodel-v2 | 0 .../cmake-build-relwithdebinfo/.cmake/api/v1/query/toolchains-v1 | 0 4 files changed, 0 insertions(+), 0 deletions(-) delete mode 100644 tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/cache-v2 delete mode 100644 tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/cmakeFiles-v1 delete mode 100644 tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/codemodel-v2 delete mode 100644 tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/toolchains-v1 diff --git a/tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/cache-v2 b/tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/cache-v2 deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/cmakeFiles-v1 b/tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/cmakeFiles-v1 deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/codemodel-v2 b/tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/codemodel-v2 deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/toolchains-v1 b/tests/src/astro/aerodynamics/cmake-build-relwithdebinfo/.cmake/api/v1/query/toolchains-v1 deleted file mode 100644 index e69de29bb2..0000000000