From cdd238da25270258919ce42b33dcee841a143168 Mon Sep 17 00:00:00 2001 From: Dominic Dirkx Date: Tue, 17 Oct 2023 08:48:34 +0200 Subject: [PATCH] Working on radiation pressure partials --- .../radiationPressureAcceleration.h | 10 ++++ .../radiationPressureTargetModel.h | 18 +++++-- .../radiationPressureAccelerationPartial.h | 18 +++++++ .../createAccelerationPartials.h | 49 ++++++++++++++++++- .../radiationPressureAcceleration.cpp | 44 +++++++++-------- .../radiationPressureTargetModel.cpp | 4 +- 6 files changed, 115 insertions(+), 28 deletions(-) diff --git a/include/tudat/astro/electromagnetism/radiationPressureAcceleration.h b/include/tudat/astro/electromagnetism/radiationPressureAcceleration.h index 9b26c7cbdf..598a0f451e 100644 --- a/include/tudat/astro/electromagnetism/radiationPressureAcceleration.h +++ b/include/tudat/astro/electromagnetism/radiationPressureAcceleration.h @@ -59,6 +59,11 @@ class RadiationPressureAcceleration: public basic_astrodynamics::AccelerationMod return receivedIrradiance; } + std::function getTargetMassFunction() + { + return targetMassFunction_; + } + protected: RadiationPressureAcceleration(const std::function& sourcePositionFunction, const std::shared_ptr& targetModel, @@ -140,8 +145,13 @@ class IsotropicPointSourceRadiationPressureAcceleration : public RadiationPressu Eigen::Vector3d calculateAcceleration() override; std::shared_ptr sourceModel_; + std::shared_ptr sourceBodyShapeModel_; + Eigen::Vector3d sourceCenterPositionInGlobalFrame_; + Eigen::Vector3d targetCenterPositionInGlobalFrame_; + Eigen::Vector3d targetCenterPositionInSourceFrame_; + // For dependent variable double sourceToTargetReceivedFraction; }; diff --git a/include/tudat/astro/electromagnetism/radiationPressureTargetModel.h b/include/tudat/astro/electromagnetism/radiationPressureTargetModel.h index 17af2627a3..0af6102e5b 100644 --- a/include/tudat/astro/electromagnetism/radiationPressureTargetModel.h +++ b/include/tudat/astro/electromagnetism/radiationPressureTargetModel.h @@ -51,17 +51,19 @@ class RadiationPressureTargetModel * Calculate radiation pressure force from incident radiation using geometrical/optical target properties. * * @param sourceIrradiance Incident irradiance magnitude [W/m²] - * @param sourceToTargetDirectionLocalFrame Direction of incoming radiation in local (i.e. target-fixed) coordinates + * @param sourceToTargetDirectionLocalFrame Direction of incoming radiation * @return Radiation pressure force vector in local (i.e. target-fixed) coordinates [N] */ virtual Eigen::Vector3d evaluateRadiationPressureForce( - const double sourceIrradiance, const Eigen::Vector3d& sourceToTargetDirectionLocalFrame) const = 0; + const double sourceIrradiance, const Eigen::Vector3d& sourceToTargetDirection) const = 0; std::map> getSourceToTargetOccultingBodies() const { return sourceToTargetOccultingBodies_; } + virtual bool forceFunctionRequiresLocalFrameInputs( ) = 0; + protected: virtual void updateMembers_(const double currentTime) {}; @@ -93,7 +95,7 @@ class CannonballRadiationPressureTargetModel : public RadiationPressureTargetMod Eigen::Vector3d evaluateRadiationPressureForce( double sourceIrradiance, - const Eigen::Vector3d& sourceToTargetDirectionLocalFrame) const override; + const Eigen::Vector3d& sourceToTargetDirection ) const override; double getArea() const { @@ -105,6 +107,11 @@ class CannonballRadiationPressureTargetModel : public RadiationPressureTargetMod return coefficient_; } + bool forceFunctionRequiresLocalFrameInputs( ) + { + return false; + } + private: double area_; double coefficient_; @@ -163,6 +170,11 @@ class PaneledRadiationPressureTargetModel : public RadiationPressureTargetModel return segmentFixedPanels_; } + bool forceFunctionRequiresLocalFrameInputs( ) + { + return true; + } + private: void updateMembers_( double currentTime ) override; diff --git a/include/tudat/astro/orbit_determination/acceleration_partials/radiationPressureAccelerationPartial.h b/include/tudat/astro/orbit_determination/acceleration_partials/radiationPressureAccelerationPartial.h index a262f404be..162b7c1b03 100644 --- a/include/tudat/astro/orbit_determination/acceleration_partials/radiationPressureAccelerationPartial.h +++ b/include/tudat/astro/orbit_determination/acceleration_partials/radiationPressureAccelerationPartial.h @@ -68,6 +68,24 @@ class CannonBallRadiationPressurePartial: public AccelerationPartial radiationPressureInterface ) ), acceleratedBodyMassFunction_( massFunction ){ } + CannonBallRadiationPressurePartial( + const std::shared_ptr< electromagnetism::CannonballRadiationPressureTargetModel > cannonballTargetModel, + const std::shared_ptr< electromagnetism::IsotropicPointSourceRadiationPressureAcceleration > sourceModel, + const std::string& acceleratedBody, const std::string& acceleratingBody ): + AccelerationPartial( acceleratedBody, acceleratingBody, + basic_astrodynamics::cannon_ball_radiation_pressure ), + sourceBodyState_( radiationPressureInterface->getSourcePositionFunction( ) ), + acceleratedBodyState_( radiationPressureInterface->getTargetPositionFunction( ) ), + areaFunction_( std::bind( &electromagnetism::CannonballRadiationPressureTargetModel::getArea, cannonballTargetModel ) ), + radiationPressureCoefficientFunction_( + std::bind( &electromagnetism::CannonballRadiationPressureTargetModel::getCoefficient, + cannonballTargetModel ) ), + radiationPressureFunction_( std::bind( &electromagnetism::RadiationPressureInterface::getCurrentRadiationPressure, + radiationPressureInterface ) ), + acceleratedBodyMassFunction_( std::bind( &electromagnetism::RadiationPressureAcceleration::getTargetMassFunction, sourceModel ) ){ } + { + + } //! Destructor. ~CannonBallRadiationPressurePartial( ){ } diff --git a/include/tudat/simulation/estimation_setup/createAccelerationPartials.h b/include/tudat/simulation/estimation_setup/createAccelerationPartials.h index b40065f55b..6b70d57b4d 100644 --- a/include/tudat/simulation/estimation_setup/createAccelerationPartials.h +++ b/include/tudat/simulation/estimation_setup/createAccelerationPartials.h @@ -177,6 +177,30 @@ std::vector< std::shared_ptr< orbit_determination::TidalLoveNumberPartialInterfa const SystemOfBodies& bodies, const std::string& acceleratingBodyName ); + +template< typename InitialStateParameterType = double > +std::shared_ptr< acceleration_partials::AccelerationPartial > createRadiationPressureAccelerationModel( + std::shared_ptr< electromagnetism::RadiationPressureAcceleration > radiationPressureAccelerationModel, + const std::pair< std::string, std::shared_ptr< simulation_setup::Body > > acceleratedBody, + const std::pair< std::string, std::shared_ptr< simulation_setup::Body > > acceleratingBody, + const simulation_setup::SystemOfBodies& bodies, + const std::shared_ptr< estimatable_parameters::EstimatableParameterSet< InitialStateParameterType > > + parametersToEstimate = + std::shared_ptr< estimatable_parameters::EstimatableParameterSet< InitialStateParameterType > >( ) ) +{ + using namespace electromagnetism; + std::shared_ptr< acceleration_partials::AccelerationPartial > accelerationPartial = nullptr; + + if( std::dynamic_pointer_cast< CannonballRadiationPressureTargetModel >( radiationPressureAccelerationModel->getTargetModel( ) ) != nullptr && + std::dynamic_pointer_cast< IsotropicPointSourceRadiationPressureAcceleration >( radiationPressureAccelerationModel ) != nullptr ) + { + accelerationPartial = std::make_shared< acceleration_partials::CannonBallRadiationPressurePartial > + ( std::dynamic_pointer_cast< CannonballRadiationPressureTargetModel >( radiationPressureAccelerationModel->getTargetModel( ) ), + std::dynamic_pointer_cast< IsotropicPointSourceRadiationPressureAcceleration >( radiationPressureAccelerationModel ), + acceleratedBody.first, acceleratingBody.first ); + } +} + //! Function to create a single acceleration partial derivative object. /*! * Function to create a single acceleration partial derivative object. @@ -216,9 +240,11 @@ std::shared_ptr< acceleration_partials::AccelerationPartial > createAnalyticalAc customPartialCalculator = createCustomAccelerationPartial( parametersToEstimate, accelerationModel, acceleratedBody, acceleratingBody, bodies ); - if( ( customPartialCalculator->getNumberOfCustomPartials( ) > 0 ) && ( accelerationType != custom_acceleration ) ) + if( ( customPartialCalculator->getNumberOfCustomPartials( ) > 0 ) && + ( accelerationType != custom_acceleration ) && + ( accelerationType != custom_acceleration ) ) { - throw std::runtime_error( "Error, custom acceleration partials only supported for custom acceleration (at present)" ); + throw std::runtime_error( "Error, custom acceleration partials only supported for custom and radiation pressure acceleration (at present)" ); } } switch( accelerationType ) @@ -702,6 +728,25 @@ std::shared_ptr< acceleration_partials::AccelerationPartial > createAnalyticalAc } break; } + case radiation_pressure: + { + // Check if identifier is consistent with type. + std::shared_ptr< RadiationPressureAcceleration > radiationPressureAccelerationModel = + std::dynamic_pointer_cast< RadiationPressureAcceleration >( accelerationModel ); + if( radiationPressureAccelerationModel == nullptr ) + { + throw std::runtime_error( + "Acceleration class type does not match acceleration type enum (radiation_pressure) set when making " + "acceleration partial." ); + } + else + { + // Create partial-calculating object. + accelerationPartial = createRadiationPressureAccelerationModel( ); + + } + break; + } case thrust_acceleration: { // Check if identifier is consistent with type. diff --git a/src/astro/electromagnetism/radiationPressureAcceleration.cpp b/src/astro/electromagnetism/radiationPressureAcceleration.cpp index 0b1a66ccc1..4664321b57 100644 --- a/src/astro/electromagnetism/radiationPressureAcceleration.cpp +++ b/src/astro/electromagnetism/radiationPressureAcceleration.cpp @@ -35,37 +35,39 @@ void RadiationPressureAcceleration::updateMembers(const double currentTime) Eigen::Vector3d IsotropicPointSourceRadiationPressureAcceleration::calculateAcceleration() { - Eigen::Vector3d sourceCenterPositionInGlobalFrame = sourcePositionFunction_(); // position of center of source (e.g. planet) - - Eigen::Vector3d targetCenterPositionInGlobalFrame = targetPositionFunction_(); - Eigen::Quaterniond targetRotationFromLocalToGlobalFrame = targetRotationFromLocalToGlobalFrameFunction_(); - Eigen::Quaterniond targetRotationFromGlobalToLocalFrame = targetRotationFromLocalToGlobalFrame.inverse(); + sourceCenterPositionInGlobalFrame_ = sourcePositionFunction_(); + targetCenterPositionInGlobalFrame_ = targetPositionFunction_(); + targetCenterPositionInSourceFrame_ = targetCenterPositionInGlobalFrame_ - sourceCenterPositionInGlobalFrame_; // Evaluate irradiances at target position in source frame // No rotation to source frame is necessary because isotropic sources are rotation-invariant - Eigen::Vector3d targetCenterPositionInSourceFrame = targetCenterPositionInGlobalFrame - sourceCenterPositionInGlobalFrame; sourceToTargetReceivedFraction = sourceToTargetOccultationModel_->evaluateReceivedFractionFromExtendedSource( - sourceCenterPositionInGlobalFrame, sourceBodyShapeModel_, targetCenterPositionInGlobalFrame); - auto sourceIrradiance = sourceModel_->evaluateIrradianceAtPosition(targetCenterPositionInSourceFrame).front().first; - auto occultedSourceIrradiance = sourceIrradiance * sourceToTargetReceivedFraction; - - // Update dependent variable - receivedIrradiance = occultedSourceIrradiance; + sourceCenterPositionInGlobalFrame_, sourceBodyShapeModel_, targetCenterPositionInGlobalFrame_ ); + receivedIrradiance = + sourceModel_->evaluateIrradianceAtPosition( targetCenterPositionInSourceFrame_).front().first * sourceToTargetReceivedFraction; - if (occultedSourceIrradiance <= 0) + if (receivedIrradiance <= 0) { // Some body is occluding source as seen from target return Eigen::Vector3d::Zero(); } - Eigen::Vector3d sourceToTargetDirectionInTargetFrame = - targetRotationFromGlobalToLocalFrame * (targetCenterPositionInGlobalFrame - sourceCenterPositionInGlobalFrame).normalized(); - // Calculate radiation pressure force due to source - Eigen::Vector3d totalForceInTargetFrame = - targetModel_->evaluateRadiationPressureForce(occultedSourceIrradiance, sourceToTargetDirectionInTargetFrame); - // Calculate acceleration due to radiation pressure in global frame - Eigen::Vector3d acceleration = targetRotationFromLocalToGlobalFrame * totalForceInTargetFrame / targetMassFunction_(); - return acceleration; + if( targetModel_->forceFunctionRequiresLocalFrameInputs( ) ) + { + Eigen::Quaterniond targetRotationFromLocalToGlobalFrame = targetRotationFromLocalToGlobalFrameFunction_(); + Eigen::Quaterniond targetRotationFromGlobalToLocalFrame = targetRotationFromLocalToGlobalFrame.inverse(); + + // Calculate acceleration due to radiation pressure in global frame + return targetRotationFromLocalToGlobalFrame * + targetModel_->evaluateRadiationPressureForce( + receivedIrradiance, targetRotationFromGlobalToLocalFrame * targetCenterPositionInSourceFrame_.normalized() ) / + targetMassFunction_(); + } + else + { + return targetModel_->evaluateRadiationPressureForce( + receivedIrradiance, targetCenterPositionInSourceFrame_.normalized() ) / targetMassFunction_(); + } } Eigen::Vector3d PaneledSourceRadiationPressureAcceleration::calculateAcceleration() diff --git a/src/astro/electromagnetism/radiationPressureTargetModel.cpp b/src/astro/electromagnetism/radiationPressureTargetModel.cpp index 33a09f2dc3..d86212f7b6 100644 --- a/src/astro/electromagnetism/radiationPressureTargetModel.cpp +++ b/src/astro/electromagnetism/radiationPressureTargetModel.cpp @@ -31,12 +31,12 @@ void RadiationPressureTargetModel::updateMembers(const double currentTime) Eigen::Vector3d CannonballRadiationPressureTargetModel::evaluateRadiationPressureForce( const double sourceIrradiance, - const Eigen::Vector3d& sourceToTargetDirectionLocalFrame) const + const Eigen::Vector3d& sourceToTargetDirection) const { // From Montenbruck (2000), Sec. 3.4 const auto radiationPressure = sourceIrradiance / physical_constants::SPEED_OF_LIGHT; const auto forceMagnitude = coefficient_ * area_ * radiationPressure; - Eigen::Vector3d force = forceMagnitude * sourceToTargetDirectionLocalFrame; + Eigen::Vector3d force = forceMagnitude * sourceToTargetDirection; return force; }