excursions.tex

\documentclass[a4paper]{scrartcl}
\usepackage{unicode-math}
\usepackage{luamplib}
\usepackage{dwmpcode}
% turn on automatic TeX labels
\mplibtextextlabel{enable}
% add at the start/end of each fig
\everymplib{input colorbrewer-rgb;input invert;interim ahangle := 30;beginfig(0);}
\everyendmplib{endfig;}
%
%\usepackage{shortvrb}\MakeShortVerb{"}
\usepackage{mflogo}
\title{Excursions in \MP}
\author{Toby Thurston}
\date{October 2022 —}
\begin{document}
\maketitle

\noindent
This example document includes geometric illustrations inspired by
\textsl{Excursions in Geometry}, C.\@ Stanley Ogilvy, OUP 1968.
The illustrations are presented roughly in the same order as the book, with notes
about how you can use \MP\ to produce similar.  The section headings also
approximately follow the book.  You might like to read the PDF of
this document side by side with the source code, so that you can see how each
illustration is done.  Each illustration is included as in-line \MP\ code, there are
no external graphics files used.

\section{A bit of background}

Ogilvie starts with a review of some circle theorems.  In fact most of the book is
about circles in one way or another.
In this first diagram, you are given the width $AB$ of the screen and the ideal viewing angle $\theta$.
The \MP\ code works out the rest from that, including a useful routine for a circle through three points.
$$
\begin{mplibcode}
numeric theta;
pair A, B, V, U;
path c;

% given: width of the screen and ideal viewing angle
A = -B = 60 left;
theta = 42;  

% use sine rule to find V
V = sind(90 - 1/2 theta) / sind(theta) * abs (A-B) * dir (1/2 theta - 90) shifted A;

vardef circle_through(expr a, b, c) = 
   save m; pair m; m = whatever * (a-b) rotated 90 shifted 1/2[a,b]
                     = whatever * (a-c) rotated 90 shifted 1/2[a,c];
   fullcircle scaled 2 abs (a-m) shifted m
enddef;

c = circle_through(A, B, V);

U = point 4.8 of c;

draw A -- V -- B -- U -- A;
draw A--B withpen pencircle scaled 1;

clip currentpicture to c; % clip the thick pen
draw c withcolor Blues 8 4;

label.top("\small Screen", 1/2[A, B]);
forsuffixes $ = U, V:
    label("$\theta$", $ + 13 * dir (1/2 theta + angle (B - $))) withcolor Reds 8 7;
endfor
forsuffixes $ = A, B, U, V:
    label(str $, $ shifted - center c scaled 1.08 shifted center c); 
endfor
\end{mplibcode}
$$
Here $\theta$ is half the angle measured by the intercepted arc, which gives us the useful corollary that 
any angle inscribed in a semicircle is a right angle, and conversely that if you can show that some angle $ACB$
    is a right angle, then the semicircle drawn with $AB$ as a diameter must pass through $C$.
$$
\begin{mplibcode}
pair A, B, C; C = origin; B = 150 dir -50; A = 80 dir -140;
% you can do a "semicircle drawn with AB as a diameter" like this
draw halfcircle scaled abs (B-A) 
                rotated angle (B-A) 
                shifted 1/2[B, A] withcolor Blues 8 4;
draw A -- B -- C -- cycle;
forsuffixes $ = A, B, C:
    label(str $, $ shifted - 1/2[A,B] scaled 1.08 shifted 1/2[A,B]); 
endfor
\end{mplibcode}
$$

\smallskip\noindent
A basic theorem: If two chords intersect, the product of the lengths of the segments of 
the one equals the product of the lengths of the segments of the other. In all three cases
below you have $PA/PD = PC/PB$ by similar triangles, hence $PA\cdot PB=PC\cdot PD$.
$$
\begin{mplibcode}
path base; base = fullcircle scaled 120;
picture P[];
P1 = image(
    numeric d; d = 0.2;
    pair A, B, C, D, P;
    A = point 5 of base;
    B = point 2.5 of base;
    C = point 4 + d of base;
    D = point 8 - d of base;

    P = whatever[A, B] = whatever[C, D];

    draw C -- A -- D -- B dashed evenly scaled 1/2 
        withpen pencircle scaled 1/4
        withcolor 3/4;
    draw C -- D;
    draw A -- B;

    forsuffixes $ = A, B, C, D:
        label("$" & str $ & "$", $ scaled 1.1); 
    endfor
    label.lrt("$P$", P);

    drawoptions(withcolor Reds 7 6);
    label("$\scriptscriptstyle\alpha$", B shifted 4 (unitvector(A-B) + unitvector (D-B)));
    label("$\scriptscriptstyle\alpha$", C shifted 4 (unitvector(A-C) + unitvector (D-C)));
    label("$\scriptscriptstyle\beta$",  P shifted 4 (unitvector(A-P) + unitvector (C-P)));
    label("$\scriptscriptstyle\beta$",  P shifted 4 (unitvector(B-P) + unitvector (D-P)));
    drawoptions();
    draw base withcolor Blues 8 4;
    label.urt("I.", lrcorner currentpicture) withcolor Purples 7 7;
);
P2 = image(
    numeric d; d = 0.9;
    pair A, B, C, D, P;
    A = point 3.8 of base;
    B = point 2 of base;
    C = point 4 + d of base;
    D = point 8 - d of base;

    P = whatever[A, B] = whatever[C, D];

    draw B -- (1 + 12 / abs(B-P))[B, P];
    draw D -- (1 + 12 / abs(D-P))[D, P];

    label.ulft("$A$", A);
    label.top ("$B$", B);
    label.llft("$C$", C);
    label.lrt ("$D$", D);
    label.ulft("$P$", P);

    draw base withcolor Blues 8 4;
    label.urt("II.", lrcorner currentpicture) withcolor Purples 7 7;
);
P3 = image(
    numeric d, t; d = 0.7; t = 2.6;
    pair A, B, C, D, P;
    A = B = point t of base;
    C = point 4 + d of base;
    D = point 8 - d of base;

    P = whatever[B, B shifted direction t of base] = whatever[C, D];

    draw B -- (1 + 14 / abs(B-P))[B, P];
    draw D -- (1 + 14 / abs(D-P))[D, P];

    label.ulft("$A,B$", A);
    label.llft("$C$", C);
    label.lrt ("$D$", D);
    label.ulft("$P$", P);
    draw base withcolor Blues 8 4;
    label.urt("III.", lrcorner currentpicture) withcolor Purples 7 7;
);

interim labeloffset := 34;
label.ulft(P1, origin);
label.urt (P2, origin);
label.bot (P3, origin);
\end{mplibcode}
$$
\centerline{This example also shows how to arrange sub-figures into one.}
\clearpage
\section{Harmonic division and Apollonian circles}

Can we find $C$ and $D$ on the line $AB$ so that $AC/CB = AD/BD$? 

\medskip\noindent
Yes:
$$
\begin{mplibcode}
pair A, B, C, D;
A = origin;
numeric a, b, c;
a = 150; b = 50; c * a / b = a + b + c;
B = (a+b, 0);
C = (a, 0);
D = (a+b+c, 0);

draw A--D;
forsuffixes $ = A, B, C, D:
    draw (up--down) scaled 2 shifted $;
    label.bot("$" & str $ & "$", $);
endfor
\end{mplibcode}
$$

\bigskip\noindent
\textit{Theorem}. The bisector of any angle of a triangle divides
the opposite side into parts proportional to the adjacent sides.

\smallskip\noindent
So given $P$, with $AP=k$ and $BP=1$:
$$ % Fig 8 in the book
\begin{mplibcode}
pair A, B, C, D, P;
A = origin;
numeric a, b, c, k;
a = 180; b = 50; c * a / b = a + b + c;
B = (a+b, 0);
C = (a, 0);
D = (a+b+c, 0);
k = a / b;
path aa, bb; 
aa = fullcircle scaled 144k shifted A;
bb = fullcircle scaled 144  shifted B;
P = aa intersectionpoint bb;

draw C -- P -- D dashed evenly scaled 1/2 
                 withpen pencircle scaled 1/4
                 withcolor 3/4;
draw A -- (1+26/abs(A-P))[A, P];
draw B -- P;
draw A--D;
forsuffixes $ = A, B, C, D:
    dotlabel.bot("$" & str $ & "$", $);
endfor
dotlabel.ulft("$P$", P);
label.ulft("$k$", 1/2[A, P]);
label.rt("$1$", 1/2[B, P]);
drawoptions(withcolor Reds 7 6);
label("$\scriptscriptstyle\alpha$", P shifted 8 (unitvector(A-P) + unitvector (C-P)));
label("$\scriptscriptstyle\alpha$", P shifted 8 (unitvector(C-P) + unitvector (B-P)));
label("$\scriptscriptstyle\beta$",  P shifted 6 (unitvector(B-P) + unitvector (D-P)));
label("$\scriptscriptstyle\beta$",  P shifted 6 (unitvector(P-A) + unitvector (D-P)));
drawoptions();
\end{mplibcode}
$$
We have $AC/CB = AP/BP = k$ from the interior angles and $AD/BD = AP/BP = k$ from the exterior pair.
The proof looks like this:
$$
\begin{mplibcode}
pair A, B, C, D, P, E, F, G, H;
A = origin;
numeric a, b, c, k;
a = 180; b = 50; c * a / b = a + b + c;
B = (a+b, 0);
C = (a, 0);
D = (a+b+c, 0);
k = a / b;
path aa, bb; 
aa = fullcircle scaled 144k shifted A;
bb = fullcircle scaled 144  shifted B;
P = aa intersectionpoint bb;

E = whatever[A, P]; E - B = whatever * (P - C);
F = whatever[B, E] = whatever[P, D];
G = whatever[A, P]; G - B = whatever * (P - D);
H = whatever[B, G] = whatever[P, C];

draw unitsquare scaled 5 rotated angle (E-F) shifted F 
                 withpen pencircle scaled 1/4
                 withcolor 3/4;
draw unitsquare scaled 5 rotated angle (P-H) shifted H 
                 withpen pencircle scaled 1/4
                 withcolor 3/4;
draw C -- P -- D dashed evenly scaled 1/2 
                 withpen pencircle scaled 1/4
                 withcolor 1/2;
draw G -- B -- E dashed evenly scaled 1/2 
                 withpen pencircle scaled 1/4
                 withcolor 1/2;
draw A -- (1+20/abs(A-E))[A, E];
draw B -- P;
draw A--D;
forsuffixes $ = A, B, C, D:
    dotlabel.bot("$" & str $ & "$", $);
endfor
forsuffixes $ = E, G, P:
    dotlabel.ulft("$" & str $ & "$", $);
endfor
label.rt("\ $F$",F);
drawoptions(withcolor Reds 7 6);
label("$\scriptscriptstyle\alpha$", P shifted 8 (unitvector(A-P) + unitvector (C-P)));
label("$\scriptscriptstyle\alpha$", P shifted 8 (unitvector(C-P) + unitvector (B-P)));
label("$\scriptscriptstyle\beta$",  P shifted 6 (unitvector(B-P) + unitvector (D-P)));
label("$\scriptscriptstyle\beta$",  P shifted 6 (unitvector(P-A) + unitvector (D-P)));
drawoptions();
\end{mplibcode}
$$
Since the lines $PC$ and $PD$ are bisectors, then we have $2\alpha + 2\beta = 180°$,
hence $\alpha+\beta$ is a right angle, so if we draw $BE$ parallel to $PC$ it will
be perpendicular to $PD$, and hence the two triangles $PFE$ and $PFB$ are congruent,
so that $PE=PB$.  But the parallel lines cut off proportional segments, so that 
$$
\frac{AC}{CB} = \frac{AP}{PE} = \frac{AP}{BP}
$$
as required.  And using the exterior angle
$$
\frac{AD}{BD} = \frac{AP}{GP} = \frac{AP}{BP}
$$

\clearpage
\subsection{Applied Apollonian example}

Given $A$, $B$, $k$, and course of the target, find intersection point $Q$.
$$
\begin{mplibcode}
pair A, B, C, D, P, Q;
A = origin; B = 180 right; k = 2;
C = (k/(k+1))[A, B];
D = (k/(k-1))[A, B];
path app; 
app = fullcircle scaled length (C-D) shifted 1/2[C,D];
P = point 2.236 of app;
Q = app intersectionpoint (B -- B + 1024 dir -60);
%Q = point 6 of app;
draw app withcolor Blues 7 4;
draw A -- D -- P -- cycle;
draw B -- P -- C;
numeric r; r = 3/8;
draw r[B, Q] -- Q -- r[A, Q] withpen pencircle scaled 1/4 dashed evenly;
drawarrow B -- r[B, Q];
drawarrow A -- r[A, Q];
% label.ulft("$" & decimal (length (P-A) / length (P-B)) & "$", 1/2[A, P]);
% numeric alpha; alpha = angle (Q-A);
% label(decimal alpha, A + 30 dir 1/2 alpha);
label.lft ("$A$", A);
label.ulft("$B$", B);
label.ulft("$C$", C);
label.rt  ("$D$", D);
label.top ("$P$", P);
label.lrt ("$Q$", Q);
    drawoptions(withcolor Reds 8 7);
       label.top("$a$", 1/2[A,C]); 
       label.top("$b$", 1/2[B,C]); 
       label.top("$c$", 1/2[B,D]); 
    drawoptions();
\end{mplibcode}
$$
But how can you find $C$ and $D$ in \MP\ given $A$, $B$, and $k$?  Using
\textcolor{textred}{$a$}, 
\textcolor{textred}{$b$}, and 
\textcolor{textred}{$c$} 
for the lengths as shown above, then using the mediation syntax,
point $C$ is \mpl{(a/(a+b))[A, B]} and $D$ is \mpl{((a+b+c)/(a+b))[A, B]}, so what
are these fractions in terms of $k$?
We know that $AC/CB=AD/BD=k$ so we have 
$a/b = (a+b+c)/c = k$ which we can develop as follows.

Starting with $a/b=k$, add 1 to each side, $a/b + b/b = k + 1$, hence $(a+b)/b = k +
1$.  We can then divide the first equation by the last to get $a/(a+b) = k/(k+1)$.

Similarly, starting with $(a+b+c)/c = k$, take 1 away from each side to get
$(a+b+c)/c - c/c = k - 1$, hence $(a+b)/c = k-1$.  And again dividing the first by
the last gives $(a+b+c)/(a+b) = k / (k-1)$. 

\vfill\noindent
So our \MP\ expressions for points $C$ and $D$ can be written like this:
\begin{code}
    C = (k/(k+1))[A, B]; D = (k/(k-1))[A, B];
\end{code}
If $k=1$ you will get a “Division by zero” error from the second expression.
This corresponds to $D$ at infinity and $C$ half way from $A$ to $B$.  The circle
through these two points can be thought of as the perpendicular bisector of $A$ and
$B$. If you have $0 < k < 1$ then the circle will be round $A$ not $B$.  If $k<0$
the positions of $C$ and $D$ will be swapped, and there is another error when
$k=-1$.

\subsection{Coaxial families}
$$
\begin{mplibcode}

pair A, B; A = -B = 42 left;
draw (left -- right) scaled 400;
dotlabel.bot("$A$", A); 
dotlabel.bot("$B$", B); 
drawoptions(withcolor Blues 7 4);
for i = -4 upto 4:
    if i = 0:
        draw (up--down) scaled 400;
    else:
        numeric k; pair C, D; path ac;
        k = 2**(i/2);
        C = (k/(k+1))[A, B]; D = (k/(k-1))[A, B];
        ac = fullcircle scaled abs(C-D) shifted 1/2[C,D];
        draw ac; 
    fi
endfor
drawoptions();
numeric wd; wd = \mpdim\textwidth;
clip currentpicture to unitsquare shifted -(1/2, 1/2) scaled wd yscaled 5/8; 
\end{mplibcode}
$$
The circles on the left are drawn with $0<k<1$ and those on the right with $k>1$.
The circles drawn with a given $k$ is the mirror image of one drawn with $1/k$; so the smallest
circle around $A$ is drawn with $k=1/4$ and the corresponding smallest circle around $B$ has $k=4$.

\vfill\noindent
The tangents at the intersections of orthogonal circles are at
right angles to each other.  
$$
\begin{mplibcode}
path c[]; numeric r[]; pair o[]; t = 7/8;
r1 = 64; r2 = 50; o1 = origin;
c1 = fullcircle scaled 2r1 shifted o1;
o2 = -r2 * unitvector (direction t of c1) shifted point t of c1;
c2 = fullcircle scaled 2r2 shifted o2;
draw o1 -- o2 -- point t of c1 -- cycle withpen pencircle scaled 1/4;
draw c1 withcolor Blues 8 6;
draw c2 withcolor Oranges 8 6;
dotlabel.bot("$O_1$", o1);
dotlabel.lrt("$O_2$", o2);
label.ulft("$r_1$", 1/2[o1, point t of c1]);
label.urt("$r_2$", 1/2[o2, point t of c1]);
\end{mplibcode}
$$
Given $O_1$, $r_1$, and $r_2$, you can find $O_2$ in \MP\ with:
\begin{code}
c1 = fullcircle scaled 2 r1 shifted o1;
o2 = unitvector (direction t of c1) scaled r2 shifted point t of c1;
\end{code}
where $t$ is the desired time on $C_1$ for a point of intersection.

\newpage
If you consider a circle with radius $r$ and diameter $AB$, where $AB$ is divided
harmonically then you have (with red braces for the labels): 
$$
\begin{mplibcode}
numeric r, k; k = 3.1415; r = 84;
pair A, B, C, D; 
-A = B = (r, 0);
C = (k/(k+1))[A, B]; 
D = (k/(k-1))[A, B];
draw subpath (-1/2, 1/2) of fullcircle scaled 2r withcolor Blues 8 3;
draw subpath (4-1/2, 4+1/2) of fullcircle scaled 2r withcolor Blues 8 3;
draw A--D;
dotlabel.bot("$A$", A); 
dotlabel.bot("$B$", B); 
dotlabel.bot("$C$", C); 
dotlabel.bot("$D$", D); 
dotlabel.bot("$O$", origin); 

vardef on_brace(expr a,b,m,r,s) =
    save d, e, n, bb; numeric d, n; pair e; path bb;
    n = 1/2 m; d = angle (b-a);
    e = up scaled m rotated d shifted r[a,b];
    bb = (
          (origin {up} .. {right} (abs n,n)) rotated d shifted a 
          --
          ((-abs n,-n){right} .. {up} origin {down} .. {right}(abs n,-n)) rotated d shifted e 
          --
          ((-abs n,n){right} .. {down} origin) rotated d shifted b
         ) shifted (up scaled n rotated d);
    draw bb withpen pencircle yscaled .6 xscaled .1666 rotated d withcolor s;
    point 3 of bb
enddef;

label.top("$r$", on_brace(A + up, up shifted 1/4 left, 4, 1/2, Reds 7 6));
label.top("$r$", on_brace(up shifted 1/4 right, up + B, 4, 1/2, Reds 7 6));

\end{mplibcode}
$$
and so you can write
$
\frac{AC}{CB}=\frac{AD}{BD}$ as $\frac{r+OC}{r-OC} = \frac{OD+r}{OD-r}$
and hence
$r^2 = OC\cdot OD$.  And this expression is therefore \textit{equivalent} to the
statement that the division is harmonic, and if a circle is drawn with $AB$ as a
diameter, \textit{any} circle drawn through $C$ and $D$ will be orthogonal.
$$
\begin{mplibcode}
vardef through@#(expr a, b) = save d; numeric d; 
    d = abs(b-a); (1+@#/d)[b,a] -- (1+@#/d)[a,b]
enddef;
vardef circle_through(expr a, b, c) = save o; pair o;
    o = whatever * (b-a) rotated 90 shifted 1/2[a,b] 
      = whatever * (c-b) rotated 90 shifted 1/2[b,c];
    fullcircle scaled 2 abs(a-o) shifted o
enddef;
numeric r, k, q; k = 3.1415; r = 84;
pair A, B, C, D, T, U, X, Y; 
-A = B = (r, 0);
C = (k/(k+1))[A, B]; 
D = (k/(k-1))[A, B]; 
path alpha, beta;
alpha = fullcircle scaled 2r rotated 54; 
T = point 0 of alpha;
beta = circle_through(T, C, D);
U = alpha intersectionpoint beta;


draw origin -- alpha withcolor Blues 8 6;
draw beta withcolor Oranges 8 6;
draw through10(A,D);
dotlabel.ulft("$A$", A); 
dotlabel.urt ("$B$", B); 
dotlabel.urt ("$C$", C); 
dotlabel.ulft("$D$", D); 
dotlabel.top("\strut $T$", T); 
dotlabel.bot("\strut\ $U$", U); 
dotlabel.bot("$O$", origin); 

label.ulft("$r$", 1/2 T);
label.ulft("$\alpha$", point 1.5 of alpha);
label.urt("$\beta$", point 1 of beta);


\end{mplibcode}
$$
Two useful routines are used in this drawing.
\begin{itemize}
    \item Return the path of a line through two points, with a small overlap at each end.
\begin{code}
vardef through(expr a, b, o) = save d; numeric d; 
    d = abs(b-a); (1+o/d)[b,a] -- (1+o/d)[a,b]
enddef;
\end{code}
\item Return the path of the (unique) circle through three points.
\begin{code}
vardef circle_through(expr a, b, c) = save o; pair o;
    o = whatever * (b-a) rotated 90 shifted 1/2[a,b] 
      = whatever * (c-b) rotated 90 shifted 1/2[b,c];
    fullcircle scaled 2 abs(a-o) shifted o
enddef;
\end{code}
\end{itemize}
\newpage
\noindent
Intersecting orthogonal circles. \count10=\inputlineno
\\\noindent
\begin{mplibcode}
pair C, D; C = -D = 42 left; 
for i=-5 upto 5:
    if i <> 0:
        pair A, B; numeric k; k = 2**(i/2);
        A = (k/(k+1))[C, D]; 
        B = (k/(k-1))[C, D]; 
        draw fullcircle scaled abs(A-B) shifted 1/2[A,B] withcolor Blues 8 6;
        pair O; O = 1/2[C,D] shifted (0, 8 i * 1.6 sqrt(abs(i)));
        draw fullcircle scaled 2 abs(C-O) shifted O withcolor Oranges 8 6;
    fi
endfor
draw (left--right) scaled 300 withcolor Oranges 8 7;
draw (down--up)    scaled 200 withcolor Blues 8 7;
draw C withpen pencircle scaled dotlabeldiam;

numeric wd; wd = \mpdim\textwidth;
clip currentpicture to unitsquare shifted -(1/2, 1/2) scaled .8 wd yscaled 3/4;
\end{mplibcode}
\count11=\inputlineno
\advance\count10 3
\advance\count11 -3\par
\mpexternal[firstline=\the\count10,lastline=\the\count11]{\jobname.tex}

\vfill
\noindent Radical axis of non-intersecting circles
$$
\begin{mplibcode}

    pair A, B, C, P, R, S, T, U;
    path a, b, c;

    z0 = origin;
    z1 = -4/3(21, 10);
    z2 = 110 dir -89;
    a = fullcircle scaled 233 shifted z0;
    b = fullcircle scaled  99 shifted z1;
    c = fullcircle scaled 164 shifted z2;

    T = a intersectionpoint c; U = reverse a intersectionpoint c;
    R = b intersectionpoint c; S = reverse b intersectionpoint c;
    P = whatever[T, U] = whatever[R, S];

    A = a intersectionpoint halfcircle rotated angle (z0-P) scaled abs(z0-P) shifted 1/2[z0, P];
    B = b intersectionpoint halfcircle rotated angle (z1-P) scaled abs(z1-P) shifted 1/2[z1, P];
    C = c intersectionpoint halfcircle rotated angle (z2-P) scaled abs(z2-P) shifted 1/2[z2, P];

    vardef through(expr a, b) = save d; numeric d; d = abs(a-b); (1+16/d)[b,a] -- (1+16/d)[a,b] enddef;

    drawoptions(withcolor Blues 8 4);
    draw quartercircle rotated -17 scaled 2 abs(A-P) shifted P;
    draw through(P, A);
    draw through(P, B);
    draw through(P, C);

    drawoptions(withcolor Oranges 8 3);
    draw c; label.urt("$c$", point 1 of c);
    draw through(P, U);
    draw through(P, S);

    drawoptions();
    draw a; label.urt("$a$", point 1 of a);
    draw b; label.urt("$b$", point 1 of b);
    draw (left--right) scaled 100 rotated (90 + angle (z1-z0)) shifted P;

    drawoptions(withcolor Oranges 8 5);
    dotlabel.llft("$T$", T);
    dotlabel.lrt ("\;$U$", U);
    dotlabel.bot ("\strut$R$", R);
    dotlabel.urt ("\vrule width 0pt depth 2pt height 0pt\,$S$", S);

    drawoptions(withcolor Blues 8 7);
    dotlabel.urt("\ $A$", A);
    dotlabel.top("\strut $B$", B);
    dotlabel.urt("\strut $C$", C);

    drawoptions();
    dotlabel.bot ("\qquad $P$", P);

    currentpicture := currentpicture scaled 0.9; 
   
\end{mplibcode}
$$

%----------------------------
\section{Inversive geometry}

\subsection{Transformations}

Given a circle of radius $r$, let each point $C$ inside the circle
be transformed to a point $D$ outside the circle in such a way that the 
distances comply with $OC\cdot OD=r^2$ and the line $OC$ passes through $D$.
$$
\begin{mplibcode}
numeric r; pair C, D, E;
r = 78;     
C = (9/16 r, 0) rotated 42;
numeric s; s = r / abs C;
D = C * s * s;
draw D -- origin -- fullcircle scaled 2r;
label.bot("$r$", 1/2(r, 0));
dotlabel.lft("$O$", origin);
dotlabel.ulft("$C$", C);
dotlabel.rt("$D$", D);
\end{mplibcode}
$$
For \MP\ the direct way to find $D$ given $O$, $C$, and $r$ is to calculate
\begin{code}
D = O + unitvector(C-O) scaled r * r / abs (C-O);
\end{code}
although if you keep $O$ at the origin this can be simplified to 
\begin{code}
D = unitvector(C) scaled r * r / abs C;
\end{code}
but with the default number system this risks an overflow if $r > \sqrt{2^{15}} \approx 181$, so some  
    safer approach is more helpful for larger circles.  Examining \mpl{plain.mp} shows that 
\mpl{unitvector} is a macro defined like this:
\begin{code}
vardef unitvector primary z = z/abs z enddef;
\end{code}
which suggests this alternative (safer) approach for inverting $C$ to find $D$:
\begin{code}
numeric s; s = r / abs C; D = C * s * s; 
\end{code}
This works well provided that $|C| > \frac1{180} r$, which is usually the case.
A more general macro to invert a point $p$ in the circle centred at $o$ with radius $r$
might be:
\begin{code}
vardef invert(expr p, o, r) = 
    save t; pair t; t = p - o; 
    save s; numeric s; s = r / abs t;
    o + t * s * s
enddef;
\end{code}
You could also consider checking that $|t|>0$ and that
$s$ was not too large.

\newpage
There are other ruler-and-compass constructions to find an inverse point, 
but they do not always lead to more efficient methods in \MP.

$$
\begin{mplibcode}
primarydef a past b = a -- (1+12/abs(a-b))[a,b] enddef;
picture M[];
M1 = image(
path c; pair P, Q, A, B, P';    
c = fullcircle scaled 120;
P = 20 left rotated 3;
A = point 1/45 angle P of c;
B = point 1/45 angle (P-A) of c;
Q = point 1.8 of c;
P' = whatever[A, B] = whatever[Q, P reflectedabout(A, Q)];

draw c withcolor Blues 8 6;
draw A -- Q -- B;
draw B past P';
draw Q past P' dashed evenly scaled 1/2 withpen pencircle scaled 1/4;
draw P -- Q;

draw origin withpen pencircle scaled dotlabeldiam;
undraw origin withpen pencircle scaled 1/2 dotlabeldiam;
dotlabel.bot ("$P$", P);
dotlabel.bot ("$P'$", P');
dotlabel.llft("$A$", A);
dotlabel.rt  ("$B$", B);
dotlabel.urt ("$Q$", Q);

drawoptions(withcolor Oranges 8 8);
label("$\scriptscriptstyle 1$", Q shifted 21 dir 1/2(angle (P-Q) + angle (A-Q)));
label("$\scriptscriptstyle 2$", Q shifted 21 dir 1/2(angle (P'-Q) + angle (A-Q)));
drawoptions());
M2 = image(
path c; pair A, B, P, P', Q, R;
c = fullcircle scaled 120 rotated 3;
A = point 4 of c; B = point 0 of c; 
numeric t; t = 1.5;
Q = point t of c;
R = point -t of c;
P = whatever[A, B] = whatever [Q, R];
P' = whatever[A, B] = whatever * direction t of c shifted point t of c;

draw c withcolor Blues 8 6;
draw A past P';
draw Q past P';
draw Q -- R;
draw A -- Q -- B dashed evenly scaled 1/2 withpen pencircle scaled 1/4;
draw origin withpen pencircle scaled dotlabeldiam;
undraw origin withpen pencircle scaled 1/2 dotlabeldiam;
dotlabel.urt("$P$", P);
dotlabel.urt ("$P'$", P');
dotlabel.llft("$A$", A);
dotlabel.lrt ("$B$", B);
dotlabel.urt ("$Q$", Q);
dotlabel.lrt ("$R$", R);
drawoptions(withcolor Oranges 8 8);
label("$\scriptscriptstyle 1$", Q shifted 17 dir 1/2(angle (P-Q) + angle (B-Q)));
label("$\scriptscriptstyle 2$", Q shifted 17 dir 1/2(angle (P'-Q) + angle (B-Q)));
drawoptions());
label.urt(M1, -(46, 21));
label.llft(M2, (46, 21));
\end{mplibcode}
$$
The first method here is to construct $\angle2$ at $Q$, such that $\angle 2 = \angle 1$, and
the dotted line cuts $AB$ in $P'$.  In \MP\ this is 
\begin{code}
P' = whatever[A, B] = whatever[Q, P reflectedabout(A, Q)];
\end{code}
The second is to construct $QR$ perpendicular to the diameter through $P$, and then
the tangent at $Q$ cuts $AB$ in $P'$ which, given the circle $c$ centred 
at the origin, and the point $P$,
translates into something like this:
\begin{code}
numeric t; (t, whatever) = c intersectiontimes 
                          P -- 200 unitvector(P) rotated 90 shifted P;
P' = whatever[A,B] = whatever * direction t of c shifted point t of c;
\end{code}
\vfill\noindent
Drawing a circle through $Q$, $P$, and $P'$ shows why this works.
The tangent to one circle is the radius of the other, so $AB$ is divided 
harmonically by $P$, and $P'$.
$$\begin{mplibcode}
primarydef a past b = a -- (1+12/abs(a-b))[a,b] enddef;
path c; pair A, B, P, P', Q, R;
c = fullcircle scaled 100 rotated 3;
A = point 4 of c; B = point 0 of c; 
numeric t; t = 1.4;
Q = point t of c;
R = point -t of c;
P = whatever[A, B] = whatever [Q, R];
P' = whatever[A, B] = whatever * direction t of c shifted point t of c;
draw fullcircle scaled abs(Q-P') shifted 1/2[Q, P'] withcolor Oranges 8 6;
draw c withcolor Blues 8 6;
draw A past P';
draw Q past P';
draw Q -- P;
draw origin withpen pencircle scaled dotlabeldiam;
undraw origin withpen pencircle scaled 1/2 dotlabeldiam;
draw   1/2[Q,P'] withpen pencircle scaled dotlabeldiam;
undraw 1/2[Q,P'] withpen pencircle scaled 1/2 dotlabeldiam;
dotlabel.urt("$P$", P);
dotlabel.urt ("$P'$", P');
dotlabel.llft("$A$", A);
dotlabel.lrt ("$B$", B);
dotlabel.urt ("$Q$", Q);
\end{mplibcode}$$

If the point is outside the circle, you need to construct the tangent
first, then the circle $b$ cuts $OP$ at the inverse point $P'$:
$$\begin{mplibcode}
path a, b, c; pair O, P, P', Q, A;
O = origin;
c = fullcircle scaled 144 rotated 3;
P = 1.732 point 0 of c;
a = halfcircle zscaled (P-O) shifted 1/2[P,O];
Q = a intersectionpoint c;
b = fullcircle zscaled (Q-P) shifted 1/2[Q,P];
%P' = whatever[origin, P]; Q-P' = whatever * (origin-P) rotated 90;
P' = b intersectionpoint (center c -- point 0 of c);

A = point 4 of c; numeric d; d = abs(P-A);
draw (1+12/d)[P, A] -- (1+20/d)[A, P];
draw origin -- Q -- P;
draw c withcolor Blues 8 6; label.ulft("$c$", point 3 of c);
draw b withcolor Oranges 8 6; label.urt("$b$", point 11/2 of b);
draw a dashed evenly scaled 1/2 withpen pencircle scaled 1/4;

dotlabel.bot("$O$", origin);
dotlabel.lrt("$P$", P);
dotlabel.llft("$P'$", P');
dotlabel.top("\strut$Q\;$", Q);
\end{mplibcode}$$
There is a neat \MP\ idiom to draw the semi-circle between $O$ and $P$
\begin{code}
path a; a = halfcircle zscaled (P-O) shifted 1/2[P,O];
\end{code}
and if $O$ is really the origin, you can abbreviate that to
\begin{code}
path a; a = halfcircle zscaled P shifted 1/2 P;
\end{code}
Then, assuming path $c$ is the blue circle above, point $Q$ is
\begin{code}
pair Q; Q = a intersectionpoint c;
\end{code}
and point $P'$ is:
\begin{code}
path b; b = fullcircle scaled abs(Q-P) shifted 1/2[Q, P];
pair P'; P' = b intersectionpoint (O--P);
\end{code}
although you could do without the circle $b$ by using the normal 
idiom to find the closest point to $Q$ on $OP$:
\begin{code}
pair P'; P' = whatever[O, P]; Q-P' = whatever * (P-O) rotated 90;
\end{code}
This is all interesting geometry, but if you just want to find the inverse point
then the \mpl{invert} macro given earlier in this section also works well for points
    outside the circle.

\vfill\noindent
\hbox to \textwidth{
    $\vcenter{\begin{mplibcode}
path c; c = fullcircle scaled 120 rotated 3;
pair P; P = 1.32 point 0 of c;

path a; a = fullcircle zscaled 2P shifted P;
numeric t, u;
(t, u) = a intersectiontimes c;
pair S; S = point u of c;
path b; b = reverse fullcircle zscaled 2S shifted S;
numeric v, w;
(v, w) = b intersectiontimes (origin -- P);
pair P'; P' = point v of b;

pair A; A = point 4 of c; numeric d; d = abs(P-A);
draw (1+12/d)[P, A] -- (1+20/d)[A, P];
draw subpath(t-1/4, 4+1/4) of a withcolor Oranges 8 4;
draw subpath(v-1/4, 4+1/4) of b withcolor Oranges 8 4;
draw P -- S -- P' withcolor Oranges 8 3;
draw S -- origin withcolor Oranges 8 3;
draw c withcolor Blues 8 6;
dotlabel.bot("$P$", P);
dotlabel.bot("$P'$", P');
dotlabel.llft("$O$", origin);
dotlabel.top("$S$", S);
label.rt("$r$", 1/2 S);

\end{mplibcode}}$ \hss
$\vcenter{\vbox{\hsize 3in\noindent
Similarly, if you are working with ruler-and-compasses, you can
    find the inverse of an external point $P$ by drawing two arcs as shown, \llap{$\leftarrow$\ }on the left.
But this is less convenient in \MP, primarily because you need to make sure 
that you find the correct intersections of the two arcs, so the approach is less
“robust” than doing the calculations. The construction works because we have 
$OP/r = r/OP'$, hence $OP\cdot OP' = r^2$ as before.
}}$}

\subsection{Invariants}

\textbf{Inversion of a straight line outside the circle of inversion}.
$$\begin{mplibcode}
input invert
path c; pair P, Q, P', Q'; numeric r; 
r = 68;
c = fullcircle scaled 2r;
P = 1.414 point 0 of c;
xpart Q = xpart P;  Q-P = .9r * up;
P' = invert(P, origin, r);
Q' = invert(Q, origin, r);

draw unitsquare scaled 5 rotated angle (Q-P) shifted P withcolor 3/4;
draw unitsquare scaled 5 rotated (angle (P'-Q')-90) shifted Q' withcolor 3/4;

draw origin -- P;
draw origin -- Q;
draw P' -- Q';
draw 1.8[Q,P] -- 1.2[P, Q] withcolor Oranges 8 5;

draw origin -- point 2.8 of c dashed evenly scaled 1/2; label.urt("$r$", 1/2 point 2.8 of c);
draw c withcolor Blues 8 6;
draw fullcircle zscaled P' shifted 1/2 P' withcolor Oranges 8 5;

dotlabel.llft("$O$", origin);
dotlabel.rt  ("$P$", P);
dotlabel.rt  ("$Q$", Q);
dotlabel.lrt ("$P'$", P');
dotlabel.top ("$Q'$", Q');
\end{mplibcode}$$
Here we have $OP\cdot OP' = r^2$ and $OQ\cdot OQ' = r^2$ so $OP/OQ = OP'/OQ'$, hence
the triangles $OPQ$ and $OQ'P'$ are similar and $\angle OQ'P'$ must always be a
right angle, since $OP$ is perpendicular to the line, but $Q$ is any point on the
line, so the line must invert to a circle through $O$.  

\smallskip\noindent
\textbf{Inversion of a line that cuts the circle of inversion}.
$$\begin{mplibcode}
path c; pair P, P', A, B; numeric r; 
r = 68;
c = fullcircle scaled 2r rotated 4;
P = 0.8 point 0 of c;
P' = invert(P, origin, r);
path line; line = (down--up) rotated angle P scaled 1.2 r shifted P;
draw c withcolor Blues 8 6;
draw fullcircle zscaled P' shifted 1/2 P' withcolor Oranges 8 5;
draw line withcolor Oranges 8 5;

dotlabel.llft("$O$", origin);
dotlabel.llft("$P$", P);
dotlabel.lrt ("$P'$", P');
numeric t, u;
(t, whatever) = c intersectiontimes subpath(1/2, 1) of line; dotlabel.urt("$A$", point t of c);
(u, whatever) = c intersectiontimes subpath(0, 1/2) of line; dotlabel.lrt("$B$", point u of c);
\end{mplibcode}$$
Note that the points $A$ and $B$ are invariant.  Using the \mpl{invert} macro
already described, you can find $P'$ directly and draw the inverted line as
a circle with $OP'$ as a diameter.  But it is also possible to find the points $A$ and $B$
as the points where the line through $P$ that is perpendicular to $OP$ cuts the circle
of inversion, and then draw a circle through $O$, $A$, and $B$ using a routine like this:
\begin{code}
vardef circle_through(expr A, B, C) = save o; pair o;
    o = whatever * (A-B) rotated 90 shifted 1/2 [A,B]
      = whatever * (B-C) rotated 90 shifted 1/2 [B,C];
    fullcircle scaled 2 abs (A-o) shifted o enddef;
\end{code}
then $P'$ is the intersection of this circle with the line through $OP$.

\smallbreak\noindent
\textbf{Inversion of a circle outside the circle of inversion}
$$\begin{mplibcode}
path c, d; 
pair P, P', Q, Q', R, R'; 
numeric r; 
r = 100; z0 = origin; z1 = 1.9 r * dir 12;
c = fullcircle scaled 2r;
d = fullcircle scaled 1.414r rotated 12 shifted z1;

numeric t, a; 
(t, whatever) = d intersectiontimes (z0 -- 8[z0,z1]);
(a, whatever) = d intersectiontimes halfcircle zscaled z1 shifted 1/2 z1;

Q = point t of d;
R = point a of d;
P = point t + 4 of d;

Q' = invert(Q, z0, r);
R' = invert(R, z0, r);
P' = invert(P, z0, r);

draw c withcolor Blues 8 6;
draw d withcolor Oranges 8 5;
draw fullcircle zscaled (P' - Q') shifted 1/2[P', Q'] withcolor Oranges 8 5;

draw P -- R -- Q withpen pencircle scaled 1/4;
draw P' -- R' -- Q' withpen pencircle scaled 1/4;
draw R -- z0 -- Q;

for i=0 step 1 until 7.5:
    draw point i of d withpen pencircle scaled 1;
    draw invert(point i of d, z0, r) withpen pencircle scaled 1;
endfor

forsuffixes $ = Q', P: dotlabel.llft("\strut$" & str $ & "$", $); endfor
forsuffixes $ = R, R': dotlabel.ulft("$" & str $ & "$", $); endfor
dotlabel.rt("$Q$", Q);
dotlabel.lft("$O$", z0);
dotlabel.lrt("$P'$", P'); 
\end{mplibcode}$$
Similar triangles again show that $\angle P'R'Q'$ is a right angle, which
proves that the inversion of the circle is also a circle.  It is sufficient
therefore to find $P$ and $Q$ as the diameter of the circle-to-be-inverted that
passes through $O$ (the centre of the circle of inversion), and then the image is the 
circle drawn with the inverted points $P'$ and $Q'$ as a diameter.
However you sould beware that the points of the source circle are not mapped
equally to the image. The small dots try to show this in the picture above.
With source circle $s$, and a circle of inversion at $z_0$ with radius $r$, 
you can invert $s$ to path $s'$ like this:
\begin{code}
s' = for i=0 upto 7: invert(point i of s, z0, r) .. endfor cycle;
\end{code}
The resulting path should look circular but will run in the opposite direction
to the source path (as you might expect since inversion is a form of reflection).
Note also that the triangle $PQR$ does \textit{not} invert to the triangle $P'Q'R'$, 
and the midpoint of $P$ and $Q$ is not the midpoint of $P'$ and $Q'$.
$$\begin{mplibcode}
path c, d; 
pair P, P', Q, Q', R, R'; 
numeric r; 
r = 100; z0 = origin; z1 = 1.9 r * dir 12;
c = fullcircle scaled 2r;
d = fullcircle scaled 1.414r rotated 12 shifted z1;

numeric t, a; 
(t, whatever) = d intersectiontimes (z0 -- z1);
(a, whatever) = d intersectiontimes halfcircle zscaled z1 shifted 1/2 z1;

P = point t of d;      P' = invert(P, z0, r);
Q = point t + 4 of d;  Q' = invert(Q, z0, r);
R = point a of d;      R' = invert(R, z0, r);

draw subpath (7, 10) of c withcolor Blues 8 6;
draw d withcolor Oranges 8 5;
draw fullcircle zscaled (P' - Q') shifted 1/2[P', Q'] withcolor Oranges 8 5;

path t[]; 
t1 = P -- z1 -- R -- cycle; 
t2 = z1 -- Q -- R -- cycle;
t3 = invert_path(t1, z0, r, 1/16);
t4 = invert_path(t2, z0, r, 1/16);

forsuffixes $=t1, t3: fill $ withcolor Greens 8 2; endfor
forsuffixes $=t2, t4: fill $ withcolor Reds 8 2; endfor
forsuffixes $=t1, t2, t3, t4: draw $ withpen pencircle scaled 1/4; endfor

forsuffixes $ = Q', P: dotlabel.llft("\strut$" & str $ & "$", $); endfor
forsuffixes $ = R, R': dotlabel.ulft("$" & str $ & "$", $); endfor
dotlabel.rt("$Q$", Q);
dotlabel.lft("$O$", z0);
dotlabel.lrt("$P'$", P'); 
\end{mplibcode}$$

And even though you have to look quite hard, you can see from the curved
image of the triangle $PQR$ that what do get preserved
by inversion, are the angles between curves.   At least in the sense of
the angles between the tangents at the intersection of two curves, as here:
$$\begin{mplibcode}
path a, b;
a = (origin {dir 30} .. (160, -20) {dir 34}) shifted 80 left rotated -21;
b = (origin {dir 50} .. 170 right {dir 24}) shifted 60 left rotated  34;
numeric t, u;
(t, u) = a intersectiontimes b;

draw (-direction t of a -- direction t of a) shifted point t of a withcolor Blues 8 4;
draw (-direction u of b -- direction u of b) shifted point u of b withcolor Blues 8 4;

draw a; 
draw b;

label("$\alpha$", point t of a - 8 (unitvector(direction t of a) + unitvector(direction u of b)))
withcolor Blues 8 6;
\end{mplibcode}$$
The \MP\ interest here is drawing the tangents at the point of intersection.  \mpl{(t, u) = a intersectiontimes b} conveniently 
gives you the “time” along each path, then the tangent is 
\begin{code}
(- direction t of a -- direction t of a) shifted point t of a
\end{code}
except that the length of the direction vectors depends on the exact placement of the control points, 
so you if you prefer precision over brevity you could try:
\begin{code}
(unitvector(-direction t of a) -- unitvector(direction t of a)) 
    scaled 64 shifted point t of a
\end{code}

\end{document}