## 20354 Reputation

16 years, 11 days

## Solution of Gardens puzzle...

In this post we present the solution with Maple to the logical problem of "Gardens Puzzle"

http://www.mathsisfun.com/puzzles/gardens-solution.html

The Puzzle:

Five friends have their gardens next to one another, where they grow three kinds of crops: fruits (apple, pear, nut, cherry), vegetables (carrot, parsley, gourd, onion) and flowers (aster, rose, tulip, lily).

1. They grow 12 different varieties.
2. Everybody grows exactly 4 different varieties
3. Each variety is at least in one garden.
4. Only one variety is in 4 gardens.
5. Only in one garden are all 3 kinds of crops.
6. Only in one garden are all 4 varieties of one kind of crops.
7. Pear is only in the two border gardens.
8. Paul's garden is in the middle with no lily.
9. Aster grower doesn't grow vegetables.
10. Rose growers don't grow parsley.
11. Nuts grower has also gourd and parsley.
12. In the first garden are apples and cherries.
13. Only in two gardens are cherries.
14. Sam has onions and cherries.
15. Luke grows exactly two kinds of fruit.
16. Tulip is only in two gardens.
17. Apple is in a single garden.
18. Only in one garden next to Zick's is parsley.
19. Sam's garden is not on the border.
20. Hank grows neither vegetables nor asters.
21. Paul has exactly three kinds of vegetable.

Who has which garden and what is grown where?

About methods of solution. At first I just wanted to generate all variations and using conditions 1 .. 21 to find all solutions. But even if we use the condition that everybody grows exactly 4 different varieties then the total number variants equals  5!^2*binomial(12,4)^5=427945522455000000

So from the very beginning using some of the conditions 1 .. 21 we maximally reduce the number of possible variants. For example from the conditions 11, 18 and 6 implies that only in one garden are all 4 varieties of flowers. Next we pass through these variants and using conditions 1 .. 21 and finally come to a unique solution:

restart;

Fruits:={apple, pear, nut, cherry}:

Vegetables:={carrot, parsley, gourd, onion}:

Flowers:={aster, rose, tulip, lily}:

Set1:=Flowers:

Garden1:={Set1}:

Set2:=Fruits union Vegetables union Flowers minus {nut,gourd,parsley} minus {apple,cherry,rose}:

Garden2:={seq({nut,gourd,parsley} union {Set2[i]}, i=1..nops(Set2))}:

Set3:=Vegetables union Flowers minus {parsley}:

Garden3:={seq({apple,cherry,pear} union {Set3[i]}, i=1..nops(Set3))}:

Set4:=combinat[choose](Fruits union Vegetables union Flowers minus {onion,cherry} minus {apple,parsley,pear,nut}, 2):

Garden4:={seq({onion,cherry} union Set4[i], i=1..nops(Set4))}:

Set5:=Fruits union Vegetables union Flowers minus {apple, cherry,nut,parsley}:

Garden5:=combinat[choose](Set5, 4):

S:=[]:

for s1 in Garden1 do

for s2 in Garden2 do

for s3 in Garden3 do

for s4 in Garden4 do

for s5 in Garden5 do

s:=[s1,s2,s3,s4,s5]: s_4:=combinat[choose](s,4): m:=0: n:=0: k:=0: p:=0: q:=0:

for i in s do

if `intersect`(i,Fruits)<>{} and `intersect`(i,Vegetables)<>{} and `intersect`(i,Flowers)<>{} then m:=m+1: fi:

if pear in i then n:=n+1: fi:

if tulip in i then k:=k+1: fi:

if aster in i and `intersect`(i,Vegetables)<>{} then p:=p+1: fi:

if i=Fruits or i=Vegetables or i=Flowers then q:=q+1: fi:

od:

if nops(`union`(op(s)))=12 and nops(`union`(seq(`intersect`(op(s_4[j])), j=1..nops(s_4))))=1 and m=1 and n=2 and k=2 and p=0 and q=1 then S:=[op(S),[s1,s2,s3,s4,s5]]: fi:

od: od: od: od: od:

L1:=[seq([[3,Paul],[2,Sam],seq([combinat[permute]([1,4,5])[i,j],[Luke,Zick,Hank][j]],j=1..3)],i=1..6)]:

