Skip to content

Commit

Permalink
Added StackModel and Legendre_Phase_Function, and a few other bug fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
@thomas-germer
thomas-germer committed Sep 29, 2017
1 parent ff980d1 commit a854b09
Show file tree
Hide file tree
Showing 88 changed files with 2,603 additions and 1,394 deletions.
1 change: 0 additions & 1 deletion code/Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,6 @@
#//** Phone: (301) 975-2876
#//** Email: [email protected]
#//**
#//** Version: 7.00 (January 2015)
#//**
#//******************************************************************************

Expand Down
15 changes: 8 additions & 7 deletions code/allrough.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
Expand All @@ -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);
}
Expand Down Expand Up @@ -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.);
Expand All @@ -118,7 +119,7 @@ namespace SCATMECH {
vector<double> a(N);
vector<double> PSDint(N);
for (int i=0; i<N; ++i) {
double argument = fourpisqr*nu*qpow*stack.get_t()[i];
double argument = fourpisqr*nu*qpow*stack->get_t()[i];
// The replication factor between the layers...
a[i] = exp(-argument);
// The intrinsic roughness of the layer...
Expand Down Expand Up @@ -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);
Expand Down
10 changes: 5 additions & 5 deletions code/allrough.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"

Expand All @@ -26,7 +26,7 @@ namespace SCATMECH {
public:

DECLARE_MODEL();
DECLARE_PARAMETER(dielectric_stack,stack);
DECLARE_PARAMETER(StackModel_Ptr,stack);

protected:

Expand All @@ -43,7 +43,7 @@ namespace SCATMECH {
public:

DECLARE_MODEL();
DECLARE_PARAMETER(dielectric_stack,stack);
DECLARE_PARAMETER(StackModel_Ptr,stack);

protected:

Expand All @@ -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);
Expand Down
4 changes: 2 additions & 2 deletions code/askuser.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
Expand Down
2 changes: 1 addition & 1 deletion code/axipart.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
18 changes: 9 additions & 9 deletions code/axipart1.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<COMPLEX> reflect_s(npts);
vector<COMPLEX> reflect_p(npts);
Expand All @@ -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...
Expand Down Expand Up @@ -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;
Expand All @@ -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;
}
Expand All @@ -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));
Expand All @@ -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;
}
Expand All @@ -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);
Expand Down
16 changes: 8 additions & 8 deletions code/axipart2.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.);

Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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.);

Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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...
Expand Down
6 changes: 3 additions & 3 deletions code/bobvlieg.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)...
Expand Down Expand Up @@ -143,7 +143,7 @@ namespace SCATMECH {
std::vector<COMPLEX>& b,std::vector<COMPLEX>& x);
};

class Gauss_Laguerre_Integration {
class Gauss_Laguerre_Integration {
public:
static double* weights[];
static double* zeros[];
Expand Down
40 changes: 20 additions & 20 deletions code/bobvlieg1.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,23 +71,23 @@ 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...
T = -(ms*psi(n,ms*ys)*(psi_(n,ys)+T*chi_(n,ys))-psi_(n,ms*ys)*(psi(n,ys)+T*chi(n,ys)))/
(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...
Expand All @@ -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;

Expand Down Expand Up @@ -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<COMPLEX> reflect_s(npts);
vector<COMPLEX> reflect_p(npts);
Expand All @@ -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...
Expand Down Expand Up @@ -479,9 +479,9 @@ namespace SCATMECH {

r0 = radius;
d = 0.;
for (int i=0; i<spherecoat.get_n(); ++i) {
r0 -= spherecoat.get_t()[i];
d += spherecoat.get_t()[i];
for (int i=0; i<spherecoat->get_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.");

Expand Down Expand Up @@ -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;
Expand All @@ -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;
}
Expand All @@ -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));
Expand All @@ -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;
}
Expand All @@ -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);
Expand Down
Loading

0 comments on commit a854b09

Please sign in to comment.