We check that the Tamm-Rubilar tensor density defined in kappaLib coincides with the Tamm-Rubilar tensor density in 

    Hehl, Obukhov: Foundations of classical electromagnetics, 2003.

SetDirectory["~/writing/WIP/KappaLib/"];
<<KappaLib.m

Define a general (2,2)-tensor

kappa=emGeneralKappa["k"];
emKappaToMatrix[kappa]//MatrixForm

Step 1: Compute Fresnel equation for kappa using kappaLib

xi={xi0,xi1,xi2,xi3};
generalFresnel=emKappaToFresnel[kappa,xi];

Step 2: Compute Fresnel equation for kappa using equations D2.45 -- D2.49 in HO:2003.

Convenience functions

{AA,BB,CC,DD}=emKappaToABCD[kappa];
Am[i_,j_]:= AA[[i]][[j]]
Bm[i_,j_]:= BB[[i]][[j]]
Cm[i_,j_]:= CC[[i]][[j]]
Dm[i_,j_]:= DD[[i]][[j]]

Ctrace = Sum[Cm[i,i],{i,1,3}];
Dtrace = Sum[Dm[i,i],{i,1,3}];

Define M

M0=Det[AA];

Define M^a

M1[a_]:=-Sum[
Signature[{b,c,d}] (Am[b,a] Am[c,e]Cm[d,e]+Am[a,b] Am[e,c]Dm[e,d]),
{b,1,3},{c,1,3},{d,1,3},{e,1,3}]

Define M^ab

preM2term1[a_,b_]:=1/2 Am[a,b] ( 
Ctrace^2 +Dtrace^2-Sum[(Cm[c,d]+Dm[d,c])(Cm[d,c]+Dm[c,d]),{c,1,3},{d,1,3}])

preM2term2[a_,b_]:=Sum[(Cm[d,c]+Dm[c,d])(Am[c,a]Cm[b,d]+Dm[d,a]Am[b,c]),{c,1,3},{d,1,3}]

preM2terms345[a_,b_]:=-Ctrace Sum[Am[c,a] Cm[b,c],{c,1,3}]-Sum[Dm[c,a]Am[b,c],{c,1,3}]Dtrace-Sum[Am[d,c]Cm[a,c]Dm[d,b],{c,1,3},{d,1,3}]

preM2term6[a_,b_]:=Sum[
(Am[a,b] Am[d,c]-Am[d,a]Am[b,c])Bm[d,c],
{c,1,3},
{d,1,3}]

preM2[a_,b_]:=preM2term1[a,b]+preM2term2[a,b]+preM2terms345[a,b]+preM2term6[a,b]

M2[a_,b_]:=(preM2[a,b]+preM2[b,a])/2

Define M^abc

M3pre[a_,b_,c_]:=Sum[Signature[{d,e,c}] (Bm[d,f] (Am[a,b] Dm[e,f]-Dm[e,a]Am[b,f])+Bm[f,d](Am[a,b]Cm[f,e]-Am[f,a]Cm[b,e])+Cm[a,f]Dm[e,b]Dm[d,f]+Dm[f,a]Cm[b,e]Cm[f,d]),
{d,1,3},
{e,1,3},
{f,1,3}]

M3[a_,b_,c_]:=Module[{res,loop,perms},

perms=Permutations[{a,b,c}];
res=0;

For[loop=1, loop≤Length[perms],loop++,
res+=M3pre[perms[[loop]][[1]],perms[[loop]][[2]],perms[[loop]][[3]]];
];

res/Length[perms]
]

Define M^abcd

M4pre[a_,b_,c_,d_]:= Sum[
Signature[{e,f,c}]Signature[{g,h,d}] Bm[h,f] (1/2 Am[a,b] Bm[g,e]-Cm[a,e]Dm[g,b]),
{e,1,3},
{f,1,3},
{g,1,3},
{h,1,3}
]

M4[a_,b_,c_,d_]:=Module[{res,loop,perms},

perms=Permutations[{a,b,c,d}];
res=0;

For[loop=1, loop≤Length[perms],loop++,
res+=M4pre[perms[[loop]][[1]],perms[[loop]][[2]],perms[[loop]][[3]],perms[[loop]][[4]]];
];

res/Length[perms]
]

Compute Fresnel equation using equation D.2.44

qx={xi1,xi2,xi3};

decompFresnel=xi0^4 M0+xi0^3 Sum[qx[[i]] M1[i],{i,1,3}]+xi0^2 Sum[qx[[i]] qx[[j]] M2[i,j],{i,1,3},{j,1,3}]+xi0 Sum[qx[[i]] qx[[j]]qx[[k]] M3[i,j,k],{i,1,3},{j,1,3},{k,1,3}]+Sum[qx[[i]] qx[[j]]qx[[k]]qx[[l]] M4[i,j,k,l],{i,1,3},{j,1,3},{k,1,3},{l,1,3}];

Check that the two expressions for the Fresnel equation match

Expand[generalFresnel-decompFresnel]

Yes, they coincide.