L2:=[seq([[3,Paul],[4,Sam],seq([combinat[permute]([1,2,5])[i,j],[Luke,Zick,Hank][j]],j=1..3)],i=1..6)]:

L0:=[op(L1),op(L2)]:

L:=[seq(op(combinat[permute](L0[i])),i=1..nops(L0))]:

Sol:=[]:

for l in L do

for s in S do

sol:=[seq([op(l[i]),s[i]], i=1..5)]:

if abs(gad1[1]-gad2[1])=1 and convert([seq(((pear in sol[i][3] implies (sol[i][1]=1 or sol[i][1]=5)) and (sol[i][2]=Paul implies (not lily in sol[i][3])) and (sol[i][1]=1 implies (apple in sol[i][3] and cherry in sol[i][3])) and (sol[i][2]=Sam implies (onion in sol[i][3] and cherry in sol[i][3])) and (sol[i][2]=Luke implies nops(`intersect`(sol[i][3],Fruits))=2) and (sol[i][2]=Hank implies (`intersect`(sol[i][3],Vegetables)={} and not aster in sol[i][3])) and (sol[i][2]=Paul implies nops(`intersect`(sol[i][3],Vegetables))=3)), i=1..5)], `and`) then Sol:=[op(Sol), sol]: fi:

od: od:

for i in Sol do

Matrix(sort(i,(x,y)->x[1]<y[1]));

od;

## Another proof of six points theorem...

In this post we present another compact proof of this remarkable theorem without using  geometry package.
The proof uses a procedure called  Cc , which for three points returns a list of the coordinates of the center and the radius of the circumscribed circle.

restart;

Cc:=proc(A,B,C)

local x1, y1, x2, y2, x3, y3, x, y;

x1,y1:=op(A);  x2,y2:=op(B);  x3,y3:=op(C);

solve({(x2-x1)*(x-(x1+x2)/2)+(y2-y1)*(y-(y1+y2)/2)=0, (x2-x3)*(x-(x2+x3)/2)+(y2-y3)*(y-(y2+y3)/2)=0},{x,y});

assign(%);

[simplify([x,y]), simplify(sqrt((x-x1)^2+(y-y1)^2))];

end proc:

Proof for arbitrary triangle:

A, B, C:=[x1,y1], [x2,y2], [x3,y3]:

A1, B1, C1, M:=(B+C)/2, (A+C)/2, (A+B)/2, (A+B+C)/3:

P1:=Cc(A,M,B1)[1]: P2:=Cc(B1,M,C)[1]: P3:=Cc(C,M,A1)[1]:

P4:=Cc(A1,M,B)[1]: P5:=Cc(B,M,C1)[1]: P6:=Cc(C1,M,A)[1]:

Cc1:=Cc(P1,P2,P3):  Cc2:=Cc(P4,P5,P6):

is(Cc1=Cc2);

true

## Congratulations on Valentine's Day...

Code of the animation:

restart;

N := 192:

A := seq(plot([[.85*sin(t)^3-2+1.25*i/N, .85*(13*cos(t)*(1/15)-(1/3)*cos(2*t)-2*cos(3*t)*(1/15)-(1/15)*cos(4*t)), t = 0 .. Pi*i/N], [-.85*sin(t)^3-2+1.25*i/N, .85*(13*cos(t)*(1/15)-(1/3)*cos(2*t)-2*cos(3*t)*(1/15)-(1/15)*cos(4*t)), t = 0 .. Pi*i/N], [sin(t)^3+2-1.25*i/N, 13*cos(t)*(1/15)-(1/3)*cos(2*t)-2*cos(3*t)*(1/15)-(1/15)*cos(4*t), t = 0 .. Pi*i/N], [-sin(t)^3+2-1.25*i/N, 13*cos(t)*(1/15)-(1/3)*cos(2*t)-2*cos(3*t)*(1/15)-(1/15)*cos(4*t), t = 0 .. Pi*i/N]], color = red, thickness = 5, view = [-3 .. 3, -1.2 .. .9]), i = 1 .. N):

plots[display](A, insequence = true, scaling = constrained, axes = none);

## Testing of two plane sets for isometry...

An isometry of the Euclidean plane is a distance-preserving transformation of the plane.  There are four types: translations,  rotations,  reflections,  and  glide reflections. See http://en.wikipedia.org/wiki/Euclidean_plane_isometry#Classification_of_Euclidean_plane_isometries

