Intersection between ellipse and circle

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!

[darova]

What about graphical solution? polyxpoly or intersections

[Ludovica Lanzieri]

I'm sorry, how do these work?

Answer 1

[Star Strider]

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) —

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’.