From 278f55926e89abf1b786c76597dcb1e1f5c8a69d Mon Sep 17 00:00:00 2001 From: "Thomas A. Germer" Date: Mon, 26 Apr 2021 15:36:18 -0400 Subject: [PATCH] Remove PartialExtinction codes from Bobbert_Vlieger_BRDF_Model, Subsurface_Bobbert_Vlieger_BRDF_Model, Axisymmetric_Particle_BRDF_Model, and Subsurface_Axisymmetric_Particle and replace them with Specular(theta). Add 90 degree phase to Bobbert_Vlieger_BRDF_Model, Subsurface_Bobbert_Vlieger_BRDF_Model, Axisymmetric_Particle_BRDF_Model, and Subsurface_Axisymmetric_Particle, so that the optical theorem works the same as for free space scatterers. Correct JonesMatrix::hermitian() and add JonesMatrix::conj(). --- code/axipart.h | 7 +- code/axipart1.cpp | 164 +++++++++++++++++------------------------- code/axipart2.cpp | 2 +- code/bobvlieg.h | 16 +---- code/bobvlieg1.cpp | 172 +++++++++++++++++++-------------------------- code/bobvlieg2.cpp | 9 +-- code/mueller.h | 14 +++- code/sphprt.cpp | 15 +++- code/sphrscat.cpp | 39 +++++++++- code/sphrscat.h | 5 ++ code/subbobvlieg.h | 19 ++--- code/vector3d.cpp | 2 +- 12 files changed, 219 insertions(+), 245 deletions(-) diff --git a/code/axipart.h b/code/axipart.h index 4402317..3879207 100644 --- a/code/axipart.h +++ b/code/axipart.h @@ -47,8 +47,10 @@ namespace SCATMECH { public: Axisymmetric_Particle_BRDF_Model(); - COMPLEX PartialExtinctionS(double theta); - COMPLEX PartialExtinctionP(double theta); + MuellerMatrix Specular(double theta); + + //COMPLEX PartialExtinctionS(double theta); + //COMPLEX PartialExtinctionP(double theta); protected: void setup(); @@ -114,7 +116,6 @@ namespace SCATMECH { void calculate_Z(double thetas); void calculate_eIP(double phis); - COMPLEX PartialExtinction(double theta,int pol); void iterative_improvement(ScatterTMatrix& Ainv, ScatterTMatrix& A, std::vector& b,std::vector& x); diff --git a/code/axipart1.cpp b/code/axipart1.cpp index a7b666e..0bc37c2 100644 --- a/code/axipart1.cpp +++ b/code/axipart1.cpp @@ -581,112 +581,78 @@ namespace SCATMECH { MMAX = 0; } - COMPLEX + MuellerMatrix Axisymmetric_Particle_BRDF_Model:: - PartialExtinctionS(double theta) - { - return PartialExtinction(theta,1); - } - - COMPLEX - Axisymmetric_Particle_BRDF_Model:: - PartialExtinctionP(double theta) - { - return PartialExtinction(theta,0); - } - - COMPLEX - Axisymmetric_Particle_BRDF_Model:: - PartialExtinction(double theta,int pol) + Specular(double theta) { SETUP(); - // This function returns the partial extinction cross section for - // an incident angle of theta and polarization pol, where the "partial" - // means the following: when type=0 or type=2, the reflectance changes - // by a factor exp(-sigma*density/cos(theta)) and for type=1 or type=3, the - // transmittance changes by a factor exp(-sigma*density/cos(theta). The total - // extinction cross section is the sum of these values for type=0 and type=1, - // for light incident downward, or the sum of these values for type=2 and type=3, - // for light incident upward. The value of the partial extinction cross - // section can be negative, which is an indication that the reflectance - // or transmittance is increased by the presence of the particle. - // - - pol = pol ? 1 : 0; - vector& W = pol ? Ws : Wp; - vector& Z = pol ? Zs : Zp; + // This function returns the Mueller matrix specular reflectance (for type==0 or + // type==2) or the regular transmittance (for type==1 or type==3) for an incident + // angle of theta (in radians). The result is derived from the optical theorem + // and is only valid when density is suffiently low that multiple scattering + // between spheres can be neglected. switch (type) { - case 0: - { - set_geometry(theta,theta,0.); - COMPLEX e = E(W,Z); - 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; - } - break; - case 1: - { - double index = substrate.n(lambda); - double sint = sin(theta)/index; - if (sint>=1. || substrate.k(lambda)!=0) return 0.; - double thetat = asin(sint); - set_geometry(theta,thetat,0.); - COMPLEX e = E(W,Z); - 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]; - double T = norm(t)*factor*index; - return 4.*pi/k/sqrt(cube(index))/factor*COMPLEX(0.,-1.)*(e/t)*T; - } - break; - case 2: - { - set_geometry(theta,theta,0.); - COMPLEX e = E(W,Z); - double index = substrate.n(lambda); - - COMPLEX sint = sin(theta)*index; - COMPLEX cost = sqrt(1.-sqr(sint)); - if (imag(cost)<0) cost = -cost; - - COMPLEX r = stack->r21i(theta,lambda,substrate,vacuum)[pol]; - double R = norm(r); - - r *= exp(-2.*cI*qq*index*cos(theta)); - - return 4.*pi/k/index*COMPLEX(0.,-1.)*(e/r)*R; - } - break; - case 3: - { - double index = substrate.n(lambda); - double sint = sin(theta)*index; - COMPLEX cost = sqrt(1.-sqr(sint)); - if (imag(cost)<0) cost = -cost; - if (sint>=1.) return 0.; - double thetat = asin(sint); - - set_geometry(theta,thetat,0.); - COMPLEX e = E(W,Z); - - double factor = cos(theta)/real(cost); - COMPLEX phase = exp(cI*qq*(cost-index*cos(theta))); - - e /= phase; - 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; - } - break; - default: - error("Invalid type = " + to_string(type)); + case 0: + { + JonesMatrix X = JonesDSC(theta, theta, 0., 0.); + JonesMatrix r = stack->r12(theta, lambda, vacuum, substrate); + r *= exp(2. * cI * qq * cos(theta)); + MuellerMatrix sigma = (4. * pi / k) * ReCrossMueller(X, r); + + return MuellerMatrix(r) - sigma * (density / cos(theta)); + } + break; + case 1: + { + double index = substrate.n(lambda); + double sint = sin(theta) / index; + if (sint >= 1. || substrate.k(lambda) != 0) return MuellerZero(); + double thetat = asin(sint); + JonesMatrix X = JonesDSC(theta, thetat, 0., 0.); + COMPLEX phase = exp(cI * qq * (cos(theta) - index * cos(thetat))); + //double factor = cos(thetat)/cos(theta); + JonesMatrix t = stack->t12(theta, lambda, vacuum, substrate); + t *= phase; + MuellerMatrix sigma = (4. * pi / k / sqrt(index)) * ReCrossMueller(X, t); + + return MuellerMatrix(t) * (cos(thetat) / cos(theta) * index) - sigma * (density / cos(theta)); + } + break; + case 2: + { + double index = substrate.n(lambda); + JonesMatrix X = JonesDSC(theta, theta, 0., 0.); + JonesMatrix r = stack->r21i(theta, lambda, substrate, vacuum); + r *= exp(-2. * cI * index * qq * cos(theta)); + MuellerMatrix sigma = (4. * pi / k / index) * ReCrossMueller(X, r); + + return MuellerMatrix(r) - sigma * (density / cos(theta)); + } + break; + case 3: + { + double index = substrate.n(lambda); + double sint = sin(theta) * index; + COMPLEX cost = sqrt(1. - sqr(sint)); + if (imag(cost) < 0) cost = -cost; + if (sint >= 1.) return MuellerZero(); + double thetat = asin(sint); + JonesMatrix X = JonesDSC(theta, thetat, 0., 0.); + COMPLEX phase = exp(cI * qq * (cost - index * cos(theta))); + //double factor = cos(thetat)/cos(theta); + JonesMatrix t = stack->t21i(theta, lambda, substrate, vacuum); + t *= phase; + MuellerMatrix sigma = (4. * pi / k / sqrt(index)) * ReCrossMueller(X, t); + + return MuellerMatrix(t) * (cos(theta) / index / cos(thetat)) - sigma * (density / cos(theta)); + } + break; + default: + error("Invalid type = " + to_string(type)); } - return 0; + return MuellerZero(); } DEFINE_MODEL(Axisymmetric_Particle_BRDF_Model,Local_BRDF_Model, diff --git a/code/axipart2.cpp b/code/axipart2.cpp index 54082fd..1f9e602 100644 --- a/code/axipart2.cpp +++ b/code/axipart2.cpp @@ -474,7 +474,7 @@ namespace SCATMECH { } } - return E; + return COMPLEX(0, -1) * E; } diff --git a/code/bobvlieg.h b/code/bobvlieg.h index 881df3e..4153530 100644 --- a/code/bobvlieg.h +++ b/code/bobvlieg.h @@ -59,17 +59,9 @@ namespace SCATMECH { public: Bobbert_Vlieger_BRDF_Model(); - - // The real part of the following are the partial extinction cross - // sections for S or P polarization, respectively, where the "partial" - // means that when type=0, the downward reflectance is reduced by the - // result, when type=1, the downward transmittance is reduced by the result, - // when type=2, the upward reflectance is reduced by the result, and - // when type=3, the upward transmittance is reduced by the result. - // The total extinction cross section can be obtained by summing these - // values for type=0 and type=1 or for type=2 and type=3. - COMPLEX PartialExtinctionS(double theta); - COMPLEX PartialExtinctionP(double theta); + + // The Mueller matrix specular reflectance or regular transmittance... + MuellerMatrix Specular(double theta); protected: void setup(); @@ -135,8 +127,6 @@ namespace SCATMECH { void calculate_Z(double thetas); void calculate_eIP(double phis); - COMPLEX PartialExtinction(double theta,int pol); - void iterative_improvement(ScatterTMatrix& Ainv, ScatterTMatrix& A, std::vector& b,std::vector& x); }; diff --git a/code/bobvlieg1.cpp b/code/bobvlieg1.cpp index 6ec9424..1a60179 100644 --- a/code/bobvlieg1.cpp +++ b/code/bobvlieg1.cpp @@ -184,8 +184,12 @@ 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); + COMPLEX rs = stack->rs12(alpha, lambda, vacuum, substrate); + COMPLEX rp = stack->rp12(alpha, lambda, vacuum, substrate); + // Added 14 Jul 2020, because reflection coefficients at large complex + // angles sometime yield numerical errors from thick dielectrics + reflect_s[j] = (rs == rs) ? rs : 0.; + reflect_p[j] = (rp == rp) ? rp : 0.; } // Calculate U, V, d+, and d- for each angle and (l,m) combination... @@ -810,112 +814,78 @@ namespace SCATMECH { old_phis=-1000; } - COMPLEX + MuellerMatrix Bobbert_Vlieger_BRDF_Model:: - PartialExtinctionS(double theta) - { - return PartialExtinction(theta,1); - } - - COMPLEX - Bobbert_Vlieger_BRDF_Model:: - PartialExtinctionP(double theta) - { - return PartialExtinction(theta,0); - } - - COMPLEX - Bobbert_Vlieger_BRDF_Model:: - PartialExtinction(double theta,int pol) + Specular(double theta) { SETUP(); - // This function returns the partial extinction cross section for - // an incident angle of theta and polarization pol, where the "partial" - // means the following: when type=0 or type=2, the reflectance changes - // by a factor exp(-sigma*density/cos(theta)) and for type=1 or type=3, the - // transmittance changes by a factor exp(-sigma*density/cos(theta)). The total - // extinction cross section is the sum of these values for type=0 and type=1, - // for light incident downward, or the sum of these values for type=2 and type=3, - // for light incident upward. The value of the partial extinction cross - // section can be negative, which is an indication that the reflectance - // or transmittance is increased by the presence of the particle. - // - - pol = pol ? 1 : 0; - vector& W = pol ? Ws : Wp; - vector& Z = pol ? Zs : Zp; + // This function returns the Mueller matrix specular reflectance (for type==0 or + // type==2) or the regular transmittance (for type==1 or type==3) for an incident + // angle of theta (in radians). The result is derived from the optical theorem + // and is only valid when density is suffiently low that multiple scattering + // between spheres can be neglected. switch (type) { - case 0: - { - set_geometry(theta,theta,0.); - COMPLEX e = E(W,Z); - 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; - } - break; - case 1: - { - double index = substrate.n(lambda); - double sint = sin(theta)/index; - if (sint>=1. || substrate.k(lambda)!=0) return 0.; - double thetat = asin(sint); - set_geometry(theta,thetat,0.); - COMPLEX e = E(W,Z); - 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]; - double T = norm(t)*factor*index; - return 4.*pi/k/sqrt(cube(index))/factor*COMPLEX(0.,-1.)*(e/t)*T; - } - break; - case 2: - { - set_geometry(theta,theta,0.); - COMPLEX e = E(W,Z); - double index = substrate.n(lambda); - - COMPLEX sint = sin(theta)*index; - COMPLEX cost = sqrt(1.-sqr(sint)); - if (imag(cost)<0) cost = -cost; - - COMPLEX r = stack->r21i(theta,lambda,substrate,vacuum)[pol]; - double R = norm(r); - - r *= exp(-2.*cI*qq*index*cos(theta)); - - return 4.*pi/k/index*COMPLEX(0.,-1.)*(e/r)*R; - } - break; - case 3: - { - double index = substrate.n(lambda); - double sint = sin(theta)*index; - COMPLEX cost = sqrt(1.-sqr(sint)); - if (imag(cost)<0) cost = -cost; - if (sint>=1.) return 0.; - double thetat = asin(sint); - - set_geometry(theta,thetat,0.); - COMPLEX e = E(W,Z); - - double factor = cos(theta)/real(cost); - COMPLEX phase = exp(cI*qq*(cost-index*cos(theta))); - - e /= phase; - 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; - } - break; - default: - error("Invalid type = " + to_string(type)); + case 0: + { + JonesMatrix X = JonesDSC(theta, theta, 0., 0.); + JonesMatrix r = stack->r12(theta, lambda, vacuum, substrate); + r *= exp(2. * cI * qq * cos(theta)); + MuellerMatrix sigma = (4. * pi / k) * ReCrossMueller(X, r); + + return MuellerMatrix(r) - sigma * (density / cos(theta)); + } + break; + case 1: + { + double index = substrate.n(lambda); + double sint = sin(theta) / index; + if (sint >= 1. || substrate.k(lambda) != 0) return MuellerZero(); + double thetat = asin(sint); + JonesMatrix X = JonesDSC(theta, thetat, 0., 0.); + COMPLEX phase = exp(cI * qq * (cos(theta) - index * cos(thetat))); + //double factor = cos(thetat)/cos(theta); + JonesMatrix t = stack->t12(theta, lambda, vacuum, substrate); + t *= phase; + MuellerMatrix sigma = (4. * pi / k / sqrt(index)) * ReCrossMueller(X, t); + + return MuellerMatrix(t) * (cos(thetat) / cos(theta) * index) - sigma * (density / cos(theta)); + } + break; + case 2: + { + double index = substrate.n(lambda); + JonesMatrix X = JonesDSC(theta, theta, 0., 0.); + JonesMatrix r = stack->r21i(theta, lambda, substrate, vacuum); + r *= exp(-2. * cI * index * qq * cos(theta)); + MuellerMatrix sigma = (4. * pi / k / index) * ReCrossMueller(X, r); + + return MuellerMatrix(r) - sigma * (density / cos(theta)); + } + break; + case 3: + { + double index = substrate.n(lambda); + double sint = sin(theta) * index; + COMPLEX cost = sqrt(1. - sqr(sint)); + if (imag(cost) < 0) cost = -cost; + if (sint >= 1.) return MuellerZero(); + double thetat = asin(sint); + JonesMatrix X = JonesDSC(theta, thetat, 0., 0.); + COMPLEX phase = exp(cI * qq * (cost - index * cos(theta))); + //double factor = cos(thetat)/cos(theta); + JonesMatrix t = stack->t21i(theta, lambda, substrate, vacuum); + t *= phase; + MuellerMatrix sigma = (4. * pi / k / sqrt(index)) * ReCrossMueller(X, t); + + return MuellerMatrix(t) * (cos(theta) / index / cos(thetat)) - sigma * (density / cos(theta)); + } + break; + default: + error("Invalid type = " + to_string(type)); } - return 0; + return MuellerZero(); } DEFINE_MODEL(Bobbert_Vlieger_BRDF_Model,Local_BRDF_Model, diff --git a/code/bobvlieg2.cpp b/code/bobvlieg2.cpp index f673f02..7693d57 100644 --- a/code/bobvlieg2.cpp +++ b/code/bobvlieg2.cpp @@ -19,13 +19,10 @@ using namespace std; - namespace SCATMECH { using namespace BobVlieg_Supp; - - // // factorial_list is a lookup table for n! // @@ -1003,7 +1000,7 @@ namespace SCATMECH { COMPLEX _cost=cos(pi-thetas); 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 phase = exp(2.*cI*qq*_cost); COMPLEX rpphase = rp*phase; COMPLEX rsphase = rs*phase; @@ -1046,7 +1043,7 @@ namespace SCATMECH { if (imag(cost)<0) cost = -cost; 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))); + 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... double factor = abs(COMPLEX(cos(pi-thetas))/cost*COMPLEX(sqrt(cube(index)))); @@ -1131,7 +1128,7 @@ namespace SCATMECH { } } - return E; + return COMPLEX(0,-1)*E; } } // namespace SCATMECH diff --git a/code/mueller.h b/code/mueller.h index 87858ab..b564b84 100644 --- a/code/mueller.h +++ b/code/mueller.h @@ -415,10 +415,20 @@ namespace SCATMECH { /// Return the Hermitian transpose /// JonesMatrix hermitian() const { - return JonesMatrix(j[0],j[1],conj(j[3]),conj(j[2])); + return JonesMatrix(std::conj(j[0]), std::conj(j[1]), std::conj(j[3]), std::conj(j[2])); } friend JonesMatrix hermitian(const JonesMatrix& j) { - return JonesMatrix(j[0],j[1],conj(j[3]),conj(j[2])); + return JonesMatrix(std::conj(j[0]), std::conj(j[1]), std::conj(j[3]), std::conj(j[2])); + } + + /// + /// Return the complex conjugate + /// + JonesMatrix conj() const { + return JonesMatrix(std::conj(j[0]), std::conj(j[1]), std::conj(j[2]), std::conj(j[3])); + } + friend JonesMatrix conj(const JonesMatrix& j) { + return JonesMatrix(std::conj(j[0]), std::conj(j[1]), std::conj(j[2]), std::conj(j[3])); } /// diff --git a/code/sphprt.cpp b/code/sphprt.cpp index 2851a51..cf6c0f6 100644 --- a/code/sphprt.cpp +++ b/code/sphprt.cpp @@ -219,8 +219,11 @@ namespace SCATMECH { // The Jacobian accounts for the solid angle changing through the interface... double jacobian = sqr(n)*cos(thetas)/cos_thetas_outside; + // This phase term accounts for the propagation from the sphere to the surface... + COMPLEX phase = exp(COMPLEX(0, 1) * k * distance * (cos_thetas_outside - n * cos(thetas))); + // The result... - return scatter*sqrt(jacobian); + return scatter*sqrt(jacobian)*phase; } } else { // is_backward() ... @@ -285,8 +288,12 @@ namespace SCATMECH { COMPLEX common=jacobian/sqr(k); + // This phase term accounts for the propagation from the sphere to the surface... + COMPLEX phasei = exp(COMPLEX(0, 1) * k * distance * (cos_thetai_at_part - n * cos(thetai))); + COMPLEX phases = exp(COMPLEX(0, 1) * k * distance * (cos_thetas_at_part - n * cos(thetas))); + // Return the whole value with factors common to all elements... - return scatter*sqrt(common); + return scatter*sqrt(common)*phasei*phases; } else { // is_transmission() @@ -368,8 +375,10 @@ namespace SCATMECH { // The transmission coefficient through the interface... JonesMatrix ti = stack->t21i(thetai,lambda,substrate,vacuum)/(COMPLEX)sqrt(n); + COMPLEX phasei = exp(COMPLEX(0, 1) * k * distance * (cos_thetai_at_part - n * cos(thetai))); + // The total scatter... - JonesMatrix scatter = (scatter_direct+scatter_indirect)*ti/k; + JonesMatrix scatter = (scatter_direct+scatter_indirect)*ti/k*phasei; return scatter; diff --git a/code/sphrscat.cpp b/code/sphrscat.cpp index 7c88bd5..4a91924 100644 --- a/code/sphrscat.cpp +++ b/code/sphrscat.cpp @@ -16,9 +16,46 @@ using namespace std; - namespace SCATMECH { + MuellerMatrix + Free_Space_Scatterer:: + extinction(const Vector& v) + { + // + // Calculates the extinction cross section from the optical theorem + // See M. A. Karam, "Polarimetric Optical Theorem," J. Opt. Soc. Am. A 15 (1) 196-201 (1998) + // And my notes of 25 March 2020 + // + SETUP(); + + JonesMatrix J = jones(v, v); + + MuellerMatrix result; + result[0][0] = real(J.SS() + J.PP()); + result[0][1] = real(J.SS() - J.PP()); + result[0][2] = real(J.PS() + J.SP()); + result[0][3] = imag(J.SP() - J.PS()); + + result[1][0] = result[0][1]; + result[1][1] = result[0][0]; + result[1][2] = real(J.PS() - J.SP()); + result[1][3] = -imag(J.PS() + J.SP()); + + result[2][0] = result[0][2]; + result[2][1] = -result[1][2]; + result[2][2] = result[0][0]; + result[2][3] = imag(J.SS() - J.PP()); + + result[3][0] = result[0][3]; + result[3][1] = -result[1][3]; + result[3][2] = -result[2][3]; + result[3][3] = result[0][0]; + + double k = 2 * pi / lambda * medium.n(lambda); + return (2 * pi / sqr(k)) * result; + } + DEFINE_VIRTUAL_MODEL(Free_Space_Scatterer,Model, "Generalized free-space scatterer"); diff --git a/code/sphrscat.h b/code/sphrscat.h index 1de5a5a..b5a2a17 100644 --- a/code/sphrscat.h +++ b/code/sphrscat.h @@ -41,6 +41,11 @@ namespace SCATMECH { // The Vectors kin and kout are unit vectors in the incident and scattering directions, respectively. // virtual JonesMatrix jones(const Vector& kin,const Vector& kout) = 0; + + // + // Exctinction cross section... + // + MuellerMatrix extinction(const Vector& k); }; typedef Model_Ptr Free_Space_Scatterer_Ptr; diff --git a/code/subbobvlieg.h b/code/subbobvlieg.h index dc138f0..a9393db 100644 --- a/code/subbobvlieg.h +++ b/code/subbobvlieg.h @@ -34,14 +34,8 @@ namespace SCATMECH { public: - COMPLEX PartialExtinctionS(double theta) { - SETUP(); - return model.PartialExtinctionS(theta); - } - - COMPLEX PartialExtinctionP(double theta) { - SETUP(); - return model.PartialExtinctionP(theta); + MuellerMatrix Specular(double theta) { + return model.Specular(theta); } protected: @@ -67,14 +61,9 @@ namespace SCATMECH { DECLARE_PARAMETER(int,improve); public: - COMPLEX PartialExtinctionS(double theta) { - SETUP(); - return model.PartialExtinctionS(theta); - } - COMPLEX PartialExtinctionP(double theta) { - SETUP(); - return model.PartialExtinctionP(theta); + MuellerMatrix Specular(double theta) { + return model.Specular(theta); } protected: diff --git a/code/vector3d.cpp b/code/vector3d.cpp index f4e80fa..bdcf313 100644 --- a/code/vector3d.cpp +++ b/code/vector3d.cpp @@ -28,7 +28,7 @@ namespace SCATMECH { // _norm to a small number, rather than zero... Vector3D _a = unit(a); Vector3D _b = unit(b); - + Vector3D temp=cross(_a,_b); double _norm=Norm(temp); if (_norm>1E-10) return temp/_norm;