(* ::Package:: *)

Print["KappaLib v1.1"];

(* 
 * KappaLib is a Mathematica package for manipulating (2,2) tensors that 
 * arise in the spacetime formulation of electromagnetics.
 *
 * Written by Matias Dahl, 2010-2011.
 *
 * Note: all routines in KappaLib start with 'em'.
 *
 * Feedback and comments are welcome: matias.dahl@tkk.fi 
 *
 * Changelog:
 *  1.0 -- Initial version
 * 
 *  1.1 -- Changed the matrix representation of a (2,2)-tensor
 *
 *)


(*
 *  Section 1:
 *  - - - - - - Routines for encoding an antisymmetric (2,2)-tensor - - - - - -
 *
 *  In local coordinates, an antisymmetric (2,2)-tensor is described as an array with
 *  4 indices kappa^ij_ab whence
 *
 *    kappa = 1/2 kappa^ij_ab d/dx^i otimes d/dx^j otimes dx^a otimes dx^b.
 *
 *  By antisymmetry, this array contains redundant information. 
 *  We therefore only store entries in the array when i<j and a<b. 
 *
 *  Suppose kap^ij_ab are components for a (2,2)-tensor kappa. In KappaLib, 
 *  it is encoded in a 4-dimensional array kappa[[i]][[j]][[a]][[b]] whose 
 *  entries are defined as.
 * 
 *    kappa[[i]][[j]][[a]][[b]] =    kap^ij_ab,     when i<j and a<b,
 *                              =    0,               otherwise.
 *
 *  In the above, a,b,i,j take values in 1,2,3,4.
 *
 *)


(*
 *  If kappa is an array that encodes a (2,2)-tensor kap with components kap^ij_ab, then
 *  emReadNormal[[i]][[j]][[a]][[b]] returns kap^ij_ab.
 *)
emReadNormal[kappa_, i_, j_ ,k_ ,l_] := Module[
    {tmpV1, tmpV2, tmpV3, tmpV4},

	tmpV1 = Min[{i,j}]; 
    tmpV2 = Max[{i,j}]; 
    tmpV3 = Min[{k,l}];
    tmpV4 = Max[{k,l}];
                    
    Signature[{i,j}] Signature[{k,l}] kappa[[tmpV1]][[tmpV2]][[tmpV3]][[tmpV4]]
                 
  ]


(*
* mediumMakeNormal[res_]
*
* Set res^ab_cd = 0 when a>=b or c>=d.
*)
emMakeNormal[res_] := Module[
	{aa, bb, cc, dd, tmp},
	
	(* make a copy that we modify and then return *)
	tmp = res;

	For[aa=1, aa<=4, aa++,
		For[bb=1, bb<=4, bb++,
			For[cc=1, cc<=4, cc++,
				For[dd=1, dd<=4, dd++,
					If[ aa>= bb || cc>=dd, 
						tmp[[aa]][[bb]][[cc]][[dd]] = 0;
					];
				];
			];
		];
	];

	tmp

]


(*
 *  Compute how many bytes are required to store an mathematical object 
 *  when converted into text form?
 *)
emSize[obj_] := StringLength[ToString[InputForm[obj]]]


(*
 * Section 2:
 *  - - - - - - Specific (2,2)-tensors - - - - - - - 
 *
 *)


(*
 *   Return an array that represents the zero (2,2)-tensor.
 *)
emZeroKappa[] := Table[0, {i,1,4},{j,1,4},{k,1,4},{l,1,4}]


(*
 *   Return (2,2)-tensor with random components
 *)
emRandomKappa[] := Module[
	{res},

	res = Table[(Random[]-1/2)*2,
                   {i,1,4},{j,1,4},{k,1,4},{l,1,4}
    ];

	emMakeNormal[res]

]


(*
 *   Return the identity (2,2)-tensor that maps a 2-form into itself
 *)
