Nonzero Tensor Elements of (X^((2, 3, 4)))

This program allows to calculate the nonzero components of nonlinear susceptibility tensors X and their connections.
With "Initialization" one has to select the rank of the tensor, the symmetry operations, under whish the crystal is invariant and choose, if the calculation should also exploit the intrinsic permutation symmetry of the last n-1 indizes that ist given by n-th harmonic generation. Please be patient, as the computer does it's calculations after the input.
This section must be evaluated before using the later sections.
To caluculate the angular dependence the program has to know the (hkl)-direction that is perpendicular to the surface. Then SS, SP, PS and PP configurations are calculated.

Initialization
(start me first)

To include further symmetry operations, place '{"description",{matrix1,matrix2,..}},' within the brackets of  'operations'.

Clear["Global`*"]
Off[General::spell1]
Needs["LinearAlgebra`MatrixManipulation`"]

(******* Choose dimension *******)
dim=Input["Dimension of Tensor:\nType `n` for (n-1)th harmonics.\nE.g. `5` for 4th-harmonics"];

(******* Choose symmetry operations *******)
Module[{z,str,operations,dz,
        dx2={{ 1,0,0},{0,-1,0},{0,0,-1}},        (*3 rotational symmetries*)
        dy2={{-1,0,0},{0, 1,0},{0,0,-1}},
        xyz={{0,0, 1},{ 1,0,0},{0, 1,0}},
        xy     ={{0, 1,0},{ 1,0,0},{0,0, 1}},    (*4 mirror symmetries*)
        zminusz={{ 1,0,0},{0, 1,0},{0,0,-1}},
        xminusx={{-1,0,0},{0, 1,0},{0,0, 1}},
        yminusy={{ 1,0,0},{0,-1,0},{0,0, 1}}},
{dz[n_Integer]:={{Cos[2 Pi/n],-Sin[2 Pi/n],0},{Sin[2 Pi/n],Cos[2 Pi/n],0},{0,0,1}};
 (*Schoenflies not.: axis like in Wyckhoff, The Structure Of Crystals, 1931*)
 operations={
      {"c1 ",{ IdentityMatrix[3]}},           (* 1 *)       (*tricinic*)
      {"ci ",{-IdentityMatrix[3]}},           (* 1- *)
      {"cs ",{yminusy}},                      (* m *)       (*monoclinic*)
      {"c2 ",{dy2}},                          (* 2 *)
      {"c2h",{dy2,yminusy}},                  (* 2/m *) 
      {"c2v",{dz[2],yminusy}},                (* 2mm *)     (*orthorombic*)
      {"d2 ",{dz[2],dx2}},                    (* 222,v *)
      {"d2h",{dz[2],yminusy,dx2}},            (* mmm,vh *)
      {"s4 ",{-dz[4]}},                       (* 4- *)      (*tetragonal*)
      {"d2d",{dz[2],dx2,xy}},                 (* 4-2m,vd *)
      {"c4 ",{dz[4]}},                        (* 4 *)
      {"c4h",{dz[4],zminusz}},                (* 4/m *)
      {"c4v",{dz[4],xminusx}},                (* 4mm *)
      {"d4 ",{dz[4],dx2}},                    (* 422 *)
      {"d4h",{dz[4],dx2,zminusz}},            (* 4/mmm *)
      {"t  ",{xyz,dz[2]}},                    (* 23 *)      (*cubic*)
      {"th ",{xyz,dz[2],zminusz}},            (* m3 *)
      {"td ",{xyz,dz[2],xy}},                 (* 4-3m *)
      {"o  ",{xyz,dz[4]}},                    (* 432 *)
      {"oh ",{xyz,dz[4],zminusz}},            (* m3m *)
      {"c3 ",{ dz[3]}},                       (* 3 *)       (*trigonal*)
      {"c3i",{-dz[3]}},                       (* 3- *)
      {"c3v",{ dz[3],xminusx}},               (* 3m *)   (*mirror pl perp to x*)
      {"d3 ",{ dz[3],dx2}},                   (* 32 *)
      {"c3d",{-dz[3],dx2}},                   (* 3-m,d3d *)
      {"c3h",{dz[3].zminusz}},                (* 6- *)      (*hexagonal*)
      {"d3h",{dz[3].zminusz,dx2}},            (* 6-m2 *)
      {"c6 ",{dz[6]}},                        (* 6 *)
      {"c6h",{dz[6],zminusz}},                (* 6/m *)
      {"c6v",{dz[6],yminusy}},                (* 6mm *)
      {"d6 ",{dz[6],dx2}},                    (* 622 *)
      {"d6h",{dz[6],dx2,zminusz}},            (* 6/mmm *)
      {"isotropic",{dz[17],xminusx}}
};
 str="Choose symmetries:\n| ";
 For[z=1,z<=Length[operations],z=z+1,
     str=str <> ToString[z] <> " : " <> operations[[z]][[1]] <> " | "]; 
 symm=operations[[Input[str]]][[2]]//N;   (*Symmetry matrices for specific problem*)
}];

