Python: пересечение сфер

Проблему лучше решить с помощью тригонометрии.

Сводя задачу к двумерным кругам, мы могли бы:

import math
import numpy

class Circle():
    def __init__(self, cx, cy, r):
        """initialise Circle and set main properties"""
        self.centre = numpy.array([cx, cy])
        self.radius = r

    def find_intercept(self, c2):
        """find the intercepts between the current Circle and a second c2"""
        #Find the distance between the circles
        s = c2.centre - self.centre
        self.dx, self.dy = s
        self.d = math.sqrt(numpy.sum(s**2))
        #Test if there is an overlap.  Note: this won't detect if one circle completly surrounds the other.
        if self.d > (self.radius + c2.radius):
            print("no interaction")
        else:
            #trigonometry
            self.theta = math.atan2(self.dy,self.dx)
            #cosine rule
            self.cosA = (c2.radius**2 - self.radius**2 + self.d**2)/(2*c2.radius*self.d)
            self.A = math.acos(self.cosA)

            self.Ia = c2.centre - [math.cos(self.A+self.theta)*c2.radius, math.sin(self.A+self.theta)*c2.radius]
            self.Ib = c2.centre - [math.cos(self.A-self.theta)*c2.radius,-math.sin(self.A-self.theta)*c2.radius]

            print("Interaction points are : ", self.Ia, " and: ", self.Ib)


#define two arbitrary circles 
c1 = Circle(2,5,5)
c2 = Circle(1,6,4)


#find the intercepts
c1.find_intercept(c2)

#test results by reversing the operation
c2.find_intercept(c1)

Используя scipy.optimize.fsolve, вы можете найти корень заданной функции, исходя из первоначального предположения, которое находится где-то в диапазоне вашего решения. Я использовал этот подход для решения вашей проблемы, и, похоже, он работает для меня. Единственным недостатком является то, что он предоставляет вам только один перекресток. Чтобы найти второй, вам придется повозиться с начальными условиями, пока fsolve не найдет второй корень.

Сначала мы определяем наши сферы, определяя (произвольные) радиусы и центры для каждой сферы:

a1 = np.array([0,0,0])
r1 = .4
a2 = np.array([.3,0,0])
r2 = .5
a3 = np.array([0,.3,0])
r3 = .5

Затем мы определяем, как преобразовать обратно в декартовы координаты, заданные углы u,v

def position(a,r,u,v):
    return a + r*np.array([np.cos(u)*np.sin(v),np.sin(u)*np.sin(v),np.cos(v)])

Теперь мы думаем, из какого уравнения нам нужно найти корень. Для любой точки пересечения считается, что для идеального u1,v1,u2,v2,u3,v3 положения position(a1,r1,u1,v1) = position(a2,r2,u2,v2) = position(a3,r3,u3,v3) равны. Таким образом, мы находим три уравнения, которые должны быть нулями, а именно разности двух векторов положения. Фактически, поскольку каждый вектор имеет 3 компонента, у нас есть 9 уравнений, которых более чем достаточно, чтобы определить наши 6 переменных.

Мы находим функцию, которую нужно минимизировать, как:

def f(args):
    u1,v1,u2,v2,u3,v3,_,_,_ = args
    pos1 = position(a1,r1,u1,v1)
    pos2 = position(a2,r2,u2,v2)
    pos3 = position(a3,r3,u3,v3)
    return np.array([pos1 - pos2, pos1 - pos3, pos2 - pos3]).flatten()

fsolve требует одинакового количества входных и выходных аргументов. Поскольку у нас есть 9 уравнений, но только 6 переменных, я просто использовал 3 фиктивные переменные, чтобы размеры совпадали. Сглаживание массива в последней строке необходимо, поскольку fsolve принимает только 1D-массивы.

Теперь пересечение можно найти с помощью fsolve и (довольно случайного) предположения:

guess = np.array([np.pi/4,np.pi/4,np.pi/4,np.pi/4,np.pi/4,np.pi/4,0,0,0])

x0 = fsolve(f,guess)
u1,v1,u2,v2,u3,v3,_,_,_ = x0

Вы можете проверить правильность результата, подключив полученные углы к функции position.