(*	twoDG

Mathematica add-on package
for enhancement of two-dimensional graphics,  
described in chapter 11 of manuscript
Mathematics and Mathematica,

Copyright © 2001  James T. Smith *)

BeginPackage["twoDG`","Graphics`Colors`"];

Off[General::spell];
Off[General::spell1];
Unprotect[Point];
Unprotect[Circle];

angle::usage=
"angle[Q,O,R]  is the angle  t = QOR  in radians, measured as
 in elementary geometry, with  0 ≤ t ≤
 π.";

angleMark::usage=
"angleMark[P,O,Q,r,dr,n]  marks angle  POQ  with  n  circular
 arcs from ray  OP  counterclockwise to ray  OQ,  starting at
 radius  r  and stepping outward at intervals of width  dr.
 For  n = 1  you can omit  dr  and  n.";
 
animate::usage=
"animate  works only if placed on the line immediately after an
 expression whose value is a list of  Graphics  objects.  It
 displays them in succession cyclically, each for  0.1  second,
 and closes the animated cell group when you stop the
 animation.";
 
arc::usage=
"arc[Z,r,a,b]  is the circular arc with center  Z  and radius
 r  from direction  a  counterclockwise to  direction  b,
 where  0 ≤ a, b < 2π.";

arrowhead::usage=
"arrowhead[V,d,s]  is an arrowhead at  V  pointed in direction
 d,  with length  s."; 

badEquation::usage=
"twoDG`  returns the undefined symbol  badEquation  in certain
 cases when you ask if nonexistent objects are equal.";

badIntersection::usage=
"twoDG`  returns the undefined symbol  badIntersection  if you
 ask for the intersection of equal or parallel lines.";

betweenQ::usage=
"betweenQ[r,s,t]  is  True  just in case  s  lies between  r
 and  t  (inclusive).";

between1Q::usage=
"between1Q[P,Q,R]  is  True  just in case  Q  lies between  P
 and  R  (inclusive).";
 
between2Q::usage=
"between2Q[P,Q,R]  is  True  just in case  Q  lies between  P
 and  R  (inclusive).";
 
c1::usage=
"c1[P]  is the 1st coordinate of point  P.
 c1[g]  is the 1st coefficient of  lineq  object  g.";

c2::usage=
"c2[P]  is the 2nd coordinate of point  P.
 c2[g]  is the 2nd coefficient of  lineq  object  g.";

c3::usage=
"c3[g]  is the 3rd coefficient of  lineq  object  g.";

Circle::usage=
"Circle[{x, y},r]  is a two-dimensional graphics primitive that
 represents a circle of radius  r  centered at the point  x,y.
 Circle[{x, y},{rx, ry}] yields an ellipse with semi-axes  rx
 and ry. Circle[{x,y},r,{theta1,theta2}] represents a circular
 arc. 
 Circle[P,r]  is the circle with center  P  and radius  r.";

cis::usage=
"cis[t]  is the point on the unit circle at angle  t  radians.";

collinearQ::usage=
"collinearQ[P,Q,R]  is  True  just in case  P, Q, R  fall on a
 single line.";

concurrentQ::usage=
"concurrentQ[g,h,j]  is  True  just in case  g, h, j  are
 concurrent or mutually parallel.";

concurrentQT::usage=
"concurrentQT[g,h,j]  is  True  just in case  g, h, j  are
 concurrent or mutually parallel.  It also prints the
 determinant of the matrix of coefficients of g, h, j.";
 
corners::usage=
"corners[G]  is the list consisting of the lower left and upper
 right corners of the  PlotRange  of  Graphics  object  G.";
 
counterclockwiseQ::usage=
"counterclockwiseQ[P,Q,R]  is  True  just in case  P, Q, R
 are noncollinear and in going from  P  to  Q  to  R  to  P
 you proceed counterclockwise.";

directedAngle::usage=
"directedAngle[Q,O,R] is the directed angle  t = QOR  in
 radians with counterclockwise sense positive and
 -Pi < t ≤ π.";

