From a854b0970dc4b75e453e8c3c85a4070802c6c5c3 Mon Sep 17 00:00:00 2001 From: germer Date: Fri, 29 Sep 2017 08:33:02 -0400 Subject: [PATCH] Added StackModel and Legendre_Phase_Function, and a few other bug fixes --- code/Makefile | 1 - code/allrough.cpp | 15 +- code/allrough.h | 10 +- code/askuser.cpp | 4 +- code/axipart.h | 2 +- code/axipart1.cpp | 18 +- code/axipart2.cpp | 16 +- code/bobvlieg.h | 6 +- code/bobvlieg1.cpp | 40 +- code/bobvlieg2.cpp | 43 +- code/brdf.cpp | 91 ++- code/crossgrating.cpp | 4 +- code/crossrcw.cpp | 12 + code/crossrcw.h | 5 +- code/crough.cpp | 7 +- code/dielfunc.cpp | 6 +- code/dielfunc.h | 3 + code/diffuse.cpp | 9 +- code/diffuse.h | 2 +- code/facet.cpp | 18 +- code/facet.h | 2 +- code/filmtran.cpp | 1370 ++++++++++++++++++----------------- code/filmtran.h | 683 +++++++++++------ code/firstdiffuse.cpp | 51 +- code/firstdiffuse.h | 16 +- code/flake.cpp | 14 +- code/flake.h | 2 +- code/grating.cpp | 119 +-- code/grating.h | 12 +- code/inherit.h | 2 +- code/jmatrix.cpp | 8 +- code/matrixmath.h | 62 +- code/matrixmath2.cpp | 4 +- code/models.cpp | 2 + code/mueller.cpp | 37 +- code/mueller.h | 6 +- code/nsphere.cpp | 10 +- code/nsphere.h | 2 +- code/rayinst.cpp | 10 +- code/rayinst.h | 2 +- code/raystack.cpp | 54 +- code/raystack.h | 2 +- code/rcw.cpp | 2 +- code/rcw.h | 2 +- code/roughnes.cpp | 12 +- code/roughnes.h | 2 +- code/scateval.cpp | 13 +- code/scattabl.cpp | 20 +- code/sphdfct.cpp | 34 +- code/sphdfct.h | 2 +- code/sphprt.cpp | 18 +- code/sphprt.h | 2 +- code/subbobvlieg.cpp | 39 +- code/subbobvlieg.h | 6 +- code/subsphere.cpp | 20 +- code/subsphere.h | 2 +- code/transmit.cpp | 10 +- code/transmit.h | 2 +- code/twoface.cpp | 7 +- code/zernikeexpansion.cpp | 180 +++++ code/zernikeexpansion.h | 57 ++ docs/axipart.htm | 7 +- docs/bobvlieg.htm | 19 +- docs/classes.htm | 2 + docs/correlated_stack.htm | 13 +- docs/diffuse.htm | 9 +- docs/download.htm | 6 +- docs/facet.htm | 14 +- docs/firstdiffuse.htm | 10 +- docs/flake.htm | 25 +- docs/grating.htm | 34 +- docs/growth_stack.htm | 11 +- docs/history.htm | 19 + docs/multilayermie.htm | 4 +- docs/phase.htm | 50 ++ docs/raydef.htm | 10 +- docs/rayleighinstrument.htm | 10 +- docs/raystack.htm | 8 +- docs/roughnes.htm | 10 +- docs/sphrpart.htm | 10 +- docs/stackmodel.htm | 231 ++++++ docs/stackmodel.htm~ | 223 ++++++ docs/subaxipart.htm | 7 +- docs/subbobvlieg.htm | 18 +- docs/subsphere.htm | 8 +- docs/torrspar.htm | 10 +- docs/transmit.htm | 7 +- docs/uncorrelated_stack.htm | 10 +- 88 files changed, 2603 insertions(+), 1394 deletions(-) create mode 100644 code/zernikeexpansion.cpp create mode 100644 code/zernikeexpansion.h create mode 100644 docs/stackmodel.htm create mode 100644 docs/stackmodel.htm~ diff --git a/code/Makefile b/code/Makefile index e8c5ed9..6c23fa9 100644 --- a/code/Makefile +++ b/code/Makefile @@ -9,7 +9,6 @@ #//** Phone: (301) 975-2876 #//** Email: thomas.germer@nist.gov #//** -#//** Version: 7.00 (January 2015) #//** #//****************************************************************************** diff --git a/code/allrough.cpp b/code/allrough.cpp index bd8f994..1ec69a9 100644 --- a/code/allrough.cpp +++ b/code/allrough.cpp @@ -38,7 +38,8 @@ namespace SCATMECH { SETUP(); JonesMatrix J=JonesZero(); - for (int layer=0; layer<=stack.get_n(); ++layer) { + + for (int layer=0; layer<=stack->get_n(); ++layer) { model.set_this_layer(layer); J = J + model.Jones(thetai,thetas,phis,rotation); } @@ -65,7 +66,7 @@ namespace SCATMECH { SETUP(); MuellerMatrix M=MuellerZero(); - for (int layer=0; layer<=stack.get_n(); ++layer) { + for (int layer=0; layer<=stack->get_n(); ++layer) { model.set_this_layer(layer); M = M + model.Mueller(thetai,thetas,phis,rotation); } @@ -104,7 +105,7 @@ namespace SCATMECH { { SETUP(); - int N = stack.get_n(); + int N = stack->get_n(); double fx,fy; Bragg_Frequency(fx,fy); double qpow = pow(4*sqr(pi)*(sqr(fx)+sqr(fy)),exponent/2.); @@ -118,7 +119,7 @@ namespace SCATMECH { vector a(N); vector PSDint(N); for (int i=0; iget_t()[i]; // The replication factor between the layers... a[i] = exp(-argument); // The intrinsic roughness of the layer... @@ -164,16 +165,16 @@ namespace SCATMECH { DEFINE_MODEL(Correlated_Roughness_Stack_BRDF_Model,Roughness_BRDF_Model, "All films in a stack equally rough and correlated"); - DEFINE_PARAMETER(Correlated_Roughness_Stack_BRDF_Model,dielectric_stack,stack,"Film stack on substrate","",0xFF); + DEFINE_PTRPARAMETER(Correlated_Roughness_Stack_BRDF_Model,StackModel_Ptr,stack,"Film stack on substrate","No_StackModel",0xFF); DEFINE_MODEL(Uncorrelated_Roughness_Stack_BRDF_Model,Roughness_BRDF_Model, "All films in a stack equally rough but uncorrelated"); - DEFINE_PARAMETER(Uncorrelated_Roughness_Stack_BRDF_Model,dielectric_stack,stack,"Film stack on substrate","",0xFF); + DEFINE_PTRPARAMETER(Uncorrelated_Roughness_Stack_BRDF_Model,StackModel_Ptr,stack,"Film stack on substrate","No_StackModel",0xFF); DEFINE_MODEL(Growth_Roughness_Stack_BRDF_Model,Roughness_BRDF_Model, "All films in a stack having intrinsic roughness, but also correlation with roughness from the stack below."); - DEFINE_PARAMETER(Growth_Roughness_Stack_BRDF_Model,dielectric_stack,stack,"Film stack on substrate","",0xFF); + DEFINE_PTRPARAMETER(Growth_Roughness_Stack_BRDF_Model,StackModel_Ptr,stack,"Film stack on substrate","No_StackModel",0xFF); DEFINE_PARAMETER(Growth_Roughness_Stack_BRDF_Model,double,relaxation,"Correlation relaxation parameter [um]","0",0xFF); DEFINE_PARAMETER(Growth_Roughness_Stack_BRDF_Model,double,exponent,"Correlation spatial frequency exponent parameter","2",0xFF); DEFINE_PTRPARAMETER(Growth_Roughness_Stack_BRDF_Model,PSD_Function_Ptr,intrinsic,"Intrinsic PSD of the film interfaces","ABC_PSD_Function",0xFF); diff --git a/code/allrough.h b/code/allrough.h index 1376be2..4f51c65 100644 --- a/code/allrough.h +++ b/code/allrough.h @@ -12,8 +12,8 @@ //** Version: 7.00 (January 2015) //** //****************************************************************************** -#ifndef SCATMECH_ALLROUGH_H -#define SCATMECH_ALLROUGH_H +#ifndef SCATMECH_GRADEDROUGH_H +#define SCATMECH_GRADEDROUGH_H #include "roughnes.h" @@ -26,7 +26,7 @@ namespace SCATMECH { public: DECLARE_MODEL(); - DECLARE_PARAMETER(dielectric_stack,stack); + DECLARE_PARAMETER(StackModel_Ptr,stack); protected: @@ -43,7 +43,7 @@ namespace SCATMECH { public: DECLARE_MODEL(); - DECLARE_PARAMETER(dielectric_stack,stack); + DECLARE_PARAMETER(StackModel_Ptr,stack); protected: @@ -60,7 +60,7 @@ namespace SCATMECH { public: DECLARE_MODEL(); - DECLARE_PARAMETER(dielectric_stack,stack); + DECLARE_PARAMETER(StackModel_Ptr,stack); DECLARE_PARAMETER(double,relaxation); DECLARE_PARAMETER(double,exponent); DECLARE_PARAMETER(PSD_Function_Ptr,intrinsic); diff --git a/code/askuser.cpp b/code/askuser.cpp index 70732d9..b2310b5 100644 --- a/code/askuser.cpp +++ b/code/askuser.cpp @@ -109,8 +109,8 @@ namespace SCATMECH { if (end==string::npos) end = pname.size(); string directory = pname.substr(0,end); if (directory[directory.size()-1] != delim) directory += delim; - string filename=directory+name; - ifstream file(filename.c_str()); + string filename=directory+name; + ifstream file(filename.c_str()); if (file) return filename; begin = end+1; } diff --git a/code/axipart.h b/code/axipart.h index 281360e..36dad7c 100644 --- a/code/axipart.h +++ b/code/axipart.h @@ -38,7 +38,7 @@ namespace SCATMECH { DECLARE_MODEL(); DECLARE_PARAMETER(Axisymmetric_Shape_Ptr,Shape); DECLARE_PARAMETER(dielectric_function,particle); - DECLARE_PARAMETER(dielectric_stack,stack); + DECLARE_PARAMETER(StackModel_Ptr,stack); DECLARE_PARAMETER(double,delta); DECLARE_PARAMETER(int,lmax); DECLARE_PARAMETER(int,mmax); diff --git a/code/axipart1.cpp b/code/axipart1.cpp index e0249e9..54c6fb1 100644 --- a/code/axipart1.cpp +++ b/code/axipart1.cpp @@ -144,8 +144,8 @@ namespace SCATMECH { } // Normal incidence reflection coefficients... - COMPLEX rsnormal = stack.rs12(0.,lambda,vacuum,substrate); - COMPLEX rpnormal = stack.rp12(0.,lambda,vacuum,substrate); + COMPLEX rsnormal = stack->rs12(0.,lambda,vacuum,substrate); + COMPLEX rpnormal = stack->rp12(0.,lambda,vacuum,substrate); vector reflect_s(npts); vector reflect_p(npts); @@ -169,8 +169,8 @@ namespace SCATMECH { reflect_s[j] = rsnormal; reflect_p[j] = rpnormal; } else { // Exact solution - reflect_s[j] = stack.rs12(alpha,lambda,vacuum,substrate); - reflect_p[j] = stack.rp12(alpha,lambda,vacuum,substrate); + reflect_s[j] = stack->rs12(alpha,lambda,vacuum,substrate); + reflect_p[j] = stack->rp12(alpha,lambda,vacuum,substrate); } // Calculate U, V, d+, and d- for each angle and (l,m) combination... @@ -624,7 +624,7 @@ namespace SCATMECH { { set_geometry(theta,theta,0.); COMPLEX e = E(W,Z); - COMPLEX r = stack.r12(theta,lambda,vacuum,substrate)[pol]; + COMPLEX r = stack->r12(theta,lambda,vacuum,substrate)[pol]; r *= exp(2.*cI*qq*cos(theta)); double R = norm(r); return 4.*pi/k*COMPLEX(0.,-1.)*(e/r)*R; @@ -641,7 +641,7 @@ namespace SCATMECH { COMPLEX phase = exp(cI*qq*(cos(theta)-index*cos(thetat))); double factor = cos(thetat)/cos(theta); e /= phase; - COMPLEX t = stack.t12(theta,lambda,vacuum,substrate)[pol]; + COMPLEX t = stack->t12(theta,lambda,vacuum,substrate)[pol]; double T = norm(t)*factor*index; return 4.*pi/k/sqrt(cube(index))/factor*COMPLEX(0.,-1.)*(e/t)*T; } @@ -656,7 +656,7 @@ namespace SCATMECH { COMPLEX cost = sqrt(1.-sqr(sint)); if (imag(cost)<0) cost = -cost; - COMPLEX r = stack.r21i(theta,lambda,substrate,vacuum)[pol]; + COMPLEX r = stack->r21i(theta,lambda,substrate,vacuum)[pol]; double R = norm(r); r *= exp(-2.*cI*qq*index*cos(theta)); @@ -680,7 +680,7 @@ namespace SCATMECH { COMPLEX phase = exp(cI*qq*(cost-index*cos(theta))); e /= phase; - COMPLEX t = stack.t21i(theta,lambda,substrate,vacuum)[pol]; + COMPLEX t = stack->t21i(theta,lambda,substrate,vacuum)[pol]; double T = norm(t)/index/factor; return 4.*pi/k*sqrt(index)*factor*COMPLEX(0.,-1.)*(e/t)*T; } @@ -700,7 +700,7 @@ namespace SCATMECH { "Particle Shape", "Ellipsoid_Axisymmetric_Shape",0xFF); DEFINE_PARAMETER(Axisymmetric_Particle_BRDF_Model,dielectric_function,particle,"Particle optical properties","(1.59,0)",0xFF); - DEFINE_PARAMETER(Axisymmetric_Particle_BRDF_Model,dielectric_stack,stack,"Substrate films","",0xFF); + DEFINE_PTRPARAMETER(Axisymmetric_Particle_BRDF_Model,StackModel_Ptr,stack,"Substrate films","No_StackModel",0xFF); DEFINE_PARAMETER(Axisymmetric_Particle_BRDF_Model,double,delta,"Separation of particle from substrate [um] (in contact: 0)","0",0xFF); DEFINE_PARAMETER(Axisymmetric_Particle_BRDF_Model,int,lmax,"Maximum polar order (lmax)","0",0xFF); DEFINE_PARAMETER(Axisymmetric_Particle_BRDF_Model,int,mmax,"Maximum azimuthal order (mmax)","0",0xFF); diff --git a/code/axipart2.cpp b/code/axipart2.cpp index 96c036b..69818af 100644 --- a/code/axipart2.cpp +++ b/code/axipart2.cpp @@ -62,7 +62,7 @@ namespace SCATMECH { COMPLEX result; COMPLEX cost = cos(thetai); COMPLEX phase = exp(2.*cI*qq*cost); - COMPLEX rp = stack.rp12(thetai,lambda,vacuum,substrate); + COMPLEX rp = stack->rp12(thetai,lambda,vacuum,substrate); COMPLEX costheta2=cos(thetai/2.); COMPLEX sintheta2=sin(thetai/2.); @@ -96,7 +96,7 @@ namespace SCATMECH { COMPLEX sintheta2=sqrt(1.-cost)/sqrt(2.); if (real(sintheta2)<0) sintheta2 = -sintheta2; COMPLEX costheta2=sint/2./sintheta2; - COMPLEX _tp = stack.tp12(_thetai,lambda,vacuum,substrate); + COMPLEX _tp = stack->tp12(_thetai,lambda,vacuum,substrate); COMPLEX phase = exp(cI*qq*(cost-substrate.n(lambda)*cos(thetai))); COMPLEX result; @@ -151,7 +151,7 @@ namespace SCATMECH { COMPLEX result; COMPLEX cost = cos(thetai); COMPLEX phase = exp(2.*cI*qq*cost); - COMPLEX rs = stack.rs12(thetai,lambda,vacuum,substrate); + COMPLEX rs = stack->rs12(thetai,lambda,vacuum,substrate); COMPLEX costheta2=cos(thetai/2.); COMPLEX sintheta2=sin(thetai/2.); @@ -185,7 +185,7 @@ namespace SCATMECH { COMPLEX sintheta2=sqrt(1.-cost)/sqrt(2.); if (real(sintheta2)<0) sintheta2 = -sintheta2; COMPLEX costheta2=sint/2./sintheta2; - COMPLEX _ts = stack.ts12(_thetai,lambda,vacuum,substrate); + COMPLEX _ts = stack->ts12(_thetai,lambda,vacuum,substrate); COMPLEX phase = exp(cI*qq*(cost-substrate.n(lambda)*cos(thetai))); COMPLEX result; @@ -342,8 +342,8 @@ namespace SCATMECH { double sint = real(sqrt(1.-sqr(cost))); COMPLEX _cost=cos(pi-thetas); - COMPLEX rp = stack.rp12(pi-thetas,lambda,vacuum,substrate); - COMPLEX rs = stack.rs12(pi-thetas,lambda,vacuum,substrate); + COMPLEX rp = stack->rp12(pi-thetas,lambda,vacuum,substrate); + COMPLEX rs = stack->rs12(pi-thetas,lambda,vacuum,substrate); COMPLEX phase = exp(2.*cI*qq*_cost); COMPLEX rpphase = rp*phase; COMPLEX rsphase = rs*phase; @@ -387,8 +387,8 @@ namespace SCATMECH { if (imag(_thetas)>0) _thetas = conj(_thetas); COMPLEX cost= sqrt(1.-sqr(sint)); if (imag(cost)<0) cost = -cost; - COMPLEX tp = stack.tp12(_thetas,lambda,vacuum,substrate); - COMPLEX ts = stack.ts12(_thetas,lambda,vacuum,substrate); + COMPLEX tp = stack->tp12(_thetas,lambda,vacuum,substrate); + COMPLEX ts = stack->ts12(_thetas,lambda,vacuum,substrate); COMPLEX phase = exp(cI*qq*(cost-substrate.n(lambda)*cos(pi-thetas))); // The following factor accounts for the transmittance across the interface and // the Jacobian as the solid angle across the interface changes... diff --git a/code/bobvlieg.h b/code/bobvlieg.h index 27f5ed8..09388ec 100644 --- a/code/bobvlieg.h +++ b/code/bobvlieg.h @@ -44,10 +44,10 @@ namespace SCATMECH { // Radius of particle... DECLARE_PARAMETER(double,radius); // Coatings on sphere... - DECLARE_PARAMETER(dielectric_stack,spherecoat); + DECLARE_PARAMETER(StackModel_Ptr,spherecoat); // Substrate coatings... - DECLARE_PARAMETER(dielectric_stack,stack); + DECLARE_PARAMETER(StackModel_Ptr,stack); // Distance between particle and substrate... DECLARE_PARAMETER(double,delta); // Scattering order (-1 means exact)... @@ -143,7 +143,7 @@ namespace SCATMECH { std::vector& b,std::vector& x); }; - class Gauss_Laguerre_Integration { + class Gauss_Laguerre_Integration { public: static double* weights[]; static double* zeros[]; diff --git a/code/bobvlieg1.cpp b/code/bobvlieg1.cpp index e51077f..73fe39c 100644 --- a/code/bobvlieg1.cpp +++ b/code/bobvlieg1.cpp @@ -71,13 +71,13 @@ namespace SCATMECH { // of C.F. Bohren and D.R.Huffman, "Absorption and Scattering of Light by Small Particles," (Wiley, New York, 1983). int n = l; - int rr = spherecoat.get_n()+1; // rr is the r in S&Q + int rr = spherecoat->get_n()+1; // rr is the r in S&Q if (f==efield) { COMPLEX T = 0.; double Rs = r0; COMPLEX ns = N2; for (int s=1; s<=rr-1; ++s) { - COMPLEX nsp1 = spherecoat.get_e()[s-1].index(lambda); + COMPLEX nsp1 = spherecoat->get_e()[s-1].index(lambda); COMPLEX ms = nsp1/ns; COMPLEX ys = 2*pi*ns/lambda*Rs; // Eq. (8a) of S&Q... @@ -85,9 +85,9 @@ namespace SCATMECH { (ms*chi(n,ms*ys)*(psi_(n,ys)+T*chi_(n,ys))-chi_(n,ms*ys)*(psi(n,ys)+T*chi(n,ys))); ns = nsp1; - Rs += spherecoat.get_t()[s-1]; + Rs += spherecoat->get_t()[s-1]; } - COMPLEX nr = spherecoat.get_e()[rr-2].index(lambda); + COMPLEX nr = spherecoat->get_e()[rr-2].index(lambda); COMPLEX mr = 1./nr; COMPLEX yr = 2*pi*nr/lambda*Rs; // Eq. (9a) of S&Q... @@ -99,16 +99,16 @@ namespace SCATMECH { double Rs = r0; COMPLEX ns = N2; for (int s=1; s<=rr-1; ++s) { - COMPLEX nsp1 = spherecoat.get_e()[s-1].index(lambda); + COMPLEX nsp1 = spherecoat->get_e()[s-1].index(lambda); COMPLEX ms = nsp1/ns; COMPLEX ys = 2*pi*ns/lambda*Rs; // Eq. (8b) of S&Q... S = -(psi(n,ms*ys)*(psi_(n,ys)+S*chi_(n,ys))-ms*psi_(n,ms*ys)*(psi(n,ys)+S*chi(n,ys)))/ (chi(n,ms*ys)*(psi_(n,ys)+S*chi_(n,ys))-ms*chi_(n,ms*ys)*(psi(n,ys)+S*chi(n,ys))); ns = nsp1; - Rs += spherecoat.get_t()[s-1]; + Rs += spherecoat->get_t()[s-1]; } - COMPLEX nr = spherecoat.get_e()[rr-2].index(lambda); + COMPLEX nr = spherecoat->get_e()[rr-2].index(lambda); COMPLEX mr = 1./nr; COMPLEX yr = 2*pi*nr/lambda*Rs; @@ -161,8 +161,8 @@ namespace SCATMECH { } // Normal incidence reflection coefficients... - COMPLEX rsnormal = stack.rs12(0.,lambda,vacuum,substrate); - COMPLEX rpnormal = stack.rp12(0.,lambda,vacuum,substrate); + COMPLEX rsnormal = stack->rs12(0.,lambda,vacuum,substrate); + COMPLEX rpnormal = stack->rp12(0.,lambda,vacuum,substrate); vector reflect_s(npts); vector reflect_p(npts); @@ -186,8 +186,8 @@ namespace SCATMECH { reflect_s[j] = rsnormal; reflect_p[j] = rpnormal; } else { // Exact solution - reflect_s[j] = stack.rs12(alpha,lambda,vacuum,substrate); - reflect_p[j] = stack.rp12(alpha,lambda,vacuum,substrate); + reflect_s[j] = stack->rs12(alpha,lambda,vacuum,substrate); + reflect_p[j] = stack->rp12(alpha,lambda,vacuum,substrate); } // Calculate U, V, d+, and d- for each angle and (l,m) combination... @@ -479,9 +479,9 @@ namespace SCATMECH { r0 = radius; d = 0.; - for (int i=0; iget_n(); ++i) { + r0 -= spherecoat->get_t()[i]; + d += spherecoat->get_t()[i]; } if (r0<=0) throw SCATMECH_exception("Radius of particle, r, less than or equal to total thickness of coatings."); @@ -853,7 +853,7 @@ namespace SCATMECH { { set_geometry(theta,theta,0.); COMPLEX e = E(W,Z); - COMPLEX r = stack.r12(theta,lambda,vacuum,substrate)[pol]; + COMPLEX r = stack->r12(theta,lambda,vacuum,substrate)[pol]; r *= exp(2.*cI*qq*cos(theta)); double R = norm(r); return 4.*pi/k*COMPLEX(0.,-1.)*(e/r)*R; @@ -870,7 +870,7 @@ namespace SCATMECH { COMPLEX phase = exp(cI*qq*(cos(theta)-index*cos(thetat))); double factor = cos(thetat)/cos(theta); e /= phase; - COMPLEX t = stack.t12(theta,lambda,vacuum,substrate)[pol]; + COMPLEX t = stack->t12(theta,lambda,vacuum,substrate)[pol]; double T = norm(t)*factor*index; return 4.*pi/k/sqrt(cube(index))/factor*COMPLEX(0.,-1.)*(e/t)*T; } @@ -885,7 +885,7 @@ namespace SCATMECH { COMPLEX cost = sqrt(1.-sqr(sint)); if (imag(cost)<0) cost = -cost; - COMPLEX r = stack.r21i(theta,lambda,substrate,vacuum)[pol]; + COMPLEX r = stack->r21i(theta,lambda,substrate,vacuum)[pol]; double R = norm(r); r *= exp(-2.*cI*qq*index*cos(theta)); @@ -909,7 +909,7 @@ namespace SCATMECH { COMPLEX phase = exp(cI*qq*(cost-index*cos(theta))); e /= phase; - COMPLEX t = stack.t21i(theta,lambda,substrate,vacuum)[pol]; + COMPLEX t = stack->t21i(theta,lambda,substrate,vacuum)[pol]; double T = norm(t)/index/factor; return 4.*pi/k*sqrt(index)*factor*COMPLEX(0.,-1.)*(e/t)*T; } @@ -925,8 +925,8 @@ namespace SCATMECH { DEFINE_PARAMETER(Bobbert_Vlieger_BRDF_Model,dielectric_function,sphere,"Sphere optical properties","(1.59,0)",0xFF); DEFINE_PARAMETER(Bobbert_Vlieger_BRDF_Model,double,radius,"Particle radius [um]","0.05",0xFF); - DEFINE_PARAMETER(Bobbert_Vlieger_BRDF_Model,dielectric_stack,spherecoat,"Coatings on the sphere","",0xFF); - DEFINE_PARAMETER(Bobbert_Vlieger_BRDF_Model,dielectric_stack,stack,"Substrate films","",0xFF); + DEFINE_PTRPARAMETER(Bobbert_Vlieger_BRDF_Model,StackModel_Ptr,spherecoat,"Coatings on the sphere","No_StackModel",0xFF); + DEFINE_PTRPARAMETER(Bobbert_Vlieger_BRDF_Model,StackModel_Ptr,stack,"Substrate films","No_StackModel",0xFF); DEFINE_PARAMETER(Bobbert_Vlieger_BRDF_Model,double,delta,"Separation of particle from substrate [um] (in contact: 0)","0",0xFF); DEFINE_PARAMETER(Bobbert_Vlieger_BRDF_Model,int,lmax,"Maximum spherical harmonic order (use Bohren & Huffman estimate: 0)","0",0xFF); DEFINE_PARAMETER(Bobbert_Vlieger_BRDF_Model,int,order,"Perturbation order (exact: -1)","-1",0xFF); diff --git a/code/bobvlieg2.cpp b/code/bobvlieg2.cpp index c943c7d..7b55df0 100644 --- a/code/bobvlieg2.cpp +++ b/code/bobvlieg2.cpp @@ -449,17 +449,25 @@ namespace SCATMECH { mutex.lock(); // First look in the tables to see if it has already been determined... lrhopairmap::iterator p = previous.find(lrhopair(l,rho)); - if (p!=previous.end()) return p->second; - + if (p!=previous.end()) { + mutex.unlock(); + return p->second; + } COMPLEX zm = sin(rho)/rho; previous[lrhopair(0,rho)]=zm; - if (l==0) return zm; + if (l==0) { + mutex.unlock(); + return zm; + } COMPLEX zmm = zm; zm = (zm-cos(rho))/rho; previous[lrhopair(1,rho)]=zm; - if (l==1) return zm; + if (l==1){ + mutex.unlock(); + return zm; + } for (int i=2; i<=l; ++i) { double nu = i+0.5; @@ -490,7 +498,10 @@ namespace SCATMECH { static mutex_t mutex; mutex.lock(); lrhopairmap::iterator p = previous.find(lrhopair(l,rho)); - if (p!=previous.end()) return p->second; + if (p!=previous.end()) { + mutex.unlock(); + return p->second; + } if (rho==0.) throw SCATMECH_exception("Encountered y(l,0.)"); @@ -561,8 +572,10 @@ namespace SCATMECH { mutex_t mutex; mutex.lock(); lrhopairmap::iterator p = previous.find(lrhopair(l,rho)); - if (p!=previous.end()) return p->second; - + if (p!=previous.end()) { + mutex.unlock(); + return p->second; + } COMPLEX result = psi(l,rho)*ContinuedFraction_A(l,rho); previous[lrhopair(l,rho)]=result; @@ -686,7 +699,7 @@ namespace SCATMECH { COMPLEX result; COMPLEX cost = cos(thetai); COMPLEX phase = exp(2.*cI*qq*cost); - COMPLEX _rp = stack.rp12(thetai,lambda,vacuum,substrate); + COMPLEX _rp = stack->rp12(thetai,lambda,vacuum,substrate); COMPLEX temp = phase*_rp; COMPLEX costheta2=cos(thetai/2.); COMPLEX sintheta2=sin(thetai/2.); @@ -721,7 +734,7 @@ namespace SCATMECH { COMPLEX sintheta2=sqrt(1.-cost)/sqrt(2.); if (real(sintheta2)<0) sintheta2 = -sintheta2; COMPLEX costheta2=sint/2./sintheta2; - COMPLEX _tp = stack.tp12(_thetai,lambda,vacuum,substrate); + COMPLEX _tp = stack->tp12(_thetai,lambda,vacuum,substrate); COMPLEX phase = exp(cI*qq*(cost-substrate.n(lambda)*cos(thetai))); COMPLEX result; @@ -776,7 +789,7 @@ namespace SCATMECH { COMPLEX result; COMPLEX cost = cos(thetai); COMPLEX phase = exp(2.*cI*qq*cost); - COMPLEX _rs = stack.rs12(thetai,lambda,vacuum,substrate); + COMPLEX _rs = stack->rs12(thetai,lambda,vacuum,substrate); COMPLEX costheta2=cos(thetai/2.); COMPLEX sintheta2=sin(thetai/2.); @@ -810,7 +823,7 @@ namespace SCATMECH { COMPLEX sintheta2=sqrt(1.-cost)/sqrt(2.); if (real(sintheta2)<0) sintheta2 = -sintheta2; COMPLEX costheta2=sint/2./sintheta2; - COMPLEX _ts = stack.ts12(_thetai,lambda,vacuum,substrate); + COMPLEX _ts = stack->ts12(_thetai,lambda,vacuum,substrate); COMPLEX phase = exp(cI*qq*(cost-substrate.n(lambda)*cos(thetai))); COMPLEX result; @@ -990,8 +1003,8 @@ namespace SCATMECH { double sint = real(sqrt(1.-sqr(cost))); COMPLEX _cost=cos(pi-thetas); - COMPLEX rp = stack.rp12(pi-thetas,lambda,vacuum,substrate); - COMPLEX rs = stack.rs12(pi-thetas,lambda,vacuum,substrate); + COMPLEX rp = stack->rp12(pi-thetas,lambda,vacuum,substrate); + COMPLEX rs = stack->rs12(pi-thetas,lambda,vacuum,substrate); COMPLEX phase = exp(2.*cI*qq*_cost); COMPLEX rpphase = rp*phase; COMPLEX rsphase = rs*phase; @@ -1033,8 +1046,8 @@ namespace SCATMECH { if (imag(_thetas)>0) _thetas = conj(_thetas); COMPLEX cost= sqrt(1.-sqr(sint)); if (imag(cost)<0) cost = -cost; - COMPLEX tp = stack.tp12(_thetas,lambda,vacuum,substrate); - COMPLEX ts = stack.ts12(_thetas,lambda,vacuum,substrate); + COMPLEX tp = stack->tp12(_thetas,lambda,vacuum,substrate); + COMPLEX ts = stack->ts12(_thetas,lambda,vacuum,substrate); COMPLEX phase = exp(cI*qq*(cost-substrate.n(lambda)*cos(pi-thetas))); // The following factor accounts for the transmittance across the interface and // the Jacobian as the solid angle across the interface changes... diff --git a/code/brdf.cpp b/code/brdf.cpp index 831959f..06bee3c 100644 --- a/code/brdf.cpp +++ b/code/brdf.cpp @@ -163,7 +163,17 @@ namespace SCATMECH { if (from==xyxy) { // Convert to BRDF_Model::psps - J = JonesRotator(phis)*J; + Vector ki = (type==0||type==1) ? Vector(sin(thetai),0.,-cos(thetai)) : Vector(sin(thetai),0.,cos(thetai)); + Vector ko = (type==0||type==3) ? Vector(sin(thetas)*cos(phis),sin(thetas)*sin(phis),cos(thetas)) : Vector(sin(thetas)*cos(phis),sin(thetas)*sin(phis),-cos(thetas)); + Vector xi,yi,xo,yo,si,pi,so,po; + GetBasisVectorsSP(ko,so,po); + GetBasisVectorsSP(ki,si,pi); + GetBasisVectorsXY(ko,xo,yo); + GetBasisVectorsXY(ki,xi,yi); + JonesMatrix ri = GetJonesRotator(yi,xi,si,pi); + JonesMatrix ro = GetJonesRotator(so,po,yo,xo); + J = ro*J*ri; + } else if (from==plane) { // Convert to BRDF_Model::psps double psii,psis; @@ -173,15 +183,38 @@ namespace SCATMECH { if (to == xyxy) { // Convert to BRDF_Model::xyxy - J = JonesRotator(-phis)*J; + Vector ki = (type==0||type==1) ? Vector(sin(thetai),0.,-cos(thetai)) : Vector(sin(thetai),0.,cos(thetai)); + Vector ko = (type==0||type==3) ? Vector(sin(thetas)*cos(phis),sin(thetas)*sin(phis),cos(thetas)) : Vector(sin(thetas)*cos(phis),sin(thetas)*sin(phis),-cos(thetas)); + Vector xi,yi,xo,yo,si,pi,so,po; + GetBasisVectorsSP(ko,so,po); + GetBasisVectorsSP(ki,si,pi); + GetBasisVectorsXY(ko,xo,yo); + GetBasisVectorsXY(ki,xi,yi); + JonesMatrix ri = GetJonesRotator(si,pi,yi,xi); + JonesMatrix ro = GetJonesRotator(yo,xo,so,po); + J = ro*J*ri; + + //if (type==0||type==2) { // Added 21 March 2016 + // J = JonesRotator(-phis)*J; + //} else { + // J = JonesRotator(phis)*J; + //} } else if (to==plane) { // Convert to BRDF_Model::plane - double psii,psis; - Coordinate_System_plane_angles(psii,psis,thetai,thetas,phis); - J = (JonesRotator(-psis)*J)*JonesRotator(-psii); + Vector ki = type==0||type==1 ? Vector(sin(thetai),0.,-cos(thetai)) : Vector(sin(thetai),0.,cos(thetai)); + Vector ko = type==0||type==3 ? Vector(sin(thetas)*cos(phis),sin(thetas)*sin(phis),cos(thetas)) : Vector(sin(thetas)*cos(phis),sin(thetas)*sin(phis),-cos(thetas)); + Vector si,pi,so,po,perp,parin,parout; + GetBasisVectorsSP(ko,so,po); + GetBasisVectorsSP(ki,si,pi); + GetBasisVectorsParPerp(ki,ko,perp,parin,parout); + + J = (GetJonesRotator(perp,parout,so,po)*J)*GetJonesRotator(si,pi,perp,parin); + //double psii,psis; + //Coordinate_System_plane_angles(psii,psis,thetai,thetas,phis); + //J = (JonesRotator(-psis)*J)*JonesRotator(-psii); } } - + void BRDF_Model:: convert(MuellerMatrix& M,Coordinate_System to) const @@ -191,7 +224,16 @@ namespace SCATMECH { if (from==xyxy) { // Convert to BRDF_Model_psps - M = ((MuellerMatrix)JonesRotator(phis))*M; + Vector ki = (type==0||type==1) ? Vector(sin(thetai),0.,-cos(thetai)) : Vector(sin(thetai),0.,cos(thetai)); + Vector ko = (type==0||type==3) ? Vector(sin(thetas)*cos(phis),sin(thetas)*sin(phis),cos(thetas)) : Vector(sin(thetas)*cos(phis),sin(thetas)*sin(phis),-cos(thetas)); + Vector xi,yi,xo,yo,si,pi,so,po; + GetBasisVectorsSP(ko,so,po); + GetBasisVectorsSP(ki,si,pi); + GetBasisVectorsXY(ko,xo,yo); + GetBasisVectorsXY(ki,xi,yi); + MuellerMatrix ri = GetJonesRotator(yi,xi,si,pi); + MuellerMatrix ro = GetJonesRotator(so,po,yo,xo); + M = ro*M*ri; } else if (from==plane) { // Convert to BRDF_Model_psps double psii,psis; @@ -202,13 +244,36 @@ namespace SCATMECH { if (to == xyxy) { // Convert to BRDF_Model_xyxy - M = ((MuellerMatrix)JonesRotator(-phis))*M; + Vector ki = (type==0||type==1) ? Vector(sin(thetai),0.,-cos(thetai)) : Vector(sin(thetai),0.,cos(thetai)); + Vector ko = (type==0||type==3) ? Vector(sin(thetas)*cos(phis),sin(thetas)*sin(phis),cos(thetas)) : Vector(sin(thetas)*cos(phis),sin(thetas)*sin(phis),-cos(thetas)); + Vector xi,yi,xo,yo,si,pi,so,po; + GetBasisVectorsSP(ko,so,po); + GetBasisVectorsSP(ki,si,pi); + GetBasisVectorsXY(ko,xo,yo); + GetBasisVectorsXY(ki,xi,yi); + MuellerMatrix ri = GetJonesRotator(si,pi,yi,xi); + MuellerMatrix ro = GetJonesRotator(yo,xo,so,po); + M = ro*M*ri; + //if (type==0 || type == 2) { + // M = ((MuellerMatrix)JonesRotator(-phis))*M; + //} else { + // M = ((MuellerMatrix)JonesRotator(phis))*M; + //} } else if (to==plane) { // Convert to BRDF_Model_plane - double psii,psis; - Coordinate_System_plane_angles(psii,psis,thetai,thetas,phis); - M= (((MuellerMatrix)JonesRotator(-psis))*M)* - ((MuellerMatrix)JonesRotator(-psii)); + //double psii,psis; + //Coordinate_System_plane_angles(psii,psis,thetai,thetas,phis); + //M= (((MuellerMatrix)JonesRotator(-psis))*M)* + // ((MuellerMatrix)JonesRotator(-psii)); + Vector ki = type==0||type==1 ? Vector(sin(thetai),0.,-cos(thetai)) : Vector(sin(thetai),0.,cos(thetai)); + Vector ko = type==0||type==3 ? Vector(sin(thetas)*cos(phis),sin(thetas)*sin(phis),cos(thetas)) : Vector(sin(thetas)*cos(phis),sin(thetas)*sin(phis),-cos(thetas)); + Vector si,pi,so,po,perp,parin,parout; + GetBasisVectorsSP(ko,so,po); + GetBasisVectorsSP(ki,si,pi); + GetBasisVectorsParPerp(ki,ko,perp,parin,parout); + + M = (MuellerMatrix(GetJonesRotator(perp,parout,so,po))*M)*MuellerMatrix(GetJonesRotator(si,pi,perp,parin)); + } } @@ -297,7 +362,7 @@ namespace SCATMECH { DEFINE_VIRTUAL_MODEL(BRDF_Model,Model,"All BRDF models."); DEFINE_PARAMETER(BRDF_Model,double,lambda,"Wavelength [um]","0.532",0xFF); DEFINE_PARAMETER(BRDF_Model,dielectric_function,substrate,"Substrate","(4.05,0.05)",0xFF); - DEFINE_PARAMETER(BRDF_Model,int,type,"(0) for Forward/Reflection, (1) for Forward/Transmission, (2) for Backward/Reflection, or (3) for Backward/Transmission","0",0xFF); + DEFINE_PARAMETER(BRDF_Model,int,type,"(0) for Forward/Reflection, (1) for Forward/Transmission, (2) for Backward/Reflection, or (3) for Backward/Transmission","0",0x80); int get_model_depth=0; diff --git a/code/crossgrating.cpp b/code/crossgrating.cpp index 0de503f..078f70e 100644 --- a/code/crossgrating.cpp +++ b/code/crossgrating.cpp @@ -18,8 +18,9 @@ #include "matrixmath.h" #include "scatmech.h" #include "fft.h" +#include "bobvlieg.h" #include - + using namespace std; namespace SCATMECH { @@ -47,6 +48,7 @@ namespace SCATMECH { void Gridded_CrossGrating::FourierFactorize() { + using BobVlieg_Supp::mpow; CrossGrating::zeta = Gridded_CrossGrating::zeta; CrossGrating::d1 = Gridded_CrossGrating::d1; CrossGrating::d2 = Gridded_CrossGrating::d2; diff --git a/code/crossrcw.cpp b/code/crossrcw.cpp index 052f4f8..a83832a 100644 --- a/code/crossrcw.cpp +++ b/code/crossrcw.cpp @@ -805,6 +805,18 @@ namespace SCATMECH { return (type&0x01) ? T(i,j) : R(i,j); } + StokesVector CrossRCW_Model::GetAbsorption() + { + SETUP(); + StokesVector result(0,0,0,0); + for (int i=1;i<=2*order1+1;++i) { + for (int j=1;j<=2*order2+1;++j) { + result += (T(i,j) + R(i,j)).transpose()*StokesVector(1,0,0,0); + } + } + return StokesVector(1,0,0,0)-result; + } + CVector CrossRCW_Model::GetEField(const JonesVector& inpol, const Vector& pos, bool incident) { SETUP(); diff --git a/code/crossrcw.h b/code/crossrcw.h index c4d3943..051ccfc 100644 --- a/code/crossrcw.h +++ b/code/crossrcw.h @@ -35,8 +35,9 @@ namespace SCATMECH { public: CrossRCW_Model() : old_type(-1) {} - JonesMatrix GetAmplitude(int i,int j); - MuellerMatrix GetIntensity(int i,int j); + JonesMatrix GetAmplitude(int i,int j); + MuellerMatrix GetIntensity(int i,int j); // Returns diffraction efficiency + StokesVector GetAbsorption(); // Returns Stokes absorption coefficient Vector GetDirection(int i,int j); CVector GetPropagationVector(int i,int j); diff --git a/code/crough.cpp b/code/crough.cpp index 91fe401..a3b7a37 100644 --- a/code/crough.cpp +++ b/code/crough.cpp @@ -31,9 +31,10 @@ namespace SCATMECH { model.set_type(type); model.set_substrate(substrate); model.set_psd(psd); - dielectric_stack stack; - stack.grow(film,thickness); - model.set_stack(stack); + SingleFilm_StackModel stack; + stack.set_material(film); + stack.set_thickness(thickness); + model.set_stack(stack.clone()); } // // Jones matrix for scattering... diff --git a/code/dielfunc.cpp b/code/dielfunc.cpp index 428a95d..bb4c88a 100644 --- a/code/dielfunc.cpp +++ b/code/dielfunc.cpp @@ -172,7 +172,8 @@ namespace SCATMECH { set_n(const Table& n) { n1 = n; - last.x=-1; + last.x=-1.; + index(1.); } void @@ -180,7 +181,8 @@ namespace SCATMECH { set_k(const Table& k) { n2 = k; - last.x=-1; + last.x=-1.; + index(1.); } void diff --git a/code/dielfunc.h b/code/dielfunc.h index bf0ae0d..062ea20 100644 --- a/code/dielfunc.h +++ b/code/dielfunc.h @@ -115,6 +115,9 @@ namespace SCATMECH { void set_n(const Table& n); void set_k(const Table& k); + const Table& get_n() {return n1;} + const Table& get_k() {return n2;} + private: Table n1; Table n2; diff --git a/code/diffuse.cpp b/code/diffuse.cpp index 43d1dbf..e559c05 100644 --- a/code/diffuse.cpp +++ b/code/diffuse.cpp @@ -34,7 +34,6 @@ namespace SCATMECH { Diffuse_Subsurface_BRDF_Model:: Diffuse_Subsurface_BRDF_Model() { - stack.wash(); } MuellerMatrix Diffuse_Subsurface_BRDF_Model:: @@ -55,8 +54,8 @@ namespace SCATMECH { MuellerMatrix m = Lambertian_BRDF_Model::mueller(); - return s*(((MuellerMatrix)(stack.t12(thetas,lambda,vacuum,substrate))*m)* - (MuellerMatrix)(stack.t12(thetai,lambda,vacuum,substrate))); + return s*(((MuellerMatrix)(stack->t12(thetas,lambda,vacuum,substrate))*m)* + (MuellerMatrix)(stack->t12(thetai,lambda,vacuum,substrate))); } void @@ -76,7 +75,7 @@ namespace SCATMECH { COMPLEX ti = ArcSin((COMPLEX)(sin(theta)*substrate.n(lambda))); // Calculate reflectance... - JonesMatrix r = stack.r21(ti,lambda,substrate,vacuum); + JonesMatrix r = stack->r21(ti,lambda,substrate,vacuum); double rr=0.5*(norm(r[0])+norm(r[1])); // Integrate the reflectance... @@ -95,7 +94,7 @@ namespace SCATMECH { DEFINE_MODEL(Diffuse_Subsurface_BRDF_Model,Lambertian_BRDF_Model, "Totally diffuse scattering under a smooth surface."); - DEFINE_PARAMETER(Diffuse_Subsurface_BRDF_Model,dielectric_stack,stack,"Film stack on substrate","",0xFF); + DEFINE_PTRPARAMETER(Diffuse_Subsurface_BRDF_Model,StackModel_Ptr,stack,"Film stack on substrate","No_StackModel",0xFF); } // namespace SCATMECH diff --git a/code/diffuse.h b/code/diffuse.h index 035aeb2..b2e90a7 100644 --- a/code/diffuse.h +++ b/code/diffuse.h @@ -32,7 +32,7 @@ namespace SCATMECH { DECLARE_MODEL(); - DECLARE_PARAMETER(dielectric_stack,stack); + DECLARE_PARAMETER(StackModel_Ptr,stack); protected: MuellerMatrix mueller(); diff --git a/code/facet.cpp b/code/facet.cpp index 8f3e9c7..442545b 100644 --- a/code/facet.cpp +++ b/code/facet.cpp @@ -138,11 +138,11 @@ namespace SCATMECH { double iota = acos(cos_iota_i); COMPLEX ts,tp; if (is_forward()) { - ts = stack.ts12(iota,lambda,vacuum,substrate); - tp = stack.tp12(iota,lambda,vacuum,substrate); + ts = stack->ts12(iota,lambda,vacuum,substrate); + tp = stack->tp12(iota,lambda,vacuum,substrate); } else { - ts = stack.ts21i(iota,lambda,substrate,vacuum); - tp = stack.tp21i(iota,lambda,substrate,vacuum); + ts = stack->ts21i(iota,lambda,substrate,vacuum); + tp = stack->tp21i(iota,lambda,substrate,vacuum); } JonesMatrix j; @@ -206,11 +206,11 @@ namespace SCATMECH { COMPLEX rsiota,rpiota; if (is_forward()) { - rsiota = stack.rs12(iota,lambda,vacuum,substrate); - rpiota = stack.rp12(iota,lambda,vacuum,substrate); + rsiota = stack->rs12(iota,lambda,vacuum,substrate); + rpiota = stack->rp12(iota,lambda,vacuum,substrate); } else { - rsiota = stack.rs21i(iota,lambda,substrate,vacuum); - rpiota = stack.rp21i(iota,lambda,substrate,vacuum); + rsiota = stack->rs21i(iota,lambda,substrate,vacuum); + rpiota = stack->rp21i(iota,lambda,substrate,vacuum); } double a1 = sqr(sin(2*iota)); @@ -434,7 +434,7 @@ namespace SCATMECH { DEFINE_PTRPARAMETER(Facet_BRDF_Model,Slope_Distribution_Function_Ptr,sdf,"Slope distribution function","Exponential_Slope_Distribution_Function",0xFF); - DEFINE_PARAMETER(Facet_BRDF_Model,dielectric_stack,stack,"Film stack on substrate","",0xFF); + DEFINE_PTRPARAMETER(Facet_BRDF_Model,StackModel_Ptr,stack,"Film stack on substrate","No_StackModel",0xFF); DEFINE_PARAMETER(Exponential_Slope_Distribution_Function,double,s,"RMS slope","0.1",0xFF); diff --git a/code/facet.h b/code/facet.h index 6920cb0..3574a29 100644 --- a/code/facet.h +++ b/code/facet.h @@ -53,7 +53,7 @@ namespace SCATMECH { DECLARE_MODEL(); DECLARE_PARAMETER(Slope_Distribution_Function_Ptr,sdf); - DECLARE_PARAMETER(dielectric_stack,stack); + DECLARE_PARAMETER(StackModel_Ptr,stack); protected: JonesMatrix jones(); diff --git a/code/filmtran.cpp b/code/filmtran.cpp index a7e7f29..2bad92e 100644 --- a/code/filmtran.cpp +++ b/code/filmtran.cpp @@ -1,662 +1,708 @@ -//****************************************************************************** -//** SCATMECH: Polarized Light Scattering C++ Class Library -//** -//** File: filmtran.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 -//** -//** Version: 7.00 (January 2015) -//** -//****************************************************************************** - -//***************************************************************************** -//* -//* Operations needed to calculate reflectance from and transmittance through -//* a film stack. -//* -//* See G. R. Fowles, _Introduction_to_Modern_Optics_, Second Edition, -//* (Holt, Rinehart and Winston, New York, 1975) -//* Section 4.4, Theory of Multilayer Films -//* Note that some of the work is left to the reader in Problem 4.10. -//* -//***************************************************************************** - -// -// In all of functions in this file, the angles were changed to be complex. -// (Version 3.02, TAG 7 AUG 2002) -// - -//***************************************************************************** -// Major changes were made in SCATMECH version 3 to the classes -// dielectric_stack and Thin_Film_Transfer_Matrix so that the wavelength -// could more easily be changed. (TAG 19 JAN 2001) -//***************************************************************************** -#include -#include -#include -#include -#include "scatmech.h" -#include "filmtran.h" -#include "fresnel.h" -#include "askuser.h" - -using namespace std; - - -namespace SCATMECH { - - - - // Largest number of layers envisioned... - // static const int MAXSTACK = 64; - - // The square root of -1... - static const COMPLEX I = COMPLEX(0,1); - - - // - // The following gives the transfer matrices for a - // layer of given thickness (normalized by 2*pi/lambda), - // given index, and given external angle of incidence. - // - Thin_Film_Transfer_Matrix - Transfer_Matrix(const optical_constant& n,double thickness,COMPLEX theta,double lambda) - { - Thin_Film_Transfer_Matrix tftm; - COMPLEX cos_thetai=cos_internal_angle(n,theta); - COMPLEX klcostheta = (COMPLEX)n*2.*pi*thickness/lambda*cos_thetai; - COMPLEX p = cos_thetai*(COMPLEX)n; - COMPLEX q = cos_thetai/(COMPLEX)n; - COMPLEX sinklcostheta = sin(klcostheta); - tftm.As = tftm.Ap = cos(klcostheta); - tftm.Bs = -I/p*sinklcostheta; - tftm.Cs = -I*p*sinklcostheta; - tftm.Bp = -I/q*sinklcostheta; - tftm.Cp = -I*q*sinklcostheta; - tftm.Ds = tftm.Dp = tftm.As; - return tftm; - } - - // - // The following gives the reflectance and transmittance coefficients for - // s- or p-polarized light with indices n0 and nt surrounding the layers. - // The light is assumed to be incident on the n0 side. That is nt is the - // index of the substrate. - // - COMPLEX - Thin_Film_Transfer_Matrix:: - rs(COMPLEX theta,const optical_constant& n0,const optical_constant& nt) - { - COMPLEX p0 = (COMPLEX)n0*cos_internal_angle(n0,theta); - COMPLEX pt = (COMPLEX)nt*cos_internal_angle(nt,theta); - return ((As+Bs*pt)*p0-(Cs+Ds*pt))/((As+Bs*pt)*p0+(Cs+Ds*pt)); - } - - COMPLEX - Thin_Film_Transfer_Matrix:: - rp(COMPLEX theta,const optical_constant& n0,const optical_constant& nt) - { - COMPLEX q0 = cos_internal_angle(n0,theta)/(COMPLEX)n0; - COMPLEX qt = cos_internal_angle(nt,theta)/(COMPLEX)nt; - return ((Ap+Bp*qt)*q0-(Cp+Dp*qt))/((Ap+Bp*qt)*q0+(Cp+Dp*qt)); - } - - COMPLEX - Thin_Film_Transfer_Matrix:: - ts(COMPLEX theta,const optical_constant& n0,const optical_constant& nt) - { - COMPLEX p0 = (COMPLEX)n0*cos_internal_angle(n0,theta); - COMPLEX pt = (COMPLEX)nt*cos_internal_angle(nt,theta); - return (2.*p0)/((As+Bs*pt)*p0+(Cs+Ds*pt)); - } - - COMPLEX - Thin_Film_Transfer_Matrix:: - tp(COMPLEX theta,const optical_constant& n0,const optical_constant& nt) - { - COMPLEX p0 = cos_internal_angle(n0,theta)/(COMPLEX)n0; - COMPLEX pt = cos_internal_angle(nt,theta)/(COMPLEX)nt; - return (2.*p0)/((Ap+Bp*pt)*p0+(Cp+Dp*pt))*(COMPLEX)n0/(COMPLEX)nt; - } - - // - // Multiple layers can be accounted for my multiplication of the transfer - // matrices for each layer... - // - Thin_Film_Transfer_Matrix - Thin_Film_Transfer_Matrix:: - operator*(Thin_Film_Transfer_Matrix a) - { - Thin_Film_Transfer_Matrix t; - t.As = As*a.As+Bs*a.Cs; - t.Ap = Ap*a.Ap+Bp*a.Cp; - t.Bs = As*a.Bs+Bs*a.Ds; - t.Bp = Ap*a.Bp+Bp*a.Dp; - t.Cs = Cs*a.As+Ds*a.Cs; - t.Cp = Cp*a.Ap+Dp*a.Cp; - t.Ds = Cs*a.Bs+Ds*a.Ds; - t.Dp = Cp*a.Bp+Dp*a.Dp; - return t; - } - - // - // The dielectric_stack class holds information about a stack of films... - // - - // - // Copy constructor... - // - /*dielectric_stack:: - dielectric_stack(const dielectric_stack& ds) - { - n=ds.n; - e.resize(0); - t.resize(0); - int i; - for (i=0;i=0; --i) - M = M*Transfer_Matrix(e[i].index(lambda),t[i],theta,lambda); - return M; - } - - // - // The net transfer matrix for light from the "substrate" side of the layer. - // - Thin_Film_Transfer_Matrix - dielectric_stack:: - matrix_backward(COMPLEX theta,double lambda) const - { - Thin_Film_Transfer_Matrix M; - int i; - for (i=0; i :"; - SCATMECH_output << prompt2; - - // Get response from user... - string response = SCATMECH_input.getstrline(); - - - if (response.size()==0) { - // If response is blank, return the default... - stack.read_stack_file(deflt); - } else if (response[0]==':') { - // If response begins with a colon, use the next word as a filename... - istringstream responsestr(response); - responsestr.get(); - string filename; - responsestr >> filename; - stack.read_stack_file(filename); - } else { - int i=1; - bool finished=false; - - string dielfunc; - double thickness=-1; - - while (!finished) { - - // Extract the dielectric function and thickness of the layer... - istringstream sresponse(response); - dielfunc.clear(); - thickness=-1; - sresponse >> dielfunc >> thickness; - - // If no dielectric function, assume finished... - if (dielfunc.size()==0) finished=true; - // Otherwise... - else { - // If thickness is less than zero, or no thickness given... - if (thickness<0.) throw SCATMECH_exception("Thickness must be greater than zero."); - - // If both the dielectric function and the thickness is good, then grow film... - stack.grow(dielectric_function(dielfunc),thickness); - - // Prompt user for next layer... - ostringstream message; - message << "Layer #" << ++i << ":"; - SCATMECH_output << message.str() << endl; - - // Get next response... - response = getstrline(sresponse); - string next; - istringstream temp(response); - temp >> next; - if (next.size()==0) { - response = SCATMECH_input.getstrline(); - } - } - } - } - ostringstream message; - stack.print(message); - SCATMECH_output << message.str(); - return stack; - } - catch(const exception& e) { - SCATMECH_output << e.what() << endl; - //instack_unwind(); - } - } - } - - void - dielectric_stack:: - read_stack_file(const string& film_file_name) - { - // If string pointer is NULL.. then assume no film... - if (film_file_name.size()==0) { - wash(); - return; - } - - // Open file stream... - string fname = find_file(film_file_name); - ifstream_with_comments film_file(fname.c_str()); - // and complain if it doesn't exist... - - if (!film_file) { - throw SCATMECH_exception("Cannot open file: " + film_file_name); - } - - do { - // Read values... - string material; - film_file >> material; - - // A blank line exits the routine... - if (material.size()==0) return; - - // Get the film thickness - double thickness=-1.; - film_file >> thickness; - if (thickness<0.) throw SCATMECH_exception("Thickness must be greater than zero."); - - // Deposit that film... - grow(dielectric_function(material),thickness); - - } while (!film_file.eof()); - } - - // - // The following prints out the status of the film... - // - void - dielectric_stack:: - print(ostream& os) const - { - os << endl << "# index thickness " << endl - << "-----------------------------------" << endl; - if (n==0) { - os << "No layers" << endl; - } else { - for (int i=0; i - void - ModelParameterSet(dielectric_stack& variable,const string& parameter,const string& value) - { - if (parameter.empty()) { - if (value.empty()) { - variable.wash(); - return; - } - istringstream iss(value); - string material; - iss >> material; - if (iss.fail()) throw SCATMECH_exception("Invalid parameter for dielectric_stack"); - double thickness; - iss >> thickness; - if (iss.fail()) { - variable.read_stack_file(material); - return; - } - variable.grow(material,thickness); - while (1) { - iss >> material; - if (iss.fail()) return; - iss >> thickness; - if (iss.fail()) throw SCATMECH_exception("Invalid thickness parameter for dielectric_stack"); - variable.grow(material,thickness); - } - } else if (parameter == "grow") { - istringstream iss(value); - string material; - double thickness; - iss >> material >> thickness; - if (iss.fail()) throw SCATMECH_exception("Invalid parameter for dielectric_stack"); - variable.grow(material,thickness); - } else if (parameter == "wash") { - variable.wash(); - } - } - - template <> - string - ModelParameterGet(dielectric_stack& variable,const string& parameter) - { - ostringstream oss; - oss << variable; - return oss.str(); - } - - template <> - void - ModelParameterAskUser(dielectric_stack& value,const string& prompt) - { - value = dielectric_stack::AskUser(prompt); - } - - -} // namespace SCATMECH - - - +//****************************************************************************** +//** SCATMECH: Polarized Light Scattering C++ Class Library +//** +//** File: filmtran.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 +//** +//** Version: 7.00 (January 2015) +//** +//****************************************************************************** + +//***************************************************************************** +//* +//* Operations needed to calculate reflectance from and transmittance through +//* a film stack. +//* +//* See G. R. Fowles, _Introduction_to_Modern_Optics_, Second Edition, +//* (Holt, Rinehart and Winston, New York, 1975) +//* Section 4.4, Theory of Multilayer Films +//* Note that some of the work is left to the reader in Problem 4.10. +//* +//***************************************************************************** + +// +// In all of functions in this file, the angles were changed to be complex. +// (Version 3.02, TAG 7 AUG 2002) +// + +//***************************************************************************** +// Major changes were made in SCATMECH version 3 to the classes +// dielectric_stack and Thin_Film_Transfer_Matrix so that the wavelength +// could more easily be changed. (TAG 19 JAN 2001) +//***************************************************************************** +#include +#include +#include +#include +#include "scatmech.h" +#include "filmtran.h" +#include "fresnel.h" +#include "askuser.h" + +using namespace std; + + +namespace SCATMECH { + + + + // Largest number of layers envisioned... + // static const int MAXSTACK = 64; + + // The square root of -1... + static const COMPLEX I = COMPLEX(0,1); + + + // + // The following gives the transfer matrices for a + // layer of given thickness (normalized by 2*pi/lambda), + // given index, and given external angle of incidence. + // + Thin_Film_Transfer_Matrix + Transfer_Matrix(const optical_constant& n,double thickness,COMPLEX theta,double lambda) + { + Thin_Film_Transfer_Matrix tftm; + COMPLEX cos_thetai=cos_internal_angle(n,theta); + COMPLEX klcostheta = (COMPLEX)n*2.*pi*thickness/lambda*cos_thetai; + COMPLEX p = cos_thetai*(COMPLEX)n; + COMPLEX q = cos_thetai/(COMPLEX)n; + COMPLEX sinklcostheta = sin(klcostheta); + tftm.As = tftm.Ap = cos(klcostheta); + tftm.Bs = -I/p*sinklcostheta; + tftm.Cs = -I*p*sinklcostheta; + tftm.Bp = -I/q*sinklcostheta; + tftm.Cp = -I*q*sinklcostheta; + tftm.Ds = tftm.Dp = tftm.As; + return tftm; + } + + // + // The following gives the reflectance and transmittance coefficients for + // s- or p-polarized light with indices n0 and nt surrounding the layers. + // The light is assumed to be incident on the n0 side. That is nt is the + // index of the substrate. + // + COMPLEX + Thin_Film_Transfer_Matrix:: + rs(COMPLEX theta,const optical_constant& n0,const optical_constant& nt) + { + COMPLEX p0 = (COMPLEX)n0*cos_internal_angle(n0,theta); + COMPLEX pt = (COMPLEX)nt*cos_internal_angle(nt,theta); + return ((As+Bs*pt)*p0-(Cs+Ds*pt))/((As+Bs*pt)*p0+(Cs+Ds*pt)); + } + + COMPLEX + Thin_Film_Transfer_Matrix:: + rp(COMPLEX theta,const optical_constant& n0,const optical_constant& nt) + { + COMPLEX q0 = cos_internal_angle(n0,theta)/(COMPLEX)n0; + COMPLEX qt = cos_internal_angle(nt,theta)/(COMPLEX)nt; + return ((Ap+Bp*qt)*q0-(Cp+Dp*qt))/((Ap+Bp*qt)*q0+(Cp+Dp*qt)); + } + + COMPLEX + Thin_Film_Transfer_Matrix:: + ts(COMPLEX theta,const optical_constant& n0,const optical_constant& nt) + { + COMPLEX p0 = (COMPLEX)n0*cos_internal_angle(n0,theta); + COMPLEX pt = (COMPLEX)nt*cos_internal_angle(nt,theta); + return (2.*p0)/((As+Bs*pt)*p0+(Cs+Ds*pt)); + } + + COMPLEX + Thin_Film_Transfer_Matrix:: + tp(COMPLEX theta,const optical_constant& n0,const optical_constant& nt) + { + COMPLEX p0 = cos_internal_angle(n0,theta)/(COMPLEX)n0; + COMPLEX pt = cos_internal_angle(nt,theta)/(COMPLEX)nt; + return (2.*p0)/((Ap+Bp*pt)*p0+(Cp+Dp*pt))*(COMPLEX)n0/(COMPLEX)nt; + } + + // + // Multiple layers can be accounted for my multiplication of the transfer + // matrices for each layer... + // + Thin_Film_Transfer_Matrix + Thin_Film_Transfer_Matrix:: + operator*(Thin_Film_Transfer_Matrix a) + { + Thin_Film_Transfer_Matrix t; + t.As = As*a.As+Bs*a.Cs; + t.Ap = Ap*a.Ap+Bp*a.Cp; + t.Bs = As*a.Bs+Bs*a.Ds; + t.Bp = Ap*a.Bp+Bp*a.Dp; + t.Cs = Cs*a.As+Ds*a.Cs; + t.Cp = Cp*a.Ap+Dp*a.Cp; + t.Ds = Cs*a.Bs+Ds*a.Ds; + t.Dp = Cp*a.Bp+Dp*a.Dp; + return t; + } + + // + // The dielectric_stack class holds information about a stack of films... + // + + // + // Copy constructor... + // + /*dielectric_stack:: + dielectric_stack(const dielectric_stack& ds) + { + n=ds.n; + e.resize(0); + t.resize(0); + int i; + for (i=0;i=0; --i) + M = M*Transfer_Matrix(e[i].index(lambda),t[i],theta,lambda); + return M; + } + + // + // The net transfer matrix for light from the "substrate" side of the layer. + // + Thin_Film_Transfer_Matrix + dielectric_stack:: + matrix_backward(COMPLEX theta,double lambda) const + { + Thin_Film_Transfer_Matrix M; + int i; + for (i=0; i :"; + SCATMECH_output << prompt2; + + // Get response from user... + string response = SCATMECH_input.getstrline(); + + + if (response.size()==0) { + // If response is blank, return the default... + stack.read_stack_file(deflt); + } else if (response[0]==':') { + // If response begins with a colon, use the next word as a filename... + istringstream responsestr(response); + responsestr.get(); + string filename; + responsestr >> filename; + stack.read_stack_file(filename); + } else { + int i=1; + bool finished=false; + + string dielfunc; + double thickness=-1; + + while (!finished) { + + // Extract the dielectric function and thickness of the layer... + istringstream sresponse(response); + dielfunc.clear(); + thickness=-1; + sresponse >> dielfunc >> thickness; + + // If no dielectric function, assume finished... + if (dielfunc.size()==0) finished=true; + // Otherwise... + else { + // If thickness is less than zero, or no thickness given... + if (thickness<0.) throw SCATMECH_exception("Thickness must be greater than zero."); + + // If both the dielectric function and the thickness is good, then grow film... + stack.grow(dielectric_function(dielfunc),thickness); + + // Prompt user for next layer... + ostringstream message; + message << "Layer #" << ++i << ":"; + SCATMECH_output << message.str() << endl; + + // Get next response... + response = getstrline(sresponse); + string next; + istringstream temp(response); + temp >> next; + if (next.size()==0) { + response = SCATMECH_input.getstrline(); + } + } + } + } + ostringstream message; + stack.print(message); + SCATMECH_output << message.str(); + return stack; + } + catch(const exception& e) { + SCATMECH_output << e.what() << endl; + //instack_unwind(); + } + } + } + + void + dielectric_stack:: + read_stack_file(const string& film_file_name) + { + // If string pointer is NULL.. then assume no film... + if (film_file_name.size()==0) { + wash(); + return; + } + + // Open file stream... + string fname = find_file(film_file_name); + ifstream_with_comments film_file(fname.c_str()); + // and complain if it doesn't exist... + + if (!film_file) { + throw SCATMECH_exception("Cannot open file: " + film_file_name); + } + + do { + // Read values... + string material; + film_file >> material; + + // A blank line exits the routine... + if (material.size()==0) return; + + // Get the film thickness + double thickness=-1.; + film_file >> thickness; + if (thickness<0.) throw SCATMECH_exception("Thickness must be greater than zero."); + + // Deposit that film... + grow(dielectric_function(material),thickness); + + } while (!film_file.eof()); + } + + // + // The following prints out the status of the film... + // + void + dielectric_stack:: + print(ostream& os) const + { + os << endl << "# index thickness " << endl + << "-----------------------------------" << endl; + if (n==0) { + os << "No layers" << endl; + } else { + for (int i=0; i + void + ModelParameterSet(dielectric_stack& variable,const string& parameter,const string& value) + { + if (parameter.empty()) { + if (value.empty()) { + variable.wash(); + return; + } + istringstream iss(value); + string material; + iss >> material; + if (iss.fail()) throw SCATMECH_exception("Invalid parameter for dielectric_stack"); + double thickness; + iss >> thickness; + if (iss.fail()) { + variable.read_stack_file(material); + return; + } + variable.grow(material,thickness); + while (1) { + iss >> material; + if (iss.fail()) return; + iss >> thickness; + if (iss.fail()) throw SCATMECH_exception("Invalid thickness parameter for dielectric_stack"); + variable.grow(material,thickness); + } + } else if (parameter == "grow") { + istringstream iss(value); + string material; + double thickness; + iss >> material >> thickness; + if (iss.fail()) throw SCATMECH_exception("Invalid parameter for dielectric_stack"); + variable.grow(material,thickness); + } else if (parameter == "wash") { + variable.wash(); + } + } + + template <> + string + ModelParameterGet(dielectric_stack& variable,const string& parameter) + { + ostringstream oss; + oss << variable; + return oss.str(); + } + + template <> + void + ModelParameterAskUser(dielectric_stack& value,const string& prompt) + { + value = dielectric_stack::AskUser(prompt); + } + + void Register(const StackModel* x) + { + static bool Models_Registered = false; + if (!Models_Registered) { + Models_Registered=true; + + Register_Model(StackModel); + Register_Model(No_StackModel); + Register_Model(Stack_StackModel); + Register_Model(SingleFilm_StackModel); + Register_Model(DoubleFilm_StackModel); + Register_Model(GradedFilm_StackModel); + + } + } + + DEFINE_VIRTUAL_MODEL(StackModel,Model,"Abstract model for film stacks, allowing parametric variation of dielectric stacks"); + + DEFINE_MODEL(No_StackModel,StackModel,"An empty film stack"); + + DEFINE_MODEL(Stack_StackModel,StackModel,"StackModel that takes a dielectric stack"); + DEFINE_PARAMETER(Stack_StackModel,dielectric_stack,stack,"Film stack","",0xFF); + + DEFINE_MODEL(SingleFilm_StackModel,StackModel,"A single film"); + DEFINE_PARAMETER(SingleFilm_StackModel,dielectric_function,material,"Material","(1.5,0)",0xFF); + DEFINE_PARAMETER(SingleFilm_StackModel,double,thickness,"Thickness [um]","0.1",0xFF); + + DEFINE_MODEL(DoubleFilm_StackModel,StackModel,"A double film"); + DEFINE_PARAMETER(DoubleFilm_StackModel,dielectric_function,material1,"Material deposited first","(1.5,0)",0xFF); + DEFINE_PARAMETER(DoubleFilm_StackModel,double,thickness1,"Thickness of first film [um]","0.1",0xFF); + DEFINE_PARAMETER(DoubleFilm_StackModel,dielectric_function,material2,"Material deposited second","(1.5,0)",0xFF); + DEFINE_PARAMETER(DoubleFilm_StackModel,double,thickness2,"Thickness of second film[um]","0.1",0xFF); + + DEFINE_MODEL(GradedFilm_StackModel,StackModel,"A graded index film"); + DEFINE_PARAMETER(GradedFilm_StackModel,dielectric_function,start,"Material closest to substrate","(1.5,0)",0xFF); + DEFINE_PARAMETER(GradedFilm_StackModel,dielectric_function,end,"Material farthest from substrate","(1,0)",0xFF); + DEFINE_PARAMETER(GradedFilm_StackModel,double,thickness,"Thickness [um]","0.1",0xFF); + DEFINE_PARAMETER(GradedFilm_StackModel,int,steps,"Number of steps","5",0xFF); + +} // namespace SCATMECH + + + diff --git a/code/filmtran.h b/code/filmtran.h index f33c509..83cb8be 100644 --- a/code/filmtran.h +++ b/code/filmtran.h @@ -1,238 +1,445 @@ -//****************************************************************************** -//** SCATMECH: Polarized Light Scattering C++ Class Library -//** -//** File: filmtran.h -//** -//** 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 -//** -//** Version: 7.00 (January 2015) -//** -//****************************************************************************** -#ifndef SCATMECH_FILMTRAN_H -#define SCATMECH_FILMTRAN_H - -// -// In all of functions in this file, the angles were changed to be complex. -// (Version 3.02, TAG 7 AUG 2002) -// - -#include "mueller.h" -#include "dielfunc.h" - - -namespace SCATMECH { - - - class Thin_Film_Transfer_Matrix; - - // - // dielectric_stack is a class which holds all the relavant information about - // a dielctric stack. - // - class dielectric_stack - { - public: - // Routines to read a file containing the film stack information... - // The following member function was removed 10 OCT 2002 - // void read_stack_file(); - void read_stack_file(const std::string& filename); - - // Routine to ask user for a stack file... - // Added 10 OCT 2002 - static dielectric_stack AskUser(const std::string& query, const std::string& deflt=""); - - // Routine to print state of film... - void print(std::ostream& os) const; - - // Constructors... - dielectric_stack() { - n=0; - e.resize(0); - t.resize(0); - } - //dielectric_stack(const dielectric_stack& ds); - - // Routine to remove all the films... - void wash(); - - // Routine to add films (one at a time)... - void grow(const dielectric_function& epsilon, double thickness); - - const std::string& get_name() const { - return name; - } - friend std::ostream& operator<<(std::ostream& os,const dielectric_stack& ds); - - int get_n() const { - return n; - } - const std::vector& get_e() const { - return e; - } - const std::vector& get_t() const { - return t; - } - - double get_total_thickness() const { - double result=0; - for (int i=0; i e; - - // An array containing the list of thicknesses ... - std::vector t; - - // A routine used by the constructors and copiers... - void init(); - - std::string name; - - private: - // The thin film transfer matrices for the stack... - Thin_Film_Transfer_Matrix matrix_forward(COMPLEX theta,double lambda) const; - Thin_Film_Transfer_Matrix matrix_backward(COMPLEX theta,double lambda) const; - }; - - // - // Thin_Film_Transfer_Matrix is a private subclass which handles the - // reflection and transmission coefficients for a - // dielectric stack of films. - // - class Thin_Film_Transfer_Matrix { - public: - // The constuctor requires the incident angle... - Thin_Film_Transfer_Matrix(); - - // The matrix multiplication operator... - Thin_Film_Transfer_Matrix operator*(Thin_Film_Transfer_Matrix a); - - // The transfer matrix for a single layer... - - // The reflection and transmission coefficients... - // n0 and nt are the indices on the incident and tranmitting - // sides of the film. - COMPLEX rs(COMPLEX theta,const optical_constant& n0,const optical_constant& nt); - COMPLEX ts(COMPLEX theta,const optical_constant& n0,const optical_constant& nt); - COMPLEX rp(COMPLEX theta,const optical_constant& n0,const optical_constant& nt); - COMPLEX tp(COMPLEX theta,const optical_constant& n0,const optical_constant& nt); - friend Thin_Film_Transfer_Matrix - Transfer_Matrix(const optical_constant& n,double thickness, - COMPLEX theta, double lambda); - private: - // There are two matrices {A,B,C,D} for s and p polarized light. - COMPLEX As; - COMPLEX Ap; - COMPLEX Bs; - COMPLEX Bp; - COMPLEX Cs; - COMPLEX Cp; - COMPLEX Ds; - COMPLEX Dp; - - }; - - inline - Thin_Film_Transfer_Matrix:: - Thin_Film_Transfer_Matrix() - { - As = 1; - Bs = 0; - Cs = 0; - Ds = 1; - Ap = 1; - Bp = 0; - Cp = 0; - Dp = 1; - } - - // - // grow() adds a new layer above the current layers... - // - inline - void - dielectric_stack:: - grow(const dielectric_function& epsilon,double thickness) - { - e.push_back(epsilon); - t.push_back(thickness); - ++n; - } - - // - // wash() removes all the films... - // - inline - void - dielectric_stack:: - wash() - { - n=0; - e.resize(0); - t.resize(0); - } - - template <> - void ModelParameterSet(dielectric_stack& variable,const std::string& subparameter,const std::string& value); - template <> - std::string ModelParameterGet(dielectric_stack& variable,const std::string& subparameter); - template <> - void ModelParameterAskUser(dielectric_stack& variable,const std::string& prompt); - -} // namespace SCATMECH - - - -#endif +//****************************************************************************** +//** SCATMECH: Polarized Light Scattering C++ Class Library +//** +//** File: filmtran.h +//** +//** 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 +//** +//** Version: 7.00 (January 2015) +//** +//****************************************************************************** +#ifndef SCATMECH_FILMTRAN_H +#define SCATMECH_FILMTRAN_H + +// +// In all of functions in this file, the angles were changed to be complex. +// (Version 3.02, TAG 7 AUG 2002) +// + +#include "mueller.h" +#include "dielfunc.h" + +namespace SCATMECH { + + + class Thin_Film_Transfer_Matrix; + class StackModel; + + // + // dielectric_stack is a class which holds all the relavant information about + // a dielctric stack. + // + class dielectric_stack + { + public: + // Routines to read a file containing the film stack information... + // The following member function was removed 10 OCT 2002 + // void read_stack_file(); + void read_stack_file(const std::string& filename); + + // Routine to ask user for a stack file... + // Added 10 OCT 2002 + static dielectric_stack AskUser(const std::string& query, const std::string& deflt=""); + + // Routine to print state of film... + void print(std::ostream& os) const; + + // Constructors... + dielectric_stack() { + n=0; + e.resize(0); + t.resize(0); + } + //dielectric_stack(const dielectric_stack& ds); + + // Routine to remove all the films... + void wash(); + + // Routine to add films (one at a time)... + void grow(const dielectric_function& epsilon, double thickness); + + const std::string& get_name() const { + return name; + } + friend std::ostream& operator<<(std::ostream& os,const dielectric_stack& ds); + + int get_n() const { + return n; + } + const std::vector& get_e() const { + return e; + } + const std::vector& get_t() const { + return t; + } + + double get_total_thickness() const { + double result=0; + for (int i=0; i e; + + // An array containing the list of thicknesses ... + std::vector t; + + // A routine used by the constructors and copiers... + void init(); + + std::string name; + + private: + // The thin film transfer matrices for the stack... + Thin_Film_Transfer_Matrix matrix_forward(COMPLEX theta,double lambda) const; + Thin_Film_Transfer_Matrix matrix_backward(COMPLEX theta,double lambda) const; + }; + + // + // Thin_Film_Transfer_Matrix is a private subclass which handles the + // reflection and transmission coefficients for a + // dielectric stack of films. + // + class Thin_Film_Transfer_Matrix { + public: + // The constuctor requires the incident angle... + Thin_Film_Transfer_Matrix(); + + // The matrix multiplication operator... + Thin_Film_Transfer_Matrix operator*(Thin_Film_Transfer_Matrix a); + + // The transfer matrix for a single layer... + + // The reflection and transmission coefficients... + // n0 and nt are the indices on the incident and tranmitting + // sides of the film. + COMPLEX rs(COMPLEX theta,const optical_constant& n0,const optical_constant& nt); + COMPLEX ts(COMPLEX theta,const optical_constant& n0,const optical_constant& nt); + COMPLEX rp(COMPLEX theta,const optical_constant& n0,const optical_constant& nt); + COMPLEX tp(COMPLEX theta,const optical_constant& n0,const optical_constant& nt); + friend Thin_Film_Transfer_Matrix + Transfer_Matrix(const optical_constant& n,double thickness, + COMPLEX theta, double lambda); + private: + // There are two matrices {A,B,C,D} for s and p polarized light. + COMPLEX As; + COMPLEX Ap; + COMPLEX Bs; + COMPLEX Bp; + COMPLEX Cs; + COMPLEX Cp; + COMPLEX Ds; + COMPLEX Dp; + + }; + + inline + Thin_Film_Transfer_Matrix:: + Thin_Film_Transfer_Matrix() + { + As = 1; + Bs = 0; + Cs = 0; + Ds = 1; + Ap = 1; + Bp = 0; + Cp = 0; + Dp = 1; + } + + // + // grow() adds a new layer above the current layers... + // + inline + void + dielectric_stack:: + grow(const dielectric_function& epsilon,double thickness) + { + e.push_back(epsilon); + t.push_back(thickness); + ++n; + } + + // + // wash() removes all the films... + // + inline + void + dielectric_stack:: + wash() + { + n=0; + e.resize(0); + t.resize(0); + } + + template <> + void ModelParameterSet(dielectric_stack& variable,const std::string& subparameter,const std::string& value); + template <> + std::string ModelParameterGet(dielectric_stack& variable,const std::string& subparameter); + template <> + void ModelParameterAskUser(dielectric_stack& variable,const std::string& prompt); + + class StackModel : public Model, protected dielectric_stack + { + DECLARE_MODEL(); + public: + + const dielectric_stack& get_stack() {SETUP(); return *this;} + + // Routine to remove all the films... + void wash() {SETUP(); dielectric_stack::wash();} + + // Routine to add films (one at a time)... + void grow(const dielectric_function& epsilon, double thickness) { + SETUP(); dielectric_stack::grow(epsilon,thickness); + } + + int get_n() { + SETUP(); + return dielectric_stack::get_n(); + } + + const std::vector& get_e() { + SETUP(); return dielectric_stack::get_e(); + } + + const std::vector& get_t() { + SETUP(); + return dielectric_stack::get_t(); + } + + double get_total_thickness() { + SETUP(); + return dielectric_stack::get_total_thickness(); + } + + // The following reverses the order of the stack... + void reverse() {SETUP(); dielectric_stack::reverse();} + + // The reflection and transmission coefficients ... + // for light propagating from the top layer to the bottom layer... + // NOTE: theta is the "vacuum" angle: theta=asin(sin(internalangle)*n0); + COMPLEX rp12(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::rp12(theta,lambda,n0,nt);} + COMPLEX rs12(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::rs12(theta,lambda,n0,nt);} + COMPLEX tp12(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::tp12(theta,lambda,n0,nt);} + COMPLEX ts12(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::ts12(theta,lambda,n0,nt);} + + // for light propagating from the bottom layer to the top layer + // NOTE: theta is the "vacuum" angle: theta=asin(sin(internalangle)*n0); + COMPLEX rp21(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::rp21(theta,lambda,n0,nt);} + COMPLEX rs21(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::rs21(theta,lambda,n0,nt);} + COMPLEX tp21(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::tp21(theta,lambda,n0,nt);} + COMPLEX ts21(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::ts21(theta,lambda,n0,nt);} + + // Jones Matrix versions of reflection coefficients... + // NOTE: theta is the "vacuum" angle: theta=asin(sin(internalangle)*n0); + JonesMatrix r12(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::r12(theta,lambda,n0,nt);} + JonesMatrix t12(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::t12(theta,lambda,n0,nt);} + JonesMatrix r21(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::r21(theta,lambda,n0,nt);} + JonesMatrix t21(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::t21(theta,lambda,n0,nt);} + + // The reflection and transmission coefficients ... + // for light propagating from the top layer to the bottom layer... + // NOTE: angle is the internal angle of incidence + COMPLEX rp12i(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::rp12i(theta,lambda,n0,nt);} + COMPLEX rs12i(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::rs12i(theta,lambda,n0,nt);} + COMPLEX tp12i(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::tp12i(theta,lambda,n0,nt);} + COMPLEX ts12i(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::ts12i(theta,lambda,n0,nt);} + + // for light propagating from the bottom layer to the top layer + // NOTE: angle is the internal angle of incidence + COMPLEX rp21i(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::rp21i(theta,lambda,n0,nt);} + COMPLEX rs21i(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::rs21i(theta,lambda,n0,nt);} + COMPLEX tp21i(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::tp21i(theta,lambda,n0,nt);} + COMPLEX ts21i(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::ts21i(theta,lambda,n0,nt);} + + // Jones Matrix versions of reflection coefficients... + // NOTE: angle is the internal angle of incidence + JonesMatrix r12i(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::r12i(theta,lambda,n0,nt);} + JonesMatrix t12i(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::t12i(theta,lambda,n0,nt);} + JonesMatrix r21i(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::r21i(theta,lambda,n0,nt);} + JonesMatrix t21i(COMPLEX theta,double lambda,const dielectric_function& n0,const dielectric_function& nt) + {SETUP(); return dielectric_stack::t21i(theta,lambda,n0,nt);} + }; + + typedef Model_Ptr StackModel_Ptr; + + void Register(const StackModel* x); + + class No_StackModel : public StackModel + { + DECLARE_MODEL(); + protected: + void setup() {StackModel::setup(); wash();} + }; + + class Stack_StackModel : public StackModel + { + DECLARE_MODEL(); + DECLARE_PARAMETER(dielectric_stack,stack); + protected: + void setup() { + StackModel::setup(); + dielectric_stack::operator=(stack); + } + }; + + class SingleFilm_StackModel : public StackModel + { + DECLARE_MODEL(); + DECLARE_PARAMETER(dielectric_function,material); + DECLARE_PARAMETER(double,thickness); + protected: + void setup() { + StackModel::setup(); + wash(); + grow(material,thickness); + } + }; + + class DoubleFilm_StackModel : public StackModel + { + DECLARE_MODEL(); + DECLARE_PARAMETER(dielectric_function,material1); + DECLARE_PARAMETER(double,thickness1); + DECLARE_PARAMETER(dielectric_function,material2); + DECLARE_PARAMETER(double,thickness2); + protected: + void setup() { + StackModel::setup(); + wash(); + grow(material1,thickness1); + grow(material2,thickness2); + } + }; + + class GradedFilm_StackModel : public StackModel + { + DECLARE_MODEL(); + DECLARE_PARAMETER(dielectric_function,start); + DECLARE_PARAMETER(dielectric_function,end); + DECLARE_PARAMETER(double,thickness); + DECLARE_PARAMETER(int,steps); + protected: + void setup() { + StackModel::setup(); + wash(); + if (steps<2) error("steps<2"); + const Table::VECTOR& startT = start.get_n().get_table(); + const Table::VECTOR& endT = end.get_n().get_table(); + + int startn = startT.size(); + int endn = endT.size(); + double start1 = startT[0].first; + double start2 = startT[startn-1].first; + double end1 = endT[0].first; + double end2 = endT[endn-1].first; + int nlambda = startn>endn ? startn : endn; + double lambdamin = start1 < end1 ? start1 : end1; + double lambdamax = start2 > end2 ? start2 : end2; + for (int i=0;its12(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; + 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)); @@ -205,6 +204,36 @@ namespace SCATMECH { return fabs((1+c)/2.*q1+(1-c)/2.*q2); } + static double LegendreP(int i,double x) { + if (i<0) return 0; + if (i==0) return 1; + if (i==1) return x; + double P1=x; + double P2=1; + for (int j=2;j<=i;++j) { + double P = ((2.*j-1.)*x*P1-(j-1.)*P2)/j; + P2=P1; + P1=P; + } + return P1; + } + + double + Legendre_Phase_Function:: + f(double theta) + { + SETUP(); + double result=0; + double costheta=cos(theta); + if (c0!=0) result += c0*LegendreP(0,costheta); + if (c1!=0) result += 3*c1*LegendreP(1,costheta); + if (c2!=0) result += 5*c2*LegendreP(2,costheta); + if (c3!=0) result += 7*c3*LegendreP(3,costheta); + if (c4!=0) result += 9*c4*LegendreP(4,costheta); + if (c5!=0) result += 11*c5*LegendreP(5,costheta); + return result; + } + void Register(const First_Diffuse_BRDF_Model* x) { @@ -229,6 +258,7 @@ namespace SCATMECH { Register_Model(Reynolds_McCormick_Phase_Function); Register_Model(Double_Reynolds_McCormick_Phase_Function); Register_Model(Isotropic_Phase_Function); + Register_Model(Legendre_Phase_Function); Register_Model(Rayleigh_Phase_Function); Register_Model(Table_Phase_Function); } @@ -261,11 +291,14 @@ namespace SCATMECH { DEFINE_MODEL(Double_Reynolds_McCormick_Phase_Function,Phase_Function, "A Reynolds-McCormick phase function with a opposition effect."); + DEFINE_MODEL(Legendre_Phase_Function,Phase_Function, + "Phase function defined by a series of up to six Legendre polynomials (l=0 to l=5)"); + DEFINE_PARAMETER(First_Diffuse_BRDF_Model,double,l_scat,"Scatter mean free path [um]","1",0xFF); DEFINE_PTRPARAMETER(First_Diffuse_BRDF_Model,Phase_Function_Ptr,phase_function,"Phase Function","Henyey_Greenstein_Phase_Function",0xFF); - DEFINE_PARAMETER(First_Diffuse_BRDF_Model,dielectric_stack,stack,"Film stack on substrate","",0xFF); + DEFINE_PTRPARAMETER(First_Diffuse_BRDF_Model,StackModel_Ptr,stack,"Film stack on substrate","No_StackModel",0xFF); DEFINE_PARAMETER(Henyey_Greenstein_Phase_Function,double,g,"Asymmetry parameter (g)","0.01",0xFF); @@ -281,6 +314,12 @@ namespace SCATMECH { DEFINE_PARAMETER(Double_Reynolds_McCormick_Phase_Function,double,alpha,"Peaking factor (alpha)","0.5",0xFF); DEFINE_PARAMETER(Double_Reynolds_McCormick_Phase_Function,double,c,"Mixing Factor (c)","0.1",0xFF); + DEFINE_PARAMETER(Legendre_Phase_Function,double,c0,"Legendre coefficient (l=0)","1",0xFF); + DEFINE_PARAMETER(Legendre_Phase_Function,double,c1,"Legendre coefficient (l=1)","0.65",0xFF); + DEFINE_PARAMETER(Legendre_Phase_Function,double,c2,"Legendre coefficient (l=2)","0.42",0xFF); + DEFINE_PARAMETER(Legendre_Phase_Function,double,c3,"Legendre coefficient (l=3)","0",0xFF); + DEFINE_PARAMETER(Legendre_Phase_Function,double,c4,"Legendre coefficient (l=4)","0",0xFF); + DEFINE_PARAMETER(Legendre_Phase_Function,double,c5,"Legendre coefficient (l=5)","0",0xFF); } // namespace SCATMECH diff --git a/code/firstdiffuse.h b/code/firstdiffuse.h index 9285a7e..d799cdc 100644 --- a/code/firstdiffuse.h +++ b/code/firstdiffuse.h @@ -40,7 +40,7 @@ namespace SCATMECH { DECLARE_MODEL(); DECLARE_PARAMETER(double,l_scat); DECLARE_PARAMETER(Phase_Function_Ptr,phase_function); - DECLARE_PARAMETER(dielectric_stack,stack); + DECLARE_PARAMETER(StackModel_Ptr,stack); protected: virtual void setup(); @@ -118,6 +118,20 @@ namespace SCATMECH { double K; }; + class Legendre_Phase_Function: public Phase_Function + { + public: + virtual double f(double theta); + DECLARE_MODEL(); + DECLARE_PARAMETER(double,c0); + DECLARE_PARAMETER(double,c1); + DECLARE_PARAMETER(double,c2); + DECLARE_PARAMETER(double,c3); + DECLARE_PARAMETER(double,c4); + DECLARE_PARAMETER(double,c5); + }; + + void Register(const First_Diffuse_BRDF_Model*); void Register(const Phase_Function* x); diff --git a/code/flake.cpp b/code/flake.cpp index 03a08cc..7ea3699 100644 --- a/code/flake.cpp +++ b/code/flake.cpp @@ -101,8 +101,8 @@ namespace SCATMECH { // The reflection coefficients are obtained from the // dielectric stack... - COMPLEX rsiota = stack.rs12(iota_external,lambda,overcoat,substrate); - COMPLEX rpiota = stack.rp12(iota_external,lambda,overcoat,substrate); + COMPLEX rsiota = stack->rs12(iota_external,lambda,overcoat,substrate); + COMPLEX rpiota = stack->rp12(iota_external,lambda,overcoat,substrate); // The local surface normal... Vector _nhat = nhat(thetai_internal,thetas_internal,phis); @@ -135,12 +135,12 @@ namespace SCATMECH { } // The transmission matrices into and out of the overcoat... - JonesMatrix jti(overcoat_films.tp12(thetai,lambda,vacuum,epsilon), - overcoat_films.ts12(thetai,lambda,vacuum,epsilon), + JonesMatrix jti(overcoat_films->tp12(thetai,lambda,vacuum,epsilon), + overcoat_films->ts12(thetai,lambda,vacuum,epsilon), (COMPLEX)(0.),(COMPLEX)(0.)); - JonesMatrix jts(overcoat_films.tp12(thetas,lambda,vacuum,epsilon), - overcoat_films.ts12(thetas,lambda,vacuum,epsilon), + JonesMatrix jts(overcoat_films->tp12(thetas,lambda,vacuum,epsilon), + overcoat_films->ts12(thetas,lambda,vacuum,epsilon), (COMPLEX)(0.),(COMPLEX)(0.)); // Total Jones matrix @@ -165,7 +165,7 @@ namespace SCATMECH { DEFINE_PARAMETER(Subsurface_Facet_BRDF_Model,dielectric_function,overcoat,"Overcoat","(1.59,0)",0xFF); - DEFINE_PARAMETER(Subsurface_Facet_BRDF_Model,dielectric_stack,overcoat_films,"Overcoat films","",0xFF); + DEFINE_PTRPARAMETER(Subsurface_Facet_BRDF_Model,StackModel_Ptr,overcoat_films,"Overcoat films","No_StackModel",0xFF); } // namespace SCATMECH diff --git a/code/flake.h b/code/flake.h index 458a930..ab6e944 100644 --- a/code/flake.h +++ b/code/flake.h @@ -35,7 +35,7 @@ namespace SCATMECH { DECLARE_PARAMETER(dielectric_function,overcoat); // Any films which exist above the overcoat... - DECLARE_PARAMETER(dielectric_stack,overcoat_films); + DECLARE_PARAMETER(StackModel_Ptr,overcoat_films); protected: diff --git a/code/grating.cpp b/code/grating.cpp index ada146a..32fbd5b 100644 --- a/code/grating.cpp +++ b/code/grating.cpp @@ -258,7 +258,10 @@ namespace SCATMECH { } return false; } - } + bool tdiff(StackModel& a, StackModel& b) { + return tdiff(a.get_stack(),b.get_stack()); + } + } void Dielectric_Stack_Grating::setup() { @@ -273,51 +276,51 @@ namespace SCATMECH { materialmuz.clear(); thickness.clear(); - int n = stackepsx.get_n(); + int n = stackepsx->get_n(); - if (stackepsx.get_n()==0 && stackepsy.get_n()==0 && stackepsz.get_n()==0 && - stackmux.get_n()==0 && stackmuy.get_n()==0 && stackmuz.get_n()==0) { // This is an empty stack + if (stackepsx->get_n()==0 && stackepsy->get_n()==0 && stackepsz->get_n()==0 && + stackmux->get_n()==0 && stackmuy->get_n()==0 && stackmuz->get_n()==0) { // This is an empty stack return; } - if (stackepsy.get_n()==0 && stackepsz.get_n()==0 && - stackmux.get_n()==0 && stackmuy.get_n()==0 && stackmuz.get_n()==0) { // This is an isotropic stack - materialx.resize(stackepsx.get_n()); - thickness.resize(stackepsx.get_n()); - position.resize(stackepsx.get_n()); - for (int i=0; iget_n()==0 && stackepsz->get_n()==0 && + stackmux->get_n()==0 && stackmuy->get_n()==0 && stackmuz->get_n()==0) { // This is an isotropic stack + materialx.resize(stackepsx->get_n()); + thickness.resize(stackepsx->get_n()); + position.resize(stackepsx->get_n()); + for (int i=0; iget_n(); ++i) { int ii = n-i-1; - thickness[i]=stackepsx.get_t()[ii]; + thickness[i]=stackepsx->get_t()[ii]; position[i].resize(2); materialx[i].resize(2); position[i][0] = 0; position[i][1] = period; - materialx[i][0] = epsilon(stackepsx.get_e()[ii]); + materialx[i][0] = epsilon(stackepsx->get_e()[ii]); materialx[i][1] = materialx[i][0]; } materialz=materialy=materialx; return; } - if (stackepsy.get_n()==0 && stackepsz.get_n()==0 && - stackmuy.get_n()==0 && stackmuz.get_n()==0) { // This is a isotropic magnetic material - if (tdiff(stackepsx,stackmux)) error("Layers inconsistent"); + if (stackepsy->get_n()==0 && stackepsz->get_n()==0 && + stackmuy->get_n()==0 && stackmuz->get_n()==0) { // This is a isotropic magnetic material + if (tdiff(*stackepsx,*stackmux)) error("Layers inconsistent"); magnetic = true; // TODO: code for isotropic magnetic material - materialx.resize(stackepsx.get_n()); - materialmux.resize(stackepsx.get_n()); - thickness.resize(stackepsx.get_n()); - position.resize(stackepsx.get_n()); - for (int i=0; iget_n()); + materialmux.resize(stackepsx->get_n()); + thickness.resize(stackepsx->get_n()); + position.resize(stackepsx->get_n()); + for (int i=0; iget_n(); ++i) { int ii = n-i-1; - thickness[i]=stackepsx.get_t()[ii]; + thickness[i]=stackepsx->get_t()[ii]; position[i].resize(2); materialx[i].resize(2); materialmux[i].resize(2); position[i][0] = 0; position[i][1] = period; - materialx[i][0] = epsilon(stackepsx.get_e()[ii]); + materialx[i][0] = epsilon(stackepsx->get_e()[ii]); materialx[i][1] = materialx[i][0]; - materialmux[i][0] = epsilon(stackmux.get_e()[ii]); + materialmux[i][0] = epsilon(stackmux->get_e()[ii]); materialmux[i][1] = materialmux[i][0]; } materialz=materialy=materialx; @@ -325,48 +328,48 @@ namespace SCATMECH { return; } - if (stackmux.get_n()==0 && stackmuy.get_n()==0 && stackmuz.get_n()==0) { // This is an anisotropic, but nonmagnetic, stack - if (tdiff(stackepsx,stackepsy) || tdiff(stackepsx,stackepsz)) error("Layers inconsistent"); + if (stackmux->get_n()==0 && stackmuy->get_n()==0 && stackmuz->get_n()==0) { // This is an anisotropic, but nonmagnetic, stack + if (tdiff(*stackepsx,*stackepsy) || tdiff(*stackepsx,*stackepsz)) error("Layers inconsistent"); anisotropic = true; - materialx.resize(stackepsx.get_n()); - materialy.resize(stackepsx.get_n()); - materialz.resize(stackepsx.get_n()); - thickness.resize(stackepsx.get_n()); - position.resize(stackepsx.get_n()); - for (int i=0; iget_n()); + materialy.resize(stackepsx->get_n()); + materialz.resize(stackepsx->get_n()); + thickness.resize(stackepsx->get_n()); + position.resize(stackepsx->get_n()); + for (int i=0; iget_n(); ++i) { int ii = n-i-1; - thickness[i]=stackepsx.get_t()[ii]; + thickness[i]=stackepsx->get_t()[ii]; position[i].resize(2); materialx[i].resize(2); materialy[i].resize(2); materialz[i].resize(2); position[i][0] = 0; position[i][1] = period; - materialx[i][0] = epsilon(stackepsx.get_e()[ii]); + materialx[i][0] = epsilon(stackepsx->get_e()[ii]); materialx[i][1] = materialx[i][0]; - materialy[i][0] = epsilon(stackepsy.get_e()[ii]); + materialy[i][0] = epsilon(stackepsy->get_e()[ii]); materialy[i][1] = materialy[i][0]; - materialz[i][0] = epsilon(stackepsz.get_e()[ii]); + materialz[i][0] = epsilon(stackepsz->get_e()[ii]); materialz[i][1] = materialz[i][0]; } return; } - if (tdiff(stackepsx,stackepsy) || tdiff(stackepsx,stackepsz) || tdiff(stackepsx,stackmux) || tdiff(stackepsx,stackmuy) || tdiff(stackepsx,stackmuz)) error("Layersinconsistent"); + if (tdiff(*stackepsx,*stackepsy) || tdiff(*stackepsx,*stackepsz) || tdiff(*stackepsx,*stackmux) || tdiff(*stackepsx,*stackmuy) || tdiff(*stackepsx,*stackmuz)) error("Layers inconsistent"); magnetic = true; anisotropic = true; - materialx.resize(stackepsx.get_n()); - materialy.resize(stackepsx.get_n()); - materialz.resize(stackepsx.get_n()); - materialmux.resize(stackepsx.get_n()); - materialmuy.resize(stackepsx.get_n()); - materialmuz.resize(stackepsx.get_n()); - thickness.resize(stackepsx.get_n()); - position.resize(stackepsx.get_n()); - for (int i=0; iget_n()); + materialy.resize(stackepsx->get_n()); + materialz.resize(stackepsx->get_n()); + materialmux.resize(stackepsx->get_n()); + materialmuy.resize(stackepsx->get_n()); + materialmuz.resize(stackepsx->get_n()); + thickness.resize(stackepsx->get_n()); + position.resize(stackepsx->get_n()); + for (int i=0; iget_n(); ++i) { int ii = n-i-1; - thickness[i]=stackepsx.get_t()[ii]; + thickness[i]=stackepsx->get_t()[ii]; position[i].resize(2); materialx[i].resize(2); materialy[i].resize(2); @@ -376,17 +379,17 @@ namespace SCATMECH { materialmuz[i].resize(2); position[i][0] = 0; position[i][1] = period; - materialx[i][0] = epsilon(stackepsx.get_e()[ii]); + materialx[i][0] = epsilon(stackepsx->get_e()[ii]); materialx[i][1] = materialx[i][0]; - materialy[i][0] = epsilon(stackepsy.get_e()[ii]); + materialy[i][0] = epsilon(stackepsy->get_e()[ii]); materialy[i][1] = materialy[i][0]; - materialz[i][0] = epsilon(stackepsz.get_e()[ii]); + materialz[i][0] = epsilon(stackepsz->get_e()[ii]); materialz[i][1] = materialz[i][0]; - materialmux[i][0] = epsilon(stackmux.get_e()[ii]); + materialmux[i][0] = epsilon(stackmux->get_e()[ii]); materialmux[i][1] = materialmux[i][0]; - materialmuy[i][0] = epsilon(stackmuy.get_e()[ii]); + materialmuy[i][0] = epsilon(stackmuy->get_e()[ii]); materialmuy[i][1] = materialmuy[i][0]; - materialmuz[i][0] = epsilon(stackmuz.get_e()[ii]); + materialmuz[i][0] = epsilon(stackmuz->get_e()[ii]); materialmuz[i][1] = materialmuz[i][0]; } return; @@ -1576,12 +1579,12 @@ namespace SCATMECH { DEFINE_PARAMETER(Generic_Grating,int,nlayers,"Approximate number of levels","20",0xFF); DEFINE_MODEL(Dielectric_Stack_Grating,Grating,"A grating with zero layers."); - DEFINE_PARAMETER(Dielectric_Stack_Grating,dielectric_stack,stackepsx,"Stack of films","",0xFF); - DEFINE_PARAMETER(Dielectric_Stack_Grating,dielectric_stack,stackepsy,"Stack of films","",0xFF); - DEFINE_PARAMETER(Dielectric_Stack_Grating,dielectric_stack,stackepsz,"Stack of films","",0xFF); - DEFINE_PARAMETER(Dielectric_Stack_Grating,dielectric_stack,stackmux,"Stack of films","",0xFF); - DEFINE_PARAMETER(Dielectric_Stack_Grating,dielectric_stack,stackmuy,"Stack of films","",0xFF); - DEFINE_PARAMETER(Dielectric_Stack_Grating,dielectric_stack,stackmuz,"Stack of films","",0xFF); + DEFINE_PTRPARAMETER(Dielectric_Stack_Grating,StackModel_Ptr,stackepsx,"Stack of films","No_StackModel",0xFF); + DEFINE_PTRPARAMETER(Dielectric_Stack_Grating,StackModel_Ptr,stackepsy,"Stack of films","No_StackModel",0xFF); + DEFINE_PTRPARAMETER(Dielectric_Stack_Grating,StackModel_Ptr,stackepsz,"Stack of films","No_StackModel",0xFF); + DEFINE_PTRPARAMETER(Dielectric_Stack_Grating,StackModel_Ptr,stackmux,"Stack of films","No_StackModel",0xFF); + DEFINE_PTRPARAMETER(Dielectric_Stack_Grating,StackModel_Ptr,stackmuy,"Stack of films","No_StackModel",0xFF); + DEFINE_PTRPARAMETER(Dielectric_Stack_Grating,StackModel_Ptr,stackmuz,"Stack of films","No_StackModel",0xFF); DEFINE_MODEL(Overlaid_Grating,Grating,"One grating on top of another, with a specified offset"); DEFINE_PTRPARAMETER(Overlaid_Grating,Grating_Ptr,top,"Top grating","Single_Line_Grating",0xFF); diff --git a/code/grating.h b/code/grating.h index ee7f19e..540ceb5 100644 --- a/code/grating.h +++ b/code/grating.h @@ -269,12 +269,12 @@ namespace SCATMECH { private: DECLARE_MODEL(); - DECLARE_PARAMETER(dielectric_stack,stackepsx); - DECLARE_PARAMETER(dielectric_stack,stackepsy); - DECLARE_PARAMETER(dielectric_stack,stackepsz); - DECLARE_PARAMETER(dielectric_stack,stackmux); - DECLARE_PARAMETER(dielectric_stack,stackmuy); - DECLARE_PARAMETER(dielectric_stack,stackmuz); + DECLARE_PARAMETER(StackModel_Ptr,stackepsx); + DECLARE_PARAMETER(StackModel_Ptr,stackepsy); + DECLARE_PARAMETER(StackModel_Ptr,stackepsz); + DECLARE_PARAMETER(StackModel_Ptr,stackmux); + DECLARE_PARAMETER(StackModel_Ptr,stackmuy); + DECLARE_PARAMETER(StackModel_Ptr,stackmuz); }; class Single_Line_Grating : public Grating { diff --git a/code/inherit.h b/code/inherit.h index 2b02ef6..a74a9bb 100644 --- a/code/inherit.h +++ b/code/inherit.h @@ -849,7 +849,7 @@ namespace SCATMECH { TYPE &_parameter = ((MODEL*)model)->*parameter; SCATMECH_output << description; _parameter = Get_Model_Ptr(); - _parameter->AskUser(); + _parameter->Model::AskUser(); model->set_recalc(recalclevel); } diff --git a/code/jmatrix.cpp b/code/jmatrix.cpp index 4d60bbf..2071afc 100644 --- a/code/jmatrix.cpp +++ b/code/jmatrix.cpp @@ -196,21 +196,21 @@ namespace SCATMECH { double norm = sqrt(sqr(k.x)+sqr(k.y)); double kx = k.x/norm; double ky = k.y/norm; - if (k.z>0) + if (k.z>0) x = -p*kx-s*ky; else x = p*kx-s*ky; } else { x = Vector(1,0,0); } - y = perpto(k,x); + y = perpto(x,k); // {y,x,k} is right-handed. } void GetBasisVectorsParPerp(const Vector& kin,const Vector& kout, Vector& perp, Vector& parin, Vector& parout) { perp = perpto(kin,kout); - parin = perpto(perp,kin); - parout = perpto(perp,kout); + parin = perpto(kin,perp); + parout = perpto(kout,perp); } JonesMatrix GetJonesRotator(const Vector& xo, const Vector& yo, const Vector& xi, const Vector& yi) diff --git a/code/matrixmath.h b/code/matrixmath.h index e7abfec..9a3046a 100644 --- a/code/matrixmath.h +++ b/code/matrixmath.h @@ -126,6 +126,13 @@ namespace SCATMECH { #endif return p[i-1]; } + const T& operator()(int i) const + { + #ifdef _DEBUG + //checkbounds(i); + #endif + return p[i-1]; + } T& operator()(int i,int j) { #ifdef _DEBUG @@ -133,6 +140,13 @@ namespace SCATMECH { #endif return p[(i-1)+step[0]*(j-1)]; } + const T& operator()(int i,int j) const + { + #ifdef _DEBUG + checkbounds(i,j); + #endif + return p[(i-1)+step[0]*(j-1)]; + } T& operator()(int i,int j,int k) { #ifdef _DEBUG @@ -140,6 +154,13 @@ namespace SCATMECH { #endif return p[(i-1)+step[0]*(j-1)+step[1]*(k-1)]; } + const T& operator()(int i,int j,int k) const + { + #ifdef _DEBUG + checkbounds(i,j,k); + #endif + return p[(i-1)+step[0]*(j-1)+step[1]*(k-1)]; + } T& operator()(int i,int j,int k,int l) { #ifdef _DEBUG @@ -147,6 +168,13 @@ namespace SCATMECH { #endif return p[(i-1)+step[0]*(j-1)+step[1]*(k-1)+step[2]*(l-1)]; } + const T& operator()(int i,int j,int k,int l) const + { + #ifdef _DEBUG + checkbounds(i,j,k,l); + #endif + return p[(i-1)+step[0]*(j-1)+step[1]*(k-1)+step[2]*(l-1)]; + } T& operator()(int i,int j,int k,int l,int I) { #ifdef _DEBUG @@ -154,6 +182,13 @@ namespace SCATMECH { #endif return p[(i-1)+step[0]*(j-1)+step[1]*(k-1)+step[2]*(l-1)+step[3]*(I-1)]; } + const T& operator()(int i,int j,int k,int l,int I) const + { + #ifdef _DEBUG + checkbounds(i,j,k,l,I); + #endif + return p[(i-1)+step[0]*(j-1)+step[1]*(k-1)+step[2]*(l-1)+step[3]*(I-1)]; + } T& operator()(int i,int j,int k,int l,int I,int J) { #ifdef _DEBUG @@ -161,6 +196,13 @@ namespace SCATMECH { #endif return p[(i-1)+step[0]*(j-1)+step[1]*(k-1)+step[2]*(l-1)+step[3]*(I-1)+step[4]*(J-1)]; } + const T& operator()(int i,int j,int k,int l,int I,int J) const + { + #ifdef _DEBUG + checkbounds(i,j,k,k,I,J); + #endif + return p[(i-1)+step[0]*(j-1)+step[1]*(k-1)+step[2]*(l-1)+step[3]*(I-1)+step[4]*(J-1)]; + } T& operator()(int i,int j,int k,int l,int I,int J,int K) { #ifdef _DEBUG @@ -168,6 +210,14 @@ namespace SCATMECH { #endif return p[(i-1)+step[0]*(j-1)+step[1]*(k-1)+step[2]*(l-1)+step[3]*(I-1)+step[4]*(J-1)+step[5]*(K-1)]; } + const T& operator()(int i,int j,int k,int l,int I,int J,int K) const + { + #ifdef _DEBUG + checkbounds(i,j,k,l,I,J,K); + #endif + return p[(i-1)+step[0]*(j-1)+step[1]*(k-1)+step[2]*(l-1)+step[3]*(I-1)+step[4]*(J-1)+step[5]*(K-1)]; + } + T& operator()(int i,int j,int k,int l,int I,int J,int K,int L) { #ifdef _DEBUG @@ -175,13 +225,23 @@ namespace SCATMECH { #endif return p[(i-1)+step[0]*(j-1)+step[1]*(k-1)+step[2]*(l-1)+step[3]*(I-1)+step[4]*(J-1)+step[5]*(K-1)+step[6]*(L-1)]; } + const T& operator()(int i,int j,int k,int l,int I,int J,int K,int L) const + { + #ifdef _DEBUG + checkbounds(i,j,k,l,I,J,K,L); + #endif + return p[(i-1)+step[0]*(j-1)+step[1]*(k-1)+step[2]*(l-1)+step[3]*(I-1)+step[4]*(J-1)+step[5]*(K-1)+step[6]*(L-1)]; + } T& operator[](int i) { return p[i]; } + const T& operator[](int i) const { + return p[i]; + } #ifdef _DEBUG - void checkbounds(int i,int j=1,int k=1,int l=1,int I=1,int J=1,int K=1,int L=1) { + void checkbounds(int i,int j=1,int k=1,int l=1,int I=1,int J=1,int K=1,int L=1) const { if (i>dims[0] || j>dims[1] || k>dims[2] || l>dims[3] || I>dims[4] || J>dims[5] || K>dims[6] || L>dims[7]) throw SCATMECH_exception("Array out of bounds in FARRAY"); } diff --git a/code/matrixmath2.cpp b/code/matrixmath2.cpp index 10c1db8..0a1ef96 100644 --- a/code/matrixmath2.cpp +++ b/code/matrixmath2.cpp @@ -190,8 +190,9 @@ namespace SCATMECH { //C //C SET THE MAXIMUM NUMBER OF ITERATIONS. //C + int MAXIT = 30; + //C***FIRST EXECUTABLE STATEMENT CSVDC - int MAXIT = 30; //C //C DETERMINE WHAT IS TO BE COMPUTED. //C @@ -391,7 +392,6 @@ namespace SCATMECH { //C MM = M; ITER = 0; - MAXIT = M; // THOM KLUDGE! //C //C QUIT IF ALL THE SINGULAR VALUES HAVE BEEN FOUND. //C diff --git a/code/models.cpp b/code/models.cpp index 1fe1bb9..c026746 100644 --- a/code/models.cpp +++ b/code/models.cpp @@ -20,6 +20,7 @@ #include "rcw.h" #include "crossrcw.h" #include "reflectance.h" +#include "filmtran.h" using namespace std; @@ -41,6 +42,7 @@ namespace SCATMECH { Register((CrossGrating*)NULL); Register_Model(RCW_Model); Register_Model(CrossRCW_Model); + Register((StackModel*)NULL); } } diff --git a/code/mueller.cpp b/code/mueller.cpp index 5727dac..15c20e4 100644 --- a/code/mueller.cpp +++ b/code/mueller.cpp @@ -1,4 +1,4 @@ -//****************************************************************************** +//********ei********************************************************************** //** SCATMECH: Polarized Light Scattering C++ Class Library //** //** File: mueller.cpp @@ -113,7 +113,7 @@ namespace SCATMECH { result[1][3] = (-imag(j2SPj1PP) + imag(j2SSj1PS) + imag(j2PPj1SP) - imag(j2PSj1SS))/2.; result[2][3] = ( imag(j2SSj1PP) + imag(j2SPj1PS) - imag(j2PSj1SP) - imag(j2PPj1SS))/2.; result[3][3] = ( real(j2SSj1PP) - real(j2SPj1PS) - real(j2PSj1SP) + real(j2PPj1SS))/2.; - + return result; } @@ -369,6 +369,17 @@ namespace SCATMECH { return m; } + MuellerMatrix + MuellerDiagonal(double m00,double m11, double m22, double m33) + { + MuellerMatrix m=MuellerZero(); + m[0][0]=m00; + m[1][1]=m11; + m[2][2]=m22; + m[3][3]=m33; + return m; + } + MuellerMatrix MuellerPartialLinearPolarizer(double tmax, double tmin, double angle) { @@ -591,15 +602,21 @@ namespace SCATMECH { } } - MuellerMatrix MuellerMatrix::Closest_NonDepolarizing() const + MuellerMatrix MuellerMatrix::Closest_NonDepolarizing(int rank) const { MuellerMatrix M1,M2,M3,M4; Cloude_Decomposition(M1,M2,M3,M4); MuellerMatrix &m1=M1,&m2=M2,&m3=M3,&m4=M4; - if (m4[0][0]>m3[0][0]) m3 = m4; - if (m3[0][0]>m2[0][0]) m2 = m3; - if (m2[0][0]>m1[0][0]) m1 = m2; - return m1; + if (m4[0][0]>m3[0][0]) swap(m3,m4); + if (m3[0][0]>m2[0][0]) swap(m2,m3); + if (m2[0][0]>m1[0][0]) swap(m1,m2); + if (rank==1) return m1; + if (m4[0][0]>m3[0][0]) swap(m3,m4); + if (m3[0][0]>m2[0][0]) swap(m2,m3); + if (rank==2) return m2; + if (m4[0][0]>m3[0][0]) swap(m3,m4); + if (rank==3) return m3; + return m4; } MuellerMatrix MuellerMatrix::log() const @@ -901,6 +918,12 @@ namespace SCATMECH { } } } + //for (int j=0;j<4;++j) { + // if (depolarizer[j][j]<0) { + // for (int i=0;i<4;++i) depolarizer[i][j] = -depolarizer[i][j]; + // for (int i=0;i<4;++i) retarder[j][i] = -retarder[j][i]; + // } + //} } MuellerMatrix diff --git a/code/mueller.h b/code/mueller.h index e10519b..de1aff1 100644 --- a/code/mueller.h +++ b/code/mueller.h @@ -637,7 +637,7 @@ namespace SCATMECH { /// /// Returns the closest Mueller matrix that is non-depolarizing /// - MuellerMatrix Closest_NonDepolarizing() const; + MuellerMatrix Closest_NonDepolarizing(int rank=1) const; /// /// The Mueller matrix normalized, M/M[0][0] @@ -1093,6 +1093,10 @@ namespace SCATMECH { double attenuation=1., ///< An overal scaling factor double depolarization=1. ///< The amount of depolarization, 0<=depolarization<=1 ); + // + // Returns a diagonal Mueller matrix + // + MuellerMatrix MuellerDiagonal(double m00,double m11, double m22, double m33); /// /// Returns a partiall linear polarizer /// diff --git a/code/nsphere.cpp b/code/nsphere.cpp index 9ab3df5..9284bd6 100644 --- a/code/nsphere.cpp +++ b/code/nsphere.cpp @@ -42,7 +42,7 @@ namespace SCATMECH { { SphericalScatterer::setup(); - int nlayers = stack.get_n(); + int nlayers = stack->get_n(); int r = nlayers+1; CFARRAY m(r,1); CFARRAY n(r+1,1); @@ -56,7 +56,7 @@ namespace SCATMECH { } else if (s==r+1) { n(s) = medium.index(lambda); } else { - n(s) = stack.get_e()[s-2].index(lambda); + n(s) = stack->get_e()[s-2].index(lambda); } } for (int s=1; s<=r; ++s) { @@ -70,11 +70,11 @@ namespace SCATMECH { if (s==1) { R(s)=SphericalScatterer::radius; for (int t=0; tget_t()[t]; } if (R(s)<0.) throw SCATMECH_exception("Thickness of stack exceeds particle radius."); } else { - R(s)=R(s-1)+stack.get_t()[s-2]; + R(s)=R(s-1)+stack->get_t()[s-2]; } y(s) = k(s)*R(s); x(s) = 2*pi/lambda*R(s); @@ -238,7 +238,7 @@ namespace SCATMECH { DEFINE_MODEL(MultilayerCoatedMieScatterer,SphericalScatterer, "Scattering from a sphere with any number of coatings."); - DEFINE_PARAMETER(MultilayerCoatedMieScatterer,dielectric_stack,stack,"Coating stack on core sphere","",0xFF); + DEFINE_PTRPARAMETER(MultilayerCoatedMieScatterer,StackModel_Ptr,stack,"Coating stack on core sphere","No_StackModel",0xFF); } // namespace SCATMECH; diff --git a/code/nsphere.h b/code/nsphere.h index c73475d..a2db534 100644 --- a/code/nsphere.h +++ b/code/nsphere.h @@ -41,7 +41,7 @@ namespace SCATMECH { public: DECLARE_MODEL(); - DECLARE_PARAMETER(dielectric_stack,stack); + DECLARE_PARAMETER(StackModel_Ptr,stack); private: diff --git a/code/rayinst.cpp b/code/rayinst.cpp index d4c4dd3..b164b83 100644 --- a/code/rayinst.cpp +++ b/code/rayinst.cpp @@ -58,10 +58,10 @@ namespace SCATMECH { double n_minus_1 = air.n(lambda)-1.; // Reflection coefficients ... - complex rss = stack.rs12(thetas,lambda,vacuum,substrate); - complex rsi = stack.rs12(thetai,lambda,vacuum,substrate); - complex rps = stack.rp12(thetas,lambda,vacuum,substrate); - complex rpi = stack.rp12(thetai,lambda,vacuum,substrate); + complex rss = stack->rs12(thetas,lambda,vacuum,substrate); + complex rsi = stack->rs12(thetai,lambda,vacuum,substrate); + complex rps = stack->rp12(thetas,lambda,vacuum,substrate); + complex rpi = stack->rp12(thetai,lambda,vacuum,substrate); // Angle subtended by viewing and scattering directions... double subtendi = subtend(-thetai,0.,thetas,phis); @@ -114,7 +114,7 @@ namespace SCATMECH { DEFINE_MODEL(Rayleigh_Instrument_BRDF_Model,Instrument_BRDF_Model, "Effective BRDF due to Rayleigh scatter by the air surrounding a smooth sample."); - DEFINE_PARAMETER(Rayleigh_Instrument_BRDF_Model,dielectric_stack,stack,"Film stack on substrate","",0xFF); + DEFINE_PTRPARAMETER(Rayleigh_Instrument_BRDF_Model,StackModel_Ptr,stack,"Film stack on substrate","No_StackModel",0xFF); DEFINE_PARAMETER(Rayleigh_Instrument_BRDF_Model,double,field_of_view,"Field of view of detector [um]","1000",0xFF); diff --git a/code/rayinst.h b/code/rayinst.h index 31841fd..c1246f0 100644 --- a/code/rayinst.h +++ b/code/rayinst.h @@ -43,7 +43,7 @@ namespace SCATMECH { DECLARE_PARAMETER(double,number_density); // Any dielectric films on the substrate... - DECLARE_PARAMETER(dielectric_stack,stack); + DECLARE_PARAMETER(StackModel_Ptr,stack); protected: // The Mueller matrix BRDF for scattering... diff --git a/code/raystack.cpp b/code/raystack.cpp index c8d3815..69c382d 100644 --- a/code/raystack.cpp +++ b/code/raystack.cpp @@ -68,9 +68,9 @@ JonesMatrix Rayleigh_Stack_BRDF_Model::jonesDSC() // in mind if one is interested in coherent scattering between two // defects... complex total_phase(0,0); - for (int i=0; i ee=stack.get_e()[i].epsilon(lambda_eff); - double tt = k*stack.get_t()[i]; + for (int i=0; iget_n(); ++i) { + complex ee=stack->get_e()[i].epsilon(lambda_eff); + double tt = k*stack->get_t()[i]; ee = 1.; total_phase += sqrt(ee-sqr(sin(thetai)))*tt; total_phase += sqrt(ee-sqr(sin(thetas)))*tt; @@ -189,9 +189,9 @@ JonesMatrix Rayleigh_Stack_BRDF_Model::jonesDSC() // The total phase is not mentioned in the text, but must be kept in mind if one // is interested in coherent scattering between two defects... complex total_phase(0,0); - for (int i=0; i ee=stack.get_e()[i].epsilon(lambda_eff); - double tt = k*stack.get_t()[i]; + for (int i=0; iget_n(); ++i) { + complex ee=stack->get_e()[i].epsilon(lambda_eff); + double tt = k*stack->get_t()[i]; total_phase += sqrt(ee-sqr(sin(thetai)))*tt; total_phase -= sqrt(ee-sqr(sinthetasext))*tt; } @@ -294,14 +294,14 @@ void Rayleigh_Stack_BRDF_Model::setup() int layer_number; - double stackthickness = stack.get_total_thickness(); + double stackthickness = stack->get_total_thickness(); if (depth>=0 && depth=0; --ii) { - double t = stack.get_t()[ii]; + for (int ii=stack->get_n()-1; ii>=0; --ii) { + double t = stack->get_t()[ii]; if (_depth<=t&&_depth>=0) { layer_number = ii; if (is_forward()) { @@ -316,37 +316,37 @@ void Rayleigh_Stack_BRDF_Model::setup() // Grow films below the defect... if (is_forward()) { for (int ii=0; iiget_e()[ii].index(lambda),stack->get_t()[ii]); } } else { for (int ii=layer_number-1; ii>=0; --ii) { - above.grow(optical_constant(COMPLEX(stack.get_e()[ii].index(lambda))/bscale),stack.get_t()[ii]); + above.grow(optical_constant(COMPLEX(stack->get_e()[ii].index(lambda))/bscale),stack->get_t()[ii]); } } // Get defect's layer information... int i=layer_number; if (is_forward()) { - material2 = stack.get_e()[i].index(lambda); + material2 = stack->get_e()[i].index(lambda); } else { - material2 = optical_constant(COMPLEX(stack.get_e()[i].index(lambda))/bscale); + material2 = optical_constant(COMPLEX(stack->get_e()[i].index(lambda))/bscale); } - tau = stack.get_t()[i]; + tau = stack->get_t()[i]; // Grow films above the defect... if (is_forward()) { - for (int i=layer_number+1; iget_n(); ++i) { + above.grow(stack->get_e()[i].index(lambda),stack->get_t()[i]); } } else { - for (int i=stack.get_n()-1; i>=layer_number+1; --i) { - below.grow(optical_constant(COMPLEX(stack.get_e()[i].index(lambda))/bscale),stack.get_t()[i]); + for (int i=stack->get_n()-1; i>=layer_number+1; --i) { + below.grow(optical_constant(COMPLEX(stack->get_e()[i].index(lambda))/bscale),stack->get_t()[i]); } } } else if (depth<0 && is_forward()) { // Grow films below the defect... - for (int ii=0; iiget_n(); ++ii) { + below.grow(stack->get_e()[ii].index(lambda),stack->get_t()[ii]); } // Get defect's layer information... @@ -359,8 +359,8 @@ void Rayleigh_Stack_BRDF_Model::setup() material2 = optical_constant(1./bscale); tau = -depth; d = 0; - for (int ii=stack.get_n()-1; ii>=0; --ii) { - above.grow(optical_constant(COMPLEX(stack.get_e()[ii].index(lambda))/bscale),stack.get_t()[ii]); + for (int ii=stack->get_n()-1; ii>=0; --ii) { + above.grow(optical_constant(COMPLEX(stack->get_e()[ii].index(lambda))/bscale),stack->get_t()[ii]); } } else if (depth>stackthickness && is_forward()) { @@ -370,8 +370,8 @@ void Rayleigh_Stack_BRDF_Model::setup() d = 0; // Grow films above the defect... - for (int i=0; iget_n(); ++i) { + above.grow(stack->get_e()[i].index(lambda),stack->get_t()[i]); } } else { @@ -381,8 +381,8 @@ void Rayleigh_Stack_BRDF_Model::setup() d = tau; // Grow films above the defect... - for (int i=stack.get_n()-1; i>=0; --i) { - below.grow(optical_constant(COMPLEX(stack.get_e()[i].index(lambda))/bscale),stack.get_t()[i]); + for (int i=stack->get_n()-1; i>=0; --i) { + below.grow(optical_constant(COMPLEX(stack->get_e()[i].index(lambda))/bscale),stack->get_t()[i]); } } @@ -400,7 +400,7 @@ void Rayleigh_Stack_BRDF_Model::setup() } DEFINE_MODEL(Rayleigh_Stack_BRDF_Model,Local_BRDF_Model,"A Rayleigh defect in a stack of layers"); -DEFINE_PARAMETER(Rayleigh_Stack_BRDF_Model,dielectric_stack,stack,"Dielectric stack","",0xFF); +DEFINE_PTRPARAMETER(Rayleigh_Stack_BRDF_Model,StackModel_Ptr,stack,"Dielectric stack","No_StackModel",0xFF); DEFINE_PARAMETER(Rayleigh_Stack_BRDF_Model,double,depth,"Depth into stack [um]","0",0xFF); DEFINE_PARAMETER(Rayleigh_Stack_BRDF_Model,double,radius,"Sphere radius [um]","0.01",0xFF); DEFINE_PARAMETER(Rayleigh_Stack_BRDF_Model,dielectric_function,sphere,"Sphere optical properties","(1,0)",0xFF); diff --git a/code/raystack.h b/code/raystack.h index e40637e..80755e1 100644 --- a/code/raystack.h +++ b/code/raystack.h @@ -34,7 +34,7 @@ namespace SCATMECH { SCATMECH::JonesMatrix jonesDSC(); DECLARE_MODEL(); - DECLARE_PARAMETER(dielectric_stack,stack); + DECLARE_PARAMETER(StackModel_Ptr,stack); DECLARE_PARAMETER(double,depth); DECLARE_PARAMETER(double,radius); DECLARE_PARAMETER(dielectric_function,sphere); diff --git a/code/rcw.cpp b/code/rcw.cpp index 617a713..bec5ea2 100644 --- a/code/rcw.cpp +++ b/code/rcw.cpp @@ -432,7 +432,7 @@ namespace SCATMECH { 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); + muinvx[i+nmat]= (i==0 ? 1.: 0); muy[i+nmat]= (i==0 ? 1.: 0); muz[i+nmat]= (i==0 ? 1.: 0); } diff --git a/code/rcw.h b/code/rcw.h index f8f77e7..c0b51eb 100644 --- a/code/rcw.h +++ b/code/rcw.h @@ -127,7 +127,7 @@ namespace SCATMECH { std::vector kxi,Kx; std::vector kIzi,kIIzi,YI,YII,ZI,ZII; - double ky; + double ky; }; class RCW_BRDF_Model: public BRDF_Model { diff --git a/code/roughnes.cpp b/code/roughnes.cpp index 695d72c..4637f26 100644 --- a/code/roughnes.cpp +++ b/code/roughnes.cpp @@ -58,7 +58,7 @@ namespace SCATMECH { const int M11=1,M12=2,M21=3,M22=0; - int L=stack.get_n(); + int L=stack->get_n(); int m = is_forward() ? this_layer : L-this_layer; int n; @@ -88,12 +88,12 @@ namespace SCATMECH { if (is_forward()) { eps[0]=(COMPLEX)(substrate.epsilon(lambda)); for (n=1; nget_e()[n-1].epsilon(lambda)); } eps[L+1]=1; d[0]=0; - for (n=1; nget_t()[n-1]/lambda; } else { // is_backward() _lambda = lambda/substrate.n(lambda); scatter_medium = is_reflection() ? 1. : 1/substrate.n(lambda); @@ -101,12 +101,12 @@ namespace SCATMECH { eps[0]=1./(COMPLEX)(substrate.epsilon(lambda)); for (n=1; nget_e()[L-n].epsilon(lambda))/(COMPLEX)(substrate.epsilon(lambda)); } eps[L+1]=1.; d[0]=0; - for (n=1; nget_t()[L-n]/_lambda; } // @@ -261,7 +261,7 @@ namespace SCATMECH { DEFINE_MODEL(Roughness_Stack_BRDF_Model,Roughness_BRDF_Model, "Scattering by a single rough interface in a stack of films."); - DEFINE_PARAMETER(Roughness_Stack_BRDF_Model,dielectric_stack,stack,"Film stack on substrate","",0xFF); + DEFINE_PTRPARAMETER(Roughness_Stack_BRDF_Model,StackModel_Ptr,stack,"Film stack on substrate","No_StackModel",0xFF); DEFINE_PARAMETER(Roughness_Stack_BRDF_Model,int,this_layer,"Rough interface","0",0xFF); diff --git a/code/roughnes.h b/code/roughnes.h index 5b54acd..c7b15bc 100644 --- a/code/roughnes.h +++ b/code/roughnes.h @@ -40,7 +40,7 @@ namespace SCATMECH { DECLARE_PARAMETER(int,this_layer); // The films - DECLARE_PARAMETER(dielectric_stack,stack); + DECLARE_PARAMETER(StackModel_Ptr,stack); protected: diff --git a/code/scateval.cpp b/code/scateval.cpp index 543d487..f865ed4 100644 --- a/code/scateval.cpp +++ b/code/scateval.cpp @@ -23,10 +23,10 @@ namespace SCATMECH { namespace { // List of binary operators... - static const char binops[] = ",&|><=-+*/^"; + static const char binops[] = ",&|><=+-*/^"; // Precedence of each of the binary operators above... - static const int precs[] = {0,1,1,2,2,2,3,3,4,4,6}; + static const int precs[] = {0,1,1,2,2,2,3,3,5,5,6}; // This function returns the position of the character c in the string p... // (used to find precedence of operator) @@ -110,8 +110,13 @@ namespace SCATMECH { char pre = input.peek(); int sign = 1; if (pre == '-') { - input.get(); - sign = -1; + if (val_stack.size()==0) { + val_stack.push(0.); + get_operator(); + } else { + input.get(); + sign = -1; + } } else if (pre == '+') { input.get(); } diff --git a/code/scattabl.cpp b/code/scattabl.cpp index b14e388..37d7fea 100644 --- a/code/scattabl.cpp +++ b/code/scattabl.cpp @@ -694,16 +694,19 @@ namespace SCATMECH { } maps_mutex.lock(); - TableFileMap::iterator t = tablefiles.find(filename2); FormulaFileMap::iterator f = formulafiles.find(filename2); - if (t!=tablefiles.end()) { + TableFileMap::iterator tend = tablefiles.end(); + FormulaFileMap::iterator fend = formulafiles.end(); + maps_mutex.unlock(); + + if (t!=tend) { x = (t->second)[0]; y = (t->second)[col-1]; icol=col; params.clear(); tablefile=&(t->second); - } else if (f!=formulafiles.end()) { + } else if (f!=fend) { x = f->second.get_column(1,params); y = f->second.get_column(col,params); params = f->second.get_parameter_map(); @@ -724,16 +727,20 @@ namespace SCATMECH { if (word.size()==0) throw SCATMECH_exception("File " + filename + " is empty"); } if (word=="PARAMETERS") { + maps_mutex.lock(); formulafiles[filename2]=FormulaFile(filename2); f = formulafiles.find(filename2); + maps_mutex.unlock(); x = f->second.get_column(1,params); y = f->second.get_column(col,params); params = f->second.get_parameter_map(); formulafile = &(f->second); icol=col; } else { - tablefiles[filename2]=TableFile(filename2); + maps_mutex.lock(); + tablefiles[filename2]=TableFile(filename2); t=tablefiles.find(filename2); + maps_mutex.unlock(); x = (t->second)[0]; y = (t->second)[col-1]; params.clear(); @@ -743,7 +750,7 @@ namespace SCATMECH { } } - if (x.size()>0) { + if (x.size()>0) { values.resize(x.size()); for (int i=0; i<(int)x.size(); ++i) { @@ -760,8 +767,6 @@ namespace SCATMECH { str << filename2 << " (col=" << col << ")"; name = str.str(); - maps_mutex.unlock(); - return 0; } @@ -828,6 +833,7 @@ namespace SCATMECH { return; } params[param]=from_string(value); + last.x=numeric_limits::max(); } string Table::get_parameter(const std::string& param) const diff --git a/code/sphdfct.cpp b/code/sphdfct.cpp index e852ee1..f2132e3 100644 --- a/code/sphdfct.cpp +++ b/code/sphdfct.cpp @@ -100,10 +100,10 @@ namespace SCATMECH { // Transmission coefficients... const dielectric_constant vacuum(1.); - 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); + 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)); @@ -167,10 +167,10 @@ namespace SCATMECH { // Reflection and transmission coefficients... - COMPLEX ts = stack.ts12(thetai,lambda,vacuum,substrate)*sqrt(n); - COMPLEX tp = stack.tp12(thetai,lambda,vacuum,substrate)*sqrt(n); - COMPLEX rs = stack.rs21i(fabs(thetas),lambda,substrate,vacuum); - COMPLEX rp = stack.rp21i(fabs(thetas),lambda,substrate,vacuum); + COMPLEX ts = stack->ts12(thetai,lambda,vacuum,substrate)*sqrt(n); + COMPLEX tp = stack->tp12(thetai,lambda,vacuum,substrate)*sqrt(n); + COMPLEX rs = stack->rs21i(fabs(thetas),lambda,substrate,vacuum); + COMPLEX rp = stack->rp21i(fabs(thetas),lambda,substrate,vacuum); COMPLEX rayleigh = (e_part-e)/(e_part+2.*e)*cube(x); @@ -219,10 +219,10 @@ namespace SCATMECH { // Reflection coefficients... - COMPLEX rsi = stack.rs21i(fabs(thetai),lambda,substrate,vacuum); - COMPLEX rpi = stack.rp21i(fabs(thetai),lambda,substrate,vacuum); - COMPLEX rss = stack.rs21i(fabs(thetas),lambda,substrate,vacuum); - COMPLEX rps = stack.rp21i(fabs(thetas),lambda,substrate,vacuum); + COMPLEX rsi = stack->rs21i(fabs(thetai),lambda,substrate,vacuum); + COMPLEX rpi = stack->rp21i(fabs(thetai),lambda,substrate,vacuum); + COMPLEX rss = stack->rs21i(fabs(thetas),lambda,substrate,vacuum); + COMPLEX rps = stack->rp21i(fabs(thetas),lambda,substrate,vacuum); COMPLEX rayleigh = (e_part-e)/(e_part+2.*e)*cube(x); @@ -271,10 +271,10 @@ namespace SCATMECH { // Reflection and transmission coefficients... - COMPLEX rsi = stack.rs21i(fabs(thetai),lambda,substrate,vacuum); - COMPLEX rpi = stack.rp21i(fabs(thetai),lambda,substrate,vacuum); - COMPLEX tss = stack.ts21(thetas,lambda,substrate,vacuum)*sqrt(cos(thetas)/cos_thetas_internal/n); - COMPLEX tps = stack.tp21(thetas,lambda,substrate,vacuum)*sqrt(cos(thetas)/cos_thetas_internal/n); + COMPLEX rsi = stack->rs21i(fabs(thetai),lambda,substrate,vacuum); + COMPLEX rpi = stack->rp21i(fabs(thetai),lambda,substrate,vacuum); + COMPLEX tss = stack->ts21(thetas,lambda,substrate,vacuum)*sqrt(cos(thetas)/cos_thetas_internal/n); + COMPLEX tps = stack->tp21(thetas,lambda,substrate,vacuum)*sqrt(cos(thetas)/cos_thetas_internal/n); COMPLEX rayleigh = (e_part-e)/(e_part+2.*e)*cube(x); @@ -313,7 +313,7 @@ namespace SCATMECH { DEFINE_MODEL(Rayleigh_Defect_BRDF_Model,Local_BRDF_Model, "Scattering by a subsurface defect in the Rayleigh limit."); - DEFINE_PARAMETER(Rayleigh_Defect_BRDF_Model,dielectric_stack,stack,"Film stack on substrate","",0xFF); + DEFINE_PTRPARAMETER(Rayleigh_Defect_BRDF_Model,StackModel_Ptr,stack,"Film stack on substrate","No_StackModel",0xFF); DEFINE_PARAMETER(Rayleigh_Defect_BRDF_Model,double,radius,"Defect radius [um]","0.001",0xFF); diff --git a/code/sphdfct.h b/code/sphdfct.h index a5a63f4..1cf0a8b 100644 --- a/code/sphdfct.h +++ b/code/sphdfct.h @@ -45,7 +45,7 @@ namespace SCATMECH { // The dielectric function of the particle DECLARE_PARAMETER(dielectric_function,defect); - DECLARE_PARAMETER(dielectric_stack,stack); + DECLARE_PARAMETER(StackModel_Ptr,stack); protected: diff --git a/code/sphprt.cpp b/code/sphprt.cpp index e64393f..6b98b23 100644 --- a/code/sphprt.cpp +++ b/code/sphprt.cpp @@ -108,8 +108,8 @@ namespace SCATMECH { Vector a4in = cross(b4in,inKr), a4out = cross(b4out,outKr); // The reflection coefficients... - JonesMatrix ri = stack.r12(thetai,lambda,vacuum,substrate); - JonesMatrix rs = stack.r12(thetas,lambda,vacuum,substrate); + JonesMatrix ri = stack->r12(thetai,lambda,vacuum,substrate); + JonesMatrix rs = stack->r12(thetas,lambda,vacuum,substrate); JonesMatrix S1 = scatterer->jones(_Euler*inK, _Euler*outK); JonesMatrix S2 = scatterer->jones(_Euler*inKr,_Euler*outK); @@ -206,13 +206,13 @@ namespace SCATMECH { // The reflected wave has a reflection and extra path length... COMPLEX phase = exp(COMPLEX(0,2)*cos(thetai)*k*distance); - JonesMatrix r = stack.r12(thetai,lambda,vacuum,substrate); + JonesMatrix r = stack->r12(thetai,lambda,vacuum,substrate); scatter_indirect = phase*matrixout*scatter*matrixin*r; } // Transmission through the interface... - JonesMatrix t = stack.t12(thetas_outside,lambda,vacuum,substrate)* + JonesMatrix t = stack->t12(thetas_outside,lambda,vacuum,substrate)* sqrt(n*cos(thetas)/cos_thetas_outside); // Total scatter... @@ -263,8 +263,8 @@ namespace SCATMECH { Vector outPi = cross(outKi,outSi); // Transmission through the interface... - JonesMatrix ti = stack.t21i(thetai,lambda,substrate,vacuum)/(COMPLEX)sqrt(n); - JonesMatrix ts = stack.t21i(thetas,lambda,substrate,vacuum)*(COMPLEX)sqrt(cos_thetas_at_part/cos(thetas)/n); + JonesMatrix ti = stack->t21i(thetai,lambda,substrate,vacuum)/(COMPLEX)sqrt(n); + JonesMatrix ts = stack->t21i(thetas,lambda,substrate,vacuum)*(COMPLEX)sqrt(cos_thetas_at_part/cos(thetas)/n); // Local polarization basis set, so that {par*,perp,k} is right handed... Vector perp = perpto(inKi,outKi); @@ -352,7 +352,7 @@ namespace SCATMECH { Vector pars = cross(perp,outKr); // Reflection coefficient from surface... - JonesMatrix r = stack.r12(thetas,lambda,vacuum,substrate); + JonesMatrix r = stack->r12(thetas,lambda,vacuum,substrate); // Rotations to/from {s,p,k} and {par*,perp,k}... JonesMatrix matrixin = GetJonesRotator(pari,perp,inS,inP); @@ -368,7 +368,7 @@ namespace SCATMECH { } // The transmission coefficient through the interface... - JonesMatrix ti = stack.t21i(thetai,lambda,substrate,vacuum)/(COMPLEX)sqrt(n); + JonesMatrix ti = stack->t21i(thetai,lambda,substrate,vacuum)/(COMPLEX)sqrt(n); // The total scatter... JonesMatrix scatter = (scatter_direct+scatter_indirect)*ti/k; @@ -398,7 +398,7 @@ namespace SCATMECH { DEFINE_MODEL(Double_Interaction_BRDF_Model,Local_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","",0xFF); + DEFINE_PTRPARAMETER(Double_Interaction_BRDF_Model,StackModel_Ptr,stack,"Film stack on substrate","No_StackModel",0xFF); DEFINE_PARAMETER(Double_Interaction_BRDF_Model,double,distance,"Distance from center to surface [um]","0.05",0xFF); DEFINE_PTRPARAMETER(Double_Interaction_BRDF_Model,Free_Space_Scatterer_Ptr,scatterer,"The scattering function","MieScatterer",0xFF); DEFINE_PARAMETER(Double_Interaction_BRDF_Model,double,alpha,"Rotation of particle about z-axis [rad]","0",0xFF); diff --git a/code/sphprt.h b/code/sphprt.h index 9ed4c7d..f8bf369 100644 --- a/code/sphprt.h +++ b/code/sphprt.h @@ -58,7 +58,7 @@ namespace SCATMECH { DECLARE_PARAMETER(Free_Space_Scatterer_Ptr,scatterer); // Any dielectric layers on the surface... - DECLARE_PARAMETER(dielectric_stack,stack); + DECLARE_PARAMETER(StackModel_Ptr,stack); // The parameters for this virtual model are the Euler angles for the rotation of the // scatterer... diff --git a/code/subbobvlieg.cpp b/code/subbobvlieg.cpp index 1e180ca..0c1db10 100644 --- a/code/subbobvlieg.cpp +++ b/code/subbobvlieg.cpp @@ -53,22 +53,25 @@ namespace SCATMECH { model.set_sphere((optical_constant)((COMPLEX)(sphere.index(lambda))/n)); model.set_radius(radius); dielectric_stack temp; - for (int i=0; iget_n(); ++i) { int j = i; - COMPLEX index = (COMPLEX)(spherecoat.get_e()[j].index(lambda))/n; - double thickness = spherecoat.get_t()[j]; + COMPLEX index = (COMPLEX)(spherecoat->get_e()[j].index(lambda))/n; + double thickness = spherecoat->get_t()[j]; temp.grow(optical_constant(index),thickness); } - model.set_spherecoat(temp); + Stack_StackModel temp0; + temp0.set_stack(temp); + model.set_spherecoat(temp0.clone()); temp.wash(); - for (int i=0; iget_n(); ++i) { // Need to reverse order of films... - int j = stack.get_n()-i-1; - COMPLEX index = (COMPLEX)(stack.get_e()[j].index(lambda))/n; - double thickness = stack.get_t()[j]; + int j = stack->get_n()-i-1; + COMPLEX index = (COMPLEX)(stack->get_e()[j].index(lambda))/n; + double thickness = stack->get_t()[j]; temp.grow(optical_constant(index),thickness); } - model.set_stack(temp); + temp0.set_stack(temp); + model.set_stack(temp0.clone()); model.set_delta(delta); model.set_lmax(lmax); model.set_order(order); @@ -120,14 +123,16 @@ namespace SCATMECH { model.set_particle((optical_constant)((COMPLEX)(particle.index(lambda))/n)); model.set_Shape(Shape); dielectric_stack temp; - for (int i=0; iget_n(); ++i) { // Need to reverse order of films... - int j = stack.get_n()-i-1; - COMPLEX index = (COMPLEX)(stack.get_e()[j].index(lambda))/n; - double thickness = stack.get_t()[j]; + int j = stack->get_n()-i-1; + COMPLEX index = (COMPLEX)(stack->get_e()[j].index(lambda))/n; + double thickness = stack->get_t()[j]; temp.grow(optical_constant(index),thickness); } - model.set_stack(temp); + Stack_StackModel temp0; + temp0.set_stack(temp); + model.set_stack(temp0.clone()); model.set_delta(delta); model.set_lmax(lmax); model.set_mmax(mmax); @@ -152,8 +157,8 @@ namespace SCATMECH { DEFINE_PARAMETER(Subsurface_Bobbert_Vlieger_BRDF_Model,dielectric_function,sphere,"Sphere optical properties","(1.59,0)",0xFF); DEFINE_PARAMETER(Subsurface_Bobbert_Vlieger_BRDF_Model,double,radius,"Particle radius [um]","0.05",0xFF); - DEFINE_PARAMETER(Subsurface_Bobbert_Vlieger_BRDF_Model,dielectric_stack,spherecoat,"Coatings on the sphere","",0xFF); - DEFINE_PARAMETER(Subsurface_Bobbert_Vlieger_BRDF_Model,dielectric_stack,stack,"Substrate films","",0xFF); + DEFINE_PTRPARAMETER(Subsurface_Bobbert_Vlieger_BRDF_Model,StackModel_Ptr,spherecoat,"Coatings on the sphere","No_StackModel",0xFF); + DEFINE_PTRPARAMETER(Subsurface_Bobbert_Vlieger_BRDF_Model,StackModel_Ptr,stack,"Substrate films","No_StackModel",0xFF); DEFINE_PARAMETER(Subsurface_Bobbert_Vlieger_BRDF_Model,double,delta,"Separation of particle from substrate [um] (in contact: 0)","0",0xFF); DEFINE_PARAMETER(Subsurface_Bobbert_Vlieger_BRDF_Model,int,lmax,"Maximum spherical harmonic order (use Bohren & Huffman estimate: 0)","0",0xFF); DEFINE_PARAMETER(Subsurface_Bobbert_Vlieger_BRDF_Model,int,order,"Perturbation order (exact: -1)","-1",0xFF); @@ -165,7 +170,7 @@ namespace SCATMECH { DEFINE_PTRPARAMETER(Subsurface_Axisymmetric_Particle_BRDF_Model,Axisymmetric_Shape_Ptr,Shape,"Particle Shape","Ellipsoid_Axisymmetric_Shape",0xFF); DEFINE_PARAMETER(Subsurface_Axisymmetric_Particle_BRDF_Model,dielectric_function,particle,"Particle optical properties","(1.59,0)",0xFF); - DEFINE_PARAMETER(Subsurface_Axisymmetric_Particle_BRDF_Model,dielectric_stack,stack,"Substrate films","",0xFF); + DEFINE_PTRPARAMETER(Subsurface_Axisymmetric_Particle_BRDF_Model,StackModel_Ptr,stack,"Substrate films","No_StackModel",0xFF); DEFINE_PARAMETER(Subsurface_Axisymmetric_Particle_BRDF_Model,double,delta,"Separation of particle from substrate [um] (in contact: 0)","0",0xFF); DEFINE_PARAMETER(Subsurface_Axisymmetric_Particle_BRDF_Model,int,lmax,"Maximum polar order (lmax)","0",0xFF); DEFINE_PARAMETER(Subsurface_Axisymmetric_Particle_BRDF_Model,int,mmax,"Maximum azimuthal order (mmax)","0",0xFF); diff --git a/code/subbobvlieg.h b/code/subbobvlieg.h index e61ec65..c3e93b1 100644 --- a/code/subbobvlieg.h +++ b/code/subbobvlieg.h @@ -26,8 +26,8 @@ namespace SCATMECH { DECLARE_MODEL(); DECLARE_PARAMETER(dielectric_function,sphere); DECLARE_PARAMETER(double,radius); - DECLARE_PARAMETER(dielectric_stack,spherecoat); - DECLARE_PARAMETER(dielectric_stack,stack); + DECLARE_PARAMETER(StackModel_Ptr,spherecoat); + DECLARE_PARAMETER(StackModel_Ptr,stack); DECLARE_PARAMETER(double,delta); DECLARE_PARAMETER(int,lmax); DECLARE_PARAMETER(int,order); @@ -60,7 +60,7 @@ namespace SCATMECH { DECLARE_MODEL(); DECLARE_PARAMETER(Axisymmetric_Shape_Ptr,Shape); DECLARE_PARAMETER(dielectric_function,particle); - DECLARE_PARAMETER(dielectric_stack,stack); + DECLARE_PARAMETER(StackModel_Ptr,stack); DECLARE_PARAMETER(double,delta); DECLARE_PARAMETER(int,lmax); DECLARE_PARAMETER(int,mmax); diff --git a/code/subsphere.cpp b/code/subsphere.cpp index 1d620e8..c6b2b25 100644 --- a/code/subsphere.cpp +++ b/code/subsphere.cpp @@ -39,7 +39,7 @@ namespace SCATMECH { if (substrate.k(lambda)!=0) error("Substrate cannot be absorbing"); if (scatterer->get_lambda()!=lambda) error("scatterer.lambda!=lambda"); - if (scatterer->get_medium().n(lambda)!=substrate.n(lambda)) error("scatterer.medium!=substrate.medium"); + if (scatterer->get_medium().n(lambda)!=substrate.n(lambda)) error("scatterer.medium!=substrate"); Euler = Matrix(cos(alpha)*cos(beta), sin(alpha), -cos(alpha)*sin(beta), -cos(beta)*sin(alpha), cos(alpha), sin(alpha)*sin(beta), @@ -99,8 +99,8 @@ namespace SCATMECH { Vector outPi = cross(outKi,outSi); // Transmission coefficients into and out of the material... - JonesMatrix ti = stack.t12(thetai,lambda,vacuum,substrate)*(COMPLEX)sqrt(n); - JonesMatrix ts = stack.t12(thetas,lambda,vacuum,substrate)*(COMPLEX)sqrt(n*cos_thetas_internal/cos(thetas)); + JonesMatrix ti = stack->t12(thetai,lambda,vacuum,substrate)*(COMPLEX)sqrt(n); + JonesMatrix ts = stack->t12(thetas,lambda,vacuum,substrate)*(COMPLEX)sqrt(n*cos_thetas_internal/cos(thetas)); // Polarization basis local to the scattering plane... {par*,perp,k} are right handed... Vector perp = perpto(inKi,outKi); @@ -180,7 +180,7 @@ namespace SCATMECH { Vector pars = cross(perp,outKr); // The reflection coefficients... - JonesMatrix r = stack.r21i(fabs(thetas),lambda,substrate,vacuum); + JonesMatrix r = stack->r21i(fabs(thetas),lambda,substrate,vacuum); // Rotation matrices... JonesMatrix matrixin = GetJonesRotator(pari,perp,inSi,inPi); @@ -197,7 +197,7 @@ namespace SCATMECH { } // Transmission into substrate... - JonesMatrix ti = stack.t12(thetai,lambda,vacuum,substrate)*(COMPLEX)sqrt(n); + JonesMatrix ti = stack->t12(thetai,lambda,vacuum,substrate)*(COMPLEX)sqrt(n); JonesMatrix scatter = (scatter_direct+scatter_indirect)*ti/k/n; @@ -248,8 +248,8 @@ namespace SCATMECH { Vector a4in = cross(b4in,inKr), a4out = cross(b4out,outKr); // The reflection coefficients... - JonesMatrix ri = stack.r21i(fabs(thetai),lambda,substrate,vacuum); - JonesMatrix rs = stack.r21i(fabs(thetas),lambda,substrate,vacuum); + JonesMatrix ri = stack->r21i(fabs(thetai),lambda,substrate,vacuum); + JonesMatrix rs = stack->r21i(fabs(thetas),lambda,substrate,vacuum); // The scattering matrices... JonesMatrix S1 = scatterer->jones(_Euler*inK,_Euler*outK); @@ -349,14 +349,14 @@ namespace SCATMECH { // Phase and reflection coefficients... COMPLEX phase = exp(COMPLEX(0,2)*cos(thetai)*k*n*depth); - JonesMatrix r = stack.r21i(fabs(thetai),lambda,substrate,vacuum); + JonesMatrix r = stack->r21i(fabs(thetai),lambda,substrate,vacuum); // Scattering matrix in global basis... scatter_indirect = phase*matrixout*scatter*matrixin*r; } // Transmission out of material... - JonesMatrix t = stack.t21(thetas,lambda,substrate,vacuum)* + JonesMatrix t = stack->t21(thetas,lambda,substrate,vacuum)* sqrt(cos(thetas)/cos(thetas_inside)/n); JonesMatrix scatter = t*(scatter_direct+scatter_indirect)/k/n; @@ -373,7 +373,7 @@ namespace SCATMECH { DEFINE_MODEL(Subsurface_Particle_BRDF_Model,Local_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","",0xFF); + DEFINE_PTRPARAMETER(Subsurface_Particle_BRDF_Model,StackModel_Ptr,stack,"Film stack on substrate","No_StackModel",0xFF); DEFINE_PARAMETER(Subsurface_Particle_BRDF_Model,double,depth,"Depth of center of particle [um]","0",0xFF); DEFINE_PTRPARAMETER(Subsurface_Particle_BRDF_Model,Free_Space_Scatterer_Ptr,scatterer,"The scattering function","MieScatterer",0xFF); DEFINE_PARAMETER(Subsurface_Particle_BRDF_Model,double,alpha,"Rotation of particle about z-axis [rad]","0",0xFF); diff --git a/code/subsphere.h b/code/subsphere.h index 1fa045e..4df3f8e 100644 --- a/code/subsphere.h +++ b/code/subsphere.h @@ -29,7 +29,7 @@ namespace SCATMECH { DECLARE_MODEL(); DECLARE_PARAMETER(double,depth); DECLARE_PARAMETER(Free_Space_Scatterer_Ptr,scatterer); - DECLARE_PARAMETER(dielectric_stack,stack); + DECLARE_PARAMETER(StackModel_Ptr,stack); // Particle can be rotated from its natural orientation... // Rotation of particle about z-axis... diff --git a/code/transmit.cpp b/code/transmit.cpp index 6049e9a..e3e187f 100644 --- a/code/transmit.cpp +++ b/code/transmit.cpp @@ -41,7 +41,7 @@ namespace SCATMECH { { double thetas_internal = asin(sin(thetas)/substrate.n(lambda)); - MuellerMatrix T = films.t12(thetas,lambda,vacuum,substrate); + MuellerMatrix T = films->t12(thetas,lambda,vacuum,substrate); double transmittancefactorout = substrate.n(lambda)*cos(thetas_internal)/cos(thetas); double jacobian = cos(thetas)/cos(thetas_internal)/sqr(substrate.n(lambda)); double ratiocosines = cos(thetas_internal)/cos(thetas); @@ -55,8 +55,8 @@ namespace SCATMECH { { double thetas_internal = asin(sin(thetas)/substrate.n(lambda)); double thetai_internal = asin(sin(thetai)/substrate.n(lambda)); - MuellerMatrix Ti = films.t12(thetai,lambda,vacuum,substrate); - MuellerMatrix Ts = films.t12(thetas,lambda,vacuum,substrate); + MuellerMatrix Ti = films->t12(thetai,lambda,vacuum,substrate); + MuellerMatrix Ts = films->t12(thetas,lambda,vacuum,substrate); double transmittancefactorin = cos(thetai_internal)/cos(thetai)*substrate.n(lambda); double transmittancefactorout = cos(thetas_internal)/cos(thetas)*substrate.n(lambda); double jacobian = cos(thetas)/cos(thetas_internal)/sqr(substrate.n(lambda)); @@ -70,7 +70,7 @@ namespace SCATMECH { case 3: { double thetai_internal = asin(sin(thetai)/substrate.n(lambda)); - MuellerMatrix Ti = films.t12(thetai,lambda,vacuum,substrate); + MuellerMatrix Ti = films->t12(thetai,lambda,vacuum,substrate); double transmittancefactorin = cos(thetai_internal)/cos(thetai)*substrate.n(lambda); double factor = transmittancefactorin; @@ -89,6 +89,6 @@ namespace SCATMECH { "Model which evaluates another BRDF_Model in transmission, outside of a material."); DEFINE_PTRPARAMETER(Transmit_BRDF_Model,BRDF_Model_Ptr,model,"The model","Microroughness_BRDF_Model",0xFF); - DEFINE_PARAMETER(Transmit_BRDF_Model,dielectric_stack,films,"Films on the bottom side","",0xFF); + DEFINE_PTRPARAMETER(Transmit_BRDF_Model,StackModel_Ptr,films,"Films on the bottom side","No_StackModel",0xFF); } // namespace SCATMECH diff --git a/code/transmit.h b/code/transmit.h index 1b51333..91494fa 100644 --- a/code/transmit.h +++ b/code/transmit.h @@ -27,7 +27,7 @@ namespace SCATMECH { DECLARE_MODEL(); DECLARE_PARAMETER(BRDF_Model_Ptr,model); - DECLARE_PARAMETER(dielectric_stack,films); + DECLARE_PARAMETER(StackModel_Ptr,films); protected: virtual MuellerMatrix mueller(); diff --git a/code/twoface.cpp b/code/twoface.cpp index 6628155..a50f5a5 100644 --- a/code/twoface.cpp +++ b/code/twoface.cpp @@ -30,9 +30,10 @@ namespace SCATMECH { model.set_type(type); model.set_substrate(substrate); model.set_psd(psd); - dielectric_stack stack; - stack.grow(film,thickness); - model.set_stack(stack); + SingleFilm_StackModel stack0; + stack0.set_material(film); + stack0.set_thickness(thickness); + model.set_stack(stack0.clone()); model.set_this_layer(face-1); } diff --git a/code/zernikeexpansion.cpp b/code/zernikeexpansion.cpp new file mode 100644 index 0000000..068f393 --- /dev/null +++ b/code/zernikeexpansion.cpp @@ -0,0 +1,180 @@ +//****************************************************************************** +//** SCATMECH: Polarized Light Scattering C++ Class Library +//** +//** File: zernikeexpansion.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 "zernikeexpansion.h" +#include "bobvlieg.h" + +using namespace std; +using namespace SCATMECH; + +namespace SCATMECH { + +/// +/// The following calculates the radial Zernike polynomial: +/// where n is the radial order, l is the azimuthal order, and rho is the argument +/// +static double Radial_Zernike_Polynomial(int n,int l,double rho) +{ + using namespace BobVlieg_Supp; // This defines mpow and Fact (mpow(s) is (-1)^s), Fact(s) is factorial of s. + + int m = abs(l); + if ((n-m)%2==1) return 0.; + + double sum = 0; + for (int s=0;s<=(n-m)/2;++s) { + sum += mpow(s)*Fact(n-s)/Fact(s)/Fact((n+m)/2-s)/Fact((n-m)/2-s)*pow(rho,(double)(n-2*s)); + } + return sum; +}; + +// The azimuthal function +static double az(int k,double phi) +{ + if (k>0) return cos(k*phi); + if (k<0) return -sin(k*phi); + else return 1./sqrt(2.0); +} + +void ZernikeExpansion_BRDF_Model::setup() +{ + // Standard procedure is for setup() to call the parent class setup() + BRDF_Model::setup(); + + // The model uses the xyxy basis, rather than the default psps + model_cs = xyxy; + + // Open the coefficient file + ifstream_with_comments cfile(coefficientfile.c_str()); + + // Throw exception if the coefficient file doesn't open + if (!cfile) error("Cannot open coefficient file."); + + // Clear the list of coefficients + coeff.clear(); + + // Read the coefficients until the end of the file + while (!cfile.eof()) { + short i,j,n,m,k,l,p; + double c; + + // Read a line + string line = cfile.getstrline(); + + // If the line isn't blank... + if (line.size()!=0) { + + // Create a stream from the line + istringstream sline(line.c_str()); + + // The following allows for comma or whitespace delimited data... + sline >> i; + while (isspace(sline.peek())||sline.peek()==',') sline.get(); + sline >> j; + while (isspace(sline.peek())||sline.peek()==',') sline.get(); + sline >> n; + while (isspace(sline.peek())||sline.peek()==',') sline.get(); + sline >> m; + while (isspace(sline.peek())||sline.peek()==',') sline.get(); + sline >> k; + while (isspace(sline.peek())||sline.peek()==',') sline.get(); + sline >> l; + while (isspace(sline.peek())||sline.peek()==',') sline.get(); + sline >> p; + while (isspace(sline.peek())||sline.peek()==',') sline.get(); + sline >> c; + + // If there was a failure, then throw an exception + if (sline.fail()) error("Error reading a line in coefficient file"); + + // Push the coefficient onto the list of coefficients + coeff.push_back(Coefficient(i,j,n,m,k,l,p,c)); + } + } + return; +} + +MuellerMatrix ZernikeExpansion_BRDF_Model::mueller() +{ + // It is standard procedure for all Models to call SETUP() at beginning of + // routines that use + SETUP(); + + // This code does not accept anything for transmission mode or for + // radiation upwelling in the material + throw_transmission(); + throw_backward(); + + // Start the sum with zero + MuellerMatrix result = MuellerZero(); + double brdf = 0; + + // Incident and scattering rho's and phi's + double rhoi = sqrt(2.)*sin(thetai/2.); + double rhor = sqrt(2.)*sin(thetas/2.); + double phii = pi-rotation; + double phir = phis-rotation; + + // Iterate through all the coefficients... + for (list::iterator q=coeff.begin();q!=coeff.end();++q) { + short i = q->i; + short j = q->j; + short n = q->n; + short m = q->m; + short k = q->k; + short l = q->l; + short p = q->p; + double coeff = q->coeff; + + // The following just evaluates lambda^p + double power = p==0?1.:p==1?lambda:p==2?sqr(lambda):p==3?sqr(lambda):pow(lambda,p); + + if (i==1&&j==1) { + // The basis set for the unpolarized BRDF... + double prefact = 0.5/pi*sqrt((n+1.)*(m+1.)/((n==0||(n==m&&l==0))?4.:((n==m||l==0)?2.:1.)))*power; + brdf += coeff * prefact * (Radial_Zernike_Polynomial(n,l,rhoi)*Radial_Zernike_Polynomial(m,l,rhor) + + Radial_Zernike_Polynomial(m,l,rhoi)*Radial_Zernike_Polynomial(n,l,rhor)) * cos(l*(phir-phii)); + } else if (i==j) { + // The basis set for the normalized diagonal elements... + double prefact = 1./sqrt(2.+2.*(n==m?1:0)*(k==l?1:0))*sqrt((n+1.)*(m+1.))/pi*power; + result[i-1][j-1] += coeff * prefact * (Radial_Zernike_Polynomial(n,k,rhoi)*Radial_Zernike_Polynomial(m,l,rhor)*az(k,phii)*az(l,phir) + + Radial_Zernike_Polynomial(m,l,rhoi)*Radial_Zernike_Polynomial(n,k,rhor)*az(k,phir)*az(l,phii)); + } else if (ij in the coefficient file"); + } + + // The normalized Mueller matrix has 1 here... + result[0][0] = 1.; + // The following three elements have sign changes upon transpose... + result[2][0] = -result[2][0]; + result[2][1] = -result[2][1]; + result[3][2] = -result[3][2]; + + // Change to un-normalized Mueller matrix... + result *= brdf*scale; + + // And return! + return result; +} + +// +// These are required for Models in SCATMECH... +// +DEFINE_MODEL(ZernikeExpansion_BRDF_Model,BRDF_Model,"Model using the model of Koenderink and van Doorn, extended by Germer to include the Mueller matrix"); +DEFINE_PARAMETER(ZernikeExpansion_BRDF_Model,string,coefficientfile,"File containing coefficients","",0xFF); +DEFINE_PARAMETER(ZernikeExpansion_BRDF_Model,double,scale,"Scale factor","1",0xFF); + + +} // namespace SCATMECH; diff --git a/code/zernikeexpansion.h b/code/zernikeexpansion.h new file mode 100644 index 0000000..0cbac6d --- /dev/null +++ b/code/zernikeexpansion.h @@ -0,0 +1,57 @@ +//****************************************************************************** +//** SCATMECH: Polarized Light Scattering C++ Class Library +//** +//** File: zernikeexpansion.h +//** +//** 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 +//** +//****************************************************************************** +#ifndef ZERNIKEEXPANSION_H +#define ZERNIKEEXPANSION_H + +#include "brdf.h" +#include + +namespace SCATMECH { + + /// + /// The class ZernikeExpansion_BRDF_Model implements the theory of + /// J.J. Koenderink and A.J. van Doorn, "Phenomenological description of bidirectional surface reflection," J. Opt. Soc. Am. A vol. 15, pp. 2903-2912 (1998) + /// as extended by + /// T.A. Germer, "Full four-dimensional and reciprocal Mueller matrix bidirectional reflectance distribution function of sintered polytetrafluoroethylene," + /// Appl. Opt., to be submitted (2017) + /// + class ZernikeExpansion_BRDF_Model : public BRDF_Model { + public: + + // List of the parameters ... + DECLARE_MODEL(); + DECLARE_PARAMETER(std::string,coefficientfile); /// File containing indices and coefficients + DECLARE_PARAMETER(double,scale); /// An overall scale factor (if the coefficients correspond to unit reflectance) + + // Function evaluating Mueller matrix + MuellerMatrix mueller(); + + protected: + + // Function performing one-time code, namely reading the coefficient file + void setup(); + + // Internal structure holding indices and coefficient + struct Coefficient { + Coefficient(short _i,short _j,short _n,short _m,short _k,short _l,short _p, double _coeff) : i(_i),j(_j),n(_n),m(_m),k(_k),l(_l),p(_p),coeff(_coeff) {} + short i,j,n,m,k,l,p; + double coeff; + }; + + // The list of coefficients + std::list coeff; + }; + +} + +#endif diff --git a/docs/axipart.htm b/docs/axipart.htm index 9165d3f..fa19022 100644 --- a/docs/axipart.htm +++ b/docs/axipart.htm @@ -149,15 +149,14 @@

