% 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,... % ) % % OR (if the first two arg's are vectors): % % SYNTAX: [R,T,anginc] = coatingRT(n,L,[],lambda,lambda0,pol,... % ) % % 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);