emIdentityKappa[] := Module[
	{res, i, j, k, l},

	res = Table[
		KroneckerDelta[i,k] KroneckerDelta[j,l]-KroneckerDelta[i,l] KroneckerDelta[j,k],
        {i,1,4}, {j,1,4}, {k,1,4}, {l,1,4}
    ];

	emMakeNormal[res]

]


(*
 *   For a 4x4 symmetric non-degenerate matrix 'metric', return a (2,2)-tensor that
 *   represents the Hodge operator of 'metric' restricted to 2-forms.
 * 
 *)
emHodge[metric_] := Module[
	{resKappa, invMetric, detSq, i,j ,m,n, a,b},

	(* compute inverse and determinant only once *)
	resKappa = emZeroKappa[];
	invMetric = Inverse[metric];
	detSq = Sqrt[Abs[Det[metric]]];
	
	For[i=1, i<=4, i++,
		For[j=i+1, j<=4, j++,
			For[m=1, m<=4, m++,
				For[n=m+1, n<=4, n++,
					resKappa[[i]][[j]][[m]][[n]] = detSq Sum[ 
							invMetric[[i]][[a]] invMetric[[j]][[b]] Signature[{a,b,m,n}],
								{a,1,4},{b,1,4}
							];
				];
			];
		];
	];

	resKappa

]


(*
 *  Section 3:
 *  - - - - - - Basic linear algebra for (2,2)-tensors - - - - - - - 
 *
 *  Any (2,2)-tensor can be treated as a linear map from 2-forms into 2-forms. 
 *  We may therefore represent (2,2)-tensors as 6x6 matrices, compute their 
 *  trace, determinant, composition, adjugate and inverse. 
 *
 *)


(*
 *  Indices for basis vectors for Omega^2(N) 
 *)
emOrd[i_]:={{1, 2}, {1, 3}, {1, 4}, {3, 4}, {4, 2}, {2, 3}}[[i]]


(* convert kappa into a 6x6 matrix with respect to basis:

    { dx^emOrd[i],    i=1,...,6 }

*)
emKappaToMatrix[kappa_] := Module[
  	{matRep, tMat, i, j},
  
  	(* allocate memory for 6x6 matrix *)
  	matRep = Table[0, {i, 1, 6}, {j, 1, 6}];
  	
  	For[i = 1, i <= 6, i++,
   		For[j = 1, j <= 6, j++,

     			matRep[[j]][[i]] = emReadNormal[kappa, 

        								emOrd[i][[1]], emOrd[i][[2]],
        								emOrd[j][[1]], emOrd[j][[2]]

									];
     		];
   	];
  
  	matRep
  
  ]


(* inverse of emlKappaToMatrix *)
emMatrixToKappa[matrix_] := Module[
  	{pList, matRep, i, j, sign, a, b, c, d},
  
  	(* allocate memory for (2,2)-tensor *)
  	res = emZeroKappa[];
  	
  	For[i = 1, i <= 6, i++,
   		For[j = 1, j <= 6, j++,

				a = Sort[emOrd[i]][[1]];
				b = Sort[emOrd[i]][[2]];
				c = Sort[emOrd[j]][[1]];
				d = Sort[emOrd[j]][[2]];

				sign = Signature[emOrd[i]] Signature[emOrd[j]];

				res[[a]][[b]][[c]][[d]] = sign matrix[[j]][[i]];

     		];
	];

(*  
	res = emMakeNormal[res]; (* this should not be necessary *)
*)

	res

]


(*
 *  Return the composition of two (2,2)-tensors.
 *)
emCompose[kappa1_, kappa2_] := Module[
	{mat1, mat2},

	mat1 = emKappaToMatrix[kappa1];
	mat2 = emKappaToMatrix[kappa2];

	emMatrixToKappa[mat1 . mat2]	

]


(* 
 *  Return the determinant of a (2,2)-tensor
 *)
emDet[kappa_] := Det[emKappaToMatrix[kappa]]


