Computations for Example 3.16

(* simplification needed for expression for the Poynting vector *)
Simplify[ Cos[x] Sin[x+Phi] - Sin[x] Cos[x + Phi]]

Compute expression for medium in local coordinates

coord={x,y,z};

Compute exterior derivative of a 1-form  in basis dy⋀dz,  dz⋀dx,  dx⋀dy

ddd[form_,i_,j_]:=D[form[[j]], coord[[i]]]-D[form[[i]], coord[[j]]]
ExtDer[oneForm_]:={
ddd[oneForm,2,3],
ddd[oneForm,3,1],
ddd[oneForm,1,2]};

(* check *)
GenOneForm={a[x,y,z],b[x,y,z],c[x,y,z]};
ExtDer[GenOneForm] //MatrixForm

(* Contract 1-form and VF. Returns scalar *)
ContractOneFormVF[OneForm_, VF_]:= OneForm.VF
(* Contract 2-form T and vector VF. Returns 1-form *)
ContractTwoFormVF[T_, VF_]:= {
T[[2]] VF[[3]]-T[[3]] VF[[2]],
T[[3]] VF[[1]]-T[[1]] VF[[3]],T[[1]] VF[[2]]-T[[2]] VF[[1]]}

Test: Show that the Reeb vector field for standard tight structure
 
     alpha = dz +/- x dy
 
 is R_+/- = d/dz

al={0,λ x, 1};  (* λ = +1 or -1 *)
GenReebVF={0, 0, 1};
Simplify[ContractOneFormVF[al,GenReebVF]==1]
Simplify[ContractTwoFormVF[ExtDer[al],GenReebVF]=={0,0,0}]

Find Reeb vector field for contact form
 
     alpha = Cos[ x] dz + λ Sin[ x] dy

alpha = {0, λ Sin[x], Cos[x]};  (* λ = +1 or -1 *)
ReebAlpha={0,λ Sin[x],Cos[x]};
Simplify[ContractOneFormVF[alpha,ReebAlpha]==1 /. λ^2 -> 1]
ContractTwoFormVF[ExtDer[alpha],ReebAlpha]/. λ^2 -> 1

Find Reeb vector field for contact form
 
     beta = Cos[ x+phi] dz + λ Sin[ x+phi] dy

beta = {0, λ Sin[x+ϕ], Cos[x+ϕ]};  (* λ = +1 or -1 *)
ReebBeta={0,λ Sin[x+ϕ], Cos[x+ϕ]};
Simplify[ContractOneFormVF[beta,ReebBeta]==1 /. λ^2 -> 1]
Simplify[ContractTwoFormVF[ExtDer[beta],ReebBeta]/. λ^2 -> 1]

Exterior product of two 1-forms (wrt. basis dy⋀dz,  dz⋀dx,  dx⋀dy)

ExtProduct[xi_, eta_]:= {
xi[[2]]eta[[3]]-xi[[3]]eta[[2]],
xi[[3]]eta[[1]]-xi[[1]]eta[[3]],
xi[[1]]eta[[2]]-xi[[2]]eta[[1]]
};
(*ExtProduct[{AA,BB,CC},{XX,YY,ZZ}]*)

Construct Le and Lm for medium in Example 3.16

ω = 1;
Le[xi_] := 1/ω ContractOneFormVF[xi,ReebAlpha]ExtDer[beta]- ExtProduct[alpha, xi]
Lm[xi_] := 1/ω ContractOneFormVF[xi,ReebBeta]ExtDer[alpha]- ExtProduct[beta, xi]

Lin2Vect[Lin_,vars_]:= {D[Lin,vars[[1]]],
D[Lin,vars[[2]]],
D[Lin,vars[[3]]]
}
Vect2Matrix[vec_,vars_]:={
Lin2Vect[vec[[1]],vars],
Lin2Vect[vec[[2]],vars],
Lin2Vect[vec[[3]],vars]
}

LeMat=Vect2Matrix[Le[{C1,C2,C3}],{C1,C2,C3}];
LmMat=Vect2Matrix[Lm[{C1,C2,C3}],{C1,C2,C3}];

(* Check that we have found the right matrix representations of Le and Lm *)
Simplify[LeMat.{C1,C2,C3}-Le[{C1,C2,C3}]]
Simplify[LmMat.{C1,C2,C3}-Lm[{C1,C2,C3}]]

LeMat //MatrixForm
LmMat //MatrixForm

Simplify[Simplify[Det[LeMat]]/. λ^2-> 1]
Simplify[Simplify[Det[LeMat]]/. λ^2 -> 1]

Define kappa tensor representing medium defined by Le and Lm.

From proof of Proposition 3.15 we know that

  κ^r0_i0 = 0,   κ^rs_ij=0
 κ^rs_i0 = -Q^rs_i        and κ^r0_ij=1/2 P^r_ij 

