三维旋转的表示

三维旋转矩阵

我们考虑主动旋转，并且通过$$\vec{r} = R\vec{r}$$定义相应的旋转矩阵$$R$$，其中$$\vec{r} = (x,y,z)^T$$，于是绕固定轴（世界轴）XYZ(extrinsic rotations)的旋转矩阵$$R$$为

$$ R(\varphi,X) = \begin{pmatrix} 1&0&0 \\ 0&\cos(\varphi) & -\sin(\varphi) \\ 0& \sin(\varphi) & \cos(\varphi) \end{pmatrix}, R(\varphi,Y) = \begin{pmatrix} \cos(\varphi)&0&\sin(\varphi) \\ 0&1&0 \\ -\sin(\varphi)&0&\cos(\varphi) \end{pmatrix}, R(\varphi,Z) = \begin{pmatrix} \cos(\varphi) & -\sin(\varphi) & 0 \\ \sin(\varphi) & \cos(\varphi) & 0 \\ 0 & 0 & 1 \end{pmatrix} $$

于是连续两次绕世界轴旋转的旋转$$(Z(\varphi)-X(\theta))$$的旋转矩阵为$$R(\theta,X)R(\varphi,Z)$$

在ROOT中，TRotation类用来描述主动转动的旋转矩阵(TRotation class specifies the rotation of the object in the original system (an active rotation))，与上文的定义一致，使用方法如下

TRotation r; 
//令r为绕轴TVector3(3,4,5)旋转Pi/3的旋转矩阵
r.Rotate(TMath::Pi()/3,TVector3(3,4,5));
//对向量v进行旋转r
TVector3 v;
TRotation r;
v.Transform(r);
v *= r;  //Attention v = r * v
//变为恒等变换
r.SetToIdentity();
//得到绕固定轴Z旋转phi,再绕固定轴X旋转theta的旋转矩阵
r.RotateZ(phi)
r.RotateX(theta)

通过欧拉角表示三维旋转

欧拉角$$Z(\varphi)-X(\theta)-Z(\psi)$$

欧拉角$$Z(\varphi)-X(\theta)-Z(\psi)$$动图

考虑刚体的旋转时，我们取初始的刚体轴作为固定的世界轴

考虑内旋(intrinsic rotations)定义下（即绕刚体轴旋转）的$$Z(\psi)-X'(\theta)-Z''(\varphi)$$欧拉角旋转，依次旋转的角度为$$\psi,\theta,\varphi$$,它与外旋(extrinsic rotations)定义下（即绕世界轴旋转）的欧拉角$$Z(\varphi)-X(\theta)-Z(\psi)$$的旋转效果相同，旋转矩阵均为 $$ R_{ZX'Z''}(\psi,\theta,\varphi) = R_{ZXZ}(\varphi,\theta,\psi) = R(\psi,Z)R(\theta,X)R(\varphi,Z)'' $$

在ROOT中构建欧拉角旋转的旋转矩阵

////////////////////////////////////////////////////////////////////////////////
/// Rotate using the x-convention.
 
TRotation & TRotation::RotateXEulerAngles(Double_t phi,
                                         Double_t theta,
                                         Double_t psi) {
   TRotation euler;
   euler.SetXEulerAngles(phi,theta,psi);
   return Transform(euler);
}

////////////////////////////////////////////////////////////////////////////////
/// Rotate using the x-convention (Landau and Lifshitz, Goldstein, &c) by
/// doing the explicit rotations.  This is slightly less efficient than
/// directly applying the rotation, but makes the code much clearer.  My
/// presumption is that this code is not going to be a speed bottle neck.
 
TRotation & TRotation::SetXEulerAngles(Double_t phi,
                                      Double_t theta,
                                      Double_t psi) {
   SetToIdentity();
   RotateZ(phi);
   RotateX(theta);
   RotateZ(psi);
 
   return *this;
}

///得到一个(Z(phi)-X(theta)-Z(psi))的外旋欧拉角旋转(extrinsic rotations)的旋转矩阵，即(Z(psi)-X'(theta)-Z''(phi))内旋欧拉角旋转(ntrinsic rotations)

TRotation r;
r.SetXEulerAngles(phi,theta,psi);

例子与ROOT代码验证