(******* intrinsic symmetry y/n *******)
createvector:=
Module[{x={},
        v=Table[3^(dim-i),{i,2,dim}],
        b={""},
        m,nulllist,newlist},
{Do[b=Flatten[Table[{"0","1","2"}[[i]]<>b[[j]],{i,3},{j,Length[b]}]],{dim-1}];
 Do[b[[i]]=ToExpression[Sort[Characters[b[[i]]]]].v,{i,1,Length[b]}];
 m=Max[b]+1;
 b=m-b;
 nulllist=Table[0,{i,3^(dim-1)}];
 While[m>0,
 {newlist=ReplacePart[nulllist,1,Position[b,m]];
  b=b-m*newlist;
  x=AppendTo[x,newlist];
  m=Max[b]}];
 BlockMatrix[Outer[Times,IdentityMatrix[3],x]]
}][[1]]

If[Input["To include intrinsic symmetry for n-th harmonics, type '1'."]==1,
   {x=Transpose[createvector];},
   x=IdentityMatrix[3^dim];];


(******* Subroutines for main program *******)

vectortomatrix[vec_]:=If[VectorQ[vec],{vec},vec]

lines[m_List]:=Dimensions[m][[1]]

rows[m_List] :=Dimensions[m][[2]]

eulerbackward[mat_List]:=
Module[{len=lines[mat],m=mat,pos},
{While[len>0,
 {pos=Flatten[Position[Sign[Abs[m[[len]]]],1]][[1]];
  If[pos>0,
     {m[[len]]=Chop[m[[len]]/m[[len,pos]]],
      Do[m[[i]]=Chop[m[[i]]-m[[i,pos]]*m[[len]]],{i,1,len-1}]}],
  len=len-1}];
 m=vectortomatrix[m]}][[1]]


pickcolumn[mat_List,vec_List]:=
Block[{m2={}},
{Do[m2=Join[m2,{Transpose[mat][[vec[[i]]]]}],{i,Length[vec]}];
 m2=vectortomatrix[m2];
 m2=Transpose[m2]}][[1]]


findcolumn[mat_List]:=
Module[{m,j=1,liste={}},
{m=Join[mat,{Table[0,{za,rows[mat]}]}];
 Do[If[m[[j,i]]==0,
       liste=Join[liste,{i}],
       j=j+1],
    {i,rows[mat]}];
 liste}][[1]]


sortrow[mat_List,col_Integer]:=
Block[{m,v,p},
{m=vectortomatrix[mat];
 v=Abs[Transpose[m][[col]]];
 p=Flatten[Position[v,Max[v]]][[1]];
 v=m[[1]];
 m[[1]]=m[[p]];
 m[[p]]=v;
 m}][[1]]


eulerforward[mat_List]:=
Block[{m,pos,mrows=rows[mat],b={}},
{m=Chop[Table[mat[[i]]/Max[Abs[mat[[i]]],1],{i,Length[mat]}]];
 m=Join[m,{Table[0,{i,mrows}]}];
 pos=1;
 While[pos<=mrows,
      {m=sortrow[m,pos];
       If[m[[1,pos]]!=0,
          {m[[1]]=Chop[m[[1]]/m[[1,pos]]],
           Do[m[[i]]=Chop[m[[i]]-m[[i,pos]]*m[[1]]],{i,2,Length[m]}],
           AppendTo[b,m[[1]]],
           m=Delete[m,1]}];
        pos=pos+1}];
 vectortomatrix[b]}][[1]]


findnullspace[mat_]:=
Block[{m,m2,lst,len,zauber},
{m=vectortomatrix[mat];
 If[Max[Chop[Abs[m]]]==0,
    m2=IdentityMatrix[rows[m]],
    {len=rows[m];
     m=eulerforward[m];
     m=eulerbackward[m];
     lst=findcolumn[m];
     If[lst=={},
        m2=Table[{0},{i,len}],
        {
         m=pickcolumn[m,lst];
         zauber=Complement[Table[i,{i,len}],lst];
         m2=Table[0,{j,len},{i,Max[Length[lst],1]}];
         Do[Do[m2[[zauber[[j]],i]]=-m[[j,i]],{i,rows[m]}],{j,lines[m]}];
         Do[m2[[lst[[i]],i]]=1,{i,Length[lst]}];
        }
       ];
    };
 ];
 m2}][[1]]
 
 
 
