//****************************************************************************** //** SCATMECH: Polarized Light Scattering C++ Class Library //** //** File: firstdiffuse.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 "firstdiffuse.h" #include "askuser.h" #include "matrix3d.h" using namespace std; namespace SCATMECH { // // Constructor... // First_Diffuse_BRDF_Model:: First_Diffuse_BRDF_Model() { l_scat = 1.; stack.wash(); phase_function = new Isotropic_Phase_Function; } int First_Diffuse_BRDF_Model:: get_recalc() { return BRDF_Model::get_recalc()||phase_function->get_recalc(); } void First_Diffuse_BRDF_Model:: setup() { BRDF_Model::setup(); } JonesMatrix First_Diffuse_BRDF_Model:: jones() { SETUP(); if (type==REFLECTION) { double n = substrate.n(lambda); // Use the absorption coefficient (k) to determine the absorption rate... double k_absorb = 4*pi/lambda*substrate.k(lambda); double atten = (k_absorb)/(1./l_scat+k_absorb); // Angles of refraction internal to the material... double sin_thetai_internal = sin(thetai)/n; double cos_thetai_internal = sqrt(1.-sqr(sin_thetai_internal)); double sin_thetas_internal = sin(thetas)/n; double cos_thetas_internal = sqrt(1.-sqr(sin_thetas_internal)); // The incident and scattering khat vectors outside the material... Vector inKo(sin(thetai),0.,-cos(thetai)); Vector outKo(sin(thetas)*cos(phis),sin(thetas)*sin(phis),cos(thetas)); // The incident and scattering khat vectors inside the material... // (one of my compilers not have a COMPLEX sin(COMPLEX)!) Vector inKi(sin_thetai_internal,0.,-cos_thetai_internal); Vector 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... Vector cnormal = normal; Vector inSi = perpto(inKi,cnormal); Vector inPi = perpto(inKi,inSi); Vector outSi = perpto(outKi,cnormal); Vector outPi = perpto(outKi,outSi); Vector inSo = perpto(inKo,cnormal); Vector inPo = perpto(inKo,inSo); Vector outSo = perpto(outKo,cnormal); Vector outPo = perpto(outKo,outSo); // The transmittance for incident light... double field_power_i = sqrt(n*cos_thetai_internal/cos(thetai)); COMPLEX tsi = stack.ts12(thetai,lambda,vacuum,substrate)*field_power_i; COMPLEX tpi = stack.tp12(thetai,lambda,vacuum,substrate)*field_power_i; // The transmittance for scattered light... double field_power_s = sqrt(n*cos_thetas_internal/cos(thetas)); COMPLEX tss = stack.ts12(thetas,lambda,vacuum,substrate)*field_power_s; COMPLEX tps = stack.tp12(thetas,lambda,vacuum,substrate)*field_power_s; // Transmission transfer matrices... CMatrix ti = (tsi*(CMatrix)outer(inSi,inSo))+(tpi*(CMatrix)outer(inPi,inPo)); CMatrix ts = (tss*(CMatrix)outer(outSo,outSi))+(tps*(CMatrix)outer(outPo,outPi)); // The Jacobian converting internal solid angle to external solid angle... double jacobian = cos(thetas)/cos_thetas_internal/sqr(n); // The cosine of the internal, local scattering angle... double cos_scatterangle = inKi*outKi; // The probability of scattering that angle per solid angle, if the scattering were scalar... // For g=0, the scatter function will be Rayleigh-like... double scatter = phase_function->f(acos(cos_scatterangle)); // Normalize that probability, because Rayleigh scattering is not isotropic... //scatter *= 1.5; scatter /= 0.5*(1.+sqr(cos_scatterangle)); // The Rayleigh scattering function... CMatrix scatterlocal = (outer(outPi,outPi)+outer(outSi,outSi)); // The probability that a ray is scattered, and that that scattered ray reaches // the surface, is given by... double first = sqr(1.-atten)*cos_thetas_internal/(cos_thetai_internal+cos_thetas_internal); // The scatter matrix globally includes transmission to/from defect... CMatrix scatterglobal = ts*(scatterlocal*ti); COMPLEX pp = (CVector)outPo*(scatterglobal*(CVector)inPo); COMPLEX sp = (CVector)outPo*(scatterglobal*(CVector)inSo); COMPLEX ps = (CVector)outSo*(scatterglobal*(CVector)inPo); COMPLEX ss = (CVector)outSo*(scatterglobal*(CVector)inSo); // Collect common factors, including the cos(thetas) required for BRDF... COMPLEX common=scatter*first*jacobian/cos(thetas); // Return the whole value with factors common to all elements... return JonesMatrix(pp,ss,ps,sp)*sqrt(common); } else { // TRANSMISSION... error("No code for TRANSMISSION mode"); return JonesMatrix(0.,0.,0.,0.); } } double Henyey_Greenstein_Phase_Function:: f(double theta) { SETUP(); return (1.-sqr(g))/pow(1.+sqr(g)-2.*g*cos(theta),1.5)/4./pi; } double Double_Henyey_Greenstein_Phase_Function:: f(double theta) { SETUP(); double q1 = (1.-sqr(g))/pow(1.+sqr(g)-2.*g*cos(theta),1.5)/4./pi; double q2 = (1.-sqr(g))/pow(1.+sqr(g)+2.*g*cos(theta),1.5)/4./pi; return fabs((1+c)/2.*q1+(1-c)/2.*q2); } double Table_Phase_Function:: f(double theta) { SETUP(); return table.value(theta/deg); } void Register(const First_Diffuse_BRDF_Model* x) { static bool Models_Registered = false; if (!Models_Registered) { Models_Registered=true; Register_Model(First_Diffuse_BRDF_Model); Register((Phase_Function*)0); } } void Register(const Phase_Function* x) { static bool Models_Registered = false; if (!Models_Registered) { Models_Registered=true; Register_Model(Phase_Function); Register_Model(Henyey_Greenstein_Phase_Function); Register_Model(Double_Henyey_Greenstein_Phase_Function); Register_Model(Isotropic_Phase_Function); Register_Model(Rayleigh_Phase_Function); Register_Model(Table_Phase_Function); } } DEFINE_MODEL(First_Diffuse_BRDF_Model,BRDF_Model, "First_Diffuse_BRDF_Model", "Model for single scattering in a diffuse material under a smooth interface."); DEFINE_VIRTUAL_MODEL(Phase_Function,Model, "Phase_Function", "General phase function for volume scattering"); DEFINE_MODEL(Henyey_Greenstein_Phase_Function,Phase_Function, "Henyey_Greenstein_Phase_Function", "Henyey-Greenstein phase function"); DEFINE_MODEL(Double_Henyey_Greenstein_Phase_Function,Phase_Function, "Double_Henyey_Greenstein_Phase_Function", "Double Henyey-Greenstein phase function with retroreflection"); DEFINE_MODEL(Isotropic_Phase_Function,Phase_Function, "Isotropic_Phase_Function", "Scattering which is isotropic."); DEFINE_MODEL(Rayleigh_Phase_Function,Phase_Function, "Rayleigh_Phase_Function", "Phase function appropriate for Rayleigh scattering"); DEFINE_MODEL(Table_Phase_Function,Phase_Function, "Table_Phase_Function", "Phase function provided by a table of values versus angle"); DEFINE_PARAMETER(First_Diffuse_BRDF_Model,double,l_scat,"Scatter mean free path [um]","1"); DEFINE_PTRPARAMETER(First_Diffuse_BRDF_Model,Phase_Function_Ptr,phase_function,"Phase Function","Henyey_Greenstein_Phase_Function"); DEFINE_PARAMETER(First_Diffuse_BRDF_Model,dielectric_stack,stack,"Film stack on substrate",""); DEFINE_PARAMETER(Henyey_Greenstein_Phase_Function,double,g,"Asymmetry parameter (g)","0.01"); DEFINE_PARAMETER(Double_Henyey_Greenstein_Phase_Function,double,g,"Asymmetry parameter (g)","0.01"); DEFINE_PARAMETER(Double_Henyey_Greenstein_Phase_Function,double,c,"Parameter c","0.1"); DEFINE_PARAMETER(Table_Phase_Function,Table,table,"Tabulated phase function","1"); } // namespace SCATMECH