diff options
Diffstat (limited to 'coatings/coatingRT.m')
| -rw-r--r-- | coatings/coatingRT.m | 229 |
1 files changed, 229 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);
|