Translating and Rotating points to measure them

Hello all,

I have 3 data points with X,Y,Z values in a co-ordinate system that I can't change. They're from another computer model and they show the centre point of a circle (A), a point on the radius of the circle (B) and the point of the peak flow (C).

My plan was to use matlab to translate these 3 points and then rotate them so points A and B line up on a plane in XY and point C is above them in the z axis. On removal of the z value for C it would create a point on the new xy plane that AB are on. I could then measure the the distance AC (peak flow radius) and AB (pipe radius).

I've code to try and achieve this and I think there's a mistake in it because it give me a AC>AB value, which shouldn't be possible.

Please may somebody help to point out where the code has gone wrong.

%Origin
O = [0 0 0]';           Ox = [O' 1]';
%Vector co-ordinates of 1 (A - Origin centreline), ...
% 2 (B - Radius of artery), and 3 (C - Radius of Peak Velocity)
A = [13.3856 -60.0377 15.8443]';           Ax = [A' 1]';
B = [26.9486 -51.0653 20.9265]';           Bx = [B' 1]';
C = [16.2240 -92.5594 40.8687]';           Cx = [C' 1]';

%Find the new i unit vector in the old co-ords (the unit vector along the new x axis)
AB = B - A;
magAB = sqrt(sum(AB.^2));
new_i=AB./magAB;

%Calculate the angle to rotate through Z when transforming
thetaZ = atan(new_i(2)/new_i(1));

%Hence, define the translation matrix (to move the origin to A) and ...
% the rotation matrixes (to align the new x axis with AB and new_i)
T = [1 0 0 -A(1) ; 0 1 0  -A(2) ; 0 0 1 -A(3) ; 0 0 0 1];
Rz = [cos(thetaZ) sin(thetaZ) 0 0 ; -sin(thetaZ) cos(thetaZ) 0 0 ; 0 0 1 0 ; 0 0 0 1];

Transform = Rz * T;

%transform Cx to the new co-ordinates by translation and rotation in Z
A_dash = round(Transform * Ax,10);
B_dash = round(Transform * Bx,10);
C_dash = round(Transform * Cx,10);
new_i_t = round(Transform * [new_i' 1]',4); new_i_t = new_i_t(1:3);
new_O = round(Transform * Ox,4); new_O = new_O(1:3);

A_dash = A_dash(1:3); B_dash = B_dash(1:3); C_dash = C_dash(1:3);

%Perform a final rotation in Y
thetaY = atan(B_dash(3)/B_dash(1));
Ry = [cos(thetaY) 0 sin(thetaY) ; 0 1 0 ; -sin(thetaY) 0 cos(thetaY)];
B_dash = Ry * B_dash;
C_dash = Ry * C_dash;

%Taking point C onto the plane of AB
C_dashflat = C_dash.*[1 1 0]';

%Find Radius of Peak Flow 
Radius_Peak_Flow = sqrt(sum(C_dashflat.^2));

%Find Radius of Artery  
Radius_Artery = magAB;
%Ratio (between 0 -1)
Ratio = Radius_Peak_Flow/Radius_Artery