I wanted to play around with generating an “elliptic pencil of circles” starting with the Hermatian matrices. An elliptic pencil of circles are all the circles (or some of them) which all intersect at two specific points. Here is a simple one, where the points of intersection are at (0,1) and (0,-1), or ±i in the complex plane.

You see there is the unit circle, plus the other circles branching off to each side:

The Hermatian matrix template for this particular pencil is elegant:

where n is the x coordinate of the center of each circle (or the complex number center, not having an imaginary component).

Here is Racket code used to generate the plot:

(require math/array)
(require math/matrix)
(require plot)
(define (matrix-circle-radius M)
(let ([A (array-ref M #[0 0])]
[d (matrix-determinant M)])
(sqrt (/ d (* -1 (* A A))))))
(define (matrix-circle-center M)
(let ([C (array-ref M #[1 0])]
[A (array-ref M #[0 0])])
(/ C (* -1 A))))
(define (circle-isoline x y r c)
(isoline
(lambda (x_ y_) (sqrt (+ (sqr (- x_ x)) (sqr (- y_ y))))) r
#:color c
#:width 2))
(define (elipt-plot-demo)
(plot
(map (lambda (C)
(circle-isoline (real-part (first C))
(imag-part (first C))
(second C)
(list 30 160 210)
))
(map (lambda (M)
(list (matrix-circle-center M)
(matrix-circle-radius M)))
(map (lambda (n)
(matrix+
(matrix-scale (array #[#[0 1] #[1 0]]) n)
(array #[#[1 0] #[0 -1]])))
(range 2.5 -3 -0.5))))
#:x-min -5
#:x-max 5
#:y-min -5
#:y-max 5
#:width 400
#:height 400))

Schwerdteger describes an interesting “invariant” relationship between any two circles:

where

, , and

with the circles represented as Hermitian matrices .

The point is, in the end you come up with this single number which represents whether the smaller circle is inside the first, overlapping the first, just touching the first, or outside the first. Here are the cases:

common invariant > 1 : smaller circle contained withing the greater circle

common invariant = 1 : touching from the inside

-1 > common invariant > 1 : overlapping at two points

common invariant = -1 : touching from the outside

common invariant < -1 : completely outside

Here is Racket code used for calculations and plotting:

I borrowed the book Geometry of Complex Numbers by Schwerdtfeger. The material has been fascinating so far, though admittedly it took me about 3 hours to comprehend the material in page 1 of chapter 1. The first subject is this intriguing idea that you can represent a circle as a matrix. More specifically, you you can represent a circle with (complex) center 𝛾 and radius 𝜌 as a matrix

Where B = – A𝛾̅, C = -A𝛾, and D = A(𝛾𝛾̅-𝜌²), so that A and D are real numbers, and B and C are complex numbers. For normal circles, you will have A=1, though you can scale the matrix to have other (non-zero) values of A and still have the same circle. (If A is zero, you end up with a straight line, or some other mysterious thing that is not a circle.)

The simple case is the unit circle:

Which is

More exotic is the circle away from the origin, centered at 10+10i, i.e., (10, 10), with radius 3:

Here is some code to generate the matrix in Racket, and to put it back:

lang racket
(require math/array)
(require math/matrix)
(define (circle-to-matrix zC r)
(let ([B (* -1 (conjugate zC))]
[mzC (magnitude zC)])
(matrix [[ 1 B ]
[ (conjugate B) (- (* mzC mzC) (* r r)) ]])))
(define (matrix-circle-radius M)
(let ([A (array-ref M #[0 0])]
[d (matrix-determinant M)])
(sqrt (/ d (* -1 (* A A))))))
(define (matrix-circle-center M)
(let ([C (array-ref M #[1 0])]
[A (array-ref M #[0 0])])
(/ C (* -1 A))))

One can generate a few slight variations on the appearance of the previous fractal by altering the scaling factoring of the circle radius a little:

It is interesting also that this fractal draws in such a way that you can paint alternate cells without painted cells touching each other: I.e., they only touch at the vertices.

This fractal idea did not originate with me, but I wrote some racket code to do the midpoint calculation using barycentric coordinates. This fractal draws a circle at the midpoint of a triangle, then subdivides the triangle and repeats:

Here is the same fractal to four iterations:

To get the midpoints, I could simple pass in the coordinates of the last triangle ABC, and then use “0.5” barycentric coordinates:

[P1 (barycentric->complex 0.5 0.5 0 A B C)]
[P2 (barycentric->complex 0.5 0.0 0.5 A B C)]
[P3 (barycentric->complex 0 0.5 0.5 A B C)]

An interesting concept I came across recently in an introductory geometry book is barycentric coordinates. The fundamental idea is defining the location of a point by its relationship to the points of a triangle. (Well, technically you can use any polygon, but we will keep things simple.)

So, each of the vertices of the triangle (A, B, and C) are given a sort of weight or pull factor. These three numbers uniquely represent the position of P. The theorem goes

If A, B, C are three non-collinear points, any vector P may be expressed in terms of the vectors A, B, C thus:

P = xA + yB +zC, where x + y + z = 1

Geometry: A Comprehensive Course, Dan Pedoe

If P is pulled all the way over to vertex A, then the coordinates become (1, 0, 0). Likewise at B it is (0, 1, 0) and for C it is (0, 0, 1). Another coordinate inside the triangle would be (⅓, ⅓, ⅓). If the point moves outside the triangle, some of those coordinates become negative.

Of course, instead of boring old vectors, you can use… complex numbers!

From the formula, it is obvious how to start with some barycentric coordinates, add a triangle, and find the point. Here is some Racket code:

(define (barycentric->complex x y z A B C)
(+ (* x A) (* y B) (* z C)))

Note, that you can use any triangle you want with the same barycentric coordinates. This could allow you to stretch the graph of some set of coordinates in various ways.

But, how do you convert in the other direction, having a point and a triangle, and wanting the barycentric coordinates? This is a little more complicated. Fortunately, it happens that the factor for each of the vertices is equal to the proportional area of the triangle opposite that vertex. So, looking at this diagram again

x = [BCP] / [ABC], y = [CAP] / [ABC], z = [ABP] / [ABC] where [BCP] is the area of triangle BCP, etc.

Now, a tricky part here is that the areas must be signed areas, or the formula doesn’t necessarily work correctly. Signed area means that you have to be consistent in assigning a negative or positive value to the area of the triangle, depending on the orientation of the vertices. So, if you happen to plug in the vertices in a counter-clockwise direction, you get a negative area, and in a clockwise direction, you get a positive area. Thankfully, there is an elegant formula for determining this, which I will represent by this Racket function

The function returns ‘ccw for counter-clockwise, ‘cw for clockwise, and ‘col for collinear (i.e., the points are on a straight line, which is not relevant in this application).

This allows us to have a signed-area function, using the common formula for determining the area of a triangle from just the length of the sides. We get the length of the sides by subtracting the complex numbers from each other and pulling the magnitudes of the resulting complex numbers.

(define (signed-area A B C)
(let* ([a (magnitude (- C B))]
[b (magnitude (- A C))]
[c (magnitude (- A B))]
[s (* (+ a b c) 0.5)]
[area (sqrt (* s (- s a) (- s b) (- s c)))]
[orient (orientation A B C)])
(if (eq? orient 'ccw)
(* -1 area)
area)))

And finally we determine x, y, and z:

(define (complex->barycentric A B C P)
(let ([ABC (signed-area A B C)])
(list (/ (signed-area B C P) ABC)
(/ (signed-area C A P) ABC)
(/ (signed-area A B P) ABC))))