(* ::Package:: *)

Print["KappaLib v1.0"];
(* 
 * 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 
 *
 *)


(*
 *  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. 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 is a (2,2)-tensor. 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. 
 *
 *)


(* convert kappa into a 6x6 matrix that represent kappa in basis

    { dx^1/\dx^2, dx^1/\dx^3, dx^1 /\ dx^4, dx^2 /\ dx^3, dx^2 /\ dx^4, dx^3 /\ dx^4 } 

*)

emKappaToMatrix[kappa_] := Module[
  	{pList, matRep, i, j},
  
  	(* list parameter pairs (i,j) with i<j and i,j=1,2,3,4 *)
  	pList = {{1, 2}, {1, 3}, {1, 4}, {2, 3}, {2, 4}, {3, 4}};
  
  	(* 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[[i]][[j]] = emReadNormal[kappa, 
        					pList[[i]][[1]], pList[[i]][[2]],
        					pList[[j]][[1]], pList[[j]][[2]]];
     		];
   	];
  
  	matRep
  
  ]


(* inverse of emlKappaToMatrix *)
emMatrixToKappa[matrix_] := Module[
  	{pList, matRep, i, j},
  
  	(* list parameter pairs (i,j) with i<j and i,j=1,2,3,4 *)
  	pList = {{1, 2}, {1, 3}, {1, 4}, {2, 3}, {2, 4}, {3, 4}};
  
  	(* allocate memory for (2,2)-tensor *)
  	res = emZeroKappa[];
  	
  	For[i = 1, i <= 6, i++,
   		For[j = 1, j <= 6, j++,
				res[[pList[[i]][[1]]]][[pList[[i]][[2]]]][[pList[[j]][[1]]]][[pList[[j]][[2]]]] = matrix[[i]][[j]];
     		];
	];
  
	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
]


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


(* 
 * 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

]
