Unfortunately Maple is weak in the transformations of trigonometric expressions. Only after several attempts I was able to prove that  a[2]=aa[2] . 

Probably the easiest way to prove that two numerical trigonometric expressions are equal (not fully correct to the position of pure mathematics) is to calculate their approximate values with high accuracy and then apply the  is  command.

Download TT.mw

Check:=proc(A::Matrix)
local m, n;
m:=op([1,1],A); n:=op([1,2],A);
{seq(`+`(convert(A[i],list)[ ]), i=1..m), seq(`+`(convert(A[..,j],list)[ ]), j=1..n)};
if nops(%)=1 then true else false fi;
end proc:

Example of use:

A:=<2,7,6; 9,5,1; 4,3,8>;
Check(A);

In fact, your functions   f₁(t₁,t₂,t₃),  f₂(t₁,t₂,t₃), f₃(t₁,t₂,t₃)  define a mapping  R^3  into  R^3. In the example  the cuboid and its image shown under the mapping  f  (a rotation and translation):

restart;
with(plottools): with(plots):
A:=<1/3,2/3,2/3; -2/3,-1/3,2/3; 2/3,-2/3,1/3>: # Matrix of a rotation
f:=unapply(convert(A.<x,y,z>+<3,3,3>, list), (x,y,z)); # A rotation and translation mapping procedure
F:=transform(f):
p:=display(cuboid([0,0,0], [1,2,2.5], color=khaki)):
display(p, F(p), scaling=constrained, axes=normal, orientation=[-15,80]);

  f := (x, y, z) -> [(1/3)*x+(2/3)*y+(2/3)*z+3, -(2/3)*x-(1/3)*y+(2/3)*z+3, (2/3)*x-(2/3)*y+(1/3)*z+3] 

a:=(-2*theta+1)*phi/(phi-1)+theta;
b:=(2*theta-1)*phi/(-phi+1)+theta;
c:=-numer(op(1,a))/``(-denom(op(1,a)))+theta;

Edit.  If you need to remove the parentheses in the denominator, you can write additional   applyop(expand, [1,3], c);

Another one-line solution (the best variant of the solutions):

c:=applyop(t->-t, {[1,1], [1,3,1]}, a);

restart;
Plane:=proc(L::listlist)
uses LinearAlgebra;
sort(Determinant(Matrix([[x,y,z]-~(L[1]),L[2]-L[1],L[3]-L[1]])));
(%=0)/igcd(coeffs(%));
end proc:

Examples of use:
Plane([[1,-1,3],[-15,-17,11],[2,1,0]]);
L := [[[-12, 2, -1], [-11, 1, -5], [-10, -2, 3]], [[-12, 2, -1], [-11, 1, -5], [-10, 6, 3]], [[-12, 2, -1], [-11, 1, -5], [-9, 5, -7]], [[-12, 2, -1], [-11, 1, -5], [-9, 8, -4]], [[-12, 2, -1], [-11, 1, -5], [-7, -6, -2]], [[-12, 2, -1], [-11, 1, -5], [-7, -2, -8]], [[-12, 2, -1], [-11, 1, -5], [-7, -2, 6]], [[-12, 2, -1], [-11, 1, -5], [-7, 3, -9]], [[-12, 2, -1], [-11, 1, -5], [-7, 3, 7]], [[-12, 2, -1], [-11, 1, -5], [-7, 6, -8]], [[-12, 2, -1], [-11, 1, -5], [-7, 9, 3]], [[-12, 2, -1], [-11, 1, -5], [-7, 10, -2]], [[-12, 2, -1], [-11, 1, -5], [-6, -4, -7]], [[-12, 2, -1], [-11, 1, -5], [-6, -4, 5]], [[-12, 2, -1], [-11, 1, -5], [-6, 8, -7]], [[-12, 2, -1], [-11, 1, -5], [-6, 8, 5]], [[-12, 2, -1], [-11, 1, -5], [-4, -6, 3]], [[-12, 2, -1], [-11, 1, -5], [-4, -2, -9]], [[-12, 2, -1], [-11, 1, -5], [-4, -2, 7]], [[-12, 2, -1], [-11, 1, -5], [-4, 6, -9]], [[-12, 2, -1], [-11, 1, -5], [-4, 6, 7]], [[-12, 2, -1], [-11, 1, -5], [-4, 10, 3]], [[-12, 2, -1], [-11, 1, -5], [-2, -6, 3]], [[-12, 2, -1], [-11, 1, -5], [-2, 6, -9]], [[-12, 2, -1], [-11, 1, -5], [-2, 6, 7]], [[-12, 2, -1], [-11, 1, -5], [-2, 10, 3]], [[-12, 2, -1], [-11, 1, -5], [3, 5, -7]], [[-12, 2, -1], [-11, 1, -5], [3, 8, -4]], [[-12, 2, -1], [-11, 1, -5], [4, -2, 3]], [[-12, 2, -1], [-11, 1, -5], [4, 6, 3]], [[-12, 2, -1], [-11, 1, -5], [5, 6, -2]], [[-12, 2, -1], [-11, 1, 3], [-10, -2, -5]], [[-12, 2, -1], [-11, 1, 3], [-10, 6, -5]], [[-12, 2, -1], [-11, 1, 3], [-9, 5, -7]], [[-12, 2, -1], [-11, 1, 3], [-9, 8, -4]], [[-12, 2, -1], [-11, 1, 3], [-7, -6, -2]], [[-12, 2, -1], [-11, 1, 3], [-7, -2, -8]], [[-12, 2, -1], [-11, 1, 3], [-7, -2, 6]], [[-12, 2, -1], [-11, 1, 3], [-7, 6, -8]], [[-12, 2, -1], [-11, 1, 3], [-7, 9, -5]]]:
Plane~(L);

