Needs["Calculus`VectorAnalysis`"];
Needs["LinearAlgebra`MatrixManipulation`"];

Πfunction[ϵ1_,  ϵ2_,  θ_,  nx_, ny_, nz_]:=
Module[{κ, Sinθ0, Cosθ0, Sinφ0, 
            Cosφ0, Sin2θ0, Cos2θ0, m, mm , ν, q, 
            e, mag, ϵ, u, Π },

m  =  {nx, ny, nz};
mm  =  
              Table[m[[i]] m[[j]], {i, 1, 3}, {j, 1, 3}];
  
Cosθ0  =  nz;
Sinθ0  =  (Sqrt[1 - Cosθ(0^2)]);
If[Sinθ0 > 0, 
              Cosφ0  =  nx/Sinθ0, Cosφ0 = 1];
If[Sinθ0 > 0, 
              Sinφ0  =  ny/Sinθ0, Sinφ0 = 0];
Sin2θ0  =  
              2 Sinθ0 Cosθ0;
Cos2θ0  =  
              Cosθ(0^2) - Sinθ(0^2);

  
κ  =  -Sin[θ  π/
                      180];

  
ϵ[
                m]  =  ϵ1*
                  IdentityMatrix[3] + (ϵ2 - ϵ1) mm;
		  
  
ν[
                1]  =  (((ϵ1 - ϵ2) κ Cosφ0  Sin2θ0)/(ϵ1 + ϵ2 + (-ϵ1 + ϵ2) Cos2θ0)) + (√(ϵ1 (ϵ1 (ϵ2 - (κ^2)) - (ϵ1 - ϵ2) (ϵ2 - (κ^2)) Cosθ(0^2) + (ϵ1 - ϵ2) (κ^2) Cosφ(0^2) Sinθ(0^2))))/(-ϵ1 + (ϵ1 - ϵ2) Cosθ(0^2)); 
ν[2] = -(Sqrt[ϵ1 - (κ^2)]);    
  
ν[
                3]  =  (((ϵ1 - ϵ2) κ Cosφ0  Sin2θ0)/(ϵ1 + ϵ2 + (-ϵ1 + ϵ2) Cos2θ0)) - (√(ϵ1 (ϵ1 (ϵ2 - (κ^2)) - (ϵ1 - ϵ2) (ϵ2 - (κ^2)) Cosθ(0^2) + (ϵ1 - ϵ2) (κ^2) Cosφ(0^2) Sinθ(0^2))))/(-ϵ1 + (ϵ1 - ϵ2) Cosθ(0^2));
ν[4] = (Sqrt[ϵ1 - (κ^2)]);   

  
q[1]  =  {κ, 0, ν[1]};   
q[2]  =  {κ, 0, ν[2]};   
q[3]  =  {κ, 0, ν[3]};   
q[4]  =  {κ, 0, ν[4]};   
  
		(* In general, one must make sure that m x k ≠ 0, otherwise the choice of eigenvector has an extra degree of freedom, and we chose the two vectors, normal to m, that together form a direct triedre *)
e[
                1]  =  -(q[1] . q[1]) m  +  (q[1] . m) 
                  q[1];
mag[e1]  =  Sqrt[e[1] . e[1]];
  
If[mag[e1]  ==  0, 
              e[1]  =  {-Cosθ0 Cosφ0, -Cosθ0 Sinφ0, Sinθ0}, e[1]  =  e[1]/mag[e1]];
e[2] =  CrossProduct[m, q[2]];
mag[e2]  =  Sqrt[e[2] . e[2]];
  
If[mag[e2]  ==  0, 
              e[2]  =  {-Sinφ0, Cosφ0, 0} , 
              e[2]  =  e[2]/mag[e2]];e[
                3]  =  (q[3] . q[3]) m  -  (q[3] . m) q[3] ;
mag[e3]  =  Sqrt[e[3] . e[3]];
  
If[mag[e3]  ==  0, 
              e[3]  =  {Cosθ0 Cosφ0, 
                  Cosθ0 Sinφ0, -Sinθ0}, 
              e[3]  =  e[3]/mag[e3]];
e[4] =  -CrossProduct[m, q[4]];
mag[e4]  =  Sqrt[e[4] . e[4]];
  
If[mag[e4]  ==  0, 
              e[4]  =  {-Sinφ0, Cosφ0, 0}, 
              e[4]  =  e[4]/mag[e4]];

		  
ϵ[-1]  =  Inverse[ϵ[m]];
Do[u[i]  =  ϵ[-1] . e[i], {i, 1, 4}];
Do[
              u[i]  =  u[i] / Sqrt[u[i] . u[i]], {i, 1, 4}];
Π[LM]  =   Transpose[
                Table[{u[m][[1]], u[m][[2]], CrossProduct[
                        q[m], u[m]][[1]], CrossProduct[q[m], 
                        u[m]][[2]]}, {m, 1, 4}]]
]