Procedure  AreIsometric  checks two plane sets  Set1  and  Set2  and  if there is an isometry of the plane mapping the first set to the second one, then the procedure returns true and false otherwise. Global variable  T  saves the type of isometry and all its parameters. For example, it returns  the rotation center and rotation angle, etc.

Each of the sets  Set1  and  Set2   is the set (or list) consisting of the following objects in any combinations:

1) The points that are defined as a list of coordinates  [х, y] .

2) Segments, which are defined as a set (or list)  of two points  {[x1, y1], [x2, y2]}  or  [[x1, y1], [x2, y2]]  .

3) Curves, which should  be defined as a list of points   [[x1, y1], [x2, y2], ..., [xn, yn]].

Of course, if  n = 2, then the curve is identical to the segment.

Code of the procedure:

AreIsometric:=proc(Set1::{set,list}, Set2::{set,list})

local n1, n2, n3, n4,s1, S, s, l1, l2, S11, f, x0, y0, phi, Sol, x, y, M1, M2, A1, A2, A3, A4, B1, B2, B3, B4, line1, line2, line3, line4, u, v, Sign, g, M, Line1, Line2, Line3, A, B, C, h, AB, CD, Eq, Eq1, T1, T2, i, S1, S2, T11;

global T;

uses combinat;

S1:={};  S2:={};  T1:={}; T2:={};

for i in Set1 do

if i[1]::realcons  then S1:={op(S1),i} else

S1:={op(i), op(S1)};  T1:={op(T1), seq({i[k],i[k+1]}, k=1..nops(i)-1)} fi;

od;

for i in Set2 do

if i[1]::realcons  then S2:={op(S2),i} else

S2:={op(i), op(S2)};  T2:={op(T2), seq({i[k],i[k+1]}, k=1..nops(i)-1)} fi;

od;

n1:=nops(S1);  n2:=nops(S2);  n3:=nops(T1); n4:=nops(T2);

if is(S1=S2) and is(T1=T2) then T:=identity;  return true fi;

if n1<>n2 or n3<>n4 then return false fi;

if n1=1 then T:=[translation, <S2[1,1]-S1[1,1], S2[1,2]-S1[1,2]>];  return true fi;

f:=(x,y,phi)->[(x-x0)*cos(phi)-(y-y0)*sin(phi)+x0, (x-x0)*sin(phi)+(y-y0)*cos(phi)+y0];  g:=(x,y)->[(B^2*x-A^2*x-2*A*B*y-2*A*C)/(A^2+B^2), (A^2*y-B^2*y-2*A*B*x-2*B*C)/(A^2+B^2)];

_Envsignum0 := 1;

s1:=[S1[1], S1[2]];  S:=select(s->is((s1[2,1]-s1[1,1])^2+(s1[2,2]-s1[1,2])^2=(s[2,1]-s[1,1])^2+(s[2,2]-s[1,2])^2),permute(S2, 2));

for s in S do

# Checking for translation

l1:=s[1]-s1[1]; l2:=s[2]-s1[2];

if is(l1=l2) then S11:=map(x->x+l1, S1);

if n3<>0 then T11:={seq(map(x->x+l1, T1[i]), i=1..nops(T1))}; fi;

if n3=0 then  if is(S11=S2) then T:=[translation, convert(l1, Vector)]; return true fi;  else

if is(S11=S2) and is(T11=T2) then T:=[translation, convert(l1, Vector)]; return true fi; fi;

fi;

# Checking for rotation

x0:='x0'; y0:='y0'; phi:='phi'; u:='u'; v:='v'; Sign:='Sign';

if  is(s1[1]-s[1]<>s1[2]-s[2]) then

M1:=[(s1[1,1]+s[1,1])/2, (s1[1,2]+s[1,2])/2]; M2:=[(s1[2,1]+s[2,1])/2, (s1[2,2]+s[2,2])/2]; A1:=s1[1,1]-s[1,1]; B1:=s1[1,2]-s[1,2]; A2:=s1[2,1]-s[2,1]; B2:=s1[2,2]-s[2,2];    line1:=A1*(x-M1[1])+B1*(y-M1[2])=0; line2:=A2*(x-M2[1])+B2*(y-M2[2])=0;

