Мне нужно интегрировать по дугам, полученным в результате пересечения круга с прямоугольником, и попасть внутрь прямоугольника. Я могу найти точки пересечения, используя пакет shapely
. Однако я не знаю, как получить интервалы интеграции. Например, на рисунке ниже мой код возвращает [-2.1562, 2.1562]
в радианах (относительно центра круга), тогда как он должен автоматически понимать, что интервалы интегрирования, попадающие в прямоугольник, равны [[2.1562, 3.1415],[-3.1415, -2.1562]]
(при условии, что число pi= 3.1415).
Вот еще один пример:
Мой код возвращается[-0.45036, -0.29576, 0.29576, 0.45036]
и ожидаемые интервалы будут [[0.29576, 0.45036], [-0.45036, -0.29576]]
.
Код также должен работать для любого другого местоположения, в котором находится круг (с любым радиусом), независимо от того, находится ли его центр снаружи или внутри прямоугольника.
Вот мой код, написанный с использованием iPython:
import matplotlib.pyplot as plt
import math
import numpy as np
from shapely.geometry import LineString, MultiPoint
from shapely.geometry import Polygon
from shapely.geometry import Point
# Utilities
def cart2pol(xy, center):
x,y = xy
x_0,y_0 = center
rho = np.sqrt((x-x_0)**2 + (y-y_0)**2)
phi = np.arctan2(y-y_0, x-x_0)
return(rho, phi)
def pol2cart(rho, phi, center):
x_0,y_0 = center
x = rho * np.cos(phi)+x_0
y = rho * np.sin(phi)+y_0
return(x, y)
def distance(A,B):
return math.sqrt((A[0]-B[0])**2+(A[1]-B[1])**2)
#######################
rad = 6
center = (-1,5)
p = Point(center)
c = p.buffer(rad).boundary
A = (10,0)
B = (0,0)
C = (0,10)
D = (10,10)
coords = [Point(A), Point(B), Point(C), Point(D)]
poly = MultiPoint(coords).convex_hull
i=c.intersection(poly)
lines = [LineString([A, D]), LineString([D, C]),
LineString([C, B]), LineString([B, A])]
points = []
for l in lines:
i = c.intersection(l)
if not i.is_empty:
if i.geom_type == 'MultiPoint':
for j in range(len(i.geoms)):
points.append(i.geoms[j].coords[0])
else:
points.append(i.coords[0])
# Repeat the tangential points
for k, point in enumerate(points.copy()):
if abs(distance(center, point)**2 + distance(point, B)**2 - distance(B, center)**2) < 1e-4:
points.insert(k+1,point)
elif abs(distance(center, point)**2 + distance(point, D)**2 -distance(D, center)**2) < 1e-4:
points.insert(k+1,point)
# Sort points in polar coordinates
phis = [cart2pol(point,center)[1] for point in points]
phis.sort()
print(phis)
# Plot the shapes
x,y = c.xy
plt.plot(*c.xy)
for l in lines:
plt.plot(*l.xy, 'b')
plt.gca().set_aspect('equal', adjustable='box')
Я попытался отсортировать точки пересечения по их углу таким образом, чтобы каждые два соседних элемента в списке точек пересечения соответствовалидуга. Проблема заключается в том, что при повороте по круговой единице произойдет скачок в углах от -pi до pi. Также я не знаю, как найти это, находится ли дуга внутри прямоугольника или нет с учетом его 2 конечных точек.