Visualizing lemniscates and Cassinian curves

Last night was fun. I found a USB stick I haven't used in years and couldn’t resist browsing its content. My excitement rose when I retrieved - among various other juicy and less juicy pieces of data, mostly dating back from the college era – an archive containing many old C programs. I stumbled upon this piece of code /rambling_begins I wrote at a time my friends and I were venerating Brian Greene et al. and our bible was the Elegant Universe. Eventually, my friends went to study physics but did not major in cosmology. Applied physics pay the bills. String theory physicists are starving to death nowadays! /rambling_over.

So the idea was to visually represent lemniscates.

Lemniscates and its derivatives are fascinating mathematical objects. (Their construction was appealing because it resembled the ellipse's.)

Remember that an ellipse is defined by its two centers C0, C1 and the set of points P such that dist(C0,P)+dist(C1,P) is constant:

Image (c) Wikipedia

A lemniscate is also defined by its two centers C0, C1 and a set of points P, but this time d=dist(C0,P)*dist(C1,P) must be constant. The most famous lemniscate is Bernouilli’s, where d is set to sqrt(dist(C0,C1)/2). This produces … the infinite sign:

Image (c) Mathcurve.com

Generalized, these lemniscates are formally called Cassini’s ovals, and we can intuitively guess that greater values of d generate what look like ellipses, whereas smaller values of d generate pairs of pseudo-circles:

Image (c) Mathcurve.com

Now the last step was to generalize that to 3, 4 and more centers. The formal name of these beauties is Cassinian curves, and they look like that:

Image (c) Mathcurve.com

Back in the days, I didn’t know these objects had real names – not to mention their use in electrostatic fields analysis. I was experimenting and simply called them generalized lemnicates. Anyway, below are a few of screenshots with 2 centers, 3 centers, 6 centers, and a couple of zoomed-in areas.

The trick to get these cool gradients is to rely on the RGB implementation used by most programming frameworks. The (Red, Green, Blue) tuple is represented by 3 bytes stored in a 32-bit dword. (The high-byte can be the alpha channel, but is irrelevant in the case of raw colors.) If we define d=dist(C0,P)*dist(C1,P)*…*dist(Cn,P) and represent d as a colored pixel whose RGB value is directly proportional to d (ie, a number between 0 and 0FFFFFFh).

The yellow-ish separations mean that the Blue byte is about to be incremented. Of course, this depends on the number of centers and how close they are to each other: these initial conditions yield the range of value d can have within the window. With more 10 centers, very close to each other, it takes a while to get to Blue=1, therefore we can see increments of Green (at least the first 4 are clearly marked out).

I’ve attached the source and binary so you can experiment yourself. (Sorry for the Linux folks out there, it’s Windows C…)

ProjectEuler problem 144

Project Euler is a website offering math challenges solvable by programming - usually few tens lines of code are enough. Last friday, after more than a year of Euler abstinence, I decided to have a look at it again, and ran across this cool problem (#144, check it out here).



So the goal is to find the number of times the laser beam will hit the ellipse-shaped cell before exiting through the input cavity.
First, let's formulate the problem with simple trigonometry. Have a look at the schema below:



The laser beam comes from A and hits the ellipse at point 0. We can calculate the slope m0: (yO-yA)/(xO-xA); we can also calculate the slope m1: -4*x/y. We need to calculate slope m2 and deduce B in order to iterate this process (find the next point of impact).

It appears clear that:
a = a0-a1
and
a = a1+PI-a2

Therefore:
tan(a) = tan(a0-a1)
= (tan(a0)-tan(a1))/(1+tan(a0)*tan(a1))
= (m0-m1)/(1+m0*m1)

and
tan(a) = tan(a1+PI-a2)
= tan(a1-a2)
= (tan(a1)-tan(a2))/(1+tan(a1)*tan(a2))
= (m1-m2)/(1+m1*m2)

Let's set X = (m0-m1)/(1+m0*m1)
X = (m1-m2)/(1+m1*m2)
-> m2 = (m1-X)/(1+X*m1)

We have the slope m2 of the reflected beam (OB). With the coordinates of O (x0,y0), we have the equation of the line:
y = a*x+b, with a=m2 and b=y0-m2*x0

The last step is calculate the coordinates of B. We have:
- the equation of the ellipse: 4*x^2+y^2 = 100
- the equation of (0B): y = m2*x+b

It boils down to solving the quadratic:
x^2*(4+a^2) + x*(2*a*b) + (b^2-100) = 0

And then we yield the coordinates of B. After a few 300-and-something iterations (nah I wouldn't spoil your Euler fun) you'll find out the reflected point falls onto the exit hole. Problem solved!