initial

The optics_toolkit code taken from
http://mercury.pr.erau.edu/~greta9a1/downloads/index.html

the older version is also available at mathwork web site
http://www.mathworks.com/matlabcentral/fileexchange/15459-basic-paraxial-optics-toolkit

3 files changed, 459 insertions, 0 deletions

diff --git a/coatings/coatingRT.m b/coatings/coatingRT.m
new file mode 100644
index 0000000..915e82b
--- /dev/null
+++ b/coatings/coatingRT.m
@@ -0,0 +1,229 @@
+% Coating reflectivity and transmittivity.

+%--------------------------------------------------------------------------

+% 

+% Front end code for "multidiel.m" by Orfanidis.

+%

+% Returns the reflectivity (fractional reflected power) and transmittivity

+% (fractional transmitted power) of a periodic, multilayer dielectric

+% coating as a function of the coating parameters and input beam properties

+% (incident wavelength, polarization, angle of incidence, etc.)

+%

+% To use coatingRT.m, the function multidiel.m by Orfanidis must be in your 

+% matlab path. The function multidiel.m is available on the mathworks website

+% in the package ewa.zip from Orfanidis:

+% http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=4456

+%

+% AUTHOR: Andri M. Gretarsson, 6/2007. 

+% LAST UPDATED: 7/2009 by AMG

+%

+% SYNTAX: [R,T,anginc] = coatingRT(nH,nL,N,lambda,lambda0,pol,...

+%             <nreflecs,Lprotect,nprotect,nb,anginc,LH,LL>)

+%

