Q: How to find the intervals that are intersection of a circle with a rectangle I need to integrate over the arcs that are resulted from the intersection of a circle with a rectangle and fall inside the rectangle. I can find the intersection points using the shapely package. However, I don't know how to obtain integration intervals. For example, in the below figure my code returns [-2.1562, 2.1562] in radians (with respect to the center of the circle), while it should be able to automatically understand that the integration intervals that falls inside the rectangle are [[2.1562, 3.1415],[-3.1415, -2.1562]] (assuming pi = 3.1415). 

Here is another example:

My code returns [-0.45036, -0.29576, 0.29576, 0.45036] and the expected intervals will be [[0.29576, 0.45036], [-0.45036, -0.29576]].
The code should also work for any other location that the circle is located (with any radius), whether its center is outside or inside the rectangle.
Here is my code, written using 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')

I tried to sort the intersection points according to their angle in a way that each two adjacent items in the list of intersection points corresponds to an arc. The problem is that there will be a jump in the angles from -pi to pi when rotating along the unit circle. Also I don't know how to find that whether an arc is inside the rectangle or not given its 2 end points.

A: Dealing with angle ranges is not straightforward.
1) select a non-ambiguous representation range, such as [-π, π) radians.
2) write a function that finds the intersections of the circle with a (h/v) half-plane and returns an angle interval. It the interval straddles the ±π border, split it in two.
3) write a function that finds the intersection between two lists of intervals (this is a modified merging problem).
4) process the four edges and intersect the resulting intervals.
5) possibly merge intervals that straddle the ±π border.

