Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix Bug in SRP Dynamic Effector Module #472

Merged
merged 2 commits into from
Oct 30, 2023
Merged

Fix Bug in SRP Dynamic Effector Module #472

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
46351ca
Remove double count of cos(theta) term
leahkiner Oct 30, 2023
5268c2d
Add the bug fix to the known issues file
leahkiner Oct 30, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 4 additions & 10 deletions ...simulation/dynamics/facetSRPDynamicEffector/_UnitTest/test_unitFacetSRPDynamicEffector.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -107,20 +107,17 @@ 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)
SRPPressure = (solarRadFlux / speedLight) * numAU * numAU

# 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,])

Expand All @@ -146,20 +143,17 @@ 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)
SRPPressure = (solarRadFlux / speedLight) * numAU * numAU

# 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, ])

Expand Down
17 changes: 5 additions & 12 deletions src/simulation/dynamics/facetSRPDynamicEffector/facetSRPDynamicEffector.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand All @@ -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
Expand Down
Loading