Description 
[center,rad,v1,v2] = circlefit3d(p1,p2,p3)
This function computes centers and radii of circles, that interpolate some given triples of points in 3d space.
The function accepts the input of matrices of trial points also, returning a matrix of corresponding centroids and a vector of corresponding radii. Corresponding points are assumed to be located in the same row in each matrix p1, p2 and p3.
Several error checks are done (collinearity, size of inputs) and return different error codes in return argument rad.
The function uses vector algebra only, is fully vectorized and should run quite fast.
Additional outputs are basis vectors v1 and v2 to describe the plane of the resulting circle.
Check:
p1=rand(4,3);p2=rand(4,3);p3=rand(4,3);
[center,rad,v1,v2] = circlefit3d(p1,p2,p3);
plot3(p1(:,1),p1(:,2),p1(:,3),'bo');hold on;plot3(p2(:,1),p2(:,2),p2(:,3),'bo');plot3(p3(:,1),p3(:,2),p3(:,3),'bo');
for i=1:361,
a = i/180*pi;
x = center(:,1)+sin(a)*rad.*v1(:,1)+cos(a)*rad.*v2(:,1);
y = center(:,2)+sin(a)*rad.*v1(:,2)+cos(a)*rad.*v2(:,2);
z = center(:,3)+sin(a)*rad.*v1(:,3)+cos(a)*rad.*v2(:,3);
plot3(x,y,z,'r.');
end
axis equal;grid on;rotate3d on;
