Working on radiation pressure partials

include/tudat/astro/electromagnetism/radiationPressureAcceleration.h

@@ -59,6 +59,11 @@ class RadiationPressureAcceleration: public basic_astrodynamics::AccelerationMod
         return receivedIrradiance;
     }
 
+    std::function<double()>  getTargetMassFunction()
+    {
+        return targetMassFunction_;
+    }
+
 protected:
     RadiationPressureAcceleration(const std::function<Eigen::Vector3d()>& sourcePositionFunction,
                                   const std::shared_ptr<RadiationPressureTargetModel>& targetModel,
@@ -140,8 +145,13 @@ class IsotropicPointSourceRadiationPressureAcceleration : public RadiationPressu
     Eigen::Vector3d calculateAcceleration() override;
 
     std::shared_ptr<IsotropicPointRadiationSourceModel> sourceModel_;
+
     std::shared_ptr<basic_astrodynamics::BodyShapeModel> sourceBodyShapeModel_;
 
+    Eigen::Vector3d sourceCenterPositionInGlobalFrame_;
+    Eigen::Vector3d targetCenterPositionInGlobalFrame_;
+    Eigen::Vector3d targetCenterPositionInSourceFrame_;
+
     // For dependent variable
     double sourceToTargetReceivedFraction;
 };

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<std::string, std::vector<std::string>> 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;
 

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( ){ }
 

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.

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()

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;
 }