(*  
 *   Compute adjugate of kappa. For invertible linear maps, the adjugate is defined as
 *
 *      Adj A = Det (A) . Inverse(A)
 *
 *)
emAdjugate[kappa_] := Module[
	{i, j, mat},

	mat = emKappaToMatrix[kappa];

	mat = Table[
			(-1)^(i+j)Det[Drop[mat,{i,i},{j,j}]],
			{j,1,6},{i,1,6}
			];

	emMatrixToKappa[mat]
]


(*
 * compute trace of kappa 
 *
 * alternative: 1/2 Sum[emlReadNormal[medium, i,j,i,j], {i,1,4},{j,1,4}]
 *
 *)
emTrace[kappa_] := Tr[emKappaToMatrix[kappa]]


(* compute inverse of (2,2) tensor kappa *)
emInverse[kappa_] := emMatrixToKappa[Inverse[emKappaToMatrix[kappa]]]


(*
 * Section 5
 * - - - - - - - -  Routines for ABCD-representation of a (2,2)-tensor - - - - - - - - - - -
 *
 *  See Section 2.5 in the paper.
 *)


(*
 *  Return A,B,C,D matrices for a (2,2)-tensor.
 *
 *)
emKappaToABCD[medium_]:=Module[
	{resA, resB, resC, resD, ii, rr, a, b, m, n},

	(* create empty return variables *)
	resA = Table[0, {i,1,3}, {j,1,3}];
	resB = Table[0, {i,1,3}, {j,1,3}];
	resC = Table[0, {i,1,3}, {j,1,3}];
	resD = Table[0, {i,1,3}, {j,1,3}];

	(* create A-matrix *)
	For[ii=1, ii<=3, ii++,
		For[rr=1, rr<=3, rr++,
			resA[[rr]][[ii]] = -1/2 Sum[ 
				Signature[{ii,a,b}] emReadNormal[medium, rr+1, 1, a+1, b+1],
				{a,1,3},{b,1,3}
			];
		];	
	];
    
	(* create B-matrix *)
	For[ii=1, ii<=3, ii++,
		For[rr=1, rr<=3, rr++,
			resB[[rr]][[ii]] = -1/2 Sum[ 
				Signature[{rr,a,b}] emReadNormal[medium,a+1,b+1,ii+1,1],
				{a,1,3},{b,1,3}
			];
		];	
	];

	(* create C-matrix *)
	For[ii=1, ii<=3, ii++,
		For[rr=1, rr<=3, rr++,
			resC[[rr]][[ii]] = emReadNormal[medium, rr+1,1, ii+1,1];
		];	
	];

	(* create D-matrix *)
	For[ii=1, ii<=3, ii++,
		For[rr=1, rr<=3, rr++,
			resD[[rr]][[ii]] = 1/4 Sum[ 
				Signature[{rr,m,n}] Signature[{ii,a,b}] emReadNormal[medium,m+1,n+1,a+1,b+1],
				{a,1,3},{b,1,3},{m,1,3},{n,1,3}
			];
		];	
	];

	{resA, resB, resC, resD}
]


