MATLAB: Intersection between ellipse and circle

circleellipseintersection

I'm trying to write a code to find the intersection between an ellipse and a circle. This is what I've got so far:
%ellipse
e=0.284576477148; %eccentricity
a=0.803468308684*150*10^6; c=e*a; b=sqrt(a^2-c^2); %parameters of ellipse
beta=linspace(0,2*pi,100);
X=a*cos(beta)-b*sin(beta)+c; %ellipse is shifted
Y=a*cos(beta)+b*sin(beta);
plot(X,Y,'r'), hold on
%circle
r=150*10^6; %radius
x=r*cos(beta);
y=r*sin(beta);
plot(x,y,'b')
axis square
z=[];
for i=1:length(beta)
  f1=a*cos(beta(i))-b*sin(beta(i))+c-r*cos(beta(i)); %X coordinates of the ellipse - x coordinates of the circle
  f2=a*cos(beta(i))+b*sin(beta(i))-r*sin(beta(i)); %Y coordinates - y coordinates
  if f1<=0.1*10^-20 && f2<=0.1*10^-20
      z=[z beta(i)];
  end
end
for k=1:length(z)
    p=r*cos(z(k));
    q=r*sin(z(k));
    plot(p,q,'*'), hold on
end
I'm trying to look for the angle beta where the curves intersect by looking for the value of beta that makes the expressions f1 and f2 close to zero. However, this code finds too many points that are actually quite far from the intersections, no matter how low I set the tolerance (if I set it to 0 no points are shown). What am I doing wrong? Thanks!

Best Answer

Try this: 
%ellipse
e=0.284576477148; %eccentricity
a=0.803468308684*150*10^6; c=e*a; b=sqrt(a^2-c^2); %parameters of ellipse
beta=linspace(0,2*pi,100);
X=a*cos(beta)-b*sin(beta)+c; %ellipse is shifted
Y=a*cos(beta)+b*sin(beta);
figure
plot(X,Y,'r')
hold on
%circle
r=150*10^6; %radius
x=r*cos(beta);
y=r*sin(beta);
plot(x,y,'b')
hold off
axis equal
[incirc,oncirc] = inpolygon(x,y, X,Y);                      % Logical Vector Of Points In Circle
[inelps,onelps] = inpolygon(X,Y, x,y);                      % Logical Vector Of Points In Ellipse
Ylix2 = find(Y == min(Y(inelps)));                          % Find Indices Of Minimum Y
Yuix2 = find(Y == max(Y(inelps)));                          % Find Indices Of Maximum Y
Max = [X(Yuix2), Y(Yuix2)];                                 % Approximate Intercept Values

Min = [X(Ylix2), Y(Ylix2)];                                 % Approximate Intercept Values
figure
plot(x, y, '-b')
hold on
plot(X, Y, '-r')
% plot(x(incirc), y(incirc), '.k')                            % In-Circle Points
% plot(X(inelps), Y(inelps), '.g')                            % In-Ellipse Points
plot(X(Yuix2), Y(Yuix2), 'pc', 'MarkerSize',10)             % Upper Intersection
plot(X(Ylix2), Y(Ylix2), 'pm', 'MarkerSize',10)             % Lower Intersection
hold off
axis equal
producing for figure(2) —
Intersection between ellipse and circle - 2019 11 11.png
You can use those indices to define a region to calculate more precise intersections. Probably the easiest way to get more precision is just to increase the resolution of ‘beta’.
Related Question