//****************************************************************************** //** SCATMECH: Polarized Light Scattering C++ Class Library //** //** File: sphprt.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 "sphprt.h" #include "vector3d.h" #include "matrix3d.h" #include "fresnel.h" #include "askuser.h" #include "rayscat.h" using namespace std; namespace SCATMECH { JonesMatrix Double_Interaction_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, and uses a generalized spherical // defect scattering function in place of Equation 5 of the paper. // // The code also uses reflection coefficients appropriate for a // dielectric layer below the sphere instead of a single interface. // This routine must call any child's setup() routine... double k = 2.*pi/lambda; double kd = k*distance; // The incident and scattering khat vectors... Vector inK = -polar(1.,-thetai,0.); Vector outK = polar(1.,thetas,phis); // The surface normal and the reflection matrix... Vector normal(0.,0.,1.); Matrix reflect(1.,0.,0., 0.,1.,0., 0.,0.,-1.); // The reflected khats... Vector inKr(inK.x,inK.y,-inK.z); // reflect*inK; Vector outKr(outK.x,outK.y,-outK.z); // reflect*outK; // Normal incidence or scattering flags bool notnormali = fabs(thetai)>1.0E-10; bool notnormals = fabs(thetas)>1.0E-10; // The polarization vectors... // IMPORTANT... // The following definitions make {s,p,k} right handed... Vector inS = notnormali ? perpto(inK,normal) : Vector(0.,-1.,0.); Vector inP = -perpto(inS,inK); Vector outS = notnormals ? perpto(outK,normal) : Vector(0.,-1.,0.); Vector outP = -perpto(outS,outK); // The reflected light polarization vectors... Vector inSr = notnormali ? perpto(inKr,normal) : Vector(0.,-1.,0.); Vector inPr = -perpto(inSr,inKr); Vector outSr = notnormals ? perpto(outKr,normal) : Vector(0.,-1.,0.); Vector outPr = -perpto(outSr,outKr); // Basis vectors appropriate for the scattering plane frame of reference // "perp" vectors... Vector b1in = perpto(inK,outK), b1out=b1in; Vector b2in = perpto(inKr,outK), b2out=b2in; Vector b3in = perpto(inK,outKr), b3out=b3in; Vector b4in = perpto(inKr,outKr), b4out=b4in; // "par" vectors... Vector a1in = cross(b1in,inK), a1out = cross(b1out,outK); Vector a2in = cross(b2in,inKr), a2out = cross(b2out,outK); Vector a3in = cross(b3in,inK), a3out = cross(b3out,outKr); Vector a4in = cross(b4in,inKr), a4out = cross(b4out,outKr); // The reflection coefficients... JonesMatrix ri = stack.r12(thetai,lambda,vacuum_function,substrate); JonesMatrix rs = stack.r12(thetas,lambda,vacuum_function,substrate); JonesMatrix S1 = scatterer->jones(Euler*rotatez(inK),Euler*rotatez(outK)); JonesMatrix S2 = scatterer->jones(Euler*rotatez(inKr),Euler*rotatez(outK)); JonesMatrix S3 = scatterer->jones(Euler*rotatez(inK),Euler*rotatez(outKr)); JonesMatrix S4 = scatterer->jones(Euler*rotatez(inKr),Euler*rotatez(outKr)); // Phase factors between the different beams... COMPLEX I(0,1); COMPLEX alpha = exp(2.*I*kd*cos(thetai)); COMPLEX beta = exp(2.*I*kd*cos(thetas)); JonesMatrix RotIn1 = JonesMatrix( b1in*inP, a1in*inS, a1in*inP, b1in*inS ); JonesMatrix RotIn2 = JonesMatrix( b2in*inPr, a2in*inSr, a2in*inPr, b2in*inSr ); JonesMatrix RotIn3 = JonesMatrix( b3in*inP, a3in*inS, a3in*inP, b3in*inS ); JonesMatrix RotIn4 = JonesMatrix( b4in*inPr, a4in*inSr, a4in*inPr, b4in*inSr ); JonesMatrix RotOut1 = JonesMatrix(b1out*outP, a1out*outS, b1out*outS, a1out*outP ); JonesMatrix RotOut2 = JonesMatrix(b2out*outP, a2out*outS, b2out*outS, a2out*outP ); JonesMatrix RotOut3 = JonesMatrix(b3out*outPr,a3out*outSr,b3out*outSr,a3out*outPr); JonesMatrix RotOut4 = JonesMatrix(b4out*outPr,a4out*outSr,b4out*outSr,a4out*outPr); S1 = RotOut1*S1*RotIn1; S2 = RotOut2*S2*RotIn2*alpha*ri; S3 = rs*RotOut3*S3*RotIn3*beta; S4 = rs*RotOut4*S4*RotIn4*alpha*beta*ri; COMPLEX common=1./k; return (S1+S2+S3+S4)*common; } else { // TRANSMISSION... double n = substrate.n(lambda); double k=2.*pi/lambda; Vector normal(0,0,1); Vector inK = polar(1.,pi-thetai,0.); Vector inKr = inK; inKr.z = -inKr.z; Vector inS = perpto(inK,normal); Vector inP = cross(inK,inS); Vector inSr = inS; Vector inPr = cross(inKr,inSr); double sin_thetas_outside = sin(thetas)*n; // This model does not handle evanescent waves at the particle... if (sin_thetas_outside>1.) return JonesZero(); double thetas_outside = asin(sin_thetas_outside); Vector outK = polar(1.,pi-thetas_outside,phis); Vector outS = perpto(outK,normal); Vector outP = cross(outK,outS); JonesMatrix scatter_direct; { Vector perp = perpto(inK,outK); // Polarization unit vectors in local Vector pari = cross(inK,perp); // coordinates. Vector pars = cross(outK,perp); JonesMatrix matrixin = JonesMatrix(perp*inP,pari*inS,pari*inP,perp*inS); JonesMatrix matrixout = JonesMatrix(outP*perp,outS*pars,outS*perp,outP*pars); JonesMatrix scatter = scatterer->jones(Euler*rotatez(inK),Euler*rotatez(outK)); scatter_direct = matrixout*scatter*matrixin; } JonesMatrix scatter_indirect; { Vector perp = perpto(inKr,outK); // Polarization unit vectors in local Vector pari = cross(inKr,perp); // coordinates. Vector pars = cross(outK,perp); JonesMatrix matrixin = JonesMatrix(perp*inPr,pari*inSr,pari*inPr,perp*inSr); JonesMatrix matrixout = JonesMatrix(outP*perp,outS*pars,outS*perp,outP*pars); JonesMatrix scatter = scatterer->jones(Euler*rotatez(inKr),Euler*rotatez(outK)); COMPLEX phase = exp(COMPLEX(0,2)*cos(thetai)*k*distance); JonesMatrix r = stack.r12(thetai,lambda,vacuum,substrate); scatter_indirect = phase*matrixout*scatter*matrixin*r; } JonesMatrix t = stack.t12(thetas_outside,lambda,vacuum,substrate)* sqrt(n*cos(thetas)/cos(thetas_outside)); JonesMatrix scatter = t*(scatter_direct+scatter_indirect)/k; double jacobian = sqr(n)*cos(thetas)/cos(thetas_outside); return scatter*sqrt(jacobian); } } // // Routine to carry out once-only operations // void Double_Interaction_BRDF_Model:: setup() { // Call parent's setup()... Local_BRDF_Model::setup(); if (scatterer->get_lambda()!=lambda) scatterer->set_lambda(lambda); if (scatterer->get_medium().n(lambda)!=1) scatterer->set_medium(vacuum); Euler = Matrix(cos(alpha), -(cos(beta)*sin(alpha)), sin(alpha)*sin(beta), sin(alpha), cos(alpha)*cos(beta), -(cos(alpha)*sin(beta)), 0, sin(beta), cos(beta)); } Vector Double_Interaction_BRDF_Model:: rotatez(const Vector& v) const { Vector result; result.x = cos(rotation)*v.x+sin(rotation)*v.y; result.y = -sin(rotation)*v.x+cos(rotation)*v.y; result.z = v.z; return result; } DEFINE_MODEL(Double_Interaction_BRDF_Model,Local_BRDF_Model, "Double_Interaction_BRDF_Model", "The double-interaction model for a spherical particle above a surface."); DEFINE_PARAMETER(Double_Interaction_BRDF_Model,dielectric_stack,stack,"Film stack on substrate",""); DEFINE_PARAMETER(Double_Interaction_BRDF_Model,double,distance,"Distance from center to surface [um]","0.05"); DEFINE_PTRPARAMETER(Double_Interaction_BRDF_Model,Free_Space_Scatterer_Ptr,scatterer,"The scattering function","MieScatterer"); DEFINE_PARAMETER(Double_Interaction_BRDF_Model,double,alpha,"Rotation of particle about z-axis [rad]","0"); DEFINE_PARAMETER(Double_Interaction_BRDF_Model,double,beta,"Rotation of particle about y-axis [rad]","0"); } // namespace SCATMECH