(* ------------------------------------- *)
(*                                       *)
(*  This file is called in for all of    *)
(*  the different related problems which *)
(*  are posed by Keith Glover.           *)
(*                                       *)
(* ------------------------------------- *)

(* -------------------- 

  Read in the basic files

   --------------------  *)

<<NCGB.m;
<<NCSR6.m;

(* -------------------- 

  Set Variables to be noncommutative

   --------------------  *)

Avariables = {A11,A12,A21,A22};
SetNonCommutative[Avariables];

(*        The following command is equivalent to:

        AvariablesTP = {tp[A11], tp[A12], tp[A21], tp[A22]};
*)
AvariablesTP = Map[tp,Avariables];

Bknowns = {B1,B2};
SetNonCommutative[Bknowns];
BknownsTP = Map[tp,Bknowns];

Bunknowns = {B12,B22,B32};
SetNonCommutative[Bunknowns];
BunknownsTP = Map[tp,Bunknowns];

Cknowns = {C1,C2}; 
SetNonCommutative[Cknowns];
CknownsTP = Map[tp,Cknowns];

Cunknowns = {C21,C22,C23}; 
SetNonCommutative[Cunknowns];
CunknownsTP = Map[tp,Cunknowns];

Dvariables = {D11,D12,D21,D22}; 
SetNonCommutative[Dvariables];
DvariablesTP = Map[tp,Dvariables];

Pvariables = {P11,P12,P13,P21,P22,P23,P31,P32,P33};
SetNonCommutative[Pvariables];
PvariablesTP = Map[tp,Pvariables];
PvariablesInv=Map[Inv,{P11,P22,P33}];
AllPiiVars={P11,Inv[P11],tp[P11],Inv[tp[P11]],
            P22,Inv[P22],tp[P22],Inv[tp[P22]],
            P33,Inv[P33],tp[P33],Inv[tp[P33]]};
AllOffPVars={P12,tp[P12],
             P13,tp[P13],
             P21,tp[P21],
             P23,tp[P23],
             P31,tp[P31],
             P32,tp[P32]};
PiiInvertible={P11**Inv[P11]-1,
               Inv[P11]**P11-1,
               P22**Inv[P22]-1,
               Inv[P22]**P22-1,
               P33**Inv[P33]-1,
               Inv[P33]**P33-1};


invPtil= {{N11,N12,N13},
          {N21,N22,N23},
          {N31,N32,N33}};
 
SetNonCommutative[N11,N12,N13,N21,N22,N23,N31,N32,N33];
AllNVars={N11,tp[N11],
          N22,tp[N22],
          N33,tp[N33],
          N12,tp[N12],
          N13,tp[N13],
          N21,tp[N21],
          N23,tp[N23],
          N31,tp[N31],
          N32,tp[N32]};

othervariables = {Sigma,Inv[Sigma]};
SetNonCommutative[othervariables];
othervariablesTP = Map[tp,othervariables];

(* -------------------- 

  Set sigma=1 to remove commutative component

   --------------------  *)

sigma = 1;

(* -------------------- 

  Set all of the tilda matrices

   --------------------  *)

matAtil = {
{ A11,    0,  A12/2},
{ 0  ,  A11, -A12/2},
{ A21, -A21, A22}
};

matBtil = {
{ B1, 0},
{ 0 , sigma Inv[Sigma]**tp[C1]},
{ B2, -tp[C2]}
};

matCtil = {
{ 0,  -2C1, C2},
{ -2 sigma tp[B1]**Inv[Sigma],  0, -tp[B2]}
};

matDtil = {
{ D11, D12},
{ D21, D22}
};

matPtil = {
{P11,P12,P13},
{P21,P22,P23},
{P31,P32,P33}
};

invPtil= {
{N11,N12,N13},
{N21,N22,N23},
{N31,N32,N33}
};

invPtiltp= Transpose[Map[tp,invPtil,{2}]];

invPtilSelfAdjoint = MapThread[Equal,
                      {Flatten[invPtil],Flatten[invPtiltp]}];

Id3 = IdentityMatrix[3];

eqinvP= Flatten[{ MM[invPtil,matPtil]- Id3, MM[matPtil,invPtil]- Id3}];
(* -------------------- 

  Conditions such that matPtil is of desired form
  as in Glover paper

   --------------------  *)

