When two spheres overlap in 3D, they create a circle, where they intersect. I would like to be able to create a (random) point on this circle.

Here is what I have so far:
Axis a through the center of the circle (long blue line)
Center of circle (where both blue lines meet)
Radius of the circle (based on 2D projection of the problem)

As a next step I would need to calculate a vector p, which, either from the circle center or from the origin points to the point on the circle.

I know that, if the circle is in the xy-plane, I could use p = a/2 +ricos(t) + rjsin(t) with i and j the unit vectors for x and y, however because the circle is not in the xy-plane I would need to create two unit vectors in the circle plane perpendicular to the axis a but I don’t know how? Any suggestions?

if a is your vector between the sphere centers, then if you cross product it with any other vector b not parallel to a, you will get a new vector, call it i, that is perpendicular to a. Cross a with i to get j. Then normalize i and j.

To come up with b, pick the shortest component of a, that is, the component whose abs() is closest to 0, set that component to 1 and the other two to 0. So, if a was ( 7, -1, -8 ), you would choose b as ( 0, 1, 0 ).

In the circle equation I wrote on top I use a/2 to point to the middle point.
If I do that, then all coordinates calculated have the origin of a as origin. That means to show them correctly I have to add the red vector (where the blue axis vector starts).

If I use the vector pointing from the coordinate system origin to the circle middle point instead of a/2, then the coordinates need no correction.