SetDirectory["~/KappaLib/"];
<<kappaLib-1.2.m
<<helper.m

Define 4-parameter medium kappa such that

1) each kappa is algebraically decomposable

2) beta^2-alpha gamma = 1

3) Fresnel surface always factorises into two linear factors 
    and a quadratic form. For suitable parameters, the quadratic 
    form can have a Lorentz signature
    
4) The nonlinear equation for D has no solution.

kappa=emMatrixToKappa[(({{tau+C1, -1, 0, 0, 0, 0}, {C1, tau, 1, 0, 0, 0}, {1, 0, tau, 0, 0, 0}, {0, 0, 1, tau, -1, 0}, {0, 1, 0, 0, tau, 1}, {0,  C3, 0, 0, C2, tau-C1}}))];

(* kappa can invertible *)
FullSimplify[emDet[kappa]]

(* kappa has axion component if and only if tau = 0*)
FullSimplify[emTrace[kappa]]

(* Since this does not simplify to {0}, kappa has a skewon component *)
Union[Flatten[kappa-emPoincare[kappa]]]

(* Since this does not simplify to {0}, kappa has a principal component *)
Union[Flatten[kappa+emPoincare[kappa]]]

Define parameters to show that kappa is algebraically decomposable

Abivector=(({{0, 1, 0, 0}, {-1, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}));
Bbivector=(({{0, 0, C1+C3, C1}, {0, 0, 0, C1-C2}, {-(C1+C3), 0, 0, -1}, {-C1, -(C1-C2), 1, 0}}));
alpha=tau^2-1;
beta=-tau;
gamma=1;
rho=1/2;

(*  Verify that kappa satisfies equation (44)   *)
RHS = alpha emIdentityKappa[] + beta(kappa + emPoincare[kappa]) + gamma emCompose[emPoincare[kappa],kappa];
LHS = emBiProduct[rho,Abivector,Bbivector]+emBiProduct[rho,Bbivector,Abivector];
eqs=Union[Flatten[RHS-LHS]];
Union[Simplify[eqs]]

Write equations that should be satisfied for D

(* contract kappa with a bivector from the left *)
contract[biv_,kappa_]:=Table[
1/2 Sum[
biv[[i]][[j]]emReadNormal[kappa,a,b,i,j]
,
{i,1,4},{j,1,4}
]
,
{a,1,4},{b,1,4}
]

There is no bivector D that satisfies 

   D (kappa + beta Id) = 1/2 trace (rho bar(D) otimes D)  A + B.

Dbivector=emMatrix["d",4,Structure->"AntiSymmetric"];
dLHS = contract[Dbivector,kappa+beta emIdentityKappa[]];
dRHS = 1/2 emTrace[emBiProduct[rho,Dbivector,Dbivector]]Abivector+Bbivector;
deqs=Union[Flatten[dLHS-dRHS]];

deqs

Fresnel surface factorises

coords={xi0,xi1,xi2,xi3};
fr=FullSimplify[emKappaToFresnel[kappa,coords]];

Mat=(({{0, C1+C3, -2 C1-C3, -C2-C1 C3}, {C1+C3, 0, -C1+C2, 1+(C1^2)-C1 C2}, {-2 C1-C3, -C1+C2, 2, -C1}, {-C2-C1 C3, 1+(C1^2)-C1 C2, -C1, -2 C2}}));
LinFactors=xi0 (xi0+xi1+xi3);
Simplify[fr-1/2 LinFactors (coords.Mat.coords)]

For suitable choices of C1,C2,C3, the quadratic factor has a Lorentz signature

sub={C1->0,C2->1/2,C3->1/2,tau->2};

FullSimplify[Det[Mat/.sub]]
FullSimplify[emDet[kappa/.sub]]