Some math
AFAIU your problem is that you have a sphere and a point on the sphere and you want to add 4 more points on the same sphere that would form a kind of a cross on the surface of the sphere around the target point.
I think it is easier to think about this problem in terms of vectors. You have a vector from the center of the sphere to your target point V of size R. All the point lying on the distance d from the target point form another sphere. The crossing of two sphere is a circle. Obviously this circle lies in a plane that is orthogonal to V. Solving a simple system of equations you can find that the distance from the target point to that plane is d^2/(2*R). So the vector from the center of the original sphere to the center of the circle:
Vc = V * (1 - d^2/(2*R^2))
and the radius of that circle is
Rc = sqrt(d^2 - (d^2/(2*R))**2)
So now to select 4 points, you need to select two orthogonal unit vectors lying in that plane D1 and D2. Then 4 points would be Vc + Rc*D1, Vc - Rc*D1, Vc + Rc*D2, and Vc - Rc*D2. To do this you may first select D1 fixing z =0 and switch x and y in Vc
D1 = (Vy/sqrt(Vx^2+Vy^2), -Vx/sqrt(Vx^2+Vy^2), 0)
and then find D2 as a result of cross-product of V and D1. This will work unless unless Vx = Vy = 0 (i.e. V goes along the z-axis) but in that case you can select
D1 = (1,0,0)
D2 = (0,1,0)
Some code
And here is some Python code that implements that math:
def cross_product(v1, v2):
return (v1[1] * v2[2] - v1[2] * v2[1],
v1[2] * v2[0] - v1[0] * v2[2],
v1[0] * v2[1] - v1[1] * v2[0])
def find_marks(sphereCenter, target, d):
lsc = list(sphereCenter)
lt0 = list(target)
lt1 = map(lambda c1, c0: (c1 - c0), lt0, lsc) # shift everything as if sphereCenter is (0,0,0)
rs2 = sum(map(lambda x: x ** 2, lt1)) # spehere radius**2
rs = rs2 ** 0.5
dv = d ** 2 / 2.0 / rs
dvf = d ** 2 / 2.0 / rs2
lcc = map(lambda c: c * (1 - dvf), lt1) # center of the circle in the orthogonal plane
rc = (d ** 2 - dv ** 2) ** 0.5 # orthogonal circle radius
relEps = 0.0001
absEps = relEps * rs
dir1 = (lt1[1], -lt1[0], 0) # select any direction orthogonal to the original vector
dl1 = (lt1[0] ** 2 + lt1[1] ** 2) ** 0.5
# if original vector is (0,0, z) then we've got dir1 = (0,0,0) but we can use (1,0,0) as our vector
if abs(dl1) < absEps:
dir1 = (rc, 0, 0)
dir2 = (0, rc, 0)
else:
dir1 = map(lambda c: rc * c / dl1, dir1)
dir2 = cross_product(lt1, dir1)
dl2 = sum(map(lambda c: c ** 2, dir2)) ** 0.5
dir2 = map(lambda c: rc * c / dl2, dir2)
p1 = map(lambda c0, c1, c2: c0 + c1 + c2, lsc, lcc, dir1)
p2 = map(lambda c0, c1, c2: c0 + c1 + c2, lsc, lcc, dir2)
p3 = map(lambda c0, c1, c2: c0 + c1 - c2, lsc, lcc, dir1)
p4 = map(lambda c0, c1, c2: c0 + c1 - c2, lsc, lcc, dir2)
return [tuple(p1), tuple(p2), tuple(p3), tuple(p4)]
For an extreme case
find_marks((0, 0, 0), (12, 5, 0), 13.0 * 2 ** 0.5)
i.e. for a circle of radius 13 with a center at (0,0,0), the target point lying on the big circle in the plane parallel to the xy-plane and d = sqrt(2)*R, the answer is
[(4.999999999999996, -12.000000000000004, 0.0),
(-5.329070518200751e-15, -2.220446049250313e-15, -13.0),
(-5.000000000000006, 12.0, 0.0),
(-5.329070518200751e-15, -2.220446049250313e-15, 13.0)]
So two points (2-nd and 4-th) are just two z-extremes and the other two are 90° rotations of the target point in the xy-plane which looks quite OK.
For a less extreme example:
find_marks((1, 2, 3), (13, 7, 3), 1)
which is the previous example with d reduced to 1 and with the original center moved to (1,2,3)
[(13.34882784191617, 6.06281317940119, 3.0),
(12.964497041420119, 6.985207100591716, 2.000739918710263),
(12.580166240924067, 7.907601021782242, 3.0),
(12.964497041420119, 6.985207100591716, 3.999260081289737)]
which also looks plausible
