//**********************************************************************************// // Copyright (C) 2009-2015 Ovidio Pena // // Copyright (C) 2013-2015 Konstantin Ladutenko // // // // This file is part of scattnlay // // // // This program is free software: you can redistribute it and/or modify // // it under the terms of the GNU General Public License as published by // // the Free Software Foundation, either version 3 of the License, or // // (at your option) any later version. // // // // This program is distributed in the hope that it will be useful, // // but WITHOUT ANY WARRANTY; without even the implied warranty of // // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // // GNU General Public License for more details. // // // // The only additional remark is that we expect that all publications // // describing work using this software, or all commercial products // // using it, cite the following reference: // // [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by // // a multilayered sphere," Computer Physics Communications, // // vol. 180, Nov. 2009, pp. 2348-2354. // // // // You should have received a copy of the GNU General Public License // // along with this program. If not, see . // //**********************************************************************************// #define VERSION "0.3.1" #include #include #include #include #include namespace nmie { int ScattCoeffs(const unsigned int L, const int pl, std::vector& x, std::vector >& m, const int nmax, std::vector > &an, std::vector > &bn); int nMie(const unsigned int L, const int pl, std::vector& x, std::vector >& m, const unsigned int nTheta, std::vector& Theta, const int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector >& S1, std::vector >& S2); int nMie(const unsigned int L, std::vector& x, std::vector >& m, const unsigned int nTheta, std::vector& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector >& S1, std::vector >& S2); int nMie(const unsigned int L, const int pl, std::vector& x, std::vector >& m, const unsigned int nTheta, std::vector& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector >& S1, std::vector >& S2); int nMie(const unsigned int L, std::vector& x, std::vector >& m, const unsigned int nTheta, std::vector& Theta, const int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector >& S1, std::vector >& S2); int nField(const unsigned int L, const int pl, const std::vector& x, const std::vector >& m, const int nmax, const unsigned int ncoord, const std::vector& Xp, const std::vector& Yp, const std::vector& Zp, std::vector > >& E, std::vector > >& H); class MultiLayerMie { public: // Run calculation void RunMieCalculation(); void RunFieldCalculation(); // Return calculation results double GetQext(); double GetQsca(); double GetQabs(); double GetQbk(); double GetQpr(); double GetAsymmetryFactor(); double GetAlbedo(); std::vector > GetS1(); std::vector > GetS2(); std::vector > GetAn(){return an_;}; std::vector > GetBn(){return bn_;}; // Problem definition // Add new layer void AddNewLayer(double layer_size, std::complex layer_index); // Modify width of the layer void SetLayerSize(std::vector layer_size, int layer_position = 0); // Modify refractive index of the layer void SetLayerIndex(std::vector< std::complex > layer_index, int layer_position = 0); // Modify size of all layers void SetLayersSize(const std::vector& layer_size); // Modify refractive index of all layers void SetLayersIndex(const std::vector< std::complex >& index); // Modify scattering (theta) angles void SetAngles(const std::vector& angles); // Modify coordinates for field calculation void SetFieldCoords(const std::vector< std::vector >& coords); // Modify PEC layer void SetPECLayer(int layer_position = 0); // Set a fixed value for the maximun number of terms void SetMaxTerms(int nmax); // Get maximun number of terms int GetMaxTerms() {return nmax_;}; // Clear layer information void ClearLayers(); // Applied units requests double GetSizeParameter(); double GetLayerWidth(int layer_position = 0); std::vector GetLayersSize(); std::vector > GetLayersIndex(); std::vector > GetFieldCoords(); std::vector > > GetFieldE(){return E_;}; // {X[], Y[], Z[]} std::vector > > GetFieldH(){return H_;}; private: void calcNstop(); void calcNmax(unsigned int first_layer); std::complex calc_an(int n, double XL, std::complex Ha, std::complex mL, std::complex PsiXL, std::complex ZetaXL, std::complex PsiXLM1, std::complex ZetaXLM1); std::complex calc_bn(int n, double XL, std::complex Hb, std::complex mL, std::complex PsiXL, std::complex ZetaXL, std::complex PsiXLM1, std::complex ZetaXLM1); void calc_an_bn_bulk(std::vector >& an, std::vector >& bn, double x, std::complex m); std::complex calc_S1(int n, std::complex an, std::complex bn, double Pi, double Tau); std::complex calc_S2(int n, std::complex an, std::complex bn, double Pi, double Tau); void calcD1D3(std::complex z, std::vector >& D1, std::vector >& D3); void calcPsiZeta(std::complex x, std::vector >& Psi, std::vector >& Zeta); void sbesjh(std::complex z, std::vector >& jn, std::vector >& jnp, std::vector >& h1n, std::vector >& h1np); void calcPiTau(const double& costheta, std::vector& Pi, std::vector& Tau); void calcSpherHarm(const double Rho, const double Theta, const double Phi, const std::complex& rn, const std::complex& Dn, const double& Pi, const double& Tau, const double& n, std::vector >& Mo1n, std::vector >& Me1n, std::vector >& No1n, std::vector >& Ne1n); void ScattCoeffs(); void ExpanCoeffs(); void ExpanCoeffsV2(); void fieldExt(const double Rho, const double Theta, const double Phi, std::vector >& E, std::vector >& H); void calcField(const double Rho, const double Theta, const double Phi, std::vector >& E, std::vector >& H); bool isExpCoeffsCalc_ = false; bool isScaCoeffsCalc_ = false; bool isMieCalculated_ = false; // Size parameter for all layers std::vector size_param_; // Refractive index for all layers std::vector< std::complex > refractive_index_; // Scattering angles for scattering pattern in radians std::vector theta_; // Should be -1 if there is no PEC. int PEC_layer_position_ = -1; // with calcNmax(int first_layer); int nmax_ = -1; int nmax_preset_ = -1; // Scattering coefficients std::vector > an_, bn_; std::vector< std::vector > coords_; // TODO: check if l index is reversed will lead to performance // boost, if $a^(L+1)_n$ stored in aln_[n][0], $a^(L)_n$ in // aln_[n][1] and so on... // at the moment order is forward! std::vector< std::vector > > aln_, bln_, cln_, dln_; /// Store result double Qsca_ = 0.0, Qext_ = 0.0, Qabs_ = 0.0, Qbk_ = 0.0, Qpr_ = 0.0, asymmetry_factor_ = 0.0, albedo_ = 0.0; std::vector > > E_, H_; // {X[], Y[], Z[]} std::vector > S1_, S2_; //Used constants const double PI_=3.14159265358979323846; // light speed [m s-1] double const cc_ = 2.99792458e8; // assume non-magnetic (MU=MU0=const) [N A-2] double const mu_ = 4.0*PI_*1.0e-7; //Temporary variables std::vector > PsiZeta_; }; // end of class MultiLayerMie } // end of namespace nmie