if is(A1*B2-A2*B1<>0) then Sol:=solve({line1, line2}); x0:=simplify(rhs(Sol[1]));   y0:=simplify(rhs(Sol[2])); u:=[s1[1,1]-x0,s1[1,2]-y0]; v:=[s[1,1]-x0,s[1,2]-y0];    else

if is(s[2]-s1[1]=s[1]-s1[2])  then   x0:=(s1[1,1]+s[1,1])/2;  y0:=(s1[1,2]+s[1,2])/2;

if is([x0,y0]<>s1[1]) then  u:=[s1[1,1]-x0,s1[1,2]-y0]; v:=[s[1,1]-x0,s[1,2]-y0]; else

u:=[s1[2,1]-x0,s1[2,2]-y0]; v:=[s[2,1]-x0,s[2,2]-y0]; fi;

else  A3:=s1[2,1]-s1[1,1];  B3:=s1[2,2]-s1[1,2]; A4:=s[2,1]-s[1,1];  B4:=s[2,2]-s[1,2];  line3:=B3*(x-s1[1,1])-A3*(y-s1[1,2])=0;  line4:=B4*(x-s[1,1])-A4*(y-s[1,2])=0;Sol:=solve({line3, line4}); x0:=simplify(rhs(Sol[1])); y0:=simplify(rhs(Sol[2]));

if is(s1[1]=s[1]) then    u:=s1[2]-[x0,y0]; v:=s[2]-[x0,y0]; else

u:=s1[1]-[x0,y0]; v:=s[1]-[x0,y0];  fi;  fi;  fi;

Sign:=signum(u[1]*v[2]-u[2]*v[1]);   phi:=Sign*arccos(expand(rationalize(simplify((u[1]*v[1]+u[2]*v[2])/sqrt(u[1]^2+u[2]^2)/sqrt(v[1]^2+v[2]^2)))));      S11:=expand(rationalize(simplify(map(x->f(op(x), phi), S1))));

if n3<>0 then T11:={seq(expand(rationalize(simplify(map(x->f(op(x), phi), T1[i])))), i=1..nops(T1))}; fi;

if n3=0 then  if is(S11=expand(rationalize(simplify(S2))))  then T:=[rotation, [x0,y0], phi]; return true fi;  else

if is(S11=expand(rationalize(simplify(S2)))) and  is(T11=expand(rationalize(simplify(T2)))) then

T:=[rotation, [x0,y0], phi]; return true fi;  fi;

fi;

od;

# Checking for reflection or glide reflection

for s in S do

AB:=s1[2]-s1[1]; CD:=s[2]-s[1];

if is(AB[1]*CD[2]-AB[2]*CD[1]=0) then  M:=(s1[2]+s[1])/2;

if  is(AB[1]*CD[1]+ AB[2]*CD[2]>0) then  A:=AB[2]; B:=-AB[1];    Line1:=A*(x-M[1])+B*(y-M[2])=0;  else

A:=AB[1]; B:=AB[2];  Line2:=A*(x-M[1])+B*(y-M[2])=0; fi;

else     u:=[AB[1]+CD[1], AB[2]+CD[2]];  A:=u[2]; B:=-u[1];     M:=[(s1[1,1]+s[1,1])/2, (s1[1,2]+s[1,2])/2]; Line3:=A*(x-M[1])+B*(y-M[2])=0;   fi;    C:=-A*M[1]-B*M[2];  h:= simplify(expand(rationalize(s[1]-g(op(s1[1])))));    S11:=expand(rationalize(simplify(map(x->g(op(x))+h, S1))));

if n3<>0 then T11:={seq(expand(rationalize(simplify(map(x->g(op(x))+h, T1[i])))), i=1..nops(T1))}; fi;

if n3=0 then   if is(S11=expand(rationalize(S2))) then

Eq:=A*x+B*y+C=0; Eq1:=`if`(is(coeff(lhs(Eq), y)<>0), y=solve(Eq, y),  x=solve(Eq, x));

if h=[0,0] then  T:=[reflection, Eq1] else T:=[glide_reflection,Eq1,convert(h, Vector)] fi; return true fi; else

if is(S11=expand(rationalize(S2))) and is(T11=expand(rationalize(T2))) then

Eq:=A*x+B*y+C=0; Eq1:=`if`(is(coeff(lhs(Eq), y)<>0), y=solve(Eq, y),  x=solve(Eq, x));

if h=[0,0] then T:=[reflection, Eq1] else