ΠfunctionIsotropic[ϵ0_, θθ0_] := 
Module[{θ}, 
θ = θθ0 π/
              180. ; 
Πmatrix = {{-((Sqrt[ϵ0 - Sin[θ]^2])/(Sqrt[ϵ]))0, 
              0, ((Sqrt[ϵ0 - Sin[θ]^2])/(Sqrt[ϵ]))0, 0}, {0, 
              1. , 0, 1. }, {0, (Sqrt[ϵ0 - Sin[θ]^2]), 
              0, -(Sqrt[ϵ0 - Sin[θ]^2])}, {(Sqrt[ϵ])0, 
              0, (Sqrt[ϵ])0, 0}}
]

Φfunction[ϵ1_,  ϵ2_,  λ_, dd_,  θ_,  nx_, ny_, nz_]:=
Module[{d, κ, k0, Sinθ0, Cosθ0, 
            Sinφ0, Cosφ0, Sin2θ0, Cos2θ0, 
            m , ν, Φ},

d  =  dd (10^(-8));
k0 = 
              2.  π/(λ  (10^(-7)));
m  =  {nx, ny, nz};
  
Cosθ0  =  nz;
Sinθ0  =  (Sqrt[1 - Cosθ(0^2)]);
If[Sinθ0 > 0, 
              Cosφ0  =  nx/Sinθ0, Cosφ0 = 1];
If[Sinθ0 > 0, 
              Sinφ0  =  ny/Sinθ0, Sinφ0 = 0];
Sin2θ0  =  
              2 Sinθ0 Cosθ0;
Cos2θ0  =  
              Cosθ(0^2) - Sinθ(0^2);

  
κ  =  -Sin[θ  π/
                      180];

  
ϵ[
                m]  =  ϵ1*
                  IdentityMatrix[3] + (ϵ2 - ϵ1) mm;
		  
  
ν[
                1]  =  (((ϵ1 - ϵ2) κ Cosφ0  Sin2θ0)/(ϵ1 + ϵ2 + (-ϵ1 + ϵ2) Cos2θ0)) + (√(ϵ1 (ϵ1 (ϵ2 - (κ^2)) - (ϵ1 - ϵ2) (ϵ2 - (κ^2)) Cosθ(0^2) + (ϵ1 - ϵ2) (κ^2) Cosφ(0^2) Sinθ(0^2))))/(-ϵ1 + (ϵ1 - ϵ2) Cosθ(0^2)); 
ν[2] = -(Sqrt[ϵ1 - (κ^2)]);    
  
ν[
                    3]  =  (((ϵ1 - ϵ2) κ Cosφ0  Sin2θ0)/(ϵ1 + ϵ2 + (-ϵ1 + ϵ2) Cos2θ0)) - (√(ϵ1 (ϵ1 (ϵ2 - (κ^2)) - (ϵ1 - ϵ2) (ϵ2 - (κ^2)) Cosθ(0^2) + (ϵ1 - ϵ2) (κ^2) Cosφ(0^2) Sinθ(0^2))))/(-ϵ1 + (ϵ1 - ϵ2) Cosθ(0^2));
ν[4] = (Sqrt[ϵ1 - (κ^2)]);   

Φ[LM]  = 
              Table[If[i == j, Exp[I k0 ν[i] d ], 0] , {i, 1, 
                  4}, {j, 1, 4}]

]