direction::usage=
"direction[O,Q]  is angle  t = POQ,  where  P  lies right of
 O,  in radians, with  0 ≤ t < 2π.";
 
distance::usage=
"distance[P,Q]  is the distance between points  P,Q.
 distance[P,g]  is the distance from point  P  to  line  g.";

draw::usage=
"draw[L]  sends  Graphics[L]  to  Show,  with  AspectRatio
 -> Automatic,  for rendering.  It returns the  Graphics
 object returned by  Show.
 draw[L,M]  does the same, only it sends option list  M  with
 AspectRatio->Automatic  appended.  That lets you use  M
 to override  AspectRatio->Automatic.";
 
filledPolygon::usage=
"filledPolygon[L]  is the closed filled polygon with vertex
 list  L  of  Point  objects.  The last vertex is considered
 joined to the first.";

horizontal::usage=
"horizontal[y0]  is the line with equation  y = y0.
 horizontal[P]   is the horizontal line through  P.";

horizontalQ::usage=
"horizontalQ[g]  is  True  just in case  g  is horizontal.";

lineClipped::usage=
"lineClipped[{LL,UR},g_]  is the segment formed by clipping  g
 to the rectangle with specified lower left and upper right
 corners.
 lineClipped[G,g]  is the segment formed by clipping  g  to the
 PlotRange  of  Graphics  object  G.";

lineq::usage=
"lineq[a,b,c]  is the line with equation  ax+by+c = 0.  twoDG`
 includes many overload definitions for operators on  lineq
 objects.";

midpoint::usage=
"midpoint[P,Q]  is the midpoint between P,Q.";

noCircle::usage=
"twoDG`  substitutes the undefined symbol  noCircle  when you
 refer to a circle with negative radius.";

norm::usage=
"norm[P]  is the distance between  P  and the origin.";

noDirection::usage=
"twoDG`  returns the undefined symbol  noDirection  if you ask
 for a nonexistent direction.";

noLine::usage=
"twoDG`  substitutes the undefined symbol  noLine  when you
 refer to  lineq[0,0,c].";
 
noSlope::usage=
"twoDG`  returns the undefined symbol  noSlope  if you ask for
 a nonexistent slope.";

onQ::usage=
"onQ[P,g]  is  True  just in case  P  lies on  g.";

origin::usage=
"origin  is the point with coordinates zero.";

parallelQ::usage=
"parallelQ[g,h]  is  True  just in case  g  and  h  are
 parallel or equal.";

perpendicularQ::usage=
"perpendicularQ[g,h]  is  True  just in case  g  and  h  are
 perpendicular.";

