(* This figure shows that the center of the unique circle through 
   three given points ("a", "b", and "c" in the figure) can be 
   constructed by finding the intersection of any two of the 
   three perpendicular bisectors to the sides of the triangle 
   defined by the three points. *)

CONST AngleSize = 15;

PROC DrawAngle(a, b, c) IS 
  IF 
    VAR 
      d ~ (0.1079, 0) REL (b, a), e ~ (0.5, 0.5) REL (d, f), 
      f ~ (0.1213, 0) REL (b, c) 
    IN 
      Geometry.Colinear(a, d, b) AND 
      Geometry.Colinear(b, f, c) AND 
      (d, b) CONG (b, f) AND 
      (a, b) PARA (e, f) AND 
      (d, e) PARA (b, c) AND 
      Geometry.Dist(b, d) = AngleSize -> 
        SAVE PS IN 
          PS.MoveTo(d); PS.LineTo(e); PS.LineTo(f); PS.Stroke() 
        END 
    END 
  FI 
END;

PROC Tri(a, b, c) IS 
  SAVE PS IN 
    PS.MoveTo(a); 
    PS.LineTo(b); 
    PS.LineTo(c); 
    PS.Close(); 
    SAVE PS IN PS.SetColor(Color.Red); PS.Fill() END; 
    PS.SetWidth(5); 
    PS.Stroke() 
  END 
END;

CONST SIZE = 4;

CONST SMALLSIZE = 2;

CONST WD = 1.5;

PRIVATE PROC Dot(center) IS 
  Circle.Draw(center, R2.Plus(center, (SIZE, 0))); PS.Fill() 
END;

PRIVATE PROC DrawCircle(origin, pt) IS 
  Dot(origin); Circle.Draw(origin, pt); PS.Stroke() 
END;

PRIVATE PROC DrawTriangle(a, b, c) IS 
  PS.SetWidth(WD); 
  PS.SetJointStyle(2); 
  PS.NewPath(); 
  PS.MoveTo(c); 
  PS.LineTo(a); 
  PS.LineTo(b); 
  PS.Close(); 
  PS.Save(); 
  PS.SetColor(Color.FromGrey(0.9)); 
  PS.Fill(); 
  PS.Restore(); 
  PS.Stroke() 
END;

PROC Doit(a, b, c) IS 
  VAR 
    len, msg1 = (20, 750), msg2 = (20, 700), 
    j ~ (0.2541, -0.1282) REL (a, c), e ~ (0.5, 0) REL (a, b), 
    f ~ (0.5, 0) REL (b, c), d ~ (0, 0), 
    h ~ (0.8677, -0.1424) REL (a, c) 
  IN 
    e = Geometry.Mid(a, b) AND 
    f = Geometry.Mid(b, c) AND 
    Angle.Right(d, e, b) AND 
    Angle.Right(d, f, c) AND 
    Geometry.Colinear(j, d, f) AND 
    Geometry.Colinear(h, d, e) -> 
      DrawTriangle(a, b, c); 
      DrawCircle(d, b); 
      PS.Type(msg1, Text.Cat("d = ", Unparse.Point(d))); 
      PS.Type(msg2, 
              Text.Cat("radius = ", 
                       Unparse.Value(Geometry.Dist(d, a)))) 
  END | SKIP 
END;

PROC Cmd0() IS 
  IF 
    VAR 
      a ~ (18.43907, -79.24909), b = (-125.14915, 47.137493), 
      c ~ (54.910183, 115.96457), d ~ (36.674625, 18.357738), 
      e ~ (-35.119484, 81.55103), f ~ (-14.627759, 27.94239) 
    IN 
      d = Geometry.Mid(a, c) AND 
      e = Geometry.Mid(b, c) AND 
      Angle.Right(b, e, f) AND 
      Angle.Right(a, d, f) -> 
        PS.SetWidth(3); 
        PS.SetJointStyle(PS.BevelJoints); 
        SAVE PS IN 
          Circle.Draw(f, a); PS.SetColor(Color.Red); PS.Stroke() 
        END; 
        Triangle.Draw(a, b, c); 
        PS.Stroke(); 
        PS.SetWidth(2); 
        PS.MoveTo(e); 
        PS.LineTo(f); 
        PS.LineTo(d); 
        PS.Stroke(); 
        DrawAngle(b, e, f); 
        DrawAngle(a, d, f) 
    END 
  FI 
END;