T:=[glide_reflection,Eq1,convert(h, Vector)] fi; return true fi;  fi;

od;

T:='T';   false;

end proc:

Three simple examples:

AreIsometric({[4, 0], [7, 4], [14, 0]}, {[4, 14], [9, 14], [10, 6]});  T;

AreIsometric({[2, 0], [2, 2], [5, 0]}, {[3, 3], [3, 6], [5, 3]});  T;

S1 := {[[5, 5], [5, 20], [10, 15], [15, 20], [15, 5]]}:
S2 := {[[21, 11], [30, 23], [31, 16], [38, 17], [29, 5]]}:
S3 := {[[50, 23], [41, 11], [51, 16], [49, 5], [58, 17]]}:
AreIsometric(S1, S2); T; AreIsometric(S1, S3);

plots[display](plottools[curve](op(S1), thickness = 2, color = green), plottools[curve](op(S2), thickness = 2, color = green), plottools[curve](op(S3), thickness = 2, color = red), scaling = constrained, view = [-1 .. 60, -1 .. 25]);  # Green sets are isometric,  green and red sets aren't

Example with animation:

S1:={[0,0],[-1,2],[2,4],[4,2]}:  S2:={[8,3],[6,6],[8,8],[10,4]}:

AreIsometric(S1, S2);  T;

with(plots): with(plottools):

#  For clarity, instead of points polygons depicted with vertices at these points

N:=50:

A:=seq(rotate(polygon([[0,0],[-1,2],[2,4],[4,2]], color=blue), (k*Pi)/(2*N), [3,7]), k=0..N):  B:=polygon([[8,3],[6,6],[8,8],[10,4]], color=green):  E:=line([3,7], [6,6], color=black, linestyle=2):

C:=seq(rotate(line([3,7], [2,4], color=black, linestyle=2), (k*Pi)/(2*N), [3,7]), k=0..N):  L:=curve([[3,7],[2,4], [-1,2], [0,0],[4,2], [2,4]], color=black, linestyle=2):  T:=textplot([3, 7.2, "Center of rotation"]):

Frames:=seq(display(A[k], B, E,T,L, C[k]), k=1..N+1):

display(seq(Frames[k],k=1..N+1), insequence=true, scaling=constrained);

Finding unique solutions to the problem of Queens (m chess queens on an n×n chessboard not attacking one another). Used the procedures  Queens  and  QueenPic  . See  http://www.mapleprimes.com/posts/200049-Queens-Problem-And-Its-Visualization-With-Maple

Queens(8, 8);  M := [ListTools[Categorize](AreIsometric, S)]:

nops(M);

seq(op(M[k])[1], k = 1 .. %);   # 12 unique solutions from total 92 solutions

QueensPic([%], 4);  #  Visualization of obtained solutions

Finding unique solutions to the problem "Polygons of matches"  (all polygons with specified perimeter and area). See http://www.mapleprimes.com/posts/142884-Polygons-Of-The-Matches

N := 12: S := 6: Polygons(N, S);

M := [ListTools[Categorize](AreIsometric, T)]:

n := nops(M);

seq([op(M[k][1])], k = 1 .. n);  #  7 unique solutions from total 35 solutions with perimeter 12 and area 6

Visual([%], 4);  #  Visualization of obtained solutions

## The largest incircle and the smallest co...

This message contains two procedures . The first procedure called  LargestIncircle  symbolically looks for the largest circle inscribed in a convex polygon . Polygon must be specified as a list of lists of its vertices . If the polygon does not have parallel sides , the problem has a unique solution , otherwise may be several solutions . It is easy to prove that there is always a maximum circle that touches three sides of the polygon . Procedure code based on this property of the optimal solution.

The second procedure called  MinCoveringCircle  symbolically seeking smallest circle that encloses a given set of points. The set of points can be anyone not necessarily being vertices of a convex polygon. Procedure code based on the following geometric properties of the optimal solution , which can easily be proved :

1) The minimum covering circle is unique.

2) The minimum covering circle of a set Set can be determined by at most three points in  Set  which lie on the boundary of the circle. If it is determined by only two points, then the line segment joining those two points must be a diameter of the minimum circle. If it is determined by three points, then the triangle consisting of those three points is not obtuse.

Code of the 1st procedure:

restart;

LargestIncircle:=proc(L::listlist)