Point::usage=
"Point[coords] is a graphics primitive that represents a point.
 Point[x,y]  is the point with coordinates  x,y.  twoDG`
 includes many overload definitions for operators on  Point
 objects.";
 
pointParallel::usage=
"pointParallel[P,g]  is the line through  P  parallel to  g.";

pointPerpendicular::usage=
"pointPerpendicular[P,g]  is the line through  P  perpendicular
 to  g.";
 
pointSlope::usage=
"pointSlope[P,m]  is the line through  P  with slope  m.";

polygon::usage=
"polygon[L]  is the polygon (not necessarily closed) with
 vertex list  L  of  Point  objects.";
 
rayClipped::usage=
"rayClipped[{LL,UR},O,P]  is the segment formed by clipping the
 ray  OP  that starts from  O  and passes through  P  to the
 rectangle with specified lower left and upper right corners.
 rayClipped[G,O,P]  is the segment formed by clipping  OP  to
 the PlotRange  of  Graphics  object  G.";
 
rightAngle::usage=
"rightAngle[P,O,Q,s]  marks right angle  POQ  with a square of
 edge length  s.";

slope::usage=
"slope[P,Q]  is the slope of segment PQ.
 slope[g]    is the slope of line  g.";
 
tick::usage=
"tick[P,Q,u,t,n]  marks segment  PQ  with  n  ticks of length
 u,  centered with an interval of width  t.
 For  n = 1  you may omit  t  and  n.";

throughQ::usage=
"throughQ[g,P]  is  True  just in case  g  passes through  P.";

unitCircle::usage=
"unitCircle  is the circle with radius 1 centered at the
 origin.";

vertical::usage=
"vertical[x0]  is the line with equation  x = x0.
 vertical[P]   is the vertical line through  P.";

verticalQ::usage=
"verticalQ[g]  is  True  just in case  g  is vertical.";

Begin["`Private`"];

Point[x_, y_] := 
        Point[{x, y}]; 

origin := Point[0, 0];

cis[θ_] := 
        Point[Cos[θ], 
          Sin[θ]];

c1[Point[{x_, y_}]] ^:= x;
c2[Point[{x_, y_}]] ^:= 
        y;

Point[P_] == Point[Q_] ^:= 
        FullSimplify[P == Q];
Point[P_] ≠  Point[Q_] ^:= 
        Not[Point[P] == 
            Point[Q]];

distance[Point[{x1_, y1_}], 
          Point[{x2_, 
              y2_}]] ^:= 
Sqrt[(x1 - x2)^2 + (y1 - 
                  y2)^2];

Point[P_] + Point[Q_] ^:= Point[P + Q];
t_*Point[P_] ^:= Point[t*P];

\(\(midpoint[P_, 
          Q_] := \((1/2)\) \((P + 
              Q)\);\)

 (*\ \(Documented\ in\ exercise\ 5\)\(:\)\ *) \)
norm[P_] := distance[origin, P];

slope[Point[{x1_, y1_}], Point[{x2_, y2_}]] ^:= 
        Module[{dx, dy}, 
          dx = x2 - x1;   dy = y2 - y1; 
If[
            FullSimplify[dx == 0], noSlope, 
            dy/dx]];

direction[O_, Q_] := 
        Module[{dx, dy, θ}, 
dx = c1[Q] - c1[O]; 
          dy = c2[Q] - c2[O]; 
If [
            dx == 0, 
If[dy == 0, 
              noDirection, 
If[dy < 0, 
                3  π/2, π/2]], 
θ = 
              ArcTan[dy/dx];                
If [
              dx < 0, θ + π, 
If[
                dy < 0, θ + 
                  2  π, θ]]]];

directedAngle[Q_, O_, R_] := 
        Module[{θ}, 
θ = 
            Mod[direction[O, R] - direction[O, Q], 
              2  π]; 
If[θ > π, θ = θ - 2  π, θ]];

angle[Q_, O_, R_] := Abs[directedAngle[Q, O, R]];

$TextStyle = "Input";

lineq[0, 0, c_] := noLine;
 
c1[lineq[a_, b_, c_]] ^:= a;
c2[lineq[a_, b_, c_]] ^:= b;
\(\(c3[lineq[a_, b_, c_]] ^:= c;\)
\)
lineq[a0_, b0_, c0_] == 
          lineq[a1_, b1_, c1_] ^:= 
FullSimplify[
          a0*b1 == a1*b0  && 
b0*c1 == b1*c0 && 
            c0*a1 == c1*a0];
lineq[a0_, b0_, c0_] ≠ lineq[a1_, b1_, c1_] ^:= 
        Not[lineq[a0, b0, c0] == lineq[a1, b1, c1]];
noLine == noLine ^= badEquation;
noLine ≠ noLine ^= badEquation;

pointSlope[P_, m_] := 
        lineq[m, -1, c2[P] - m*c1[P]];
 
horizontal[y_] :=  lineq[0, 1, -y];
horizontal[Point[x_, y_]] ^:=  horizontal[y];
vertical[x_] := lineq[1, 0, -x];
vertical[Point[x_, y_]] ^:= 
        vertical[x];

pointParallel[P_, g_] := 
        lineq[c1[g], c2[g], -c1[g] c1[P] - c2[g] c2[P]];
pointPerpendicular[P_, g_] := 
        lineq[-c2[g], c1[g], c2[g] c1[P] - c1[g] c2[P]];

Point[x1_, y1_]⋁Point[x2_, y2_] ^:= 
        FullSimplify[
          lineq[y2 - y1, x1 - x2, x2*y1 - x1*y2]];

lineq[a1_, b1_, c1_] ⋂ lineq[a2_, b2_, c2_] ^:= 
        Module[{d}, 
          d = a1*b2 - a2*b1; 
If[
            FullSimplify[d == 0], 
            badIntersection,  (1/d)  Point[b1*c2 - b2*c1, 
                c1*a2 - c2*a1]]];

\(\(concurrentQ[
          lineq[a1_, b1_, c1_], 
		\ \ lineq[a2_, b2_, 
            c2_], 
		\ \ lineq[a3_, b3_, c3_]] ^:= 
        FullSimplify[
          Det[{{a1, b1, c1}, 
\ \ \ \ \ {a2, b2, 
                  c2}, 
\ \ \ \ \ {a3, b3, c3}}] == 
            0];\)\ 
\)
\(\(concurrentQT[
          lineq[a1_, b1_, c1_], 
		\ \ \ \ \ lineq[a2_, 
            b2_, c2_], 
		\ \ \ \ \ lineq[a3_, b3_, 
            c3_]] ^:= Module[{d}, 
          d = Det[{{a1, b1, c1}, 
\ \ \ \ \ \ \ {a2, b2, 
                  c2}, 
\ \ \ \ \ \ \ {a3, b3, 
                  c3}}]; 
Print["Det = ", 
            d]; 
FullSimplify[d == 0]];\)\)

distance[Point[x_, y_], lineq[a_, b_, c_]] ^:= 
        Abs[a*x + b*y + c]/Sqrt[a^2 + b^2]; 

slope[lineq[a_, b_, c_]] ^:= 
If[
          FullSimplify[b == 0], 
          noSlope, -a/
            b];

  
horizontalQ[g_] := c1[g] == 0;
verticalQ[g_] := c2[g] == 0;
onQ[P_, g_] := 
        FullSimplify[
          c1[P] c1[g] + c2[P] c2[g] + c3[g] == 0];
throughQ[g_, P_] := 
        onQ[P, g];

  

between1Q[P_, Q_, R_] := 
        FullSimplify[
          distance[P, Q] + distance[Q, R] == 
            distance[P, 
              R]];

  
perpendicularQ[g_, h_] := 
        FullSimplify[c1[g] c1[h] + c2[g] c2[h] == 0];
parallelQ[g_, h_] := 
        FullSimplify[
          c1[g] c2[h] == 
            c1[h] c2[
                g]];

  
\(\(collinearQ[
          Point[x1_, 
            y1_], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Point[x2_, 
            y2_], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Point[x3_, y3_]] := 
        FullSimplify[
          Det[{{x1, y1, 1}, 
\ \ \ \ \ {x2, y2, 
                  1}, 
\ \ \ \ \ {x3, y3, 1}}] == 
            0];\)

 (*\ \(Documented\ in\ exercise\ 6\)\(:\)\ *) \)
\(\(counterclockwiseQ[
          Point[x1_, 
            y1_], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Point[x2_, 
            y2_], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Point[x3_, y3_]] := 
        FullSimplify[
          Det[{{x1, y1, 1}, 
\ \ \ \ \ {x2, y2, 
                  1}, 
\ \ \ \ \ {x3, y3, 1}}] > 
            0];\)

 (*\ \(Documented\ in\ exercise\ 7\)\(:\)\ *) \)
\(\(betweenQ[s_, t_, u_] := 
        FullSimplify[\((s ≤ t ≤ 
                u)\) ∨ \((u ≤ t ≤ 
                s)\)];\)

 (*\ \(Documented\ in\ exercise\ 8\)\(:\)\ *) \)
between2Q[P_, Q_, 
          R_] := 
collinearQ[P, Q, 
            R] && ((FullSimplify[c2[P] == c2[Q]] && 
                  betweenQ[c1[P], c1[Q], c1[R]]) ∨     betweenQ[
                c2[P], c2[Q], c2[R]]);

draw[L_] := 
        Show[Graphics[L], AspectRatio -> Automatic];
draw[L_, M_] := 
        Show[Graphics[L], 
          Append[{M}, AspectRatio -> Automatic]];

Point[{x0_, y0_}] Point[{x1_, y1_}] ^:= 
        Line[{{x0, y0}, {x1, y1}}];
polygon[L_] := 
        Line[Table[{c1[L[[k]]], c2[L[[k]]]}, {k, 1, 
              Length[L]}]];
filledPolygon[L_] := 
        Polygon[Table[{c1[L[[k]]], c2[L[[k]]]}, {k, 1, 
              Length[L]}]];

corners[G_] := 
        Module[{xMin, xMax, yMin, yMax}, {{xMin, xMax}, {yMin, yMax}} = 
            PlotRange /. AbsoluteOptions[G]; 
{Point[xMin, 
              yMin], Point[xMax, yMax]}];
lineClipped[{LL_, UR_}, g_] := 
        Module[{xleft, xright, ybottom, ytop, left, right, bottom, top, L, R, 
            B, T, xB, xT}, 
          If[g == noLine, Return[{}]]; 
xleft = 
            c1[LL];  xright = c1[UR]; 
ybottom = c2[LL];  
          ytop = c2[UR]; 
left = vertical[LL];  
          right = vertical[UR]; 
bottom = horizontal[LL];  
          top = horizontal[UR]; 
L = 
            g ⋂ left;   R = 
            g ⋂ right; 
B = 
            g ⋂ bottom;   T = 
            g ⋂ top; 
xB = c1[B];   xT = 
            c1[T]; 
If[horizontalQ[g], 
            Return[If[ybottom ≤ c2[L] ≤ ytop, 
                L R, {}]]]; 
If[verticalQ[g], 
            Return[If[xleft ≤ xB ≤ xright, 
                B T, {}]]]; 
If[xT < xleft, 
            Return[If[xB < xleft, {}, 
                If[xB ≤ xright, L B, 
                  L R]]]]; 
If[xT ≤ xright, 
            Return[If[xB < xleft, T L, 
                If[xB ≤ xright, T B, 
                  T R]]]]; 
If[xB < xleft, L R, 
            If[xB ≤ xright, B R, {}]]];
lineClipped[G_, g_] := 
        lineClipped[corners[G], 
          g];

  
rayClipped[{LL_, UR_}, O_, P_] := 
        Module[{g, Q, xO, yO, xP, yP, left, right, bottom, top, xB, xT, 
            xleft, xright}, 
If[O == P, 
            Return[O P]];  
          g = O⋁P; 
{xO, yO} = {c1[O], 
              c2[O]}; 
{xP, yP} = {c1[P], 
              c2[P]}; 
left = vertical[LL];  
          right = vertical[UR]; 
bottom = horizontal[LL];  
          top = horizontal[UR]; 
xleft = c1[LL];  
          xright = c1[UR]; 
If[horizontalQ[g], 
            Q = g ⋂ If[xO > xP, left, right];   Return[
              O Q]]; 
If[verticalQ[g], 
            Q = g ⋂ If[yO > yP, bottom, top];   Return[
              O Q]]; 
xB = 
            c1[g ⋂ bottom];   xT = 
            c1[g ⋂ 
                top]; 
If[(xB < xleft && 
                  xT > xright) ∨ (xT < xleft && xB > xright), 
            Q = g ⋂ If[xO > xP, left, right];   Return[
              O Q]]; 
If[
            xleft ≤ xB ≤ xright && 
              xleft ≤ xT ≤ xright, 
            Q = g ⋂ If[yO > yP, bottom, top];   Return[
              O Q]]; 
If[
            xleft ≤ xB ≤ xright, 
            Q = g ⋂ 
                If[yO > yP, bottom, 
                  If[xO > xP, left, right]]; 
Return[
              O Q]]; 
Q = 
            g ⋂ 
              If[yO > yP, If[xO > xP, left, right], 
                top]; 
Return[O Q]];
rayClipped[G_, O_, P_] := rayClipped[corners[G], O, P];

Circle[Point[x_, y_], r_] ^:= 
If[r < 0, noCircle, 
          Circle[{x, y}, r]];
unitCircle = Circle[origin, 1];
arc[Z_, r_, α_, β_] := 
        If[r < 0, noCircle, 
          Circle[{c1[Z], c2[Z]}, 
            r, {α, 
              If[α ≤ β, β, β + 
                  2  π]}]];

  
Circle[{x_, y_}, r_] ⋂ lineq[a_, b_, c_] ^:= 
        Module[{Z, g, d, P, V, s, Q}, 
Z = 
            Point[x, y];   g = lineq[a, b, c]; 
d = 
            distance[Z, g]; 
If[d > r, 
            Return[{}]]; 
P = 
            pointPerpendicular[Z, g] ⋂ g; 
If[
            FullSimplify[d == r], Return[{P}]]; 
V = 
            Point[b, -a];   s = 
            norm[V]; 
Q = (Sqrt[r^2 - d^2]/s) 
              V; 
{P + Q, P - Q}];
Circle[{x1_, y1_}, r1_] ⋂ 
          Circle[{x2_, y2_}, 
            r2_] ^:= 
If[{x1, y1} == {x2, y2}, 
          If[r1 == r2, Circle[{x1, y1}, r1], {}], 
          Circle[{x1, y1}, r1] ⋂ 
lineq[
              2 (x1 - x2), 2 (y1 - y2), 
              x1^2 + x2^2 - y1^2 + y2^2 + r1^2 - r2^2]];

angleMark[P_, O_, Q_, r_, dr_, n_] := 
        Module[{α, β}, α = direction[O, P];   β = 
            direction[O, Q]; 
Table[
            arc[O, r + k*dr, α, β], {k, 0, n - 1}]];
angleMark[P_, O_, Q_, r_] := angleMark[P, O, Q, r, 0, 1];
rightAngle[P_, O_, Q_, s_] := 
        Module[{t, u, T, U, V}, 
          t = s/distance[O, P]; 
u = 
            s/distance[O, Q]; 
T = 
            t (P - O); 
U = 
            u (Q - O); 
polygon[{O + T, O + T + U, 
              O + U}]];

tick[P_, Q_, u_, t_, n_] := 
        Module[{M, d, T, U, dT, dU}, 
          M = midpoint[P, Q];               d = 
            distance[P, Q]; 
T = 
            Q - P;                                      dT = (t/(2  d)) T;  U = Point[-c2[T], c1[T]];  
          dU = (u/(2  d)) 
              U; 
Table[(M + k dT + dU) (M + k dT - 
                  dU), {k, -(n - 1), n - 1, 2}]];
tick[P_, Q_, u_] := tick[P, Q, u, 0, 1];

arrowhead[V_, θ_, d_] := 
        Module[{P, Q}, 
P = 
            V + d*cis[θ + 5  π/4];   Q = 
            V + d*cis[θ + 3  π/4]; 
{V P, 
            V Q}];

\(\(animate := {SelectionMove[SelectedNotebook[], All, 
            GeneratedCell]; 
FrontEndTokenExecute["SelectionCloseAllGroups"]; 
SelectionAnimate[SelectedNotebook[], 
            AnimationDisplayTime -> 0.1]};\)\)

Protect["twoDG`Private`*"];
End[];
Protect["twoDG`*"];
On[General::spell1];
On[General::spell];
Protect[Circle];
Protect[Point];
EndPackage[];

(* 
      Appended to work around some bizarre glitch: *) 
Unprotect[Point];
Point[x1_, y1_]⋁Point[x2_, y2_] ^:= 
        FullSimplify[
          lineq[y2 - y1, x1 - x2, 
            x2*y1 - x1*y2]]; 
Protect[Point];