//****************************************************************************** //** SCATMECH: Polarized Light Scattering C++ Class Library //** //** File: rcw.cpp //** //** Thomas A. Germer //** Sensor Science Division, National Institute of Standards and Technology //** 100 Bureau Dr. Stop 8443; Gaithersburg, MD 20899-8443 //** Phone: (301) 975-2876 //** Email: thomas.germer@nist.gov //** //****************************************************************************** #include "rcw.h" #include "matrixmath.h" using namespace std; // // Equation numbers refer to Moharam et al. J. Opt. Soc. Am. A 12, 1068 (1995) or // Moharam et al. J. Opt. Soc. Am. A 12, 1077 (1995) namespace SCATMECH { using namespace CMLIB; namespace { template const T& min(const T& a,const T& b) { return (a const T& max(const T& a,const T& b) { return (a>b) ? a : b; } static const COMPLEX cI(0,1); } CVector RCW_Model::GetEField(const JonesVector& inpol, const Vector& pos, bool incident) { SETUP(); CVector result(0.,0.,0.); if (type==0 && pos.z<0) return result; if (type==1 && pos.z>-totalthickness) return result; if (type==2 && pos.z>-totalthickness) return result; if (type==3 && pos.z<0) return result; Vector R = pos; if (type==1) R.z += totalthickness; if (type==2) R.z += totalthickness; if (type==0 && incident) { Vector kin(nI*k0*sin(thetai*deg),0.,-nI*k0*cos(thetai*deg)); Vector s(0.,1.,0.); Vector p = cross(kin,s)/(k0*nI); COMPLEX phase = exp(cI*COMPLEX(kin*R)); result = (inpol.S()*CVector(s)+inpol.P()*CVector(p))*phase; } if (type==2 && incident) { CVector kin(conj(nII)*k0*sin(thetai*deg),0.,conj(nII)*k0*cos(thetai*deg)); CVector s(0.,1.,0.); CVector p = cross(kin,s)/(k0*conj(nII)); COMPLEX phase = exp(cI*COMPLEX(kin*(CVector)R)); result = (inpol.S()*CVector(s)+inpol.P()*CVector(p))*phase; } for (int i=0; i-totalthickness) return result; if (type==2 && pos.z>-totalthickness) return result; if (type==3 && pos.z<0) return result; Vector R = pos; if (type==1) R.z += totalthickness; if (type==2) R.z += totalthickness; if (type==0 && incident) { Vector kin(nI*k0*sin(thetai*deg),0.,-nI*k0*cos(thetai*deg)); Vector p(0.,-1,0.); Vector s = cross(p,kin)/(k0*nI); COMPLEX phase = nI*exp(cI*COMPLEX(kin*R)); result = (inpol.S()*CVector(s)+inpol.P()*CVector(p))*phase; } if (type==2 && incident) { CVector kin(conj(nII)*k0*sin(thetai*deg),0.,conj(nII)*k0*cos(thetai*deg)); CVector p(0.,-1.,0.); CVector s = cross(p,kin)/(k0*conj(nII)); COMPLEX phase = conj(nII)*exp(cI*COMPLEX(kin*(CVector)R)); result = (inpol.S()*CVector(s)+inpol.P()*CVector(p))*phase; } for (int i=0; iget_lambda()!=lambda) grating->set_lambda(lambda); // Exchange some information with grating... period = grating->get_period(); totalthickness=0; for (i=0; iget_levels(); ++i) { totalthickness+=grating->get_thickness(i); } dielectric_function medium_i = grating->get_medium_i(); dielectric_function medium_t = grating->get_medium_t(); k0 = 2*pi/lambda; n = order; nmat = 2*n+1; if (medium_i.k(lambda)!=0.) error("medium_i must be non-absorbing"); nI = medium_i.n(lambda); eI = sqr(nI); // Note: Moharam, et al., use n-ik convention... nII = conj((COMPLEX)medium_t.index(lambda)); eII = sqr(nII); kxi.resize(nmat); Kx.resize(nmat); kIzi.resize(nmat); kIIzi.resize(nmat); YI.resize(nmat); YII.resize(nmat); ZI.resize(nmat); ZII.resize(nmat); V.resize(nmat); kvector.resize(nmat); r.resize(nmat); R.resize(nmat); if (type==1 || type==2 || type==3) { if (imag(nII)!=0) { error("medium_t must be non-absorbing for type == 1, 2, or 3"); } } // Incident wavevector, with coordinates aligned along grating. // Note that SCATMECH generally assumes the z axis points out // of the grating. rotation represents a counterclockwise rotation // of the grating, looking down at the surface. if (type==0 || type==1) { inckx = k0*nI*sin(thetai*deg)*cos(rotation*deg); incky = -k0*nI*sin(thetai*deg)*sin(rotation*deg); inckz = -k0*nI*cos(thetai*deg); // This program has a problem if kx of incident light is zero... if (fabs(inckx/k0/nI)<1E-7) { inckx = k0*nI*1E-7; inckz = -k0*sqrt(eI - sqr(inckx/k0) - sqr(incky/k0)); } } else { inckx = k0*real(nII)*sin(thetai*deg)*cos(rotation*deg); incky = -k0*real(nII)*sin(thetai*deg)*sin(rotation*deg); inckz = -k0*real(nII)*cos(thetai*deg); // This program has a problem if kx of incident light is zero... if (fabs(inckx/k0/real(nII))<1E-7) { inckx = k0*real(nII)*1E-7; inckz = -k0*sqrt(real(eII) - sqr(inckx/k0) - sqr(incky/k0)); } } // Eq. (51) ... // Moharam et al. assume z axis points into grating, so y(Moharam) = -y(SCATMECH) ky = -incky; // Kvector (magnitude) associated with grating... double Kgrating = 2.*pi/period; for (i=0; i0) aa = conj(aa); kIzi[i] = aa; COMPLEX bb = k0*sqrt(eII*mu - kk - sqr(ky/k0)); if (imag(bb)>0) bb = conj(bb); kIIzi[i] = bb; // Defined after Eq. (24)... YI[i] = kIzi[i]/k0/mu; YII[i] = kIIzi[i]/k0/mu; // Defined after Eq. (44)... ZI[i] = kIzi[i]/k0/sqr(nI); ZII[i] = kIIzi[i]/k0/sqr(nII); } // Min and max orders that propagate... min_order=order+1; max_order=-order-1; // Rotation matrix elements to rotate from sample coordinates to // incident wave coordinates... double cost = cos(rotation*deg); double sint = sin(rotation*deg); for (i=-order; i<=order; ++i) { int j=i+n; R[j] = JonesZero(); V[j] = Vector(0,0,0); double kx = inckx-i*Kgrating; double ky = incky; if (type==0 || type==3) { COMPLEX eps = medium_i.epsilon(lambda); COMPLEX kz = sqrt(sqr(k0)*eps-sqr(kx)-sqr(ky)); // Diffracted wavevector in incident wave coordinates... kvector[j] = CVector(kx*cost-ky*sint,kx*sint+ky*cost,kz); // If it is propagating... if (imag(kz)==0.) { V[j] = Vector(real(kvector[j].x),real(kvector[j].y),real(kvector[j].z))/(k0*nI); if (imax_order) max_order=i; } else { V[j] = Vector(0.,0.,0.); } } else { // type==1 || type==2 COMPLEX eps = medium_t.epsilon(lambda); COMPLEX kz = -sqrt(sqr(k0)*eps-sqr(kx)-sqr(ky)); // Diffracted wavevector in incident wave coordinates... kvector[j] = CVector(kx*cost-ky*sint,kx*sint+ky*cost,kz); // If it is propagating... if (imag(kz)==0) { V[j] = Vector(real(kvector[j].x),real(kvector[j].y),-real(kvector[j].z))/(k0*real(nII)); if (imax_order) max_order=i; } else { V[j] = Vector(0.,0.,0.); } } } if (rotation==0) { switch (type) { case 0: InPlaneReflection(false); break; case 1: InPlaneTransmission(false); break; case 2: InPlaneReflection(true); break; case 3: InPlaneTransmission(true); break; default: error("Invalid type (must be 0, 1, 2, or 3)"); } } else { switch (type) { case 0: if (grating->is_anisotropic()) { ConicalAnisoReflection(false); } else { ConicalReflection(false); } break; case 1: if (grating->is_anisotropic()) { ConicalAnisoTransmission(false); } else { ConicalTransmission(false); } break; case 2: if (grating->is_anisotropic()) { ConicalAnisoReflection(true); } else { ConicalReflection(true); } break; case 3: if (grating->is_anisotropic()) { ConicalAnisoTransmission(true); } else { ConicalTransmission(true); } break; default: error("Invalid type (must be 0, 1, 2, or 3)"); } } kxi.resize(0); Kx.resize(0); kIzi.resize(0); YI.resize(0); YII.resize(0); ZI.resize(0); ZII.resize(0); } //*********************************************************************** //** ** //** In-Plane Reflection ** //** ** //*********************************************************************** void RCW_Model::InPlaneReflection(bool backward) { int i,j,k; int nnmat = 2*nmat; vector epsinvx(2*nmat+1),epsy(2*nmat+1),epsz(2*nmat+1); vector muinvx(2*nmat+1),muy(2*nmat+1),muz(2*nmat+1); vector Ex(sqr(nmat)),Ey(sqr(nmat)),Ez(sqr(nmat)); vector Exinv(sqr(nmat)),Eyinv(sqr(nmat)),Ezinv(sqr(nmat)); vector Mx(sqr(nmat)),My(sqr(nmat)),Mz(sqr(nmat)); vector Mxinv(sqr(nmat)),Myinv(sqr(nmat)),Mzinv(sqr(nmat)); vector A(sqr(nmat)); vector B(sqr(nmat)); vector EB(sqr(nmat)); vector MA(sqr(nmat)); vector Ws(sqr(nmat)),Qs(nmat),Xs(nmat),Vs(sqr(nmat)); vector Wp(sqr(nmat)),Qp(nmat),Xp(nmat),Vp(sqr(nmat)); vector Fs(sqr(nnmat)),Fp(sqr(nnmat)); vector indexs(nnmat),indexp(nnmat); vector WXVXs(nnmat); vector WXVXp(nnmat); vector as(2*nmat*nmat),ap(2*nmat*nmat); vector fs(nmat*nmat),gs(nmat*nmat); vector fp(nmat*nmat),gp(nmat*nmat); for (i=0; iget_levels()-1,backward ? +1 : -1 ); for (forloopint::iterator layer = layerloop.begin(); layer!=layerloop.end(); ++layer) { string layerstr = format("(layer = %d)",(int)layer); double d = grating->get_thickness(layer); if (d<0.) error("Thickness of layer " + to_string((int)layer) + " = " + to_string(d) + " less than zero"); // // Get the Fourier components for this layer... // message(layerstr + "Creating E and M matrices"); for (i=-nmat; i<=nmat; ++i) { if (grating->is_anisotropic()) { epsinvx[i+nmat]=conj(grating->fourierx(i,layer,1)); epsy[i+nmat]=conj(grating->fouriery(i,layer,0)); epsz[i+nmat]=conj(grating->fourierz(i,layer,0)); } else { epsinvx[i+nmat]=conj(grating->fourierx(i,layer,1)); epsy[i+nmat]=conj(grating->fourierx(i,layer,0)); epsz[i+nmat]=epsy[i+nmat]; } if (grating->is_magnetic()) { muinvx[i+nmat]=conj(grating->fouriermux(i,layer,1)); muy[i+nmat]=conj(grating->fouriermuy(i,layer,0)); muz[i+nmat]=conj(grating->fouriermuz(i,layer,0)); } else { muinvx[i+nmat]= (i==0 ? 1.: 0); muy[i+nmat]= (i==0 ? 1.: 0); muz[i+nmat]= (i==0 ? 1.: 0); } } for (i=0; iis_magnetic()) { Inverse(Mx,nmat); Inverse(Myinv,nmat); Inverse(Mzinv,nmat); } // // Building matrix A... // for (i=0; i Gs(nnmat*nnmat),Gp(nnmat*nnmat); for (i=0; i RRs(nnmat,0.); vector RRp(nnmat,0.); RRs[n] = -1.; //RRs[n+nmat] = COMPLEX(0.,-cos(thetai*deg)*nI); RRs[n+nmat] = COMPLEX(0.,inckz/k0); RRp[n] = -1.; if (backward) { RRp[n+nmat] = COMPLEX(0.,inckz/k0/sqr(real(nII))); } else { RRp[n+nmat] = COMPLEX(0.,inckz/k0/sqr(nI)); } message("Solving for reflected fields"); { vector indexs(nnmat),indexp(nnmat); message("LU decomposition of Gs"); LUdecompose(Gs,nnmat,indexs); message("LU decomposition of Gs"); LUbacksubstitute(Gs,nnmat,indexs,RRs); message("LU decomposition of Gp"); LUdecompose(Gp,nnmat,indexp); message("LU decomposition of Gp"); LUbacksubstitute(Gp,nnmat,indexp,RRp); } message("Finishing"); for (i=-order; i<=order; ++i) { j= n+i; r[j].PS() = 0.; r[j].PP() = conj(RRp[j]); r[j].SP() = 0.; r[j].SS() = conj(RRs[j]); if (backward) { R[j]=MuellerMatrix(r[j])*fabs(real(kIIzi[j])/inckz); } else { R[j]=MuellerMatrix(r[j])*fabs(real(kIzi[j])/inckz); } } } //*********************************************************************** //** ** //** In-Plane Transmission ** //** ** //*********************************************************************** void RCW_Model::InPlaneTransmission(bool backward) { int i,j,k; int nnmat = 2*nmat; vector epsinvx(2*nmat+1),epsy(2*nmat+1),epsz(2*nmat+1); vector muinvx(2*nmat+1),muy(2*nmat+1),muz(2*nmat+1); vector Ex(sqr(nmat)),Ey(sqr(nmat)),Ez(sqr(nmat)); vector Exinv(sqr(nmat)),Eyinv(sqr(nmat)),Ezinv(sqr(nmat)); vector Mx(sqr(nmat)),My(sqr(nmat)),Mz(sqr(nmat)); vector Mxinv(sqr(nmat)),Myinv(sqr(nmat)),Mzinv(sqr(nmat)); vector A(sqr(nmat)); vector B(sqr(nmat)); vector EB(sqr(nmat)); vector MA(sqr(nmat)); vector Ws(sqr(nmat)),Qs(nmat),Xs(nmat),Vs(sqr(nmat)); vector Wp(sqr(nmat)),Qp(nmat),Xp(nmat),Vp(sqr(nmat)); vector Fs(sqr(nnmat)),Fp(sqr(nnmat)); vector indexs(nnmat),indexp(nnmat); vector WXVXs(nnmat); vector WXVXp(nnmat); vector as(2*nmat*nmat),ap(2*nmat*nmat),as_(2*nmat*nmat),ap_(2*nmat*nmat); vector fs(nmat*nmat),gs(nmat*nmat),ss(nmat*nmat),ts(nmat*nmat); vector fp(nmat*nmat),gp(nmat*nmat),sp(nmat*nmat),tp(nmat*nmat); double costhetai = backward ? -inckz/k0/real(nII) : -inckz/k0/nI; for (i=0; iget_levels()-1,backward ? -1 : +1 ); for (forloopint::iterator layer = layerloop.begin(); layer!=layerloop.end(); ++layer) { string layerstr = format("(layer = %d)",(int)layer); double d = grating->get_thickness(layer); if (d<0.) error("Thickness of layer " + to_string((int)layer) + " = " + to_string(d) + " less than zero"); // // Get the Fourier components for this layer... // message(layerstr + "Creating E and M matrices"); for (i=-nmat; i<=nmat; ++i) { if (grating->is_anisotropic()) { epsinvx[i+nmat]=conj(grating->fourierx(i,layer,1)); epsy[i+nmat]=conj(grating->fouriery(i,layer,0)); epsz[i+nmat]=conj(grating->fourierz(i,layer,0)); } else { epsinvx[i+nmat]=conj(grating->fourierx(i,layer,1)); epsy[i+nmat]=conj(grating->fourierx(i,layer,0)); epsz[i+nmat]=epsy[i+nmat]; } if (grating->is_magnetic()) { muinvx[i+nmat]=conj(grating->fouriermux(i,layer,1)); muy[i+nmat]=conj(grating->fouriermuy(i,layer,0)); muz[i+nmat]=conj(grating->fouriermuz(i,layer,0)); } else { muinvx[i+nmat]=(i==0 ? 1.:0.); muy[i+nmat]=(i==0 ? 1.:0.); muz[i+nmat]=(i==0 ? 1.:0.); } } for (i=0; iis_magnetic()) { Inverse(Mx,nmat); Inverse(Myinv,nmat); Inverse(Mzinv,nmat); } // // Building matrix A... // for (i=0; i Gs(nnmat*nnmat),Gp(nnmat*nnmat); for (i=0; i Ts(nnmat,0.); vector Tp(nnmat,0.); for (i=0; i indexs(nnmat),indexp(nnmat); LUdecompose(Gs,nnmat,indexs); LUbacksubstitute(Gs,nnmat,indexs,Ts); LUdecompose(Gp,nnmat,indexp); LUbacksubstitute(Gp,nnmat,indexp,Tp); } message("Finishing"); for (i=-order; i<=order; ++i) { j= n+i; if (backward) { r[j].PS() = 0; r[j].PP() = conj(Tp[j+nmat]/nI*nII); r[j].SP() = 0.; r[j].SS() = conj(Ts[j+nmat]); R[j]=MuellerMatrix(r[j])*fabs(real(kIzi[j])/inckz); } else { r[j].PP() = conj(Tp[j+nmat]/nII*nI); r[j].SP() = 0.; r[j].SS() = conj(Ts[j+nmat]); R[j]=MuellerMatrix(r[j])*fabs(real(kIIzi[j])/inckz); } } } //*********************************************************************** //** ** //** Conical Reflection for Isotropic Media ** //** ** //*********************************************************************** void RCW_Model::ConicalReflection(bool backward) { int i,j,k; int nnmat = 4*nmat; int nmat2 = 2*nmat; int nmat4 = 4*nmat; vector EEE(2*nmat+1),EEEE(2*nmat+1); vector E(sqr(nmat)), Einv(sqr(nmat)); vector EE(sqr(nmat)), EEinv(sqr(nmat)); vector MMM(2*nmat+1),MMMM(2*nmat+1); vector M(sqr(nmat)), Minv(sqr(nmat)); vector MM(sqr(nmat)), MMinv(sqr(nmat)); vector A(sqr(nmat)), Ainv(sqr(nmat)); vector B(sqr(nmat)), Binv(sqr(nmat)); vector Matrix1(sqr(nmat)),Matrix2(sqr(nmat)); vector W1(sqr(nmat)),Q1(nmat),X1(nmat); vector W2(sqr(nmat)),Q2(nmat),X2(nmat); vector V11(sqr(nmat)),V12(sqr(nmat)),V21(sqr(nmat)),V22(sqr(nmat)),v21(sqr(nmat)),v12(sqr(nmat)); vector Vsp(sqr(nmat)),Wsp(sqr(nmat)),Wpp(sqr(nmat)),Vpp(sqr(nmat)), Vss(sqr(nmat)),Wss(sqr(nmat)),Wps(sqr(nmat)),Vps(sqr(nmat)); vector BigMat(sqr(nnmat)),BigMatLU(sqr(nnmat)); vector index(nnmat); vector WXVXl(nnmat),WXVXr(nnmat),WXVXlx(nnmat),WXVXrx(nnmat); vector ab(nmat2*nmat4), fg(nmat2*nmat4); vector Fc(nmat),Fs(nmat); // // Set the initial values for the fg matrix... // for (i=0; iget_levels()-1,backward ? +1 : -1 ); for (forloopint::iterator layer = layerloop.begin(); layer!=layerloop.end(); ++layer) { string layerstr = format("(layer = %d)",(int)layer); double d = grating->get_thickness(layer); if (d<0.) error("Thickness of layer " + to_string((int)layer) + " = " + to_string(d) + " less than zero"); // // Get the Fourier components for this layer... // for (i=-nmat; i<=nmat; ++i) { EEE[i+nmat]=conj(grating->fourierx(i,layer,0)); EEEE[i+nmat]=conj(grating->fourierx(i,layer,1)); MMM[i+nmat]=conj(grating->fouriermux(i,layer,0)); MMMM[i+nmat]=conj(grating->fouriermux(i,layer,1)); } for (i=0; iis_magnetic()) { Inverse(Minv,nmat); Inverse(MMinv,nmat); } // // Building matrix A... // message(layerstr + "Building and inverting A matrix"); for (i=0; i RRs(nnmat,0.); vector RRp(nnmat,0.); if (backward) { double costhetai = -inckz/k0/real(nII); RRs[ n ] = -1.; RRs[ n+nmat ] = -cI*nII*costhetai; RRp[n+2*nmat] = -cI*nII; RRp[n+3*nmat] = costhetai; } else { double costhetai = -inckz/k0/nI; RRs[ n ] = -1.; RRs[ n+nmat ] = -cI*nI*costhetai; RRp[n+2*nmat] = -cI*nI; RRp[n+3*nmat] = costhetai; } message("Solving for fields"); LUbacksubstitute(BigMat,nnmat,index,RRs); LUbacksubstitute(BigMat,nnmat,index,RRp); message("Finishing"); for (i=-order; i<=order; ++i) { j= n+i; // Conventions for incident polarizations are handled above. // Conventions for reflected polarizations need factor of -1 for S // and sqrt(-1) for P. Conjugates are needed for conversion to -iwt from iwt. if (backward) { r[j].PS() = conj( RRp[j] ); r[j].PP() = conj( cI*RRp[j+nmat]/nII ); r[j].SP() = conj( -cI*RRs[j+nmat]/nII ); r[j].SS() = conj( -RRs[j] ); R[j]=MuellerMatrix(r[j])*fabs(real(kIIzi[j])/inckz); } else { r[j].PS() = conj( -RRp[j] ); r[j].PP() = conj( cI*RRp[j+nmat]/nI ); r[j].SP() = conj( cI*RRs[j+nmat]/nI ); r[j].SS() = conj( -RRs[j] ); R[j]=MuellerMatrix(r[j])*fabs(real(kIzi[j])/inckz); } } } //*********************************************************************** //** ** //** Conical Transmission for Isotropic Media ** //** ** //*********************************************************************** void RCW_Model::ConicalTransmission(bool backward) { int i,j,k; int nnmat = 4*nmat; int nmat2 = 2*nmat; int nmat4 = 4*nmat; vector EEE(2*nmat+1),EEEE(2*nmat+1); vector E(sqr(nmat)), Einv(sqr(nmat)); vector EE(sqr(nmat)), EEinv(sqr(nmat)); vector MMM(2*nmat+1),MMMM(2*nmat+1); vector M(sqr(nmat)), Minv(sqr(nmat)); vector MM(sqr(nmat)), MMinv(sqr(nmat)); vector A(sqr(nmat)), Ainv(sqr(nmat)); vector B(sqr(nmat)), Binv(sqr(nmat)); vector Matrix1(sqr(nmat)),Matrix2(sqr(nmat)); vector W1(sqr(nmat)),Q1(nmat),X1(nmat); vector W2(sqr(nmat)),Q2(nmat),X2(nmat); vector V11(sqr(nmat)),V12(sqr(nmat)),V21(sqr(nmat)),V22(sqr(nmat)),v21(sqr(nmat)),v12(sqr(nmat)); vector Vsp(sqr(nmat)),Wsp(sqr(nmat)),Wpp(sqr(nmat)),Vpp(sqr(nmat)), Vss(sqr(nmat)),Wss(sqr(nmat)),Wps(sqr(nmat)),Vps(sqr(nmat)); vector BigMat(sqr(nnmat)),BigMat2(sqr(nnmat)); vector index(nnmat); vector WXVX(nnmat); vector ab(nmat4*nmat4), fg(nmat2*nmat4),st(nmat2*nmat4); vector Fc(nmat),Fs(nmat); // // Set the initial values for the fg matrix... // for (i=0; iget_levels()-1,backward ? -1 : +1 ); for (forloopint::iterator layer = layerloop.begin(); layer!=layerloop.end(); ++layer) { string layerstr = format("(layer = %d)",(int)layer); double d = grating->get_thickness(layer); if (d<0.) error("Thickness of layer " + to_string((int)layer) + " = " + to_string(d) + " less than zero"); // // Get the Fourier components for this layer... // for (i=-nmat; i<=nmat; ++i) { EEE[i+nmat]=conj(grating->fourierx(i,layer,0)); EEEE[i+nmat]=conj(grating->fourierx(i,layer,1)); } if (grating->is_magnetic()) { for (i=-nmat; i<=nmat; ++i) { MMM[i+nmat]=conj(grating->fouriermux(i,layer,0)); MMMM[i+nmat]=conj(grating->fouriermux(i,layer,1)); } } else { for (i=-nmat; i<=nmat; ++i) { MMM[i+nmat]= (i==0 ? 1. : 0.); MMMM[i+nmat]= (i==0 ? 1. : 0.); } } for (i=0; iis_magnetic()) { Inverse(Minv,nmat); Inverse(MMinv,nmat); } // // Building matrix A... // message(layerstr + "Building and inverting A matrix"); for (i=0; i Ts(nnmat); vector Tp(nnmat); for (i=0; i EEEx(2*nmat+1),EEEEx(2*nmat+1); vector Ex(sqr(nmat)), Exinv(sqr(nmat)); vector EEx(sqr(nmat)), EExinv(sqr(nmat)); vector EEEy(2*nmat+1),EEEEy(2*nmat+1); vector Ey(sqr(nmat)), Eyinv(sqr(nmat)); vector EEy(sqr(nmat)), EEyinv(sqr(nmat)); vector EEEz(2*nmat+1),EEEEz(2*nmat+1); vector Ez(sqr(nmat)), Ezinv(sqr(nmat)); vector EEz(sqr(nmat)), EEzinv(sqr(nmat)); vector MMMx(2*nmat+1),MMMMx(2*nmat+1); vector Mx(sqr(nmat)), Mxinv(sqr(nmat)); vector MMx(sqr(nmat)), MMxinv(sqr(nmat)); vector MMMy(2*nmat+1),MMMMy(2*nmat+1); vector My(sqr(nmat)), Myinv(sqr(nmat)); vector MMy(sqr(nmat)), MMyinv(sqr(nmat)); vector MMMz(2*nmat+1),MMMMz(2*nmat+1); vector Mz(sqr(nmat)), Mzinv(sqr(nmat)); vector MMz(sqr(nmat)), MMzinv(sqr(nmat)); vector F(sqr(nmat2)), G(sqr(nmat2)); vector FG(sqr(nmat2)); vector W(sqr(nmat2)),Q(nmat2),v(sqr(nmat2)); vector X(nmat2); vector Ws(nmat*nmat2); vector Wp(nmat*nmat2); vector Vs(nmat*nmat2); vector Vp(nmat*nmat2); vector X1(nmat); vector X2(nmat); vector Vs1(sqr(nmat)),Vs2(sqr(nmat)),Vp1(sqr(nmat)),Vp2(sqr(nmat)), Ws1(sqr(nmat)),Ws2(sqr(nmat)),Wp1(sqr(nmat)),Wp2(sqr(nmat)); vector BigMat(sqr(nnmat)),BigMatLU(sqr(nnmat)); vector index(nnmat); vector WXVXl(nnmat),WXVXr(nnmat),WXVXlx(nnmat),WXVXrx(nnmat); vector ab(nmat2*nmat4), fg(nmat2*nmat4); vector Fc(nmat),Fs(nmat); // // Set the initial values for the fg matrix... // for (i=0; iget_levels()-1,backward ? +1 : -1 ); for (forloopint::iterator layer = layerloop.begin(); layer!=layerloop.end(); ++layer) { string layerstr = format("(layer = %d)",(int)layer); double d = grating->get_thickness(layer); if (d<0.) error("Thickness of layer " + to_string((int)layer) + " = " + to_string(d) + " less than zero"); // // Get the Fourier components for this layer... // for (i=-nmat; i<=nmat; ++i) { if (grating->is_anisotropic()) { EEEx[i+nmat]=conj(grating->fourierx(i,layer,0)); EEEEx[i+nmat]=conj(grating->fourierx(i,layer,1)); EEEy[i+nmat]=conj(grating->fouriery(i,layer,0)); EEEEy[i+nmat]=conj(grating->fouriery(i,layer,1)); EEEz[i+nmat]=conj(grating->fourierz(i,layer,0)); EEEEz[i+nmat]=conj(grating->fourierz(i,layer,1)); } else { EEEx[i+nmat]=conj(grating->fourierx(i,layer,0)); EEEEx[i+nmat]=conj(grating->fourierx(i,layer,1)); EEEy[i+nmat]=EEEx[i+nmat]; EEEEy[i+nmat]=EEEEx[i+nmat]; EEEz[i+nmat]=EEEx[i+nmat]; EEEEz[i+nmat]=EEEEx[i+nmat]; } if (grating->is_magnetic()) { MMMx[i+nmat]=conj(grating->fouriermux(i,layer,0)); MMMMx[i+nmat]=conj(grating->fouriermux(i,layer,1)); MMMy[i+nmat]=conj(grating->fouriermuy(i,layer,0)); MMMMy[i+nmat]=conj(grating->fouriermuy(i,layer,1)); MMMz[i+nmat]=conj(grating->fouriermuz(i,layer,0)); MMMMz[i+nmat]=conj(grating->fouriermuz(i,layer,1)); } else { MMMx[i+nmat]=(i==0 ? 1.: 0.); MMMMx[i+nmat]=(i==0 ? 1.: 0.); MMMy[i+nmat]=(i==0 ? 1.: 0.); MMMMy[i+nmat]=(i==0 ? 1.: 0.); MMMz[i+nmat]=(i==0 ? 1.: 0.); MMMMz[i+nmat]=(i==0 ? 1.: 0.); } } for (i=0; iis_magnetic()) { Inverse(MMxinv,nmat); Inverse(Myinv,nmat); Inverse(Mzinv,nmat); } // // Building matrix R and S (the upper right and lower left off-diagonals of Eq. (57) of M&G)... // message(layerstr + "Building matrices F and G "); for (i=0; i RRs(nnmat,0.); vector RRp(nnmat,0.); if (backward) { double costhetai = -inckz/k0/real(nII); RRs[ n ] = -1.; RRs[ n+nmat ] = -cI*nII*costhetai; RRp[n+2*nmat] = -cI*nII; RRp[n+3*nmat] = costhetai; } else { double costhetai = -inckz/k0/nI; RRs[ n ] = -1.; RRs[ n+nmat ] = -cI*nI*costhetai; RRp[n+2*nmat] = -cI*nI; RRp[n+3*nmat] = costhetai; } message("Solving for fields"); LUbacksubstitute(BigMat,nnmat,index,RRs); LUbacksubstitute(BigMat,nnmat,index,RRp); message("Finishing"); for (i=-order; i<=order; ++i) { j= n+i; // Conventions for incident polarizations are handled above. // Conventions for reflected polarizations need factor of -1 for S // and sqrt(-1) for P. Conjugates are needed for conversion to -iwt from iwt. if (backward) { r[j].PS() = conj( RRp[j] ); r[j].PP() = conj( cI*RRp[j+nmat]/nII ); r[j].SP() = conj(-cI*RRs[j+nmat]/nII ); r[j].SS() = conj( -RRs[j] ); R[j]=MuellerMatrix(r[j])*fabs(real(kIIzi[j])/inckz); } else { r[j].PS() = conj( -RRp[j] ); r[j].PP() = conj( cI*RRp[j+nmat]/nI ); r[j].SP() = conj( cI*RRs[j+nmat]/nI ); r[j].SS() = conj( -RRs[j] ); R[j]=MuellerMatrix(r[j])*fabs(real(kIzi[j])/inckz); } } } //*********************************************************************** //** ** //** Conical Transmission for Anisotropic Media ** //** ** //*********************************************************************** void RCW_Model::ConicalAnisoTransmission(bool backward) { int i,j,k; int nnmat = 4*nmat; int nmat2 = 2*nmat; int nmat4 = 4*nmat; vector EEEx(2*nmat+1),EEEEx(2*nmat+1); vector Ex(sqr(nmat)), Exinv(sqr(nmat)); vector EEx(sqr(nmat)), EExinv(sqr(nmat)); vector EEEy(2*nmat+1),EEEEy(2*nmat+1); vector Ey(sqr(nmat)), Eyinv(sqr(nmat)); vector EEy(sqr(nmat)), EEyinv(sqr(nmat)); vector EEEz(2*nmat+1),EEEEz(2*nmat+1); vector Ez(sqr(nmat)), Ezinv(sqr(nmat)); vector EEz(sqr(nmat)), EEzinv(sqr(nmat)); vector MMMx(2*nmat+1),MMMMx(2*nmat+1); vector Mx(sqr(nmat)), Mxinv(sqr(nmat)); vector MMx(sqr(nmat)), MMxinv(sqr(nmat)); vector MMMy(2*nmat+1),MMMMy(2*nmat+1); vector My(sqr(nmat)), Myinv(sqr(nmat)); vector MMy(sqr(nmat)), MMyinv(sqr(nmat)); vector MMMz(2*nmat+1),MMMMz(2*nmat+1); vector Mz(sqr(nmat)), Mzinv(sqr(nmat)); vector MMz(sqr(nmat)), MMzinv(sqr(nmat)); vector F(sqr(nmat2)), G(sqr(nmat2)); vector FG(sqr(nmat2)); vector W(sqr(nmat2)),Q(nmat2),v(sqr(nmat2)); vector X(nmat2); vector Ws(nmat*nmat2); vector Wp(nmat*nmat2); vector Vs(nmat*nmat2); vector Vp(nmat*nmat2); vector X1(nmat); vector X2(nmat); vector Vs1(sqr(nmat)),Vs2(sqr(nmat)),Vp1(sqr(nmat)),Vp2(sqr(nmat)), Ws1(sqr(nmat)),Ws2(sqr(nmat)),Wp1(sqr(nmat)),Wp2(sqr(nmat)); //vector Matrix1(sqr(nmat)),Matrix2(sqr(nmat)); vector BigMat(sqr(nnmat)),BigMat2(sqr(nnmat)); vector index(nnmat); vector WXVX(nnmat); vector ab(nmat4*nmat4), fg(nmat2*nmat4),st(nmat2*nmat4); vector Fc(nmat),Fs(nmat); // // Set the initial values for the fg matrix... // for (i=0; iget_levels()-1,backward ? -1 : +1 ); for (forloopint::iterator layer = layerloop.begin(); layer!=layerloop.end(); ++layer) { string layerstr = format("(layer = %d)",(int)layer); double d = grating->get_thickness(layer); if (d<0.) error("Thickness of layer " + to_string((int)layer) + " = " + to_string(d) + " less than zero"); // // Get the Fourier components for this layer... // for (i=-nmat; i<=nmat; ++i) { if (grating->is_anisotropic()) { EEEx[i+nmat]=conj(grating->fourierx(i,layer,0)); EEEEx[i+nmat]=conj(grating->fourierx(i,layer,1)); EEEy[i+nmat]=conj(grating->fouriery(i,layer,0)); EEEEy[i+nmat]=conj(grating->fouriery(i,layer,1)); EEEz[i+nmat]=conj(grating->fourierz(i,layer,0)); EEEEz[i+nmat]=conj(grating->fourierz(i,layer,1)); } else { EEEx[i+nmat]=conj(grating->fourierx(i,layer,0)); EEEEx[i+nmat]=conj(grating->fourierx(i,layer,1)); EEEy[i+nmat]=EEEx[i+nmat]; EEEEy[i+nmat]=EEEEx[i+nmat]; EEEz[i+nmat]=EEEx[i+nmat]; EEEEz[i+nmat]=EEEEx[i+nmat]; } if (grating->is_magnetic()) { MMMx[i+nmat]=conj(grating->fouriermux(i,layer,0)); MMMMx[i+nmat]=conj(grating->fouriermux(i,layer,1)); MMMy[i+nmat]=conj(grating->fouriermuy(i,layer,0)); MMMMy[i+nmat]=conj(grating->fouriermuy(i,layer,1)); MMMz[i+nmat]=conj(grating->fouriermuz(i,layer,0)); MMMMz[i+nmat]=conj(grating->fouriermuz(i,layer,1)); } else { MMMx[i+nmat]=(i==0 ? 1.: 0.); MMMMx[i+nmat]=(i==0 ? 1.: 0.); MMMy[i+nmat]=(i==0 ? 1.: 0.); MMMMy[i+nmat]=(i==0 ? 1.: 0.); MMMz[i+nmat]=(i==0 ? 1.: 0.); MMMMz[i+nmat]=(i==0 ? 1.: 0.); } } for (i=0; iis_magnetic()) { Inverse(MMxinv,nmat); Inverse(Myinv,nmat); Inverse(Mzinv,nmat); } // // Building matrix R and S (the upper right and lower left off-diagonals of Eq. (57) of M&G)... // message(layerstr + "Building matrices F and G "); for (i=0; i Ts(nnmat); vector Tp(nnmat); for (i=0; iget_recalc()) { grating->SETUP(); RCW.set_grating(grating); } if (order<0) error("order < 0"); if (grating->get_medium_t().index(lambda)!=substrate.index(lambda)) error("grating->medium_t != substrate"); if (order!=RCW.get_order()) RCW.set_order(order); if (type!=RCW.get_type()) RCW.set_type(type); if (lambda!=RCW.get_lambda()) RCW.set_lambda(lambda); if (thetai/deg!=RCW.get_thetai()) RCW.set_thetai(thetai/deg); if (rotation/deg!=RCW.get_rotation()) RCW.set_rotation(rotation/deg); V.resize(2*order+1); R.resize(2*order+1); for (int i=-order; i<=order; ++i) { int j=i+order; V[j]=RCW.GetDirection(i); R[j]=RCW.GetIntensity(i); } cosalpha = cos(alpha); Omega = pi*sqr(alpha); low = RCW.GetMinimumPropagatingOrder(); high = RCW.GetMaximumPropagatingOrder(); } MuellerMatrix RCW_BRDF_Model::mueller() { SETUP(); Vector v = polar(1.,thetas,phis); MuellerMatrix result = MuellerZero(); for (int i=low; i<=high; ++i) { if (V[i+order]*v>=cosalpha) { result += R[i+order]/(Omega*v.z); } } return result; } DEFINE_MODEL(RCW_Model,Model,"Rigorous coupled wave theory for a grating"); DEFINE_PARAMETER(RCW_Model,int,order,"Order of calculation","25",0xFF); DEFINE_PARAMETER(RCW_Model,int,type,"(0) for Forward/Reflection, (1) for Forward/Transmission, (2) for Backward/Reflection, or (3) for Backward/Transmission","0",0xFF); DEFINE_PARAMETER(RCW_Model,double,lambda,"Wavelength (vacuum)","0.532",0xFF); DEFINE_PARAMETER(RCW_Model,double,thetai,"Incident angle [deg]","0",0xFF); DEFINE_PARAMETER(RCW_Model,double,rotation,"Azimuthal rotation of sample [deg]","0",0xFF); DEFINE_PTRPARAMETER(RCW_Model,Grating_Ptr,grating,"Grating","Single_Line_Grating",0xFF); DEFINE_MODEL(RCW_BRDF_Model,BRDF_Model,"Rigorous coupled wave theory for a grating, form fitted into a BRDF_Model"); DEFINE_PARAMETER(RCW_BRDF_Model,double,alpha,"Half angle of diffraction cone [rad]","0.0175",0x01); DEFINE_PARAMETER(RCW_BRDF_Model,int,order,"Maximum order","25",0x02); DEFINE_PTRPARAMETER(RCW_BRDF_Model,Grating_Ptr,grating,"Grating","Single_Line_Grating",0x02); }