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

Metaclass I:

vars={x0,x1,x2,x3};
kappa=emMatrixToKappa[(({{a1, 0, 0, -b1, 0, 0}, {0, a2, 0, 0, -b2, 0}, {0, 0, a3, 0, 0, -b3}, {b1, 0, 0, a1, 0, 0}, {0, b2, 0, 0, a2, 0}, {0, 0, b3, 0, 0, a3}}))];
fr=emKappaToFresnel[kappa,vars];

FullSimplify[emDet[kappa]]

Write Fresnel polynomial using D0, D1, D2, D3.

subSym={
D1 ->((((((a2-a3))^2)+(b2^2)+(b3^2)))/(b2 b3)),
D3->((((((a1-a2))^2)+(b1^2)+(b2^2)))/(b1 b2)),
D2->((((((a1-a3))^2)+(b1^2)+(b3^2)))/(b1 b3)),
D0->2((a1 (b2^2-b3^2 )+a2 ((b3^2)-(b1^2))+a3 (b1^2-b2^2))-(a1-a2) (a1-a3) (a2-a3))/(b1 b2 b3)
};

frSym=(x0^4)+(x1^4)+(x2^4)+(x3^4)-D0 x0 x1 x2 x3+D1 ((x2^2) (x3^2)-(x0^2) (x1^2))+D2 ((x1^2) (x3^2)-(x0^2) (x2^2))+D3 ((x1^2) (x2^2)-(x0^2) (x3^2))

Simplify[fr-b1 b2 b3 frSym/.subSym]

Verify: If D1 = D2 = D3 = 2 then the Fresnel surface is a single light cone

FullSimplify[frSym/.{D1->2,D2->2,D3->2,D0->0}]

Implicit equation for D0

Simplify[(D0^2)+4 (-4+(D1^2)+(D2^2)+(D3^2)-D1 D2 D3)/.subSym]

We assume that the Fresnel polynomial factorises:

A=Table[ToExpression["A"<>ToString[Min[{i,j}]]<>ToString[Max[{i,j}]]],{i,0,3},{j,0,3}];
B=Table[ToExpression["B"<>ToString[Min[{i,j}]]<>ToString[Max[{i,j}]]],{i,0,3},{j,0,3}];
A//MatrixForm
B//MatrixForm
factorised=(vars.A.vars) (vars.B.vars);

cons=Union[Flatten[CoefficientList[factorised-frSym,vars]]];
cons=simp[cons];
show[cons]

Equation (2): By rescaling we may assume that A00 = 1.

sub={A00->1,B00->1};
cons=simp[cons//.sub];
show[cons]

tmp=Take[cons,{5,7}]

Solve[toEqs[tmp],{B01,B02,B03}]

sub=Join[sub,%[[1]]]

tmp=Join[Take[cons,{17,19}],Take[cons,{23,25}]];
tmp//MatrixForm

Solve[toEqs[tmp],{B33,B22,B11,B12,B13,B23}]

sub=Join[sub,%[[1]]]

cons=simp[cons//.sub];
show[cons]

Solve[cons[[20]]==0,D0]

sub=Join[sub,%[[1]]]

cons=simp[cons//.sub];
show[FullSimplify[cons]]

Variables[cons]

elimVars=Variables[A]
condVars={D1,D2,D3}

gb=simp[GroebnerBasis[cons,condVars,elimVars]];//Timing

show[gb]

Since we assume that 2 <= D_1 <= D_2 <= D_3 and D_2, D_3>2 the first equation implies that D_1 = 2

show[simp[gb/.D1->2]]

.. and since D_2, D_3>2, equation (2) implies that D_2 = D_3. We have shown that D1 = 2 and D2=D3>2.

D0 vanishes

Simplify[4 (-4+(D1^2)+(D2^2)+(D3^2)-D1 D2 D3)/.{D1->2, D2->D3}]

Verify decomposition

AA= (({{1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, (1/2) (-D3+(Sqrt[-4+(D3^2)])), 0}, {0, 0, 0, (1/2) (-D3+(Sqrt[-4+(D3^2)]))}}));
BB = (({{1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, (1/2) (-D3-(Sqrt[-4+(D3^2)])), 0}, {0, 0, 0, (1/2) (-D3-(Sqrt[-4+(D3^2)]))}}));

verify=(vars.AA.vars) (vars.BB.vars);

Simplify[verify-frSym/.{D0->0,D1->2,D2->D3}]