function fresnelplot(n1,n2,k2,plotopt)
%
% fresnelplot(n1,n2,k2,plotopt) 
% 
% plots the Fresnel reflectivities and/or absorptivities (s-pol, p-pol 
% and circ.) for light incident from a dielectric medium with 
% refractive index n1 on an opaque absorbing medium with complex refractive 
% index n2+i*k2.
%
% Input:  n1      - refractive index (medium 1)
%         n2      - refractive index (medium 2)
%         k2      - extinction coefficient (medium 2)
%         plotopt - plotting option (type 'R' for plotting reflectivity,
%                   'A' for absorptivity or 'RA' for both)
%
% Output: plots of reflectivities and/or absorptivities
%
% Last updated: 2011-10-26 (David Bergström)
%

thetadeg = (0:0.1:90);
theta = thetadeg*pi/180;

[a,b,c] = intrc(n1,n2,k2,theta);

if nargin<4 || isempty(plotopt)
    fprintf('No plotting option specified. Type \''help fresnelplot\'' for further details.\n');
    return;
end;
if nargin==4 
    if ischar(plotopt)
        switch plotopt
            case {'R','r'}
                plot(thetadeg,a,thetadeg,b,thetadeg,c)
                legend('Rs','Rp','R','Location','northwest');
                xlabel('Angle of incidence \theta');
                ylabel('Reflectivity');
            case {'A','a'}
                plot(thetadeg,1-a,thetadeg,1-b,thetadeg,1-c)
                legend('As','Ap','A','Location','northwest');
                xlabel('Angle of incidence \theta');
                ylabel('Absorptivity');
            case {'RA','ra','AR','ar'}
                subplot(2,1,1)
                plot(thetadeg,a,thetadeg,b,thetadeg,c)
                legend('Rs','Rp','R','Location','northwest');
                xlabel('Angle of incidence \theta');
                ylabel('Reflectivity');
                subplot(2,1,2)
                plot(thetadeg,1-a,thetadeg,1-b,thetadeg,1-c)
                legend('As','Ap','A','Location','northwest');
                xlabel('Angle of incidence \theta');
                ylabel('Absorptivity');
            otherwise 
                fprintf('%s is not a valid option. Type \''help fresnelplot\'' for further details.\n',plotopt);
        end;
    else fprintf('Option must be a string. Type \''help fresnelplot\'' for further details.\n');
    end;
end;

function [Rs,Rp,R] = intrc(n1,n2,k2,theta)
%
% [Rs,Rp,R] = intrc(n1,n2,k2,theta) 
%
% calculates the reflection coefficients for s- and p-polarized components 
% and circular polarized (or unpolarized) reflectivity for light incident 
% from a perfect dielectric medium with refractive index n1 to an absorbing 
% medium with complex refractive index n2+i*k2 with angle of incidence 
% theta (given in radians). 
%
% Input:  n1    - refractive index (medium 1)
%         n2    - refractive index (medium 2)
%         k2    - extinction coefficient (medium 2)
%         theta - angle of incidence
%
% Output: Rs    - reflection coefficient of s-polarized component 
%         Rp    - reflection coefficient of p-polarized component
%         R     - circular polarized (or unpolarized) reflectivity
%
% Last updated: 2009-02-24 (David Bergström)
%

p = sqrt(1/2*(sqrt((n2.^2-k2.^2-n1^2*sin(theta).^2).^2+4*n2.^2.*k2.^2)+(n2.^2-k2.^2-n1^2*sin(theta).^2)));
q = sqrt(1/2*(sqrt((n2.^2-k2.^2-n1^2*sin(theta).^2).^2+4*n2.^2.*k2.^2)-(n2.^2-k2.^2-n1^2*sin(theta).^2)));

Rs = ((n1*cos(theta)-p).^2+q.^2)./((n1*cos(theta)+p).^2+q.^2);
Rp = ((p-n1*sin(theta).*tan(theta)).^2+q.^2)./((p+n1*sin(theta).*tan(theta)).^2+q.^2).*Rs;
R = (Rs+Rp)./2;
