Radiation pressure partial correction

tudat-team · Oct 22, 2024 · 5243b11 · 5243b11
1 parent 61c2a19
commit 5243b11
Show file tree

Hide file tree

Showing 4 changed files with 71 additions and 52 deletions.
diff --git a/include/tudat/astro/electromagnetism/radiationPressureTargetModel.h b/include/tudat/astro/electromagnetism/radiationPressureTargetModel.h
@@ -516,7 +516,7 @@ class PaneledRadiationPressureTargetModel : public RadiationPressureTargetModel
         {
             panelTorques_.at( i ).setZero( );
         }
-        panelForcesPerSource_[ sourceName ] = panelTorques_;
+        panelTorquesPerSource_[ sourceName ] = panelTorques_;
     }
 
 

diff --git a/src/astro/electromagnetism/radiationPressureTargetModel.cpp b/src/astro/electromagnetism/radiationPressureTargetModel.cpp
@@ -195,7 +195,6 @@ Eigen::Vector3d PaneledRadiationPressureTargetModel::evaluateRadiationPressureFo
     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++ )
     {

diff --git a/...rbit_determination/acceleration_partials/panelledRadiationPressureAccelerationPartial.cpp b/...rbit_determination/acceleration_partials/panelledRadiationPressureAccelerationPartial.cpp
@@ -33,6 +33,9 @@ void PanelledRadiationPressurePartial::update( const double currentTime )
                     2.0 * currentRadiationPressure * bodyFixedUnitVectorToSource.transpose( ) / ( distanceToSource );
 
             currentCosineAnglePartial_ = Eigen::Matrix< double, 1, 3 >::Zero( );
+
+            Eigen::Matrix3d rotationToInertialFrame =
+                Eigen::Matrix3d( radiationPressureAcceleration_->getTargetRotationFromLocalToGlobalFrame( ) );
             Eigen::Vector3d currentPanelReactionVector = Eigen::Vector3d::Zero( );
             Eigen::Vector3d currentPanelNormal = Eigen::Vector3d::Zero( );
             Eigen::Matrix3d currentPanelPartialContribution = Eigen::Matrix3d::Zero( );
@@ -42,31 +45,39 @@ void PanelledRadiationPressurePartial::update( const double currentTime )
 
             for( int i = 0; i < panelledTargetModel_->getTotalNumberOfPanels( ); i++ )
             {
-                currentPanelNormal = panelledTargetModel_->getSurfaceNormals( ).at( i );
-                cosineOfPanelInclination = panelledTargetModel_->getSurfacePanelCosines( acceleratingBody_ ).at( i );
-
                 currentPanelPartialContribution.setZero( );
-                if( cosineOfPanelInclination > 0.0 )
+                if( panelledTargetModel_->getBodyFixedPanels( ).at( i )->getTrackedBody( ) != acceleratingBody_ )
                 {
-                    currentCosineAnglePartial_ = currentPanelNormal.transpose( ) * currentSourceUnitVectorPartial_;
-
-                    currentPanelArea = panelledTargetModel_->getBodyFixedPanels( ).at( i )->getPanelArea( );
-                    currentPanelReactionVector = panelledTargetModel_->getPanelForces( acceleratingBody_ ).at( i ) / ( currentRadiationPressure * currentPanelArea );
+                    currentPanelNormal = panelledTargetModel_->getSurfaceNormals( ).at( i );
+                    cosineOfPanelInclination = panelledTargetModel_->getSurfacePanelCosines( acceleratingBody_ ).at( i );
 
-                    currentPanelPartialContribution += panelledTargetModel_->getFullPanels( ).at( i )->getReflectionLaw( )->
-                        evaluateReactionVectorDerivativeWrtTargetPosition(
-                            currentPanelNormal, -bodyFixedUnitVectorToSource, cosineOfPanelInclination, currentPanelReactionVector,
-                            currentSourceUnitVectorPartial_, currentCosineAnglePartial_ );
+                    currentPanelPartialContribution.setZero( );
+                    if( cosineOfPanelInclination > 0.0 )
+                    {
+                        currentCosineAnglePartial_ = currentPanelNormal.transpose( ) * currentSourceUnitVectorPartial_;
 
-                    currentPanelPartialContribution *= currentRadiationPressure * currentPanelArea;
+                        currentPanelArea = panelledTargetModel_->getBodyFixedPanels( ).at( i )->getPanelArea( );
+                        currentPanelReactionVector = panelledTargetModel_->getPanelForces( acceleratingBody_ ).at( i ) / ( currentRadiationPressure * currentPanelArea );
 
+                        currentPanelPartialContribution += panelledTargetModel_->getFullPanels( ).at( i )->getReflectionLaw( )->
+                            evaluateReactionVectorDerivativeWrtTargetPosition(
+                                currentPanelNormal, -bodyFixedUnitVectorToSource, cosineOfPanelInclination, currentPanelReactionVector,
+                                currentSourceUnitVectorPartial_, currentCosineAnglePartial_ );
 
+                        currentPanelPartialContribution *= currentRadiationPressure * currentPanelArea;
+                    }
+                }
+                else
+                {
+                    currentPanelPartialContribution =
+                        panelledTargetModel_->getPanelForces( acceleratingBody_ ).at( i ).norm( ) *
+                        rotationToInertialFrame.transpose( ) *
+                        ( Eigen::Matrix3d::Identity( ) / distanceToSource - inertialVectorFromSource * inertialVectorFromSource.transpose( ) / ( distanceToSource * distanceToSource * distanceToSource ) ) *
+                        rotationToInertialFrame;
                 }
                 currentPartialWrtPosition_ += currentPanelPartialContribution;
             }
 
-            Eigen::Matrix3d rotationToInertialFrame =
-                Eigen::Matrix3d( radiationPressureAcceleration_->getTargetRotationFromLocalToGlobalFrame( ) );
             currentPartialWrtPosition_ = rotationToInertialFrame * currentPartialWrtPosition_ * rotationToInertialFrame.transpose( );
             currentPartialWrtPosition_ /= currentMass;
 

diff --git a/tests/src/astro/orbit_determination/acceleration_partials/unitTestAccelerationPartials.cpp b/tests/src/astro/orbit_determination/acceleration_partials/unitTestAccelerationPartials.cpp
@@ -74,7 +74,7 @@ BOOST_AUTO_TEST_CASE( testPanelledRadiationPressureAccelerationPartials )
 {
     tudat::spice_interface::loadStandardSpiceKernels( );
 
-    for( int testIndex = 0; testIndex < 6; testIndex++ )
+    for( int testIndex = 0; testIndex < 7; testIndex++ )
     {
         // Create empty bodies, earth and sun.
         std::shared_ptr< Body > vehicle = std::make_shared< Body >( );
@@ -114,43 +114,58 @@ BOOST_AUTO_TEST_CASE( testPanelledRadiationPressureAccelerationPartials )
         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 )
+
+        if( testIndex < 6 )
         {
+            if( testIndex % 2 == 1 )
+            {
+                useInstantaneousReradiation = true;
+            }
+            std::string panelName = "SolarPanel";
+            double scalingValue = 1.0;
+            if( testIndex > 3 )
+            {
+                scalingValue = 0.5;
+                panelName = "PanelName";
+            }
             panelSettingsList.push_back( std::make_shared< BodyPanelSettings >(
-                std::make_shared< FrameFixedBodyPanelGeometrySettings >( -Eigen::Vector3d::UnitY( ), 3.254 ),
+                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 ) );
+            }
+        }
+        else
+        {
             panelSettingsList.push_back( std::make_shared< BodyPanelSettings >(
-                std::make_shared< FrameFixedBodyPanelGeometrySettings >( -Eigen::Vector3d::UnitZ( ), 8.654 ),
+                bodyTrackingPanelGeometry( "Sun", true , 1.0 ),
                 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 ) );
+                bodyTrackingPanelGeometry( "Sun", false , 1.346 ),
+                std::make_shared< SpecularDiffuseBodyPanelReflectionLawSettings >( specularReflectivity, diffuseReflectivity, useInstantaneousReradiation ), "SolarPanel" ) );
+
         }
 
+
         addBodyExteriorPanelledShape(
             std::make_shared< FullPanelledBodySettings >( panelSettingsList ), "Vehicle", bodies );
 
@@ -250,12 +265,6 @@ BOOST_AUTO_TEST_CASE( testPanelledRadiationPressureAccelerationPartials )
         testPartialWrtSpecularReflectivity = calculateAccelerationWrtParameterPartials(
             specularReflectivityParameter, accelerationModel, 0.1, updateFunction );
 
-        std::cout<<partialWrtSpecularReflectivity<<std::endl<<std::endl;
-        std::cout<<testPartialWrtSpecularReflectivity<<std::endl<<std::endl;
-
-        std::cout<<partialWrtDiffuseReflectivity<<std::endl<<std::endl;
-        std::cout<<testPartialWrtDiffuseReflectivity<<std::endl<<std::endl;
-
         // Compare numerical and analytical results.
         TUDAT_CHECK_MATRIX_CLOSE_FRACTION( testPartialWrtSunPosition,
                                            partialWrtSunPosition, 1.0e-6 );