//****************************************************************************** //** SCATMECH: Polarized Light Scattering C++ Class Library //** //** File: subsphere.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 "subsphere.h" #include "askuser.h" #include "rayscat.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))); } void Subsurface_Particle_BRDF_Model:: setup() { Local_BRDF_Model::setup(); if (substrate.k(lambda)!=0) error("Substrate cannot be absorbing"); if (scatterer->get_lambda()!=lambda) scatterer->set_lambda(lambda); if (scatterer->get_medium().n(lambda)!=substrate.n(lambda)) scatterer->set_medium(substrate); 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 Subsurface_Particle_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; } JonesMatrix Subsurface_Particle_BRDF_Model:: jonesDSC() { SETUP(); double n = substrate.n(lambda); double k=2.*pi/lambda; // 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); if (type==REFLECTION) { double sin_thetas_internal = sin(thetas)/n; double cos_thetas_internal = sqrt(1.-sqr(sin_thetas_internal)); Vector outKi(sin_thetas_internal*cos(phis), sin_thetas_internal*sin(phis), cos_thetas_internal); // The following definitions make {s,p,k} right handed... Vector outSi = perpto(outKi,normal); Vector outPi = perpto(outKi,outSi); JonesMatrix ti = stack.t12(thetai,lambda,vacuum,substrate)*(COMPLEX)sqrt(n*cos_thetai_internal/cos(thetai)); JonesMatrix ts = stack.t12(thetas,lambda,vacuum,substrate)*(COMPLEX)sqrt(n*cos_thetas_internal/cos(thetas)); 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 = scatterer->jones(Euler*rotatez(inKi),Euler*rotatez(outKi)); scatter=matrixout*scatter*matrixin; double j1 = cos(thetas)/cos_thetas_internal/sqr(n); double j2 = cos(thetai)/cos_thetai_internal; // The scatter matrix globally includes transmission to/from defect... scatter = ts*scatter*ti; COMPLEX common=j1*j2/sqr(k*n); // Return the whole value with factors common to all elements... return scatter*sqrt(common); } else { // TRANSMISSION... 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); 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 = scatterer->jones(Euler*rotatez(inKi),Euler*rotatez(outKi)); 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 = scatterer->jones(Euler*rotatez(inKi),Euler*rotatez(outKr)); JonesMatrix r = stack.r21(thetas_external,lambda,substrate,vacuum); COMPLEX phase = exp(COMPLEX(0,2)*cos_thetas*k*depth*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; } } DEFINE_MODEL(Subsurface_Particle_BRDF_Model,Local_BRDF_Model, "Subsurface_Particle_BRDF_Model", "Scattering by a particle beneath an interface of a nonabsorbing media"); DEFINE_PARAMETER(Subsurface_Particle_BRDF_Model,dielectric_stack,stack,"Film stack on substrate",""); DEFINE_PARAMETER(Subsurface_Particle_BRDF_Model,double,depth,"Depth of center of particle [um]","0"); DEFINE_PTRPARAMETER(Subsurface_Particle_BRDF_Model,Free_Space_Scatterer_Ptr,scatterer,"The scattering function","MieScatterer"); DEFINE_PARAMETER(Subsurface_Particle_BRDF_Model,double,alpha,"Rotation of particle about z-axis [rad]","0"); DEFINE_PARAMETER(Subsurface_Particle_BRDF_Model,double,beta,"Rotation of particle about y-axis [rad]","0"); } //namespace SCATMECH