when P represent Le and Q represent Lm^-1

Perform decomposition of medium in Proposition 3.15

(* start with zero (2,2)-tensor *)
kappa = Table[0,
{i,1,4},{j,1,4},
{k,1,4},{l,1,4}
];
(* we only keep kappa^ij_mn updated for i<j and m<n. *)

kappaRead[a_,b_,i_,j_]:= Signature[{a,b}] Signature[{i,j}] kappa[[Min[{a,b}]]][[Max[{a,b}]]][[Min[{i,j}]]][[Max[{i,j}]]]

(* fill in P^r_ij-entries *)
For[i=2, i≤4,i++,
For[j=i, j≤4,j++,
For[r=2, r≤4,r++,
kappa[[1]][[r]][[i]][[j]] = -1/2 Sum[LeMat[[s]][[r-1]]Signature[{s,i-1,j-1}],{s,1,3}]
]
]
]

LmInv = Simplify[Simplify[Inverse[LmMat]]/. λ^2 -> 1];
LmInv//MatrixForm

(* fill in Q^rs_i-entries *)
For[r=2, r≤4,r++,
For[s=r, s≤4,s++,
For[i=2, i≤4,i++,
kappa[[r]][[s]][[1]][[i]] =1/2 Sum[LmInv[[i-1]][[u]]Signature[{u,r-1,s-1}],{u,1,3}]
]
]
]

kappa //MatrixForm

Check that we have defined medium κ correctly. That is, check that 
    G=κ F 
corresponds to 
    B = L_m H,   D = L_e E.

EE={e1,e2,e3};
(* Here B is defined such that BB_ab dx^a tensor dx^b 
= D12 dx^1 ⋀ dx^2 + D23 dx^2 ⋀ dx^3 + D31 dx^3 ⋀ dx^1 *) 
BB={
{0,B12, -B31},
{-B12,0,B23},
{B31,-B23,0}
};
(* B (wrt.basis dy⋀dz,dz⋀dx,dx⋀dy) *)
BBvec = {B23, B31, B12};

(* define 2-forms F and G. See Section 2 *)(* start with arbitrary F=(E,B) *)
Fform = Table[0,{i,1,4},{j,1,4}];
For[i=1, i≤ 3, i++,
Fform[[i+1]][[1]] = EE[[i]];
Fform[[1]][[i+1]] = -EE[[i]];
For[j=1, j≤ 3, j++,
Fform[[i+1]][[j+1]] = BB[[i]][[j]]
]
]
Fform //MatrixForm
kappaF = Table[
Sum[ kappaRead[i,j,a,b] Fform[[i]][[j]],
{i,1,4},{j,1,4}
],
{a,1,4},{b,1,4}
];

(* extract D and H from G=kappa(F) *)
DD = Table[0,{i,1,3},{j,1,3}];
HH = Table[0,{i,1,3}];
For[i=1, i≤ 3, i++,
HH[[i]] = -kappaF[[i+1]][[1]];
For[j=1, j≤ 3, j++,
DD[[i]][[j]] = kappaF[[i+1]][[j+1]];
]
]

Compare H extrated from G = kappa F and
 B = L_m H

(* HH only depends on B23, B31, B12 *)
HH //MatrixForm

tt=Vect2Matrix[HH,{B23,B31,B12}];
ttInv = Simplify[Simplify[Inverse[tt]]/.λ^2 -> 1];
Simplify[tt.{B23,B31,B12}==HH]

(* The last computation shows that tt is a 3x3 matrix that represents inverse of Lm 
from basis {dx^23, dx^31, dx^12} to basis {dx^1, dx^2, dx^3}.
If we have defined κ correctly we should then have tt^-1 = LmMat
*)
Simplify[ttInv == LmMat /. λ^2 -> 1]

Show that G = kappa F contains the constitutive equation D = L_e E

(* Check that DD is anti-symmetric *)
DD==-Transpose[DD]

(* Check that DD only depends on {e1,e2,e3} *)
DD //MatrixForm

(* Claim: Suppose D is a 2-form  
   D = D_ij dx^i tensor dx^j. 
Then 
    D^a = 1/2 ϵ^abc D_bc 
is representation of D in basis {dx^23, dx^31, dx^12} 
*)
Dtmp={
{0,dt12, -dt31},
{-dt12,0,dt23},
{dt31,-dt23,0}
};
Table[ 1/2 Sum[
Signature[{b,u,v}] Dtmp[[u]][[v]],
{u,1,3},{v,1,3}
],
{b,1,3}
]

(* transform DD_ij into basis
{dx^2⋀dx^3,dx^3⋀dx^1,dx^1⋀dx^2} *) DDrightBasis=Table[ 1/2 Sum[Signature[{b,u,v}] DD[[u]][[v]],{u,1,3},{v,1,3}],{b,1,3}]

