How to derive the formula in class FGInitialCondition

[emptyxiao]

How to derive the formula in function "calcThetaBeta"

jsbsim/src/initialization/FGInitialCondition.cpp

Lines 321 to 390 in 59c0bd9

	// When the AoA is modified, we need to update the angles theta and beta to
	// keep the true airspeed amplitude, the climb rate and the heading unchanged.
	// Beta will be modified if the aircraft roll angle is not null.
	void FGInitialCondition::calcThetaBeta(double alfa, const FGColumnVector3& _vt_NED)
	{
	FGColumnVector3 vOrient = orientation.GetEuler();
	double calpha = cos(alfa), salpha = sin(alfa);
	double cpsi = orientation.GetCosEuler(ePsi), spsi = orientation.GetSinEuler(ePsi);
	double cphi = orientation.GetCosEuler(ePhi), sphi = orientation.GetSinEuler(ePhi);
	FGMatrix33 Tpsi( cpsi, spsi, 0.,
	-spsi, cpsi, 0.,
	0., 0., 1.);
	FGMatrix33 Tphi(1., 0., 0.,
	0., cphi, sphi,
	0.,-sphi, cphi);
	FGMatrix33 Talpha( calpha, 0., salpha,
	0., 1., 0.,
	-salpha, 0., calpha);
	FGColumnVector3 v0 = Tpsi * _vt_NED;
	FGColumnVector3 n = (Talpha * Tphi).Transposed() * FGColumnVector3(0., 0., 1.);
	FGColumnVector3 y = {0., 1., 0.};
	FGColumnVector3 u = y - DotProduct(y, n) * n;
	FGColumnVector3 p = y * n;
	if (DotProduct(p, v0) < 0) p *= -1.0;
	p.Normalize();
	u *= DotProduct(v0, y) / DotProduct(u, y);
	// There are situations where the desired alpha angle cannot be obtained. This
	// is not a limitation of the algorithm but is due to the mathematical problem
	// not having a solution. This can only be cured by limiting the alpha angle
	// or by modifying an additional angle (psi ?). Since this is anticipated to
	// be a pathological case (mainly when a high roll angle is required) this
	// situation is not addressed below. However if there are complaints about the
	// following error being raised too often, we might need to reconsider this
	// position.
	if (DotProduct(v0, v0) < DotProduct(u, u)) {
	cerr << "Cannot modify angle 'alpha' from " << alpha << " to " << alfa << endl;
	return;
	}
	FGColumnVector3 v1 = u + sqrt(DotProduct(v0, v0) - DotProduct(u, u))*p;
	FGColumnVector3 v0xz(v0(eU), 0., v0(eW));
	FGColumnVector3 v1xz(v1(eU), 0., v1(eW));
	v0xz.Normalize();
	v1xz.Normalize();
	double sinTheta = (v1xz * v0xz)(eY);
	vOrient(eTht) = asin(sinTheta);
	orientation = FGQuaternion(vOrient);
	const FGMatrix33& Tl2b = orientation.GetT();
	FGColumnVector3 v2 = Talpha * Tl2b * _vt_NED;
	alpha = alfa;
	beta = atan2(v2(eV), v2(eU));
	double cbeta=1.0, sbeta=0.0;
	if (vt != 0.0) {
	cbeta = v2(eU) / vt;
	sbeta = v2(eV) / vt;
	}
	Tw2b = { calpha*cbeta, -calpha*sbeta, -salpha,
	sbeta, cbeta, 0.0,
	salpha*cbeta, -salpha*sbeta, calpha };
	Tb2w = Tw2b.Transposed();
	}

[seanmcleod]

So in the simple case assuming 0 bank angle ( $\phi$ ) and the supplied velocity vector $V$ and a requested AoA ( $\alpha$ ) the calculation is simple:

$$ \theta = \gamma + \alpha $$

However the calcThetaBeta(alpha, vt_NED) function needs to handle the more general case in terms a velocity vector ( $V$ ) supplied in the NED coordinate frame, a requested $\alpha$ and orientation/attitude constraints that have been set separately, i.e. the bank angle ( $\phi$ ) and heading angle ( $\psi$ ).

Given the constraints, in particular the orientation constraints this means that the function isn't always able to calculate a solution. For example say the bank angle is 90 degrees and the heading angle is aligned with the velocity vector's horizontal angle and a non-zero alpha is requested then it's impossible to achieve the requested alpha angle. In this case the only way to achieve the requested alpha angle would be to relax the heading angle constraint.

Here are some images drawn by @agodemar which shows all the relevant angles for the general case. Maybe they'll make it a bit clearer in terms of how the function goes about it's calculation.

https://jsbsim-team.github.io/jsbsim-reference-manual/mypages/user-manual-frames-of-reference/

[emptyxiao]

Thank you for your reply! But I also want to know the mathematical derivation process. It really bothers me.
Can you demonstrate the mathematical derivation process, or provide relevant mathematical derivation process papers?
Thank you very much!

[seanmcleod]

@emptyxiao are you looking for a detailed explanation, i.e. virtually line by line explaining the maths/geometry behind every line in the function, or do you follow say 80% of the code and you just need a detailed explanation for some 20% of the code, i.e. particular steps you don't follow?

I don't have a detailed line by line explanation handy, I would need to spend a fair bit of time going through the code.

Going through the commit history I've figured out when this particular function was added and by whom.

"Complete rewrite of the FGInitialCondition class" - da0b947, in 2010, method originally called SetAlphaRadIC().

The good news is the code author, @bcoconni is still around and very active in terms of JSBSim.

[agodemar]

@emptyxiao

But I also want to know the mathematical derivation process. It really bothers me.
Can you demonstrate the mathematical derivation process, or provide relevant mathematical derivation process papers?
Thank you very much!

I can feel your thirst for learning the source of these formulas.

I think that this paper De Marco A., Duke E. L., Berndt, J. S., "A General Solution to the Aircraft Trim Problem." AIAA Paper 2007-6703, AIAA Modeling and Simulation Technologies Conference and Exhibit 20-23 August 2007, Hilton Head, South Carolina, USA might be of some help to you, pages from 15 to 24.

You coud help us to double check if the formulas we use in code are correct when compared to those in the theory recalled in the paper.

[bcoconni]

Going through the commit history I've figured out when this particular function was added and by whom.

"Complete rewrite of the FGInitialCondition class" - da0b947, in 2010, method originally called SetAlphaRadIC().

The good news is the code author, @bcoconni is still around and very active in terms of JSBSim.

Oh man, that was many moons ago 🤔. I've found some notes I wrote at the time on this topic but there are not many details involved: I have most likely dumped the sheets on which I derived these equations. I'll need (quite) some time to wrap my head around these formulas and give a detailed explanation.

@emptyxiao, you might have found there are some errors in the results that these formulas give ? If that is so, I'd be interested in the cases you found that break the computations

[emptyxiao]

First of all, thank you for your reply！Recently, I was deriving the formulas in JSBSim. At present, this process stays at FGInitialCondition class .Learning FGTrim class is the next stage.
I have found a problem in function "UpdateCustomTransformMatrix()".I think it should use "mT.Transposed()"instead of "mT".I'm not sure I'm right.If you are also interesting in this problem.You may check this function.
@bcoconni If you can give me a detailed explanation about ‘’calcThetaBeta‘’,I would appreciate it very much.