Parameters:

stack - dielectric_stack + StackModel_Ptr - Description of any stack of coatings on the particle, usually specified by a file. - See dielectric_stack for more information. + Description of any stack of coatings on the particle. - none + No_StackModel diff --git a/docs/bobvlieg.htm b/docs/bobvlieg.htm index b2afd37..83da6a2 100644 --- a/docs/bobvlieg.htm +++ b/docs/bobvlieg.htm @@ -147,29 +147,24 @@

Parameters:

spherecoat - dielectric_stack + StackModel_Ptr Description of any stack of - coatings on the particle, usually specified by a file.
- See dielectric_stack, for - more information.
+ coatings on the particle. + - none + No_StackModel stack - dielectric_stack + StackModel_Ptr Description of any stack - of films deposited on the substrate, usually specified - by a file.
- See dielectric_stack, for - more information.
+ of films deposited on the substrate. - no - films + No_StackModel diff --git a/docs/classes.htm b/docs/classes.htm index ac25e6c..9258499 100644 --- a/docs/classes.htm +++ b/docs/classes.htm @@ -82,6 +82,8 @@

Optical Properties of Materials

  • class dielectric_stack
    • +
  • class StackModel
    • +
  • Reflection and transmission coefficients
    • diff --git a/docs/correlated_stack.htm b/docs/correlated_stack.htm index b32f4fc..ca90637 100644 --- a/docs/correlated_stack.htm +++ b/docs/correlated_stack.htm @@ -118,17 +118,13 @@

    Parameters:

    stack - dielectric_stack + StackModel_Ptr Description of the stack - of films deposited on the substrate, usually specified - by a file.
    - See dielectric_stack, for - more information. + of films deposited on the substrate. - no - films + No_StackModel @@ -141,8 +137,7 @@

    See also:

    "BRDF_Model.htm">BRDF_Model,   Roughness_BRDF_Model,   Roughness_Stack_BRDF_Model,   - dielectric_stack

    + StackModel

    J. M. Elson, "Multilayer-coated optics: guided-wave coupling and scattering by means of interface random diff --git a/docs/diffuse.htm b/docs/diffuse.htm index a452e6c..e0232b8 100644 --- a/docs/diffuse.htm +++ b/docs/diffuse.htm @@ -120,15 +120,12 @@

    Parameters:

    stack - dielectric_stack + StackModel_Ptr Description of any stack - of films deposited on the substrate, usually specified - by a file.
    - See dielectric_stack, for - more information. + of films deposited on the substrate. - no films + No_StackModel diff --git a/docs/download.htm b/docs/download.htm index e88b48a..d2067ee 100644 --- a/docs/download.htm +++ b/docs/download.htm @@ -22,7 +22,11 @@

    Download

    Download:

    -

    To download the library, go to GitHub SCATMECH repository and choose Clone or Download. +

    Download release Version 7.00 as + SCATMECH-7.00.zip or + SCATMECH-7.00.tar.gz

    + +

    To download the latest pre-release version of the library, go to GitHub SCATMECH repository and choose Clone or Download.

    diff --git a/docs/facet.htm b/docs/facet.htm index 6bd4aef..7272da6 100644 --- a/docs/facet.htm +++ b/docs/facet.htm @@ -106,20 +106,14 @@

    Parameters:

    stack - dielectric_stack + StackModel_Ptr Description of any stack - of films deposited on the substrate, usually specified - by a file. The films are assumed to be conformal and of + of films deposited on the substrate. The films are assumed to be conformal and of total thickness much less than the horizontal scale of - the roughness.
    - See dielectric_stack, for - more information.
    - (Inherited from Facet_BRDF_Model). + the roughness. - no - films + No_StackModel diff --git a/docs/firstdiffuse.htm b/docs/firstdiffuse.htm index 4395b62..9528840 100644 --- a/docs/firstdiffuse.htm +++ b/docs/firstdiffuse.htm @@ -136,16 +136,12 @@

    Parameters:

    stack - dielectric_stack + StackModel_Ptr Description of any stack - of films deposited on the substrate, usually specified - by a file. See dielectric_stack, for more - information.
    + of films deposited on the substrate. - no - films + No_StackModel diff --git a/docs/flake.htm b/docs/flake.htm index e47f379..ce05065 100644 --- a/docs/flake.htm +++ b/docs/flake.htm @@ -120,20 +120,16 @@

    Parameters:

    stack - dielectric_stack + StackModel_Ptr Description of any stack - of films deposited on the substrate, usually specified - by a file. The films are assumed to be conformal and of + of films deposited on the substrate. The films are assumed to be conformal and of total thickness much less than the horizontal scale of the roughness.
    - See dielectric_stack, for - more information.
    (Inherited from Facet_BRDF_Model). - no - films + No_StackModel @@ -153,17 +149,14 @@

    Parameters:

    overcoat_stack - dielectric_stack + StackModel_Ptr Description of any stack - of films deposited on top of the coating, usually - specified by a file. The films are assumed to be - conformal and to not be rough.
    - See dielectric_stack, for - more information.
    - - no - films + of films deposited on top of the coating. The films are assumed to be + conformal and to not be rough. + + + No_StackModel diff --git a/docs/grating.htm b/docs/grating.htm index 6f66d38..5189b0d 100644 --- a/docs/grating.htm +++ b/docs/grating.htm @@ -819,71 +819,71 @@

    Instantiable Models and Their Parameters:

    stackepsx - dielectric_stack + StackModel_Ptr The description of any film stacks on the substrate. This particular stack corresponds to the dielectric function in the x direction (the grating direction). - none + No_StackModel stackepsy - dielectric_stack + StackModel_Ptr The description of any film stacks on the substrate. This particular stack corresponds to the dielectric function in the y direction (the grating direction). If this parameter is set, it must have the same number of layers and the same layer thicknesses as stackepsx and stackepsz must be set, too. - If this parameter is not set, then the materials are taken to be isotropic and the dielectric functions taken from stackepsx. - none + If this parameter is not set (No_StackModel), then the materials are taken to be isotropic and the dielectric functions taken from stackepsx. + No_StackModel stackepsz - dielectric_stack + StackModel_Ptr The description of any film stacks on the substrate. This particular stack corresponds to the dielectric function in the z direction (the grating direction). If this parameter is set, it must have the same number of layers and the same layer thicknesses as stackepsx and stackepsy must be set, too. - If this parameter is not set, then the materials are taken to be isotropic and the dielectric functions taken from stackepsx. - none + If this parameter is not set (No_StackModel), then the materials are taken to be isotropic and the dielectric functions taken from stackepsx. + No_StackModel stackmux - dielectric_stack + StackModel_Ptr The description of any film stacks on the substrate. This particular stack corresponds to the magnetic function in the x direction (the grating direction). If this parameter is set, it must have the same number of layers and the same layer thicknesses as stackepsx. - If this parameter is not set, any layers are taken to be nonmagnetic. - none + If this parameter is not set (No_StackModel), any layers are taken to be nonmagnetic. + No_StackModel stackmuy - dielectric_stack + StackModel_Ptr The description of any film stacks on the substrate. This particular stack corresponds to the magnetic function in the y direction (the grating direction). If this parameter is set, it must have the same number of layers and the same layer thicknesses as stackepsx, and stackmuz must be set, too. - If this parameter is not set, any layers are taken to be isotropic. - none + If this parameter is not set (No_StackModel), any layers are taken to be isotropic. + No_StackModel stackmuz - dielectric_stack + StackModel_Ptr The description of any film stacks on the substrate. This particular stack corresponds to the magnetic function in the z direction (the grating direction). If this parameter is set, it must have the same number of layers and the same layer thicknesses as stackepsx, and stackmuy must be set, too. - If this parameter is not set, any layers are taken to be isotropic. - none + If this parameter is not set (No_StackModel), any layers are taken to be isotropic. + No_StackModel diff --git a/docs/growth_stack.htm b/docs/growth_stack.htm index 9078428..623a0ef 100644 --- a/docs/growth_stack.htm +++ b/docs/growth_stack.htm @@ -128,16 +128,13 @@

    Parameters:

    stack - dielectric_stack + StackModel_Ptr Description of the stack - of films deposited on the substrate, usually specified - by a file.
    - See dielectric_stack, for - more information. + of films deposited on the substrate. + - no - films + No_StackModel diff --git a/docs/history.htm b/docs/history.htm index c3695fa..d25d03e 100644 --- a/docs/history.htm +++ b/docs/history.htm @@ -18,6 +18,8 @@

    Version History

    + +

    Version 7.01 (September 2017)

    + + +

    Abstract class StackModel was created to give + dielectric_stack the properties of a Model. Specific classes + No_StackModel, Stack_StackModel, SingleFilm_StackModel, + DoubleFilm_StackModel, and GradedFilm_StackModel implement specific film stack models. These + models inherit Model, so they make it easier to + parametrically vary parameters associated with a film + stack. All scattering models that allowed for a dielectric_stack have been upgraded to use StackModel instead.

    + +

    Added Legendre_Phase_Function to the list of classes inheriting Phase_Function.

    + +

    Corrected some errors in the thread safety using mutexes.

    + +

    Corrected some errors in coordinate systems for polarization in BRDF_Model.

    Version 7.00 (January 2015)

    diff --git a/docs/multilayermie.htm b/docs/multilayermie.htm index 1718559..1712381 100644 --- a/docs/multilayermie.htm +++ b/docs/multilayermie.htm @@ -96,7 +96,7 @@

    Parameters:

    stack - dielectric_stack + StackModel_Ptr The dielectric stack describing the thicknesses and optical properties of the coatings on the sphere (everything but the core). The sum @@ -104,7 +104,7 @@

    Parameters:

    corresponds to that closest to the core. - no films + No_StackModel diff --git a/docs/phase.htm b/docs/phase.htm index fbb0b84..016288e 100644 --- a/docs/phase.htm +++ b/docs/phase.htm @@ -140,6 +140,56 @@

    Instantiable Models and Their Parameters:

    has no adjustable parameters. + + + Legendre_Phase_Function:
    + The Legendre phase function is parameterized by a sum of up to six Legendre + polynomials (orders l=0 to l=5). + + + + + c0 + double + The 0-coefficient of the Legendre polynomial expansion. + 1 + + + + c1 + double + The 1-coefficient of the Legendre polynomial expansion. + 0.65 + + + + c2 + double + The 2-coefficient of the Legendre polynomial expansion. + 0.42 + + + + c3 + double + The 3-coefficient of the Legendre polynomial expansion. + 0 + + + + c4 + double + The 4-coefficient of the Legendre polynomial expansion. + 0 + + + + c5 + double + The 5-coefficient of the Legendre polynomial expansion. + 0 + + Table_Phase_Function:
    The class Table_Phase_Function allows the user to diff --git a/docs/raydef.htm b/docs/raydef.htm index ce8e130..e6a8db2 100644 --- a/docs/raydef.htm +++ b/docs/raydef.htm @@ -112,16 +112,12 @@

    Parameters:

    stack - dielectric_stack + StackModel_Ptr Description of any stack - of films deposited on the substrate, usually specified - by a file. See dielectric_stack, for more - information.
    + of films deposited on the substrate. - no - films + No_StackModel diff --git a/docs/rayleighinstrument.htm b/docs/rayleighinstrument.htm index b329f2a..c7fb0ca 100644 --- a/docs/rayleighinstrument.htm +++ b/docs/rayleighinstrument.htm @@ -102,16 +102,12 @@

    Parameters:

    stack - dielectric_stack + StackModel_Ptr Description of any stack - of films deposited on the substrate, usually specified - by a file. See dielectric_stack, for more - information.
    + of films deposited on the substrate. - no - films + No_StackModel diff --git a/docs/raystack.htm b/docs/raystack.htm index e22c730..dfd90f2 100644 --- a/docs/raystack.htm +++ b/docs/raystack.htm @@ -86,11 +86,9 @@

    Parameters:

    stack - dielectric_stack - Description of any stack of films deposited on the substrate, usually specified - by a file. See dielectric_stack, for more - information.
    - no films + StackModel_Ptr + Description of any stack of films deposited on the substrate. + No_StackModel diff --git a/docs/roughnes.htm b/docs/roughnes.htm index be7e50c..fc914e2 100644 --- a/docs/roughnes.htm +++ b/docs/roughnes.htm @@ -113,16 +113,12 @@

    Parameters:

    stack - dielectric_stack + StackModel_Ptr Description of the stack - of films deposited on the substrate, usually specified - by a file.
    - See dielectric_stack, for - more information. + of films deposited on the substrate. - no - films + No_StackModel diff --git a/docs/sphrpart.htm b/docs/sphrpart.htm index a2ff658..3ae97f3 100644 --- a/docs/sphrpart.htm +++ b/docs/sphrpart.htm @@ -114,16 +114,12 @@

    Parameters:

    stack - dielectric_stack + StackModel_Ptr Description of any stack - of films deposited on the substrate, usually specified - by a file. See dielectric_stack, for more - information.
    + of films deposited on the substrate. - no - films + No_StackModel diff --git a/docs/stackmodel.htm b/docs/stackmodel.htm new file mode 100644 index 0000000..dfc4255 --- /dev/null +++ b/docs/stackmodel.htm @@ -0,0 +1,231 @@ + + + + + + SCATMECH: StackModel + + + + + + + + + + +
    +

    SCATMECH > Classes and Functions > + Optical Properties of Materials + > StackModel

    + +

    Abstract class StackModel

    + +
    + +

    The abstract class StackModel is used to give + dielectric_stack the properties of + a Model. + The provided StackModelss are: + +

    Instantiable Models and Their Parameters:

    + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
    ParameterData + TypeDescriptionDefault
    No_StackModel:
    + This model represents an interface without any films. +
    No_StackModel has no adjustable parameters.
    Stack_StackModel:
    + This model is a model that uses a dielectric_stack. +
    stackdielectric_stackThe description of the stack.none
    SingleFilm_StackModel:
    + This model represents a single film. +
    materialdielectric_functionThe optical properties (n and k) of the film.(1.5,0)
    thicknessdoubleThe thickness of the film [µm].0.1
    DoubleFilm_StackModel:
    + This model represents two films. +
    material1dielectric_functionThe optical properties (n and k) of the film closer to the substrate.(1.5,0)
    thickness1doubleThe thickness of the film closer to the substrate [µm].0.1
    material2dielectric_functionThe optical properties (n and k) of the film farther from the substrate.(1.5,0)
    thickness2doubleThe thickness of the film [µm] farther from the substrate.0.1
    GradedFilm_StackModel:
    + A film with a graded index. +
    startdielectric_functionThe optical properties (n and k) of the film closer to the substrate.(1.5,0)
    enddielectric_functionThe optical properties (n and k) of the film farther from the substrate.(1,0)
    thicknessdoubleThe thickness of the film [µm].0.1
    stepsintThe number of steps with which to approximate the film.5
    + +

    See also:

    + +

    SCATMECH Home,   + dielectric_stack +  

    + +

    Include file:

    + + 
    +#include "filmtrans.h"