我们考虑沿`X`轴的单位矢量$$\hat{\vec{x}}$$，即TVector3(1,0,0)，经过欧拉角$$(Z(30^{\circ})-X'(45^{\circ})-Z''(60^{\circ}))$$旋转后的坐标,即上文中令$$\psi=\frac{\pi}{6},\theta=\frac{\pi}{4},\varphi=\frac{\pi}{3}$$。

几何计算可得坐标为$$(\frac{\sqrt{3}}{4}-\frac{\sqrt{3}}{4\sqrt{2}},\frac{1}{4}+\frac{3}{4\sqrt{2}},\frac{\sqrt{3}}{2\sqrt{2}})\approx(0.13,0.78,0.61)$$

使用ROOT代码验证（两种方法+一种错误使用） 注：在TRotation的文档里有这么一句 With a minor caveat, the Euler angles of the rotation may be set using SetXEulerAngles() or individually set with SetXPhi(), SetXTheta(), and SetXPsi(). These routines set the Euler angles using the X-convention which is defined by a rotation about the Z-axis, about the new X-axis, and about the new Z-axis. This is the convention used in Landau and Lifshitz, Goldstein and other common physics texts. 这很容易让人困惑，因为SetXEulerAngles(phi,theta,psi)并不是如文档所说的按新的轴一次转phi,theta,psi,而是按固定轴依次转phi,theta,psi；另外，通过依次SetXPhi(), SetXTheta(), and SetXPsi()并不能正确的实现欧拉角旋转。下面验证这两件事。

TRotation R;
R.SetXEulerAngles(TMath::Pi()/3,TMath::Pi()/4,TMath::Pi()/6);
r = TVector3(1,0,0);
(R*r).Print()//并不会改变r的值
//输出TVector3 A 3D physics vector (x,y,z)=(0.126826,0.780330,0.612372)

//错误的示范
TRotation R2;
R2.SetXEulerAngles(TMath::Pi()/6,TMath::Pi()/4,TMath::Pi()/3);
(R2*r).Print()
//输出TVector3 A 3D physics vector (x,y,z)=(0.126826,0.926777,0.353553)

//另一种方式
TVector3 Z(0,0,1);
TVector3 X(1,0,0);
r = TVector3(1,0,0);
r.Rotate(TMath::Pi()/6,Z);//会改变r的值
X.Rotate(TMath::Pi()/6,Z);
r.Rotate(TMath::Pi()/4,X);
Z.Rotate(TMath::Pi()/4,X);
r.Rotate(TMath::Pi()/3,Z);
r.Print()
//输出TVector3 A 3D physics vector (x,y,z)=(0.126826,0.780330,0.612372) 

//通过文档所说的`SetXPhi(), SetXTheta(), and SetXPsi()`来进行欧拉角旋转，会得到错误的结果
root [13] TRotation R;
root [14] R.SetToIdentity();R.SetXPhi(TMath::Pi()/6);R.SetXTheta(TMath::Pi()/4);R.SetXPsi(TMath::Pi()/3);(R*TVector3(1,0,0)).Print()
TVector3 A 3D physics vector (x,y,z)=(0.324469,0.928023,0.183013) (rho,theta,phi)=(1.000000,79.454709,70.728583)
root [17] R.SetToIdentity();R.SetXPhi(TMath::Pi()/3);R.SetXTheta(TMath::Pi()/4);R.SetXPsi(TMath::Pi()/6);(R*TVector3(1,0,0)).Print()
TVector3 A 3D physics vector (x,y,z)=(0.573223,0.739199,0.353553) (rho,theta,phi)=(1.000000,69.295189,52.207654)
//均不能得到满意的结果，检查SetXPhi()的源代码可以发现，它调用了GetXPhi()等函数，但GetXPhi()并不是一个定义良好的函数，比如：
root [0] TRotation R
(TRotation &) Name: TRotation Title: Rotations of TVector3 objects
root [1] R.SetXPhi(TMath::Pi()/3)
root [2] R.GetXPhi()
(double) 0.52359878
root [3] R.SetXTheta(0)
root [4] R.SetXPsi(0)
root [5] R.GetXPhi()
(double) 0.26179939
root [6] R.GetXTheta()
(double) 0.0000000
root [7] R.GetXPhi()
(double) 0.26179939
//这也许是SetXPhi(), SetXTheta(), and SetXPsi()不能正常使用的原因