local n, x, y, P, S, T, Eqs, IsInside, OC, Pt, E, R, t, Sol, P0, R1, Circle;

uses combinat, geometry;

n:=nops(L); _EnvHorizontalName := x; _EnvVerticalName := y;

P:=[seq(point(A||i, L[i]), i=1..n)];

if nops(convexhull(P))<n then error "L should be a convex polygon" fi;

S:={seq(line(l||i,[A||i, A||(i+1)]), i=1..n-1), line(l||n,[A||n, A||1])};

Eqs:=map(lhs,{seq(Equation(l||i), i=1..n)});

T:=choose({seq(l||i, i=1..n)}, 3);

IsInside:=proc(Point, OC, Eqs)

convert([seq(is(eval(Eqs[i], {x=Point[1], y=Point[2]})*eval(Eqs[i], {x=OC[1], y=OC[2]})>0), i=1..n)], `and`);

end proc;

R:=0; point(Pt,[x,y]);

for t in T do

solve([distance(Pt,t[1])=distance(Pt,t[2]),      distance(Pt,t[2])=distance(Pt,t[3])],[x,y]);

Sol:=map(a->[rhs(a[1]), rhs(a[2])], %);

P0:=op(select(b->IsInside(b, OC, Eqs), Sol)); if P0<>NULL then R1:=distance(point(E,P0),      t[1]);

if convert([seq(is(R1<=distance(E, i)), i=S minus t)], `and`) and is(R1>R) then R:=R1;      Circle:=[P0, R] fi; fi;

od;

Circle;

end proc:

Code the 2nd procedure:

MinCoveringCircle:=proc(Set::{set, list}(list))

local S, T, t, d, A, B, C, P, R, R1, O1, Coord, L;

uses combinat, geometry;

if type(Set, list) then S:=convert(Set, set) else S:=Set fi;

T:=choose(S, 2);

for t in T do

d:=simplify(distance(point(A, t[1]), point(B, t[2]))); midpoint(C, A, B);

if convert([seq(is(distance(point(P, p), C)<=d/2), p in S minus {t[1], t[2]})], `and`)   then return            [coordinates(C), d/2] fi;

od;

T:=choose(S, 3);

R:=infinity;

for t in T do

point(A, t[1]); point(B, t[2]); point(C, t[3]);

if is(distance(A, B)^2<=distance(B, C)^2+distance(A, C)^2 and        distance(A,C)^2<=distance(B,        C)^2+distance(A, B)^2 and distance(B, C)^2 <= distance(A, C)^2 + distance(A, B)^2) then

if convert([seq(is(distance(point(P, p), O1)<=R1), p in S minus {t[1], t[2], t[3]})], `and`) and is(R1<R) then    R:=R1; L:=[Coord, R] fi; fi;

od;

L;

end proc:

Examples of using the first procedure:

L:=[[0,2],[1,4],[4,4],[5,1],[3,0]]:

C:=LargestIncircle(L);

A:=plottools[disk](op(C), color=green):

B:=plottools[polygon](L, style=line, thickness=2):

plots[display](A,B, scaling=constrained);

The procedure can be used to find the largest circle inscribed in a convex shape, bounded by the curves, if these curves are approximated by broken lines:

L:=[seq([k,k^2], k=-1..2, 0.1)]:

C:=LargestIncircle(L);

A:=plottools[disk](op(C), color=green):

B:=plottools[polygon](L, style=line, thickness=2):

plots[display](A,B, scaling=constrained);

Examples of using the second procedure:

L:=[[0,2],[1,4],[4,4],[5,1],[3,0]]:

C:=MinCoveringCircle(L);

A:=plottools[circle](op(C), color=red, thickness=2):

B:=plottools[polygon](L, color=cyan):

plots[display](A,B);

The smallest circle covering 30 random points:

> roll:=rand(1..10):

L:=[seq([roll(), roll()], i=1..30)]:

C:=MinCoveringCircle(L);

A:=plottools[circle](op(C), color=red, thickness=2):

B:=seq(plottools[disk](L[i],0.07,color=blue), i=1..nops(L)):

C1:=plottools[disk](C[1], 0.07, color=red):

T:=plots[textplot]([C[1,1]-0.2, C[1,2]-0.2, "C"], font=[TIMES,ROMAN,16]):

plots[display](A, B, C1, T);

 First 6 7 8 9 10 11 Page 8 of 11
﻿