+
    +

    Source code:

    + + 
    +filmtran.cpp
+
    + +

    Top of Page

    + +
    +

    For More Information

    + +

    +SCATMECH Technical Information and Questions
    +Sensor Science Division Home Page
    +Sensor Science Division Inquiries
    +Website Comments + +

    +Current SCATMECH version: 7.00 (January 2015)
    +This page first online: Version 4.00 (July 2003)
    +This page last modified: Version 7.00 (January 2015)
    + +



















    +



















    + +

    + + + diff --git a/docs/stackmodel.htm~ b/docs/stackmodel.htm~ new file mode 100644 index 0000000..bf0f2f5 --- /dev/null +++ b/docs/stackmodel.htm~ @@ -0,0 +1,223 @@ + + + + + + SCATMECH: StackModel + + + + + + + + + + +
    +

    SCATMECH > Classes and Functions > + Optical Properties of Materials + > StackModel

    + +

    Abstract class StackModel

    + +
    + +

    The abstract class StackModel is used to give + dielectric_stack the properties of + a Model. + The provided StackModelss are: + +

    Instantiable Models and Their Parameters:

    + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
    ParameterData + TypeDescriptionDefault
    No_StackModel:
    + This model represents an interface without any films. +
    Stack_StackModel:
    + This model is a model that uses a dielectric_stack. +
    stackdielectric_stackThe description of the stack.none
    SingleFilm_StackModel:
    + The isotropic phase function corresponds to scattering + which is isotropic. It has no parameters. +
    Isotropic_Phase_Function + has no adjustable parameters.
    DoubleFilm_StackModel:
    + The Rayleigh phase function corresponds to scattering + which follows that expected for unpolarized Rayleigh + scattering. It has no parameters. +
    Rayleigh_Phase_Function + has no adjustable parameters.
    GradedFilm_StackModel:
    + The class Table_Phase_Function allows the user to + provide a interpolation table of values for the phase + function. The ordinate is the angle in degrees, while the + absissa is the phase function, which is assumed to be + properly normalized. The parameter table is the + table of values (see Table). +
    tableTableA look-up + table of phase function versus angle in + degrees.1
    + +

    See also:

    + +

    SCATMECH Home,   + First_Diffuse_BRDF_Model +  

    + +

    H.C. van de Hulst, Multiple Light Scattering: Tables, + Formulas, and Applications, (Academic Press, New York, + 1980), Volume 2, Chap. 10.
    + L.C. Henyey and J.L. Greenstein, "Diffuse radiation in the galaxy," Astrophys. J. 93, 70-83 (1941). +

    + +

    Include file:

    + + 
    +#include "firstdiffuse.h"
