#ifndef SRC_NMIE_NMIE_WRAPPER_H_ #define SRC_NMIE_NMIE_WRAPPER_H_ /// /// @file nmie-wrapper.h /// @author Ladutenko Konstantin /// @date Tue Sep 3 00:40:47 2013 /// @copyright 2013 Ladutenko Konstantin /// /// nmie-wrapper 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. /// /// nmie-wrapper 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. /// /// You should have received a copy of the GNU General Public License /// along with nmie-wrapper. If not, see . /// /// nmie-wrapper uses nmie.c from scattnlay by Ovidio Pena /// as a linked library. He has an additional condition to /// his library: // The only additional condition 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. // /// /// @brief Wrapper class around nMie function for ease of use /// /// #include #include #include #include #include #ifndef NDEBUG # define ASSERT(condition, message) \ do { \ if (! (condition)) { \ std::cerr << "Assertion `" #condition "` failed in " << __FILE__ \ << " line " << __LINE__ << ": " << message << std::endl; \ std::exit(EXIT_FAILURE); \ } \ } while (false) #else # define ASSERT(condition, message) do { } while (false) #endif namespace nmie { int nMie_wrapper(int L, const std::vector& x, const std::vector >& m, int nTheta, const std::vector& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector >& S1, std::vector >& S2); class MultiLayerMie { // Will throw for any error! // SP stands for size parameter units. public: // Set parameters in applied units void SetWavelength(double wavelength) {wavelength_ = wavelength;}; // It is possible to set only a multilayer target to run calculaitons. // For many runs it can be convenient to separate target and coating layers. // Per layer void AddTargetLayer(double layer_width, std::complex layer_index); void AddCoatingLayer(double layer_width, std::complex layer_index); // For all layers void SetTargetWidth(std::vector width); void SetTargetIndex(std::vector< std::complex > index); void SetCoatingWidth(std::vector width); void SetCoatingIndex(std::vector< std::complex > index); void SetFieldPoints(std::vector< std::array > coords); //Set parameters in size parameter units void SetWidthSP(const std::vector& width); void SetIndexSP(const std::vector< std::complex >& index); void SetFieldPointsSP(std::vector< std::array > coords); // Set common parameters void SetAnglesForPattern(double from_angle, double to_angle, int samples); void SetAngles(const std::vector& angles); std::vector GetAngles(); void SetPEC(int layer_position = 0); // By default set PEC layer to be the first one void SetMaxTermsNumber(int nmax); void ClearTarget(); void ClearCoating(); void ClearLayers(); // Applied units requests double GetTotalRadius(); double GetTargetRadius(); double GetCoatingWidth(); std::vector GetTargetLayersWidth(); std::vector< std::complex > GetTargetLayersIndex(); std::vector GetCoatingLayersWidth(); std::vector< std::complex > GetCoatingLayersIndex(); std::vector< std::array > GetFieldPoints(); std::vector,3 > > GetFieldE(); std::vector,3 > > GetFieldH(); std::vector< std::array > GetSpectra(double from_WL, double to_WL, int samples); // ext, sca, abs, bk double GetRCSext(); double GetRCSsca(); double GetRCSabs(); double GetRCSbk(); std::vector GetPatternEk(); std::vector GetPatternHk(); std::vector GetPatternUnpolarized(); // Size parameter units std::vector GetLayerWidthSP(); // Same as to get target and coating index std::vector< std::complex > GetLayerIndex(); std::vector< std::array > GetFieldPointsSP(); // Do we need normalize field to size parameter? /* std::vector > > GetFieldESP(); */ /* std::vector > > GetFieldHSP(); */ std::vector< std::array > GetSpectraSP(double from_SP, double to_SP, int samples); // WL,ext, sca, abs, bk double GetQext(); double GetQsca(); double GetQabs(); double GetQbk(); double GetQpr(); double GetAsymmetryFactor(); double GetAlbedo(); std::vector > GetS1(); std::vector > GetS2(); std::vector GetPatternEkSP(); std::vector GetPatternHkSP(); std::vector GetPatternUnpolarizedSP(); // Run calculation void RunMieCalculations(); void RunFieldCalculations(); // Output results (data file + python script to plot it with matplotlib) void PlotSpectra(); void PlotSpectraSP(); void PlotField(); void PlotFieldSP(); void PlotPattern(); void PlotPatternSP(); private: void GenerateSizeParameter(); void GenerateIndex(); void InitMieCalculations(); 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(double x, std::vector > D1, std::vector > D3, std::vector >& Psi, std::vector >& Zeta); void calcD1D3(std::complex z, std::vector >& D1, std::vector >& D3); void calcPiTau( std::vector< std::vector >& Pi, std::vector< std::vector >& Tau); void ScattCoeffs(std::vector >& an, std::vector >& bn); void fieldExt( double Rho, double Phi, double Theta, std::vector Pi, std::vector Tau, std::vector > an, std::vector > bn, std::vector >& E, std::vector >& H); bool isMieCalculated_ = false; double wavelength_ = 1.0; double total_radius_ = 0.0; /// Width and index for each layer of the structure std::vector target_width_, coating_width_; std::vector< std::complex > target_index_, coating_index_; /// Size parameters for all layers std::vector size_parameter_; /// Complex index values for each layers. std::vector< std::complex > index_; /// Scattering angles for RCS pattern in radians std::vector theta_; // Should be -1 if there is no PEC. int PEC_layer_position_ = -1; // Set nmax_ manualy with SetMaxTermsNumber(int nmax) or in ScattCoeffs(..) // with Nmax(int first_layer); int nmax_ = -1; /// 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 > 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 #endif // SRC_NMIE_NMIE_WRAPPER_H_