r1:=[[1,2,3,4,5,6,7,8,9,10,11,12]];
r2:=[[1,13,14,15,16,17,7,18,19,20,21,22]];
r3:=[[4,30,29,20,28,27,10,26,25,15,24,23]];

perminteger:=proc(L::list,n::integer)
#L should be a list representing a cyclic permutation
#this returns the image of the integer under the
#normal cyclic action of L on the integers
#(it will have no effect if n does not occur in L)
 local i,j;
 for i from 1 to nops(L)-1 do
  if L[i]=n then RETURN(L[i+1]); fi;
 od; ## for i
 if L[nops(L)]=n then RETURN(L[1]); fi;
RETURN(n);
end; 

permlist:=proc(L,v::list)
#L should be a list of lists representing a permutation in
#disjoint cycle notation 
#this returns the image of the list under the
#normal (left-to-right) permutation action of L on the integers
 local i,j,vv;
 vv:=v;
 if nops(L)=1 then 
  for i from 1 to nops(v) do vv[i]:=perminteger(op(L),v[i]); 
  od;
 fi;
 if nops(L)>1 then
  for j from 1 to nops(L) do
   vv:=permlist([L[j]],vv);
  od; ## for j
 fi;
RETURN(vv);
end;

permultiply:=proc(L::list)
#L is a list of permutations and/or 
#their integral powers, each permutation written
#in disjoint cycle notation,
#eg, L:=[perm1^n1,perm2^n2,...,permk^nk];
#This returns the product 
 local i,j,g,gtemp,gtemp1,gtemp2,newbase,base,power,temp_power;
 g:=L[1];
 if (nops(L)=1 and type(g,list)) then
#print(`1`,L);
   RETURN(g); 
 fi;
 if nops(L)=1 then
  power:=op(2,L[1]);
  base:=op(1,L[1]);
  if power=0 then RETURN([]); fi;
  if power>0 then g:=base; fi;
  if power<0 then 
     temp_power:=-power;
     newbase:=group[invperm](base);
     g:=newbase;
     base:=newbase; 
     power:=temp_power; 
  fi; 
#print(`2`,power,base);
  for i from 1 to power-1 do
   gtemp:=g;
   g:=group[mulperms](base,gtemp);
  od;
 RETURN(g);
 fi;
 if nops(L)>1 then
#print(`3`,L);
  g:=[]; 
  for j from 1 to nops(L) do
   gtemp1:=g;  
   gtemp2:=permultiply([L[j]]);
   g:=group[mulperms](gtemp1,gtemp2);
  od;
 RETURN(g);
 fi; 
end;

viewmatrix:=proc(v::list)
local A;
A:=array([[v[1],v[2],v[3],v[4],v[5],v[6],v[7],v[8],v[9],v[10],v[11],v[12]],
[v[1],v[13],v[14],v[15],v[16],v[17],v[7],v[18],v[19],v[20],v[21],v[22]],
[v[4],v[23],v[24],v[15],v[25],v[26],v[10],v[27],v[28],v[20],v[29],v[30]]]);
RETURN(A);
end;

P0:=vector([1,0]);
C0:=vector([0,0]);
A0:=matrix([[1,2,3,4,5,6,7,8,9,10,11,12],
 [a11,13,14,15,16,17,a17,18,19,20,21,22],
 [a14,23,24,a24,25,26,a1_10,27,28,a2_10,29,30]]);
 a11:=A0[1,1];
 a17:=A0[1,7];
 a14:=A0[1,4];
 a24:=A0[2,4];
 a1_10:=A0[1,10];
 a2_10:=A0[2,10];
A0:=matrix([[1,2,3,4,5,6,7,8,9,10,11,12],
 [a11,13,14,15,16,17,a17,18,19,20,21,22],
 [a14,23,24,a24,25,26,a1_10,27,28,a2_10,29,30]]);

rotation:=proc(j,N)
local A;
A:=matrix([[cos(2*j*Pi/N),-sin(2*j*Pi/N)],
    [sin(2*j*Pi/N),cos(2*j*Pi/N)]]);
RETURN(A);
end;

for j from 1 to 12 do 
 Pa.j:=evalm(rotation(j,12)&*(P0-C0)+C0):
 od:
for j from 1 to 12 do 
 Pta.j:=[10*evalf(Pa.j[1]),10*evalf(Pa.j[2]),0]; od:
for j from 1 to 12 do
 Ptb.j:=[0,10*evalf(Pa.j[1]),10*evalf(Pa.j[2])]: 
 od:
for j from 1 to 12 do
 Ptc.j:=[10*evalf(Pa.j[1]),0,10*evalf(Pa.j[2])]: 
 od:

mod12:=proc(n::integer)
 local i;
 i:=`mod`(n,12);
 if i=0 then i:=12; fi;
 RETURN(i);
 end;

move1:=proc(A::array,power::integer)
 local i,j,AA;
 AA:=A;
 for i from 1 to 12 do AA[1,i]:=A[1,mod12(i+power)]; od;
 for i from 1 to 12 do AA[2,i]:=A[2,i]; od;
 for i from 1 to 12 do AA[3,i]:=A[3,i]; od;
 AA[2,1]:=AA[1,1];
 AA[2,7]:=AA[1,7];
 AA[3,1]:=AA[1,4];
 AA[3,7]:=AA[1,10];
 RETURN(AA);
end:

move2:=proc(A::array,power::integer)
 local i,j,AA;
 AA:=A;
 for i from 1 to 12 do AA[2,i]:=A[2,mod12(i+power)]; od;
 for i from 1 to 12 do AA[1,i]:=A[1,i]; od;
 for i from 1 to 12 do AA[3,i]:=A[3,i]; od;
 AA[1,1]:=AA[2,1];
 AA[1,7]:=AA[2,7];
 AA[3,4]:=AA[2,4];
 AA[3,10]:=AA[2,10];
 RETURN(AA);
end:

move3:=proc(A::array,power::integer)
 local i,j,AA;
 AA:=A;
 for i from 1 to 12 do AA[3,i]:=A[3,mod12(i+power)]; od;
 for i from 1 to 12 do AA[1,i]:=A[1,i]; od;
 for i from 1 to 12 do AA[2,i]:=A[2,i]; od;
 AA[2,4]:=AA[3,4];
 AA[1,4]:=AA[3,1];
 AA[2,10]:=AA[3,10];
 AA[1,10]:=AA[3,7];
 RETURN(AA);
end:

#The following procedure plots the globe in 3d 
showsphere:=proc(M::array)
 local label,plota,plotb,plotc,i,j,spher;
 for i from 1 to 12 do
 label:=convert(M[1,i],string); 
 plota.i:=textplot3d([op(Pta.i),label],axes=none,color=magenta):od:
 for i from 1 to 12 do
 label:=convert(M[2,i],string); plotb.i:=textplot3d([op(Ptb.i),label],axes=none,color=blue):
 od:
 for i from 1 to 12 do
 label:=convert(M[3,i],string); plotc.i:=textplot3d([op(Ptc.i),label],axes=none,color=red):
 od:
 spher:=sphereplot(9,theta=0..2*Pi,phi=0..Pi,
 style=wireframe,numpoints=50,color=grey):
 display3d([plota1,plota2,plota3,plota4,plota5,
 plota6,plota7,plota8,plota9,plota10,plota11,plota12,
 plotb1,plotb2,plotb3,plotb4,plotb5,
 plotb6,plotb7,plotb8,plotb9,plotb10,plotb11,plotb12,
 plotc1,plotc2,plotc3,plotc4,plotc5,
 plotc6,plotc7,plotc8,plotc9,plotc10,plotc11,plotc12,spher]);
end:

puzzlemap:=proc(v0::list)
 local i,Asia,NA1,NA2,Pac,Eu,eq,label;
 for i from 1 to 7 do
  label:=convert(v0[i],string): 
  NA1.i:=textplot([1,8-i,label],axes=none):
 od:
 label:=convert(v0[1],string): 
 Eu.1:=textplot([4,7,label],axes=none):
 label:=convert(v0[1],string): 
 Asia.1:=textplot([7,7,label],axes=none):
 label:=convert(v0[1],string): 
 Pac.1:=textplot([10,7,label],axes=none):
 for i from 2 to 6 do
  label:=convert(v0[i+11],string): 
  Eu.i:=textplot([4,8-i,label],axes=none):
 od:
 label:=convert(v0[7],string): 
 Eu.7:=textplot([4,1,label],axes=none):
 Asia.7:=textplot([7,1,label],axes=none):
 Pac.7:=textplot([10,1,label],axes=none):
 for i from 2 to 6 do
  label:=convert(v0[14-i],string): 
  Asia.i:=textplot([7,8-i,label],axes=none):
 od:
 for i from 2 to 6 do
  label:=convert(v0[24-i],string): 
  Pac.i:=textplot([10,8-i,label],axes=none):
 od:
 label:=convert(v0[10],string): 
 Asia.4:=textplot([7,4,label],axes=none):
 for i from 1 to 7 do
  label:=convert(v0[i],string): 
  NA2.i:=textplot([13,8-i,label],axes=none):
 od:
 label:=convert(v0[30],string): 
 eq.1:=textplot([0,4,label],axes=none):
 for i from 2 to 3 do
  label:=convert(v0[i+21],string): 
  eq.i:=textplot([i,4,label],axes=none):
 od:
 for i from 4 to 5 do
  label:=convert(v0[i+21],string): 
  eq.i:=textplot([i+1,4,label],axes=none):
 od:
 for i from 6 to 7 do
  label:=convert(v0[i+21],string): 
  eq.i:=textplot([i+2,4,label],axes=none):
 od:
 for i from 8 to 9 do
  label:=convert(v0[i+21],string): 
  eq.i:=textplot([i+3,4,label],axes=none):
 od:
 label:=convert(v0[23],string): 
 eq.10:=textplot([14,4,label],axes=none):
 display([NA1.1,NA1.2,NA1.3,NA1.4,NA1.5,NA1.6,NA1.7,
 Eu1,Eu2,Eu3,Eu4,Eu5,Eu6,Eu7,
 Asia1,Asia2,Asia3,Asia4,Asia5,Asia6,Asia7,
 Pac1,Pac2,Pac3,Pac4,Pac5,Pac6,Pac7,
 NA2.1,NA2.2,NA2.3,NA2.4,NA2.5,NA2.6,NA2.7,
 eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10]);
end;

print(`"circles" is now loaded`);
print(` `);
print(`This creates global variables: A0, P0, C0, r1, r2, r3, Pa1 Pa2, ...`);
print(`Pta1 Pta2, ...,Ptb1 Ptb2, ...,Ptc1 Ptc2, ...`);
print(` `);
print(`To use "showsphere" you must load linalg and plots`);
print(`To use "permultiply" you must load group`);
print(`To use "puzzlemap" you must load plots`);