I do not understand your code, and wrote a new recursive procedure that solves your problem:

bit:=proc(n::posint,k::nonnegint)
if k=0 then return [[0$n]] fi;
if k=n then return [[1$n]] fi;
if k>0 and k<n then {seq(seq(`if`(i[j]=0,subsop(j=1,i),NULL), j=1..n), i=bit(n,k-1))} fi;
[op(%)];
end proc:

Example of use: 

for k from 0 to 5 do
bit(5, k);
od;

Edit.  Of course, the easiest way to solve this problem, is the usage  combinat:-permute  command, for example:

combinat:-permute([1,1,0,0,0]);

[[1, 1, 0, 0, 0], [1, 0, 1, 0, 0], [1, 0, 0, 1, 0], [1, 0, 0, 0, 1], [0, 1, 1, 0, 0], [0, 1, 0, 1, 0], [0, 1, 0, 0, 1], [0, 0, 1, 1, 0], [0, 0, 1, 0, 1], [0, 0, 0, 1, 1]]

Here is another version without recursion:

bit1:=proc(n::posint,k::nonnegint)
local NestedSeq, L;
uses ListTools;
if k=0 then return cat(0$n) fi;
NestedSeq:=proc(Expr::uneval, L::list)
local S;
eval(subs(S=seq, foldl(S, Expr, op(L))));
end proc:
[NestedSeq(cat(seq(`if`(`or`(seq(s=i||m,m=1..k)),1,0), s=1..n)),[Reverse([seq(i||p=i||(p-1)+1..n-k+p,p=2..k)])[],i1=1..n-k+1])];
end proc:

Example of use:

bit1(5,2);

Zeineb, in your example, you missed the square  [22, 33,  42, 31] .

The procedure  Anew  solves the problem for any matrix, in which both the dimensions of at least 3.

Anew:=proc(A::Matrix)
local m, n;
m:=op([1,1],A);  n:=op([1,2],A);
if m<3 or n<3 then return NULL fi;
if n::even then if m::even then
return Matrix([seq(op([seq([A[i,j],A[i+1,j+1],A[i+2,j],A[i+1,j-1]], j=2..n-2,2), seq([A[i+1,j],A[i+2,j+1],A[i+3,j],A[i+2,j-1]], j=3..n-1,2)]), i=1..m-2,2)]) else
return Matrix([seq(op([seq([A[i,j],A[i+1,j+1],A[i+2,j],A[i+1,j-1]], j=2..n-2,2), seq([A[i+1,j],A[i+2,j+1],A[i+3,j],A[i+2,j-1]], j=3..n-1,2)]), i=1..m-4,2), seq([A[m-2,j],A[m-2+1,j+1],A[m-2+2,j],A[m-2+1,j-1]], j=2..n-2,2)]) fi; fi;
if n::odd then if m::even then
return Matrix([seq(op([seq([A[i,j],A[i+1,j+1],A[i+2,j],A[i+1,j-1]], j=2..n-1,2), seq([A[i+1,j],A[i+2,j+1],A[i+3,j],A[i+2,j-1]], j=3..n-2,2)]), i=1..m-2,2)]) else
Matrix([seq(op([seq([A[i,j],A[i+1,j+1],A[i+2,j],A[i+1,j-1]], j=2..n-1,2), seq([A[i+1,j],A[i+2,j+1],A[i+3,j],A[i+2,j-1]], j=3..n-2,2)]), i=1..m-4,2), seq([A[m-2,j],A[m-2+1,j+1],A[m-2+2,j],A[m-2+1,j-1]], j=2..n-1,2)]) fi; fi;
end proc:

Example of use:

interface(rtablesize=infinity):
A := Matrix(5,10, [$ 1..50]); 
Anew(A);

For a given triangle it is not a very difficult task. We only need to set the ranges for parameterization:

