Computes the reflection/refraction problem at a multilayer interface for TM or TE polarizations The problem is defined by j=1..nNbSurf+1 regions separated by i=1..nNbSurf interfaces at z=Z_i (see Fig. F.3). The first interface is at z=0. Region 1 is z<0 and contains the incident wave (incident from z=-infty) The plane of incidence is by convention xOz. Parameters: - sPol: string sPol='TE' for TE or sPol='TM' for TM (default value) polarizations - nNbSurf: scalar integer [1x1] number of interfaces (i.e. there are nNbSurf+1 regions) see Sec. F.4.1 and Fig. F.3 for more details. - lambda: column vector [Lx1] wavelengths in NANOMETERS (nm) - Cepsilon: cell of nNbSurf+1 column vectors (same size as lambda): i.e.: {(nNbSurf+1) * [Lx1]} lambda-dependent epsilons in each region. - CL: cell of nNbSurf scalar real numbers: {nNbSurf * [1x1]} region thicknesses (L_j=Z_{j+1}-Z_j) in NANOMETERS (i.e. distance between two interfaces). Note that CL{1}=L_1 must be zero by convention (see Fig. F.3) - aideg: row vector [1xA] angles of incidence in DEGREES. - d: optional scalar real number [1x1] (d=0 if not supplied) distance in NANOMETERS from interface where local field EFs are estimated Returns: a stResMulti structure containing fields as listed below: Note that EM fields are defined as in Eq. F.11 (TM) and Eq. F.18 (TE) For TM waves, we have: ** The Fresnel coefficients, which can be used for general properties: - rP: matrix [LxA] Fresnel reflection coefficients (Eqs. F.45, F.79) for each wavelength and angle of incidence - tP: matrix [LxA] Fresnel transmission coefficients (Eqs. F.46, F.80) for each wavelength and angle of incidence ** The field amplitudes, which can be used to calculate the EM field at any points: - HiyOvH1y: Cell of nNbSurf+1 matrices, {(nNbsurf+1)*[LxA]} field amplitude of first wave in each region i characterized by the ratio \bar{H}_{iy}/H_{1y} (see Eq. F.74 for a definition of \bar{H}_{iy}) given for each wavelength and angle of incidence - HipyOvH1y: Cell of nNbSurf+1 matrices, {(nNbsurf+1)*[LxA]} same as HiyOvH1y but for second wave: \bar{H}'_{iy}/H_{1y} - EixOvE1: Cell of nNbSurf+1 matrices, {(nNbsurf+1)*[LxA]} \bar{E}_{ix}/E_1 where E_1=H_{1y}/(\epsilon_0 c \sqrt{\epsilon_1}) - EizOvE1: {(nNbsurf+1)*[LxA]} \bar{E}_{iz}/E_1 - EipxOvE1: {(nNbsurf+1)*[LxA]} \bar{E}'_{ix}/E_1 - EipzOvE1: {(nNbsurf+1)*[LxA]} \bar{E}'_{iz}/E_1 ** The local field intensity enhancement factors (LFIEF) at a distance d on either side of each of the nNbSurf interface. For interface i, "outside" means towards region i and inside means towards region i+1. - MoutPerp: {nNbSurf * [LxA]} perpendicular LFIEF, i.e. |E_z|^2/|E_1|^2 at a distance d on the "outside" of each interface - MoutPara: {nNbSurf * [LxA]} parallel LFIEF, i.e. |E_x|^2/|E_1|^2 at a distance d on the "outside" of each interface - MinPerp: {nNbSurf * [LxA]} perpendicular LFIEF, i.e. |E_z|^2/|E_1|^2 at a distance d on the "inside" of each interface - MinPara: {nNbSurf * [LxA]} parallel LFIEF, i.e. |E_x|^2/|E_1|^2 at a distance d on the "inside" of each interface For TE waves, the output is similar but the name are changed to reflect the different nature of the fiels: - rS: matrix [LxA] Fresnel reflection coefficients (Eqs. F.62, F.79) for each wavelength and angle of incidence - tS: matrix [LxA] Fresnel transmission coefficients (Eqs. F.63, F.80) for each wavelength and angle of incidence - EiyOvE1y: {(nNbsurf+1)*[LxA]} field amplitude of first wave in each region i characterized by the ratio \bar{E}_{iy}/E_{1y} (see Eq. F.74 for a definition of \bar{E}_{iy}) given for each wavelength and angle of incidence - EipyOvE1y: {(nNbsurf+1)*[LxA]} same as EiyOvE1y but for second wave: \bar{E}'_{iy}/E_{1y} - HixOvH1: {(nNbsurf+1)*[LxA]} \bar{H}_{ix}/H_1 where H_1=E_{1y}*(\epsilon_0 c \sqrt{\epsilon_1}) - HizOvH1: {(nNbsurf+1)*[LxA]} \bar{H}_{iz}/H_1 - HipxOvH1: {(nNbsurf+1)*[LxA]} \bar{H}'_{ix}/H_1 - HipzOvH1: {(nNbsurf+1)*[LxA]} \bar{H}'_{iz}/H_1 - MoutPerp: {nNbSurf * [LxA]} perpendicular LFIEF, is zero for TE waves - MoutPara: {nNbSurf * [LxA]} parallel LFIEF, i.e. |E_x|^2/|E_1|^2 at a distance d on the "outside" of each interface - MinPerp: {nNbSurf * [LxA]} perpendicular LFIEF is zero for TE waves - MinPara: {nNbSurf * [LxA]} parallel LFIEF, i.e. |E_x|^2/|E_1|^2 at a distance d on the "inside" of each interface %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% This file is part of the SPlaC v1.0 package (copyright 2008) Check the README file for further information
0001 0002 %%%%%%%%%%%%%%%%%%%%%%% 0003 % Check parameters 0004 %%%%%%%%%%%%%%%%%%%%%%% 0005 % optional parameters 0006 if nargin < 7, d = 0; end 0007 0008 if ~((strcmp(sPol,'TM')) || (strcmp(sPol,'TE'))) 0009 disp 'Warning (PlnMultiRef): Polarization not properly specified. Assuming TM polarization.'; 0010 sPol='TM'; 0011 end 0012 0013 if (CL{1}~=0) 0014 disp 'Warning (PlnMultiRef): L_1 must be zero. Check CL parameter. Changing to CL{1}=0'; 0015 % enforce L_1=0 (by convention) 0016 CL{1}=0; 0017 end 0018 0019 %%%%%%%%%%%%%%%%%%%%%%% 0020 % Computations 0021 %%%%%%%%%%%%%%%%%%%%%%% 0022 0023 nNbAngle=length(aideg); % nb of columns (angles), i.e. A 0024 nNbLambda=length(lambda); % nb of rows (wavelengths), i.e. L 0025 0026 k0=2*pi./lambda; % column [Lx1] 0027 0028 % angle of incidence in radian 0029 ai=aideg/180*pi; % row [1xA] 0030 0031 % calculate k_x (same for all regions), Eq. F.42 0032 k1=k0.*sqrt(Cepsilon{1}); % column [Lx1] 0033 kx=k1 * sin(ai); % matrix [LxA] (lambda in rows and angles in columns) 0034 0035 % calculates k_{iz} in each region (matrices) and place in a cell 0036 Ckiz=cell(1,nNbSurf+1); 0037 for ii=1:(nNbSurf+1) 0038 % each CKiz{ii} is a matrix [LxA] 0039 % Computed from Eq. F.71 0040 Ckiz{ii}=sqrt(repmat(Cepsilon{ii}.*k0.^2,1,nNbAngle)-kx.^2); 0041 % Note that k_{iz} is assumed non-zero 0042 % except possibly in last region i=nNbSurf+1 0043 % Check this and issue a warning if not the case. 0044 % Assume k_{iz} is zero if k_{iz}/k0 smaller than 1e-6 0045 if (ii<nNbSurf+1) 0046 indPb=find(abs(Ckiz{ii}./repmat(k0,1,nNbAngle))<1e-6); 0047 for jj=1:length(indPb) 0048 % issue warning for each problem 0049 disp (['Warning: kiz is close to zero in region ' ... 0050 num2str(ii) ' at wavelength ' num2str(lambda(mod(indPb(jj),nNbLambda))) ... 0051 ' nm and angle ' num2str(aideg(floor(indPb(jj)-1)/nNbLambda)) ' degrees.']); 0052 disp (['kiz = ' num2str(Ckiz{indPb(jj)}/k0(mod(indPb(jj),nNbLambda))) ' * k0']); 0053 end 0054 end 0055 end 0056 0057 0058 % initialize M matrix to identity 0059 M11=1; 0060 M21=0; 0061 M12=0; 0062 M22=1; 0063 0064 % loop through the interfaces to compute Mi (Eq. F.76) and M (Eq. F.78) 0065 Mi11=cell(1,nNbSurf); 0066 Mi21=cell(1,nNbSurf); 0067 Mi12=cell(1,nNbSurf); 0068 Mi22=cell(1,nNbSurf); 0069 for ii=1:nNbSurf 0070 % calculate Ki [LxA] 0071 if (strcmp(sPol,'TM')) % TM wave, see Eq. F.16 and Eq. F.72 0072 Ki=repmat(Cepsilon{ii}./Cepsilon{ii+1},1,nNbAngle).*Ckiz{ii+1}./Ckiz{ii}; 0073 else % TE wave, see Eq. F.22 0074 Ki=Ckiz{ii+1}./Ckiz{ii}; 0075 end 0076 0077 % calculate matrix Mi components (all are matrices [LxA]) 0078 % using Eq. F.76 0079 AuxMi=exp(i*CL{ii}*Ckiz{ii}); 0080 Mi11{ii}=0.5*(1+Ki)./ AuxMi; 0081 Mi21{ii}=0.5*(1-Ki).* AuxMi; 0082 Mi12{ii}=0.5*(1-Ki)./ AuxMi; 0083 Mi22{ii}=0.5*(1+Ki).* AuxMi; 0084 0085 % update M matrix (Eq. F.78) 0086 M11new=M11.*Mi11{ii} + M12.*Mi21{ii}; 0087 M21new=M21.*Mi11{ii} + M22.*Mi21{ii}; 0088 M12new=M11.*Mi12{ii} + M12.*Mi22{ii}; 0089 M22new=M21.*Mi12{ii} + M22.*Mi22{ii}; 0090 M11=M11new; 0091 M21=M21new; 0092 M12=M12new; 0093 M22=M22new; 0094 end 0095 0096 %%%%%%%%%%%%%%%%%%%%%%% 0097 % Outputs 0098 %%%%%%%%%%%%%%%%%%%%%%% 0099 0100 % calculate Fresnel coefficients 0101 % see Eq. F.79 and F.80 (for TM) 0102 % same expressions for TE but relate to electric field 0103 tF=1./M11; 0104 rF=M21.*tF; % Eq. F.80 (for TM) 0105 0106 % initialize local field intensity EFs at each interface 0107 MoutPerp=cell(1,nNbSurf); 0108 MoutPara=cell(1,nNbSurf); 0109 MinPerp=cell(1,nNbSurf); 0110 MinPara=cell(1,nNbSurf); 0111 0112 0113 % treat separately TM and TE from now on 0114 % All field outputs are cells of matrices [LxA] 0115 if (strcmp(sPol,'TM')) 0116 %%%%%%%%%%%%%%%%% 0117 % TM waves 0118 %%%%%%%%%%%%%%%%% 0119 % initialize field amplitudes in each region 0120 HiyOvH1y=cell(1,nNbSurf+1); 0121 HpiyOvH1y=cell(1,nNbSurf+1); 0122 EixOvE1=cell(1,nNbSurf+1); 0123 EpixOvE1=cell(1,nNbSurf+1); 0124 EizOvE1=cell(1,nNbSurf+1); 0125 EpizOvE1=cell(1,nNbSurf+1); 0126 0127 % Fields in last region (nNbSurf+1) 0128 HiyOvH1y{nNbSurf+1}=tF; % H_{iy} / H_{1y} 0129 HpiyOvH1y{nNbSurf+1}=0; % H'_{iy} / H_{1y} 0130 % use for reference (possibly complex) E field: E1= H_{1y}/|H_{1y}| |E_1| 0131 % E1=H_{1y}/ [(epsilon_0 c) sqrt{epsilon_1}] 0132 AuxE1=repmat(sqrt(Cepsilon{1})./k0./Cepsilon{nNbSurf+1},1,nNbAngle); 0133 EixOvE1{nNbSurf+1}=HiyOvH1y{nNbSurf+1}.* Ckiz{nNbSurf+1} .* AuxE1; 0134 EpixOvE1{nNbSurf+1}=0; 0135 EizOvE1{nNbSurf+1}=-HiyOvH1y{nNbSurf+1}.* ( (Cepsilon{1}./Cepsilon{nNbSurf+1}) * sin(ai) ); 0136 EpizOvE1{nNbSurf+1}=0; 0137 0138 % loop (downwards) to compute field amplitudes in other regions 0139 for ii=nNbSurf:-1:1 0140 HiyOvH1y{ii}=Mi11{ii}.*HiyOvH1y{ii+1}+Mi12{ii}.*HpiyOvH1y{ii+1}; 0141 HpiyOvH1y{ii}=Mi21{ii}.*HiyOvH1y{ii+1}+Mi22{ii}.*HpiyOvH1y{ii+1}; 0142 0143 AuxE1=repmat(sqrt(Cepsilon{1})./k0./Cepsilon{ii},1,nNbAngle); 0144 EixOvE1{ii}=HiyOvH1y{ii}.* Ckiz{ii} .* AuxE1; 0145 EpixOvE1{ii}=-HpiyOvH1y{ii}.* Ckiz{ii} .* AuxE1; 0146 AuxE2=(Cepsilon{1}./Cepsilon{ii}) * sin(ai) ; % matrix product 0147 EizOvE1{ii}=-HiyOvH1y{ii}.* AuxE2; 0148 EpizOvE1{ii}=-HpiyOvH1y{ii}.* AuxE2; 0149 end 0150 0151 % loop again to compute local field EFs on each of the nNbSurf 0152 % interfaces 0153 for ii=1:nNbSurf 0154 MoutPerp{ii}=abs( EizOvE1{ii}.*exp(i*(CL{ii}-d)*Ckiz{ii})+ ... 0155 EpizOvE1{ii}.*exp(-i*(CL{ii}-d)*Ckiz{ii}) ).^2; 0156 MoutPara{ii}=abs( EixOvE1{ii}.*exp(i*(CL{ii}-d)*Ckiz{ii})+ ... 0157 EpixOvE1{ii}.*exp(-i*(CL{ii}-d)*Ckiz{ii}) ).^2; 0158 0159 MinPerp{ii}=abs( EizOvE1{ii+1}.*exp(i*d*Ckiz{ii+1})+ ... 0160 EpizOvE1{ii+1}.*exp(-i*d*Ckiz{ii+1}) ).^2; 0161 MinPara{ii}=abs( EixOvE1{ii+1}.*exp(i*d*Ckiz{ii+1})+ ... 0162 EpixOvE1{ii+1}.*exp(-i*d*Ckiz{ii+1}) ).^2; 0163 end 0164 0165 0166 else 0167 %%%%%%%%%%%%%%%%% 0168 % TE waves 0169 %%%%%%%%%%%%%%%%% 0170 0171 % initialize field amplitudes in each region 0172 EiyOvE1y=cell(1,nNbSurf+1); 0173 EpiyOvE1y=cell(1,nNbSurf+1); 0174 HixOvH1=cell(1,nNbSurf+1); 0175 HpixOvH1=cell(1,nNbSurf+1); 0176 HizOvH1=cell(1,nNbSurf+1); 0177 HpizOvH1=cell(1,nNbSurf+1); 0178 0179 % Fields in last region (nNbSurf+1) 0180 EiyOvE1y{nNbSurf+1}=tF; % E_{iy} / E_{1y} 0181 EpiyOvE1y{nNbSurf+1}=0; % E'_{iy} / E_{1y} 0182 % use for reference (possibly complex) H field: H1= E_{1y}/|E_{1y}| |H_1| 0183 % H1=E_{1y}* [(epsilon_0 c) sqrt{epsilon_1}] 0184 AuxH1=repmat(1./(sqrt(Cepsilon{1}).*k0),1,nNbAngle); 0185 HixOvH1{nNbSurf+1}=-EiyOvE1y{nNbSurf+1}.* Ckiz{nNbSurf+1} .* AuxH1; 0186 HpixOvH1{nNbSurf+1}=0; 0187 AuxH2=repmat(sin(ai),nNbLambda,1); 0188 HizOvH1{nNbSurf+1}=EiyOvE1y{nNbSurf+1}.* AuxH2; 0189 HpizOvH1{nNbSurf+1}=0; 0190 0191 % loop (downwards) to compute field amplitudes in other regions 0192 for ii=nNbSurf:-1:1 0193 EiyOvE1y{ii}=Mi11{ii}.*EiyOvE1y{ii+1}+Mi12{ii}.*EpiyOvE1y{ii+1}; 0194 EpiyOvE1y{ii}=Mi21{ii}.*EiyOvE1y{ii+1}+Mi22{ii}.*EpiyOvE1y{ii+1}; 0195 0196 HixOvH1{ii}=-EiyOvE1y{ii}.* Ckiz{ii} .* AuxH1; 0197 HpixOvH1{ii}=EpiyOvE1y{ii}.* Ckiz{ii} .* AuxH1; 0198 HizOvH1{ii}=EiyOvE1y{ii}.* AuxH2; 0199 HpizOvH1{ii}=EpiyOvE1y{ii}.* AuxH2; 0200 end 0201 0202 % loop again to compute local field EFs on each of the nNbSurf 0203 % interfaces 0204 for ii=1:nNbSurf 0205 MoutPara{ii}=abs( EiyOvE1y{ii}.*exp(i*(CL{ii}-d)*Ckiz{ii})+ ... 0206 EpiyOvE1y{ii}.*exp(-i*(CL{ii}-d)*Ckiz{ii}) ).^2; 0207 0208 MinPara{ii}=abs( EiyOvE1y{ii+1}.*exp(i*d*Ckiz{ii+1})+ ... 0209 EpiyOvE1y{ii+1}.*exp(-i*d*Ckiz{ii+1}) ).^2; 0210 MoutPerp{ii}=zeros(nNbLambda,nNbAngle); 0211 MinPerp{ii}=zeros(nNbLambda,nNbAngle); 0212 end 0213 end 0214 0215 % Finally, 0216 % create structure for results 0217 0218 % results common to TM and TE waves 0219 stResMulti.sPol=sPol; 0220 stResMulti.nNbSurf=nNbSurf; 0221 stResMulti.lambda=lambda; 0222 stResMulti.Cepsilon=Cepsilon; 0223 stResMulti.CL=CL; 0224 stResMulti.aideg=aideg; 0225 stResMulti.d=d; 0226 0227 stResMulti.MoutPerp=MoutPerp; 0228 stResMulti.MoutPara=MoutPara; 0229 stResMulti.MinPerp=MinPerp; 0230 stResMulti.MinPara=MinPara; 0231 0232 0233 if (strcmp(sPol,'TM')) % TM waves 0234 stResMulti.rP=rF; 0235 stResMulti.tP=tF; 0236 stResMulti.HiyOvH1y=HiyOvH1y; 0237 stResMulti.HpiyOvH1y=HpiyOvH1y; 0238 stResMulti.EixOvE1=EixOvE1; 0239 stResMulti.EizOvE1=EizOvE1; 0240 stResMulti.EpixOvE1=EpixOvE1; 0241 stResMulti.EpizOvE1=EpizOvE1; 0242 else % TE waves 0243 stResMulti.rS=rF; 0244 stResMulti.tS=tF; 0245 stResMulti.EiyOvE1y=EiyOvE1y; 0246 stResMulti.EpiyOvE1y=EpiyOvE1y; 0247 stResMulti.HixOvH1=HixOvH1; 0248 stResMulti.HizOvH1=HizOvH1; 0249 stResMulti.HpixOvH1=HpixOvH1; 0250 stResMulti.HpizOvH1=HpizOvH1; 0251 end 0252 0253 0254 0255 0256 0257 0258