Skip to content

Commit

Permalink
Bug fix for issue #553
Browse files Browse the repository at this point in the history
Improper use of the local frame rotation rate leads to the divergence vehicle rotation rates.
  • Loading branch information
bcoconni committed Jan 16, 2022
1 parent de2c2e8 commit 7115d78
Showing 1 changed file with 10 additions and 19 deletions.
29 changes: 10 additions & 19 deletions src/initialization/FGInitialCondition.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1187,15 +1187,7 @@ bool FGInitialCondition::Load_v1(Element* document)
if (document->FindElement("trim"))
SetTrimRequest(document->FindElementValue("trim"));

// Refer to Stevens and Lewis, 1.5-14a, pg. 49.
// This is the rotation rate of the "Local" frame, expressed in the local frame.
const FGMatrix33& Tl2b = orientation.GetT();
double radInv = 1.0 / position.GetRadius();
FGColumnVector3 vOmegaLocal = {radInv*vUVW_NED(eEast),
-radInv*vUVW_NED(eNorth),
-radInv*vUVW_NED(eEast)*tan(position.GetLatitude())};

vPQR_body = Tl2b * vOmegaLocal;
vPQR_body.InitMatrix();

return result;
}
Expand Down Expand Up @@ -1408,19 +1400,12 @@ bool FGInitialCondition::Load_v2(Element* document)
// - Body

Element* attrate_el = document->FindElement("attitude_rate");
const FGMatrix33& Tl2b = orientation.GetT();

// Refer to Stevens and Lewis, 1.5-14a, pg. 49.
// This is the rotation rate of the "Local" frame, expressed in the local frame.
double radInv = 1.0 / position.GetRadius();
FGColumnVector3 vOmegaLocal = { radInv*vUVW_NED(eEast),
-radInv*vUVW_NED(eNorth),
-radInv*vUVW_NED(eEast)*tan(position.GetLatitude())};

if (attrate_el) {

string frame = attrate_el->GetAttributeValue("frame");
frame = to_lower(frame);
const FGMatrix33& Tl2b = orientation.GetT();
FGColumnVector3 vAttRate = attrate_el->FindElementTripletConvertTo("RAD/SEC");

if (frame == "eci") {
Expand All @@ -1429,6 +1414,12 @@ bool FGInitialCondition::Load_v2(Element* document)
} else if (frame == "ecef") {
vPQR_body = Tl2b * position.GetTec2l() * vAttRate;
} else if (frame == "local") {
// Refer to Stevens and Lewis, 1.5-14a, pg. 49.
// This is the rotation rate of the "Local" frame, expressed in the local frame.
double radInv = 1.0 / position.GetRadius();
FGColumnVector3 vOmegaLocal = {radInv*vUVW_NED(eEast),
-radInv*vUVW_NED(eNorth),
-radInv*vUVW_NED(eEast)*tan(position.GetLatitude())};

This comment has been minimized.

Copy link
@seanmcleod

seanmcleod Jan 17, 2022

Member

But we'll still have the issue of tan(90) blowing up in this case right?

This comment has been minimized.

Copy link
@bcoconni

bcoconni Jan 18, 2022

Author Member

You are correct: in that particular case, the rotational rate can still blow up when evaluated at the Poles. Given that in that particular case (frame == "local"), the rotational rate AttRate is expressed with respect to the local frame, we have no choice but to use the Stevens & Lewis formula. However we might include some safeguards to this formula. Either by capping the max value vOmegaLocal components, by issuing a warning message or maybe throwing an exception ?

vPQR_body = Tl2b * (vAttRate + vOmegaLocal);
} else if (frame == "body") {
vPQR_body = vAttRate;
Expand All @@ -1439,11 +1430,11 @@ bool FGInitialCondition::Load_v2(Element* document)
result = false;

} else if (frame.empty()) {
vPQR_body = Tl2b * vOmegaLocal;
vPQR_body.InitMatrix();
}

} else { // Body frame attitude rate assumed 0 relative to local.
vPQR_body = Tl2b * vOmegaLocal;
vPQR_body.InitMatrix();
}

return result;
Expand Down

0 comments on commit 7115d78

Please sign in to comment.