diff --git a/examples/tudat/satellite_propagation/CMakeLists.txt b/examples/tudat/satellite_propagation/CMakeLists.txt index e5847b9992..7567301da2 100755 --- a/examples/tudat/satellite_propagation/CMakeLists.txt +++ b/examples/tudat/satellite_propagation/CMakeLists.txt @@ -13,20 +13,20 @@ # ${Tudat_PROPAGATION_LIBRARIES} # ) -TUDAT_ADD_EXECUTABLE(application_SinglePerturbedSatellitePropagator - "singlePerturbedSatellitePropagator.cpp" - ${Tudat_PROPAGATION_LIBRARIES} - ) +#TUDAT_ADD_EXECUTABLE(application_SinglePerturbedSatellitePropagator +# "singlePerturbedSatellitePropagator.cpp" +# ${Tudat_PROPAGATION_LIBRARIES} +# ) -TUDAT_ADD_EXECUTABLE(application_InnerSolarSystemPropagation - "innerSolarSystemPropagation.cpp" - ${Tudat_PROPAGATION_LIBRARIES} - ) +#TUDAT_ADD_EXECUTABLE(application_InnerSolarSystemPropagation +# "innerSolarSystemPropagation.cpp" +# ${Tudat_PROPAGATION_LIBRARIES} +# ) -TUDAT_ADD_EXECUTABLE(application_lageosRadiationPressureAcceleration - "lageosRadiationPressureAcceleration.cpp" - ${Tudat_PROPAGATION_LIBRARIES} - ) +#TUDAT_ADD_EXECUTABLE(application_lageosRadiationPressureAcceleration +# "lageosRadiationPressureAcceleration.cpp" +# ${Tudat_PROPAGATION_LIBRARIES} +# ) TUDAT_ADD_EXECUTABLE(application_EarthOrbiterBasicStateEstimation "earthOrbiterBasicStateEstimation.cpp" diff --git a/examples/tudat/satellite_propagation/earthOrbiterBasicStateEstimation.cpp b/examples/tudat/satellite_propagation/earthOrbiterBasicStateEstimation.cpp index a392bca57d..5394442483 100644 --- a/examples/tudat/satellite_propagation/earthOrbiterBasicStateEstimation.cpp +++ b/examples/tudat/satellite_propagation/earthOrbiterBasicStateEstimation.cpp @@ -9,7 +9,10 @@ */ #include - +#include "tudat/astro/ephemerides/constantRotationalEphemeris.h" +#include "tudat/simulation/environment_setup/createRadiationPressureInterface.h" +#include "tudat/astro/electromagnetism/radiationPressureAcceleration.h" +#include "tudat/astro/electromagnetism/radiationPressureTargetModel.h" #include //! Execute propagation of orbits of Vehicle and Obelix around the Earth. @@ -34,6 +37,8 @@ int main( ) using namespace tudat::coordinate_conversions; using namespace tudat::ground_stations; using namespace tudat::observation_models; + using namespace tudat::electromagnetism; + /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /////////////////////// CREATE ENVIRONMENT AND VEHICLE ////////////////////////////////////////////////////// @@ -46,12 +51,10 @@ int main( ) std::vector< std::string > bodyNames; bodyNames.push_back( "Earth" ); bodyNames.push_back( "Sun" ); - bodyNames.push_back( "Moon" ); - bodyNames.push_back( "Mars" ); // Specify initial time double initialEphemerisTime = double( 1.0E7 ); - double finalEphemerisTime = initialEphemerisTime + 3.0 * 86400.0; + double finalEphemerisTime = initialEphemerisTime + 0.5 * 86400.0; // Create bodies needed in simulation BodyListSettings bodySettings = getDefaultBodySettings( bodyNames, initialEphemerisTime - 3600.0,finalEphemerisTime + 3600.0 ); @@ -62,25 +65,36 @@ int main( ) bodies.createEmptyBody( "Vehicle" ); bodies.at( "Vehicle" )->setConstantBodyMass( 400.0 ); - // Create and add aerodynamic coefficient interface - double referenceArea = 4.0; - double aerodynamicDragCoefficient = 1.2; - std::shared_ptr< AerodynamicCoefficientSettings > aerodynamicCoefficientSettings = - constantAerodynamicCoefficientSettings( - referenceArea, ( Eigen::Vector3d( ) << aerodynamicDragCoefficient, -0.01, 0.1 ).finished( ), 1, 1 ); - addAerodynamicCoefficientInterface( - bodies, "Vehicle", aerodynamicCoefficientSettings ); - - // Create and add radiation pressure interace - double referenceAreaRadiation = 4.0; - double radiationPressureCoefficient = 1.2; - std::vector< std::string > occultingBodies = { "Earth" }; - std::shared_ptr< RadiationPressureInterfaceSettings > vehicleRadiationPressureSettings = - cannonBallRadiationPressureSettings( - "Sun", referenceAreaRadiation, radiationPressureCoefficient, occultingBodies ); - addRadiationPressureInterface( - bodies, "Vehicle", vehicleRadiationPressureSettings ); - addEmptyTabulateEphemeris( bodies, "Vehicle", "Earth" ); + // Define constant rotational ephemeris + //bodies.at( "Vehicle" )->setRotationalEphemeris( + // createRotationModel( + // orbitalStateBasedRotationSettings( "Sun", false, false, "ECLIPJ2000", "VehicleFixed" ), + // "Vehicle", bodies )); + Eigen::Vector7d rotationalStateVehicle; + rotationalStateVehicle.segment( 0, 4 ) = linear_algebra::convertQuaternionToVectorFormat( Eigen::Quaterniond( Eigen::Matrix3d::Identity() )); + rotationalStateVehicle.segment( 4, 3 ) = Eigen::Vector3d::Zero(); + bodies.at( "Vehicle" )->setRotationalEphemeris( + std::make_shared< ConstantRotationalEphemeris >( + rotationalStateVehicle, "ECLIPJ2000", "VehicleFixed" ) ); + //bodySettings.at( "Vehicle" )->rotationModelSettings = std::make_shared< SynchronousRotationModelSettings >( "Jupiter", "ECLIPJ2000", "VehicleFixed" ); + + std::vector < std::shared_ptr< system_models::VehicleExteriorPanel > > panels; + panels = { + std::make_shared< system_models::VehicleExteriorPanel >(9.9, Eigen::Vector3d::UnitX(), + reflectionLawFromSpecularAndDiffuseReflectivity(0.35, 0.20)), + std::make_shared< system_models::VehicleExteriorPanel >(9.9, Eigen::Vector3d::UnitY(), + reflectionLawFromSpecularAndDiffuseReflectivity(0.35, 0.25)), + }; + const std::string panelTypeId = "SolarPanel"; + panels.at(0)->setPanelTypeId(panelTypeId); + panels.at(1)->setPanelTypeId(panelTypeId); + + bodies.at( "Vehicle" )->setRadiationPressureTargetModels( + { std::make_shared(panels) } ); + + + //const auto bodyShape = std::make_shared< FullPanelledBodySettings > bodyWingPanelledGeometry(2,2,2,2,0.1,0.2,0.1,0.2); + /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /////////////////////// CREATE GROUND STATIONS ////////////////////////////////////////////////////// @@ -108,15 +122,11 @@ int main( ) std::map< std::string, std::vector< std::shared_ptr< AccelerationSettings > > > accelerationsOfVehicle; accelerationsOfVehicle[ "Earth" ] = { - sphericalHarmonicAcceleration( 32, 32 ), - aerodynamicAcceleration( ) }; + sphericalHarmonicAcceleration( 32, 32 )}; accelerationsOfVehicle[ "Sun" ] = { pointMassGravityAcceleration( ), - cannonBallRadiationPressureAcceleration( ) }; - accelerationsOfVehicle[ "Mars" ] = { - pointMassGravityAcceleration( ) }; - accelerationsOfVehicle[ "Moon" ] = { - pointMassGravityAcceleration( ) }; + radiationPressureAcceleration( ) }; + accelerationSettingsList[ "Vehicle" ] = accelerationsOfVehicle; // Set bodies for which initial state is to be estimated and integrated. @@ -161,25 +171,25 @@ int main( ) /////////////////////// DEFINE LINK ENDS FOR OBSERVATIONS ////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - // Create list of link ends in which station is receiver and in which station is transmitter (with spacecraft other link end). - std::vector< LinkEnds > stationReceiverLinkEnds; - std::vector< LinkEnds > stationTransmitterLinkEnds; + + // Define link ends. + std::vector< LinkDefinition > stationReceiverLinkEnds; + std::vector< LinkDefinition > stationTransmitterLinkEnds; for( unsigned int i = 0; i < groundStationNames.size( ); i++ ) { - LinkEnds linkEnds; - linkEnds[ transmitter ] = std::make_pair( "Earth", groundStationNames.at( i ) ); - linkEnds[ receiver ] = std::make_pair( "Vehicle", "" ); + LinkDefinition linkEnds; + linkEnds[ transmitter ] = std::pair< std::string, std::string >( std::make_pair( "Earth", groundStationNames.at( i ) ) ); + linkEnds[ receiver ] = std::make_pair< std::string, std::string >( "Vehicle", "" ); stationTransmitterLinkEnds.push_back( linkEnds ); - linkEnds.clear( ); - linkEnds[ receiver ] = std::make_pair( "Earth", groundStationNames.at( i ) ); - linkEnds[ transmitter ] = std::make_pair( "Vehicle", "" ); + linkEnds[ receiver ] = std::pair< std::string, std::string >( std::make_pair( "Earth", groundStationNames.at( i ) ) ); + linkEnds[ transmitter ] = std::make_pair< std::string, std::string >( "Vehicle", "" ); stationReceiverLinkEnds.push_back( linkEnds ); } // Define (arbitrarily) link ends to be used for 1-way range, 1-way doppler and angular position observables - std::map< ObservableType, std::vector< LinkEnds > > linkEndsPerObservable; + std::map< ObservableType, std::vector< LinkDefinition > > linkEndsPerObservable; linkEndsPerObservable[ one_way_range ].push_back( stationReceiverLinkEnds[ 0 ] ); linkEndsPerObservable[ one_way_range ].push_back( stationTransmitterLinkEnds[ 0 ] ); linkEndsPerObservable[ one_way_range ].push_back( stationReceiverLinkEnds[ 1 ] ); @@ -187,27 +197,6 @@ int main( ) linkEndsPerObservable[ one_way_doppler ].push_back( stationReceiverLinkEnds[ 1 ] ); linkEndsPerObservable[ one_way_doppler ].push_back( stationTransmitterLinkEnds[ 2 ] ); - linkEndsPerObservable[ angular_position ].push_back( stationReceiverLinkEnds[ 2 ] ); - linkEndsPerObservable[ angular_position ].push_back( stationTransmitterLinkEnds[ 1 ] ); - - /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - /////////////////////// CREATE OBSERVATION SETTINGS //////////////////////////////////////////// - /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - - observation_models::ObservationSettingsMap observationSettingsMap; - for( std::map< ObservableType, std::vector< LinkEnds > >::iterator linkEndIterator = linkEndsPerObservable.begin( ); - linkEndIterator != linkEndsPerObservable.end( ); linkEndIterator++ ) - { - ObservableType currentObservable = linkEndIterator->first; - - std::vector< LinkEnds > currentLinkEndsList = linkEndIterator->second; - for( unsigned int i = 0; i < currentLinkEndsList.size( ); i++ ) - { - // Define settings for observable, no light-time corrections, and biases for selected 1-way range links - observationSettingsMap.insert( std::make_pair( currentLinkEndsList.at( i ), - std::make_shared< ObservationSettings >( currentObservable ) ) ); - } - } /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /////////////////////// DEFINE PARAMETERS THAT ARE TO BE ESTIMATED //////////////////////////////////////////// @@ -216,12 +205,9 @@ int main( ) // Define list of parameters to estimate. std::vector< std::shared_ptr< EstimatableParameterSettings > > parameterNames; parameterNames = getInitialStateParameterSettings< double >( propagatorSettings, bodies ); - parameterNames.push_back( std::make_shared< EstimatableParameterSettings >( "Vehicle", radiation_pressure_coefficient ) ); - parameterNames.push_back( std::make_shared< EstimatableParameterSettings >( "Vehicle", constant_drag_coefficient ) ); - parameterNames.push_back( std::make_shared< SphericalHarmonicEstimatableParameterSettings >( - 2, 0, 2, 2, "Earth", spherical_harmonics_cosine_coefficient_block ) ); - parameterNames.push_back( std::make_shared< SphericalHarmonicEstimatableParameterSettings >( - 2, 1, 2, 2, "Earth", spherical_harmonics_sine_coefficient_block ) ); + const std::string dummy = "SolarPanel"; + parameterNames.push_back(std::make_shared("Vehicle", specular_reflectivity, dummy)); + parameterNames.push_back(std::make_shared("Vehicle", diffuse_reflectivity, dummy)); // Create parameters std::shared_ptr< estimatable_parameters::EstimatableParameterSet< double > > parametersToEstimate = @@ -230,166 +216,111 @@ int main( ) // Print identifiers and indices of parameters to terminal. printEstimatableParameterEntries( parametersToEstimate ); - /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - /////////////////////// INITIALIZE ORBIT DETERMINATION OBJECT //////////////////////////////////////////// - /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - - // Create orbit determination object (propagate orbit, create observation models) - OrbitDeterminationManager< double, double > orbitDeterminationManager = - OrbitDeterminationManager< double, double >( - bodies, parametersToEstimate, observationSettingsMap, - integratorSettings, propagatorSettings ); + std::vector< std::shared_ptr< ObservationModelSettings > > observationSettingsList; + for( std::map< ObservableType, std::vector< LinkDefinition > >::iterator linkEndIterator = linkEndsPerObservable.begin( ); + linkEndIterator != linkEndsPerObservable.end( ); linkEndIterator++ ) + { + ObservableType currentObservable = linkEndIterator->first; - /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - /////////////////////// SIMULATE OBSERVATIONS //////////////////////////////////////////// - /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - // Define time of first observation - double observationTimeStart = initialEphemerisTime + 1000.0; + std::vector< LinkDefinition > currentLinkEndsList = linkEndIterator->second; + for( unsigned int i = 0; i < currentLinkEndsList.size( ); i++ ) + { + observationSettingsList.push_back( + std::make_shared< ObservationModelSettings >( + currentObservable, currentLinkEndsList.at( i ), std::shared_ptr< LightTimeCorrectionSettings >( ) ) ); + } + } - // Define time between two observations - double observationInterval = 20.0; + // Create orbit determination object. + OrbitDeterminationManager< double, double > orbitDeterminationManager = OrbitDeterminationManager< double, double >( + bodies, parametersToEstimate, observationSettingsList, integratorSettings, propagatorSettings ); - // Simulate observations for 3 days + // Compute list of observation times. std::vector< double > baseTimeList; - for( unsigned int i = 0; i < 3; i++ ) + double observationTime = initialEphemerisTime + 1000.0; + double observationInterval = 60.0; + while(observationTime < finalEphemerisTime - 1000.0) { - // Simulate 500 observations per day (observationInterval apart) - for( unsigned int j = 0; j < 500; j++ ) - { - baseTimeList.push_back( observationTimeStart + static_cast< double >( i ) * 86400.0 + - static_cast< double >( j ) * observationInterval ); - } + baseTimeList.push_back( observationTime); + observationTime += observationInterval; } - // Create measureement simulation input - std::map< ObservableType, std::map< LinkEnds, std::pair< std::vector< double >, LinkEndType > > > measurementSimulationInput; - for( std::map< ObservableType, std::vector< LinkEnds > >::iterator linkEndIterator = linkEndsPerObservable.begin( ); + std::vector< std::shared_ptr< ObservationSimulationSettings< double > > > measurementSimulationInput; + for( std::map< ObservableType, std::vector< LinkDefinition > >::iterator linkEndIterator = linkEndsPerObservable.begin( ); linkEndIterator != linkEndsPerObservable.end( ); linkEndIterator++ ) { - // Define observable type and link ends ObservableType currentObservable = linkEndIterator->first; - std::vector< LinkEnds > currentLinkEndsList = linkEndIterator->second; - - // Define observation times and reference link ends + std::vector< LinkDefinition > currentLinkEndsList = linkEndIterator->second; for( unsigned int i = 0; i < currentLinkEndsList.size( ); i++ ) { - measurementSimulationInput[ currentObservable ][ currentLinkEndsList.at( i ) ] = - std::make_pair( baseTimeList, receiver ); + measurementSimulationInput.push_back( + std::make_shared< TabulatedObservationSimulationSettings< > >( + currentObservable, currentLinkEndsList[ i ], baseTimeList, receiver ) ); } } + // Simulate observations. + std::shared_ptr< ObservationCollection< > > observationsAndTimes = simulateObservations< double, double >( + measurementSimulationInput, orbitDeterminationManager.getObservationSimulators( ), bodies ); + + // Set typedefs for POD input (observation types, observation link ends, observation values, associated times with // reference link ends. typedef Eigen::Matrix< double, Eigen::Dynamic, 1 > ObservationVectorType; typedef std::map< LinkEnds, std::pair< ObservationVectorType, std::pair< std::vector< double >, LinkEndType > > > SingleObservablePodInputType; - typedef std::map< ObservableType, SingleObservablePodInputType > PodInputDataType; - - // Simulate observations - PodInputDataType observationsAndTimes = simulateObservations< double, double >( - measurementSimulationInput, orbitDeterminationManager.getObservationSimulators( ) ); /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////// PERTURB PARAMETER VECTOR AND ESTIMATE PARAMETERS //////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - // Perturb parameter estimate + // Perturb parameter estimate. Eigen::Matrix< double, Eigen::Dynamic, 1 > initialParameterEstimate = parametersToEstimate->template getFullParameterValues< double >( ); Eigen::Matrix< double, Eigen::Dynamic, 1 > truthParameters = initialParameterEstimate; Eigen::Matrix< double, Eigen::Dynamic, 1 > parameterPerturbation = Eigen::Matrix< double, Eigen::Dynamic, 1 >::Zero( truthParameters.rows( ) ); - parameterPerturbation.segment( 0, 3 ) = Eigen::Vector3d::Constant( 10.0 ); - parameterPerturbation.segment( 3, 3 ) = Eigen::Vector3d::Constant( 1.0E-2 ); - parameterPerturbation( 6 ) = 0.01; - parameterPerturbation( 7 ) = 0.01; + + // Perturbe initial state estimate. + parameterPerturbation.segment( 0, 3 ) = Eigen::Vector3d::Constant( 1.0 ); + parameterPerturbation.segment( 3, 3 ) = Eigen::Vector3d::Constant( 1.E-3 ); + parameterPerturbation.segment( 6, 1 ) = Eigen::Vector1d::Constant( 0.05 ); + parameterPerturbation.segment( 7, 1 ) = Eigen::Vector1d::Constant( 0.05 ); + initialParameterEstimate += parameterPerturbation; + parametersToEstimate->resetParameterValues( initialParameterEstimate ); + // Define estimation input - std::shared_ptr< PodInput< double, double > > podInput = - std::make_shared< PodInput< double, double > >( - observationsAndTimes, initialParameterEstimate.rows( ), - Eigen::MatrixXd::Zero( truthParameters.rows( ), truthParameters.rows( ) ), - initialParameterEstimate - truthParameters ); - podInput->defineEstimationSettings( true, true, false, true ); - - // Define observation weights (constant per observable type) + std::shared_ptr< EstimationInput< double, double > > estimationInput = std::make_shared< EstimationInput< double, double > >( + observationsAndTimes ); + std::map< observation_models::ObservableType, double > weightPerObservable; weightPerObservable[ one_way_range ] = 1.0 / ( 1.0 * 1.0 ); - weightPerObservable[ angular_position ] = 1.0 / ( 1.0E-5 * 1.0E-5 ); - weightPerObservable[ one_way_doppler ] = 1.0 / ( 1.0E-11 * 1.0E-11 ); - podInput->setConstantPerObservableWeightsMatrix( weightPerObservable ); + weightPerObservable[ one_way_doppler ] = 1.0 / ( 1.0E-11 * 1.0E-11 * physical_constants::SPEED_OF_LIGHT * physical_constants::SPEED_OF_LIGHT ); + + estimationInput->setConstantPerObservableWeightsMatrix( weightPerObservable ); + estimationInput->defineEstimationSettings( true, true, true, true, true ); + estimationInput->setConvergenceChecker( + std::make_shared< EstimationConvergenceChecker >( 10 ) ); // Perform estimation - std::shared_ptr< PodOutput< double > > podOutput = orbitDeterminationManager.estimateParameters( - podInput, std::make_shared< EstimationConvergenceChecker >( 4 ) ); + std::shared_ptr< EstimationOutput< double > > estimationOutput = orbitDeterminationManager.estimateParameters( + estimationInput ); - /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - /////////////////////// PROVIDE OUTPUT TO CONSOLE AND FILES ////////////////////////////////////////// - /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + input_output::writeMatrixToFile( estimationOutput->getParameterHistoryMatrix( ), + "MichaelearthOrbitParameterHistory.dat", 16, + tudat_applications::getOutputPath( ) ); + + Eigen::VectorXd parameterEstimate = estimationOutput->parameterEstimate_; + std::cout <<"parameter estimate: "<< ( parameterEstimate ).transpose( ) << std::endl; + + Eigen::VectorXd estimationError = parameterEstimate - truthParameters; + std::cout <<"estimation error: "<< ( estimationError ).transpose( ) << std::endl; - std::string outputSubFolder = "EarthOrbiterStateEstimationExample/"; - - // Print true estimation error, limited mostly by numerical error - Eigen::VectorXd estimationError = podOutput->parameterEstimate_ - truthParameters; - std::cout << "True estimation error is: " << std::endl << ( estimationError ).transpose( ) << std::endl; - std::cout << "Formal estimation error is: " << std::endl << podOutput->getFormalErrorVector( ).transpose( ) << std::endl; - - std::map< double, Eigen::VectorXd > propagatedErrors; - - propagateFormalErrors( - propagatedErrors, podOutput->getUnnormalizedCovarianceMatrix( ), - orbitDeterminationManager.getStateTransitionAndSensitivityMatrixInterface( ), - 60.0, initialEphemerisTime + 3600.0, finalEphemerisTime - 3600.0 ); - input_output::writeDataMapToTextFile( propagatedErrors, - "earthOrbitBasicEstimationPropagatedErrors.dat", - tudat_applications::getOutputPath( ) + outputSubFolder ); - - input_output::writeMatrixToFile( podOutput->getUnnormalizedCovarianceMatrix( ), - "earthOrbitBasicEstimationCovariance.dat", 16, - tudat_applications::getOutputPath( ) + outputSubFolder ); - input_output::writeMatrixToFile( podOutput->normalizedInformationMatrix_, - "earthOrbitBasicEstimationInformationMatrix.dat", 16, - tudat_applications::getOutputPath( ) + outputSubFolder ); - input_output::writeMatrixToFile( podOutput->weightsMatrixDiagonal_, - "earthOrbitBasicEstimationWeightsDiagonal.dat", 16, - tudat_applications::getOutputPath( ) + outputSubFolder ); - input_output::writeMatrixToFile( podOutput->residuals_, - "earthOrbitBasicEstimationResiduals.dat", 16, - tudat_applications::getOutputPath( ) + outputSubFolder ); - input_output::writeMatrixToFile( podOutput->getCorrelationMatrix( ), - "earthOrbitBasicEstimationCorrelations.dat", 16, - tudat_applications::getOutputPath( ) + outputSubFolder ); - input_output::writeMatrixToFile( podOutput->getResidualHistoryMatrix( ), - "earthOrbitBasicResidualHistory.dat", 16, - tudat_applications::getOutputPath( ) + outputSubFolder ); - input_output::writeMatrixToFile( podOutput->getParameterHistoryMatrix( ), - "earthOrbitBasicParameterHistory.dat", 16, - tudat_applications::getOutputPath( ) + outputSubFolder ); - input_output::writeMatrixToFile( getConcatenatedMeasurementVector( podInput->getObservationsAndTimes( ) ), - "earthOrbitBasicObservationMeasurements.dat", 16, - tudat_applications::getOutputPath( ) + outputSubFolder ); - input_output::writeMatrixToFile( utilities::convertStlVectorToEigenVector( - getConcatenatedTimeVector( podInput->getObservationsAndTimes( ) ) ), - "earthOrbitBasicObservationTimes.dat", 16, - tudat_applications::getOutputPath( ) + outputSubFolder ); - input_output::writeMatrixToFile( utilities::convertStlVectorToEigenVector( - getConcatenatedGroundStationIndex( podInput->getObservationsAndTimes( ) ).first ), - "earthOrbitBasicObservationLinkEnds.dat", 16, - tudat_applications::getOutputPath( ) + outputSubFolder ); - input_output::writeMatrixToFile( utilities::convertStlVectorToEigenVector( - getConcatenatedObservableTypes( podInput->getObservationsAndTimes( ) ) ), - "earthOrbitBasicObservationObservableTypes.dat", 16, - tudat_applications::getOutputPath( ) + outputSubFolder ); - input_output::writeMatrixToFile( estimationError, - "earthOrbitBasicObservationTrueEstimationError.dat", 16, - tudat_applications::getOutputPath( ) + outputSubFolder ); - input_output::writeMatrixToFile( podOutput->getFormalErrorVector( ), - "earthOrbitBasicObservationFormalEstimationError.dat", 16, - tudat_applications::getOutputPath( ) + outputSubFolder ); // Final statement. // The exit code EXIT_SUCCESS indicates that the program was successfully executed. return EXIT_SUCCESS; -} +} \ No newline at end of file diff --git a/include/tudat/astro/electromagnetism/radiationPressureTargetModel.h b/include/tudat/astro/electromagnetism/radiationPressureTargetModel.h index 1621e7646b..870aae130a 100644 --- a/include/tudat/astro/electromagnetism/radiationPressureTargetModel.h +++ b/include/tudat/astro/electromagnetism/radiationPressureTargetModel.h @@ -295,6 +295,14 @@ class PaneledRadiationPressureTargetModel : public RadiationPressureTargetModel const bool resetForces, const std::string sourceName = "" ) override; + Eigen::Vector3d evaluateRadiationPressureForcePartialWrtDiffuseReflectivity( + double sourceIrradiance, + const Eigen::Vector3d &sourceToTargetDirectionLocalFrame); + + Eigen::Vector3d evaluateRadiationPressureForcePartialWrtSpecularReflectivity( + double sourceIrradiance, + const Eigen::Vector3d &sourceToTargetDirectionLocalFrame); + std::vector< std::shared_ptr< system_models::VehicleExteriorPanel > >& getBodyFixedPanels( ) { @@ -347,6 +355,152 @@ class PaneledRadiationPressureTargetModel : public RadiationPressureTargetModel void saveLocalComputations( const std::string sourceName, const bool saveCosines ) override ; + std::vector< std::shared_ptr< system_models::VehicleExteriorPanel > > getPanelsFromId( + const std::string& panelTypeId) + { + std::vector< std::shared_ptr< system_models::VehicleExteriorPanel > > panels; + for( unsigned int i = 0; i < static_cast(totalNumberOfPanels_); i++ ) + { + if(fullPanels_.at(i)->getPanelTypeId() == panelTypeId) + { + panels.push_back(fullPanels_.at(i)); + } + } + return panels; + } + + double getAverageDiffuseReflectivity(const std::string& panelTypeId) + { + std::vector coefficients; + std::vector< std::shared_ptr< system_models::VehicleExteriorPanel > > panelsFromId; + panelsFromId = this->getPanelsFromId(panelTypeId); + for( unsigned int i = 0; i < panelsFromId.size(); i++ ){ + { + auto reflectionLaw = std::dynamic_pointer_cast(panelsFromId.at(i)->getReflectionLaw()); + coefficients.push_back(reflectionLaw->getDiffuseReflectivity()); + } + } + + // Check if coefficients are consistent + if (!coefficients.empty()) { + double firstCoefficient = coefficients[0]; + for (size_t i = 1; i < coefficients.size(); ++i) { + if (coefficients[i] != firstCoefficient) { + std::cerr << "Warning: Diffuse reflectivity coefficients for panel group " + << panelTypeId << " are not consistent. Returning average value" << std::endl; + break; } } } + + // Calculate average coefficient + double averageCoefficient = 0.0; + if (!coefficients.empty()) + { + averageCoefficient = std::accumulate(coefficients.begin(), coefficients.end(), 0.0) / coefficients.size(); + } + + return averageCoefficient; + }; + + double getAverageSpecularReflectivity(const std::string& panelTypeId) + { + std::vector coefficients; + std::vector< std::shared_ptr< system_models::VehicleExteriorPanel > > panelsFromId; + panelsFromId = this->getPanelsFromId(panelTypeId); + for( unsigned int i = 0; i < panelsFromId.size(); i++ ){ + { + auto reflectionLaw = std::dynamic_pointer_cast(panelsFromId.at(i)->getReflectionLaw()); + coefficients.push_back(reflectionLaw->getSpecularReflectivity()); + } + } + + // Check if coefficients are consistent + if (!coefficients.empty()) { + double firstCoefficient = coefficients[0]; + for (size_t i = 1; i < coefficients.size(); ++i) { + if (coefficients[i] != firstCoefficient) { + std::cerr << "Warning: Specular reflectivity coefficients for panel group " + << panelTypeId << " are not consistent. Returning average value" << std::endl; + break; } } } + + // Calculate average coefficient + double averageCoefficient = 0.0; + if (!coefficients.empty()) + { + averageCoefficient = std::accumulate(coefficients.begin(), coefficients.end(), 0.0) / coefficients.size(); + } + + return averageCoefficient; + }; + + std::vector getSpecularReflectivityForPanelTypeId(const std::string& panelTypeId) + { + std::vector coefficients; + std::vector< std::shared_ptr< system_models::VehicleExteriorPanel > > panelsFromId; + + panelsFromId = this->getPanelsFromId(panelTypeId); + for (unsigned int i = 0; i < panelsFromId.size(); i++) { + auto reflectionLaw = std::dynamic_pointer_cast(panelsFromId.at(i)->getReflectionLaw()); + coefficients.push_back(reflectionLaw->getSpecularReflectivity()); + } + + return coefficients; + } + + std::vector getDiffuseReflectivityForPanelTypeId(const std::string& panelTypeId) + { + std::vector coefficients; + std::vector< std::shared_ptr< system_models::VehicleExteriorPanel > > panelsFromId; + + panelsFromId = this->getPanelsFromId(panelTypeId); + for (unsigned int i = 0; i < panelsFromId.size(); i++) { + auto reflectionLaw = std::dynamic_pointer_cast(panelsFromId.at(i)->getReflectionLaw()); + coefficients.push_back(reflectionLaw->getDiffuseReflectivity()); + } + + return coefficients; + } + + + void setGroupSpecularReflectivity(const std::string& panelTypeId, double specularReflectivity) + { + std::vector< std::shared_ptr< system_models::VehicleExteriorPanel > > panelsFromId; + panelsFromId = this->getPanelsFromId(panelTypeId); + for( unsigned int i = 0; i < panelsFromId.size(); i++ ) + { + if(panelsFromId.at(i)->getReflectionLaw() != nullptr ) + { + auto reflectionLaw = std::dynamic_pointer_cast(fullPanels_.at(i)->getReflectionLaw()); + reflectionLaw->setSpecularReflectivity(specularReflectivity); + } + } + }; + void setGroupDiffuseReflectivity(const std::string& panelTypeId, double diffuseReflectivity) + { + std::vector< std::shared_ptr< system_models::VehicleExteriorPanel > > panelsFromId; + panelsFromId = this->getPanelsFromId(panelTypeId); + for( unsigned int i = 0; i < panelsFromId.size(); i++ ) + { + if(panelsFromId.at(i)->getReflectionLaw() != nullptr ) + { + auto reflectionLaw = std::dynamic_pointer_cast(fullPanels_.at(i)->getReflectionLaw()); + reflectionLaw->setDiffuseReflectivity(diffuseReflectivity); + } + } + }; + + void setGroupAbsorptivity(const std::string& panelTypeId, double absorbtivity) + { + std::vector< std::shared_ptr< system_models::VehicleExteriorPanel > > panelsFromId; + panelsFromId = this->getPanelsFromId(panelTypeId); + for( unsigned int i = 0; i < panelsFromId.size(); i++ ) + { + if(panelsFromId.at(i)->getReflectionLaw() != nullptr ) + { + auto reflectionLaw = std::dynamic_pointer_cast(fullPanels_.at(i)->getReflectionLaw()); + reflectionLaw->setAbsorptivity(absorbtivity); + } + } + }; + private: void updateMembers_( double currentTime ) override; diff --git a/include/tudat/astro/electromagnetism/reflectionLaw.h b/include/tudat/astro/electromagnetism/reflectionLaw.h index f76dbe39c7..f5ac50c300 100644 --- a/include/tudat/astro/electromagnetism/reflectionLaw.h +++ b/include/tudat/astro/electromagnetism/reflectionLaw.h @@ -71,6 +71,14 @@ class ReflectionLaw const Eigen::Vector3d& surfaceNormal, const Eigen::Vector3d& incomingDirection) const = 0; + virtual Eigen::Vector3d evaluateReactionVectorPartialWrtSpecularReflectivity( + const Eigen::Vector3d& surfaceNormal, + const Eigen::Vector3d& incomingDirection) const = 0; + + virtual Eigen::Vector3d evaluateReactionVectorPartialWrtDiffuseReflectivity( + const Eigen::Vector3d& surfaceNormal, + const Eigen::Vector3d& incomingDirection) const = 0; + virtual Eigen::Matrix3d evaluateReactionVectorDerivativeWrtTargetPosition( const Eigen::Vector3d& surfaceNormal, const Eigen::Vector3d& incomingDirection, @@ -124,6 +132,14 @@ class SpecularDiffuseMixReflectionLaw : public ReflectionLaw const Eigen::Vector3d& surfaceNormal, const Eigen::Vector3d& incomingDirection) const override; + Eigen::Vector3d evaluateReactionVectorPartialWrtSpecularReflectivity( + const Eigen::Vector3d& surfaceNormal, + const Eigen::Vector3d& incomingDirection) const override; + + Eigen::Vector3d evaluateReactionVectorPartialWrtDiffuseReflectivity( + const Eigen::Vector3d& surfaceNormal, + const Eigen::Vector3d& incomingDirection) const override; + Eigen::Matrix3d evaluateReactionVectorDerivativeWrtTargetPosition( const Eigen::Vector3d& surfaceNormal, const Eigen::Vector3d& incomingDirection, @@ -137,16 +153,31 @@ class SpecularDiffuseMixReflectionLaw : public ReflectionLaw return absorptivity_; } + void setAbsorptivity( const double absorptivity ) + { + absorptivity_ = absorptivity; + } + double getSpecularReflectivity() const { return specularReflectivity_; } + void setSpecularReflectivity( const double specularReflectivity) + { + specularReflectivity_ = specularReflectivity; + } + double getDiffuseReflectivity() const { return diffuseReflectivity_; } + void setDiffuseReflectivity( const double diffuseReflectivity) + { + diffuseReflectivity_ = diffuseReflectivity; + } + bool isWithInstantaneousReradiation() const { return withInstantaneousReradiation_; diff --git a/include/tudat/astro/orbit_determination.h b/include/tudat/astro/orbit_determination.h index 016244816c..9c3facaa9e 100644 --- a/include/tudat/astro/orbit_determination.h +++ b/include/tudat/astro/orbit_determination.h @@ -34,6 +34,8 @@ #include "orbit_determination/estimatable_parameters/directTidalTimeLag.h" #include "orbit_determination/estimatable_parameters/inverseTidalQualityFactor.h" #include "orbit_determination/estimatable_parameters/empiricalAccelerationCoefficients.h" +#include "orbit_determination/estimatable_parameters/specularReflectivity.h" +#include "orbit_determination/estimatable_parameters/diffuseReflectivity.h" #include "orbit_determination/estimatable_parameters/equivalencePrincipleViolationParameter.h" #include "orbit_determination/estimatable_parameters/estimatableParameter.h" #include "orbit_determination/estimatable_parameters/freeCoreNutationRate.h" diff --git a/include/tudat/astro/orbit_determination/acceleration_partials/fullRadiationPressureAccelerationPartial.h b/include/tudat/astro/orbit_determination/acceleration_partials/fullRadiationPressureAccelerationPartial.h index 75fd009dcc..977873dbd6 100644 --- a/include/tudat/astro/orbit_determination/acceleration_partials/fullRadiationPressureAccelerationPartial.h +++ b/include/tudat/astro/orbit_determination/acceleration_partials/fullRadiationPressureAccelerationPartial.h @@ -185,6 +185,46 @@ class RadiationPressureAccelerationPartial: public AccelerationPartial parameterSize = 1; } } + else if( parameter->getParameterName( ).first == estimatable_parameters::specular_reflectivity && + parameter->getParameterName( ).second.first == acceleratedBody_) + { + if(std::dynamic_pointer_cast( + radiationPressureAcceleration_->getTargetModel( ) ) != nullptr){ + throw std::runtime_error( "Error when creating specular reflectivity partial, PaneledRadiationPressureTargetModel not specified" ); + } + if(parameter->getParameterName( ).second.second == ""){ + throw std::runtime_error( "Error when creating specular reflectivity partial, panel group name not specified" ); + } + else{ + partialFunction = std::bind( &RadiationPressureAccelerationPartial::wrtSpecularReflectivity, + this, + std::placeholders::_1, + std::dynamic_pointer_cast( + radiationPressureAcceleration_->getTargetModel( ) ), + parameter->getParameterName( ).second.second ); + parameterSize = 1; + }; + } + else if( parameter->getParameterName( ).first == estimatable_parameters::diffuse_reflectivity && + parameter->getParameterName( ).second.first == acceleratedBody_) + { + if(std::dynamic_pointer_cast( + radiationPressureAcceleration_->getTargetModel( ) ) != nullptr){ + throw std::runtime_error( "Error when creating diffuse reflectivity partial, PaneledRadiationPressureTargetModel not specified" ); + } + if(parameter->getParameterName( ).second.second == ""){ + throw std::runtime_error( "Error when creating diffuse reflectivity partial, panel group name not specified" ); + } + else{ + partialFunction = std::bind( &RadiationPressureAccelerationPartial::wrtDiffuseReflectivity, + this, + std::placeholders::_1, + std::dynamic_pointer_cast( + radiationPressureAcceleration_->getTargetModel( ) ), + parameter->getParameterName( ).second.second ); + parameterSize = 1; + }; + } // Check if parameter dependency exists. else if( parameter->getParameterName( ).second.first == acceleratedBody_ && parameter->getParameterName( ).second.second == acceleratingBody_ ) { @@ -249,6 +289,16 @@ class RadiationPressureAccelerationPartial: public AccelerationPartial void wrtRadiationPressureCoefficient( Eigen::MatrixXd& partial, std::shared_ptr< electromagnetism::CannonballRadiationPressureTargetModel > targetModel ); + void wrtSpecularReflectivity( + Eigen::MatrixXd& partial, + std::shared_ptr< electromagnetism::PaneledRadiationPressureTargetModel > targetModel, + const std::string& panelTypeId); + + void wrtDiffuseReflectivity( + Eigen::MatrixXd& partial, + std::shared_ptr< electromagnetism::PaneledRadiationPressureTargetModel > targetModel, + const std::string& panelTypeId); + std::shared_ptr< electromagnetism::PaneledSourceRadiationPressureAcceleration > radiationPressureAcceleration_; std::shared_ptr< estimatable_parameters::CustomSingleAccelerationPartialCalculatorSet > customAccelerationPartialSet_; diff --git a/include/tudat/astro/orbit_determination/acceleration_partials/panelledRadiationPressureAccelerationPartial.h b/include/tudat/astro/orbit_determination/acceleration_partials/panelledRadiationPressureAccelerationPartial.h index 95e9482c73..d2ad05456c 100755 --- a/include/tudat/astro/orbit_determination/acceleration_partials/panelledRadiationPressureAccelerationPartial.h +++ b/include/tudat/astro/orbit_determination/acceleration_partials/panelledRadiationPressureAccelerationPartial.h @@ -116,6 +116,10 @@ class PanelledRadiationPressurePartial: public AccelerationPartial } } + void wrtSpecularReflectivity( Eigen::MatrixXd& partial, const std::string& panelTypeId ); + + void wrtDiffuseReflectivity( Eigen::MatrixXd& partial, const std::string& panelTypeId ); + //! Function for updating partial w.r.t. the bodies' positions /*! * Function for updating common blocks of partial to current state. For the panelled radiation @@ -182,10 +186,28 @@ class PanelledRadiationPressurePartial: public AccelerationPartial numberOfRows = 1; break; + default: break; } } + switch( parameter->getParameterName( ).first ) + { + case estimatable_parameters::specular_reflectivity: + partialFunction = std::bind( &PanelledRadiationPressurePartial::wrtSpecularReflectivity, + this, std::placeholders::_1, parameter->getParameterName( ).second.second ); + numberOfRows = 1; + break; + + case estimatable_parameters::diffuse_reflectivity: + partialFunction = std::bind( &PanelledRadiationPressurePartial::wrtDiffuseReflectivity, + this, std::placeholders::_1, parameter->getParameterName( ).second.second ); + numberOfRows = 1; + break; + default: + break; + } + } return std::make_pair( partialFunction, numberOfRows ); } diff --git a/include/tudat/astro/orbit_determination/estimatable_parameters/diffuseReflectivity.h b/include/tudat/astro/orbit_determination/estimatable_parameters/diffuseReflectivity.h new file mode 100644 index 0000000000..c5fd15af86 --- /dev/null +++ b/include/tudat/astro/orbit_determination/estimatable_parameters/diffuseReflectivity.h @@ -0,0 +1,137 @@ +/* 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_DIFFUSEREFLECTIVITY_H +#define TUDAT_DIFFUSEREFLECTIVITY_H + + + + +#include + +#include + +#include "tudat/basics/utilities.h" +#include "tudat/astro/orbit_determination/estimatable_parameters/estimatableParameter.h" +#include "tudat/astro/electromagnetism/radiationPressureInterface.h" + +namespace tudat +{ + +namespace estimatable_parameters +{ + +//! Interface class for the estimation of a diffuse reflectivity coefficient per group of panels +class DiffuseReflectivity: public EstimatableParameter< double > +{ + public: + //! Constructor. + /*! + * Constructor DiffuseReflectivity + * \param radiationPressureInterface Object containing the radiation pressure coefficient to be estimated. + * \param associatedBody Name of body containing the radiationPressureInterface object + * \param panelTypeID Name of panel group for which to estimate the coefficient + */ + DiffuseReflectivity( + const std::shared_ptr< electromagnetism::PaneledRadiationPressureTargetModel > radiationPressureInterface, + const std::string& associatedBody, + const std::string panelTypeId): + EstimatableParameter< double >( diffuse_reflectivity, associatedBody, panelTypeId ), + radiationPressureInterface_( radiationPressureInterface ), + panelTypeId_( panelTypeId ) + { + // check if the panelTypeID exists among the panels + std::vector< std::shared_ptr< system_models::VehicleExteriorPanel > > panelsFromId = + radiationPressureInterface_->getPanelsFromId( panelTypeId_); + if(panelsFromId.size() < 1){ + throw std::runtime_error( "Error when creating estimated diffuse reflectivity for " + + panelTypeId + " of " + associatedBody + ", no corresponding panels defined" ); + } + + } + //! Destructor. + ~DiffuseReflectivity( ) { } + + double getParameterValue( ) + { + // Retrieve all diffuse reflectivity values for the panels corresponding to the given panelTypeId + std::vector diffuseReflectivities = radiationPressureInterface_->getDiffuseReflectivityForPanelTypeId(panelTypeId_); + + // Calculate the average diffuse reflectivity + double averageDiffuseReflectivity = std::accumulate(diffuseReflectivities.begin(), diffuseReflectivities.end(), 0.0) / diffuseReflectivities.size(); + + // Check if all values are the same + bool allValuesSame = std::all_of(diffuseReflectivities.begin(), diffuseReflectivities.end(), + [&](double value) { return value == diffuseReflectivities[0]; }); + + // If not all values are the same, print a warning and reset the values to the average + if (!allValuesSame) + { + std::cerr << "Warning: Diffuse reflectivity values for panel group " + << panelTypeId_ << " are not consistent. Resetting all to the average value." << std::endl; + + // Set all panels' diffuse reflectivity to the average value + radiationPressureInterface_->setGroupDiffuseReflectivity(panelTypeId_, averageDiffuseReflectivity); + + // Adjust the absorptivity + double averageSpecularReflectivity = radiationPressureInterface_->getAverageSpecularReflectivity(panelTypeId_); + double apsorptivity = 1 - averageDiffuseReflectivity - averageSpecularReflectivity; + radiationPressureInterface_->setGroupAbsorptivity(panelTypeId_, apsorptivity); + } + + // Return the average diffuse reflectivity in all cases + return averageDiffuseReflectivity; + } + + + void setParameterValue( double parameterValue ) + { + // Set the diffuse reflectivity, even if it's greater than 1 (non-physical case) + radiationPressureInterface_->setGroupDiffuseReflectivity(panelTypeId_, parameterValue); + + // Get the current specular reflectivity + double specularReflectivity = radiationPressureInterface_->getAverageSpecularReflectivity( panelTypeId_ ); + + // Compute the absorptivity such that the sum of specular + diffuse + absorptivity = 1 + double absorptivity = 1.0 - parameterValue - specularReflectivity; + + // Check if the absorptivity is negative, which indicates non-physical behavior + if ( absorptivity < 0.0 ) + { + // Print a warning with both the non-physical specular reflectivity and the resulting negative absorptivity + std::cerr << "Warning: Non-physical behavior detected for panel group " + << panelTypeId_ << ". Diffuse reflectivity = " << parameterValue + << ", resulting absorptivity = " << absorptivity << "." + << std::endl; + } + + // Set the absorptivity for the panel group (even if it's negative) + radiationPressureInterface_->setGroupAbsorptivity(panelTypeId_, absorptivity); + + } + + int getParameterSize( ){ return 1; } + + protected: + + private: + + std::shared_ptr< electromagnetism::PaneledRadiationPressureTargetModel > radiationPressureInterface_; + + std::string panelTypeId_; + +}; + + +} // namespace estimatable_parameters + +} // namespace tudat + +#endif // TUDAT_DIFFUSEREFLECTIVITY_H \ No newline at end of file diff --git a/include/tudat/astro/orbit_determination/estimatable_parameters/estimatableParameter.h b/include/tudat/astro/orbit_determination/estimatable_parameters/estimatableParameter.h index 4485c01fb9..32ff89b808 100644 --- a/include/tudat/astro/orbit_determination/estimatable_parameters/estimatableParameter.h +++ b/include/tudat/astro/orbit_determination/estimatable_parameters/estimatableParameter.h @@ -84,6 +84,8 @@ enum EstimatebleParametersEnum periodic_gravity_field_variation_amplitudes, source_direction_radiation_pressure_scaling_factor, source_perpendicular_direction_radiation_pressure_scaling_factor, + specular_reflectivity, + diffuse_reflectivity, mode_coupled_tidal_love_numbers }; diff --git a/include/tudat/astro/orbit_determination/estimatable_parameters/specularReflectivity.h b/include/tudat/astro/orbit_determination/estimatable_parameters/specularReflectivity.h new file mode 100644 index 0000000000..ae18dbda6d --- /dev/null +++ b/include/tudat/astro/orbit_determination/estimatable_parameters/specularReflectivity.h @@ -0,0 +1,135 @@ +/* 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_SPECULARREFLECTIVITY_H +#define TUDAT_SPECULARREFLECTIVITY_H + + + + +#include + +#include + +#include "tudat/basics/utilities.h" +#include "tudat/astro/orbit_determination/estimatable_parameters/estimatableParameter.h" +#include "tudat/astro/electromagnetism/radiationPressureInterface.h" + +namespace tudat +{ + +namespace estimatable_parameters +{ + +//! Interface class for the estimation of a specular reflectivity coefficient per group of panels +class SpecularReflectivity: public EstimatableParameter< double > +{ + public: + //! Constructor. + /*! + * ConstructorSpecularReflectivity + * \param radiationPressureInterface Object containing the radiation pressure coefficient to be estimated. + * \param associatedBody Name of body containing the radiationPressureInterface object + * \param panelTypeID Name of panel group for which to estimate the coefficient + */ + SpecularReflectivity( + const std::shared_ptr< electromagnetism::PaneledRadiationPressureTargetModel > radiationPressureInterface, + const std::string& associatedBody, + const std::string panelTypeId): + EstimatableParameter< double >( specular_reflectivity, associatedBody, panelTypeId ), + radiationPressureInterface_( radiationPressureInterface ), + panelTypeId_( panelTypeId ) + { + // check if the panelTypeID exists among the panels + std::vector< std::shared_ptr< system_models::VehicleExteriorPanel > > panelsFromId = + radiationPressureInterface_->getPanelsFromId( panelTypeId_); + if(panelsFromId.size() < 1){ + throw std::runtime_error( "Error when creating estimated specular reflectivity for " + + panelTypeId + " of " + associatedBody + ", no corresponding panels defined" ); + } + + } + //! Destructor. + ~SpecularReflectivity( ) { } + + double getParameterValue( ) + { + // Retrieve all specular reflectivity values for the panels corresponding to the given panelTypeId + std::vector specularReflectivities = radiationPressureInterface_->getSpecularReflectivityForPanelTypeId(panelTypeId_); + + // Calculate the average diffuse reflectivity + double averageSpecularReflectivity = std::accumulate(specularReflectivities.begin(), specularReflectivities.end(), 0.0) / specularReflectivities.size(); + + // Check if all values are the same + bool allValuesSame = std::all_of(specularReflectivities.begin(), specularReflectivities.end(), + [&](double value) { return value == specularReflectivities[0]; }); + + // If not all values are the same, print a warning and reset the values to the average + if (!allValuesSame) + { + std::cerr << "Warning: Specular reflectivity values for panel group " + << panelTypeId_ << " are not consistent. Resetting all to the average value." << std::endl; + + // Set all panels' specular reflectivity to the average value + radiationPressureInterface_->setGroupSpecularReflectivity(panelTypeId_, averageSpecularReflectivity); + + // Adjust the absorptivity + double averageDiffuseReflectivity = radiationPressureInterface_->getAverageDiffuseReflectivity(panelTypeId_); + double apsorptivity = 1 - averageDiffuseReflectivity - averageSpecularReflectivity; + radiationPressureInterface_->setGroupAbsorptivity(panelTypeId_, apsorptivity); + } + + // Return the average diffuse reflectivity in all cases + return averageSpecularReflectivity; + } + + void setParameterValue( double parameterValue ) + { + // Set the specular reflectivity, even if it's greater than 1 (non-physical case) + radiationPressureInterface_->setGroupSpecularReflectivity(panelTypeId_, parameterValue); + + // Get the current diffuse reflectivity + double diffuseReflectivity = radiationPressureInterface_->getAverageDiffuseReflectivity( panelTypeId_ ); + + // Compute the absorptivity such that the sum of specular + diffuse + absorptivity = 1 + double absorptivity = 1.0 - parameterValue - diffuseReflectivity; + + // Check if the absorptivity is negative, which indicates non-physical behavior + if ( absorptivity < 0.0 ) + { + // Print a warning with both the non-physical specular reflectivity and the resulting negative absorptivity + std::cerr << "Warning: Non-physical behavior detected for panel group " + << panelTypeId_ << ". Specular reflectivity = " << parameterValue + << ", resulting absorptivity = " << absorptivity << "." + << std::endl; + } + + // Set the absorptivity for the panel group (even if it's negative) + radiationPressureInterface_->setGroupAbsorptivity(panelTypeId_, absorptivity); + } + + int getParameterSize( ){ return 1; } + + protected: + + private: + + std::shared_ptr< electromagnetism::PaneledRadiationPressureTargetModel > radiationPressureInterface_; + + std::string panelTypeId_; + +}; + + +} // namespace estimatable_parameters + +} // namespace tudat + +#endif // TUDAT_SPECULARREFLECTIVITY_H \ No newline at end of file diff --git a/include/tudat/astro/system_models/vehicleExteriorPanels.h b/include/tudat/astro/system_models/vehicleExteriorPanels.h index e2a0f2a512..56c5f72a66 100644 --- a/include/tudat/astro/system_models/vehicleExteriorPanels.h +++ b/include/tudat/astro/system_models/vehicleExteriorPanels.h @@ -141,6 +141,16 @@ class VehicleExteriorPanel return trackedBody_; } + std::string getPanelTypeId( ) const + { + return panelTypeId_; + } + + void setPanelTypeId( const std::string panelTypeId ) + { + panelTypeId_ = panelTypeId; + } + protected: std::function< Eigen::Vector3d( ) > frameFixedSurfaceNormal_; @@ -155,6 +165,8 @@ class VehicleExteriorPanel std::shared_ptr< electromagnetism::ReflectionLaw > reflectionLaw_; + std::string panelTypeId_; + }; } // namespace system_models diff --git a/include/tudat/simulation/estimation_setup/createEstimatableParameters.h b/include/tudat/simulation/estimation_setup/createEstimatableParameters.h index 173b7efeb1..23da5f8f08 100644 --- a/include/tudat/simulation/estimation_setup/createEstimatableParameters.h +++ b/include/tudat/simulation/estimation_setup/createEstimatableParameters.h @@ -49,6 +49,8 @@ #include "tudat/simulation/estimation_setup/estimatableParameterSettings.h" #include "tudat/simulation/propagation_setup/dynamicsSimulator.h" #include "tudat/simulation/environment_setup/body.h" +#include "tudat/astro/orbit_determination/estimatable_parameters/specularReflectivity.h" +#include "tudat/astro/orbit_determination/estimatable_parameters/diffuseReflectivity.h" namespace tudat { @@ -1310,6 +1312,38 @@ std::shared_ptr< estimatable_parameters::EstimatableParameter< double > > create currentBodyName, doubleParameterName->parameterType_.second.second ); break; } + case specular_reflectivity: + { + if( std::dynamic_pointer_cast< electromagnetism::PaneledRadiationPressureTargetModel >( currentBody->getRadiationPressureTargetModel( ) ) == nullptr) + { + std::string errorMessage = "Error, no panelled radiation pressure target model found in body " + + currentBodyName + " when making specular reflectivity parameter."; + throw std::runtime_error( errorMessage ); + } + else + { + doubleParameterToEstimate = std::make_shared< SpecularReflectivity >( + std::dynamic_pointer_cast< electromagnetism::PaneledRadiationPressureTargetModel >( currentBody->getRadiationPressureTargetModel( ) ), + currentBodyName, doubleParameterName->parameterType_.second.second); + } + break; + } + case diffuse_reflectivity: + { + if( std::dynamic_pointer_cast< electromagnetism::PaneledRadiationPressureTargetModel >( currentBody->getRadiationPressureTargetModel( ) ) == nullptr) + { + std::string errorMessage = "Error, no panelled radiation pressure target model found in body " + + currentBodyName + " when making diffuse reflectivity parameter."; + throw std::runtime_error( errorMessage ); + } + else + { + doubleParameterToEstimate = std::make_shared< DiffuseReflectivity >( + std::dynamic_pointer_cast< electromagnetism::PaneledRadiationPressureTargetModel >( currentBody->getRadiationPressureTargetModel( ) ), + currentBodyName, doubleParameterName->parameterType_.second.second); + } + break; + } default: throw std::runtime_error( "Warning, this double parameter has not yet been implemented when making parameters" ); break; @@ -1788,6 +1822,7 @@ std::shared_ptr< estimatable_parameters::EstimatableParameter< Eigen::VectorXd > } break; } + case arc_wise_constant_drag_coefficient: { // Check input consistency diff --git a/include/tudat/simulation/estimation_setup/estimatableParameterSettings.h b/include/tudat/simulation/estimation_setup/estimatableParameterSettings.h index 3474394726..02d37136c6 100644 --- a/include/tudat/simulation/estimation_setup/estimatableParameterSettings.h +++ b/include/tudat/simulation/estimation_setup/estimatableParameterSettings.h @@ -1604,6 +1604,19 @@ inline std::shared_ptr< EstimatableParameterSettings > radiationPressureTargetPe return std::make_shared< EstimatableParameterSettings >( targetName, source_perpendicular_direction_radiation_pressure_scaling_factor, sourceName ); } +inline std::shared_ptr< EstimatableParameterSettings > diffuseReflectivity( + const std::string bodyName, + const std::string panel_group_id ) +{ + return std::make_shared< EstimatableParameterSettings >( bodyName, diffuse_reflectivity, panel_group_id ); +} + +inline std::shared_ptr< EstimatableParameterSettings > specularReflectivity( + const std::string bodyName, + const std::string panel_group_id ) +{ + return std::make_shared< EstimatableParameterSettings >( bodyName, specular_reflectivity, panel_group_id ); +} inline std::shared_ptr< EstimatableParameterSettings > polynomialGravityFieldVariationParameter( diff --git a/src/astro/electromagnetism/radiationPressureTargetModel.cpp b/src/astro/electromagnetism/radiationPressureTargetModel.cpp index 1091223cf5..74718da534 100644 --- a/src/astro/electromagnetism/radiationPressureTargetModel.cpp +++ b/src/astro/electromagnetism/radiationPressureTargetModel.cpp @@ -141,6 +141,95 @@ void PaneledRadiationPressureTargetModel::saveLocalComputations( const std::stri panelTorquesPerSource_[ sourceName ] = panelTorques_; } +Eigen::Vector3d PaneledRadiationPressureTargetModel::evaluateRadiationPressureForcePartialWrtDiffuseReflectivity( + double sourceIrradiance, + const Eigen::Vector3d& sourceToTargetDirectionLocalFrame) +{ + Eigen::Vector3d forcePartialWrtDiffuseReflectivity = Eigen::Vector3d::Zero(); + + double radiationPressure = sourceIrradiance / physical_constants::SPEED_OF_LIGHT; + auto segmentFixedPanelsIterator = segmentFixedPanels_.begin( ); + int counter = 0; + Eigen::Quaterniond currentOrientation; + + for( unsigned int i = 0; i < segmentFixedPanels_.size( ) + 1; i++ ) + { + currentOrientation = Eigen::Quaterniond( Eigen::Matrix3d::Identity( ) ); + if( i > 0 ) + { + currentOrientation = segmentFixedToBodyFixedRotations_.at( segmentFixedPanelsIterator->first )( ); + + } + + const std::vector< std::shared_ptr< system_models::VehicleExteriorPanel > >& currentPanels_ = + ( i == 0 ) ? bodyFixedPanels_ : segmentFixedPanels_.at( segmentFixedPanelsIterator->first ); + for( unsigned int j = 0; j < currentPanels_.size( ); j++ ) + { + surfaceNormals_[ counter ] = currentOrientation * currentPanels_.at( j )->getFrameFixedSurfaceNormal( )( ); + surfacePanelCosines_[ counter ] = (-sourceToTargetDirectionLocalFrame).dot(surfaceNormals_[ counter ]); + + if (surfacePanelCosines_[ counter ] > 0) + { + Eigen::Vector3d panelForce = radiationPressure * currentPanels_.at( j )->getPanelArea() * surfacePanelCosines_[ counter ] * + currentPanels_.at( j )->getReflectionLaw()->evaluateReactionVectorPartialWrtDiffuseReflectivity(surfaceNormals_[ counter ], sourceToTargetDirectionLocalFrame ); + forcePartialWrtDiffuseReflectivity += panelForce; + } + + counter++; + } + if( i > 0 ) + { + segmentFixedPanelsIterator++; + } + } + return forcePartialWrtDiffuseReflectivity; +} + +Eigen::Vector3d PaneledRadiationPressureTargetModel::evaluateRadiationPressureForcePartialWrtSpecularReflectivity( + double sourceIrradiance, + const Eigen::Vector3d& sourceToTargetDirectionLocalFrame) +{ + Eigen::Vector3d forcePartialWrtSpecularReflectivity = Eigen::Vector3d::Zero(); + + double radiationPressure = sourceIrradiance / physical_constants::SPEED_OF_LIGHT; + auto segmentFixedPanelsIterator = segmentFixedPanels_.begin( ); + int counter = 0; + Eigen::Quaterniond currentOrientation; + Eigen::Vector3d currentCenterOfMass = Eigen::Vector3d::Constant( TUDAT_NAN ); + + for( unsigned int i = 0; i < segmentFixedPanels_.size( ) + 1; i++ ) + { + currentOrientation = Eigen::Quaterniond( Eigen::Matrix3d::Identity( ) ); + if( i > 0 ) + { + currentOrientation = segmentFixedToBodyFixedRotations_.at( segmentFixedPanelsIterator->first )( ); + + } + + const std::vector< std::shared_ptr< system_models::VehicleExteriorPanel > >& currentPanels_ = + ( i == 0 ) ? bodyFixedPanels_ : segmentFixedPanels_.at( segmentFixedPanelsIterator->first ); + for( unsigned int j = 0; j < currentPanels_.size( ); j++ ) + { + surfaceNormals_[ counter ] = currentOrientation * currentPanels_.at( j )->getFrameFixedSurfaceNormal( )( ); + surfacePanelCosines_[ counter ] = (-sourceToTargetDirectionLocalFrame).dot(surfaceNormals_[ counter ]); + + if (surfacePanelCosines_[ counter ] > 0) + { + Eigen::Vector3d panelForce = radiationPressure * currentPanels_.at( j )->getPanelArea() * surfacePanelCosines_[ counter ] * + currentPanels_.at( j )->getReflectionLaw()->evaluateReactionVectorPartialWrtSpecularReflectivity(surfaceNormals_[ counter ], sourceToTargetDirectionLocalFrame ); + forcePartialWrtSpecularReflectivity += panelForce; + } + + counter++; + } + if( i > 0 ) + { + segmentFixedPanelsIterator++; + } + } + return forcePartialWrtSpecularReflectivity; +} + void PaneledRadiationPressureTargetModel::updateMembers_(double currentTime) { diff --git a/src/astro/electromagnetism/reflectionLaw.cpp b/src/astro/electromagnetism/reflectionLaw.cpp index 08518a6248..b3b57c2931 100644 --- a/src/astro/electromagnetism/reflectionLaw.cpp +++ b/src/astro/electromagnetism/reflectionLaw.cpp @@ -127,5 +127,59 @@ Eigen::Vector3d computeMirrorlikeReflection( } } +Eigen::Vector3d SpecularDiffuseMixReflectionLaw::evaluateReactionVectorPartialWrtSpecularReflectivity( + const Eigen::Vector3d& surfaceNormal, + const Eigen::Vector3d& incomingDirection) const +{ + const double cosBetweenNormalAndIncoming = surfaceNormal.dot(-incomingDirection); + if (cosBetweenNormalAndIncoming <= 0) + { + // Radiation is incident on backside of surface + return Eigen::Vector3d::Zero(); + } + Eigen::Vector3d pureSpecularPartial = - 2 * cosBetweenNormalAndIncoming * surfaceNormal; + Eigen::Vector3d reactionFromInstantaneousReradiation; + if (withInstantaneousReradiation_) + { + reactionFromInstantaneousReradiation = -2. / 3 * surfaceNormal; + } + else + { + reactionFromInstantaneousReradiation = Eigen::Vector3d::Zero(); + } + Eigen::Vector3d adjustmentForAbsorptivity = incomingDirection + reactionFromInstantaneousReradiation; + + Eigen::Vector3d partial = pureSpecularPartial - adjustmentForAbsorptivity; + return partial; + +}; + +Eigen::Vector3d SpecularDiffuseMixReflectionLaw::evaluateReactionVectorPartialWrtDiffuseReflectivity( + const Eigen::Vector3d& surfaceNormal, + const Eigen::Vector3d& incomingDirection) const +{ + const double cosBetweenNormalAndIncoming = surfaceNormal.dot(-incomingDirection); + if (cosBetweenNormalAndIncoming <= 0) + { + // Radiation is incident on backside of surface + return Eigen::Vector3d::Zero(); + } + Eigen::Vector3d pureDiffusePartial = incomingDirection - 2. / 3 * surfaceNormal; + Eigen::Vector3d reactionFromInstantaneousReradiation; + if (withInstantaneousReradiation_) + { + reactionFromInstantaneousReradiation = -2. / 3 * surfaceNormal; + } + else + { + reactionFromInstantaneousReradiation = Eigen::Vector3d::Zero(); + } + Eigen::Vector3d adjustmentForAbsorptivity = incomingDirection + reactionFromInstantaneousReradiation; + Eigen::Vector3d partial = pureDiffusePartial - adjustmentForAbsorptivity; + + return partial; + +}; + } // tudat } // electromagnetism diff --git a/src/astro/orbit_determination/acceleration_partials/fullRadiationPressureAccelerationPartial.cpp b/src/astro/orbit_determination/acceleration_partials/fullRadiationPressureAccelerationPartial.cpp index 062113a76d..1a0da6281d 100755 --- a/src/astro/orbit_determination/acceleration_partials/fullRadiationPressureAccelerationPartial.cpp +++ b/src/astro/orbit_determination/acceleration_partials/fullRadiationPressureAccelerationPartial.cpp @@ -75,6 +75,57 @@ void RadiationPressureAccelerationPartial::wrtRadiationPressureCoefficient( partial = radiationPressureAcceleration_->getAcceleration( ) / targetModel->getCoefficient( ); } + +void RadiationPressureAccelerationPartial::wrtSpecularReflectivity( + Eigen::MatrixXd& partial, + std::shared_ptr< electromagnetism::PaneledRadiationPressureTargetModel > targetModel, + const std::string& panelTypeId) +{ + Eigen::Vector3d partial_temp = Eigen::Vector3d::Zero(); + partial = partial_temp; + + std::function targetMassFunction = radiationPressureAcceleration_->getTargetMassFunction(); + double spacecraftMass = targetMassFunction(); + + std::function targetRotationFromLocalToGlobalFrameFunction = radiationPressureAcceleration_->getTargetRotationFromLocalToGlobalFrameFunction(); + Eigen::Quaterniond targetRotationFromGlobalToLocalFrame = targetRotationFromLocalToGlobalFrameFunction().inverse( ); + std::function sourceCenterPositionInGlobalFrameFunction = radiationPressureAcceleration_->getSourcePositionFunction(); + std::function targetCenterPositionInGlobalFrameFunction = radiationPressureAcceleration_->getTargetPositionFunction(); + Eigen::Vector3d targetCenterPositionInGlobalFrame = targetCenterPositionInGlobalFrameFunction() - sourceCenterPositionInGlobalFrameFunction(); + Eigen::Vector3d sourceToTargetDirectionLocalFrame = targetRotationFromGlobalToLocalFrame * targetCenterPositionInGlobalFrame.normalized( ); + double receivedIrradiance = radiationPressureAcceleration_->getReceivedIrradiance(); + if (receivedIrradiance >= 0) + { + Eigen::Vector3d forcePartialWrtSpecularReflectivity = targetRotationFromLocalToGlobalFrameFunction() * + targetModel->evaluateRadiationPressureForcePartialWrtSpecularReflectivity(receivedIrradiance, sourceToTargetDirectionLocalFrame); + partial += forcePartialWrtSpecularReflectivity/spacecraftMass; + } + +} +void RadiationPressureAccelerationPartial::wrtDiffuseReflectivity( + Eigen::MatrixXd& partial, + std::shared_ptr< electromagnetism::PaneledRadiationPressureTargetModel > targetModel, + const std::string& panelTypeId) +{ + std::function targetMassFunction = radiationPressureAcceleration_->getTargetMassFunction(); + double spacecraftMass = targetMassFunction(); + std::function targetRotationFromLocalToGlobalFrameFunction = radiationPressureAcceleration_->getTargetRotationFromLocalToGlobalFrameFunction(); + Eigen::Quaterniond targetRotationFromGlobalToLocalFrame = targetRotationFromLocalToGlobalFrameFunction().inverse( ); + std::function sourceCenterPositionInGlobalFrameFunction = radiationPressureAcceleration_->getSourcePositionFunction(); + std::function targetCenterPositionInGlobalFrameFunction = radiationPressureAcceleration_->getTargetPositionFunction(); + Eigen::Vector3d targetCenterPositionInGlobalFrame = targetCenterPositionInGlobalFrameFunction() - sourceCenterPositionInGlobalFrameFunction(); + Eigen::Vector3d sourceToTargetDirectionLocalFrame = targetRotationFromGlobalToLocalFrame * targetCenterPositionInGlobalFrame.normalized( ); + + double receivedIrradiance = radiationPressureAcceleration_->getReceivedIrradiance(); + if (receivedIrradiance >= 0) + { + Eigen::Vector3d forcePartialWrtDiffuseReflectivity = targetRotationFromLocalToGlobalFrameFunction() * + targetModel->evaluateRadiationPressureForcePartialWrtDiffuseReflectivity(receivedIrradiance, sourceToTargetDirectionLocalFrame); + partial += forcePartialWrtDiffuseReflectivity/spacecraftMass; + } +} + + } } diff --git a/src/astro/orbit_determination/acceleration_partials/panelledRadiationPressureAccelerationPartial.cpp b/src/astro/orbit_determination/acceleration_partials/panelledRadiationPressureAccelerationPartial.cpp index cf807d0f6c..016b9fcaef 100755 --- a/src/astro/orbit_determination/acceleration_partials/panelledRadiationPressureAccelerationPartial.cpp +++ b/src/astro/orbit_determination/acceleration_partials/panelledRadiationPressureAccelerationPartial.cpp @@ -80,6 +80,59 @@ void PanelledRadiationPressurePartial::update( const double currentTime ) } +void PanelledRadiationPressurePartial::wrtDiffuseReflectivity( + Eigen::MatrixXd& partial, + const std::string& panelTypeId) +{ + Eigen::Vector3d partial_temp = Eigen::Vector3d::Zero(); + partial = partial_temp; + + std::function targetMassFunction = radiationPressureAcceleration_->getTargetMassFunction(); + double spacecraftMass = targetMassFunction(); + + std::function targetRotationFromLocalToGlobalFrameFunction = radiationPressureAcceleration_->getTargetRotationFromLocalToGlobalFrameFunction(); + Eigen::Quaterniond targetRotationFromGlobalToLocalFrame = targetRotationFromLocalToGlobalFrameFunction().inverse( ); + std::function sourceCenterPositionInGlobalFrameFunction = radiationPressureAcceleration_->getSourcePositionFunction(); + std::function targetCenterPositionInGlobalFrameFunction = radiationPressureAcceleration_->getTargetPositionFunction(); + Eigen::Vector3d targetCenterPositionInGlobalFrame = targetCenterPositionInGlobalFrameFunction() - sourceCenterPositionInGlobalFrameFunction(); + Eigen::Vector3d sourceToTargetDirectionLocalFrame = targetRotationFromGlobalToLocalFrame * targetCenterPositionInGlobalFrame.normalized( ); + double receivedIrradiance = radiationPressureAcceleration_->getReceivedIrradiance(); + if (receivedIrradiance >= 0) + { + Eigen::Vector3d forcePartialWrtDiffuseReflectivity = targetRotationFromLocalToGlobalFrameFunction() * + panelledTargetModel_->evaluateRadiationPressureForcePartialWrtDiffuseReflectivity(receivedIrradiance, sourceToTargetDirectionLocalFrame); + partial += forcePartialWrtDiffuseReflectivity/spacecraftMass; + } + +} + +void PanelledRadiationPressurePartial::wrtSpecularReflectivity( + Eigen::MatrixXd& partial, + const std::string& panelTypeId) +{ + Eigen::Vector3d partial_temp = Eigen::Vector3d::Zero(); + partial = partial_temp; + + std::function targetMassFunction = radiationPressureAcceleration_->getTargetMassFunction(); + double spacecraftMass = targetMassFunction(); + + std::function targetRotationFromLocalToGlobalFrameFunction = radiationPressureAcceleration_->getTargetRotationFromLocalToGlobalFrameFunction(); + Eigen::Quaterniond targetRotationFromGlobalToLocalFrame = targetRotationFromLocalToGlobalFrameFunction().inverse( ); + std::function sourceCenterPositionInGlobalFrameFunction = radiationPressureAcceleration_->getSourcePositionFunction(); + std::function targetCenterPositionInGlobalFrameFunction = radiationPressureAcceleration_->getTargetPositionFunction(); + Eigen::Vector3d targetCenterPositionInGlobalFrame = targetCenterPositionInGlobalFrameFunction() - sourceCenterPositionInGlobalFrameFunction(); + Eigen::Vector3d sourceToTargetDirectionLocalFrame = targetRotationFromGlobalToLocalFrame * targetCenterPositionInGlobalFrame.normalized( ); + double receivedIrradiance = radiationPressureAcceleration_->getReceivedIrradiance(); + if (receivedIrradiance >= 0) + { + Eigen::Vector3d forcePartialWrtSpecularReflectivity = targetRotationFromLocalToGlobalFrameFunction() * + panelledTargetModel_->evaluateRadiationPressureForcePartialWrtSpecularReflectivity(receivedIrradiance, sourceToTargetDirectionLocalFrame); + partial += forcePartialWrtSpecularReflectivity/spacecraftMass; + } + +} + + } } diff --git a/src/astro/orbit_determination/estimatable_parameters/CMakeLists.txt b/src/astro/orbit_determination/estimatable_parameters/CMakeLists.txt index 05287e3716..972bc496bd 100644 --- a/src/astro/orbit_determination/estimatable_parameters/CMakeLists.txt +++ b/src/astro/orbit_determination/estimatable_parameters/CMakeLists.txt @@ -43,6 +43,8 @@ set(estimatable_parameters_HEADERS "polarMotionAmplitude.h" "coreFactor.h" "freeCoreNutationRate.h" + "specularReflectivity.h" + "diffuseReflectivity.h" "polynomialClockCorrections.h" ) diff --git a/src/astro/orbit_determination/estimatable_parameters/estimatableParameter.cpp b/src/astro/orbit_determination/estimatable_parameters/estimatableParameter.cpp index 38dbb3b966..49475ed5c4 100644 --- a/src/astro/orbit_determination/estimatable_parameters/estimatableParameter.cpp +++ b/src/astro/orbit_determination/estimatable_parameters/estimatableParameter.cpp @@ -165,6 +165,12 @@ std::string getParameterTypeString( const EstimatebleParametersEnum parameterTyp case source_perpendicular_direction_radiation_pressure_scaling_factor: parameterDescription = " Radiation pressure acceleration scaling factor perpendicular to source "; break; + case specular_reflectivity: + parameterDescription = " specular reflectivity for panel group "; + break; + case diffuse_reflectivity: + parameterDescription = " diffuse reflectivity for panel group "; + break; case mode_coupled_tidal_love_numbers: parameterDescription = " Mode-coupled tidal Love numbers"; break; @@ -342,6 +348,12 @@ bool isDoubleParameter( const EstimatebleParametersEnum parameterType ) case source_perpendicular_direction_radiation_pressure_scaling_factor: isDoubleParameter = true; break; + case specular_reflectivity: + isDoubleParameter = true; + break; + case diffuse_reflectivity: + isDoubleParameter = true; + break; case mode_coupled_tidal_love_numbers: isDoubleParameter = false; break; diff --git a/src/simulation/environment_setup/createSystemModel.cpp b/src/simulation/environment_setup/createSystemModel.cpp index 2583324100..efc1510d3e 100644 --- a/src/simulation/environment_setup/createSystemModel.cpp +++ b/src/simulation/environment_setup/createSystemModel.cpp @@ -147,6 +147,10 @@ std::pair< std::shared_ptr< system_models::VehicleExteriorPanel >, std::string > std::shared_ptr< electromagnetism::ReflectionLaw > reflectionLaw = createReflectionLaw( panelSettings->reflectionLawSettings_ ); exteriorPanel->setReflectionLaw( reflectionLaw ); + if( panelSettings->panelTypeId_ != "" ) + { + exteriorPanel->setPanelTypeId( panelSettings->panelTypeId_ ); + } } if( bodies.at( bodyName )->getRotationalEphemeris( ) == nullptr ) diff --git a/tests/src/astro/orbit_determination/acceleration_partials/unitTestAccelerationPartials.cpp b/tests/src/astro/orbit_determination/acceleration_partials/unitTestAccelerationPartials.cpp index 37fa80e2d2..23e033eccd 100644 --- a/tests/src/astro/orbit_determination/acceleration_partials/unitTestAccelerationPartials.cpp +++ b/tests/src/astro/orbit_determination/acceleration_partials/unitTestAccelerationPartials.cpp @@ -70,6 +70,222 @@ using namespace tudat::propulsion; BOOST_AUTO_TEST_SUITE( test_acceleration_partials ) +BOOST_AUTO_TEST_CASE( testPanelledRadiationPressureAccelerationPartials ) +{ + tudat::spice_interface::loadStandardSpiceKernels( ); + + for( int testIndex = 0; testIndex < 6; testIndex++ ) + { + // Create empty bodies, earth and sun. + std::shared_ptr< Body > vehicle = std::make_shared< Body >( ); + double vehicleMass = 400.0; + vehicle->setConstantBodyMass( vehicleMass ); + std::shared_ptr< Body > sun = std::make_shared< Body >( ); + SystemOfBodies bodies; + bodies.addBody( vehicle, "Vehicle" ); + bodies.addBody( sun, "Sun" ); + + // Load spice kernels. + + // Set current state of sun and earth. + sun->setState( Eigen::Vector6d::Zero( ) );//getBodyCartesianStateAtEpoch( "Sun", "SSB", "J2000", "NONE", 1.0E6 ) ); + sun->setRadiationSourceModel( + createRadiationSourceModel( + getDefaultRadiationSourceModelSettings( "Sun", TUDAT_NAN, TUDAT_NAN ), + "Sun", bodies ) ); + vehicle->setState( + ( Eigen::Vector6d( ) << 1.4E11, 1.0E11, 1.1E11, 0.0, 0.0, 0.0 ).finished( ) );//getBodyCartesianStateAtEpoch( "Earth", "SSB", "J2000", "NONE", 1.0E6 ) ); + vehicle->setRotationalEphemeris( + // std::make_shared< ephemerides::CustomRotationalEphemeris >( + // [=](const double){return Eigen::Quaterniond( Eigen::Matrix3d::Identity( ) ); }, "ECLIPJ2000", "VehicleFixed" ) ); + std::make_shared< tudat::ephemerides::SimpleRotationalEphemeris >( 0.2, 0.4, -0.2, 1.0E-5, 0.0, "ECLIPJ2000", "VehicleFixed" ) ); + vehicle->setCurrentRotationalStateToLocalFrameFromEphemeris( 0.0 ); + + // Create links to set and get state functions of bodies. + std::function< void( Eigen::Vector6d ) > sunStateSetFunction = + std::bind( &Body::setState, sun, std::placeholders::_1 ); + std::function< void( Eigen::Vector6d ) > vehicleStateSetFunction = + std::bind( &Body::setState, vehicle, std::placeholders::_1 ); + std::function< Eigen::Vector6d( ) > sunStateGetFunction = + std::bind( &Body::getState, sun ); + std::function< Eigen::Vector6d( ) > vehicleStateGetFunction = + std::bind( &Body::getState, vehicle ); + + double specularReflectivity = 0.5; + double diffuseReflectivity = 0.1; + bool useInstantaneousReradiation = false; + if( testIndex % 2 == 1 ) + { + useInstantaneousReradiation = true; + } + std::string panelName = "SolarPanel"; + double scalingValue = 1.0; + if( testIndex > 3 ) + { + scalingValue = 0.5; + panelName = "PanelName"; + } + std::vector< std::shared_ptr< BodyPanelSettings > > panelSettingsList; + panelSettingsList.push_back( std::make_shared< BodyPanelSettings >( + std::make_shared< FrameFixedBodyPanelGeometrySettings >( -Eigen::Vector3d::UnitX( ), 1.0 ), + std::make_shared< SpecularDiffuseBodyPanelReflectionLawSettings >( specularReflectivity, diffuseReflectivity, useInstantaneousReradiation ), "SolarPanel" ) ); + if( testIndex > 2 ) + { + panelSettingsList.push_back( std::make_shared< BodyPanelSettings >( + std::make_shared< FrameFixedBodyPanelGeometrySettings >( -Eigen::Vector3d::UnitY( ), 3.254 ), + std::make_shared< SpecularDiffuseBodyPanelReflectionLawSettings >( specularReflectivity, diffuseReflectivity, useInstantaneousReradiation ), "SolarPanel" ) ); + panelSettingsList.push_back( std::make_shared< BodyPanelSettings >( + std::make_shared< FrameFixedBodyPanelGeometrySettings >( -Eigen::Vector3d::UnitZ( ), 8.654 ), + std::make_shared< SpecularDiffuseBodyPanelReflectionLawSettings >( specularReflectivity, diffuseReflectivity, useInstantaneousReradiation ), "SolarPanel" ) ); + panelSettingsList.push_back( std::make_shared< BodyPanelSettings >( + std::make_shared< FrameFixedBodyPanelGeometrySettings >( Eigen::Vector3d::UnitX( ), 1.346 ), + std::make_shared< SpecularDiffuseBodyPanelReflectionLawSettings >( + scalingValue * specularReflectivity, scalingValue * diffuseReflectivity, useInstantaneousReradiation ), panelName ) ); + panelSettingsList.push_back( std::make_shared< BodyPanelSettings >( + std::make_shared< FrameFixedBodyPanelGeometrySettings >( Eigen::Vector3d::UnitY( ), 10.4783 ), + std::make_shared< SpecularDiffuseBodyPanelReflectionLawSettings >( + scalingValue * specularReflectivity, scalingValue * diffuseReflectivity, useInstantaneousReradiation ), panelName ) ); + panelSettingsList.push_back( std::make_shared< BodyPanelSettings >( + std::make_shared< FrameFixedBodyPanelGeometrySettings >( Eigen::Vector3d::UnitZ( ), 6.4235 ), + std::make_shared< SpecularDiffuseBodyPanelReflectionLawSettings >( + scalingValue * specularReflectivity, scalingValue * diffuseReflectivity, useInstantaneousReradiation ), panelName ) ); + } + + addBodyExteriorPanelledShape( + std::make_shared< FullPanelledBodySettings >( panelSettingsList ), "Vehicle", bodies ); + + auto paneledRadiationPressureTargetSettings = + std::make_shared< RadiationPressureTargetModelSettings >( paneled_target ); + + std::shared_ptr radiationPressureInterface = + std::dynamic_pointer_cast( + createRadiationPressureTargetModel( + paneledRadiationPressureTargetSettings, "Vehicle", bodies ).at( 0 ) ); + vehicle->addRadiationPressureTargetModel( radiationPressureInterface ); + + // Create acceleration model. + std::shared_ptr< RadiationPressureAcceleration > accelerationModel = + std::dynamic_pointer_cast< RadiationPressureAcceleration >( + createRadiationPressureAccelerationModel( + vehicle, sun, "Vehicle", "Sun", bodies ) ); + accelerationModel->updateMembers( 0.0 ); + // Create partial-calculating object. + std::shared_ptr< AccelerationPartial > accelerationPartial = + createAnalyticalAccelerationPartial( accelerationModel, { "Vehicle", vehicle }, { "Sun", sun}, bodies ); + + std::shared_ptr< EstimatableParameter< double > > parallelScalingFactor = + std::make_shared< RadiationPressureScalingFactor >( accelerationModel, source_direction_radiation_pressure_scaling_factor, "Vehicle", "Sun" ); + std::shared_ptr< EstimatableParameter< double > > perpendicularScalingFactor = + std::make_shared< RadiationPressureScalingFactor >( accelerationModel, source_perpendicular_direction_radiation_pressure_scaling_factor, "Vehicle", "Sun" ); + std::shared_ptr< EstimatableParameter< double > > diffuseReflectivityParameter = + std::make_shared< DiffuseReflectivity >( radiationPressureInterface, "Vehicle", "SolarPanel" ); + std::shared_ptr< EstimatableParameter< double > > specularReflectivityParameter = + std::make_shared< SpecularReflectivity >( radiationPressureInterface, "Vehicle", "SolarPanel" ); + + // Calculate analytical partials. + accelerationPartial->update( 0.0 ); + Eigen::MatrixXd partialWrtSunPosition = Eigen::Matrix3d::Zero( ); + accelerationPartial->wrtPositionOfAcceleratingBody( partialWrtSunPosition.block( 0, 0, 3, 3 ) ); + + Eigen::MatrixXd partialWrtSunVelocity = Eigen::Matrix3d::Zero( ); + accelerationPartial->wrtVelocityOfAcceleratingBody( partialWrtSunVelocity.block( 0, 0, 3, 3 ) ); + + Eigen::MatrixXd partialWrtVehiclePosition = Eigen::Matrix3d::Zero( ); + accelerationPartial->wrtPositionOfAcceleratedBody( partialWrtVehiclePosition.block( 0, 0, 3, 3 ) ); + + Eigen::MatrixXd partialWrtVehicleVelocity = Eigen::Matrix3d::Zero( ); + accelerationPartial->wrtVelocityOfAcceleratedBody( partialWrtVehicleVelocity.block( 0, 0, 3, 3 ) ); + + //Eigen::MatrixXd partialWrtEmissivities = accelerationPartial->wrtParameter( panelEmissivitiesParameter ); + Eigen::MatrixXd partialWrtParallelScaling = accelerationPartial->wrtParameter( parallelScalingFactor ); + Eigen::MatrixXd partialWrtPerpendicularScaling = accelerationPartial->wrtParameter( perpendicularScalingFactor ); + Eigen::Vector3d partialWrtSpecularReflectivity = accelerationPartial->wrtParameter( specularReflectivityParameter ); + Eigen::Vector3d partialWrtDiffuseReflectivity = accelerationPartial->wrtParameter( diffuseReflectivityParameter ); + + // Declare numerical partials. + Eigen::Matrix3d testPartialWrtVehiclePosition = Eigen::Matrix3d::Zero( ); + Eigen::Matrix3d testPartialWrtVehicleVelocity = Eigen::Matrix3d::Zero( ); + Eigen::Matrix3d testPartialWrtSunPosition = Eigen::Matrix3d::Zero( ); + Eigen::Matrix3d testPartialWrtSunVelocity = Eigen::Matrix3d::Zero( ); + Eigen::MatrixXd testPartialWrtEmissivities = Eigen::MatrixXd::Zero( 3, 2 ); + Eigen::MatrixXd testPartialWrtParallelScaling = Eigen::MatrixXd::Zero( 3, 1 ); + Eigen::MatrixXd testPartialWrtPerpendicularScaling = Eigen::MatrixXd::Zero( 3, 1 ); + Eigen::Vector3d testPartialWrtSpecularReflectivity = Eigen::Vector3d::Zero( ); + Eigen::Vector3d testPartialWrtDiffuseReflectivity = Eigen::Vector3d::Zero( ); + + // Declare perturbations in position for numerical partial/ + Eigen::Vector3d positionPerturbation; + positionPerturbation<< 10000.0, 10000.0, 10000.0; + Eigen::Vector3d velocityPerturbation; + velocityPerturbation<< 0.1, 0.1, 0.1; + + // Calculate numerical partials. + std::function< void( const double ) > updateFunction1 = std::bind( + &RadiationPressureTargetModel::updateMembers, vehicle->getRadiationPressureTargetModel( ), std::placeholders::_1 ); + std::function< void( const double ) > updateFunction2 = std::bind( + &RadiationSourceModel::updateMembers, sun->getRadiationSourceModel( ), std::placeholders::_1 ); + std::function< void( ) > updateFunction = [=]( ){ + updateFunction2( 0.0 ); + updateFunction1( 0.0 ); + }; + + testPartialWrtSunPosition = calculateAccelerationWrtStatePartials( + sunStateSetFunction, accelerationModel, sun->getState( ), positionPerturbation, 0, updateFunction ); + testPartialWrtVehiclePosition = calculateAccelerationWrtStatePartials( + vehicleStateSetFunction, accelerationModel, vehicle->getState( ), positionPerturbation, 0, updateFunction ); + testPartialWrtSunVelocity = calculateAccelerationWrtStatePartials( + sunStateSetFunction, accelerationModel, sun->getState( ),velocityPerturbation, 3, updateFunction ); + testPartialWrtVehicleVelocity = calculateAccelerationWrtStatePartials( + vehicleStateSetFunction, accelerationModel, vehicle->getState( ), velocityPerturbation, 3, updateFunction ); + // testPartialWrtEmissivities = calculateAccelerationWrtParameterPartials( + // panelEmissivitiesParameter, accelerationModel, emissivityPerturbations ); + + testPartialWrtParallelScaling = calculateAccelerationWrtParameterPartials( + parallelScalingFactor, accelerationModel, 10.0, updateFunction ); + testPartialWrtPerpendicularScaling = calculateAccelerationWrtParameterPartials( + perpendicularScalingFactor, accelerationModel, 10.0, updateFunction ); + + testPartialWrtDiffuseReflectivity = calculateAccelerationWrtParameterPartials( + diffuseReflectivityParameter, accelerationModel, 0.1, updateFunction ); + testPartialWrtSpecularReflectivity = calculateAccelerationWrtParameterPartials( + specularReflectivityParameter, accelerationModel, 0.1, updateFunction ); + + std::cout<::epsilon( ) ); + TUDAT_CHECK_MATRIX_CLOSE_FRACTION( testPartialWrtVehiclePosition, + partialWrtVehiclePosition, 1.0e-6 ); + TUDAT_CHECK_MATRIX_CLOSE_FRACTION( testPartialWrtVehicleVelocity, + partialWrtVehicleVelocity, std::numeric_limits< double >::epsilon( ) ); + TUDAT_CHECK_MATRIX_CLOSE_FRACTION( testPartialWrtParallelScaling, + partialWrtParallelScaling, 1.0E-13 ); + TUDAT_CHECK_MATRIX_CLOSE_FRACTION( testPartialWrtPerpendicularScaling, + partialWrtPerpendicularScaling, 1.0E-13 ); + TUDAT_CHECK_MATRIX_CLOSE_FRACTION( partialWrtSpecularReflectivity, + testPartialWrtSpecularReflectivity, 1.0E-8 ); + if( testIndex % 2 == 0 ) + { + TUDAT_CHECK_MATRIX_CLOSE_FRACTION( partialWrtDiffuseReflectivity, + testPartialWrtDiffuseReflectivity, 1.0E-13 ); + } + else + { + BOOST_CHECK_EQUAL( testPartialWrtDiffuseReflectivity.norm( ), 0.0 ); + BOOST_CHECK_EQUAL( partialWrtDiffuseReflectivity.norm( ), 0.0 ); + } + } + +} + + BOOST_AUTO_TEST_CASE( testCentralGravityPartials ) { // Create empty bodies, earth and sun. @@ -1202,188 +1418,6 @@ BOOST_AUTO_TEST_CASE( testDirectDissipationAccelerationPartial ) } -BOOST_AUTO_TEST_CASE( testPanelledRadiationPressureAccelerationPartials ) -{ - // Create empty bodies, earth and sun. - std::shared_ptr< Body > vehicle = std::make_shared< Body >( ); - double vehicleMass = 400.0; - vehicle->setConstantBodyMass( vehicleMass ); - std::shared_ptr< Body > sun = std::make_shared< Body >( ); - SystemOfBodies bodies; - bodies.addBody( vehicle, "Vehicle" ); - bodies.addBody( sun, "Sun" ); - - // Load spice kernels. - tudat::spice_interface::loadStandardSpiceKernels( ); - - // Set current state of sun and earth. - sun->setState( Eigen::Vector6d::Zero( ) );//getBodyCartesianStateAtEpoch( "Sun", "SSB", "J2000", "NONE", 1.0E6 ) ); - sun->setRadiationSourceModel( - createRadiationSourceModel( - getDefaultRadiationSourceModelSettings( "Sun", TUDAT_NAN, TUDAT_NAN ), - "Sun", bodies ) ); - vehicle->setState( - ( Eigen::Vector6d( ) << 1.4E11, 1.0E11, 1.1E11, 0.0, 0.0, 0.0 ).finished( ) );//getBodyCartesianStateAtEpoch( "Earth", "SSB", "J2000", "NONE", 1.0E6 ) ); - vehicle->setRotationalEphemeris( -// std::make_shared< ephemerides::CustomRotationalEphemeris >( -// [=](const double){return Eigen::Quaterniond( Eigen::Matrix3d::Identity( ) ); }, "ECLIPJ2000", "VehicleFixed" ) ); - std::make_shared< tudat::ephemerides::SimpleRotationalEphemeris >( 0.2, 0.4, -0.2, 1.0E-5, 0.0, "ECLIPJ2000", "VehicleFixed" ) ); - vehicle->setCurrentRotationalStateToLocalFrameFromEphemeris( 0.0 ); - - // Create links to set and get state functions of bodies. - std::function< void( Eigen::Vector6d ) > sunStateSetFunction = - std::bind( &Body::setState, sun, std::placeholders::_1 ); - std::function< void( Eigen::Vector6d ) > vehicleStateSetFunction = - std::bind( &Body::setState, vehicle, std::placeholders::_1 ); - std::function< Eigen::Vector6d( ) > sunStateGetFunction = - std::bind( &Body::getState, sun ); - std::function< Eigen::Vector6d( ) > vehicleStateGetFunction = - std::bind( &Body::getState, vehicle ); - - std::vector< std::shared_ptr< BodyPanelSettings > > panelSettingsList; - panelSettingsList.push_back( std::make_shared< BodyPanelSettings >( - std::make_shared< FrameFixedBodyPanelGeometrySettings >( -Eigen::Vector3d::UnitX( ), 1.0 ), - std::make_shared< SpecularDiffuseBodyPanelReflectionLawSettings >( 0.5, 0.1, true ) ) ); - panelSettingsList.push_back( std::make_shared< BodyPanelSettings >( - std::make_shared< FrameFixedBodyPanelGeometrySettings >( -Eigen::Vector3d::UnitY( ), 3.254 ), - std::make_shared< SpecularDiffuseBodyPanelReflectionLawSettings >( 0.4, 0.2, true ) ) ); - panelSettingsList.push_back( std::make_shared< BodyPanelSettings >( - std::make_shared< FrameFixedBodyPanelGeometrySettings >( -Eigen::Vector3d::UnitZ( ), 8.654 ), - std::make_shared< SpecularDiffuseBodyPanelReflectionLawSettings >( 0.2, 0.3, true ) ) ); - panelSettingsList.push_back( std::make_shared< BodyPanelSettings >( - std::make_shared< FrameFixedBodyPanelGeometrySettings >( Eigen::Vector3d::UnitX( ), 1.346 ), - std::make_shared< SpecularDiffuseBodyPanelReflectionLawSettings >( 0.25, 0.15, true ) ) ); - panelSettingsList.push_back( std::make_shared< BodyPanelSettings >( - std::make_shared< FrameFixedBodyPanelGeometrySettings >( Eigen::Vector3d::UnitY( ), 10.4783 ), - std::make_shared< SpecularDiffuseBodyPanelReflectionLawSettings >( 0.44, 0.51, true ) ) ); - panelSettingsList.push_back( std::make_shared< BodyPanelSettings >( - std::make_shared< FrameFixedBodyPanelGeometrySettings >( Eigen::Vector3d::UnitZ( ), 6.4235 ), - std::make_shared< SpecularDiffuseBodyPanelReflectionLawSettings >( 0.16, 0.34, true ) ) ); - - addBodyExteriorPanelledShape( - std::make_shared< FullPanelledBodySettings >( panelSettingsList ), "Vehicle", bodies ); - vehicle->addRadiationPressureTargetModel( - createRadiationPressureTargetModel( - std::make_shared< RadiationPressureTargetModelSettings >( paneled_target ), "Vehicle", bodies ).at( 0 ) ); - - // Create acceleration model. - std::shared_ptr< RadiationPressureAcceleration > accelerationModel = - std::dynamic_pointer_cast< RadiationPressureAcceleration >( - createRadiationPressureAccelerationModel( - vehicle, sun, "Vehicle", "Sun", bodies ) ); - accelerationModel->updateMembers( 0.0 ); - // Create partial-calculating object. - std::shared_ptr< AccelerationPartial > accelerationPartial = - createAnalyticalAccelerationPartial( accelerationModel, { "Vehicle", vehicle }, { "Sun", sun}, bodies ); - - // std::vector< int > panelIndices1 = boost::assign::list_of( 0 )( 6 ); - // std::vector< int > panelIndices2 = boost::assign::list_of( 2 ); - - // std::vector< std::vector< int > > panelIndices; - // panelIndices.push_back( panelIndices2 ); - // panelIndices.push_back( panelIndices1 ); - - // std::shared_ptr< EstimatableParameterSettings > panelEmissivitiesSettings = - // std::make_shared< PanelRadiationEmissivitiesParameterSettings >( "Vehicle", panelIndices ); - // std::vector< std::shared_ptr< EstimatableParameterSettings > > parameterSettingsVector; - // parameterSettingsVector.push_back( panelEmissivitiesSettings ); - // std::shared_ptr< estimatable_parameters::EstimatableParameterSet< double > > estimatableParameters = createParametersToEstimate( - // parameterSettingsVector, bodies ); - // std::shared_ptr< EstimatableParameter< Eigen::VectorXd > > panelEmissivitiesParameter = - // estimatableParameters->getVectorParameters( ).begin( )->second; - - - std::shared_ptr< EstimatableParameter< double > > parallelScalingFactor = - std::make_shared< RadiationPressureScalingFactor >( accelerationModel, source_direction_radiation_pressure_scaling_factor, "Vehicle", "Sun" ); - std::shared_ptr< EstimatableParameter< double > > perpendicularScalingFactor = - std::make_shared< RadiationPressureScalingFactor >( accelerationModel, source_perpendicular_direction_radiation_pressure_scaling_factor, "Vehicle", "Sun" ); - - - // Calculate analytical partials. - accelerationPartial->update( 0.0 ); - Eigen::MatrixXd partialWrtSunPosition = Eigen::Matrix3d::Zero( ); - accelerationPartial->wrtPositionOfAcceleratingBody( partialWrtSunPosition.block( 0, 0, 3, 3 ) ); - - Eigen::MatrixXd partialWrtSunVelocity = Eigen::Matrix3d::Zero( ); - accelerationPartial->wrtVelocityOfAcceleratingBody( partialWrtSunVelocity.block( 0, 0, 3, 3 ) ); - - Eigen::MatrixXd partialWrtVehiclePosition = Eigen::Matrix3d::Zero( ); - accelerationPartial->wrtPositionOfAcceleratedBody( partialWrtVehiclePosition.block( 0, 0, 3, 3 ) ); - - Eigen::MatrixXd partialWrtVehicleVelocity = Eigen::Matrix3d::Zero( ); - accelerationPartial->wrtVelocityOfAcceleratedBody( partialWrtVehicleVelocity.block( 0, 0, 3, 3 ) ); - - //Eigen::MatrixXd partialWrtEmissivities = accelerationPartial->wrtParameter( panelEmissivitiesParameter ); - Eigen::MatrixXd partialWrtParallelScaling = accelerationPartial->wrtParameter( parallelScalingFactor ); - Eigen::MatrixXd partialWrtPerpendicularScaling = accelerationPartial->wrtParameter( perpendicularScalingFactor ); - - // Declare numerical partials. - Eigen::Matrix3d testPartialWrtVehiclePosition = Eigen::Matrix3d::Zero( ); - Eigen::Matrix3d testPartialWrtVehicleVelocity = Eigen::Matrix3d::Zero( ); - Eigen::Matrix3d testPartialWrtSunPosition = Eigen::Matrix3d::Zero( ); - Eigen::Matrix3d testPartialWrtSunVelocity = Eigen::Matrix3d::Zero( ); - Eigen::MatrixXd testPartialWrtEmissivities = Eigen::MatrixXd::Zero( 3, 2 ); - Eigen::MatrixXd testPartialWrtParallelScaling = Eigen::MatrixXd::Zero( 3, 1 ); - Eigen::MatrixXd testPartialWrtPerpendicularScaling = Eigen::MatrixXd::Zero( 3, 1 ); - - // Declare perturbations in position for numerical partial/ - Eigen::Vector3d positionPerturbation; - positionPerturbation<< 10000.0, 10000.0, 10000.0; - Eigen::Vector3d velocityPerturbation; - velocityPerturbation<< 0.1, 0.1, 0.1; - // Eigen::VectorXd emissivityPerturbations = Eigen::VectorXd::Zero( 2 ); - // emissivityPerturbations( 0 ) = 1.0; - // emissivityPerturbations( 1 ) = 1.0; - - // Calculate numerical partials. - std::function< void( const double ) > updateFunction1 = std::bind( - &RadiationPressureTargetModel::updateMembers, vehicle->getRadiationPressureTargetModel( ), std::placeholders::_1 ); - std::function< void( const double ) > updateFunction2 = std::bind( - &RadiationSourceModel::updateMembers, sun->getRadiationSourceModel( ), std::placeholders::_1 ); - std::function< void( ) > updateFunction = [=]( ){ - updateFunction2( 0.0 ); - updateFunction1( 0.0 ); - }; - - testPartialWrtSunPosition = calculateAccelerationWrtStatePartials( - sunStateSetFunction, accelerationModel, sun->getState( ), positionPerturbation, 0, updateFunction ); - testPartialWrtVehiclePosition = calculateAccelerationWrtStatePartials( - vehicleStateSetFunction, accelerationModel, vehicle->getState( ), positionPerturbation, 0, updateFunction ); - testPartialWrtSunVelocity = calculateAccelerationWrtStatePartials( - sunStateSetFunction, accelerationModel, sun->getState( ),velocityPerturbation, 3, updateFunction ); - testPartialWrtVehicleVelocity = calculateAccelerationWrtStatePartials( - vehicleStateSetFunction, accelerationModel, vehicle->getState( ), velocityPerturbation, 3, updateFunction ); - // testPartialWrtEmissivities = calculateAccelerationWrtParameterPartials( - // panelEmissivitiesParameter, accelerationModel, emissivityPerturbations ); - - testPartialWrtParallelScaling = calculateAccelerationWrtParameterPartials( - parallelScalingFactor, accelerationModel, 10.0, updateFunction ); - testPartialWrtPerpendicularScaling = calculateAccelerationWrtParameterPartials( - perpendicularScalingFactor, accelerationModel, 10.0, updateFunction ); - -// std::cout<::epsilon( ) ); - TUDAT_CHECK_MATRIX_CLOSE_FRACTION( testPartialWrtVehiclePosition, - partialWrtVehiclePosition, 1.0e-6 ); - TUDAT_CHECK_MATRIX_CLOSE_FRACTION( testPartialWrtVehicleVelocity, - partialWrtVehicleVelocity, std::numeric_limits< double >::epsilon( ) ); - // TUDAT_CHECK_MATRIX_CLOSE_FRACTION( testPartialWrtEmissivities, - // partialWrtEmissivities, 1.0e-14 ); - TUDAT_CHECK_MATRIX_CLOSE_FRACTION( testPartialWrtParallelScaling, - partialWrtParallelScaling, 1.0E-13 ); - TUDAT_CHECK_MATRIX_CLOSE_FRACTION( testPartialWrtPerpendicularScaling, - partialWrtPerpendicularScaling, 1.0E-13 ); - -} BOOST_AUTO_TEST_CASE( testPanelledSurfaceRadiationPressureAccelerationPartials ) @@ -1774,6 +1808,7 @@ BOOST_AUTO_TEST_CASE( testYarkovskyPartials ) partialWrtSunYarkovskyParameter, 1.0E-8 ); } + BOOST_AUTO_TEST_SUITE_END( ) } // namespace unit_tests