diff --git a/docs/source/Support/bskKnownIssues.rst b/docs/source/Support/bskKnownIssues.rst index 311b0da7f4..c24e85f797 100644 --- a/docs/source/Support/bskKnownIssues.rst +++ b/docs/source/Support/bskKnownIssues.rst @@ -17,6 +17,8 @@ Version |release| install script is updated this is corrected in the current release. - We found a slow memory leak if messages with arrays or vectors were accessed from python. The ``swig`` issue has now been fixed in the current release. +- The :ref:`facetSRPDynamicEffector` module was double counting a cosine term in the SRP force calculation. This is + corrected in the current release. Version 2.2.0 ------------- diff --git a/src/simulation/dynamics/facetSRPDynamicEffector/_UnitTest/test_unitFacetSRPDynamicEffector.py b/src/simulation/dynamics/facetSRPDynamicEffector/_UnitTest/test_unitFacetSRPDynamicEffector.py index b6eb725f05..ddc2cb435a 100644 --- a/src/simulation/dynamics/facetSRPDynamicEffector/_UnitTest/test_unitFacetSRPDynamicEffector.py +++ b/src/simulation/dynamics/facetSRPDynamicEffector/_UnitTest/test_unitFacetSRPDynamicEffector.py @@ -107,12 +107,9 @@ def checkFacetSRPForce(index, area, specCoeff, diffCoeff, normal_B, sigma_BN, sc # Calculate the incidence angle theta between the facet normal vector and the Sun-direction vector cosTheta = np.dot(sHat, normal_B) - intermediate = np.cross(sHat, normal_B) - sinTheta = np.linalg.norm(intermediate) - theta = np.arctan2(sinTheta, cosTheta) # Calculate the facet projected area onto the plane whose normal vector is the Sun-direction vector - projArea = area * np.cos(theta) + projArea = area * cosTheta # Calculate the solar radiation pressure acting on the facet numAU = AstU / np.linalg.norm(r_SB_B) @@ -120,7 +117,7 @@ def checkFacetSRPForce(index, area, specCoeff, diffCoeff, normal_B, sigma_BN, sc # Compute the SRP force acting on the facet only if the facet is illuminated by the Sun if projArea > 0: - srp_force = -SRPPressure * projArea * np.cos(theta) * ( (1-specCoeff) * sHat + 2 * ( (diffCoeff / 3) + specCoeff * np.cos(theta)) * normal_B ) + srp_force = -SRPPressure * projArea * ( (1-specCoeff) * sHat + 2 * ( (diffCoeff / 3) + specCoeff * cosTheta) * normal_B ) else: srp_force = np.zeros([3,]) @@ -146,12 +143,9 @@ def checkFacetSRPTorque(index, area, specCoeff, diffCoeff, normal_B, locationPnt # Calculate the incidence angle theta between the facet normal vector and the Sun-direction vector cosTheta = np.dot(sHat, normal_B) - intermediate = np.cross(sHat, normal_B) - sinTheta = np.linalg.norm(intermediate) - theta = np.arctan2(sinTheta, cosTheta) # Calculate the facet projected area onto the plane whose normal vector is the Sun-direction vector - projArea = area * np.cos(theta) + projArea = area * cosTheta # Calculate the solar radiation pressure acting on the facet numAU = AstU / np.linalg.norm(r_SB_B) @@ -159,7 +153,7 @@ def checkFacetSRPTorque(index, area, specCoeff, diffCoeff, normal_B, locationPnt # Compute the SRP force contribution from the facet only if the facet is illuminated by the Sun if projArea > 0: - srp_force = -SRPPressure * projArea * np.cos(theta) * ( (1-specCoeff) * sHat + 2 * ( (diffCoeff / 3) + specCoeff * np.cos(theta)) * normal_B ) + srp_force = -SRPPressure * projArea * ( (1-specCoeff) * sHat + 2 * ( (diffCoeff / 3) + specCoeff * cosTheta) * normal_B ) else: srp_force = np.zeros([3, ]) diff --git a/src/simulation/dynamics/facetSRPDynamicEffector/facetSRPDynamicEffector.cpp b/src/simulation/dynamics/facetSRPDynamicEffector/facetSRPDynamicEffector.cpp index 31eb3f763a..fcc2c5b3f7 100644 --- a/src/simulation/dynamics/facetSRPDynamicEffector/facetSRPDynamicEffector.cpp +++ b/src/simulation/dynamics/facetSRPDynamicEffector/facetSRPDynamicEffector.cpp @@ -126,7 +126,7 @@ void FacetSRPDynamicEffector::computeForceTorque(double integTime, double timeSt // Zero storage information double projectedArea = 0.0; - double projectionTerm = 0.0; + double cosTheta = 0.0; facetSRPForcePntB_B.setZero(); facetSRPTorquePntB_B.setZero(); totalSRPForcePntB_B.setZero(); @@ -141,22 +141,15 @@ void FacetSRPDynamicEffector::computeForceTorque(double integTime, double timeSt // Loop through the facets and calculate the SRP force and torque acting on point B for(int i = 0; i < this->numFacets; i++) { - projectionTerm = this->scGeometry.facetNormals_B[i].dot(sHat); - projectedArea = this->scGeometry.facetAreas[i] * projectionTerm; + cosTheta = this->scGeometry.facetNormals_B[i].dot(sHat); + projectedArea = this->scGeometry.facetAreas[i] * cosTheta; if(projectedArea > 0.0){ - - // Calculate the incidence angle theta between the facet normal vector and the Sun-direction vector - double cosTheta = projectionTerm; - Eigen::Vector3d intermediate = sHat.cross(this->scGeometry.facetNormals_B[i]); - double sinTheta = intermediate.norm(); - double theta = atan2(sinTheta, cosTheta); - // Compute the SRP force acting on the ith facet - facetSRPForcePntB_B = -SRPPressure * projectedArea * cos(theta) + facetSRPForcePntB_B = -SRPPressure * projectedArea * ( (1-this->scGeometry.facetSpecCoeffs[i]) * sHat + 2 * ( (this->scGeometry.facetDiffCoeffs[i] / 3) - + this->scGeometry.facetSpecCoeffs[i] * cos(theta)) + + this->scGeometry.facetSpecCoeffs[i] * cosTheta) * this->scGeometry.facetNormals_B[i] ); // Compute the SRP torque acting on the ith facet