单晶衍射仪系列2-Kappa测角仪旋转矩阵

casjxm

kappa测角仪有三个角可以旋转（o: omega, k: kappa, p: phi）,暂不考虑(t: 2sita)。
如果我们去查3D旋转矩阵时，通常采用的是三个欧拉角（分别绕坐标轴x,y,z旋转的3个角）的旋转矩阵来对3D空间中的矢量或坐标系进行旋转，陀螺仪也采用的是这种：
https://blog.csdn.net/qq_44092699/article/details/136036998
这里不做细节描述，欧拉角旋转方式有个不足的地方是旋转后的结果与3个欧拉角的旋转顺序是有关的，也就是，如果3个欧拉角不变，旋转顺序不同，XYZ, YZX, ...，等等，最终的旋转结果是不同的，这在单晶衍射测试时是非常棘手的，因为晶体在测试过程中，旋转的顺利是很难固定的。
kappa测角仪有个神奇之处是（o,k,p）3个角旋转结果与顺序无关，也就是，无论旋转的过程是怎样的，只要o,k,p固定，最终结果都不会变。
测角仪：


（o,k,p）设置：


在推导kappa测角仪旋转矩阵时，需要用到罗德里格斯公式：

https://blog.csdn.net/u014265928/article/details/117170481

有个特别需要注意的地方，旋转矩阵必须要明确是用来左乘矢量，还是右乘；是顺时针旋转，还是逆时针旋转；是旋转坐标系，还是旋转矢量；

规律：（顺时针，逆时针）（左乘，右乘）（旋转坐标系，旋转矢量），3种设置同时置换两种，旋转矩阵不变；否则转置；

直接上代码（matlab）：

均采用：顺时针旋转，旋转矢量，左乘（v'=R*V, 这里的v是列向量）

绕w旋转矩阵：

function result=rot_omega_mat(w)

%[nx ny nz]=[0 0 1];

n=[0 0 1];

nx=n(1,1);ny=n(1,2);nz=n(1,3);

k=w;

res=[(1-cos(k))*nx^2+cos(k) (1-cos(k))*nx*ny-nz*sin(k) (1-cos(k))*nx*nz+ny*sin(k) ;

(1-cos(k))*nx*ny+nz*sin(k) (1-cos(k))*ny^2+cos(k) (1-cos(k))*ny*nz-nx*sin(k) ;

(1-cos(k))*nx*nz-ny*sin(k) (1-cos(k))*ny*nz+nx*sin(k) (1-cos(k))*nz^2+cos(k) ];

%anticlockwise, vector rotation, left multiply.

result=res';

%change res to closewise

绕kappa旋转矩阵：

%rotation matrix around the kappa axis
%clockwise, vector rotation, left multiply.
%v'=R(n,t)v
% ok: angle between kappa axis and z
% k: kappa of rotation
% t: sita of rotation
% w2=pi/2-w,w:omega of rotation
% x,y,z, projection of kappa axis on the x y and z ordinate system
function result=rot_kap_mat(w,k,ok)
w2=pi-w;
r=1.0;
nx=r*sin(ok)*cos(w2);
ny=r*sin(ok)*sin(w2);
nz=r*cos(ok);
% [nx ny nz] is the kappa axis vector after closewise rotating w around z axis
res=[(1-cos(k))*nx^2+cos(k) (1-cos(k))*nx*ny-nz*sin(k) (1-cos(k))*nx*nz+ny*sin(k) ;
(1-cos(k))*nx*ny+nz*sin(k) (1-cos(k))*ny^2+cos(k) (1-cos(k))*ny*nz-nx*sin(k) ;
(1-cos(k))*nx*nz-ny*sin(k) (1-cos(k))*ny*nz+nx*sin(k) (1-cos(k))*nz^2+cos(k) ];
% res:anticlockwise, vector rotation, left multiply.
result=res';
%change res to closewise

