(* 
 *  TRproof.in
 *
 *  This Mathematica script shows that if kappa is a (2,2)-tensor, then 
 * 
 *       (det kappa)^2 TR(kappa) + TR(adjugate(kappa)) = 0
 *
 *  where TR(kappa) is the Tamm-Rubilar tensor densities for kappa 
 *  and adjugate(kappa) is the adjugate of kappa.
 *
 *  To run script: type:
 *
 *     nohup math -noprompt -run "<<TRproof.in" > TRproof.out 2> TRproof.err &
 *
 *)

scriptStartTimer = TimeUsed[];
<<../kappaLib.m

(*
 *  Define most general (2,2)-tensor kappa using four 3x3 matrices
 *)
Ax = emGeneral3x3["a"];
Bx = emGeneral3x3["b"];
Cx = emGeneral3x3["c"];
Dx = emGeneral3x3["d"];
kappa = emABCDToKappa[Ax, Bx, Cx, Dx];


(* 
 *  The next command computes the TR-tensor of kappa, and Expand simplifies the
 *  expression by canceling/combining common terms. However, it does not try to 
 *  factor the expression. 
 *)
Print["Computing Tamm-Rubilar tensor for kappa..."];
TRkappa = Expand[emKappaToTR[kappa]];  


(* 
 *  Compute Det[kappa]^2
 *
 *  Notes: 
 *    - Note that the power is not expanded yet.  
 *    - Adding Expand[..] here will slow down to program in the main loop.
 *)
kappaDetSq = emDet[kappa]^2;


(*
 * Next compute the Tamm-Rubilar tensor density of Adj(kappa). That is, we  
 * compute the composite function 
 *
 *  TR o Adj (kappa).
 *
 * To simplify the end result we first expand the expression for 
 * TR operation before substituting Adj(kappa)
 *
 *)
Print["Computing Tamm-Rubilar tensor for adj(Kappa)..."];
coEntry[i_, j_] := ToExpression["co" <> ToString[i] <> ToString[j]];
MatAdj = Table[coEntry[i,j], {i, 1, 6}, {j, 1, 6}];
TRkappaAdj = Expand[emKappaToTR[emMatrixToKappa[MatAdj]]]; 
MatAdjX = emKappaToMatrix[emAdjugate[kappa]];
Sub = Flatten[
           Table[coEntry[i, j] -> MatAdjX[[i]][[j]], 
           {i, 1, 6}, {j, 1, 6}]
      ];

TRkappaAdj = TRkappaAdj /. Sub;


SimplifyRecursive[f_]:= Module[
   {ll, simp, t0, var, i, bound, fList},

   (* 
    *  If argument is zero, do nothing. 
    *)
   If[f==0, Return[]];

   ll = emSize[f];

   If[ Length[Variables[f]] < 27,

      (* 
       *   Size of f is "small". Try to simplify 
       *)

	t0 = TimeUsed[];	
	simp = Expand[f];
	If[simp==0, 
	     Print["Verified formula of size "<>ToString[ll]<>" in " <>ToString[TimeUsed[]-t0] <> " CPU seconds"];
	     res = 0;
	     ,
	     Print["Fail: (a)"];
	     Quit[];
	     ,
	     Print["Fail: (b)"];
	     Quit[];
        ];
   
   ,

     (* 
      *   Size of f is "large". Eliminate one variable. 
      *) 
     var = Variables[f];
     var = var[[1]];
     fList = {f};

     bound = 5;     
     For[i=1, i<= bound, i++,
     	fList = Append[fList, D[fList[[Length[fList]]], var] ];
     ];

     If[fList[[Length[fList]]] ==0,
         (* true. do nothing *)
	 ,
         (* false: abort *)
	 Print["Fail: bound too small (a)."];
	 Quit[];
	 ,
	 (* otherwise: abort *)
	 Print["Fail: bound too small (b)."];
	 Quit[];
     ];

     fList = fList /. var->0;

     Print["Did Taylor expansion on variable "<>ToString[var]
            <> ".  Spawing " <> ToString[Length[fList]] <> " new calls..."];

     For[i=1, i<= Length[fList], i++,
     	SimplifyRecursive[ fList[[i]] ];
     ];

   ];

]


Print["Making list of polynomials to prove..."];
polys = {};
totSize = 0;
For[i = 1, i <= 4, i++,
 	For[j = i, j <= 4, j++,
   		For[k = j, k <= 4, k++,
     			For[l = k, l <= 4, l++,

				tmpA = emReadTR[TRkappaAdj, i, j, k, l]; 	
 				tmpB = kappaDetSq emReadTR[TRkappa, i, j, k, l];

				tmp = tmpA+tmpB;

				Print[""];
			       	Print[" ***  Length of polynomial at index " 
					<> ToString[i]
		     	     	        <> ToString[j]
 				        <> ToString[k]
				        <> ToString[l]
				        <> ": "
					<> ToString[emSize[tmp]]
				];

				Print["Simplifying ..."];
				time0 = TimeUsed[];
				res = SimplifyRecursive[tmp];
				time1 = TimeUsed[];
				Print[""];
				Print["Done in computation time : " <> ToString[time1-time0] <> "CPU seconds"];
				Print[""];

				totSize += emSize[tmp];

	     		];
     		];
   	];
];

Print[""];
Print["Total length of polynomials "<>ToString[totSize]];


Print["Main loop done. Total CPU time:"<> ToString[N[TimeUsed[]-scriptStartTimer, 2]]<> "s. Exiting."];

Exit[];
