|
| 1 | +/* |
| 2 | + ISC License |
| 3 | +
|
| 4 | + Copyright (c) 2026, Autonomous Vehicle Systems Lab, University of Colorado at Boulder |
| 5 | +
|
| 6 | + Permission to use, copy, modify, and/or distribute this software for any |
| 7 | + purpose with or without fee is hereby granted, provided that the above |
| 8 | + copyright notice and this permission notice appear in all copies. |
| 9 | +
|
| 10 | + THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES |
| 11 | + WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF |
| 12 | + MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR |
| 13 | + ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES |
| 14 | + WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN |
| 15 | + ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF |
| 16 | + OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. |
| 17 | +
|
| 18 | +*/ |
| 19 | + |
| 20 | +#include "facetedSRPEffector.h" |
| 21 | +#include "architecture/utilities/avsEigenSupport.h" |
| 22 | + |
| 23 | +const double speedLight = 299792458.0; // [m/s] Speed of light |
| 24 | +const double AstU = 149597870700.0; // [m] Astronomical unit |
| 25 | +const double solarRadFlux = 1368.0; // [W/m^2] Solar radiation flux at 1 AU |
| 26 | + |
| 27 | +/*! This method resets required module variables and checks the input messages to ensure they are linked. |
| 28 | +
|
| 29 | + @param currentSimNanos [ns] Time the method is called |
| 30 | +*/ |
| 31 | +void FacetedSRPEffector::Reset(uint64_t currentSimNanos) { |
| 32 | + // Check Sun state input message is linked |
| 33 | + if (!this->sunStateInMsg.isLinked()) { |
| 34 | + bskLogger.bskLog(BSK_ERROR, "FacetedSRPEffector.sunStateInMsg was not linked."); |
| 35 | + return; |
| 36 | + } |
| 37 | + |
| 38 | + // Check all facet input messages are linked |
| 39 | + for (uint64_t idx = 0; idx < this->facetElementBodyInMsgs.size(); idx++) { |
| 40 | + if (!this->facetElementBodyInMsgs[idx].isLinked()) { |
| 41 | + this->bskLogger.bskLog(BSK_ERROR, |
| 42 | + "FacetedSpacecraftProjectedArea: Input message error. facetElementBodyInMsgs is not linked."); |
| 43 | + return; |
| 44 | + } |
| 45 | + if (!this->facetProjectedAreaInMsgs[idx].isLinked()) { |
| 46 | + this->bskLogger.bskLog(BSK_ERROR, |
| 47 | + "FacetedSpacecraftProjectedArea: Input message error. facetProjectedAreaInMsgs is not linked."); |
| 48 | + return; |
| 49 | + } |
| 50 | + } |
| 51 | + |
| 52 | + // Clear, allocate and initialize data lists |
| 53 | + this->facetAreaList.clear(); |
| 54 | + this->facetProjectedAreaList.clear(); |
| 55 | + this->facetR_CopB_BList.clear(); |
| 56 | + this->facetNHat_BList.clear(); |
| 57 | + this->facetDiffuseCoeffList.clear(); |
| 58 | + this->facetSpecularCoeffList.clear(); |
| 59 | + this->facetAreaList.reserve(this->numFacets); |
| 60 | + this->facetProjectedAreaList.reserve(this->numFacets); |
| 61 | + this->facetR_CopB_BList.reserve(this->numFacets); |
| 62 | + this->facetNHat_BList.reserve(this->numFacets); |
| 63 | + this->facetDiffuseCoeffList.reserve(this->numFacets); |
| 64 | + this->facetSpecularCoeffList.reserve(this->numFacets); |
| 65 | + for (uint64_t idx = 0; idx < this->numFacets; idx++) { |
| 66 | + this->facetAreaList.push_back(0.0); |
| 67 | + this->facetProjectedAreaList.push_back(0.0); |
| 68 | + this->facetR_CopB_BList.push_back(Eigen::Vector3d::Zero()); |
| 69 | + this->facetNHat_BList.push_back(Eigen::Vector3d::Zero()); |
| 70 | + this->facetDiffuseCoeffList.push_back(0.0); |
| 71 | + this->facetSpecularCoeffList.push_back(0.0); |
| 72 | + } |
| 73 | +} |
| 74 | + |
| 75 | +/*! This method gives the module access to the hub inertial attitude and position. |
| 76 | +
|
| 77 | + @param states Dynamic parameter states |
| 78 | +*/ |
| 79 | +void FacetedSRPEffector::linkInStates(DynParamManager& states) { |
| 80 | + this->hubSigma = states.getStateObject(this->stateNameOfSigma); |
| 81 | + this->hubPosition = states.getStateObject(this->stateNameOfPosition); |
| 82 | +} |
| 83 | + |
| 84 | +/*! This method computes the srp force and torque acting about the hub point B in B frame components. |
| 85 | +
|
| 86 | + @param callTime [s] Time the method is called |
| 87 | + @param timeStep [s] Simulation time step |
| 88 | +*/ |
| 89 | +void FacetedSRPEffector::computeForceTorque(double callTime, double timeStep) { |
| 90 | + |
| 91 | + this->forceExternal_B.setZero(); |
| 92 | + this->torqueExternalPntB_B.setZero(); |
| 93 | + if (this->facetAreaList.size() != this->numFacets) { |
| 94 | + return; |
| 95 | + } |
| 96 | + |
| 97 | + // Read the sun state input message |
| 98 | + SpicePlanetStateMsgPayload sunStateIn{}; |
| 99 | + Eigen::Vector3d r_SN_N = Eigen::Vector3d::Zero(); |
| 100 | + if (this->sunStateInMsg.isWritten()) { |
| 101 | + sunStateIn = this->sunStateInMsg(); |
| 102 | + r_SN_N = cArray2EigenVector3d(sunStateIn.PositionVector); |
| 103 | + } |
| 104 | + |
| 105 | + // Compute the sun direction unit vector in spacecraft body frame components |
| 106 | + Eigen::MRPd sigma_BN; |
| 107 | + sigma_BN = (Eigen::Vector3d)this->hubSigma->getState(); |
| 108 | + Eigen::Matrix3d dcm_BN = sigma_BN.toRotationMatrix().transpose(); |
| 109 | + Eigen::Vector3d r_BN_N = this->hubPosition->getState(); |
| 110 | + Eigen::Vector3d r_SB_B = dcm_BN * (r_SN_N - r_BN_N); |
| 111 | + Eigen::Vector3d sunDirHat_B = r_SB_B.normalized(); |
| 112 | + |
| 113 | + // Read the facet geometry and projected area input messages |
| 114 | + FacetElementBodyMsgPayload facetElementBodyIn{}; |
| 115 | + ProjectedAreaMsgPayload facetProjectedAreaIn{}; |
| 116 | + for (uint64_t idx = 0; idx < this->numFacets; idx++) { |
| 117 | + if (this->facetElementBodyInMsgs[idx].isWritten() && this->facetProjectedAreaInMsgs[idx].isWritten()) { |
| 118 | + facetElementBodyIn = this->facetElementBodyInMsgs[idx](); |
| 119 | + this->facetAreaList[idx] = facetElementBodyIn.area; |
| 120 | + this->facetR_CopB_BList[idx] = cArray2EigenVector3d(facetElementBodyIn.r_CopB_B); |
| 121 | + this->facetNHat_BList[idx] = cArray2EigenVector3d(facetElementBodyIn.nHat_B); |
| 122 | + this->facetDiffuseCoeffList[idx] = facetElementBodyIn.c_diffuse; |
| 123 | + this->facetSpecularCoeffList[idx] = facetElementBodyIn.c_specular; |
| 124 | + |
| 125 | + facetProjectedAreaIn = this->facetProjectedAreaInMsgs[idx](); |
| 126 | + this->facetProjectedAreaList[idx] = facetProjectedAreaIn.area; |
| 127 | + } |
| 128 | + } |
| 129 | + |
| 130 | + // Calculate the SRP pressure acting at the current spacecraft location |
| 131 | + double numAU = AstU / r_SB_B.norm(); |
| 132 | + double SRPPressure = (solarRadFlux / speedLight) * numAU * numAU; // [Pa] |
| 133 | + |
| 134 | + // Compute the per-facet and total SRP force and torque |
| 135 | + Eigen::Vector3d facetSRPForcePntB_B = Eigen::Vector3d::Zero(); |
| 136 | + Eigen::Vector3d facetSRPTorquePntB_B = Eigen::Vector3d::Zero(); |
| 137 | + Eigen::Vector3d totalSRPForcePntB_B = Eigen::Vector3d::Zero(); |
| 138 | + Eigen::Vector3d totalSRPTorquePntB_B = Eigen::Vector3d::Zero(); |
| 139 | + |
| 140 | + for (uint64_t idx = 0; idx < this->numFacets; idx++) { |
| 141 | + if (this->facetProjectedAreaList[idx] > 0.0) { |
| 142 | + double cosTheta = facetProjectedAreaList[idx] / facetAreaList[idx]; |
| 143 | + |
| 144 | + facetSRPForcePntB_B = -SRPPressure * this->facetProjectedAreaList[idx] |
| 145 | + * ((1 - this->facetSpecularCoeffList[idx]) |
| 146 | + * sunDirHat_B + 2 * ((this->facetDiffuseCoeffList[idx] / 3) |
| 147 | + + this->facetSpecularCoeffList[idx] * cosTheta) |
| 148 | + * this->facetNHat_BList[idx]); // [N] |
| 149 | + facetSRPTorquePntB_B = this->facetR_CopB_BList[idx].cross(facetSRPForcePntB_B); // [Nm] |
| 150 | + |
| 151 | + // Add the facet contribution to the total SRP force and torque acting on the spacecraft |
| 152 | + totalSRPForcePntB_B += facetSRPForcePntB_B; |
| 153 | + totalSRPTorquePntB_B += facetSRPTorquePntB_B; |
| 154 | + } |
| 155 | + } |
| 156 | + |
| 157 | + // Update the force and torque vectors in the dynamic effector base class |
| 158 | + this->forceExternal_B = totalSRPForcePntB_B; |
| 159 | + this->torqueExternalPntB_B = totalSRPTorquePntB_B; |
| 160 | +} |
| 161 | + |
| 162 | +/*! Setter method for the total number of spacecraft facets. |
| 163 | + @param numFacets [-] |
| 164 | +*/ |
| 165 | +void FacetedSRPEffector::setNumFacets(const uint64_t numFacets) { |
| 166 | + this->numFacets = numFacets; |
| 167 | + |
| 168 | + this->facetElementBodyInMsgs.clear(); |
| 169 | + this->facetProjectedAreaInMsgs.clear(); |
| 170 | + this->facetElementBodyInMsgs.reserve(this->numFacets); |
| 171 | + this->facetProjectedAreaInMsgs.reserve(this->numFacets); |
| 172 | + |
| 173 | + // Push back facet message vectors |
| 174 | + for (uint64_t idx = 0; idx < this->numFacets; ++idx) { |
| 175 | + this->facetElementBodyInMsgs.push_back(ReadFunctor<FacetElementBodyMsgPayload>{}); |
| 176 | + this->facetProjectedAreaInMsgs.push_back(ReadFunctor<ProjectedAreaMsgPayload>()); |
| 177 | + } |
| 178 | +} |
| 179 | + |
| 180 | +/*! Getter method for the total number of spacecraft facets. |
| 181 | + @return const uint64_t |
| 182 | +*/ |
| 183 | +const uint64_t FacetedSRPEffector::getNumFacets() const { return this->numFacets; } |
0 commit comments