restart;
with(geom3d):
point(A,1,0,0), point(B,1/2,0,sqrt(3)/2),point(C,0,sqrt(3)/2,1/2), point(OO,0,0,0):
plane(P1,[A,C,OO]), plane(P2,[B,C,OO]), plane(P3,y=tan(phi)*x, [x,y,z]):
Plane1:=Equation(P1,[x,y,z]):
Plane2:=Equation(P2,[x,y,z]):
Tr:={x=cos(phi)*sin(theta),y=sin(phi)*sin(theta),z=cos(theta)}:
eval({Plane1,y=tan(phi)*x}, Tr):
theta1:=solve(eval({Plane1,y=tan(phi)*x}, Tr), theta,explicit):
eval({Plane2,y=tan(phi)*x}, Tr):
theta2:=solve(eval({Plane2,y=tan(phi)*x}, Tr), theta,explicit):
plots:-display(plottools:-sphere(color="LightBlue", style=patch),plot3d([1.01*cos(phi)*sin(theta),1.01*sin(phi)*sin(theta),1.01*cos(theta)], phi=0..Pi/2, theta=eval(theta,theta2)..eval(theta,theta1),color="Green", style=surface), axes=none);
                     

plots:-display(plottools:-circle(), plots:-textplot([0, 0.9, infinity], font=[times, bold, 22]), axes=none);

I do not quite understand your problem. If you want to avoid premature calculation, you can write down all the constants as symbols.
 

Example:

`12`*x^`-4`*y^`2`/`3`/x^`6`/y^`-5`;
evalindets(%, 'symbol', parse);

phephoi := proc (A, B)
local check, i, j, C, E;
uses ListTools;
C := A;
for i from 1 to nops(A) do
   check := 0; 
   for j from 1 to nops(B) do
      if B[j] = A[i] then check := 1; break; end if;
   end do; 
   if check = 0 then C := subsop(Search(A[i],C) = NULL, C); end if;
end do;
C;
end proc:

Example of use:

A := [1, 2, 3, 5, 6, 7];
B := [2, 4, 7, 11, 8];
phephoi(A, B);

If you are interested in the image of a clock face with Maple, look at this post

http://www.mapleprimes.com/posts/143906-Visualization-Of-Clock-And-Degrees-Problem

vv, this is a very elegant solution, vote up!

Here is another solution, based on the procedure named  ArcOnSphere . The procedure draws the shortest arc on a sphere, connecting the specified points (the points do not have to be opposites). Required parameters:  A  and  B  are two points on the sphere, C is the center of the sphere (by default  C = [0, 0, 0] ), G are the graphics options for the space curve  (by default  G=[color=red, thickness=3, scaling=constrained] ).

restart;
ArcOnSphere:=proc(A::list, B::list, C::list:=[0,0,0], G::list(equation):=[color=red, thickness=3, scaling=constrained])
local Av, Bv, Cv, CA, CB, V, M, Phi;
uses Student[LinearAlgebra], plots;
Av, Bv, Cv:= op(convert~([A, B, C], Vector));
CA:=Av-Cv; CB:=Bv-Cv;
V:=CA &x CB;
M:=RotationMatrix(alpha, V);
Phi:=arccos(CA.CB/Norm(CA)/Norm(CB));
spacecurve(M.CA+Cv, alpha=0..Phi, op(G));
end proc:

Example of use:

L:=[[sqrt(3)/2,0,1/2], [0,1/2,sqrt(3)/2], [0,1,0], [1/sqrt(2),1/sqrt(2),0]]:  # 4 points on the unit sphere
L1, L2, L3, L4:=seq(ArcOnSphere(op(L[i..i+1])), i=1..nops(L)-1), ArcOnSphere(L[nops(L)], L[1]):  # The arcs connecting these points
Labels:=["A", "B", "C", "D"]:
T:=plots[textplot3d]([seq([op(L[i]+0.05*L[i]), Labels[i]], i=1..nops(L))], font=[times,bold,16]):
plots[display](L1, L2, L3, L4, T, plottools[sphere](style=surface));

interface(rtablesize=infinity):  # It is necessary that the matrices be displayed explicitly
A, B := LinearAlgebra:-RandomMatrix(15),  LinearAlgebra:-RandomMatrix(15);  # Specifying of two matrices
A+B,  3*A,  A^2,  A.B,  A^(-1);

Please note that for multiplying matrices, we use a dot rather than a sign  "*" 

DuplicationMatrix:=proc(n)
local A, a, L1, L2;
uses ListTools;
A:=Matrix(n, symbol=a, shape=symmetric);
L1:=[seq(seq(A[i,j], j=1..n), i=1..n)];
L2:=convert(L1, set);
Matrix(n^2,n*(n+1)/2, {seq(seq((i,j)=`if`(j=Search(L1[i], L2),1, 0), j=1..n*(n+1)/2), i=1..n^2)});
end proc:

Example of use:
C:=DuplicationMatrix(3);

Check.

A, V1, V2:=<a,b,c; b,d,e; c,e,f>, <a,b,c,d,e,f>, <seq(seq(A[i,j], j=1..3), i=1..3)>;
C.V1=V2;