emABCDToKappa[matA_, matB_, matC_, matD_] := Module[
	{res, ii, jj, kk, ll, rr, ss, pp, qq},

	(* create empty return variable *)
    res = emZeroKappa[];

    (* Note: 
     *   We only need to fill in res[[i]][[j]][[k]][[l]] for which i<j and k<l. 
     *)

	(* import C-matrix *)
	For[ii=1, ii<=3, ii++,
		For[rr=1, rr<=3, rr++,
			res[[1]][[rr+1]][[1]][[ii+1]] = matC[[rr]][[ii]];
		];	
	];
    
	(* import B-matrix *)
	For[ii=1, ii<=3, ii++,
		For[jj=ii+1, jj<=3, jj++,
			For[rr=1, rr<=3, rr++,
				res[[ii+1]][[jj+1]][[1]][[rr+1]] = Sum[ 
					Signature[{kk,ii,jj}] matB[[kk]][[rr]],
					{kk, 1, 3}
				];
			];	
		];
	];
    
	(* import A-matrix *)
	For[ii=1, ii<=3, ii++,
		For[rr=1, rr<=3, rr++,
			For[ss=rr+1, ss<=3, ss++,
				res[[1]][[ii+1]][[rr+1]][[ss+1]] = Sum[ 
					Signature[{kk,rr,ss}] matA[[ii]][[kk]],
					{kk, 1, 3}
				];
			];
		];	
	];

	(* import D-matrix *)
	For[rr=1, rr<=3, rr++,
		For[ss=rr+1, ss<=3, ss++,
			For[ii=1, ii<=3, ii++,
				For[jj=ii+1, jj<=3, jj++,
					res[[ii+1]][[jj+1]][[rr+1]][[ss+1]]= Sum[ 
						Signature[{kk,rr,ss}] Signature[{ll,ii,jj}] matD[[ll]][[kk]],
						{kk, 1, 3}, {ll, 1, 3}
					];
				];	
			];
		];	
	];

	res

]


(*
 * - - - - - - -  General objects - - - - - 
 *  
 * When proving theorems it is useful to start with a general object.
 * The below routines help to create such objects.
 *)


(*
 * Return a (2,2)-tensor that depends on A,B,C,D metrices with arbitrary entries
 *)
emGeneralABCD[aName_, bName_, cName_, dName_]:=Module[

	{Ax, Bx, Cx, Dx},

	Ax = Table[ToExpression[aName<>ToString[i]<>ToString[j]],{i,1,3},{j,1,3}];
	Bx = Table[ToExpression[bName<>ToString[i]<>ToString[j]],{i,1,3},{j,1,3}];
	Cx = Table[ToExpression[cName<>ToString[i]<>ToString[j]],{i,1,3},{j,1,3}];
	Dx = Table[ToExpression[dName<>ToString[i]<>ToString[j]],{i,1,3},{j,1,3}];

	emABCDToKappa[Ax,Bx,Cx,Dx]

]


(*
 * Construct a general (2,2)-tensor that depends a 6x6 matrix with entries 
 *
 *    {varName_ij  : i,j = 1,..., 6} 
 *
 *)
emGeneralKappa[varName_] := Module[

	{mat},

	mat = Table[ToExpression[varName<>ToString[i]<>ToString[j]], {i,1,6}, {j,1,6}];
	
	emMatrixToKappa[mat]

]


(*
 *  return a general 3x3 matrix with entries
 * 
 *    {{v11,v12,v13},{v12,v22,v23},{v13,v23,v33}} 
 *)
emGeneralSymmetric3x3[v_] := Table[

        ToExpression[v<>ToString[Min[{i,j}]]
		    				<>ToString[Max[{i,j}]]],
         {i,1,3},{j,1,3}
]


(*
 *  return a general symmetric 3x3 matrix with entries
 * 
 *   {{v11,v12,v13},{v21,v22,v23},{v31,v32,v33}}
 *)
emGeneral3x3[varName_] := Table[

        ToExpression[varName<>ToString[i]
		    				<>ToString[j]],
         {i,1,3},{j,1,3}
]


(*
* - - - - - -- - - - - - Routines for computing the Tamm-Rubilar tensor density
*)


