|
@@ -45,6 +45,8 @@
|
|
|
#include <stdexcept>
|
|
|
#include <vector>
|
|
|
|
|
|
+using namespace std;
|
|
|
+
|
|
|
namespace nmie {
|
|
|
//helpers
|
|
|
template<class T> inline T pow2(const T value) {return value*value;}
|
|
@@ -72,10 +74,10 @@ namespace nmie {
|
|
|
// Return value: //
|
|
|
// Number of multipolar expansion terms used for the calculations //
|
|
|
//**********************************************************************************//
|
|
|
- int ScattCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const int nmax, std::vector<std::complex<double> >& an, std::vector<std::complex<double> >& bn) {
|
|
|
+ int ScattCoeffs(const unsigned int L, const int pl, vector<double>& x, vector<complex<double> >& m, const int nmax, vector<complex<double> >& an, vector<complex<double> >& bn) {
|
|
|
|
|
|
if (x.size() != L || m.size() != L)
|
|
|
- throw std::invalid_argument("Declared number of layers do not fit x and m!");
|
|
|
+ throw invalid_argument("Declared number of layers do not fit x and m!");
|
|
|
try {
|
|
|
MultiLayerMie ml_mie;
|
|
|
ml_mie.SetAllLayersSize(x);
|
|
@@ -89,10 +91,10 @@ namespace nmie {
|
|
|
bn = ml_mie.GetBn();
|
|
|
|
|
|
return ml_mie.GetMaxTerms();
|
|
|
- } catch(const std::invalid_argument& ia) {
|
|
|
+ } catch(const invalid_argument& ia) {
|
|
|
// Will catch if ml_mie fails or other errors.
|
|
|
- std::cerr << "Invalid argument: " << ia.what() << std::endl;
|
|
|
- throw std::invalid_argument(ia);
|
|
|
+ cerr << "Invalid argument: " << ia.what() << endl;
|
|
|
+ throw invalid_argument(ia);
|
|
|
return -1;
|
|
|
}
|
|
|
return 0;
|
|
@@ -127,12 +129,12 @@ namespace nmie {
|
|
|
// Return value: //
|
|
|
// Number of multipolar expansion terms used for the calculations //
|
|
|
//**********************************************************************************//
|
|
|
- int nMie(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const unsigned int nTheta, std::vector<double>& Theta, const int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
|
|
|
+ int nMie(const unsigned int L, const int pl, vector<double>& x, vector<complex<double> >& m, const unsigned int nTheta, vector<double>& Theta, const int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, vector<complex<double> >& S1, vector<complex<double> >& S2) {
|
|
|
|
|
|
if (x.size() != L || m.size() != L)
|
|
|
- throw std::invalid_argument("Declared number of layers do not fit x and m!");
|
|
|
+ throw invalid_argument("Declared number of layers do not fit x and m!");
|
|
|
if (Theta.size() != nTheta)
|
|
|
- throw std::invalid_argument("Declared number of sample for Theta is not correct!");
|
|
|
+ throw invalid_argument("Declared number of sample for Theta is not correct!");
|
|
|
try {
|
|
|
MultiLayerMie ml_mie;
|
|
|
ml_mie.SetAllLayersSize(x);
|
|
@@ -154,10 +156,10 @@ namespace nmie {
|
|
|
S2 = ml_mie.GetS2();
|
|
|
|
|
|
return ml_mie.GetMaxTerms();
|
|
|
- } catch(const std::invalid_argument& ia) {
|
|
|
+ } catch(const invalid_argument& ia) {
|
|
|
// Will catch if ml_mie fails or other errors.
|
|
|
- std::cerr << "Invalid argument: " << ia.what() << std::endl;
|
|
|
- throw std::invalid_argument(ia);
|
|
|
+ cerr << "Invalid argument: " << ia.what() << endl;
|
|
|
+ throw invalid_argument(ia);
|
|
|
return -1;
|
|
|
}
|
|
|
return 0;
|
|
@@ -191,7 +193,7 @@ namespace nmie {
|
|
|
// Return value: //
|
|
|
// Number of multipolar expansion terms used for the calculations //
|
|
|
//**********************************************************************************//
|
|
|
- int nMie(const unsigned int L, std::vector<double>& x, std::vector<std::complex<double> >& m, const unsigned int nTheta, std::vector<double>& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
|
|
|
+ int nMie(const unsigned int L, vector<double>& x, vector<complex<double> >& m, const unsigned int nTheta, vector<double>& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, vector<complex<double> >& S1, vector<complex<double> >& S2) {
|
|
|
return nmie::nMie(L, -1, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
|
|
|
}
|
|
|
|
|
@@ -223,7 +225,7 @@ namespace nmie {
|
|
|
// Return value: //
|
|
|
// Number of multipolar expansion terms used for the calculations //
|
|
|
//**********************************************************************************//
|
|
|
- int nMie(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const unsigned int nTheta, std::vector<double>& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
|
|
|
+ int nMie(const unsigned int L, const int pl, vector<double>& x, vector<complex<double> >& m, const unsigned int nTheta, vector<double>& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, vector<complex<double> >& S1, vector<complex<double> >& S2) {
|
|
|
return nmie::nMie(L, pl, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
|
|
|
}
|
|
|
|
|
@@ -257,7 +259,7 @@ namespace nmie {
|
|
|
// Return value: //
|
|
|
// Number of multipolar expansion terms used for the calculations //
|
|
|
//**********************************************************************************//
|
|
|
- int nMie(const unsigned int L, std::vector<double>& x, std::vector<std::complex<double> >& m, const unsigned int nTheta, std::vector<double>& Theta, const int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
|
|
|
+ int nMie(const unsigned int L, vector<double>& x, vector<complex<double> >& m, const unsigned int nTheta, vector<double>& Theta, const int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, vector<complex<double> >& S1, vector<complex<double> >& S2) {
|
|
|
return nmie::nMie(L, -1, x, m, nTheta, Theta, nmax, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
|
|
|
}
|
|
|
|
|
@@ -284,18 +286,18 @@ namespace nmie {
|
|
|
// Return value: //
|
|
|
// Number of multipolar expansion terms used for the calculations //
|
|
|
//**********************************************************************************//
|
|
|
- int nField(const unsigned int L, const int pl, const std::vector<double>& x, const std::vector<std::complex<double> >& m, const int nmax, const unsigned int ncoord, const std::vector<double>& Xp_vec, const std::vector<double>& Yp_vec, const std::vector<double>& Zp_vec, std::vector<std::vector<std::complex<double> > >& E, std::vector<std::vector<std::complex<double> > >& H) {
|
|
|
+ int nField(const unsigned int L, const int pl, const vector<double>& x, const vector<complex<double> >& m, const int nmax, const unsigned int ncoord, const vector<double>& Xp_vec, const vector<double>& Yp_vec, const vector<double>& Zp_vec, vector<vector<complex<double> > >& E, vector<vector<complex<double> > >& H) {
|
|
|
if (x.size() != L || m.size() != L)
|
|
|
- throw std::invalid_argument("Declared number of layers do not fit x and m!");
|
|
|
+ throw invalid_argument("Declared number of layers do not fit x and m!");
|
|
|
if (Xp_vec.size() != ncoord || Yp_vec.size() != ncoord || Zp_vec.size() != ncoord
|
|
|
|| E.size() != ncoord || H.size() != ncoord)
|
|
|
- throw std::invalid_argument("Declared number of coords do not fit Xp, Yp, Zp, E, or H!");
|
|
|
+ throw invalid_argument("Declared number of coords do not fit Xp, Yp, Zp, E, or H!");
|
|
|
for (auto f:E)
|
|
|
if (f.size() != 3)
|
|
|
- throw std::invalid_argument("Field E is not 3D!");
|
|
|
+ throw invalid_argument("Field E is not 3D!");
|
|
|
for (auto f:H)
|
|
|
if (f.size() != 3)
|
|
|
- throw std::invalid_argument("Field H is not 3D!");
|
|
|
+ throw invalid_argument("Field H is not 3D!");
|
|
|
try {
|
|
|
MultiLayerMie ml_mie;
|
|
|
//ml_mie.SetPECLayer(pl); // TODO add PEC layer to field plotting
|
|
@@ -307,10 +309,10 @@ namespace nmie {
|
|
|
H = ml_mie.GetFieldH();
|
|
|
|
|
|
return ml_mie.GetMaxTerms();
|
|
|
- } catch(const std::invalid_argument& ia) {
|
|
|
+ } catch(const invalid_argument& ia) {
|
|
|
// Will catch if ml_mie fails or other errors.
|
|
|
- std::cerr << "Invalid argument: " << ia.what() << std::endl;
|
|
|
- throw std::invalid_argument(ia);
|
|
|
+ cerr << "Invalid argument: " << ia.what() << endl;
|
|
|
+ throw invalid_argument(ia);
|
|
|
return - 1;
|
|
|
}
|
|
|
return 0;
|
|
@@ -322,7 +324,7 @@ namespace nmie {
|
|
|
// ********************************************************************** //
|
|
|
double MultiLayerMie::GetQext() {
|
|
|
if (!isMieCalculated_)
|
|
|
- throw std::invalid_argument("You should run calculations before result request!");
|
|
|
+ throw invalid_argument("You should run calculations before result request!");
|
|
|
return Qext_;
|
|
|
}
|
|
|
|
|
@@ -332,7 +334,7 @@ namespace nmie {
|
|
|
// ********************************************************************** //
|
|
|
double MultiLayerMie::GetQabs() {
|
|
|
if (!isMieCalculated_)
|
|
|
- throw std::invalid_argument("You should run calculations before result request!");
|
|
|
+ throw invalid_argument("You should run calculations before result request!");
|
|
|
return Qabs_;
|
|
|
}
|
|
|
|
|
@@ -342,7 +344,7 @@ namespace nmie {
|
|
|
// ********************************************************************** //
|
|
|
double MultiLayerMie::GetQsca() {
|
|
|
if (!isMieCalculated_)
|
|
|
- throw std::invalid_argument("You should run calculations before result request!");
|
|
|
+ throw invalid_argument("You should run calculations before result request!");
|
|
|
return Qsca_;
|
|
|
}
|
|
|
|
|
@@ -352,7 +354,7 @@ namespace nmie {
|
|
|
// ********************************************************************** //
|
|
|
double MultiLayerMie::GetQbk() {
|
|
|
if (!isMieCalculated_)
|
|
|
- throw std::invalid_argument("You should run calculations before result request!");
|
|
|
+ throw invalid_argument("You should run calculations before result request!");
|
|
|
return Qbk_;
|
|
|
}
|
|
|
|
|
@@ -362,7 +364,7 @@ namespace nmie {
|
|
|
// ********************************************************************** //
|
|
|
double MultiLayerMie::GetQpr() {
|
|
|
if (!isMieCalculated_)
|
|
|
- throw std::invalid_argument("You should run calculations before result request!");
|
|
|
+ throw invalid_argument("You should run calculations before result request!");
|
|
|
return Qpr_;
|
|
|
}
|
|
|
|
|
@@ -372,7 +374,7 @@ namespace nmie {
|
|
|
// ********************************************************************** //
|
|
|
double MultiLayerMie::GetAsymmetryFactor() {
|
|
|
if (!isMieCalculated_)
|
|
|
- throw std::invalid_argument("You should run calculations before result request!");
|
|
|
+ throw invalid_argument("You should run calculations before result request!");
|
|
|
return asymmetry_factor_;
|
|
|
}
|
|
|
|
|
@@ -382,7 +384,7 @@ namespace nmie {
|
|
|
// ********************************************************************** //
|
|
|
double MultiLayerMie::GetAlbedo() {
|
|
|
if (!isMieCalculated_)
|
|
|
- throw std::invalid_argument("You should run calculations before result request!");
|
|
|
+ throw invalid_argument("You should run calculations before result request!");
|
|
|
return albedo_;
|
|
|
}
|
|
|
|
|
@@ -390,9 +392,9 @@ namespace nmie {
|
|
|
// ********************************************************************** //
|
|
|
// Returns previously calculated S1 //
|
|
|
// ********************************************************************** //
|
|
|
- std::vector<std::complex<double> > MultiLayerMie::GetS1() {
|
|
|
+ vector<complex<double> > MultiLayerMie::GetS1() {
|
|
|
if (!isMieCalculated_)
|
|
|
- throw std::invalid_argument("You should run calculations before result request!");
|
|
|
+ throw invalid_argument("You should run calculations before result request!");
|
|
|
return S1_;
|
|
|
}
|
|
|
|
|
@@ -400,9 +402,9 @@ namespace nmie {
|
|
|
// ********************************************************************** //
|
|
|
// Returns previously calculated S2 //
|
|
|
// ********************************************************************** //
|
|
|
- std::vector<std::complex<double> > MultiLayerMie::GetS2() {
|
|
|
+ vector<complex<double> > MultiLayerMie::GetS2() {
|
|
|
if (!isMieCalculated_)
|
|
|
- throw std::invalid_argument("You should run calculations before result request!");
|
|
|
+ throw invalid_argument("You should run calculations before result request!");
|
|
|
return S2_;
|
|
|
}
|
|
|
|
|
@@ -410,7 +412,7 @@ namespace nmie {
|
|
|
// ********************************************************************** //
|
|
|
// Modify scattering (theta) angles //
|
|
|
// ********************************************************************** //
|
|
|
- void MultiLayerMie::SetAngles(const std::vector<double>& angles) {
|
|
|
+ void MultiLayerMie::SetAngles(const vector<double>& angles) {
|
|
|
isExpCoeffsCalc_ = false;
|
|
|
isScaCoeffsCalc_ = false;
|
|
|
isMieCalculated_ = false;
|
|
@@ -421,7 +423,7 @@ namespace nmie {
|
|
|
// ********************************************************************** //
|
|
|
// Modify size of all layers //
|
|
|
// ********************************************************************** //
|
|
|
- void MultiLayerMie::SetAllLayersSize(const std::vector<double>& layer_size) {
|
|
|
+ void MultiLayerMie::SetAllLayersSize(const vector<double>& layer_size) {
|
|
|
isExpCoeffsCalc_ = false;
|
|
|
isScaCoeffsCalc_ = false;
|
|
|
isMieCalculated_ = false;
|
|
@@ -429,9 +431,9 @@ namespace nmie {
|
|
|
double prev_layer_size = 0.0;
|
|
|
for (auto curr_layer_size : layer_size) {
|
|
|
if (curr_layer_size <= 0.0)
|
|
|
- throw std::invalid_argument("Size parameter should be positive!");
|
|
|
+ throw invalid_argument("Size parameter should be positive!");
|
|
|
if (prev_layer_size > curr_layer_size)
|
|
|
- throw std::invalid_argument
|
|
|
+ throw invalid_argument
|
|
|
("Size parameter for next layer should be larger than the previous one!");
|
|
|
prev_layer_size = curr_layer_size;
|
|
|
size_param_.push_back(curr_layer_size);
|
|
@@ -442,7 +444,7 @@ namespace nmie {
|
|
|
// ********************************************************************** //
|
|
|
// Modify refractive index of all layers //
|
|
|
// ********************************************************************** //
|
|
|
- void MultiLayerMie::SetAllLayersIndex(const std::vector< std::complex<double> >& index) {
|
|
|
+ void MultiLayerMie::SetAllLayersIndex(const vector< complex<double> >& index) {
|
|
|
isExpCoeffsCalc_ = false;
|
|
|
isScaCoeffsCalc_ = false;
|
|
|
isMieCalculated_ = false;
|
|
@@ -453,11 +455,11 @@ namespace nmie {
|
|
|
// ********************************************************************** //
|
|
|
// Modify coordinates for field calculation //
|
|
|
// ********************************************************************** //
|
|
|
- void MultiLayerMie::SetFieldCoords(const std::vector< std::vector<double> >& coords) {
|
|
|
+ void MultiLayerMie::SetFieldCoords(const vector< vector<double> >& coords) {
|
|
|
if (coords.size() != 3)
|
|
|
- throw std::invalid_argument("Error! Wrong dimension of field monitor points!");
|
|
|
+ throw invalid_argument("Error! Wrong dimension of field monitor points!");
|
|
|
if (coords[0].size() != coords[1].size() || coords[0].size() != coords[2].size())
|
|
|
- throw std::invalid_argument("Error! Missing coordinates for field monitor points!");
|
|
|
+ throw invalid_argument("Error! Missing coordinates for field monitor points!");
|
|
|
coords_ = coords;
|
|
|
}
|
|
|
|
|
@@ -470,7 +472,7 @@ namespace nmie {
|
|
|
isScaCoeffsCalc_ = false;
|
|
|
isMieCalculated_ = false;
|
|
|
if (layer_position < 0 && layer_position != -1)
|
|
|
- throw std::invalid_argument("Error! Layers are numbered from 0!");
|
|
|
+ throw invalid_argument("Error! Layers are numbered from 0!");
|
|
|
PEC_layer_position_ = layer_position;
|
|
|
}
|
|
|
|
|
@@ -538,21 +540,21 @@ namespace nmie {
|
|
|
// ********************************************************************** //
|
|
|
void MultiLayerMie::calcNmax(unsigned int first_layer) {
|
|
|
int ri, riM1;
|
|
|
- const std::vector<double>& x = size_param_;
|
|
|
- const std::vector<std::complex<double> >& m = refractive_index_;
|
|
|
+ const vector<double>& x = size_param_;
|
|
|
+ const vector<complex<double> >& m = refractive_index_;
|
|
|
calcNstop(); // Set initial nmax_ value
|
|
|
for (unsigned int i = first_layer; i < x.size(); i++) {
|
|
|
if (static_cast<int>(i) > PEC_layer_position_) // static_cast used to avoid warning
|
|
|
- ri = round(std::abs(x[i]*m[i]));
|
|
|
+ ri = round(abs(x[i]*m[i]));
|
|
|
else
|
|
|
ri = 0;
|
|
|
- nmax_ = std::max(nmax_, ri);
|
|
|
+ nmax_ = max(nmax_, ri);
|
|
|
// first layer is pec, if pec is present
|
|
|
if ((i > first_layer) && (static_cast<int>(i - 1) > PEC_layer_position_))
|
|
|
- riM1 = round(std::abs(x[i - 1]* m[i]));
|
|
|
+ riM1 = round(abs(x[i - 1]* m[i]));
|
|
|
else
|
|
|
riM1 = 0;
|
|
|
- nmax_ = std::max(nmax_, riM1);
|
|
|
+ nmax_ = max(nmax_, riM1);
|
|
|
}
|
|
|
nmax_ += 15; // Final nmax_ value
|
|
|
}
|
|
@@ -561,12 +563,12 @@ namespace nmie {
|
|
|
// ********************************************************************** //
|
|
|
// Calculate an - equation (5) //
|
|
|
// ********************************************************************** //
|
|
|
- std::complex<double> MultiLayerMie::calc_an(int n, double XL, std::complex<double> Ha, std::complex<double> mL,
|
|
|
- std::complex<double> PsiXL, std::complex<double> ZetaXL,
|
|
|
- std::complex<double> PsiXLM1, std::complex<double> ZetaXLM1) {
|
|
|
+ complex<double> MultiLayerMie::calc_an(int n, double XL, complex<double> Ha, complex<double> mL,
|
|
|
+ complex<double> PsiXL, complex<double> ZetaXL,
|
|
|
+ complex<double> PsiXLM1, complex<double> ZetaXLM1) {
|
|
|
|
|
|
- std::complex<double> Num = (Ha/mL + n/XL)*PsiXL - PsiXLM1;
|
|
|
- std::complex<double> Denom = (Ha/mL + n/XL)*ZetaXL - ZetaXLM1;
|
|
|
+ complex<double> Num = (Ha/mL + n/XL)*PsiXL - PsiXLM1;
|
|
|
+ complex<double> Denom = (Ha/mL + n/XL)*ZetaXL - ZetaXLM1;
|
|
|
|
|
|
return Num/Denom;
|
|
|
}
|
|
@@ -575,12 +577,12 @@ namespace nmie {
|
|
|
// ********************************************************************** //
|
|
|
// Calculate bn - equation (6) //
|
|
|
// ********************************************************************** //
|
|
|
- std::complex<double> MultiLayerMie::calc_bn(int n, double XL, std::complex<double> Hb, std::complex<double> mL,
|
|
|
- std::complex<double> PsiXL, std::complex<double> ZetaXL,
|
|
|
- std::complex<double> PsiXLM1, std::complex<double> ZetaXLM1) {
|
|
|
+ complex<double> MultiLayerMie::calc_bn(int n, double XL, complex<double> Hb, complex<double> mL,
|
|
|
+ complex<double> PsiXL, complex<double> ZetaXL,
|
|
|
+ complex<double> PsiXLM1, complex<double> ZetaXLM1) {
|
|
|
|
|
|
- std::complex<double> Num = (mL*Hb + n/XL)*PsiXL - PsiXLM1;
|
|
|
- std::complex<double> Denom = (mL*Hb + n/XL)*ZetaXL - ZetaXLM1;
|
|
|
+ complex<double> Num = (mL*Hb + n/XL)*PsiXL - PsiXLM1;
|
|
|
+ complex<double> Denom = (mL*Hb + n/XL)*ZetaXL - ZetaXLM1;
|
|
|
|
|
|
return Num/Denom;
|
|
|
}
|
|
@@ -589,8 +591,9 @@ namespace nmie {
|
|
|
// ********************************************************************** //
|
|
|
// Calculates S1 - equation (25a) //
|
|
|
// ********************************************************************** //
|
|
|
- std::complex<double> MultiLayerMie::calc_S1(int n, std::complex<double> an, std::complex<double> bn,
|
|
|
- double Pi, double Tau) {
|
|
|
+ complex<double> MultiLayerMie::calc_S1(int n, complex<double> an, complex<double> bn,
|
|
|
+ double Pi, double Tau) {
|
|
|
+
|
|
|
return double(n + n + 1)*(Pi*an + Tau*bn)/double(n*n + n);
|
|
|
}
|
|
|
|
|
@@ -599,8 +602,9 @@ namespace nmie {
|
|
|
// Calculates S2 - equation (25b) (it's the same as (25a), just switches //
|
|
|
// Pi and Tau) //
|
|
|
// ********************************************************************** //
|
|
|
- std::complex<double> MultiLayerMie::calc_S2(int n, std::complex<double> an, std::complex<double> bn,
|
|
|
- double Pi, double Tau) {
|
|
|
+ complex<double> MultiLayerMie::calc_S2(int n, complex<double> an, complex<double> bn,
|
|
|
+ double Pi, double Tau) {
|
|
|
+
|
|
|
return calc_S1(n, an, bn, Tau, Pi);
|
|
|
}
|
|
|
|
|
@@ -617,31 +621,31 @@ namespace nmie {
|
|
|
// Output parameters: //
|
|
|
// D1, D3: Logarithmic derivatives of the Riccati-Bessel functions //
|
|
|
//**********************************************************************************//
|
|
|
- void MultiLayerMie::calcD1D3(const std::complex<double> z,
|
|
|
- std::vector<std::complex<double> >& D1,
|
|
|
- std::vector<std::complex<double> >& D3) {
|
|
|
+ void MultiLayerMie::calcD1D3(const complex<double> z,
|
|
|
+ vector<complex<double> >& D1,
|
|
|
+ vector<complex<double> >& D3) {
|
|
|
|
|
|
// Downward recurrence for D1 - equations (16a) and (16b)
|
|
|
- D1[nmax_] = std::complex<double>(0.0, 0.0);
|
|
|
- const std::complex<double> zinv = std::complex<double>(1.0, 0.0)/z;
|
|
|
+ D1[nmax_] = complex<double>(0.0, 0.0);
|
|
|
+ const complex<double> zinv = complex<double>(1.0, 0.0)/z;
|
|
|
|
|
|
for (int n = nmax_; n > 0; n--) {
|
|
|
D1[n - 1] = static_cast<double>(n)*zinv - 1.0/(D1[n] + static_cast<double>(n)*zinv);
|
|
|
}
|
|
|
|
|
|
- if (std::abs(D1[0]) > 1.0e15) {
|
|
|
- throw std::invalid_argument("Unstable D1! Please, try to change input parameters!\n");
|
|
|
+ if (abs(D1[0]) > 1.0e15) {
|
|
|
+ throw invalid_argument("Unstable D1! Please, try to change input parameters!\n");
|
|
|
//printf("Warning: Potentially unstable D1! Please, try to change input parameters!\n");
|
|
|
}
|
|
|
|
|
|
// Upward recurrence for PsiZeta and D3 - equations (18a) - (18d)
|
|
|
- PsiZeta_[0] = 0.5*(1.0 - std::complex<double>(std::cos(2.0*z.real()), std::sin(2.0*z.real()))
|
|
|
- *std::exp(-2.0*z.imag()));
|
|
|
- D3[0] = std::complex<double>(0.0, 1.0);
|
|
|
+ PsiZeta_[0] = 0.5*(1.0 - complex<double>(cos(2.0*z.real()), sin(2.0*z.real()))
|
|
|
+ *exp(-2.0*z.imag()));
|
|
|
+ D3[0] = complex<double>(0.0, 1.0);
|
|
|
for (int n = 1; n <= nmax_; n++) {
|
|
|
PsiZeta_[n] = PsiZeta_[n - 1]*(static_cast<double>(n)*zinv - D1[n - 1])
|
|
|
*(static_cast<double>(n)*zinv - D3[n - 1]);
|
|
|
- D3[n] = D1[n] + std::complex<double>(0.0, 1.0)/PsiZeta_[n];
|
|
|
+ D3[n] = D1[n] + complex<double>(0.0, 1.0)/PsiZeta_[n];
|
|
|
}
|
|
|
}
|
|
|
|
|
@@ -658,19 +662,19 @@ namespace nmie {
|
|
|
// Output parameters: //
|
|
|
// Psi, Zeta: Riccati-Bessel functions //
|
|
|
//**********************************************************************************//
|
|
|
- void MultiLayerMie::calcPsiZeta(std::complex<double> z,
|
|
|
- std::vector<std::complex<double> >& Psi,
|
|
|
- std::vector<std::complex<double> >& Zeta) {
|
|
|
+ void MultiLayerMie::calcPsiZeta(complex<double> z,
|
|
|
+ vector<complex<double> >& Psi,
|
|
|
+ vector<complex<double> >& Zeta) {
|
|
|
|
|
|
- std::complex<double> c_i(0.0, 1.0);
|
|
|
- std::vector<std::complex<double> > D1(nmax_ + 1), D3(nmax_ + 1);
|
|
|
+ complex<double> c_i(0.0, 1.0);
|
|
|
+ vector<complex<double> > D1(nmax_ + 1), D3(nmax_ + 1);
|
|
|
|
|
|
// First, calculate the logarithmic derivatives
|
|
|
calcD1D3(z, D1, D3);
|
|
|
|
|
|
// Now, use the upward recurrence to calculate Psi and Zeta - equations (20a) - (21b)
|
|
|
- Psi[0] = std::sin(z);
|
|
|
- Zeta[0] = std::sin(z) - c_i*std::cos(z);
|
|
|
+ Psi[0] = sin(z);
|
|
|
+ Zeta[0] = sin(z) - c_i*cos(z);
|
|
|
for (int n = 1; n <= nmax_; n++) {
|
|
|
Psi[n] = Psi[n - 1]*(static_cast<double>(n)/z - D1[n - 1]);
|
|
|
Zeta[n] = Zeta[n - 1]*(static_cast<double>(n)/z - D3[n - 1]);
|
|
@@ -692,7 +696,7 @@ namespace nmie {
|
|
|
// Pi, Tau: Angular functions Pi and Tau, as defined in equations (26a) - (26c) //
|
|
|
//**********************************************************************************//
|
|
|
void MultiLayerMie::calcPiTau(const double& costheta,
|
|
|
- std::vector<double>& Pi, std::vector<double>& Tau) {
|
|
|
+ vector<double>& Pi, vector<double>& Tau) {
|
|
|
|
|
|
int i;
|
|
|
//****************************************************//
|
|
@@ -729,17 +733,15 @@ namespace nmie {
|
|
|
// Output parameters: //
|
|
|
// Mo1n, Me1n, No1n, Ne1n: Complex vector spherical harmonics //
|
|
|
//**********************************************************************************//
|
|
|
- void MultiLayerMie::calcSpherHarm(const std::complex<double> Rho, const double Theta, const double Phi,
|
|
|
- const std::complex<double>& rn, const std::complex<double>& Dn,
|
|
|
+ void MultiLayerMie::calcSpherHarm(const complex<double> Rho, const double Theta, const double Phi,
|
|
|
+ const complex<double>& rn, const complex<double>& Dn,
|
|
|
const double& Pi, const double& Tau, const double& n,
|
|
|
- std::vector<std::complex<double> >& Mo1n, std::vector<std::complex<double> >& Me1n,
|
|
|
- std::vector<std::complex<double> >& No1n, std::vector<std::complex<double> >& Ne1n) {
|
|
|
+ vector<complex<double> >& Mo1n, vector<complex<double> >& Me1n,
|
|
|
+ vector<complex<double> >& No1n, vector<complex<double> >& Ne1n) {
|
|
|
|
|
|
// using eq 4.50 in BH
|
|
|
- std::complex<double> c_zero(0.0, 0.0);
|
|
|
+ complex<double> c_zero(0.0, 0.0);
|
|
|
|
|
|
- using std::sin;
|
|
|
- using std::cos;
|
|
|
Mo1n[0] = c_zero;
|
|
|
Mo1n[1] = cos(Phi)*Pi*rn/Rho;
|
|
|
Mo1n[2] = -sin(Phi)*Tau*rn/Rho;
|
|
@@ -778,8 +780,8 @@ namespace nmie {
|
|
|
|
|
|
isScaCoeffsCalc_ = false;
|
|
|
|
|
|
- const std::vector<double>& x = size_param_;
|
|
|
- const std::vector<std::complex<double> >& m = refractive_index_;
|
|
|
+ const vector<double>& x = size_param_;
|
|
|
+ const vector<complex<double> >& m = refractive_index_;
|
|
|
const int& pl = PEC_layer_position_;
|
|
|
const int L = refractive_index_.size();
|
|
|
|
|
@@ -792,7 +794,7 @@ namespace nmie {
|
|
|
if (nmax_preset_ <= 0) calcNmax(fl);
|
|
|
else nmax_ = nmax_preset_;
|
|
|
|
|
|
- std::complex<double> z1, z2;
|
|
|
+ complex<double> z1, z2;
|
|
|
//**************************************************************************//
|
|
|
// Note that since Fri, Nov 14, 2014 all arrays start from 0 (zero), which //
|
|
|
// means that index = layer number - 1 or index = n - 1. The only exception //
|
|
@@ -801,10 +803,10 @@ namespace nmie {
|
|
|
// between different arrays. The change was done to optimize memory usage. //
|
|
|
//**************************************************************************//
|
|
|
// Allocate memory to the arrays
|
|
|
- std::vector<std::complex<double> > D1_mlxl(nmax_ + 1), D1_mlxlM1(nmax_ + 1),
|
|
|
+ vector<complex<double> > D1_mlxl(nmax_ + 1), D1_mlxlM1(nmax_ + 1),
|
|
|
D3_mlxl(nmax_ + 1), D3_mlxlM1(nmax_ + 1);
|
|
|
|
|
|
- std::vector<std::vector<std::complex<double> > > Q(L), Ha(L), Hb(L);
|
|
|
+ vector<vector<complex<double> > > Q(L), Ha(L), Hb(L);
|
|
|
|
|
|
for (int l = 0; l < L; l++) {
|
|
|
Q[l].resize(nmax_ + 1);
|
|
@@ -816,15 +818,15 @@ namespace nmie {
|
|
|
bn_.resize(nmax_);
|
|
|
PsiZeta_.resize(nmax_ + 1);
|
|
|
|
|
|
- std::vector<std::complex<double> > PsiXL(nmax_ + 1), ZetaXL(nmax_ + 1);
|
|
|
+ vector<complex<double> > PsiXL(nmax_ + 1), ZetaXL(nmax_ + 1);
|
|
|
|
|
|
//*************************************************//
|
|
|
// Calculate D1 and D3 for z1 in the first layer //
|
|
|
//*************************************************//
|
|
|
if (fl == pl) { // PEC layer
|
|
|
for (int n = 0; n <= nmax_; n++) {
|
|
|
- D1_mlxl[n] = std::complex<double>(0.0, - 1.0);
|
|
|
- D3_mlxl[n] = std::complex<double>(0.0, 1.0);
|
|
|
+ D1_mlxl[n] = complex<double>(0.0, - 1.0);
|
|
|
+ D3_mlxl[n] = complex<double>(0.0, 1.0);
|
|
|
}
|
|
|
} else { // Regular layer
|
|
|
z1 = x[fl]* m[fl];
|
|
@@ -842,8 +844,8 @@ namespace nmie {
|
|
|
//*****************************************************//
|
|
|
// Iteration from the second layer to the last one (L) //
|
|
|
//*****************************************************//
|
|
|
- std::complex<double> Temp, Num, Denom;
|
|
|
- std::complex<double> G1, G2;
|
|
|
+ complex<double> Temp, Num, Denom;
|
|
|
+ complex<double> G1, G2;
|
|
|
for (int l = fl + 1; l < L; l++) {
|
|
|
//************************************************************//
|
|
|
//Calculate D1 and D3 for z1 and z2 in the layers fl + 1..L //
|
|
@@ -859,9 +861,9 @@ namespace nmie {
|
|
|
//Calculate Q, Ha and Hb in the layers fl + 1..L //
|
|
|
//*************************************************//
|
|
|
// Upward recurrence for Q - equations (19a) and (19b)
|
|
|
- Num = std::exp(-2.0*(z1.imag() - z2.imag()))
|
|
|
- *std::complex<double>(std::cos(-2.0*z2.real()) - std::exp(-2.0*z2.imag()), std::sin(-2.0*z2.real()));
|
|
|
- Denom = std::complex<double>(std::cos(-2.0*z1.real()) - std::exp(-2.0*z1.imag()), std::sin(-2.0*z1.real()));
|
|
|
+ Num = exp(-2.0*(z1.imag() - z2.imag()))
|
|
|
+ *complex<double>(cos(-2.0*z2.real()) - exp(-2.0*z2.imag()), sin(-2.0*z2.real()));
|
|
|
+ Denom = complex<double>(cos(-2.0*z1.real()) - exp(-2.0*z1.imag()), sin(-2.0*z1.real()));
|
|
|
Q[l][0] = Num/Denom;
|
|
|
for (int n = 1; n <= nmax_; n++) {
|
|
|
Num = (z1*D1_mlxl[n] + double(n))*(double(n) - z1*D3_mlxl[n - 1]);
|
|
@@ -919,7 +921,7 @@ namespace nmie {
|
|
|
an_[n] = calc_an(n + 1, x[L - 1], Ha[L - 1][n], m[L - 1], PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
|
|
|
bn_[n] = calc_bn(n + 1, x[L - 1], Hb[L - 1][n], m[L - 1], PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
|
|
|
} else {
|
|
|
- an_[n] = calc_an(n + 1, x[L - 1], std::complex<double>(0.0, 0.0), std::complex<double>(1.0, 0.0), PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
|
|
|
+ an_[n] = calc_an(n + 1, x[L - 1], complex<double>(0.0, 0.0), complex<double>(1.0, 0.0), PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
|
|
|
bn_[n] = PsiXL[n + 1]/ZetaXL[n + 1];
|
|
|
}
|
|
|
} // end of for an and bn terms
|
|
@@ -957,11 +959,11 @@ namespace nmie {
|
|
|
//**********************************************************************************//
|
|
|
void MultiLayerMie::RunMieCalculation() {
|
|
|
if (size_param_.size() != refractive_index_.size())
|
|
|
- throw std::invalid_argument("Each size parameter should have only one index!");
|
|
|
+ throw invalid_argument("Each size parameter should have only one index!");
|
|
|
if (size_param_.size() == 0)
|
|
|
- throw std::invalid_argument("Initialize model first!");
|
|
|
+ throw invalid_argument("Initialize model first!");
|
|
|
|
|
|
- const std::vector<double>& x = size_param_;
|
|
|
+ const vector<double>& x = size_param_;
|
|
|
|
|
|
isExpCoeffsCalc_ = false;
|
|
|
isScaCoeffsCalc_ = false;
|
|
@@ -971,7 +973,7 @@ namespace nmie {
|
|
|
calcScattCoeffs();
|
|
|
|
|
|
if (!isScaCoeffsCalc_) // TODO seems to be unreachable
|
|
|
- throw std::invalid_argument("Calculation of scattering coefficients failed!");
|
|
|
+ throw invalid_argument("Calculation of scattering coefficients failed!");
|
|
|
|
|
|
// Initialize the scattering parameters
|
|
|
Qext_ = 0.0;
|
|
@@ -983,14 +985,14 @@ namespace nmie {
|
|
|
albedo_ = 0.0;
|
|
|
|
|
|
// Initialize the scattering amplitudes
|
|
|
- std::vector<std::complex<double> > tmp1(theta_.size(),std::complex<double>(0.0, 0.0));
|
|
|
+ vector<complex<double> > tmp1(theta_.size(),complex<double>(0.0, 0.0));
|
|
|
S1_.swap(tmp1);
|
|
|
S2_ = S1_;
|
|
|
|
|
|
- std::vector<double> Pi(nmax_), Tau(nmax_);
|
|
|
+ vector<double> Pi(nmax_), Tau(nmax_);
|
|
|
|
|
|
- std::complex<double> Qbktmp(0.0, 0.0);
|
|
|
- std::vector< std::complex<double> > Qbktmp_ch(nmax_ - 1, Qbktmp);
|
|
|
+ complex<double> Qbktmp(0.0, 0.0);
|
|
|
+ vector< complex<double> > Qbktmp_ch(nmax_ - 1, Qbktmp);
|
|
|
// By using downward recurrence we avoid loss of precision due to float rounding errors
|
|
|
// See: https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html
|
|
|
// http://en.wikipedia.org/wiki/Loss_of_significance
|
|
@@ -1002,14 +1004,14 @@ namespace nmie {
|
|
|
Qsca_ += (n + n + 1.0)*(an_[i].real()*an_[i].real() + an_[i].imag()*an_[i].imag()
|
|
|
+ bn_[i].real()*bn_[i].real() + bn_[i].imag()*bn_[i].imag());
|
|
|
// Equation (29)
|
|
|
- Qpr_ += ((n*(n + 2)/(n + 1))*((an_[i]*std::conj(an_[n]) + bn_[i]*std::conj(bn_[n])).real())
|
|
|
- + ((double)(n + n + 1)/(n*(n + 1)))*(an_[i]*std::conj(bn_[i])).real());
|
|
|
+ Qpr_ += ((n*(n + 2)/(n + 1))*((an_[i]*conj(an_[n]) + bn_[i]*conj(bn_[n])).real())
|
|
|
+ + ((double)(n + n + 1)/(n*(n + 1)))*(an_[i]*conj(bn_[i])).real());
|
|
|
// Equation (33)
|
|
|
Qbktmp += (double)(n + n + 1)*(1 - 2*(n % 2))*(an_[i]- bn_[i]);
|
|
|
// Calculate the scattering amplitudes (S1 and S2) //
|
|
|
// Equations (25a) - (25b) //
|
|
|
for (unsigned int t = 0; t < theta_.size(); t++) {
|
|
|
- calcPiTau(std::cos(theta_[t]), Pi, Tau);
|
|
|
+ calcPiTau(cos(theta_[t]), Pi, Tau);
|
|
|
|
|
|
S1_[t] += calc_S1(n, an_[i], bn_[i], Pi[i], Tau[i]);
|
|
|
S2_[t] += calc_S2(n, an_[i], bn_[i], Pi[i], Tau[i]);
|
|
@@ -1049,11 +1051,11 @@ namespace nmie {
|
|
|
//**********************************************************************************//
|
|
|
void MultiLayerMie::calcExpanCoeffs() {
|
|
|
if (!isScaCoeffsCalc_)
|
|
|
- throw std::invalid_argument("(ExpanCoeffs) You should calculate external coefficients first!");
|
|
|
+ throw invalid_argument("(ExpanCoeffs) You should calculate external coefficients first!");
|
|
|
|
|
|
isExpCoeffsCalc_ = false;
|
|
|
|
|
|
- std::complex<double> c_one(1.0, 0.0), c_zero(0.0, 0.0);
|
|
|
+ complex<double> c_one(1.0, 0.0), c_zero(0.0, 0.0);
|
|
|
|
|
|
const int L = refractive_index_.size();
|
|
|
|
|
@@ -1077,17 +1079,17 @@ namespace nmie {
|
|
|
dln_[L][n] = c_one;
|
|
|
}
|
|
|
|
|
|
- std::vector<std::complex<double> > D1z(nmax_ + 1), D1z1(nmax_ + 1), D3z(nmax_ + 1), D3z1(nmax_ + 1);
|
|
|
- std::vector<std::complex<double> > Psiz(nmax_ + 1), Psiz1(nmax_ + 1), Zetaz(nmax_ + 1), Zetaz1(nmax_ + 1);
|
|
|
- std::complex<double> denomZeta, denomPsi, T1, T2, T3, T4;
|
|
|
+ vector<complex<double> > D1z(nmax_ + 1), D1z1(nmax_ + 1), D3z(nmax_ + 1), D3z1(nmax_ + 1);
|
|
|
+ vector<complex<double> > Psiz(nmax_ + 1), Psiz1(nmax_ + 1), Zetaz(nmax_ + 1), Zetaz1(nmax_ + 1);
|
|
|
+ complex<double> denomZeta, denomPsi, T1, T2, T3, T4;
|
|
|
|
|
|
auto& m = refractive_index_;
|
|
|
- std::vector< std::complex<double> > m1(L);
|
|
|
+ vector< complex<double> > m1(L);
|
|
|
|
|
|
for (int l = 0; l < L - 1; l++) m1[l] = m[l + 1];
|
|
|
- m1[L - 1] = std::complex<double> (1.0, 0.0);
|
|
|
+ m1[L - 1] = complex<double> (1.0, 0.0);
|
|
|
|
|
|
- std::complex<double> z, z1;
|
|
|
+ complex<double> z, z1;
|
|
|
for (int l = L - 1; l >= 0; l--) {
|
|
|
z = size_param_[l]*m[l];
|
|
|
z1 = size_param_[l]*m1[l];
|
|
@@ -1122,15 +1124,15 @@ namespace nmie {
|
|
|
|
|
|
// Check the result and change aln_[0][n] and aln_[0][n] for exact zero
|
|
|
for (int n = 0; n < nmax_; ++n) {
|
|
|
- if (std::abs(aln_[0][n]) < 1e-10) aln_[0][n] = 0.0;
|
|
|
+ if (abs(aln_[0][n]) < 1e-10) aln_[0][n] = 0.0;
|
|
|
else {
|
|
|
- //throw std::invalid_argument("Unstable calculation of aln_[0][n]!");
|
|
|
+ //throw invalid_argument("Unstable calculation of aln_[0][n]!");
|
|
|
printf("Warning: Potentially unstable calculation of aln (aln[0][%i] = %g, %gi)\n", n, aln_[0][n].real(), aln_[0][n].imag());
|
|
|
aln_[0][n] = 0.0;
|
|
|
}
|
|
|
- if (std::abs(bln_[0][n]) < 1e-10) bln_[0][n] = 0.0;
|
|
|
+ if (abs(bln_[0][n]) < 1e-10) bln_[0][n] = 0.0;
|
|
|
else {
|
|
|
- //throw std::invalid_argument("Unstable calculation of bln_[0][n]!");
|
|
|
+ //throw invalid_argument("Unstable calculation of bln_[0][n]!");
|
|
|
printf("Warning: Potentially unstable calculation of bln (bln[0][%i] = %g, %gi)\n", n, bln_[0][n].real(), bln_[0][n].imag());
|
|
|
bln_[0][n] = 0.0;
|
|
|
}
|
|
@@ -1153,17 +1155,17 @@ namespace nmie {
|
|
|
// E, H: Complex electric and magnetic fields //
|
|
|
//**********************************************************************************//
|
|
|
void MultiLayerMie::calcField(const double Rho, const double Theta, const double Phi,
|
|
|
- std::vector<std::complex<double> >& E, std::vector<std::complex<double> >& H) {
|
|
|
+ vector<complex<double> >& E, vector<complex<double> >& H) {
|
|
|
|
|
|
- std::complex<double> c_zero(0.0, 0.0), c_i(0.0, 1.0), c_one(1.0, 0.0);
|
|
|
- std::vector<std::complex<double> > ipow = {c_one, c_i, -c_one, -c_i}; // Vector containing precomputed integer powers of i to avoid computation
|
|
|
- std::vector<std::complex<double> > M3o1n(3), M3e1n(3), N3o1n(3), N3e1n(3);
|
|
|
- std::vector<std::complex<double> > M1o1n(3), M1e1n(3), N1o1n(3), N1e1n(3);
|
|
|
- std::vector<std::complex<double> > Psi(nmax_ + 1), D1n(nmax_ + 1), Zeta(nmax_ + 1), D3n(nmax_ + 1);
|
|
|
- std::vector<double> Pi(nmax_), Tau(nmax_);
|
|
|
+ complex<double> c_zero(0.0, 0.0), c_i(0.0, 1.0), c_one(1.0, 0.0);
|
|
|
+ vector<complex<double> > ipow = {c_one, c_i, -c_one, -c_i}; // Vector containing precomputed integer powers of i to avoid computation
|
|
|
+ vector<complex<double> > M3o1n(3), M3e1n(3), N3o1n(3), N3e1n(3);
|
|
|
+ vector<complex<double> > M1o1n(3), M1e1n(3), N1o1n(3), N1e1n(3);
|
|
|
+ vector<complex<double> > Psi(nmax_ + 1), D1n(nmax_ + 1), Zeta(nmax_ + 1), D3n(nmax_ + 1);
|
|
|
+ vector<double> Pi(nmax_), Tau(nmax_);
|
|
|
|
|
|
int l = 0; // Layer number
|
|
|
- std::complex<double> ml;
|
|
|
+ complex<double> ml;
|
|
|
|
|
|
// Initialize E and H
|
|
|
for (int i = 0; i < 3; i++) {
|
|
@@ -1189,7 +1191,7 @@ namespace nmie {
|
|
|
calcPsiZeta(Rho*ml, Psi, Zeta);
|
|
|
|
|
|
// Calculate angular functions Pi and Tau
|
|
|
- calcPiTau(std::cos(Theta), Pi, Tau);
|
|
|
+ calcPiTau(cos(Theta), Pi, Tau);
|
|
|
|
|
|
for (int n = nmax_ - 2; n >= 0; n--) {
|
|
|
int n1 = n + 1;
|
|
@@ -1200,7 +1202,7 @@ namespace nmie {
|
|
|
calcSpherHarm(Rho*ml, Theta, Phi, Zeta[n1], D3n[n1], Pi[n], Tau[n], rn, M3o1n, M3e1n, N3o1n, N3e1n);
|
|
|
|
|
|
// Total field in the lth layer: eqs. (1) and (2) in Yang, Appl. Opt., 42 (2003) 1710-1720
|
|
|
- std::complex<double> En = ipow[n1 % 4]*(rn + rn + 1.0)/(rn*rn + rn);
|
|
|
+ complex<double> En = ipow[n1 % 4]*(rn + rn + 1.0)/(rn*rn + rn);
|
|
|
for (int i = 0; i < 3; i++) {
|
|
|
// electric field E [V m - 1] = EF*E0
|
|
|
E[i] += En*(cln_[l][n]*M1o1n[i] - c_i*dln_[l][n]*N1e1n[i]
|
|
@@ -1212,7 +1214,7 @@ namespace nmie {
|
|
|
} // end of for all n
|
|
|
|
|
|
// magnetic field
|
|
|
- std::complex<double> hffact = ml/(cc_*mu_);
|
|
|
+ complex<double> hffact = ml/(cc_*mu_);
|
|
|
for (int i = 0; i < 3; i++) {
|
|
|
H[i] = hffact*H[i];
|
|
|
}
|
|
@@ -1262,16 +1264,16 @@ namespace nmie {
|
|
|
const double& Zp = coords_[2][point];
|
|
|
|
|
|
// Convert to spherical coordinates
|
|
|
- Rho = std::sqrt(pow2(Xp) + pow2(Yp) + pow2(Zp));
|
|
|
+ Rho = sqrt(pow2(Xp) + pow2(Yp) + pow2(Zp));
|
|
|
|
|
|
// If Rho=0 then Theta is undefined. Just set it to zero to avoid problems
|
|
|
- Theta = (Rho > 0.0) ? std::acos(Zp/Rho) : 0.0;
|
|
|
+ Theta = (Rho > 0.0) ? acos(Zp/Rho) : 0.0;
|
|
|
|
|
|
// If Xp=Yp=0 then Phi is undefined. Just set it to zero to avoid problems
|
|
|
if (Xp == 0.0)
|
|
|
- Phi = (Yp != 0.0) ? std::asin(Yp/std::sqrt(pow2(Xp) + pow2(Yp))) : 0.0;
|
|
|
+ Phi = (Yp != 0.0) ? asin(Yp/sqrt(pow2(Xp) + pow2(Yp))) : 0.0;
|
|
|
else
|
|
|
- Phi = std::acos(Xp/std::sqrt(pow2(Xp) + pow2(Yp)));
|
|
|
+ Phi = acos(Xp/sqrt(pow2(Xp) + pow2(Yp)));
|
|
|
|
|
|
// Avoid convergence problems due to Rho too small
|
|
|
if (Rho < 1e-5) Rho = 1e-5;
|
|
@@ -1283,14 +1285,12 @@ namespace nmie {
|
|
|
//*******************************************************//
|
|
|
|
|
|
// This array contains the fields in spherical coordinates
|
|
|
- std::vector<std::complex<double> > Es(3), Hs(3);
|
|
|
+ vector<complex<double> > Es(3), Hs(3);
|
|
|
|
|
|
// Do the actual calculation of electric and magnetic field
|
|
|
calcField(Rho, Theta, Phi, Es, Hs);
|
|
|
|
|
|
{ //Now, convert the fields back to cartesian coordinates
|
|
|
- using std::sin;
|
|
|
- using std::cos;
|
|
|
E_[point][0] = sin(Theta)*cos(Phi)*Es[0] + cos(Theta)*cos(Phi)*Es[1] - sin(Phi)*Es[2];
|
|
|
E_[point][1] = sin(Theta)*sin(Phi)*Es[0] + cos(Theta)*sin(Phi)*Es[1] + cos(Phi)*Es[2];
|
|
|
E_[point][2] = cos(Theta)*Es[0] - sin(Theta)*Es[1];
|