Remove["Global`*"]
(*This skript exemplifies in real space that the number of*)
(*n-fold rotational symmetries that can be distinguished*)
(*is equal to n,all other n-fold rotational symmetries*)
(*coincide with a continuous rotational symmetry (transverse*)
(*isotropy)*)

(*This is of course not a full proof.For this purpose,one can*)
(*leave the roation angle unspecified and show that the rank*)
(*of the coefficient matrix only lowers for specific angles,*)
(*which happen to be fractions of 2 Pi/m with 1≤m≤n,*)
(*which is ultimatly a consequence of the fact that any power n*)
(*of a sin or cos function can be written as a linear comb.*)
(*of sin and cos functions with the argument divided by m.*)
(*Hermanns proof is much more elegant,but does not capture*)
(*the equivalence of the action of C3,D3 and C4,D4.*)

order=4;(*tensorial order*)comps=3^order;(*no.of indep.comps.in absence of symms.*)(*create a symbolic array of all components*)CC=Array[C,ConstantArray[3,order]];

(*Q=RotationMatrix[kappa,{0,0,1}]//TrigToExp*)
(*definition of the product Q_im Q_jn... Q_lp C_ij..l*)
(*for any tensorial order*)

Rayleigh[rotmat_,arg_,depth_]:=Module[{tmp},tmp=arg;
For[i=1,i≤depth,i++,tmp=Transpose[tmp,1<->i];
tmp=rotmat.tmp;
tmp=Transpose[tmp,1<->i];];
tmp]

(*The following creates generators of the symm.group*)
(*and defines the invariances.Uncomment for the group*)
(*of interest.*)

(*Triclinic*)
Q=IdentityMatrix[3];
Q2=Q;

(*Monoclinic*)
Q=RotationMatrix[2 Pi/2,{0,0,1}];
Q2=IdentityMatrix[3];

(*Orthotropic*)
Q=RotationMatrix[2 Pi/2,{0,0,1}];
Q2=RotationMatrix[2 Pi/2,{0,1,0}];

(*Tetragonal (D4)*)
Q=RotationMatrix[2 Pi/4,{0,0,1}];
Q2=IdentityMatrix[3];

(*Piezo-Tetragonal (C4)*)
Q=RotationMatrix[2 Pi/4,{0,0,1}];
Q2=RotationMatrix[Pi,{0,1,0}];

(*Trigonal (D3)*)
Q=RotationMatrix[2 Pi/3,{0,0,1}];
Q2=IdentityMatrix[3];

(*Piezo-Trigonal (C3)*)
Q=RotationMatrix[2 Pi/3,{0,0,1}];
Q2=RotationMatrix[Pi,{0,1,0}];

(*Hexagonal*)
Q=RotationMatrix[2 Pi/6,{0,0,1}];
Q2=IdentityMatrix[3];

(*Cubic*)
Q=RotationMatrix[2 Pi/3,{1,1,1}];
Q2=RotationMatrix[Pi,{0,1,0}];

(*Isotropic*)
Q=RotationMatrix[2 Pi/6,{0,0,1}];
Q2=RotationMatrix[2 Pi/6,{0,1,0}];

rotsymmetries=Join[Flatten[Rayleigh[Q,CC,4]-CC],Flatten[Rayleigh[Q2,CC,4]-CC]];

symmetries=Join[rotsymmetries,Flatten[CC-Transpose[CC,{3,4,1,2}]],Flatten[Transpose[CC,1<->2]-CC]];

(*generate hom.lin.sys.of the symmetry transformation...*)
matrix=CoefficientArrays[{symmetries},Flatten[CC]][[2,1]];
(*... and check its rank*)
rank=MatrixRank[matrix]
(*show coefficient matrix*)
(*MatrixForm[matrix]*)
(*The no.of indep.comps.equals the no.of comps.*)
(*minus the no.of constraints (above matrix rank)*)
Print["CC has ",comps-rank," independent components. 
The constraints can be imposed as following in the componentes:"]
(*Give a symbolic solution*)
result=Solve[symmetries==0,Flatten[CC]]
Print["The resulting tensor is:"]
CC/.result[[1]]//MatrixForm