(*	twoDG

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

Copyright \[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 \[LessEqual] t \[LessEqual]
 \[Pi].";

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 \[LessEqual] a, b < 2\[Pi].";

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 \[LessEqual] \[Pi].";

direction::usage=
"direction[O,Q]  is angle  t = POQ,  where  P  lies right of
 O,  in radians, with  0 \[LessEqual] t < 2\[Pi].";
 
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
 \[Rule] 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\[Rule]Automatic  appended.  That lets you use  M
 to override  AspectRatio\[Rule]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`"];

(* 11.2 definitions *)
(* Heading :  Point  ADT *)

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

origin := Point[0, 0];

cis[\[Theta]_] := Point[Cos[\[Theta]], Sin[\[Theta]]];

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];

(* 11.2 definitions *)
(* Heading :  Vectors *)

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];
  
(* 11.2 definitions *)
(* Heading :  Measuring angles *)

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, \[Theta]},
      dx = c1[Q] - c1[O]; dy = c2[Q] - c2[O];
      If [dx == 0,
        If[dy == 0, noDirection,
          If[dy < 0, 3\[Pi]/2, \[Pi]/2]],
        \[Theta] = 
          ArcTan[dy/dx];             (* -\[Pi]/2 < \[Theta] < \[Pi]/2 *)
        If [dx < 0, \[Theta] + \[Pi],
          If[dy < 0, \[Theta] + 2\[Pi], \[Theta]]]]];

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

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

(* 11.3 definitions *)
(* Heading :  Text  primitive *)

$TextStyle = "Input";

(* 11.7 definitions *)
(* Heading :  The  lineq  ADT *)

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;

(* 11.7 definitions *)
(* Heading :  Constructing  lineq  objects *)

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]];

(* 11.7 definitions *)
(* Heading :  Join and intersection *)

Point[x1_, y1_]\[Vee]Point[x2_, y2_] ^:= 
    FullSimplify[lineq[y2 - y1, x1 - x2, x2*y1 - x1*y2]]; 
lineq[a1_, b1_, c1_] \[Intersection] 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]]]; 

(* 11.7 definitions *)
(* Heading :  Testing for concurrence *)

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]];

(* 11.7 definitions *)
(* Heading :  Additional operations *)

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];

(* Documented in exercise 1 : *)
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];

(* Documented in exercise 2 : *)

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

(* Documented in exercise 3 : *)
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]];

(* Documented in exercise 5 : *)
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) \[Or] (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]]) \[Or]    
          betweenQ[c2[P], c2[Q], c2[R]]);

(* 11.9 definitions *)
(* Heading :  Function  draw *)

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

(* 11.9 definitions *)
(* Heading :  Drawing segments and polygons *)

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]}]];

(* 11.9 definitions *)
(* Heading :  Clipping *)

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 \[Intersection] left;  R = g \[Intersection] right;
      B = g \[Intersection] bottom;  T = g \[Intersection] 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];

(* Documented in exercise 5 : *)
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\[Vee]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 \[Intersection] If[xO > xP, left, right];  
        Return[O Q]];
      If[verticalQ[g], Q = g \[Intersection] If[yO > yP, bottom, top];  
        Return[O Q]];
      xB = c1[g \[Intersection] bottom];  xT = c1[g \[Intersection] top];
      If[(xB < xleft && xT > xright) \[Or] (xT < xleft && xB > xright), 
        Q = g \[Intersection] If[xO > xP, left, right];  Return[O Q]];
      If[xleft <= xB <= xright && xleft <= xT <= xright, 
        Q = g \[Intersection] If[yO > yP, bottom, top];  Return[O Q]];
      If[xleft <= xB <= xright, 
        Q = g \[Intersection] If[yO > yP, bottom, If[xO > xP, left, right]];
        Return[O Q]];
      Q = g \[Intersection] If[yO > yP, If[xO > xP, left, right], top];
      Return[O Q]];
rayClipped[G_, O_, P_] := rayClipped[corners[G], O, P];

(* 11.9 definitions *)
(* Heading :  Arcs and circles *)

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

(* Documented in exercise 7 : *)
Circle[{x_, y_}, r_] \[Intersection] 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] \[Intersection] 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_] \[Intersection] Circle[{x2_, y2_}, r2_] ^:=
    If[{x1, y1} == {x2, y2}, If[r1 == r2, Circle[{x1, y1}, r1], {}], 
      Circle[{x1, y1}, r1] \[Intersection]
        lineq[2(x1 - x2), 2(y1 - y2), 
          x1^2 + x2^2 - y1^2 + y2^2 + r1^2 - r2^2]];

(* 11.9 definitions *)
(* Heading :  Angle marks *)

angleMark[P_, O_, Q_, r_, dr_, n_] := 
    Module[{\[Alpha], \[Beta]}, \[Alpha] = direction[O, P];  \[Beta] = 
        direction[O, Q];
      Table[arc[O, r + k*dr, \[Alpha], \[Beta]], {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}]];

(* 11.9 definitions *)
(* Heading :  Ticks *)

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/(2d)) T; 
      U = Point[-c2[T], c1[T]]; dU = (u/(2d))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];

(* 11.9 definitions *)
(* Heading :  Arrows *)

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

(* 11.10 definitions *)
(* Heading :  Animation *)

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_]\[Vee]Point[x2_, y2_] ^:= 
    FullSimplify[lineq[y2 - y1, x1 - x2, x2*y1 - x1*y2]]; 
Protect[Point];