AbsRef[p1_, p2_, p3_, p4_, p5_, p6_, p7_, p8_]:=
Module[{κ, θ, φ, ϵ1, ϵ2, 
            d, δi, δr, ℰ, 
            nn, Π, Φ, M, TM, u, air, water, LM, 
            reflected, inwave, refwave, α, 
            analyzer, ℛ},

θ = p1 π/180;ϵ1 = p2;ϵ2 = p3;d = p4;λ = p5;δi = p6 ;δr = p7 ;φ = p8 π/180;

  
nn  =  {Sin[θ] 
                  Cos[φ], Sin[θ] Sin[φ], 
                Cos[θ]};

Π[
                air]  =  ΠfunctionIsotropic[1.0, 53.0];
Π[
                water]  =  ΠfunctionIsotropic[1.77, 53.0];
Π[
                LM] = Πfunction[ϵ1, ϵ2, 53.0, 
                nn[[1]], nn[[2]], nn[[3]]];
Φ[
                LM]  =  Φfunction[ϵ1, ϵ2, λ, d, 53.0, nn[[1]], nn[[2]], nn[[3]]];

M[air, LM]  =  
              Inverse[Π[air]] . Π[LM];
M[LM, water]  =  
              Inverse[Π[LM]] . Π[water];
		
		  
    TM  =  
              M[air, LM] . Conjugate[Φ[LM]] . M[LM, water];
	
	  ;

u[air][
                1]  =  {-Cos[53.  π/180], 0, Sin[53.  π/180]};   
    u[air][2]  =  {0, 1, 0};   
    u[air][3]  =  { Cos[53.  π/180], 0,  
                Sin[53.  π/180]};   
    u[air][4]  =  {0, 1, 0};   	

reflected[air]  =  
              SubMatrix[TM, {3, 1}, {2, 2}] . 
                Inverse[
                  SubMatrix[
                    TM, {1, 1}, {2, 2}]] . {ℰ[air][
                    1], ℰ[air][2]};
		  
refwave[
                air]  =  reflected[air][[1]] u[air][
                    3]  +  reflected[air][[2]] u[air][4];
inwave[
                air]  =  ℰ[air][1] u[air][
                    1]  +  ℰ[air][2] u[air][2] ;
		  
α[
                analyzer] =  u[air][3] 
                  Cos[δr π/180] + u[air][4] 
                  Sin[δr π/180] ;
		  
ℛ = (Abs[refwave[air] . α[analyzer]]^2/
                      Abs[inwave[air] . inwave[air]])/.{ℰ[air][1] ->  
                      Cos[δi π/180], ℰ[air][
                        2] ->  Sin[δi  π/180]}

	]

\(\(Reflectivity\  = \ 
	
        Compile[{p1, p2, p3, p4, p5, p6, p7, p8}, 
          AbsRef[p1, p2, p3, p4, p5, p6, p7, p8]];\)\)

BrewsterAngle  =  
        53;
ϵ[water]  =  1.77;
ϵ[air]  =  1.0;
ϵ[LM][perp]  =  2.16;
ϵ[LM][paral] = 2.37;
λ  =  640.0;  
length  =  28.0   ; 
polarizer  =  0.0;
analyzer  =  60.0;
θ0  =  20.0 ;
φ0  =  60.0 ;

RefAzimuthTable  =  
        Table[{φ0, 
            Reflectivity[θ0, ϵ[LM][
                perp], ϵ[LM][paral], 
              length*Cos[(π/180) θ0], λ, polarizer, 
              analyzer, φ0]}, {φ0, 0, 360, 
            30}];

RefAzimuth  =  Interpolation[RefAzimuthTable];

refl[φ_]  =  
        RefAzimuth[Mod[φ, 360]];

Plot[refl[x], {x, 0, 360}, Axes -> True]