(****** Main program: determination of relationship between tensor elements *****)
For[z=1,z<=Length[symm],z=z+1,
    {mt=Transpose[symm[[z]]],

     m =Inverse[symm[[z]]],
     Do[m=BlockMatrix[Outer[Times,m,mt]],{dim-1}],
     m=m.x-x,
     inv=findnullspace[m];
     x=Chop[x.inv]}]

x=Rationalize[x,10^(-10)];
x=Transpose[x];
If[Max[x]!=0,{x=eulerforward[x],x=eulerbackward[x]}];

chi=Transpose[x].Table["X"<>ToString[i],{i,Length[x]}];
Do[chi=Partition[chi,3],{dim}]

Nonzero elements and their relationship 
(run Initialization first)

readout={""};
Do[readout=Flatten[Table[{"x","y","z"}[[i]]<>readout[[j]],{i,3},
                         {j,Length[readout]}]],{dim}]

(*Non zero elements*)
   y=Apply[Plus,Sign[Abs[x]]];
   output1=readout[[Flatten[Position[Sign[y],1]]]];
   Print[Length[Position[Sign[y],1]]," non zero tensor elements: (dependend on ",
         Length[x]," parameters)"]
   Print[output1]
   Print[""]

(*Relationship between non zero elements*)
   output2=" ";
   If[Apply[Plus,Chop[y^2-y]]!=0,
      {Print["Tensor elements are not separable."],
       Print["'readout' = Flatten[chi]"]},
      {Do[{maximum=Max[x[[i]]],
           Do[If[x[[i,j]]!=0,
                 output2=output2<>ToString[readout[[j]]*maximum/x[[i,j]]]<>" = "],
              {j,3^dim}],
            output2=output2 <> ToString[maximum*X] <> ToString[i]<>" \n "},
          {i,Length[x]}],
       Print[Length[x]"Independent tensor elements:"],
       Print[output2]}];

Angular dependence of the signal in reflection
(run Initialization first)

a=IdentityMatrix[3];

a[[1,3]]=Input["Insert (h)"];
a[[2,3]]=Input["Insert (k)"];
a[[3,3]]=Input["Insert (l)"];
a=a/Sqrt[a[[1,3]]^2+a[[2,3]]^2+a[[3,3]]^2];

If[N[a[[2,3]]]==0.,
   a[[2,2]]=1,
   a[[2,2]]=Sqrt[ a[[3,3]]^2 / (a[[3,3]]^2+a[[2,3]]^2)]];

If[N[a[[3,3]]]==0.,
   a[[3,2]]=-Sign[a[[2,3]]],
   a[[3,2]]=-Sign[a[[2,3]] a[[3,3]]] Sqrt[1-a[[2,2]]^2]];

a[[1,1]]= a[[2,2]] a[[3,3]]-a[[3,2]] a[[2,3]];
a[[2,1]]= a[[3,2]] a[[1,3]];
a[[3,1]]=-a[[2,2]] a[[1,3]];
a=Chop[a];

picture=Show[Graphics3D[{Polygon[{{-0.3,-0.3,0},{-0.3,0.3,0},{0.3,0.3,0},{0.3,-0.3,0}}],
                         Line[{{0,0,0},a[[1]]}],Text["x",0.9 a[[1]]],
                         Line[{{0,0,0},a[[2]]}],Text["y",0.9 a[[2]]],
                         Line[{{0,0,0},a[[3]]}],Text["z",0.9 a[[3]]],

                         Line[{{0,0,0},{1,0,0}}],Text["e1",{1.1,0,0}],
                         Line[{{0,0,0},{0,1,0}}],Text["e2",{0,1.1,0}],
                         Line[{{0,0,0},{0,0,1}}],Text["e3",{0,0,1.1}]},
                        Boxed->False]];

(* delete linear fresnel factors *)
Clear[lfx,lfy,lfz];

Print["Angular dependence of signal:"]

Print["\nSS configuration:"]
ss=chi;
Do[ss =ss.(a.{-Sin[α] lfy,Cos[α] lfy,0}),{dim-1}];
ss=(ss.(a.{-Sin[α] fy,Cos[α] fy,0}))[[1]]//Expand//Chop//TrigReduce//Factor;
Print[ss]