+% OR (if the first two arg's are vectors):

+%

+% SYNTAX: [R,T,anginc] = coatingRT(n,L,[],lambda,lambda0,pol,...

+%             <nreflecs,Lprotect,nprotect,nb,anginc>)

+%

+% where n is a vector of [incident medium index, the coating layer indices of refraction, substrate index]

+% and L is a vector of length length(n)-2 containing the coating layer OPTICAL thicknesses 

+% (in units of lambda0). In other words the arguments are in the same format as the first two argument

+% of multidiel.m

+%

+% INPUT VARIABLES

+% ---------------

+% nH        = index of refraction of the high-index coating layers (e.g. Ta2O5) OR vector

+%             of layer indices

+% nL        = index of refraction of the high-index coating layers (e.g. SiO2) OR vector

+%             of layer OPTICAL thicknesses

+% N         = number of coating layer pairs

+% lambda    = (1xN vector) incident wavelength in nanometers. This is the wavelength

+%             of the light that is actually hitting the coating.  Note: Units are

+%             nanometers.

+% lambda0   = center wavelength in nm.  This is the wavelength of

+%             light in terms of which the layer thicknesses are specified. For

+%             a high-reflective coating with quarter-wave layers, the center

+%             wavelength is usually also the wavelength of maximum reflectivity. (Hence

+%             the term "center wavelength.") Note: Units are nanometers.

+% pol       = polarization of the incident light, should be 'te' for

+%             s-polarization or 'tm' for p-polarization. ('te' stands for

+%             "transverse-electric." In other words, the electric field is

+%             transverse to the plane of incidence.) It's also OK to specify 'p' or

+%             's' as these are translated to 'tm' and 'te' respectively.

+% nreflecs  = This should usually be 0.  Setting this to a higher number

+%             causes the function to add in that number of _secondary_

+%             reflections from the back surface of the sample (assumed to be

+%             uncoated).

+% Lprotect  = OPTICAL thickness (in terms of lambda0) of any protective coating layer laid 

+%             down on top of the coating.  Usually, the topmost periodic layer

+%             is a low index layer (usually SiO2).  However, the topmost layer is 

+%             usually a diffeent thickness layer (e.g. 1/2 wavelength thick) rather than a 1/4 

+%             wavelength layer like the others. 

+%             This is accounted for in this code by setting Lprotect to 1/4,

+%             and nprotect=nL.  In other words the 1/2 wavelength topmost 

+%             layer is constructed as two 1/4 wavelength layers (the periodic

+%             layer and the protective layer) each with the same index of 

+%             refraction. And these are indeed the default values for Lprotect

+%             and nprotect.  Default: Lprotect=1/4.

+% nprotect  = Index of refraction of any additional coating layer laid down

+%             to protect the periodic coating.  Default: nprotect=nL.

+% nb        = Index of refraction of the substrate.  Default: nb=1.46 (silica).

+% anginc    = 1xM vector of incidence angles in radians. Default is 500

+%             evenly spaced angles between zero and pi/2.

+% LH        = OPTICAL thickness in terms of lambda0 of the high index layers if not 1/4.

+%             Default: LH=1/4.

+% LL        = OPTICAL thickness in terms of lambda0 of the low index layers if not 1/4.

+%             Default:  LL=1/4.

+%

+% OUTPUT VARIABLES

+% ----------------

+% R, T      = MxN matrix of Reflectivity and Transmissivity (fractional reflected or 

+%             transmitted _power_).  R and T are vectors of size(anginc).

+% lambda_mesh = This is the lambda mesh needed to plot R or T as a surface

+% anginc_mesh = This is the anginc mesh needed to plot R or T as a surface

+% r = The amplitude reflectance  (R = r x r*).

+%

+%

+% EXAMPLE 1  (R as a function of incident angle and wavelength)

+% ---------

+% lambda  = [199:3:1201]*1e-9;

+% lambda0 = 1064e-9;

+% anginc  = [0:0.5:90]*pi/180;

+% [R,T,lambda_mesh,anginc_mesh,r] = coatingRT(2.02,1.46,15,lambda,lambda0,'tm',...

+%                                 0,0,0,1.46,anginc);

+% h=pcolor(lambda_mesh*1e9,anginc_mesh*180/pi,R); 

+% set(h,'edgecolor','none');

+% colorbar;

+% caxis([0 1]);

+%

+%

+% EXAMPLE 2  (R as a function of incident angle at lambda=900 nm)

+% ---------

+% lambda  = 900e-9;

+% lambda0 = 1064e-9;

+% anginc  = [0:0.1:90]*pi/180;

+% [R,T,lambda,anginc,r] = coatingRT(2.02,1.46,15,lambda,lambda0,'tm',...

+%                                 0,0,0,1.46,anginc);

+% plot(anginc*180/pi,R); 

+%

+%

+% EXAMPLE 3 (R as a function of wavelength at an incident angle of 70 degrees)

+% ---------

+% lambda_in  = [500:0.1:1300]*1e-9;

+% lambda0 = 1064e-9;

+% anginc_in  = 70*pi/180;

+% [Rp,Tp,lambda,anginc] = coatingRT(2.02,1.46,15,lambda_in,lambda0,'tm',...

+%                                 0,0,0,1.46,anginc_in);

+% [Rs,Ts,lambda,anginc,r] = coatingRT(2.02,1.46,15,lambda_in,lambda0,'te',...

+%                                 0,0,0,1.46,anginc_in);

+% plot(lambda*1e9,Rp,lambda*1e9,Rs); 

+% axis([500 1300 0 1.1]); 

+%

+%

+% EXAMPLE 4 (R as a function of wavelength at normal incidence with an succession of irregular coatings)

+% ---------

+% for s=1:25

+% n=[1,repmat([2.0+0.2*rand(1,1),1.5+0.2*rand(1,1)],1,15),1.46];

+% L=0.25+[0.1*rand(1,30)];

+% lambda_in  = [500:0.1:1300]*1e-9;

+% lambda0 = 900e-9;

+% anginc_in  = 0;

+% [Rp,Tp,lambda,anginc] = coatingRT(n,L,[],lambda_in,lambda0,'tm',5,[],[],[],anginc_in);

+% plot(lambda*1e9,Rp);

+% axis([500 1300 0 1.1]);

+% drawnow

+% pause(0.1)

+% end 

+%

+% % EXAMPLE 5  (R from a simple uncoated glass surface as a function of incident angle at lambda=632 nm)

+% ---------

+% lambda  = 632e-9;

+% lambda0 = 1064e-9;

+% anginc  = [0:0.1:90]*pi/180;

+% [R,T,lambda,anginc,r] = coatingRT([1 1.5],[],[],lambda,lambda0,'p',...

+%                                 [],[],[],[],anginc);

+% plot(anginc*180/pi,R); set(gca,'xlim',[0 90]); shg;

+%

+%--------------------------------------------------------------------------

+% SYNTAX: [R,T,lambda_mesh,anginc_mesh,r] = coatingRT(nH,nL,N,lambda,...

+%             lambda0,pol <,nreflecs,Lprotect,nprotect,nb,anginc,LH,LL>)

+%--------------------------------------------------------------------------

+

+

+function [R,T,lambda_mesh,anginc_mesh,r] = coatingRT(nH,nL,N,lambda,lambda0,pol,varargin)

+

+m=6; n=1;

+if nargin>=m+n, nreflecs=varargin{n}; else nreflecs=0; end; n=n+1;

+if nargin>=m+n, Lprotect=varargin{n}; else Lprotect=0; end; n=n+1;

+if nargin>=m+n, nprotect=varargin{n}; else nprotect=nL; end; n=n+1;

+if nargin>=m+n, nb=varargin{n}; else nb=1.46; end; n=n+1;

+if nargin>=m+n, anginc=varargin{n}; else anginc=linspace(0,90,500)*pi/180; end; n=n+1;

+if nargin>=m+n, LH=varargin{n}; else LH=1/4; end; n=n+1;    % OPTICAL thicknesses of the layers in units of lambda0,

+if nargin>=m+n, LL=varargin{n}; else LL=1/4; end; n=n+1;    % where lambda0 is the design wavelength in units of nm

+

+na = 1;         % refractive index of medium outside the sample (usually air or vacuum);

+

+if strcmp(pol,'p')|strcmp(pol,'P')

+    pol='tm';

+elseif strcmp(pol,'s')|strcmp(pol,'S')

+    pol='te';

+elseif strcmp(pol,'tm')|strcmp(pol,'te')

+else

+    error('Polarization string was not understood. Should be ''te'' or ''tm''.');

+end

+

+if length(nH)>1 

+    n_in=nH; 

+    n_out=reverse(n_in);

+    L_in=nL;

+    L_out=reverse(L_in);

+    na=n_in(1);

+    nb=n_in(end);

+else

+    coating=repmat([nL,nH], 1, N);

+    L_in = (coating==nL)*LL + (coating==nH)*LH; % lengths of the layers in order as seen entering sample

+    if Lprotect~=0;

+        L_in=[Lprotect,L_in];

+        coating=[nprotect,coating];

+    end

+    L_out= reverse(L_in); % lengths of the layers in order as seen exiting sample

+    n_in = [na,coating,nb]; % indices for the layers in order as seen entering sample

+    n_out = reverse(n_in); % indices for the layers in order as seen exiting sample

+end

+

+angincg=asin(na*sin(anginc)/nb);                %Snell's law

+r_in=zeros(length(anginc),length(lambda));

+r_out=zeros(length(anginc),length(lambda));

+R_in=zeros(length(anginc),length(lambda));

+R_out=zeros(length(anginc),length(lambda));

+R_g=zeros(length(anginc),1);                    %Fresnel Eqn's don't depend on lambda explicitly, so neither does R_g

+R=zeros(length(anginc),length(lambda));

+[r_gs r_gp]=fresnel(nb,na,angincg*180/pi);      % nb is the "left hand" medium in the call to fresnel because the beam is inside the sample

+R_gs=r_gs.*conj(r_gs);

+R_gp=r_gp.*conj(r_gp);

+

+

+for s=1:length(anginc)

+    r_in(s,:)=multidiel(n_in,L_in,lambda/lambda0,anginc(s)*180/pi,pol);

+    r_out(s,:)=multidiel(n_out,L_out,lambda/lambda0,angincg(s)*180/pi,pol);

+    R_in(s,:)=(r_in(s,:).*conj(r_in(s,:)));

+    R_out(s,:)=(r_out(s,:).*conj(r_out(s,:)));

+    if strcmp(pol,'te'), R_g(s)=R_gs(s); else R_g(s)=R_gp(s); end

+    if strcmp(pol,'te'), r_g(s)=r_gs(s); else r_g(s)=r_gp(s); end

+    if nreflecs~=0

+        k=R_out(s,:)*R_g(s);  

+        % next line is sum from 0 to nreflecs internally reflecting light beams

+        % NOTE: nreflecs = 2 is equivalent the having 3 beams in total reflected

+        % from the sample (primary reflection + 2 secondary reflections).

+        R(s,:)=R_in(s,:) + ( 1-R_in(s,:) ) * R_g(s) .* (1-k.^nreflecs)./( 1-k ) .* ( 1-R_out(s,:) );

+        if nargout==5, 

+            warning('The amplitude reflection coefficient, r, does not include the effect of internal bounces. (Not yet implemented).'); 

+        end

+        r(s,:)=r_in(s,:);

+    else

+        R(s,:)=R_in(s,:);

+        r(s,:)=r_in(s,:);

+    end

+end

+

+T=1-R;

+[lambda_mesh,anginc_mesh]=meshgrid(lambda,anginc);

diff --git a/coatings/coatingRTnew.m b/coatings/coatingRTnew.m
new file mode 100644
index 0000000..4ce14ef
--- /dev/null
+++ b/coatings/coatingRTnew.m
@@ -0,0 +1,225 @@
+% Coating reflectivity and transmittivity.

+%--------------------------------------------------------------------------

+% 

+% Front end code for "multidiel.m" by Orfanidis.

+%

+% Returns the reflectivity (fractional reflected power) and transmittivity

+% (fractional transmitted power) of a periodic, multilayer dielectric

+% coating as a function of the coating parameters and input beam properties

+% (incident wavelength, polarization, angle of incidence, etc.)

+%

+% To use coatingRT.m, the function multidiel.m by Orfanidis must be in your 

+% matlab path. The function multidiel.m is available on the mathworks website

+% in the package ewa.zip from Orfanidis:

+% http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=4456

+%

+% AUTHOR: Andri M. Gretarsson, 6/2007. 

+% LAST UPDATED: 7/2009 by AMG

+%

+% SYNTAX: [R,T,anginc] = coatingRT(nH,nL,N,lambda,lambda0,pol,...

+%             <nreflecs,Lprotect,nprotect,nb,anginc,LH,LL>)

+%

+% OR (if the first two arg's are vectors):

+%

+% SYNTAX: [R,T,anginc] = coatingRT(n,L,[],lambda,lambda0,pol,...

+%             <nreflecs,Lprotect,nprotect,nb,anginc>)

+%

+% where n is a vector of [incident medium index, the coating layer indices of refraction, substrate index]

+% and L is a vector of length length(n)-2 containing the coating layer OPTICAL thicknesses 

+% (in units of lambda0). In other words the arguments are in the same format as the first two argument

+% of multidiel.m

+%

+% INPUT VARIABLES

+% ---------------

+% nH        = index of refraction of the high-index coating layers (e.g. Ta2O5) OR vector

+%             of layer indices

+% nL        = index of refraction of the high-index coating layers (e.g. SiO2) OR vector

+%             of physical layer thicknesses

+% N         = number of coating layers

+% lambda    = (1xN vector) incident wavelength in nanometers. This is the wavelength

+%             of the light that is actually hitting the coating.  Note: Units

+%             must be the same as those specifying the layer thicknesses (e.g. either

+%             meters or nanometers but not a mixture of the two!)

+% pol       = polarization of the incident light, should be 'te' for

+%             s-polarization or 'tm' for p-polarization. ('te' stands for

+%             "transverse-electric." In other words, the electric field is

+%             transverse to the plane of incidence.) It's also OK to specify 'p' or

+%             's' as these are translated to 'tm' and 'te' respectively.

+% nreflecs  = This should usually be 0.  Setting this to a higher number

+%             causes the function to add in that number of _secondary_

+%             reflections from the back surface of the sample (assumed to be

+%             uncoated).

+% Lprotect  = OPTICAL thickness (in terms of lambda0) of any protective coating layer laid 

+%             down on top of the coating.  Usually, the topmost periodic layer

+%             is a low index layer (usually SiO2).  However, the topmost layer is 

+%             usually a diffeent thickness layer (e.g. 1/2 wavelength thick) rather than a 1/4 

+%             wavelength layer like the others. 

+%             This is accounted for in this code by setting Lprotect to 1/4,

+%             and nprotect=nL.  In other words the 1/2 wavelength topmost 

+%             layer is constructed as two 1/4 wavelength layers (the periodic

+%             layer and the protective layer) each with the same index of 

+%             refraction. And these are indeed the default values for Lprotect

+%             and nprotect.  Default: Lprotect=1/4.

+% nprotect  = Index of refraction of any additional coating layer laid down

+%             to protect the periodic coating.  Default: nprotect=nL.

+% nb        = Index of refraction of the substrate.  Default: nb=1.46 (silica).

+% anginc    = 1xM vector of incidence angles in radians. Default is 500

+%             evenly spaced angles between zero and pi/2.

+% LH        = Actual physical thicknesses of the high index layers. Default is

+%             1064e-9/4/2.2 nanometers. (Quarter wavelength at 1064 nm of tantala.)

+% LL        = Actual physical thicknesses of the low index layers.

+%             1064e-9/4/1.46 nanometers. (Quarter wavelength at 1064 nm of silica.)

+%

+% OUTPUT VARIABLES

+% ----------------

+% R, T      = MxN matrix of Reflectivity and Transmissivity (fractional reflected or 

+%             transmitted _power_).  R and T are vectors of size(anginc).

+% lambda_mesh = This is the lambda mesh needed to plot R or T as a surface

+% anginc_mesh = This is the anginc mesh needed to plot R or T as a surface

+% r = The amplitude reflectance  (R = r x r*).

+%

+%

+% EXAMPLE 1  (R as a function of incident angle and wavelength for a tantala/silica 15 layer coating)

+% ---------

+% lambda  = [199:3:1201]*1e-9;

+% anginc  = [0:0.5:90]*pi/180;

+% [R,T,lambda_mesh,anginc_mesh,r] = coatingRTnew(2.2,1.46,15,lambda,'tm',...

+%                                 0,0,0,1.46,anginc,1.21e-7,1.82e-7);

+% h=pcolor(lambda_mesh*1e9,anginc_mesh*180/pi,R); 

+% set(h,'edgecolor','none');

+% colorbar;

+% caxis([0 1]);

+%

+%

+% EXAMPLE 2  (R as a function of incident angle at lambda=900 nm)

+% ---------

+% lambda  = 900e-9;

+% lambda0 = 1064e-9;

+% anginc  = [0:0.1:90]*pi/180;

+% [R,T,lambda,anginc,r] = coatingRTnew(2.2,1.46,15,lambda,'tm',...

+%                                 0,0,0,1.46,anginc,1.21e-7,1.82e-7);

+% plot(anginc*180/pi,R); 

+%

+%

+% EXAMPLE 3 (R as a function of wavelength at an incident angle of 70 degrees)

+% ---------

+% lambda_in  = [500:0.1:1300]*1e-9;

+% anginc_in  = 70*pi/180;

+% [Rp,Tp,lambda,anginc] = coatingRTnew(2.2,1.46,15,lambda_in,'tm',...

+%                                 0,0,0,1.46,anginc_in,1.21e-7,1.82e-7);

+% [Rs,Ts,lambda,anginc] = coatingRTnew(2.2,1.46,15,lambda_in,'te',...

+%                                 0,0,0,1.46,anginc_in,1.21e-7,1.82e-7);

+% plot(lambda*1e9,Rp,lambda*1e9,Rs); 

+% axis([500 1300 0 1.1]); 

+%

+%

+% EXAMPLE 4 (R as a function of wavelength at normal incidence from an irregular coating)

+% ---------

+% for s=1:25

+% n=[1,repmat([2.0+0.2*rand(1,1),1.5+0.2*rand(1,1)],1,15),1.46];

+% L=repmat([1/2.2 1/1.46]*0.25*850e-9,1,15)+[0.1e-7*randn(1,30)];

+% anginc_in  = 0;

+% [Rp,Tp,lambda,anginc] = coatingRTnew(n,L,[],lambda_in,'tm',5,[],[],[],anginc_in);

+% plot(lambda*1e9,Rp);

+% axis([500 1300 0 1.1]);

+% drawnow

+% pause(0.1)

+% end 

+%

+% % EXAMPLE 5  (R from a simple uncoated glass surface as a function of incident angle at lambda=632 nm)

+% ---------

+% lambda  = 632e-9;

+% anginc  = [0:0.1:90]*pi/180;

+% [R,T,lambda,anginc,r] = coatingRTnew([1 1.5],[],[],lambda,'p',...

+%                                 [],[],[],[],anginc);

+% plot(anginc*180/pi,R); set(gca,'xlim',[0 90]); shg;

+%

+%--------------------------------------------------------------------------

+% SYNTAX: [R,T,lambda_mesh,anginc_mesh,r] = coatingRT(nH,nL,N,lambda,...

+%             pol <,nreflecs,Lprotect,nprotect,nb,anginc,LH,LL>)

+%--------------------------------------------------------------------------

+

+

+function [R,T,lambda_mesh,anginc_mesh,r] = coatingRTnew(nH,nL,N,lambda,pol,varargin)

+lambda0=1064e-9;

+m=5; n=1;

+if nargin>=m+n, nreflecs=varargin{n}; else nreflecs=0; end; n=n+1;

+if nargin>=m+n, Lprotect=varargin{n}; else Lprotect=0; end; n=n+1;

+if nargin>=m+n, nprotect=varargin{n}; else nprotect=nL; end; n=n+1;

+if nargin>=m+n, nb=varargin{n}; else nb=1.46; end; n=n+1;

+if nargin>=m+n, anginc=varargin{n}; else anginc=linspace(0,90,500)*pi/180; end; n=n+1;

+if nargin>=m+n, LH=varargin{n}; else LH=lambda0/4/2.2; end; n=n+1;     % Actual physical thicknesses of the high index layers.

+if nargin>=m+n, LL=varargin{n}; else LL=lambda0/4/1.46; end; n=n+1;    % Actual physical thicknesses of the low index layers.

+na = 1;         % refractive index of medium outside the sample (usually air or vacuum);

+

+

+if strcmp(pol,'p')|strcmp(pol,'P')

+    pol='tm';

+elseif strcmp(pol,'s')|strcmp(pol,'S')

+    pol='te';

+elseif strcmp(pol,'tm')|strcmp(pol,'te')

+else

+    error('Polarization string was not understood. Should be ''te'',''tm'',''p'', or ''s''.');

+end

+

+if length(nH)>1 

+    n_in=nH; 

+    n_out=reverse(n_in);

+    L_in=nL;

+    if length(n_in)>2

+        L_in=L_in.*n_in(2:end-1)/lambda0;    % convert actual thicknesses to optical thickensses in terms of lambda0;

+    end

+    L_out=reverse(L_in);

+    na=n_in(1);

+    nb=n_in(end);

+else

+    coating=repmat([nL,nH], 1, N);

+    L_in = (coating==nL)*LL + (coating==nH)*LH; % lengths of the layers in order as seen entering sample   

+    if Lprotect~=0;

+        L_in=[Lprotect,L_in];

+        coating=[nprotect,coating];

+    end

+    n_in = [na,coating,nb]; % indices for the layers in order as seen entering sample

+    n_out = reverse(n_in); % indices for the layers in order as seen exiting sample    

+    L_in=L_in.*n_in(2:end-1)/lambda0;    % convert actual thicknesses to optical thickensses in terms of lambda0;

+    

+    L_out= reverse(L_in); % lengths of the layers in order as seen exiting sample

+end

+

+angincg=asin(na*sin(anginc)/nb);                %Snell's law

+r_in=zeros(length(anginc),length(lambda));

+r_out=zeros(length(anginc),length(lambda));

+R_in=zeros(length(anginc),length(lambda));

+R_out=zeros(length(anginc),length(lambda));

+R_g=zeros(length(anginc),1);                    %Fresnel Eqn's don't depend on lambda explicitly, so neither does R_g

+R=zeros(length(anginc),length(lambda));

+[r_gs r_gp]=fresnel(nb,na,angincg*180/pi);      % nb is the "left hand" medium in the call to fresnel because the beam is inside the sample

+R_gs=r_gs.*conj(r_gs);

+R_gp=r_gp.*conj(r_gp);

+

+

+for s=1:length(anginc)

+    r_in(s,:)=multidiel(n_in,L_in,lambda/lambda0,anginc(s)*180/pi,pol);

+    r_out(s,:)=multidiel(n_out,L_out,lambda/lambda0,angincg(s)*180/pi,pol);

+    R_in(s,:)=(r_in(s,:).*conj(r_in(s,:)));

+    R_out(s,:)=(r_out(s,:).*conj(r_out(s,:)));

+    if strcmp(pol,'te'), R_g(s)=R_gs(s); else R_g(s)=R_gp(s); end

+    if strcmp(pol,'te'), r_g(s)=r_gs(s); else r_g(s)=r_gp(s); end

+    if nreflecs~=0

+        k=R_out(s,:)*R_g(s);  

+        % next line is sum from 0 to nreflecs internally reflecting light beams

+        % NOTE: nreflecs = 2 is equivalent the having 3 beams in total reflected

+        % from the sample (primary reflection + 2 secondary reflections).

+        R(s,:)=R_in(s,:) + ( 1-R_in(s,:) ) * R_g(s) .* (1-k.^nreflecs)./( 1-k ) .* ( 1-R_out(s,:) );

+        if nargout==5, 

+            warning('The amplitude reflection coefficient, r, does not include the effect of internal bounces. (Not yet implemented).'); 

+        end

+        r(s,:)=r_in(s,:);

+    else

+        R(s,:)=R_in(s,:);

+        r(s,:)=r_in(s,:);

+    end

+end

+

+T=1-R;

+[lambda_mesh,anginc_mesh]=meshgrid(lambda,anginc);

diff --git a/coatings/prerequisits.txt b/coatings/prerequisits.txt
new file mode 100644
index 0000000..1900f1b
--- /dev/null
+++ b/coatings/prerequisits.txt
@@ -0,0 +1,5 @@
+To use coatingRT.m, the function multidiel.m by Orfanidis must be in your matlab path.
+
+multidiel.m is available on the mathworks website in the package ewa.zip from Orfanidis:
+
+http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=4456
\ No newline at end of file