+
    +

    Source code:

    + + 
    +firstdiffuse.cpp
+
    +

    Definition of public elements:

    + + 
    +class Phase_Function : public Model {
+public:
+    Phase_Function(int ask=0);
+    virtual double f(double theta)=0;
+    virtual Model* clone() const =0;
+    static Inheritance inheritance;
+    virtual void set_parameter(const std::string& name,const std::string& value);
+};
+
+typedef Model_Ptr<Phase_Function> Phase_Function_Ptr;
+
    + +
    + +

    double f(double + theta)

    + +

    + This + function evaluates the phase function for angle + theta. The angle theta is measured in + radians from the forward scattering direction. It should + be normalized so that the integral over the hemisphere is + one. +

    + +

    Top of Page

    + +

    typedef Model_Ptr<Phase_Function> + Phase_Function_Ptr

    + +

    + The typedef Phase_Function_Ptr behaves like a + pointer to an instance of class Phase_Function. + The following statement will query the user for an + instance of class Phase_Function: + + 

    +Phase_Function_Ptr model = Get_Model_Ptr();
+
    + +

    The + next statement will also create an instance of class + Phase_Function: + + 

    +Phase_Function_Ptr model = "Henyey_Greenstein_Phase_Function";
+
    +

    See + Model_Ptr<model>. +

    + +

    Top of Page

    + +
    +

    For More Information

    + +

    +SCATMECH Technical Information and Questions
    +Sensor Science Division Home Page
    +Sensor Science Division Inquiries
    +Website Comments + +

    +Current SCATMECH version: 7.00 (January 2015)
    +This page first online: Version 4.00 (July 2003)
    +This page last modified: Version 7.00 (January 2015)
    + +



















    +



















    + +

    + + + diff --git a/docs/subaxipart.htm b/docs/subaxipart.htm index feb7f49..f387b16 100644 --- a/docs/subaxipart.htm +++ b/docs/subaxipart.htm @@ -147,15 +147,14 @@

    Parameters:

    stack - dielectric_stack + StackModel_Ptr - Description of any stack of coatings on the particle, usually specified by a file. - See dielectric_stack for more information. + Description of any stack of coatings on the particle. - none + No_StackModel diff --git a/docs/subbobvlieg.htm b/docs/subbobvlieg.htm index 1ed8224..6cd13f5 100644 --- a/docs/subbobvlieg.htm +++ b/docs/subbobvlieg.htm @@ -143,29 +143,23 @@

    Parameters:

    spherecoat - dielectric_stack + StackModel_Ptr Description of any stack of - coatings on the particle, usually specified by a file.
    - See dielectric_stack, for - more information.
    + coatings on the particle. - none + No_StackModel stack - dielectric_stack + StackModel_Ptr Description of any stack - of films deposited on the substrate, usually specified - by a file.
    - See dielectric_stack, for - more information.
    + of films deposited on the substrate. - no - films + No_StackModel diff --git a/docs/subsphere.htm b/docs/subsphere.htm index bee03ed..c6fe536 100644 --- a/docs/subsphere.htm +++ b/docs/subsphere.htm @@ -97,11 +97,9 @@

    Parameters:

    stack - dielectric_stack - Description of any stack of films deposited on the substrate, usually specified - by a file. See dielectric_stack, for more - information.
    - no films + StackModel_Ptr + Description of any stack of films deposited on the substrate. + No_StackModel diff --git a/docs/torrspar.htm b/docs/torrspar.htm index 01f76e7..086439c 100644 --- a/docs/torrspar.htm +++ b/docs/torrspar.htm @@ -117,20 +117,16 @@

    Parameters:

    stack - dielectric_stack + StackModel_Ptr Description of any stack - of films deposited on the substrate, usually specified - by a file. The films are assumed to be conformal and of + of films deposited on the substrate.The films are assumed to be conformal and of total thickness much less than the horizontal scale of the roughness.
    - See dielectric_stack, for - more information.
    (Inherited from Facet_BRDF_Model). - no - films + No_StackModel diff --git a/docs/transmit.htm b/docs/transmit.htm index 06def30..f69f20f 100644 --- a/docs/transmit.htm +++ b/docs/transmit.htm @@ -126,12 +126,11 @@ films - dielectric_stack + StackModel_Ptr - Description of the stack of films deposited on the bottom surface, usually specified by a file. - See dielectric_stack, for more information.
    + Description of the stack of films deposited on the bottom surface - no films + No_StackModel
    diff --git a/docs/uncorrelated_stack.htm b/docs/uncorrelated_stack.htm index abfc6b5..dd14f6e 100644 --- a/docs/uncorrelated_stack.htm +++ b/docs/uncorrelated_stack.htm @@ -118,16 +118,12 @@

    Parameters:

    stack - dielectric_stack + StackModel_Ptr Description of the stack - of films deposited on the substrate, usually specified - by a file.
    - See dielectric_stack, for - more information. + of films deposited on the substrate. - no - films + No_StackModel