(* fast implementation of mediumPreTammRubilar[kappa_, i_,j_,k_,l_] *)
emPreTammRubilar[kappa_, i_,j_,k_,l_] := Module[

	{pList, pSig, NN, 
	AA, BB, CC, DD,
	a1, a2, a3, a4, 
	b1, b2, b3, b4, b5, b6,
	xx, yy,
	res},

	(* precompute permutations of 1,2,3,4 and their signatures *)
	pList = Permutations[{1,2,3,4}];
	NN = Length[pList];
    pSig = Table[
			Signature[{pList[[AA]][[1]], pList[[AA]][[2]], pList[[AA]][[3]], pList[[AA]][[4]]}],
			{AA, 1, NN}
		];

	(* init return value *)
	res = 0;
	
	For[AA=1, AA <= NN, AA++,
		{b1,b2,b5,xx} = pList[[AA]];

		If[xx == k && b1<b2,

			For[BB=1, BB <= NN, BB++,

				{b3,b4,b6,yy} = pList[[BB]];

				If[yy == l && b3<b4,

					For[CC=1, CC <= NN, CC++,
						{a1,a2,a3,a4} = pList[[CC]];

						If[a1<a2,

						res += 2^3 /(48)  pSig[[AA]] pSig[[BB]] pSig[[CC]] emReadNormal[kappa, a1, a2, b1, b2] emReadNormal[kappa, a3, i, b3, b4] emReadNormal[kappa, a4, j, b5, b6];

						];
						
					];
				];			
			];
		];
	];

	res

]


emKappaToTR[kappa_]:=Module[

	{preTR, TR, i, j, k, l, loop, perms},

	preTR = Table[
				emPreTammRubilar[kappa, i, j, k, l],
				{i,1,4},{j,1,4},{k,1,4},{l,1,4}];

	(* symmetrize *)
	TR = Table[0, {i,1,4}, {j, 1, 4}, {k,1,4}, {l,1,4}];

	(* Note: We only compute TR[[i]][[j]][[k]][[l]] for which i<=j<=k<=l. *)

	For[i=1, i<=4, i++,
		For[j=i, j<=4, j++,
			For[k=j, k<=4, k++,
				For[l=k, l<=4, l++,
					perms = Permutations[{i,j,k,l}];
					(* note: Permutations[{1,1,1,3}] will have 4 entries,
                     * and Permutations[{1,1,1,1}] will have only 1. *)
						
					TR[[i]][[j]][[k]][[l]] = Sum[
						preTR[[perms[[loop]][[1]]]][[perms[[loop]][[2]]]][[perms[[loop]][[3]]]][[perms[[loop]][[4]]]]
						, {loop, 1, Length[perms]}
						]/Length[perms];
				];
			];
		];
	];

	TR

]


emReadTR[tr_, i_, j_, k_, l_] := Module[

	{si, sj, sk, sl},

	{si, sj, sk, sl} = Sort[{i,j,k,l}];

	tr[[si]][[sj]][[sk]][[sl]]

]


(*
*
*  - - - - - - - Routines for Fresnel equation
*
*)


emKappaToFresnel[kappa_, xi_] := Module[

	{TR},

	TR = emKappaToTR[kappa];

	emTRToFresnel[TR, xi]
]


emTRToFresnel[TRtensor_, xi_] := Module[
	{i,j,k,l},

	Sum[

		emReadTR[TRtensor, i, j, k, l] xi[[i]] xi[[j]] xi[[k]] xi[[l]],

		{i,1,4}, {j,1,4}, {k,1,4}, {l,1,4}
	]

]


emFresnelToTR[fresnel_, xi_] := Module[
	{i,j,k,l, res},

   res = Table[0, {i,1,4},{j,1,4},{k,1,4},{l,1,4}];

	For[i=1, i<=4, i++,
		For[j=i, j<=4, j++,
			For[k=j, k<=4, k++,
				For[l=k, l<=4, l++,
                                    res[[i]][[j]][[k]][[l]] = 1/4! 
                                             D[D[D[D[fresnel, xi[[i]]],xi[[j]]],xi[[k]]],xi[[l]]]
                                    ];
                                ];
                          ];
       ];

       res
]


(* 
 *   - - - - - - - - 2-forms
 *
 *  Suppose F = F_ij dx^i otimes dx^j where F_ij is antisymmetric in ij.
 *  Then we store elements (F_ij)_ {i<j} in F as a 4x4 matrix
 *
 *)


(*
 *   
 *)
