function reflection = reflection_grapheneStack(efev,lambda,d)
% the Fermi energy (efev) is input in units of electron volts
% lambda (wavelength) and d (oxide thickness) are input in meters

c=3e8;          % speed of light (m/s)
n1 = 3.5;       % refractive index silicon
n2 = 1.47;      % refractive index silicon dioxide
n3 = 1;         % refractive index air
L = lambda;     %wavelength (m)
f = c./L;       %frequency (hz)
w = 2.*pi.*f;   % angular frequency
z0 = 377;       % free space impedance (ohms)

% define the fresnel coefficients and the propagation term
t12 = 2*n1 ./ (n1 + n2);
t21 = 2*n2 ./ (n1+n2);
r12 = ((n1 - n2) ./ (n1+n2));
r21 = (n2 - n1) ./ (n1+n2);
pd = exp(-1i.*w.*n2.*d./c);

% conductivity of graphene
sigma = GrapheneConductivity(efev,lambda);

% reflection term at graphene interface
ref = (n2 - n3 - z0.*sigma) ./ (n2+n3+z0.*sigma);

% calculate stack reflection (infinite paths) from taylor series
R = r12 + (t12*t21/r21).*(1./(1 - pd.*pd.*ref.*r21)-1);

% find the time averaged intensity
reflection = real(R.*conj(R));

end