Print["\nSP configuration:"]
spx=chi;
Do[spx=spx.(a.{-Sin[α] lfy,Cos[α] lfy,0}),{dim-1}];
spx=(spx.(a.{Cos[α] fx, Sin[α] fx, 0}))[[1]]//Expand//Chop//TrigReduce//Factor;
spz=chi;
Do[spz=spz.(a.{-Sin[α] lfy,Cos[α] lfy,0}),{dim-1}];
spz=(spz.(a.{0,0,fz}))[[1]]//Expand//Chop//TrigReduce//Factor;
Print[spx," +\n",spz]

Print["\nPS configuration:"]
ps=chi;
Do[ps =ps.(a.{Cos[α] lfx,Sin[α] lfx,lfz}),{dim-1}];
ps=(ps.(a.{-Sin[α] fy,Cos[α] fy,0}))[[1]]//Expand//Chop//TrigReduce//Factor;
Print[ps]

Print["\nPP configuration:"]
ppx=chi;
Do[ppx=ppx.(a.{Cos[α] lfx,Sin[α] lfx,lfz}),{dim-1}];
ppx=(ppx.(a.{Cos[α] fx, Sin[α] fx, 0}))[[1]]//Expand//Chop//TrigReduce//Factor;
ppz=chi;
Do[ppz=ppz.(a.{Cos[α] lfx,Sin[α] lfx,lfz}),{dim-1}];
ppz=(ppz.(a.{0,0,fz}))[[1]]//Expand//Chop//TrigReduce//Factor;
Print[ppx," +\n",ppz]

Graphs
(run Initialization and Angular dependence first)

(* refracted index of fundamental and harmonic signal *)
ϵfundamental=Sqrt[3.692];
ϵharmonic   =Sqrt[1.831];

(* linear fresnel factors *)
lfy=2 Cos[θ]/(Cos[θ]+Sqrt[ϵfundamental-Sin[θ]^2]);
lfx=2 Cos[θ]/(ϵfundamental Cos[θ]+Sqrt[ϵfundamental-Sin[θ]^2]);
lfz=2 Sqrt[ϵfundamental-Sin[θ]^2]/(ϵfundamental Cos[θ]+Sqrt[ϵfundamental-Sin[θ]^2]);

(* nonlinear fresnel factors *)
fresnels [e_,e4_,t_]=
   4 Pi/(e-e4) (Sqrt[e4-Sin[t]^2]-Sqrt[e-Sin[t]^2])/(Sqrt[e4-Sin[t]^2]+Cos[t]);
fresnelp [e_,e4_,t_]=4 Pi /(Sqrt[e-Sin[t]^2]+Sqrt[e4-Sin[t]^2])/
                    (Sqrt[e4-Sin[t]^2]+e4 Cos[t]);
fresnelpx[e_,e4_,t_]=fresnelp[e,e4,t] Sqrt[e4-Sin[t]^2];
fresnelpz[e_,e4_,t_]=fresnelp[e,e4,t] Sin[t];
fy=fresnels [ϵfundamental,ϵharmonic,θ];
fx=fresnelpx[ϵfundamental,ϵharmonic,θ];
fz=fresnelpz[ϵfundamental,ϵharmonic,θ];

(* various plots *)
Plot[ss^2       /.{θ->Pi/4.0,"X1"->5,"X2"->1},{α,0,2 Pi}];
Plot[(spx+spz)^2/.{θ->Pi/4.0,"X1"->5,"X2"->1},{α,0,2 Pi}];
Plot[ps^2       /.{θ->Pi/4.0,"X1"->5,"X2"->1},{α,0,2 Pi}];
Plot[(ppx+ppz)^2/.{θ->Pi/4.0,"X1"->5,"X2"->1},{α,0,2 Pi}];

Plot3D[ss^2/.{"X1"->5,"X2"->1}
       ,{α,0,2 Pi},{θ,0,Pi/2},PlotPoints->50];
Plot3D[(spx+spz)^2/.{"X1"->5,"X2"->1}
       ,{α,0,2 Pi},{θ,0,Pi/2},PlotPoints->50];
Plot3D[ps^2/.{"X1"->5,"X2"->1}
      ,{α,0,2 Pi},{θ,0,Pi/2},PlotPoints->50];
Plot3D[(ppx+ppz)^2/.{"X1"->5,"X2"->1}
       ,{α,0,2 Pi},{θ,0,Pi/2},PlotPoints->50];