[principle_axes,new_mass_weight_matrix]=align_inertia(mass_weight_matrix)

function [principle_axes,new_mass_weight_matrix]=align_inertia(mass_weight_matrix)

inertia_tensor=inertia(mass_weight_matrix);

[principle_axes_unsort,inertia_unsort]=eig(inertia_tensor);

[inertia_sort,I]=sort(diag(inertia_unsort));

principle_axes=principle_axes_unsort(:,4-I);

right_handed_axes_test=cross(principle_axes(:,1),principle_axes(:,2));

angle=acos(dot(right_handed_axes_test,principle_axes(:,3))/(norm(right_handed_axes_test)*norm(principle_axes(:,3))));

if angle>pi/2
    
    principle_axes(:,3)=-principle_axes(:,3);
        
end

new_mass_weight_matrix=mass_weight_matrix*principle_axes;