Ellipse Intersect Ellipse

[STxAxTIC]

This one was a lot to chew. I think I worked it out properly on paper for the general case of any ellipse touching any other ellipse, but I had to do a mini-problem along the way to sanity-check my equations. It was pretty fun, so I thought I'd not only start this post to fill in the answer eventually, but also throw a problem.

Consider two copies of the same ellipse, except the second is rotated 90 degrees and centered on the original as shown in the screenshot. If the long length is "a", and the short length is "b", Calculate the intersection points represented by "?". To be definite: if a=10 and b=5, Calculate the exact number for "?". 100 digits of precision will do.

[_vince]

Code: [Select]

deflng a-z
const sw = 640
const sh = 480
dim shared pi as double
pi = 4*atn(1)
declare sub circlef(x, y, r)
declare sub ellipse(x, y, rx, ry)

'screenres sw, sh, 32
screen 12


a = 280
b = 80

ellipse sw/2, sh/2, a, b
ellipse sw/2, sh/2, b, a

color 12'rgb(255,0,0)

dim r as double
dim p as double
p = a*tan(pi/4)/b
r = (a/(cos(pi/4))/sqr(1 + p*p))

circlef sw/2 + r*cos(pi/4), sh/2 - r*sin(pi/4), 10
circlef sw/2 - r*cos(pi/4), sh/2 - r*sin(pi/4), 10
circlef sw/2 - r*cos(pi/4), sh/2 + r*sin(pi/4), 10
circlef sw/2 + r*cos(pi/4), sh/2 + r*sin(pi/4), 10

sleep
system

sub circlef(x, y, r)
        x0 = r
        y0 = 0
        e = -r

        do while y0 < x0
                if e <=0 then
                        y0 = y0 + 1
                        line (x - x0, y + y0)-(x + x0, y + y0),, bf
                        line (x - x0, y - y0)-(x + x0, y - y0),, bf
                        e = e + 2*y0
                else
                        line (x - y0, y - x0)-(x + y0, y - x0),, bf
                        line (x - y0, y + x0)-(x + y0, y + x0),, bf
                        x0 = x0 - 1
                        e = e - 2*x0
                end if
        loop
        line (x - r, y)-(x + r, y),, bf
end sub

sub ellipse(x, y, rx, ry)
        a = 2*rx*rx
        b = 2*ry*ry

        x0 = rx
        y0 = 0
        xx = ry*ry*(1 - 2*rx)
        yy = rx*rx

        e = 0
        sx = b*rx
        sy = 0

        do while sx >= sy
                pset (x - x0, y + y0)
                pset (x + x0, y + y0)
                pset (x - x0, y - y0)
                pset (x + x0, y - y0)

                y0 = y0 + 1
                sy = sy + a
                e = e + yy
                yy = yy + a

                if 2*e + xx > 0 then
                        x0 = x0 - 1
                        sx = sx - b
                        e = e + xx
                        xx = xx + b
                end if
        loop

        x0 = 0
        y0 = ry
        xx = rx*ry
        yy = rx*rx*(1 - 2*ry)

        e = 0
        sx = 0
        sy = a*ry

        do while sx <= sy
                pset (x - x0, y - y0)
                pset (x + x0, y - y0)
                pset (x - x0, y + y0)
                pset (x + x0, y + y0)

                x0 = x0 + 1
                sx = sx + b
                e = e + xx
                xx = xx + b

                if 2*e + yy > 0 then
                        y0 = y0 - 1
                        sy = sy - a
                        e = e + yy
                        yy = yy + a
                end if
        loop
end sub


[STxAxTIC]

Perfect _vince.

Here's my writeup for the solution to any ellipse touching any other ellipse in the most general 2D case - that math is explosive but when we reduce it to the example on hand, things are pretty tight:

[_vince]

a cleaner expression for r

[STxAxTIC]

Isn't it even cleaner when you erase tan(pi/4)?

On that note, I thought mine was cleaner anyway! I still like it though.

[_vince]

I was mocking your dislike for trig!

[STxAxTIC]

LOL! You so got me.