(* Compute transformation from E to D with respect to bases {dx^1, dx^2, dx^3} and {dx^23, dx^31, dx^12} *)
tt = Vect2Matrix[DDrightBasis,{e1,e2,e3}];
Simplify[tt==LeMat/. λ^2-> 1]//MatrixForm

Extra:

(*Lemma: For any antisymmetric a_ij we have
A_uv = 1/2 ϵ_uva ϵ^abc A_bc 
*)

(* proof: *)
anti={
{0,a12, -a31},
{-a12,0,a23},
{a31,-a23,0}
};
anti - 1/2Table[Sum[ Signature[{uu,vv,aa}] Signature[{aa,bb,cc}] anti[[bb]][[cc]],{aa,1,3},{bb,1,3},{cc,1,3}] ,
{uu,1,3},{vv,1,3}
]

Decompose medium kappa into principal, skewon, and axion part

(* Verify that κ has no axion part. See Proposition 3.15 *)
TraceKappa = Sum[ kappaRead[i,j,i,j],{i,1,4},{j,1,4}]

KappaStrike[i_,j_]:= Sum[ kappaRead[i,m,j,m],{m,1,4}]
Table[KappaStrike[i,j],{i,1,4},{j,1,4}] //MatrixForm

PreSkewon[i_,j_,l_,m_]:= 2 KappaStrike[i,l] KroneckerDelta[j,m];
PreSkewonII[i_,j_,l_,m_]:= 1/2 (PreSkewon[i,j,l,m]-PreSkewon[j,i,l,m])
SkewonPart[i_,j_,l_,m_]:= 1/2 (PreSkewonII[i,j,l,m]-PreSkewonII[i,j,m,l])

PrincipalPart[i_,j_,l_,m_]:=kappaRead[i,j,l,m]-SkewonPart[i,j,l,m]

Extract ϵ-matrix for skewon and principal medium

EE={e1,e2,e3};
(* Here B is defined such that BB_ab dx^a tensor dx^b 
= D12 dx^1 ⋀ dx^2 + D23 dx^2 ⋀ dx^3 + D31 dx^3 ⋀ dx^1 *)

(* define 2-forms F and G. See Section 2 *)(* start with arbitrary F=(E,B) *)
Fform = Table[0,{i,1,4},{j,1,4}];
For[i=1, i≤ 3, i++,
Fform[[i+1]][[1]] = EE[[i]];
Fform[[1]][[i+1]] = -EE[[i]];
]
kappaFsk = Table[
Sum[ SkewonPart[i,j,a,b] Fform[[i]][[j]],
{i,1,4},{j,1,4}
],
{a,1,4},{b,1,4}
];
kappaFpr= Table[
Sum[ PrincipalPart[i,j,a,b] Fform[[i]][[j]],
{i,1,4},{j,1,4}
],
{a,1,4},{b,1,4}
];

(* extract D from G=kappa(F) *)
DDsk = Table[0,{i,1,3},{j,1,3}];
DDpr = Table[0,{i,1,3},{j,1,3}];
For[i=1, i≤ 3, i++,
For[j=1, j≤ 3, j++,
DDsk[[i]][[j]] = kappaFsk[[i+1]][[j+1]];
DDpr[[i]][[j]] = kappaFpr[[i+1]][[j+1]];
]
]

DDskRightBasis=Table[ 1/2 Sum[Signature[{b,u,v}] DDsk[[u]][[v]],{u,1,3},{v,1,3}],{b,1,3}] /. λ^2-> 1;
DDprRightBasis=Table[ 1/2 Sum[Signature[{b,u,v}] DDpr[[u]][[v]],{u,1,3},{v,1,3}],{b,1,3}] /. λ^2-> 1;

(* Compute transformation from E to D with respect to bases {dx^1, dx^2, dx^3} and {dx^23, dx^31, dx^12} *)
ttsk = Simplify[Vect2Matrix[DDskRightBasis,{e1,e2,e3}]];
ttpr = Simplify[Vect2Matrix[DDprRightBasis,{e1,e2,e3}]];

ttsk //MatrixForm

ttpr //MatrixForm

{Det[ttpr],Det[ttsk],Simplify[Simplify[Det[tt]]/.λ^2-> 1]}

Simplify[ttsk + ttpr == tt /. λ^2 -> 1]

LeMatSym = Simplify[1/2 (LeMat + Transpose[LeMat]) /. λ^2 -> 1];
LeMatASym = Simplify[1/2 (LeMat - Transpose[LeMat]) /. λ^2 -> 1];

Simplify[LeMatSym == ttpr]
Simplify[LeMatASym == ttsk]