绕phi旋转矩阵：
function result=rot_phi_mat(w,k,p,ok)
w2=pi-w;
r=1.0;
nx=r*sin(ok)*cos(w2);
ny=r*sin(ok)*sin(w2);
nz=r*cos(ok);
% [nx ny nz] is the kappa axis vector after closewise rotating w around z axis
rotkapmat=[(1-cos(k))*nx^2+cos(k) (1-cos(k))*nx*ny-nz*sin(k) (1-cos(k))*nx*nz+ny*sin(k) ;
(1-cos(k))*nx*ny+nz*sin(k) (1-cos(k))*ny^2+cos(k) (1-cos(k))*ny*nz-nx*sin(k) ;
(1-cos(k))*nx*nz-ny*sin(k) (1-cos(k))*ny*nz+nx*sin(k) (1-cos(k))*nz^2+cos(k) ];
%[nx ny nz]'=rotkapmat'*[0 0 1]';
n=rotkapmat'*[0 0 1]';
nx=n(1,1);ny=n(2,1);nz=n(3,1);
% [nx ny nz] is the z axis vector after closewise rotating k around kappa axis
k=p;
res=[(1-cos(k))*nx^2+cos(k) (1-cos(k))*nx*ny-nz*sin(k) (1-cos(k))*nx*nz+ny*sin(k) ;
(1-cos(k))*nx*ny+nz*sin(k) (1-cos(k))*ny^2+cos(k) (1-cos(k))*ny*nz-nx*sin(k) ;
(1-cos(k))*nx*nz-ny*sin(k) (1-cos(k))*ny*nz+nx*sin(k) (1-cos(k))*nz^2+cos(k) ];
% res:anticlockwise, vector rotation, left multiply.
result=res';
%change res to closewise
%closewise rotating p around z axis after rotating k around kappa axis after rotating w around z axis

绕w,k,p旋转矩阵：
这里给出了6种旋转方式，会发现旋转矩阵不会因为旋转过程的不同而不同。
%rotation matrix of vector around omega,kappa and phi axes
%clockwise, vector rotation, left multiply.
function result=rot_all_mat(w,k,p,ok)
res1=rot_phi_mat(w,k,p,ok)*rot_kap_mat(w,k,ok)*rot_omega_mat(w);
res2=rot_phi_mat(w,k,p,ok)*rot_omega_mat(w)*rot_kap_mat(0,k,ok);
res3=rot_kap_mat(w,k,ok)*rot_phi_mat(w,0,p,ok)*rot_omega_mat(w);
res4=rot_kap_mat(w,k,ok)*rot_omega_mat(w)*rot_phi_mat(0,0,p,ok);
res5=rot_omega_mat(w)*rot_phi_mat(0,k,p,ok)*rot_kap_mat(0,k,ok);
res6=rot_omega_mat(w)*rot_kap_mat(0,k,ok)*rot_phi_mat(0,0,p,ok);
diff=res2+res3+res4+res5+res6-5*res1;
diff2=sum(sum(diff))
result=res1;

测试代码，随机测试了一些不同的w,k,p：
rot_all_mat(10/180*pi,30/180*pi,50/180*pi,100/180*pi);
rot_all_mat(-10/180*pi,30/180*pi,50/180*pi,100/180*pi);
rot_all_mat(10/180*pi,-30/180*pi,50/180*pi,100/180*pi);
rot_all_mat(10/180*pi,30/180*pi,-50/180*pi,100/180*pi);
rot_all_mat(10/180*pi,30/180*pi,50/180*pi,-100/180*pi);
rot_all_mat(100/180*pi,100/180*pi,100/180*pi,30/180*pi);
rot_all_mat(-100/180*pi,-100/180*pi,100/180*pi,30/180*pi);
rot_all_mat(100/180*pi,-100/180*pi,-100/180*pi,30/180*pi);
rot_all_mat(100/180*pi,100/180*pi,-100/180*pi,-30/180*pi);
rot_all_mat(-100/180*pi,100/180*pi,100/180*pi,-30/180*pi);

运行结果：
>> test
diff2 =
  -1.8319e-15
diff2 =
   2.3315e-15
diff2 =
  -1.6653e-15
diff2 =
   1.6653e-15
diff2 =
  -8.3267e-16
diff2 =
   4.4409e-16
diff2 =
   1.7764e-15
diff2 =
  -4.5519e-15
diff2 =
  -8.8818e-16
diff2 =
  -8.8818e-16

测试通过！