matPtilSelfAdjoint =  {
tp[P11]==P11, tp[P12]==P21, tp[P13]==P31,
tp[P21]==P12, tp[P22]==P22, tp[P23]==P32,
tp[P31]==P13, tp[P32]==P23, tp[P33]==P33};

matAtilPartC = {
{ A11, A12,   0},
{ A21, A22,   0},
{   0,   0, A11}};

matBtilPartC = {
{B1,B12},
{B2,B22},
{B1,B32}};


(* B12, B22, B32 unknowns *)

matCtilPartC = {
{ C1, C2, -C1},
{ C21, C22, C23}
};


matPtilDIAGONAL = { 
P12==0, P13==0, P21==0, P23==0, P31==0, P32==0,
tp[P12]==0, tp[P13]==0, tp[P21]==0, tp[P23]==0, tp[P31]==0, tp[P32]==0
};

matPtilSPECIFIC = {
P11==Sigma, P22 == sigma Inv[Sigma], P33 == 2 sigma,
tp[P11]==Sigma, tp[P22] == sigma Inv[Sigma], tp[P33] == 2 sigma
};

matDtilSPECIFIC = {
D11==0, D12==2 sigma, D21==2 sigma ,D22==0 ,
tp[D11]==0, tp[D12]==2 sigma, tp[D21]==2 sigma ,tp[D22]==0 
};

(* -------------------- 

  First equation from Glover result

   --------------------  *)

Equation1 = MatMult[matAtil,matPtil] + 
MatMult[matPtil,tpMat[matAtil]] + 
MatMult[matBtil,tpMat[matBtil]];

(* -------------------- 

  Second equation from Glover result

   --------------------  *)

Equation2 = MatMult[matPtil,tpMat[matCtil]] + 
MatMult[matBtil,tpMat[matDtil]]; 

(* --------------------

  Transpose of second equation

   --------------------  *)

Equation2Prime = MatMult[matCtil,matPtil] + 
MatMult[matDtil,tpMat[matBtil]]; 

(* -------------------- 

  Third equation from Glover result

   --------------------  *)

Equation3 = MatMult[matDtil,tpMat[matDtil]] -
                         4 sigma^2 IdentityMatrix[2];

(* -------------------- 

  Transpose of the third equation from 

   --------------------  *)

Equation3Prime = MatMult[tpMat[matDtil],matDtil] - 4 sigma^2 IdentityMatrix[2];

(* --------------------

  This equation is a consequence of equations 1 and 2.

   --------------------  *)

Equation4 = MatMult[invPtil,matAtil] + 
            MatMult[tpMat[matAtil],invPtil] + 
            MatMult[tpMat[matCtil],matCtil]/(4 sigma^2);


(* -------------------- 

  Conditions such that the matrices of equations
  1 2 and 3 are all zero.

   --------------------  *)

AllThreeEquations = Flatten[{Equation1,Equation2,Equation3}];
AllThreeEquations = Map[(#==0)&,AllThreeEquations];

AllSixEquations = Flatten[{Equation1,Equation2,Equation3,
                           Equation2Prime,Equation3Prime,Equation4}];
AllSixEquations = Map[(#==0)&,AllSixEquations];

(* -------------------- 

  Set the A, B and C matrices
  Set a 2 by 2 matrix based on Sigma

   --------------------  *)

matA = {{ A11,A12},{A21,A22}};

matSigma = {{ Sigma,0},{0,1}};

matB = {{B1},{B2}};

matC = {{C1,C2}};

(* -------------------- 

  Record the matrices in the alpha and beta equations

   --------------------  *)

alphaEquation = 
MM[matA,matSigma] + MM[matSigma,tpMat[matA]] + MM[matB,tpMat[matB]];

betaEquation = 
MM[tpMat[matA],matSigma] + MM[matSigma,matA] + MM[tpMat[matC],matC];

(* -------------------- 

  Record coniditons such that the above
  two matrices are zero

   --------------------  *)

greekEquations = Flatten[{alphaEquation,betaEquation}];

(* -------------------- 

  Record assumption made throughout

   --------------------  *)

ASSUMPTIONS = {
Sigma**Inv[Sigma]==1,
Inv[Sigma]**Sigma==1,
tp[Sigma]**tp[Inv[Sigma]]==1,
tp[Inv[Sigma]]**tp[Sigma]==1,
tp[Sigma]==Sigma,
tp[Inv[Sigma]]==Inv[Sigma]
};
