Finding every point of intersection between any number of circles and ellipses centered at the origin

I'm trying to find each possible coordinate point of intersection between any number of circles and any number of ellipses, all centered at the origin. This originated from a problem of finding where the intersection points of a conical light distribution and an angled cylindrical light distribution at some plane would be located. I'm starting with the intersection between one circle and one ellipse and then progressively decreasing the size of the ellipses with each decrease in the radius of the circles. Basically, I'm trying to find each point of the following graph (of course with any number of ellipses and any number of circles):

I expect that the first two points would be located where the minor axis of the ellipse meets the radius of the circle on the (in this case) y-axis. For each circle size, multiple ellipses could intersect, depending on the ratio of the major and minor axis of the ellipse. (Of note: here we are assuming the major axis is always greater than the minor axis).

I know the equations for each point as well as the conditions necessary for a non-imaginary point. In a 3-dimensional coordinate axis, with intersections located at the xy-plane at z=0, I should have that each coordinate obey the following equation:

Moreover, the conditions to return a non-imaginary point are:

and

.

I have attempted to write the code that would run through and find each point and plot them on the xy-coordinate axis. Basically, I will run though each iteration of smaller and smaller ellipse for a circle size, beginning with 1 (m = multiplicative factor of circle size and n = multiplicative factor of ellipse size) and then repeating the whole process for each progressively smaller circle.

rho = 16; %distance from the plane on intersection
rad = 2; %radius of outer circle
PhiMin = 0;
PhiMax = atand(rad/rho);
Phi = linspace(PhiMin, PhiMax, 5);

Theta = 50; %angle of cone that results in the circle at the plane

b = rho*tand(PhiMax); %minor axis of the ellipse

a = (rho*tand(PhiMax))/cos(Theta); %major axis of the ellipse

r = rad; % radius of the circle

dummy = 0:0.1:1;
x1 = ones(length(dummy),length(dummy)); 
y1 = x1;
for j = 1:1:11 %step size for ellipse 
    n = (j-1)/10;
    for i = 1:1:11 %step size for circle 
        m = (i-1)/10;
        if ((m^2 >=n^2) && (a^2/b^2 >= m.^2/n.^2)) %conditions to get a real coordinate point of intersection
            x1(i,j) = sqrt(m^2*r^2-(b^2*((n^2*a^2-m^2*r^2)/(a^2-b^2)))); %x-coordinate equation
            y1(i,j) = b*sqrt((n^2*a^2-m^2*r^2)/(a^2-b^2)); %y-coordinate equation 
        else
        x1(i,j) = NaN; %if the above is not satisfied, we don't want imaginary numbers so instead we have that this is outputted instead
        y1(i,j) = NaN;
        end
    end
end

x = [-x1, x1];
y = [-y1, y1];

hold on
figure(1)
plot(x1,y1, 'o') %stitchiing together all four plots (I assume there is a much simpler way to go about this...) 
plot(-x1,y1, '- o')
plot(-x1,-y1,'- o')
plot(x1,-y1, '- o')

For some reason it's not giving me all the points and I wonder if there is a simpler solution (I'm sure there is).

Any help would be greatly appreciated. Thank you in advance!

[Guillaume]

I've already fixed many typos in your code, but I can't fix the undefined m. Please make sure that we can run the code you post.

[Gaëtan Poirier]

The program should run now, accidentally left that out when I copied it over here. Sorry about that.

Answer 1

[Guillaume]

Your initialisation variables are a bit strange. You create PhiMax and Phi but never use them. Your b calculation just reduces to b=rad and I don't see the point in creating a variable rad to then rename it to r. Choose one name and stick to it.

There's certainly no need for the loops:

rho = 16;         %distance from the plane on intersection
r = 2;            %radius of outer circle
Theta = 50;       %angle of cone that results in the circle at the plane
b = r;            %minor axis of the ellipse
a = b/cos(Theta); %major axis of the ellipse

m = ((1:11)' - 1)/10;  %COLUMN vector of m values
n = ((1:11) - 1)/10;   %ROW vector of n values
x = sqrt(m.^2 * r^2 - b^2 * (n.^2 .* a^2 - m.^2 * r^2) / (a^2 + b^2));  %See note below
y = b * sqrt((n.^2 .* a^2 - m.^2 * r^2) / (a^2 + b^2));                 %See note below
notreal = x ~= real(x) | y ~= real(y);   %simpler than your test
x(notreal) = NaN;
y(notreal) = NaN;

plot([x, x, -x, -x], [y, -y, y, -y], 'o');

Note that the denominator in your formulae is a^2 + b^2 but you've used a^2 - b^2 in your code. I don't know if it's another typo. I've used a^2 + b^2 here.

Also, I'm not bothering with your m^2/n^2 <= a^2/b^2 and m^2 >= n^2 tests. I'm just seeing which results are not pure real and setting them to NaN.