Fresnel equation solver for multi-layer thin films

By Steven Byrnes. Last updated March 6, 2012. This program is available at my website, http://sjbyrnes.com. I'm releasing the entire contents to the public domain, but I'd appreciate it if you would let me know by email if you find any mistakes.

I recently ported this program into Python. The Python version is, I think, better written and easier to use, and also more powerful and general. For example, the Python version has a module for calculating structures that have both thin layers and "thick" layers (where the light's coherence is lost). You can find the Python program at http://sjbyrnes.com/ . This program has a "manual":  http://sjbyrnes.com/fresnel_manual.pdf . Please refer to the manual for "theoretical background" -- the equations used and their derivations. Actually, it is a manual for the Python version, but the differences are minor, such as the fact that Mathematica lists start with index 1 but Python lists start with index 0.

This program will use the Fresnel equations (in the "transfer matrix" technique) to calculate the percentages of light reflected, absorbed, and transmitted through a planar structure with multiple layers. Each layer must have a well-defined complex index of refraction, although it can vary with wavelength (I put in an example of this below). There's no scattering, no incoherence, and everything has to be perfectly planar (this is a one-dimensional simulation). The calculation has two varieties. The "Simple program" does just what I said: Total reflection percentage, and total transmission percentage. The "Spatially-resolved program" can also calculate, at a given location within a given layer, what the Poynting vector is, and what the absorbed energy density is. This is nice for solar cell modeling, for example, where you need to know precisely where the light is absorbed. On the other hand, the spatially-resolved program is a bit slower.

The "Simple program" also calculates the parameters "ψ" and "Δ" measured in ellipsometry.

WARNING: Some program limitations:
->The light must be incident from a non-dissipative medium, i.e. a real-valued index of refraction and a real-valued incident angle. (The final medium can be dissipative.) [The Python program does not have this restriction.]
->The real parts of all refractive indices must be positive (no "left-handed materials")
->The imaginary parts of all refractive indices must be nonnegative (no stimulated emission).
(Why? Snell's law requires taking arcsines that may yield complex-valued "angles". The program is therefore sensitive to the "branch cuts" that Mathematica puts into its ArcSin function. If you don't follow the above constraints, Mathematica will sometimes end up choosing the wrong (unphysical) arcsine value, and output a wrong answer. It should be possible to relax those constraints with some re-programming. See manual for further discussion and details.)

HOW TO USE THIS PROGRAM:
You don't need to understand all the steps of the calculation. Just execute those cells once, so that Mathematica knows the definitions of the functions, and then you can start making your output plots modeled on the examples I provide. If the code doesn't run properly, try deleting all my interspersed text comments, Mathematica sometimes mistakenly parses them as code.

Simple program

CLEAR VARIABLES BEFORE STARTING

ClearAll["Global`*"];

Subroutines

The four single-interface Fresnel amplitudes. "ts1" = amplitude for transmitted s-polarization at 1 interface. rp1 = amplitude for reflected p-polarization at 1 interface, etc.

ts1[nIncident_,nFinal_,θIncident_,θFinal_]:=((2nIncident Cos[θIncident])/(nIncident Cos[θIncident]+nFinal Cos[θFinal]));
tp1[nIncident_,nFinal_,θIncident_,θFinal_]:=((2nIncident Cos[θIncident])/(nFinal Cos[θIncident]+nIncident Cos[θFinal]));
rs1[nIncident_,nFinal_,θIncident_,θFinal_]:=((nIncident Cos[θIncident]-nFinal Cos[θFinal])/(nIncident Cos[θIncident]+nFinal Cos[θFinal]));
rp1[nIncident_,nFinal_,θIncident_,θFinal_]:=((nFinal Cos[θIncident]-nIncident Cos[θFinal])/(nFinal Cos[θIncident]+nIncident Cos[θFinal]));

Take the dot-product of a list of matrices:

RepeatedDotProduct[ListOfMatrices_]:=Block[{n=Length[ListOfMatrices],index,SoFar},
If[n==0,Return[{{1,0},{0,1}}]];
SoFar=ListOfMatrices[[1]];
For[index=2,index≤n,index++,SoFar=SoFar . ListOfMatrices[[index]]];
Return[SoFar];];

Actual program

INPUTS:
⟹"nList" is the list of refractive indices, in the order that the light would pass through them. The first element of the list should be the semi-infinite medium from which the light enters, the last element should be the semi-infinite medium to which the light exits (if any exits). The first element must be real, otherwise there are program errors associated with the branch cut discontinuities in ArcSin.
⟹"θ0" is the angle of incidence, θ0=0 is normal incidence, θ0=π/2 for glancing. θ0 should be real.
⟹"dList" is the list of layer thicknesses (front to back). It should correspond one-to-one to the elements of "nList".
⟹"λvacuum" is the vacuum wavelength of the light.

FresnelCalc[nList_,dList_,θ0_,λvacuum_]:=
Block[{NumMedia,δ,θList,tsList,tpList,rsList,rpList,MsList,MpList,MsTemp,MpTemp,Ms,Mp,rstotal,rptotal,tstotal,tptotal,Rs,Rp,Ts,Tp,RUnpol,TUnpol,temp,ψ,Δ},
NumMedia=Length[nList];
If[(Length[dList]≠NumMedia) || (dList[[1]]≠∞)||(dList[[NumMedia]]≠∞),Return["Error in d or n!"]];
If[Im[nList[[1]]]≠0 || Im[θ0]≠ 0,Return["Sorry, program does not work correctly for a dissipative incoming medium."]]; 
(* θList is a list with, for each layer, the angle that the light travels through the layer. Computed with Snell's law. Note that the "angles" may be complex!*)
θList=ArcSin[((nList[[1]] Sin[θ0])/nList)];
(* δ is a list with, for each layer (except the semi-infinite top and bottom "layers"),the phase accrued by traveling through that layer. δ[[i]] is what Sernelius would call (Subscript[δ, i]), apart from the renumbering.*)
δ=Table[If[i==1||i==NumMedia,Null,((2π dList[[i]])/(λvacuum))nList[[i]] Cos[θList[[i]]]],{i,1,NumMedia}];
(* tsList, tpList, rsList, rpList are the transmission and reflection amplitudes at each of the interfaces. For example,tsList[[i,j]] corresponds to what Sernelius would call (Subscript[t, i,j]) for s-polarization (apart from the renumbering),and is the (complex) ratio of (Subscript[E, transmitted])/(Subscript[E, original]), where the "original" comes from the i side and the "transmitted" is on the j side. We're only interested in the case that j=i+1; the rest don't enter the calculation.*)
tsList=Table[If[j-i==1,ts1[nList[[i]],nList[[j]],θList[[i]],θList[[j]]],Null],{i,1,NumMedia-1},{j,1,NumMedia}];
tpList=Table[If[j-i==1,tp1[nList[[i]],nList[[j]],θList[[i]],θList[[j]]],Null],{i,1,NumMedia-1},{j,1,NumMedia}];
rsList=Table[If[j-i==1,rs1[nList[[i]],nList[[j]],θList[[i]],θList[[j]]],Null],{i,1,NumMedia-1},{j,1,NumMedia}];
rpList=Table[If[j-i==1,rp1[nList[[i]],nList[[j]],θList[[i]],θList[[j]]],Null],{i,1,NumMedia-1},{j,1,NumMedia}];
(* At the interface between the (n-1)st and nth material, let (Subscript[x, n]) be the amplitude of the wave on the (n-1)st side heading towards the boundary, and let (Subscript[y, n]) be the amplitude of the wave on the (n-1)st side heading away from the boundary. Then (({{(Subscript[x, n])}, {(Subscript[y, n])}}))=(Subscript[M, n])(({{(Subscript[x, n+1])}, {(Subscript[y, n+1])}})), where (Subscript[M, n])=(1/(Subscript[t, n-1,n]))(({{1, (Subscript[r, n-1,n])}, {(Subscript[r, n-1,n]), 1}}))(({{Exp[-I (Subscript[δ, n])], 0}, {0, Exp[I (Subscript[δ, n])]}})),from Sernelius. (Subscript[M, 1]) and (Subscript[M, NumMedia]) are not defined. MsList[[n]] will be (Subscript[M, n]) for s-polarization, likewise MpList for p-polarization. Ms and Mp will be the matrix Sernelius calls "tilde-M", which you get by dot-producting all the (Subscript[M, i]), plus an extra factor at the end.*)
MsList=Table[If[i==1 || i== NumMedia,Null,(1/(tsList[[i-1,i]]))(({{1, rsList[[i-1,i]]}, {rsList[[i-1,i]], 1}})) . (({{Exp[-I δ[[i]]], 0}, {0, Exp[I δ[[i]]]}}))],{i,1,NumMedia}];
MpList=Table[If[i==1 || i==NumMedia,Null,(1/(tpList[[i-1,i]]))(({{1, rpList[[i-1,i]]}, {rpList[[i-1,i]], 1}})) . (({{Exp[-I δ[[i]]], 0}, {0, Exp[I δ[[i]]]}}))],{i,1,NumMedia}];
(* Ms and Mp will be the matrix Sernelius calls "tilde-M", which you get by dot-producting all the (Subscript[M, i]), plus an extra factor at the end.*)
MsTemp=RepeatedDotProduct[Take[MsList,{2,NumMedia-1}]];Ms = MsTemp . ((1/(tsList[[NumMedia-1,NumMedia]]))(({{1, rsList[[NumMedia-1,NumMedia]]}, {rsList[[NumMedia-1,NumMedia]], 1}})) );
MpTemp=RepeatedDotProduct[Take[MpList,{2,NumMedia-1}]];
Mp = MpTemp . ((1/(tpList[[NumMedia-1,NumMedia]]))(({{1, rpList[[NumMedia-1,NumMedia]]}, {rpList[[NumMedia-1,NumMedia]], 1}})) );
(*Net complex transmission and reflection amplitudes:*)
rstotal=((Ms[[2,1]])/(Ms[[1,1]]));
rptotal=((Mp[[2,1]])/(Mp[[1,1]]));
tstotal=(1/(Ms[[1,1]]));
tptotal=(1/(Mp[[1,1]]));
(*Net transmitted and reflected power,as a proportion of the incoming light power,for polarized and unpolarized incoming light. Don't forget,Ts,Tp,TUnpol represent the light power that enters the final, semi-infinite medium. If that medium has a complex index of refraction,the light will further attenuate as it travels through that medium to your detector,and you should take that into account separately.  See manual for discussion of formulas.*)
Rs=Abs[(rstotal^2)];
Rp=Abs[(rptotal^2)];
Ts=Abs[(tstotal^2)]×((Re[nList[[NumMedia]]Cos[θList[[NumMedia]]]])/(Re[nList[[1]]Cos[θList[[1]]]]));
Tp=Abs[(tptotal^2)]×((Re[nList[[NumMedia]]Conjugate[Cos[θList[[NumMedia]]]]])/(Re[nList[[1]]Cos[θList[[1]]]]));
RUnpol=(1/2)(Rs+Rp);
TUnpol=(1/2)(Ts+Tp);
(* Next, ψ and Δ are the parameters measured in ellipsometry. Note that both are in RADIANS, not degrees. To get the value in degrees, divide by the degree symbol, e.g. " ψ/° ". (In Mathematica, the degree symbol is just treated as the number π/180 ~ 0.0175.) ψ is defined to be between 0 and +π/2, and Δ  between -π and +π.  Warning: There are different conventions for Δ in use, you may need to use "π-Δ" or "2π-Δ" or "Δ-π", etc., to get agreement with your ellipsometer output.  *)
ψ=ArcTan[Abs[rptotal/rstotal]];
Δ=Arg[-rptotal/rstotal];
Return[{Rs,Ts,Rp,Tp,RUnpol,TUnpol,ψ,Δ}];]

Define separate functions for the reflection and transmission intensities for s-polarization, p-polarization, and unpolarized light.

Rs[nList_,dList_,θ0_,λvacuum_]:=FresnelCalc[nList,dList,θ0,λvacuum][[1]];
Ts[nList_,dList_,θ0_,λvacuum_]:=FresnelCalc[nList,dList,θ0,λvacuum][[2]];
Rp[nList_,dList_,θ0_,λvacuum_]:=FresnelCalc[nList,dList,θ0,λvacuum][[3]];
Tp[nList_,dList_,θ0_,λvacuum_]:=FresnelCalc[nList,dList,θ0,λvacuum][[4]];
RUnpol[nList_,dList_,θ0_,λvacuum_]:=FresnelCalc[nList,dList,θ0,λvacuum][[5]];
TUnpol[nList_,dList_,θ0_,λvacuum_]:=FresnelCalc[nList,dList,θ0,λvacuum][[6]];
ψ[nList_,dList_,θ0_,λvacuum_]:=FresnelCalc[nList,dList,θ0,λvacuum][[7]];
Δ[nList_,dList_,θ0_,λvacuum_]:=FresnelCalc[nList,dList,θ0,λvacuum][[8]];

Output -- Example plots

For these examples, I'm using using the SteveUnits package (Downloadable from http://sjbyrnes.com/). See package documentation for use. To check dimensional analysis, re-run the program from the start and confirm that the final answers don't change.

<<SteveUnits`

Example 1

Here's a thin non-absorbing layer, on top of a thick absorbing layer, with air on both sides. Plotting reflected intensity versus wavenumber, at two different incident angles.

dListTemp = {∞,100nm,300nm,∞};
nListTemp = {1,2.2,3.3+0.3I,1};
Plot[{RUnpol[nListTemp,dListTemp,0,(1/(k (cm^(-1))))],RUnpol[nListTemp,dListTemp,π/4,(1/(k (cm^(-1))))]},{k,0,100000},PlotRange->All,Frame->True,FrameLabel->{"Wavenumber ((cm^(-1)))","Fraction reflected"},PlotLabel->"Reflection of unpolarized light at normal incidence (blue), π/4 incidence (purple)",ImageSize->500]

Example 2

Here's the transmitted intensity versus wavelength through a single-layer film which has some complicated wavelength-dependent index of refraction. (I made these numbers up, but in real life they could be read out of a graph / table published in the literature.) Air is on both sides of the film, and the light is normally incident.

FilmRefractiveIndexList=
{{200nm,2.1+0.1I},
{300nm,2.4+0.3I},
{400nm,2.3+0.4I},
{500nm,2.2+0.4I},
{750nm,2.2+0.5I}};
FilmIndexFunction=Interpolation[FilmRefractiveIndexList];
dListTemp2={∞,300nm,∞};
nListTemp2[λ_]={1,FilmIndexFunction[λ],1};
Plot[TUnpol[nListTemp2[λ nm],dListTemp2,0,λ nm],{λ,200,750},PlotRange->All,Frame->True,FrameLabel->{"Vacuum wavelength (nm)","Fraction transmitted"},PlotLabel->"Transmission of unpolarized light at normal incidence",ImageSize->500]

Example 3

Here is a calculation of the ψ and Δ parameters measured in ellipsometry. This reproduces Fig. 1.14 in Handbook of Ellipsometry by Tompkins, 2005.

dListTemp3[dinnm_]:= {∞,dinnm nm,∞};
nListTemp3 = {1,1.46,3.87+0.02I};
Plot[{ψ[nListTemp3,dListTemp3[d],70 °,633 nm]/°,Δ[nListTemp3,dListTemp3[d],70°,633 nm]/°},{d,0,1000},PlotRange->All,Frame->True,FrameLabel->{"SiO2 thickness (nm)","ψ (blue) or Δ (purple) in degrees"},PlotLabel->"Ellipsometric parameters for air/SiO2/Si, varying SiO2 thickness.\nIncident light is: 70°, 633nm.\nShould agree with Handbook of Ellipsometry Fig. 1.14",ImageSize->500]

Spatially-resolved program

CLEAR VARIABLES BEFORE STARTING

ClearAll["Global`*"];

Subroutines

The four single-interface Fresnel amplitudes. "ts1" = amplitude for transmitted s-polarization for 1 interface. rp1 = amplitude for reflected p-polarization for 1 interface, etc.

ts1[nIncident_,nFinal_,θIncident_,θFinal_]:=((2nIncident Cos[θIncident])/(nIncident Cos[θIncident]+nFinal Cos[θFinal]));
tp1[nIncident_,nFinal_,θIncident_,θFinal_]:=((2nIncident Cos[θIncident])/(nFinal Cos[θIncident]+nIncident Cos[θFinal]));
rs1[nIncident_,nFinal_,θIncident_,θFinal_]:=((nIncident Cos[θIncident]-nFinal Cos[θFinal])/(nIncident Cos[θIncident]+nFinal Cos[θFinal]));
rp1[nIncident_,nFinal_,θIncident_,θFinal_]:=((nFinal Cos[θIncident]-nIncident Cos[θFinal])/(nFinal Cos[θIncident]+nIncident Cos[θFinal]));

Take the dot-product of a list of matrices:

RepeatedDotProduct[ListOfMatrices_]:=Block[{n=Length[ListOfMatrices],index,SoFar},
If[n==0,Return[{{1,0},{0,1}}]];
SoFar=ListOfMatrices[[1]];
For[index=2,index≤n,index++,SoFar=SoFar . ListOfMatrices[[index]]];
Return[SoFar];];

Actual programs

INPUTS:
⟹"nList" is the list of refractive indices, in the order that the light would pass through them. The first element of the list should be the semi-infinite medium from which the light enters, the last element should be the semi-infinite medium to which the light exits (if any exits). The first element must be real, otherwise there are program errors associated with the branch cut discontinuities in ArcSin.
⟹"θ0" is the angle of incidence, θ0=0 is normal incidence, θ0=π/2 for glancing. θ0 should be real.
⟹"dList" is the list of layer thicknesses (front to back). It should correspond one-to-one to the elements of "nList".
⟹"λvacuum" is the vacuum wavelength of the light.
⟹"LayerNumber" and "DistanceIntoLayer" give the exact point at which you want to calculate absorbed energy density, Poynting vector, etc. "LayerNumber" is the number of the layer  (e.g. LayerNumber=1 would be the semi-infinite medium from which the light enters), and "DistanceIntoLayer" is the orthogonal distance from the interface with the previous layer. LayerNumber must be more than 1, but it *can* refer to the final, semi-infinite medium.

The computation is split into two sub-programs, "SRFresnelCalc1" and "SRFresnelCalc2".

SRFresnelCalc1[nList_,dList_,θ0_,λvacuum_]:=
Block[{NumMedia,δ,θList,tsList,tpList,rsList,rpList,MsList,MpList,MsTemp,MpTemp,Ms,Mp,rstotal,rptotal,tstotal,tptotal,Rs,Rp,Ts,Tp,RUnpol,TUnpol,temp,xysList,xypList,vwsList,vwpList,i},
NumMedia=Length[nList];
If[(Length[dList]≠NumMedia) || (dList[[1]]≠∞)||(dList[[NumMedia]]≠∞),Print["Error in d or n!"];Return[]];
If[Im[nList[[1]]]≠0 || Im[θ0]≠ 0,Print["Sorry, program does not work correctly for a dissipative incoming medium."];Return[]]; 
(* θList is a list with, for each layer, the angle that the light travels through the layer. Computed with Snell's law. Note that the "angles" may be complex!*)
θList=ArcSin[((nList[[1]] Sin[θ0])/nList)];
(* δ is a list with, for each layer (except the semi-infinite top and bottom "layers"),the phase accrued by traveling through that layer. δ[[i]] is what Sernelius would call (Subscript[δ, i]), apart from the renumbering.*)
δ=Table[If[i==1||i==NumMedia,Null,((2π dList[[i]])/(λvacuum))nList[[i]] Cos[θList[[i]]]],{i,1,NumMedia}];
(* tsList, tpList, rsList, rpList are the transmission and reflection amplitudes at each of the interfaces. For example,tsList[[i,j]] corresponds to what Sernelius would call (Subscript[t, i,j]) for s-polarization (apart from the renumbering),and is the (complex) ratio of (Subscript[E, transmitted])/(Subscript[E, original]), where the \ comes from the i side and the \ is on the j side. We're only interested in the case that j=i+1; the rest don't enter the calculation.*)
tsList=Table[If[j-i==1,ts1[nList[[i]],nList[[j]],θList[[i]],θList[[j]]],Null],{i,1,NumMedia-1},{j,1,NumMedia}];
tpList=Table[If[j-i==1,tp1[nList[[i]],nList[[j]],θList[[i]],θList[[j]]],Null],{i,1,NumMedia-1},{j,1,NumMedia}];
rsList=Table[If[j-i==1,rs1[nList[[i]],nList[[j]],θList[[i]],θList[[j]]],Null],{i,1,NumMedia-1},{j,1,NumMedia}];
rpList=Table[If[j-i==1,rp1[nList[[i]],nList[[j]],θList[[i]],θList[[j]]],Null],{i,1,NumMedia-1},{j,1,NumMedia}];
(* At the interface between the (n-1)st and nth material, let (Subscript[x, n]) be the amplitude of the wave on the (n-1)st side heading towards the boundary, and let (Subscript[y, n]) be the amplitude of the wave on the (n-1)st side heading away from the boundary. Then (({{(Subscript[x, n])}, {(Subscript[y, n])}}))=(Subscript[M, n])(({{(Subscript[x, n+1])}, {(Subscript[y, n+1])}})), where (Subscript[M, n])=(1/(Subscript[t, n-1,n]))(({{1, (Subscript[r, n-1,n])}, {(Subscript[r, n-1,n]), 1}}))(({{Exp[-I (Subscript[δ, n])], 0}, {0, Exp[I (Subscript[δ, n])]}})),from Sernelius. (Subscript[M, 1]) and (Subscript[M, NumMedia]) are not defined. MsList[[n]] will be (Subscript[M, n]) for s-polarization, likewise MpList for p-polarization. Ms and Mp will be the matrix Sernelius calls "tilde-M", which you get by dot-producting all the (Subscript[M, i]), plus an extra factor at the end.*)
MsList=Table[If[i==1 || i== NumMedia,Null,(1/(tsList[[i-1,i]]))(({{1, rsList[[i-1,i]]}, {rsList[[i-1,i]], 1}})) . (({{Exp[-I δ[[i]]], 0}, {0, Exp[I δ[[i]]]}}))],{i,1,NumMedia}];
MpList=Table[If[i==1 || i==NumMedia,Null,(1/(tpList[[i-1,i]]))(({{1, rpList[[i-1,i]]}, {rpList[[i-1,i]], 1}})) . (({{Exp[-I δ[[i]]], 0}, {0, Exp[I δ[[i]]]}}))],{i,1,NumMedia}];
(* Ms and Mp will be the matrix Sernelius calls "tilde-M", which you get by dot-producting all the (Subscript[M, i]), plus an extra factor at the end.*)
MsTemp=RepeatedDotProduct[Take[MsList,{2,NumMedia-1}]];Ms = MsTemp . ((1/(tsList[[NumMedia-1,NumMedia]]))(({{1, rsList[[NumMedia-1,NumMedia]]}, {rsList[[NumMedia-1,NumMedia]], 1}})) );
MpTemp=RepeatedDotProduct[Take[MpList,{2,NumMedia-1}]];
Mp = MpTemp . ((1/(tpList[[NumMedia-1,NumMedia]]))(({{1, rpList[[NumMedia-1,NumMedia]]}, {rpList[[NumMedia-1,NumMedia]], 1}})) );
(*Net complex transmission and reflection amplitudes:*)
rstotal=((Ms[[2,1]])/(Ms[[1,1]]));
rptotal=((Mp[[2,1]])/(Mp[[1,1]]));
tstotal=(1/(Ms[[1,1]]));
tptotal=(1/(Mp[[1,1]]));
(* Beyond here isn't in Sernelius. We want to set up a calculation of absorption as a function of depth. We now know the full (({{(Subscript[x, n])}, {(Subscript[y, n])}})) vector for n=2 (the initial, incident state) -- it's  (({{1}, {rstotal}})) (for s-polarization). Note that we've set the incident amplitude to 1, which is just a nice convention. We'll now calculate (({{(Subscript[x, n])}, {(Subscript[y, n])}})) for the other values of n, using the variable names (({{(Subscript[x, n])}, {(Subscript[y, n])}})) = xysList[[n]] for s-polarization, or xypList[[n]] for p-polarization. Remember, (({{(Subscript[x, n+1])}, {(Subscript[y, n+1])}}))=((Subscript[M, n])^(-1))(({{(Subscript[x, n])}, {(Subscript[y, n])}})) , and we already computed the (Subscript[M, n]) earlier. *)
xysList=Table[Null,{i,1,NumMedia}];
xysList[[2]]=(({{1}, {rstotal}}));
For[i=3,i≤NumMedia,i++,xysList[[i]]=Inverse[MsList[[i-1]]] . xysList[[i-1]]];
xypList=Table[Null,{i,1,NumMedia}];
xypList[[2]]=(({{1}, {rptotal}}));
For[i=3,i≤NumMedia,i++,xypList[[i]]=Inverse[MpList[[i-1]]] . xypList[[i-1]]];
(* At the interface between the (n-1)st and nth material, let (Subscript[v, n]) be the amplitude of the wave on the n'th side heading away from the boundary, and let (Subscript[w, n]) be the amplitude of the wave on the n'th side heading towards the boundary. We can calculate (Subscript[v, n]) and (Subscript[w, n]) from (Subscript[x, n]) and (Subscript[y, n]), using a transformation matrix similar to (Subscript[M, n]). We use the variable names  (({{(Subscript[v, n])}, {(Subscript[w, n])}})) = vwsList[[n]] for s-polarization, or vwpList[[n]] for p-polarization. *)
vwsList=Table[Null,{i,1,NumMedia}];
For[i=2,i≤NumMedia,i++,vwsList[[i]]= Inverse[((1/(tsList[[i-1,i]]))(({{1, rsList[[i-1,i]]}, {rsList[[i-1,i]], 1}})) )].xysList[[i]]];
vwpList=Table[Null,{i,1,NumMedia}];
For[i=2,i≤NumMedia,i++,vwpList[[i]]= Inverse[((1/(tpList[[i-1,i]]))(({{1, rpList[[i-1,i]]}, {rpList[[i-1,i]], 1}})) )].xypList[[i]]];
(*Net complex transmission and reflection amplitudes:*)
rstotal=((Ms[[2,1]])/(Ms[[1,1]]));
rptotal=((Mp[[2,1]])/(Mp[[1,1]]));
tstotal=(1/(Ms[[1,1]]));
tptotal=(1/(Mp[[1,1]]));
(*Net transmitted and reflected power,as a proportion of the incoming light power,for polarized and unpolarized incoming light. See manual for discussion of formulas. *)
Rs=Abs[(rstotal^2)];
Rp=Abs[(rptotal^2)];
Ts=Abs[(tstotal^2)]×((Re[nList[[NumMedia]]Cos[θList[[NumMedia]]]])/(Re[nList[[1]]Cos[θList[[1]]]]));
Tp=Abs[(tptotal^2)]×((Re[nList[[NumMedia]]Conjugate[Cos[θList[[NumMedia]]]]])/(Re[nList[[1]]Cos[θList[[1]]]]));
RUnpol=(1/2)(Rs+Rp);
TUnpol=(1/2)(Ts+Tp);
Return[{Rs,Ts,Rp,Tp,RUnpol,TUnpol,vwsList,vwpList,θList}];]

SRFresnelCalc2[nList_,dList_,θ0_,λvacuum_,LayerNumber_,DistanceIntoLayer_]:=Block[{fc,vwsList,vwpList,θList,PhaseRate,AmplitudeVectorS,AmplitudeVectorP,PoyntingS,PoyntingP,AbsorbedS,AbsorbedP,fs,bs,fp,bp},
If[LayerNumber≤1 || LayerNumber > Length[nList] || LayerNumber∉Integers,Print["Bad LayerNumber"];Return[]];
If[DistanceIntoLayer>dList[[LayerNumber]],Print["Unphysical distance"];Return[]];
fc=SRFresnelCalc1[nList,dList,θ0,λvacuum];
vwsList=fc[[7]];
vwpList=fc[[8]];
θList=fc[[9]];
(*PhaseRate is the (phase accrued per distance) for an EM wave traveling forward within the layer in question.*) 
PhaseRate=((2π)/(λvacuum))nList[[LayerNumber]] Cos[θList[[LayerNumber]]];
(*AmplitudeVectorS and AmplitudeVectorP are two-element column vectors whose top entry is the amplitude of the forward-moving EM wave and whose bottom entry is the amplitude of the backwards-moving EM wave. S and P are the two polarizations.*)
AmplitudeVectorS= (({{Exp[I PhaseRate DistanceIntoLayer], 0}, {0, Exp[-I PhaseRate DistanceIntoLayer]}})) . vwsList[[LayerNumber]];
AmplitudeVectorP= (({{Exp[I PhaseRate DistanceIntoLayer], 0}, {0, Exp[-I PhaseRate DistanceIntoLayer]}})) . vwpList[[LayerNumber]];
fs=AmplitudeVectorS[[1,1]];
bs=AmplitudeVectorS[[2,1]];
fp=AmplitudeVectorP[[1,1]];
bp=AmplitudeVectorP[[2,1]];
(*PoyntingS is the Poynting vector (for s-polarization), projected along the "forward" direction (normal to the interfaces). It's given as a proportion of the (projected) Poynting vector of the incoming light. See manual for derivation of formulas. *)
PoyntingS=((Re[nList[[LayerNumber]]Cos[θList[[LayerNumber]]]Conjugate[fs+bs](fs-bs)])/(Re[nList[[1]]Cos[θList[[1]]]]));
(* p-polarization:*)
PoyntingP=((Re[nList[[LayerNumber]]Conjugate[Cos[θList[[LayerNumber]]]](fp+bp)Conjugate[fp-bp]])/(Re[nList[[1]]Cos[θList[[1]]]]));
(*Next, absorbed energy density. In principle this has units of [power]/[volume], but we'll output it as a multiple of incoming (projected) light power density (units [power]/[area]), so the output has units ((of [length])^(-1)). This is the negative derivative (with respect to distance) of the "Poynting" expressions above: *)
AbsorbedS=-Re[nList[[LayerNumber]]Cos[θList[[LayerNumber]]](I PhaseRate ((Abs[fs+bs])^2)-I Conjugate[PhaseRate]((Abs[fs-bs])^2))]/Re[nList[[1]]Cos[θList[[1]]]];
AbsorbedP=-(1/(Re[nList[[1]]Cos[θList[[1]]]]))Re[nList[[LayerNumber]]Conjugate[Cos[θList[[LayerNumber]]]](-I Conjugate[PhaseRate]((Abs[fp+bp])^2)+I PhaseRate ((Abs[fp-bp])^2))];
Return[{PoyntingS,PoyntingP,((PoyntingS+PoyntingP)/2),AbsorbedS,AbsorbedP,((AbsorbedS+AbsorbedP)/2)}];];

Define separate functions for the reflection and transmission intensities for s-polarization, p-polarization, and unpolarized light, and likewise the Poynting vectors and absorbed energy densities.

Rs[nList_,dList_,θ0_,λvacuum_]:=SRFresnelCalc1[nList,dList,θ0,λvacuum][[1]];
Ts[nList_,dList_,θ0_,λvacuum_]:=SRFresnelCalc1[nList,dList,θ0,λvacuum][[2]];
Rp[nList_,dList_,θ0_,λvacuum_]:=SRFresnelCalc1[nList,dList,θ0,λvacuum][[3]];
Tp[nList_,dList_,θ0_,λvacuum_]:=SRFresnelCalc1[nList,dList,θ0,λvacuum][[4]];
RUnpol[nList_,dList_,θ0_,λvacuum_]:=SRFresnelCalc1[nList,dList,θ0,λvacuum][[5]];
TUnpol[nList_,dList_,θ0_,λvacuum_]:=SRFresnelCalc1[nList,dList,θ0,λvacuum][[6]];
PoyntingS[nList_,dList_,θ0_,λvacuum_,LayerNumber_,DistanceIntoLayer_]:=SRFresnelCalc2[nList,dList,θ0,λvacuum,LayerNumber,DistanceIntoLayer][[1]];
PoyntingP[nList_,dList_,θ0_,λvacuum_,LayerNumber_,DistanceIntoLayer_]:=SRFresnelCalc2[nList,dList,θ0,λvacuum,LayerNumber,DistanceIntoLayer][[2]];
PoyntingUnpol[nList_,dList_,θ0_,λvacuum_,LayerNumber_,DistanceIntoLayer_]:=SRFresnelCalc2[nList,dList,θ0,λvacuum,LayerNumber,DistanceIntoLayer][[3]];
AbsorbedS[nList_,dList_,θ0_,λvacuum_,LayerNumber_,DistanceIntoLayer_]:=SRFresnelCalc2[nList,dList,θ0,λvacuum,LayerNumber,DistanceIntoLayer][[4]];
AbsorbedP[nList_,dList_,θ0_,λvacuum_,LayerNumber_,DistanceIntoLayer_]:=SRFresnelCalc2[nList,dList,θ0,λvacuum,LayerNumber,DistanceIntoLayer][[5]];
AbsorbedUnpol[nList_,dList_,θ0_,λvacuum_,LayerNumber_,DistanceIntoLayer_]:=SRFresnelCalc2[nList,dList,θ0,λvacuum,LayerNumber,DistanceIntoLayer][[6]];

Output -- Example plots

For these examples, I'm using using the SteveUnits package (Downloadable from http://sjbyrnes.com/). See package documentation for use. To check dimensional analysis, re-run the program from the start and confirm that the final answers don't change.

<<SteveUnits`

Example

Here's a device with two absorbing layers (the first 100nm thick, the second 300nm thick), and air on both sides. The light is λ=400nm violet light, s-polarized, incident at π/4 angle. We'll plot absorbed energy density as a function of position within the absorbing layer, in units of ((μm)^(-1)) (remember, this is absorbed energy density as a fraction of incoming light power density).

dListTemp = {∞,100nm,300nm,∞};
nListTemp = {1,2.2+0.2I,3.3+0.3I,1};
θ0=π/4;
λvac=400nm;
temp[d_]:=Piecewise[{{AbsorbedS[nListTemp,dListTemp,θ0,λvac,2,d nm]/((μm)^(-1)),d<100},{AbsorbedS[nListTemp,dListTemp,θ0,λvac,3,(d-100) nm]/((μm)^(-1)),d≥100}}]
Plot[temp[d],{d,0,400},PlotRange->All,Frame->True,FrameLabel->{"Distance into device (nm)","Absorbed energy (((μm)^(-1)))"},PlotLabel->"Absorbed light energy as a function of depth",ImageSize->500]

Note the discontinuity at 100nm -- that's the interface between the two layers.

Here's a cross-check...The integral of this should equal the total absorbed power. So the integral, plus reflected power, plus transmitted power, should be 1...

NIntegrate[AbsorbedS[nListTemp,dListTemp,θ0,λvac,2,d],{d,0,100nm}]+NIntegrate[AbsorbedS[nListTemp,dListTemp,θ0,λvac,3,d],{d,0,300nm}]+Rs[nListTemp,dListTemp,θ0,λvac]+Ts[nListTemp,dListTemp,θ0,λvac]

It works! Next, let's look at the Poynting vector for the same device:

dListTemp = {∞,100nm,300nm,∞};
nListTemp = {1,2.2+0.2I,3.3+0.3I,1};
θ0=π/4;
λvac=400nm;
temp[d_]:=Piecewise[{{PoyntingS[nListTemp,dListTemp,θ0,λvac,2,d nm],d<100},{PoyntingS[nListTemp,dListTemp,θ0,λvac,3,(d-100) nm],d≥100}}]
Plot[temp[d],{d,0,400},PlotRange->All,Frame->True,FrameLabel->{"Distance into device (nm)","Poynting vector"},PlotLabel->"Poynting vector as a function of depth",ImageSize->500]

There's a kink at the film interface (100nm), but the Poynting vector is a continuous function, as is required by physics. Good. Next, make sure that it agrees with the transmitted intensity and reflected intensity calculations:

{PoyntingS[nListTemp,dListTemp,θ0,λvac,4,0nm],Ts[nListTemp,dListTemp,θ0,λvac]}

{PoyntingS[nListTemp,dListTemp,θ0,λvac,2,0nm],1-Rs[nListTemp,dListTemp,θ0,λvac]}

Good!