//**********************************************************************************// // Copyright (C) 2009-2015 Ovidio Pena // // // // 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 nMie(int L, std::vector& x, std::vector >& m, 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 nField(const int L, const int pl, const std::vector& x, const std::vector >& m, const int nmax, const 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 RunMieCalculations(); void RunFieldCalculations(); // 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_width, std::complex layer_index); // Modify width of the layer void SetLayerWidth(std::vector layer_width, int layer_position = 0); // Modify refractive index of the layer void SetLayerIndex(std::vector< std::complex > layer_index, int layer_position = 0); // Modify width of all layers void SetLayersWidth(std::vector layer_width); // Modify refractive index of all layers void SetLayerIndex(std::vector< std::complex > layer_index); // Set PEC layer void SetPECLayer(int layer_position = 0); // Set maximun number of terms to be used void SetMaxTerms(int nmax); // Get maximun number of terms int GetMaxTermsUsed() {return nmax_used_;}; // Clear layer information void ClearLayers(); // Applied units requests double GetTotalRadius(); double GetCoreRadius(); double GetLayerWidth(int layer_position = 0); std::vector GetLayersWidth(); std::vector > GetLayersIndex(); std::vector > GetFieldCoords(); std::vector > > GetFieldE(){return E_field_;}; // {X[], Y[], Z[]} std::vector > > GetFieldH(){return H_field_;}; private: void Nstop(); void Nmax(int first_layer); void sbesjh(std::complex z, std::vector >& jn, std::vector >& jnp, std::vector >& h1n, std::vector >& h1np); void sphericalBessel(std::complex z, std::vector >& bj, std::vector >& by, std::vector >& bd); 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); 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 calcPsiZeta(std::complex x, std::vector > D1, std::vector > D3, std::vector >& Psi, std::vector >& Zeta); std::complex calcD1confra(int N, const std::complex z); void calcD1D3(std::complex z, std::vector >& D1, std::vector >& D3); void calcSinglePiTau(const double& costheta, std::vector& Pi, std::vector& Tau); void calcAllPiTau(std::vector< std::vector >& Pi, std::vector< std::vector >& Tau); void ExtScattCoeffs(std::vector >& an, std::vector >& bn); void IntScattCoeffs(); void IntScattCoeffsInit(); void fieldExt(const double Rho, const double Phi, const double Theta, const std::vector& Pi, const std::vector& Tau, std::vector >& E, std::vector >& H); void fieldInt(const double Rho, const double Phi, const double Theta, const std::vector& Pi, const std::vector& Tau, std::vector >& E, std::vector >& H); bool areIntCoeffsCalc_ = false; bool areExtCoeffsCalc_ = false; bool isMieCalculated_ = false; double wavelength_ = 1.0; double total_radius_ = 0.0; // Size parameter for all layers std::vector width_; // Refractive index for all layers std::vector< std::complex > 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 Nmax(int first_layer); int nmax_ = -1; int nmax_used_ = -1; int nmax_preset_ = -1; // Scattering coefficients std::vector > an_, bn_; std::vector< std::vector > coords_sp_; // TODO: check if l index is reversed will lead to performance // boost, if $a^(L+1)_n$ stored in al_n_[n][0], $a^(L)_n$ in // al_n_[n][1] and so on... // at the moment order is forward! std::vector< std::vector > > al_n_, bl_n_, cl_n_, dl_n_; /// 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_field_, H_field_; // {X[], Y[], Z[]} // Mie efficinecy from each multipole channel. std::vector Qsca_ch_, Qext_ch_, Qabs_ch_, Qbk_ch_, Qpr_ch_; std::vector Qsca_ch_norm_, Qext_ch_norm_, Qabs_ch_norm_, Qbk_ch_norm_, Qpr_ch_norm_; 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