emReadForm[form_, i_, j_] := Module[
    {tmpV1, tmpV2},

	tmpV1 = Min[{i,j}]; 
    tmpV2 = Max[{i,j}]; 
                      
    Signature[{i,j}] form[[tmpV1]][[tmpV2]]
                 
  ]


(*
 *  Apply (2,2)-tensor 'tensor' to 2-form 'form'.
 *)
emTensorForm[tensor_, form_] := Module[
	
	{res, i, j, a, b},
	
	res = Table[0, {i,1,4}, {j,1,4}];

	For[i = 1, i <= 4, i++,
		For[j = i+1, j <= 4, j++,
			
			res[[i]][[j]] = 1/2 Sum[emReadNormal[tensor, a,b,i,j] emReadForm[form, a,b], {a,1,4}, {b,1,4}];

		];
	];

	res

]


(* convert 2-form 'form' into a 6-vector in basis

    { dx^emOrd[i]  :  i = 1, ..., 6 }

 *)
emFormToVector[form_] := Module[
  	{vec, i, j},
  
  	(* allocate memory for 6-vector *)
  	vec = Table[0, {i, 1, 6}];
  	
  	For[i = 1, i <= 6, i++,

		vec[[i]] = emReadForm[form, emOrd[i][[1]],emOrd[i][[2]]];

   	];
  
  	vec
  
  ]


(* inverse of emFormToVector *)
emVectorToForm[vec_] := Module[

  	{res, i, sign, a, b},
  
  	(* allocate memory for (2,2)-tensor *)
  	res = Table[0, {i,1,4}, {j,1,4}];
  	
  	For[i = 1, i <= 6, i++,
   		
		a = Min[emOrd[i]];
		b = Max[emOrd[i]];

		sign = Signature[emOrd[i]];

		res[[a]][[b]] = sign vec[[i]];
	];

	res

]


emGeneralForm[sym_] := emVectorToForm[Table[ToExpression[sym<>ToString[i]],{i,1,6}]]


emPrintForm[form_]:=Module[

	{i, basis, vec},

	basis = Table["dx^"<>ToString[emOrd[i][[1]]-1]<>" \[And] dx^"<>ToString[emOrd[i][[2]]-1], {i,1,6}];
	vec = emFormToVector[form];

	For[i=1, i<= 6, i++,
		Print[basis[[i]] <> ": " <>ToString[InputForm[vec[[i]]]]];
	];

]


(* 
*   - - - - - - - - Coordinate change for (2,2)-tensor and 2-forms
* 
*)


emFormCoordinateChange[form_, L_]:=Module[


]


(* 
 * Perform a coordinate transformation on components of kappa
 *
 * Suppose L is the matrix 
 *
 *   L[[i]][[j]] = d tildex^i / dx^j
 *
 * Then L^-1 is the matrix L^-1[[i]][[j]] = d x^i / d tildex^j 
 * since then 
 *  (L o L^-1)[[i]][[j]] = L_ik L^-1_kj = d tildex^i/dx^k dx^k/d tildex^j = delta^i_j.
 *
 * If input kappa describe a (2,2)-tensor in x^i-coordinates, then 
 * the output of this routine describe kappa in tilde x^i coordinates.
 *
 *)
emCoordinateChange[kappa_, L_]:=Module[
	{i,j,m,n,p,q,r,s, kappaTransform, Linv},

	kappaTransform = emZeroKappa[];

	Linv = Inverse[L];
	
	For[i=1, i<=4, i++,
		For[j=i+1, j<=4, j++,
			For[m=1, m<=4, m++,
				For[n=m+1, n<=4, n++,
					kappaTransform[[i]][[j]][[m]][[n]] =  Sum[

						emReadNormal[kappa, p,q,r,s] *
							L[[i]][[p]] L[[j]][[q]] Linv[[r]][[m]] Linv[[s]][[n]],
							{p,1,4},{q,1,4},{r,1,4},{s,1,4}

					];
				];
			];
		];
	];

	kappaTransform

]
