//****************************************************************************** //** SCATMECH: Polarized Light Scattering C++ Class Library //** //** File: sphdfct.cpp //** //** Thomas A. Germer //** Optical Technology Division, National Institute of Standards and Technology //** 100 Bureau Dr. Stop 8443; Gaithersburg, MD 20899-8443 //** Phone: (301) 975-2876; FAX: (301) 975-6991 //** Email: thomas.germer@nist.gov //** //** Version: 6.00 (February 2008) //** //****************************************************************************** #include "scatmech.h" #include "sphdfct.h" #include "vector3d.h" #include "matrix3d.h" #include "fresnel.h" #include "askuser.h" using namespace std; namespace SCATMECH { static COMPLEX ArcSin(const COMPLEX& a) { COMPLEX temp = sqrt(1. - sqr(a)) + COMPLEX(0,1)*a; return COMPLEX(arg(temp),-log(abs(temp))); } // // Routine for calculating jones matrix from arbitrary spherical defect... // JonesMatrix Rayleigh_Defect_BRDF_Model:: jonesDSC() { SETUP(); if (type==REFLECTION) { // This routine uses the methodology outlined in Thomas A. Germer, // "Angular dependence and polarization of out-of-plane optical // scattering from particulate contamination, subsurface defects, // and surface microroughness," Appl. Opt. 36, 8798--8805 (1997). // // The code follows the derivation, except that the transmission // coefficients are chosen to be appropriate for a film. // // Also, there was an error in the paper regarding the Eq. (8) // and Appendix A. // COMPLEX n = substrate.index(lambda); COMPLEX e = substrate.epsilon(lambda); COMPLEX e_part=defect.epsilon(lambda); dielectric_constant epsilon=substrate.epsilon(lambda); COMPLEX x = kr*n; double k=2.*pi/lambda; // Angles of refraction internal to the material... COMPLEX sin_thetai_internal = sin_internal_angle(epsilon,thetai); COMPLEX cos_thetai_internal = cos_internal_angle(epsilon,thetai); COMPLEX sin_thetas_internal = sin_internal_angle(epsilon,thetas); COMPLEX cos_thetas_internal = cos_internal_angle(epsilon,thetas); // The incident and scattering khat vectors outside the material... CVector inKo((COMPLEX)(sin(thetai)),(COMPLEX)(0.),(COMPLEX)(-cos(thetai))); CVector outKo((COMPLEX)(sin(thetas)*cos(phis)), (COMPLEX)(sin(thetas)*sin(phis)),(COMPLEX)(cos(thetas))); // The incident and scattering khat vectors inside the material... // (one of my compilers not have a COMPLEX sin(COMPLEX)!) CVector inKi(sin_thetai_internal,(COMPLEX)(0.),-cos_thetai_internal); CVector outKi(sin_thetas_internal*cos(phis), sin_thetas_internal*sin(phis), cos_thetas_internal); // The surface normal... Vector normal(0.,0.,1.); // The following definitions make {s,p,k} right handed... CVector cnormal = normal; CVector inSi = perpto(inKi,cnormal); CVector inPi = perpto(inKi,inSi); CVector outSi = perpto(outKi,cnormal); CVector outPi = perpto(outKi,outSi); CVector inSo = perpto(inKo,cnormal); CVector inPo = perpto(inKo,inSo); CVector outSo = perpto(outKo,cnormal); CVector outPo = perpto(outKo,outSo); // Transmission coefficients... const dielectric_constant vacuum(1.); // The following four lines were removed by TAG 31.11.2000: //COMPLEX tsi = stack.ts12(thetai,vacuum,epsilon); //COMPLEX tpi = stack.tp12(thetai,vacuum,epsilon); //COMPLEX tss = stack.ts21(thetas,epsilon,vacuum); //COMPLEX tps = stack.tp21(thetas,epsilon,vacuum); // The following four lines were added by TAG 03.04.2001: COMPLEX tsi = stack.ts12(thetai,lambda,vacuum,substrate); COMPLEX tpi = stack.tp12(thetai,lambda,vacuum,substrate); COMPLEX tss = stack.ts21(thetas,lambda,substrate,vacuum); COMPLEX tps = stack.tp21(thetas,lambda,substrate,vacuum); // Transmission transfer matrices... CMatrix ti = (tsi*outer(inSi,inSo))+(tpi*outer(inPi,inPo)); // Eq. 8 of the paper is corrected here... // Factor of n corrected...TAG 31 DEC 2002... CMatrix ts = cos(thetas)/cos_thetas_internal/n* ((tss*outer(outSo,outSi))+(tps*outer(outPo,outPi))); COMPLEX rayleigh = (e_part-e)/(e_part+2.*e)*cube(x); CMatrix scatterlocal = (outer(outPi,outPi)+outer(outSi,outSi))*rayleigh; // Loss and phase factors... COMPLEX I(0,1); // Error found in following two lines in SCATMECH Version 1.0 // COMPLEX alpha = exp(2.*I*kd*cos(thetai)); // COMPLEX beta = exp(2.*I*kd*cos(thetas)); // Replaced with following two lines (TAG 11 April 2000): COMPLEX alpha = exp(I*n*kd*cos_thetai_internal); COMPLEX beta = exp(I*n*kd*cos_thetas_internal); // The scatter matrix globally includes transmission to/from defect... CMatrix scatterglobal = ts*(scatterlocal*ti)*alpha*beta; COMPLEX pp = outPo*(scatterglobal*inPo); COMPLEX sp = outPo*(scatterglobal*inSo); COMPLEX ps = outSo*(scatterglobal*inPo); COMPLEX ss = outSo*(scatterglobal*inSo); COMPLEX common=1./(k*n); // Return the whole value with factors common to all elements... return JonesMatrix(pp,ss,ps,sp)*common; } else { // TRANSMISSION... double n = substrate.n(lambda); double e = substrate.e1(lambda); COMPLEX e_part=defect.epsilon(lambda); dielectric_constant epsilon=substrate.epsilon(lambda); double x = kr*n; double k = 2.*pi/lambda; double sin_thetas = sin(thetas); double cos_thetas = sqrt(1.-sqr(sin_thetas)); double sin_thetas_external = sin_thetas*n; COMPLEX thetas_external = ArcSin((COMPLEX)sin_thetas_external); COMPLEX n2 = e_part/e; COMPLEX polarizability = (n2-1.)/(n2+2.)*cube(x); // Angles of refraction internal to the material... double sin_thetai_internal = sin(thetai)/n; double cos_thetai_internal = sqrt(1.-sqr(sin_thetai_internal)); Vector inKi(sin_thetai_internal,0.,-cos_thetai_internal); // The surface normal... Vector normal(0.,0.,1.); Vector inSi = perpto(inKi,normal); Vector inPi = perpto(inKi,inSi); JonesMatrix scatter_direct; { Vector outKi(sin_thetas*cos(phis),sin_thetas*sin(phis),-cos_thetas); Vector outSi = perpto(outKi,normal); Vector outPi = cross(outKi,outSi); Vector perp = perpto(inKi,outKi); // Polarization unit vectors in local Vector pari = cross(inKi,perp); // coordinates. Vector pars = cross(outKi,perp); JonesMatrix matrixin = JonesMatrix(perp*inPi,pari*inSi,pari*inPi,perp*inSi); JonesMatrix matrixout = JonesMatrix(outPi*perp,outSi*pars,outSi*perp,outPi*pars); double scatterangle = acos(inKi*outKi); JonesMatrix scatter = JonesMatrix(polarizability, polarizability*cos(scatterangle), 0.,0.); scatter_direct=matrixout*scatter*matrixin; } JonesMatrix scatter_indirect; { Vector outKr(sin_thetas*cos(phis),sin_thetas*sin(phis),cos_thetas); Vector outSr = perpto(outKr,normal); Vector outPr = cross(outKr,outSr); Vector perp = perpto(inKi,outKr); // Polarization unit vectors in local Vector pari = cross(inKi,perp); // coordinates. Vector pars = cross(outKr,perp); JonesMatrix matrixin = JonesMatrix(perp*inPi,pari*inSi,pari*inPi,perp*inSi); JonesMatrix matrixout = JonesMatrix(outPr*perp,outSr*pars,outSr*perp,outPr*pars); double scatterangle = acos(inKi*outKr); JonesMatrix scatter = JonesMatrix(polarizability, polarizability*cos(scatterangle), 0.,0.); JonesMatrix r = stack.r21(thetas_external,lambda,substrate,vacuum); COMPLEX phase = exp(COMPLEX(0,2)*cos_thetas*kd*n); scatter_indirect=matrixout*scatter*matrixin*r*phase; } JonesMatrix ti = stack.t12(thetai,lambda,vacuum,substrate)*(COMPLEX)sqrt(n); JonesMatrix scatter = (scatter_direct+scatter_indirect)*ti/k/n; return scatter; } } // // Routines to carry out one-time calculations... // void Rayleigh_Defect_BRDF_Model:: setup() { // Call parent's setup()... Local_BRDF_Model::setup(); double k = 2.*pi/lambda; kd = k*distance; kr = k*radius; } DEFINE_MODEL(Rayleigh_Defect_BRDF_Model,Local_BRDF_Model, "Rayleigh_Defect_BRDF_Model", "Scattering by a subsurface defect in the Rayleigh limit."); DEFINE_PARAMETER(Rayleigh_Defect_BRDF_Model,dielectric_stack,stack,"Film stack on substrate",""); DEFINE_PARAMETER(Rayleigh_Defect_BRDF_Model,double,radius,"Defect radius [um]","0.001"); DEFINE_PARAMETER(Rayleigh_Defect_BRDF_Model,double,distance,"Distance from center to surface [um]","0"); DEFINE_PARAMETER(Rayleigh_Defect_BRDF_Model,dielectric_function,defect,"Defect","(1,0)"); } // namespace SCATMECH