directional cosine matrix for coordinate transformation

Good Afternoon All,

I know this is highly engineering based but the photo shown below is a reference paper to calculate the Torque Roll Axis. It sounds all crazy fancy but to my understanding its a coordinate transformation to make the new x axis (x') parallel to the x moment of inertia. I have applied his formulas to create this new DCM for transformation. I have called his directional cosine vector for x' Tx, the directional cosine vector for y' Ty and the directional cosine vector for z' Tz. I am not sure if I am producing his equations correctly in matlab. i.e maybe using the wrong functions. For example when applying an arbitrary value it shouldn't matter what value I choose correct? However when changing the arbitrary values it does change my DCM. I have also checked to see if the TRA_Transform_Matrix is orthogonal and it does meet that requirement.

If anyone could take a look at the code and tell me if I am on the right track I would greatly appreciate it. I know its a lot to ask but I am getting a crazy DCM.

Thanks so much.

Here is just some sample test data:

CG=[1.429925 .01714458 .7631665]; 
Mount=[2.32062 -.10742 .51072];
I=[17.24533 66.20312 62.17456 -.5598926 13.70156 1.363275]
%I is the only variable that serves an actual importance.

function[TRA_Transform]= TRA_Transformation_Matrix(CG,Mount,I)

%This program calculates the Transformation Matrix for the TRA
%
% User Input
%This is a subfunction of the main program of PT_Wizard_NF_Calculator.
%Please see the main routine document for further information of inputs.
%  
%   CG: Center of Gravity Position Matrix expressed as [xc,yc,zc]
%   Mount: Any Mount Position Matrix expressed as [x,y,z]
%   I: Moment of Inertia Matrix expressed as [Ixx,Iyy,Izz,Ixy,Ixz,Iyz]
%

%Arbitrary Values in TRA_Transform (distances)
ry=CG(2)-Mount(2);
ryz=CG(2)-Mount(3);
rz=CG(3)-Mount(3);

%Mass Inertia Matrix 
M_Int=zeros(3,3);
M_Int(1,1)=I(1);
M_Int(1,2)=-I(4);
M_Int(1,3)=-I(5);
M_Int(2,1)=-I(4);
M_Int(2,2)=I(2);
M_Int(2,3)=-I(6);
M_Int(3,1)=-I(5);
M_Int(3,2)=-I(6);
M_Int(3,3)=I(3);

%Inverse of Mass-Inertia Matrix to find TRA Direction
Inv_M_Int=inv(M_Int);
%Formulating the norm vector used for finding the normalization constant
Norm_Vec=[Inv_M_Int(1,1) Inv_M_Int(2,1) Inv_M_Int(3,1)];
%Calculating the normalization constant (a)
Norm_Constant_a=norm(Norm_Vec);
%Applying normalization to M_Rot_Inv
Norm_Inv_M_Int=Inv_M_Int/Norm_Constant_a;
%TRA direction (assuming torque applied to x-axis)
TRA_Direction=[0 0 0 Norm_Inv_M_Int(1,1) Norm_Inv_M_Int(2,1) Norm_Inv_M_Int(3,1)];

%TRA DCM for x
TRA_Xx=TRA_Direction(4);
TRA_Xy=TRA_Direction(5);
TRA_Xz=TRA_Direction(6);

%TRA DCM for y
TRA_Yxi=(-(TRA_Xy*ry+TRA_Xz*ryz)/TRA_Xx);
TRA_Yyi=ry;
TRA_Yzi=ryz;
%Creating TRA Transform Y vector for norm
TRA_Y_Vec=[TRA_Yxi TRA_Yyi TRA_Yzi];
%TRA Transform for y
TRA_Yx=TRA_Yxi/norm(TRA_Y_Vec);
TRA_Yy=TRA_Yyi/norm(TRA_Y_Vec);
TRA_Yz=TRA_Yzi/norm(TRA_Y_Vec);

%TRA DCM for z
%Initial Matrix of TRA Transform z
Int_TRA_z=[TRA_Xx TRA_Xy; TRA_Yx TRA_Yy];
%Solving for the z elements for TRA Transform
TRA_z=-inv(Int_TRA_z)*[TRA_Xz; TRA_Yz]*(rz);
TRA_Zxi=TRA_z(1);
TRA_Zyi=TRA_z(2);
TRA_Zzi=rz;
%Creating TRA Transform Y vector for norm
TRA_Z_Vec=[TRA_Zxi TRA_Zyi TRA_Zzi];
%TRA Transform for z
TRA_Zx=TRA_Zxi/norm(TRA_Z_Vec);
TRA_Zy=TRA_Zyi/norm(TRA_Z_Vec);
TRA_Zz=TRA_Zzi/norm(TRA_Z_Vec);

%TRA DCM Formulation
TRA=zeros(3,3);
TRA(1,1)=TRA_Xx;
TRA(1,2)=TRA_Xy;
TRA(1,3)=TRA_Xz;
TRA(2,1)=TRA_Yx;
TRA(2,2)=TRA_Yy;
TRA(2,3)=TRA_Yz;
TRA(3,1)=TRA_Zx;
TRA(3,2)=TRA_Zy;
TRA(3,3)=TRA_Zz;

%Completed TRA Transformation Matrix
zeros_temp=zeros(3,3);
TRA_Transform=[TRA zeros_temp; zeros_temp TRA];