#ifndef SRC_NMIE_NMIE_H_ #define SRC_NMIE_NMIE_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); 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 { // Will throw for any error! // SP stands for size parameter units. public: void GetFailed(); long iformat = 0; bool output = true; void prn(double var) { do { if (!output) break; ++iformat; printf("%23.13e",var); if (iformat%4 == 0) printf("\n"); } while (false); } // 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 SetTargetPEC(double radius); 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(const std::vector< std::vector >& coords_sp); // 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); int GetMaxTermsUsed() {return nmax_used_;}; void ClearTarget(); void ClearCoating(); void ClearLayers(); void ClearAllDesign(); //Layers + SP + index_ // 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::vector > GetFieldPoints(); std::vector > > GetFieldE(); // {X[], Y[], Z[]} std::vector > > GetFieldH(); std::vector< std::vector > 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(); std::vector GetQsca_channel(); std::vector GetQabs_channel(); std::vector GetQsca_channel_normalized(); std::vector GetQabs_channel_normalized(); std::vector > GetAn(){return an_;}; std::vector > GetBn(){return bn_;}; 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 ConvertToSP(); 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); 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 ScattCoeffs(std::vector >& an, std::vector >& bn); void ScattCoeffsLayerd(); void ScattCoeffsLayerdInit(); 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); 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; 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 #endif // SRC_NMIE_NMIE_H_