Generate random point on the ellipsoid in 3d.

Question

I want to generate the trajectory of point M in space, such that the sum of its distances to two other points, P and Q in space, remains constant, that is, |PM| + |MQ| = constant. Therefore, I believe the trajectory of this point should be an ellipsoid with PQ as its foci and rotating about the PQ axis. However, I wrote some code for this but it failed to generate the correct trajectory. I noticed that the code only succeeds when the midpoint of P and Q is at the origin (0,0,0). In other cases, the generated trajectory is incorrect. There seems to be a logical issue somewhere in the code. Can anyone help?

clear;
clc;
d  = 50;
f1 = [0.5,7,1.0];
f2 = [3.5,6.0,0.0];
origin = [(f1(1) + f2(1))/2 ; (f1(2) + f2(2))/2; (f1(3) + f2(3))/2];
cc = sqrt((f1(1)-f2(1))^2 + (f1(2)-f2(2))^2 + (f1(3)-f2(3))^2)/2;
a = d / 2;
b = sqrt(a^2-cc^2);
c = b;
N = 99;
u = rand(1,N);
v = rand(1,N);
theta = u * 2 * pi;
phi = acos(v * 2 - 1);
sinTheta = sin(theta);
cosTheta = cos(theta);
sinPhi = sin(phi);
cosPhi = cos(phi);
rx = a .* sinPhi .* cosTheta;
ry = b .* sinPhi .* sinTheta;
rz = c .* cosPhi;
Q = [rx;ry;rz];
P = [cc,0,0];
P1 = f2;
cP = cross(P, P1);
cP = cP / norm(cP);
dP = dot(P, P1);
theta = acos(dP/(norm(P)*norm(P1)));
crossProductMatrix = [0, -cP(3), cP(2);
                     cP(3), 0, -cP(1);
                     -cP(2), cP(1), 0];
R = eye(3) + sin(theta) * crossProductMatrix + (1-cos(theta))*crossProductMatrix*crossProductMatrix;
Q1 = R * Q;
rx = Q1(1,:) + origin(1);
ry = Q1(2,:) + origin(2);
rz = Q1(3,:) + origin(3);
distance_f1 = sqrt((rx - f1(1)).^2 + (ry - f1(2)).^2 + (rz - f1(3)).^2);
distance_f2 = sqrt((rx - f2(1)).^2 + (ry - f2(2)).^2 + (rz - f2(3)).^2);
distance_sum = distance_f1 + distance_f2;
disp(distance_sum);

Answer 1

P = [cc,0,0]; 
P1 = f2-origin'; 
cP = cross(P, P1);

I've identified the mistake I made, which is that when calculating the angle of rotation, I should compute it relative to the center point. Therefore, P1 should first subtract the origin.