Specular and diffuse

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"

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<unsigned int>(totalNumberOfPanels_); i++ )
+        {
+            if(fullPanels_.at(i)->getPanelTypeId() == panelTypeId)
+            {
+                panels.push_back(fullPanels_.at(i));
+            }
+        }
+        return panels;
+    }
+
+    double getAverageDiffuseReflectivity(const std::string& panelTypeId)
+    {
+        std::vector<double> 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<SpecularDiffuseMixReflectionLaw>(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<double> 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<SpecularDiffuseMixReflectionLaw>(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<double> getSpecularReflectivityForPanelTypeId(const std::string& panelTypeId)
+    {
+        std::vector<double> 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<SpecularDiffuseMixReflectionLaw>(panelsFromId.at(i)->getReflectionLaw());
+            coefficients.push_back(reflectionLaw->getSpecularReflectivity());
+        }
+
+        return coefficients;
+    }
+
+    std::vector<double> getDiffuseReflectivityForPanelTypeId(const std::string& panelTypeId)
+    {
+        std::vector<double> 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<SpecularDiffuseMixReflectionLaw>(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<SpecularDiffuseMixReflectionLaw>(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<SpecularDiffuseMixReflectionLaw>(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<SpecularDiffuseMixReflectionLaw>(fullPanels_.at(i)->getReflectionLaw());
+                reflectionLaw->setAbsorptivity(absorbtivity);
+            }
+        }
+    };
+
 private:
     void updateMembers_( double currentTime ) override;
 

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

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"

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<electromagnetism::PaneledRadiationPressureTargetModel>(
+                    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<electromagnetism::PaneledRadiationPressureTargetModel>(
+                                                     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<electromagnetism::PaneledRadiationPressureTargetModel>(
+                    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<electromagnetism::PaneledRadiationPressureTargetModel>(
+                                                     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_;

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