[FIXED] How to find the points of intersection between a circle and an ellipseM

Issue

So I have been trying to solve a question lately, given an ellipse and a circle such that the center of the circle lies on the ellipse. We need to find the area between the two curves. As for inputs we have a,b (The axes of the ellipse). R(The radius of the circle) and Theta(The angle that the circles center makes with the horizontal axis) Also the ellipse is centred at the origin.

I figured out that we can get the co ordinates for the center of the circle and then solve the two equations giving me the necessary points which can be then used to get the correct area. However finding those points actually leads to a fourth order equation, and I am stuck at finding the necessary intersection points. I tried solving the equation using numpy but that also gave a huge overhead.

Another approach that I was thinking of was trying to fit a convex polygon within the region and then return it's area however I am clueless with regards to the accuracy and complexity of this method. Thank you.

Solution

I think the easiest way to do this is simply by numerical integration.

You can find the lower and upper curves for both circle and ellipse and hence the top and bottom parts of the area at any x value. You can also readily find the limits of integration.

Note also that, for the purposes of finding the relevant area, you can rotate and reflect such that the circle centre is in the first quadrant of the ellipse.

I think the below code is OK, but you'll have to count squares on the graph as a rough check. Note that the angle is specified in degrees, between 0 and 360.

from math import cos, sin, sqrt, pi
import matplotlib.pyplot as plt
import numpy as np


def getOverlap( a, b, R, phi_deg, N ):
    phi_rad = phi_deg * pi / 180.0

    # Identify point of intersection (= centre of circle)
    scale = a / sqrt( cos( phi_rad ) ** 2 + ( a * sin( phi_rad ) / b ) ** 2 )
    x0, y0 = scale * cos( phi_rad ), scale * sin( phi_rad )
    # WLOG, work in first quadrant
    xc, yc = abs( x0 ), abs( y0 )

    # Set maximum possible limits of numerical integration
    # Note: height of intersection may still be 0 between part of these bounds
    left, right = max( xc - R, -a ), min( xc + R, a )

    # Numerical integration (mid-ordinate rule)
    dx = ( right - left ) / N
    area = 0.0
    for i in range( N ):
        x = left + ( i + 0.5 ) * dx
        dy = sqrt( R ** 2 - ( x - xc ) ** 2 )
        circle_top    = yc + dy
        circle_bottom = yc - dy
        dy = b * sqrt( 1.0 - ( x / a ) ** 2 )
        ellipse_top    =  dy
        ellipse_bottom = -dy
        top    = min( circle_top   , ellipse_top    )
        bottom = max( circle_bottom, ellipse_bottom )
        if top > bottom: area += ( top - bottom ) * dx

    return area, x0, y0



a, b, R, phi_deg = 3.0, 1.0, 2.5, 30.0
N = 1000000
overlap, x0, y0 = getOverlap( a, b, R, phi_deg, N )
print( "Area of overlap = ", overlap )


npts = 1000
t = np.linspace( 0.0, 2 * pi, npts )
xellipse, yellipse = a * np.cos( t ), b * np.sin( t )
xcircle , ycircle  = x0 + R * np.cos( t ), y0 + R * np.sin( t )
plt.plot( xcircle, ycircle )
plt.plot( xellipse, yellipse )
plt.axis( "square" )
plt.grid()
plt.show()

Output:

Area of overlap =  6.213925590052305

Graph:

An alternative case with 4 points of intersection instead of 2 would be a, b, R, phi_deg = 5.0, 1.0, 3.0, 90.0 whence

Area of overlap =  10.49361094712482

Answered By - lastchance