|
|
@@ -1,897 +1,954 @@
|
|
|
#ifndef SRC_NMIE_BASIC_HPP_
|
|
|
#define SRC_NMIE_BASIC_HPP_
|
|
|
-//**********************************************************************************//
|
|
|
-// Copyright (C) 2009-2018 Ovidio Pena <ovidio@bytesfall.com> //
|
|
|
-// Copyright (C) 2013-2018 Konstantin Ladutenko <kostyfisik@gmail.com> //
|
|
|
-// //
|
|
|
-// 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 at least one of the following references: //
|
|
|
-// [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by //
|
|
|
-// a multilayered sphere," Computer Physics Communications, //
|
|
|
-// vol. 180, Nov. 2009, pp. 2348-2354. //
|
|
|
-// [2] K. Ladutenko, U. Pal, A. Rivera, and O. Pena-Rodriguez, "Mie //
|
|
|
-// calculation of electromagnetic near-field for a multilayered //
|
|
|
-// sphere," Computer Physics Communications, vol. 214, May 2017, //
|
|
|
-// pp. 225-230. //
|
|
|
-// //
|
|
|
-// You should have received a copy of the GNU General Public License //
|
|
|
-// along with this program. If not, see <http://www.gnu.org/licenses/>. //
|
|
|
-//**********************************************************************************//
|
|
|
-
|
|
|
-//**********************************************************************************//
|
|
|
-// This class implements the algorithm for a multilayered sphere described by: //
|
|
|
-// [1] W. Yang, "Improved recursive algorithm for light scattering by a //
|
|
|
-// multilayered sphere,” Applied Optics, vol. 42, Mar. 2003, pp. 1710-1720. //
|
|
|
-// //
|
|
|
-// You can find the description of all the used equations in: //
|
|
|
-// [2] O. Pena and U. Pal, "Scattering of electromagnetic radiation by //
|
|
|
-// a multilayered sphere," Computer Physics Communications, //
|
|
|
-// vol. 180, Nov. 2009, pp. 2348-2354. //
|
|
|
-// [3] K. Ladutenko, U. Pal, A. Rivera, and O. Pena-Rodriguez, "Mie //
|
|
|
-// calculation of electromagnetic near-field for a multilayered //
|
|
|
-// sphere," Computer Physics Communications, vol. 214, May 2017, //
|
|
|
-// pp. 225-230. //
|
|
|
-// //
|
|
|
-// Hereinafter all equations numbers refer to [2] //
|
|
|
-//**********************************************************************************//
|
|
|
-#include <iostream>
|
|
|
+//***************************************************************************//
|
|
|
+// Copyright (C) 2009-2022 Ovidio Pena <ovidio@bytesfall.com> //
|
|
|
+// Copyright (C) 2013-202 Konstantin Ladutenko <kostyfisik@gmail.com> //
|
|
|
+// //
|
|
|
+// 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 at least one of the following references: //
|
|
|
+// [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by //
|
|
|
+// a multilayered sphere," Computer Physics Communications, //
|
|
|
+// vol. 180, Nov. 2009, pp. 2348-2354. //
|
|
|
+// [2] K. Ladutenko, U. Pal, A. Rivera, and O. Pena-Rodriguez, "Mie //
|
|
|
+// calculation of electromagnetic near-field for a multilayered //
|
|
|
+// sphere," Computer Physics Communications, vol. 214, May 2017, //
|
|
|
+// pp. 225-230. //
|
|
|
+// //
|
|
|
+// You should have received a copy of the GNU General Public License //
|
|
|
+// along with this program. If not, see <http://www.gnu.org/licenses/>. //
|
|
|
+//***************************************************************************//
|
|
|
+//***************************************************************************//
|
|
|
+// This class implements the algorithm for a multilayered sphere described
|
|
|
+// by:
|
|
|
+// [1] W. Yang, "Improved recursive algorithm for light scattering by a
|
|
|
+// multilayered sphere,” Applied Optics, vol. 42, Mar. 2003, pp.
|
|
|
+// 1710-1720.
|
|
|
+//
|
|
|
+// You can find the description of all the used equations in:
|
|
|
+// [2] O. Pena and U. Pal, "Scattering of electromagnetic radiation by
|
|
|
+// a multilayered sphere," Computer Physics Communications,
|
|
|
+// vol. 180, Nov. 2009, pp. 2348-2354.
|
|
|
+// [3] K. Ladutenko, U. Pal, A. Rivera, and O. Pena-Rodriguez, "Mie
|
|
|
+// calculation of electromagnetic near-field for a multilayered
|
|
|
+// sphere," Computer Physics Communications, vol. 214, May 2017,
|
|
|
+// pp. 225-230.
|
|
|
+//
|
|
|
+// Hereinafter all equations numbers refer to [2]
|
|
|
+//*****************************************************************************//
|
|
|
#include <iomanip>
|
|
|
+#include <iostream>
|
|
|
#include <stdexcept>
|
|
|
#include <vector>
|
|
|
|
|
|
-#include "special-functions-impl.hpp"
|
|
|
#include "nmie.hpp"
|
|
|
+#include "special-functions-impl.hpp"
|
|
|
|
|
|
namespace nmie {
|
|
|
|
|
|
+// class implementation
|
|
|
|
|
|
- //class implementation
|
|
|
-
|
|
|
- // ********************************************************************** //
|
|
|
- // Returns previously calculated Qext //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- template <typename outputType>
|
|
|
- outputType MultiLayerMie<FloatType>::GetQext() {
|
|
|
- if (!isMieCalculated_)
|
|
|
- throw std::invalid_argument("You should run calculations before result request!");
|
|
|
- return static_cast<outputType>(Qext_);
|
|
|
- }
|
|
|
-
|
|
|
- // ********************************************************************** //
|
|
|
- // Returns previously calculated Qabs //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- template <typename outputType>
|
|
|
- outputType MultiLayerMie<FloatType>::GetQabs() {
|
|
|
- if (!isMieCalculated_)
|
|
|
- throw std::invalid_argument("You should run calculations before result request!");
|
|
|
- return static_cast<outputType>(Qabs_);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- // ********************************************************************** //
|
|
|
- // Returns previously calculated Qsca //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- template <typename outputType>
|
|
|
- outputType MultiLayerMie<FloatType>::GetQsca() {
|
|
|
- if (!isMieCalculated_)
|
|
|
- throw std::invalid_argument("You should run calculations before result request!");
|
|
|
- return static_cast<outputType>(Qsca_);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- // ********************************************************************** //
|
|
|
- // Returns previously calculated Qbk //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- template <typename outputType>
|
|
|
- outputType MultiLayerMie<FloatType>::GetQbk() {
|
|
|
- if (!isMieCalculated_)
|
|
|
- throw std::invalid_argument("You should run calculations before result request!");
|
|
|
- return static_cast<outputType>(Qbk_);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- // ********************************************************************** //
|
|
|
- // Returns previously calculated Qpr //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- template <typename outputType>
|
|
|
- outputType MultiLayerMie<FloatType>::GetQpr() {
|
|
|
- if (!isMieCalculated_)
|
|
|
- throw std::invalid_argument("You should run calculations before result request!");
|
|
|
- return static_cast<outputType>(Qpr_);
|
|
|
- }
|
|
|
-
|
|
|
+// ********************************************************************** //
|
|
|
+// Returns previously calculated Qext //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+template <typename outputType>
|
|
|
+outputType MultiLayerMie<FloatType>::GetQext() {
|
|
|
+ if (!isMieCalculated_)
|
|
|
+ throw std::invalid_argument(
|
|
|
+ "You should run calculations before result request!");
|
|
|
+ return static_cast<outputType>(Qext_);
|
|
|
+}
|
|
|
|
|
|
- // ********************************************************************** //
|
|
|
- // Returns previously calculated assymetry factor //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- template <typename outputType>
|
|
|
- outputType MultiLayerMie<FloatType>::GetAsymmetryFactor() {
|
|
|
- if (!isMieCalculated_)
|
|
|
- throw std::invalid_argument("You should run calculations before result request!");
|
|
|
- return static_cast<outputType>(asymmetry_factor_);
|
|
|
- }
|
|
|
+// ********************************************************************** //
|
|
|
+// Returns previously calculated Qabs //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+template <typename outputType>
|
|
|
+outputType MultiLayerMie<FloatType>::GetQabs() {
|
|
|
+ if (!isMieCalculated_)
|
|
|
+ throw std::invalid_argument(
|
|
|
+ "You should run calculations before result request!");
|
|
|
+ return static_cast<outputType>(Qabs_);
|
|
|
+}
|
|
|
|
|
|
+// ********************************************************************** //
|
|
|
+// Returns previously calculated Qsca //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+template <typename outputType>
|
|
|
+outputType MultiLayerMie<FloatType>::GetQsca() {
|
|
|
+ if (!isMieCalculated_)
|
|
|
+ throw std::invalid_argument(
|
|
|
+ "You should run calculations before result request!");
|
|
|
+ return static_cast<outputType>(Qsca_);
|
|
|
+}
|
|
|
|
|
|
- // ********************************************************************** //
|
|
|
- // Returns previously calculated Albedo //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- template <typename outputType>
|
|
|
- outputType MultiLayerMie<FloatType>::GetAlbedo() {
|
|
|
- if (!isMieCalculated_)
|
|
|
- throw std::invalid_argument("You should run calculations before result request!");
|
|
|
- return static_cast<outputType>(albedo_);
|
|
|
- }
|
|
|
+// ********************************************************************** //
|
|
|
+// Returns previously calculated Qbk //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+template <typename outputType>
|
|
|
+outputType MultiLayerMie<FloatType>::GetQbk() {
|
|
|
+ if (!isMieCalculated_)
|
|
|
+ throw std::invalid_argument(
|
|
|
+ "You should run calculations before result request!");
|
|
|
+ return static_cast<outputType>(Qbk_);
|
|
|
+}
|
|
|
|
|
|
+// ********************************************************************** //
|
|
|
+// Returns previously calculated Qpr //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+template <typename outputType>
|
|
|
+outputType MultiLayerMie<FloatType>::GetQpr() {
|
|
|
+ if (!isMieCalculated_)
|
|
|
+ throw std::invalid_argument(
|
|
|
+ "You should run calculations before result request!");
|
|
|
+ return static_cast<outputType>(Qpr_);
|
|
|
+}
|
|
|
|
|
|
- // ********************************************************************** //
|
|
|
- // Returns previously calculated S1 //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- std::vector<std::complex<FloatType> > MultiLayerMie<FloatType>::GetS1() {
|
|
|
- if (!isMieCalculated_)
|
|
|
- throw std::invalid_argument("You should run calculations before result request!");
|
|
|
- return S1_;
|
|
|
- }
|
|
|
+// ********************************************************************** //
|
|
|
+// Returns previously calculated asymmetry factor //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+template <typename outputType>
|
|
|
+outputType MultiLayerMie<FloatType>::GetAsymmetryFactor() {
|
|
|
+ if (!isMieCalculated_)
|
|
|
+ throw std::invalid_argument(
|
|
|
+ "You should run calculations before result request!");
|
|
|
+ return static_cast<outputType>(asymmetry_factor_);
|
|
|
+}
|
|
|
|
|
|
+// ********************************************************************** //
|
|
|
+// Returns previously calculated Albedo //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+template <typename outputType>
|
|
|
+outputType MultiLayerMie<FloatType>::GetAlbedo() {
|
|
|
+ if (!isMieCalculated_)
|
|
|
+ throw std::invalid_argument(
|
|
|
+ "You should run calculations before result request!");
|
|
|
+ return static_cast<outputType>(albedo_);
|
|
|
+}
|
|
|
|
|
|
- // ********************************************************************** //
|
|
|
- // Returns previously calculated S2 //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- std::vector<std::complex<FloatType> > MultiLayerMie<FloatType>::GetS2() {
|
|
|
- if (!isMieCalculated_)
|
|
|
- throw std::invalid_argument("You should run calculations before result request!");
|
|
|
- return S2_;
|
|
|
- }
|
|
|
+// ********************************************************************** //
|
|
|
+// Returns previously calculated S1 //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+std::vector<std::complex<FloatType>> MultiLayerMie<FloatType>::GetS1() {
|
|
|
+ if (!isMieCalculated_)
|
|
|
+ throw std::invalid_argument(
|
|
|
+ "You should run calculations before result request!");
|
|
|
+ return S1_;
|
|
|
+}
|
|
|
|
|
|
+// ********************************************************************** //
|
|
|
+// Returns previously calculated S2 //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+std::vector<std::complex<FloatType>> MultiLayerMie<FloatType>::GetS2() {
|
|
|
+ if (!isMieCalculated_)
|
|
|
+ throw std::invalid_argument(
|
|
|
+ "You should run calculations before result request!");
|
|
|
+ return S2_;
|
|
|
+}
|
|
|
|
|
|
+// ********************************************************************** //
|
|
|
+// Modify scattering (theta) angles //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+void MultiLayerMie<FloatType>::SetAngles(const std::vector<FloatType>& angles) {
|
|
|
+ MarkUncalculated();
|
|
|
+ theta_ = angles;
|
|
|
+}
|
|
|
|
|
|
// ********************************************************************** //
|
|
|
- // Modify scattering (theta) angles //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- void MultiLayerMie<FloatType>::SetAngles(const std::vector<FloatType> &angles) {
|
|
|
- MarkUncalculated();
|
|
|
- theta_ = angles;
|
|
|
+// Modify size of all layers //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+void MultiLayerMie<FloatType>::SetLayersSize(
|
|
|
+ const std::vector<FloatType>& layer_size) {
|
|
|
+ MarkUncalculated();
|
|
|
+ size_param_.clear();
|
|
|
+ FloatType 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!");
|
|
|
+ if (prev_layer_size > curr_layer_size)
|
|
|
+ throw std::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);
|
|
|
}
|
|
|
+}
|
|
|
|
|
|
+// ********************************************************************** //
|
|
|
+// Modify refractive index of all layers //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+void MultiLayerMie<FloatType>::SetLayersIndex(
|
|
|
+ const std::vector<std::complex<FloatType>>& index) {
|
|
|
+ MarkUncalculated();
|
|
|
+ refractive_index_ = index;
|
|
|
+}
|
|
|
|
|
|
- // ********************************************************************** //
|
|
|
- // Modify size of all layers //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- void MultiLayerMie<FloatType>::SetLayersSize(const std::vector<FloatType> &layer_size) {
|
|
|
- MarkUncalculated();
|
|
|
- size_param_.clear();
|
|
|
- FloatType 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!");
|
|
|
- if (prev_layer_size > curr_layer_size)
|
|
|
- throw std::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);
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
+// ********************************************************************** //
|
|
|
+// Modify coordinates for field calculation //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+void MultiLayerMie<FloatType>::SetFieldCoords(
|
|
|
+ const std::vector<std::vector<FloatType>>& coords) {
|
|
|
+ if (coords.size() != 3)
|
|
|
+ throw std::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!");
|
|
|
+ coords_ = coords;
|
|
|
+}
|
|
|
|
|
|
+// ********************************************************************** //
|
|
|
+// Modify index of PEC layer //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+void MultiLayerMie<FloatType>::SetPECLayer(int layer_position) {
|
|
|
+ MarkUncalculated();
|
|
|
+ if (layer_position < 0 && layer_position != -1)
|
|
|
+ throw std::invalid_argument("Error! Layers are numbered from 0!");
|
|
|
+ PEC_layer_position_ = layer_position;
|
|
|
+}
|
|
|
|
|
|
+// ********************************************************************** //
|
|
|
+// Set maximun number of terms to be used //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+void MultiLayerMie<FloatType>::SetMaxTerms(int nmax) {
|
|
|
+ MarkUncalculated();
|
|
|
+ nmax_preset_ = nmax;
|
|
|
+}
|
|
|
|
|
|
// ********************************************************************** //
|
|
|
- // Modify refractive index of all layers //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- void MultiLayerMie<FloatType>::SetLayersIndex(const std::vector< std::complex<FloatType> > &index) {
|
|
|
- MarkUncalculated();
|
|
|
- refractive_index_ = index;
|
|
|
- }
|
|
|
+// Get total size parameter of particle //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+FloatType MultiLayerMie<FloatType>::GetSizeParameter() {
|
|
|
+ if (size_param_.size() > 0)
|
|
|
+ return size_param_.back();
|
|
|
+ else
|
|
|
+ return 0;
|
|
|
+}
|
|
|
|
|
|
- // ********************************************************************** //
|
|
|
- // Modify coordinates for field calculation //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- void MultiLayerMie<FloatType>::SetFieldCoords(const std::vector< std::vector<FloatType> > &coords) {
|
|
|
- if (coords.size() != 3)
|
|
|
- throw std::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!");
|
|
|
- coords_ = coords;
|
|
|
- }
|
|
|
+// ********************************************************************** //
|
|
|
+// Mark uncalculated //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+void MultiLayerMie<FloatType>::MarkUncalculated() {
|
|
|
+ isExpCoeffsCalc_ = false;
|
|
|
+ isScaCoeffsCalc_ = false;
|
|
|
|
|
|
+ isMieCalculated_ = false;
|
|
|
+}
|
|
|
+// ********************************************************************** //
|
|
|
+// Clear layer information //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+void MultiLayerMie<FloatType>::ClearLayers() {
|
|
|
+ MarkUncalculated();
|
|
|
+ size_param_.clear();
|
|
|
+ refractive_index_.clear();
|
|
|
+}
|
|
|
|
|
|
- // ********************************************************************** //
|
|
|
- // Modify index of PEC layer //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- void MultiLayerMie<FloatType>::SetPECLayer(int layer_position) {
|
|
|
- MarkUncalculated();
|
|
|
- if (layer_position < 0 && layer_position != -1)
|
|
|
- throw std::invalid_argument("Error! Layers are numbered from 0!");
|
|
|
- PEC_layer_position_ = layer_position;
|
|
|
- }
|
|
|
+// ********************************************************************** //
|
|
|
+// ********************************************************************** //
|
|
|
+// ********************************************************************** //
|
|
|
+// Computational core
|
|
|
+// ********************************************************************** //
|
|
|
+// ********************************************************************** //
|
|
|
+// ********************************************************************** //
|
|
|
|
|
|
+template <typename FloatType>
|
|
|
+unsigned int LeRu_near_field_cutoff(const std::complex<FloatType> zz) {
|
|
|
+ std::complex<double> z = ConvertComplex<double>(zz);
|
|
|
+ auto x = std::abs(z);
|
|
|
+ return std::round(x + 11 * std::pow(x, (1.0 / 3.0)) + 1);
|
|
|
+ // return 10000;
|
|
|
+}
|
|
|
|
|
|
- // ********************************************************************** //
|
|
|
- // Set maximun number of terms to be used //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- void MultiLayerMie<FloatType>::SetMaxTerms(int nmax) {
|
|
|
- MarkUncalculated();
|
|
|
- nmax_preset_ = nmax;
|
|
|
+// ********************************************************************** //
|
|
|
+// Calculate calcNstop - equation (17) //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+unsigned int MultiLayerMie<FloatType>::calcNstop(FloatType xL) {
|
|
|
+ unsigned int nmax = 0;
|
|
|
+ // Wiscombe
|
|
|
+ if (xL < size_param_.back())
|
|
|
+ xL = size_param_.back();
|
|
|
+ if (xL <= 8) {
|
|
|
+ nmax = newround(xL + 4.0 * pow(xL, 1.0 / 3.0) + 1);
|
|
|
+ } else if (xL <= 4200) {
|
|
|
+ nmax = newround(xL + 4.05 * pow(xL, 1.0 / 3.0) + 2);
|
|
|
+ } else {
|
|
|
+ nmax = newround(xL + 4.0 * pow(xL, 1.0 / 3.0) + 2);
|
|
|
}
|
|
|
+ // Use Le Ru cutoff for near field, as a universal one.
|
|
|
+ auto Nstop = nmie::LeRu_near_field_cutoff(std::complex<FloatType>(xL, 0)) + 1;
|
|
|
+ if (Nstop > nmax)
|
|
|
+ nmax = Nstop;
|
|
|
+ return nmax;
|
|
|
+}
|
|
|
|
|
|
// ********************************************************************** //
|
|
|
- // Get total size parameter of particle //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- FloatType MultiLayerMie<FloatType>::GetSizeParameter() {
|
|
|
- if (size_param_.size() > 0)
|
|
|
- return size_param_.back();
|
|
|
+// Maximum number of terms required for the calculation //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+unsigned int MultiLayerMie<FloatType>::calcNmax(FloatType xL) {
|
|
|
+ const int pl = PEC_layer_position_;
|
|
|
+ const unsigned int first_layer = (pl > 0) ? pl : 0;
|
|
|
+ unsigned int ri, riM1, nmax = 0;
|
|
|
+ const std::vector<FloatType>& x = size_param_;
|
|
|
+ const std::vector<std::complex<FloatType>>& m = refractive_index_;
|
|
|
+ nmax = calcNstop(xL);
|
|
|
+ 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 = newround(cabs(x[i] * m[i]));
|
|
|
else
|
|
|
- return 0;
|
|
|
- }
|
|
|
-
|
|
|
- // ********************************************************************** //
|
|
|
- // Mark uncalculated //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- void MultiLayerMie<FloatType>::MarkUncalculated() {
|
|
|
- isExpCoeffsCalc_ = false;
|
|
|
- isScaCoeffsCalc_ = false;
|
|
|
-
|
|
|
- isMieCalculated_ = false;
|
|
|
- }
|
|
|
- // ********************************************************************** //
|
|
|
- // Clear layer information //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- void MultiLayerMie<FloatType>::ClearLayers() {
|
|
|
- MarkUncalculated();
|
|
|
- size_param_.clear();
|
|
|
- refractive_index_.clear();
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- // ********************************************************************** //
|
|
|
- // ********************************************************************** //
|
|
|
- // ********************************************************************** //
|
|
|
- // Computational core
|
|
|
- // ********************************************************************** //
|
|
|
- // ********************************************************************** //
|
|
|
- // ********************************************************************** //
|
|
|
-
|
|
|
- template <typename FloatType>
|
|
|
- unsigned int LeRu_near_field_cutoff(const std::complex<FloatType> zz) {
|
|
|
- std::complex<double> z = ConvertComplex<double>(zz);
|
|
|
- auto x = std::abs(z);
|
|
|
- return std::round(x + 11 * std::pow(x, (1.0 / 3.0)) + 1);
|
|
|
-// return 10000;
|
|
|
- }
|
|
|
-
|
|
|
- // ********************************************************************** //
|
|
|
- // Calculate calcNstop - equation (17) //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- unsigned int MultiLayerMie<FloatType>::calcNstop(FloatType xL) {
|
|
|
- unsigned int nmax = 0;
|
|
|
- //Wiscombe
|
|
|
- if (xL < size_param_.back()) xL = size_param_.back();
|
|
|
- if (xL <= 8) {
|
|
|
- nmax = newround(xL + 4.0*pow(xL, 1.0/3.0) + 1);
|
|
|
- } else if (xL <= 4200) {
|
|
|
- nmax = newround(xL + 4.05*pow(xL, 1.0/3.0) + 2);
|
|
|
- } else {
|
|
|
- nmax = newround(xL + 4.0*pow(xL, 1.0/3.0) + 2);
|
|
|
- }
|
|
|
- //Use Le Ru cutoff for near field, as a universal one.
|
|
|
- auto Nstop = nmie::LeRu_near_field_cutoff(std::complex<FloatType>(xL, 0))+1;
|
|
|
- if (Nstop > nmax) nmax = Nstop;
|
|
|
- return nmax;
|
|
|
+ ri = 0;
|
|
|
+ nmax = std::max(nmax, ri);
|
|
|
+ // first layer is pec, if pec is present
|
|
|
+ if ((i > first_layer) && (static_cast<int>(i - 1) > PEC_layer_position_))
|
|
|
+ riM1 = newround(cabs(x[i - 1] * m[i]));
|
|
|
+ else
|
|
|
+ riM1 = 0;
|
|
|
+ nmax = std::max(nmax, riM1);
|
|
|
}
|
|
|
-
|
|
|
-
|
|
|
- // ********************************************************************** //
|
|
|
- // Maximum number of terms required for the calculation //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- unsigned int MultiLayerMie<FloatType>::calcNmax(FloatType xL) {
|
|
|
- const int pl = PEC_layer_position_;
|
|
|
- const unsigned int first_layer = (pl > 0) ? pl : 0;
|
|
|
- unsigned int ri, riM1, nmax = 0;
|
|
|
- const std::vector<FloatType> &x = size_param_;
|
|
|
- const std::vector<std::complex<FloatType> > &m = refractive_index_;
|
|
|
- nmax = calcNstop(xL);
|
|
|
- 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 = newround(cabs(x[i]*m[i]));
|
|
|
- else
|
|
|
- ri = 0;
|
|
|
- nmax = std::max(nmax, ri);
|
|
|
- // first layer is pec, if pec is present
|
|
|
- if ((i > first_layer) && (static_cast<int>(i - 1) > PEC_layer_position_))
|
|
|
- riM1 = newround(cabs(x[i - 1]* m[i]));
|
|
|
- else
|
|
|
- riM1 = 0;
|
|
|
- nmax = std::max(nmax, riM1);
|
|
|
- }
|
|
|
- nmax += 15; // Final nmax value
|
|
|
+ nmax += 15; // Final nmax value
|
|
|
#ifdef MULTI_PRECISION
|
|
|
- nmax += MULTI_PRECISION; //TODO we may need to use more terms that this for MP computations.
|
|
|
+ nmax += MULTI_PRECISION; // TODO we may need to use more terms that this for
|
|
|
+ // MP computations.
|
|
|
#endif
|
|
|
- // nmax *= nmax;
|
|
|
- // printf("using nmax %i\n", nmax);
|
|
|
- return nmax;
|
|
|
- }
|
|
|
-
|
|
|
+ // nmax *= nmax;
|
|
|
+ // printf("using nmax %i\n", nmax);
|
|
|
+ return nmax;
|
|
|
+}
|
|
|
|
|
|
- // ********************************************************************** //
|
|
|
- // Calculate an - equation (5) //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- std::complex<FloatType> MultiLayerMie<FloatType>::
|
|
|
- calc_an(int n, FloatType XL, std::complex<FloatType> Ha, std::complex<FloatType> mL,
|
|
|
- std::complex<FloatType> PsiXL, std::complex<FloatType> ZetaXL,
|
|
|
- std::complex<FloatType> PsiXLM1, std::complex<FloatType> ZetaXLM1) {
|
|
|
+// ********************************************************************** //
|
|
|
+// Calculate an - equation (5) //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+std::complex<FloatType> MultiLayerMie<FloatType>::calc_an(
|
|
|
+ int n,
|
|
|
+ FloatType XL,
|
|
|
+ std::complex<FloatType> Ha,
|
|
|
+ std::complex<FloatType> mL,
|
|
|
+ std::complex<FloatType> PsiXL,
|
|
|
+ std::complex<FloatType> ZetaXL,
|
|
|
+ std::complex<FloatType> PsiXLM1,
|
|
|
+ std::complex<FloatType> ZetaXLM1) {
|
|
|
+ std::complex<FloatType> Num = (Ha / mL + n / XL) * PsiXL - PsiXLM1;
|
|
|
+ std::complex<FloatType> Denom = (Ha / mL + n / XL) * ZetaXL - ZetaXLM1;
|
|
|
+ // std::cout<< std::setprecision(100)
|
|
|
+ // << "Ql " << PsiXL
|
|
|
+ // << std::endl;
|
|
|
+
|
|
|
+ return Num / Denom;
|
|
|
+}
|
|
|
|
|
|
- std::complex<FloatType> Num = (Ha/mL + n/XL)*PsiXL - PsiXLM1;
|
|
|
- std::complex<FloatType> Denom = (Ha/mL + n/XL)*ZetaXL - ZetaXLM1;
|
|
|
- // std::cout<< std::setprecision(100)
|
|
|
- // << "Ql " << PsiXL
|
|
|
- // << std::endl;
|
|
|
+// ********************************************************************** //
|
|
|
+// Calculate bn - equation (6) //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+std::complex<FloatType> MultiLayerMie<FloatType>::calc_bn(
|
|
|
+ int n,
|
|
|
+ FloatType XL,
|
|
|
+ std::complex<FloatType> Hb,
|
|
|
+ std::complex<FloatType> mL,
|
|
|
+ std::complex<FloatType> PsiXL,
|
|
|
+ std::complex<FloatType> ZetaXL,
|
|
|
+ std::complex<FloatType> PsiXLM1,
|
|
|
+ std::complex<FloatType> ZetaXLM1) {
|
|
|
+ std::complex<FloatType> Num = (mL * Hb + n / XL) * PsiXL - PsiXLM1;
|
|
|
+ std::complex<FloatType> Denom = (mL * Hb + n / XL) * ZetaXL - ZetaXLM1;
|
|
|
+
|
|
|
+ return Num / Denom;
|
|
|
+}
|
|
|
|
|
|
+// ********************************************************************** //
|
|
|
+// Calculates S1 - equation (25a) //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+std::complex<FloatType> MultiLayerMie<FloatType>::calc_S1(
|
|
|
+ int n,
|
|
|
+ std::complex<FloatType> an,
|
|
|
+ std::complex<FloatType> bn,
|
|
|
+ FloatType Pi,
|
|
|
+ FloatType Tau) {
|
|
|
+ return FloatType(n + n + 1) * (Pi * an + Tau * bn) / FloatType(n * n + n);
|
|
|
+}
|
|
|
|
|
|
- return Num/Denom;
|
|
|
+// ********************************************************************** //
|
|
|
+// Calculates S2 - equation (25b) (it's the same as (25a), just switches //
|
|
|
+// Pi and Tau) //
|
|
|
+// ********************************************************************** //
|
|
|
+template <typename FloatType>
|
|
|
+std::complex<FloatType> MultiLayerMie<FloatType>::calc_S2(
|
|
|
+ int n,
|
|
|
+ std::complex<FloatType> an,
|
|
|
+ std::complex<FloatType> bn,
|
|
|
+ FloatType Pi,
|
|
|
+ FloatType Tau) {
|
|
|
+ return calc_S1(n, an, bn, Tau, Pi);
|
|
|
+}
|
|
|
+
|
|
|
+//****************************************************************************
|
|
|
+// This function calculates the logarithmic derivatives of the Riccati-Bessel
|
|
|
+// functions (D1 and D3) for a complex argument (z).
|
|
|
+// Equations (16a), (16b) and (18a) - (18d)
|
|
|
+//
|
|
|
+// Input parameters:
|
|
|
+// z: Complex argument to evaluate D1 and D3
|
|
|
+// nmax_: Maximum number of terms to calculate D1 and D3
|
|
|
+//
|
|
|
+// Output parameters:
|
|
|
+// D1, D3: Logarithmic derivatives of the Riccati-Bessel functions
|
|
|
+//****************************************************************************
|
|
|
+template <typename FloatType>
|
|
|
+void MultiLayerMie<FloatType>::calcD1D3(
|
|
|
+ const std::complex<FloatType> z,
|
|
|
+ std::vector<std::complex<FloatType>>& D1,
|
|
|
+ std::vector<std::complex<FloatType>>& D3) {
|
|
|
+ std::vector<std::complex<FloatType>> PsiZeta(nmax_ + 1);
|
|
|
+ evalDownwardD1(z, D1);
|
|
|
+ evalUpwardD3(z, D1, D3, PsiZeta);
|
|
|
+}
|
|
|
+
|
|
|
+//*****************************************************************************
|
|
|
+// This function calculates the Riccati-Bessel functions (Psi and Zeta) for a
|
|
|
+// complex argument (z).
|
|
|
+// Equations (20a) - (21b)
|
|
|
+//
|
|
|
+// Input parameters:
|
|
|
+// z: Complex argument to evaluate Psi and Zeta
|
|
|
+// nmax: Maximum number of terms to calculate Psi and Zeta
|
|
|
+//
|
|
|
+// Output parameters:
|
|
|
+// Psi, Zeta: Riccati-Bessel functions
|
|
|
+//*****************************************************************************
|
|
|
+template <typename FloatType>
|
|
|
+void MultiLayerMie<FloatType>::calcPsiZeta(
|
|
|
+ std::complex<FloatType> z,
|
|
|
+ std::vector<std::complex<FloatType>>& Psi,
|
|
|
+ std::vector<std::complex<FloatType>>& Zeta) {
|
|
|
+ std::vector<std::complex<FloatType>> D1(nmax_ + 1), D3(nmax_ + 1),
|
|
|
+ PsiZeta(nmax_ + 1);
|
|
|
+ // First, calculate the logarithmic derivatives
|
|
|
+ evalDownwardD1(z, D1);
|
|
|
+ // Now, use the upward recurrence to calculate Psi equations (20ab)
|
|
|
+ evalUpwardPsi(z, D1, Psi);
|
|
|
+ // Now, use the upward recurrence to calculate Psi*Zeta equations (18ad)
|
|
|
+ evalUpwardD3(z, D1, D3, PsiZeta);
|
|
|
+ for (unsigned int i = 0; i < Zeta.size(); i++) {
|
|
|
+ Zeta[i] = PsiZeta[i] / Psi[i];
|
|
|
}
|
|
|
-
|
|
|
-
|
|
|
- // ********************************************************************** //
|
|
|
- // Calculate bn - equation (6) //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- std::complex<FloatType> MultiLayerMie<FloatType>::calc_bn(int n, FloatType XL, std::complex<FloatType> Hb, std::complex<FloatType> mL,
|
|
|
- std::complex<FloatType> PsiXL, std::complex<FloatType> ZetaXL,
|
|
|
- std::complex<FloatType> PsiXLM1, std::complex<FloatType> ZetaXLM1) {
|
|
|
-
|
|
|
- std::complex<FloatType> Num = (mL*Hb + n/XL)*PsiXL - PsiXLM1;
|
|
|
- std::complex<FloatType> Denom = (mL*Hb + n/XL)*ZetaXL - ZetaXLM1;
|
|
|
-
|
|
|
- return Num/Denom;
|
|
|
+ // evalUpwardZeta(z, D3, Zeta);
|
|
|
+}
|
|
|
+
|
|
|
+template <typename FloatType>
|
|
|
+void MultiLayerMie<FloatType>::calcPiTauAllTheta(
|
|
|
+ const double from_Theta,
|
|
|
+ const double to_Theta,
|
|
|
+ std::vector<std::vector<FloatType>>& Pi,
|
|
|
+ std::vector<std::vector<FloatType>>& Tau) {
|
|
|
+ const unsigned int perimeter_points = Pi.size();
|
|
|
+ for (auto& val : Pi)
|
|
|
+ val.resize(available_maximal_nmax_, static_cast<FloatType>(0.0));
|
|
|
+ for (auto& val : Tau)
|
|
|
+ val.resize(available_maximal_nmax_, static_cast<FloatType>(0.0));
|
|
|
+ double delta_Theta =
|
|
|
+ eval_delta<double>(perimeter_points, from_Theta, to_Theta);
|
|
|
+ for (unsigned int i = 0; i < perimeter_points; i++) {
|
|
|
+ auto Theta = static_cast<FloatType>(from_Theta + i * delta_Theta);
|
|
|
+ // Calculate angular functions Pi and Tau
|
|
|
+ calcPiTau(nmm::cos(Theta), Pi[i], Tau[i]);
|
|
|
}
|
|
|
-
|
|
|
-
|
|
|
- // ********************************************************************** //
|
|
|
- // Calculates S1 - equation (25a) //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- std::complex<FloatType> MultiLayerMie<FloatType>::calc_S1(int n, std::complex<FloatType> an, std::complex<FloatType> bn,
|
|
|
- FloatType Pi, FloatType Tau) {
|
|
|
- return FloatType(n + n + 1)*(Pi*an + Tau*bn)/FloatType(n*n + n);
|
|
|
+}
|
|
|
+
|
|
|
+//*******************************************************************************
|
|
|
+// This function calculates Pi and Tau for a given value of cos(Theta).
|
|
|
+// Equations (26a) - (26c)
|
|
|
+//
|
|
|
+// Input parameters:
|
|
|
+// nmax_: Maximum number of terms to calculate Pi and Tau
|
|
|
+// nTheta: Number of scattering angles
|
|
|
+// Theta: Array containing all the scattering angles where the scattering
|
|
|
+// amplitudes will be calculated
|
|
|
+//
|
|
|
+// Output parameters:
|
|
|
+// Pi, Tau: Angular functions Pi and Tau, as defined in equations (26a) -
|
|
|
+// (26c)
|
|
|
+//*******************************************************************************
|
|
|
+template <typename FloatType>
|
|
|
+void MultiLayerMie<FloatType>::calcPiTau(const FloatType& costheta,
|
|
|
+ std::vector<FloatType>& Pi,
|
|
|
+ std::vector<FloatType>& Tau) {
|
|
|
+ int nmax = Pi.size();
|
|
|
+ if (Pi.size() != Tau.size())
|
|
|
+ throw std::invalid_argument(
|
|
|
+ "Error! Pi and Tau vectors should have the same size!");
|
|
|
+
|
|
|
+ //****************************************************//
|
|
|
+ // Equations (26a) - (26c) //
|
|
|
+ //****************************************************//
|
|
|
+ // Initialize Pi and Tau
|
|
|
+ Pi[0] = 1.0; // n=1
|
|
|
+ Tau[0] = costheta;
|
|
|
+ // Calculate the actual values
|
|
|
+ if (nmax > 1) {
|
|
|
+ Pi[1] = 3 * costheta * Pi[0]; // n=2
|
|
|
+ Tau[1] = 2 * costheta * Pi[1] - 3 * Pi[0];
|
|
|
+ for (int i = 2; i < nmax; i++) { // n=[3..nmax_]
|
|
|
+ Pi[i] = ((i + i + 1) * costheta * Pi[i - 1] - (i + 1) * Pi[i - 2]) / i;
|
|
|
+ Tau[i] = (i + 1) * costheta * Pi[i] - (i + 2) * Pi[i - 1];
|
|
|
+ }
|
|
|
}
|
|
|
-
|
|
|
-
|
|
|
- // ********************************************************************** //
|
|
|
- // Calculates S2 - equation (25b) (it's the same as (25a), just switches //
|
|
|
- // Pi and Tau) //
|
|
|
- // ********************************************************************** //
|
|
|
- template <typename FloatType>
|
|
|
- std::complex<FloatType> MultiLayerMie<FloatType>::calc_S2(int n, std::complex<FloatType> an, std::complex<FloatType> bn,
|
|
|
- FloatType Pi, FloatType Tau) {
|
|
|
- return calc_S1(n, an, bn, Tau, Pi);
|
|
|
+} // end of MultiLayerMie::calcPiTau(...)
|
|
|
+
|
|
|
+//*****************************************************************************
|
|
|
+// This function calculates vector spherical harmonics (eq. 4.50, p. 95 BH),
|
|
|
+// required to calculate the near-field parameters.
|
|
|
+//
|
|
|
+// Input parameters:
|
|
|
+// Rho: Radial distance
|
|
|
+// Phi: Azimuthal angle
|
|
|
+// Theta: Polar angle
|
|
|
+// rn: Either the spherical Ricatti-Bessel function of first or third kind
|
|
|
+// Dn: Logarithmic derivative of rn
|
|
|
+// Pi, Tau: Angular functions Pi and Tau
|
|
|
+// n: Order of vector spherical harmonics
|
|
|
+//
|
|
|
+// Output parameters:
|
|
|
+// Mo1n, Me1n, No1n, Ne1n: Complex vector spherical harmonics
|
|
|
+//*****************************************************************************
|
|
|
+template <typename FloatType>
|
|
|
+template <typename evalType>
|
|
|
+void MultiLayerMie<FloatType>::calcSpherHarm(
|
|
|
+ const std::complex<evalType> Rho,
|
|
|
+ const evalType Theta,
|
|
|
+ const evalType Phi,
|
|
|
+ const std::complex<evalType>& rn,
|
|
|
+ const std::complex<evalType>& Dn,
|
|
|
+ const evalType& Pi,
|
|
|
+ const evalType& Tau,
|
|
|
+ const evalType& n,
|
|
|
+ std::vector<std::complex<evalType>>& Mo1n,
|
|
|
+ std::vector<std::complex<evalType>>& Me1n,
|
|
|
+ std::vector<std::complex<evalType>>& No1n,
|
|
|
+ std::vector<std::complex<evalType>>& Ne1n) {
|
|
|
+ // using eq 4.50 in BH
|
|
|
+ std::complex<evalType> c_zero(0.0, 0.0);
|
|
|
+
|
|
|
+ // using nmm::sin;
|
|
|
+ // using nmm::cos;
|
|
|
+ auto sin_Phi = sin_t(Phi);
|
|
|
+ auto cos_Phi = cos_t(Phi);
|
|
|
+ auto sin_Theta = sin(Theta);
|
|
|
+ Mo1n[0] = c_zero;
|
|
|
+ Mo1n[1] = cos_Phi * Pi * rn / Rho;
|
|
|
+ Mo1n[2] = -sin_Phi * Tau * rn / Rho;
|
|
|
+
|
|
|
+ Me1n[0] = c_zero;
|
|
|
+ Me1n[1] = -sin_Phi * Pi * rn / Rho;
|
|
|
+ Me1n[2] = -cos_Phi * Tau * rn / Rho;
|
|
|
+
|
|
|
+ No1n[0] = sin_Phi * (n * n + n) * sin_Theta * Pi * rn / Rho / Rho;
|
|
|
+ No1n[1] = sin_Phi * Tau * Dn * rn / Rho;
|
|
|
+ No1n[2] = cos_Phi * Pi * Dn * rn / Rho;
|
|
|
+
|
|
|
+ Ne1n[0] = cos_Phi * (n * n + n) * sin_Theta * Pi * rn / Rho / Rho;
|
|
|
+ Ne1n[1] = cos_Phi * Tau * Dn * rn / Rho;
|
|
|
+ Ne1n[2] = -sin_Phi * Pi * Dn * rn / Rho;
|
|
|
+} // end of MultiLayerMie::calcSpherHarm(...)
|
|
|
+
|
|
|
+//********************************************************************************
|
|
|
+// This function calculates the scattering coefficients required to calculate
|
|
|
+// both the near- and far-field parameters.
|
|
|
+//
|
|
|
+// Input parameters:
|
|
|
+// L: Number of layers
|
|
|
+// pl: Index of PEC layer. If there is none just send -1
|
|
|
+// x: Array containing the size parameters of the layers [0..L-1]
|
|
|
+// m: Array containing the relative refractive indexes of the layers [0..L-1]
|
|
|
+// nmax: Maximum number of multipolar expansion terms to be used for the
|
|
|
+// calculations. Only use it if you know what you are doing, otherwise
|
|
|
+// set this parameter to -1 and the function will calculate it.
|
|
|
+//
|
|
|
+// Output parameters:
|
|
|
+// an, bn: Complex scattering amplitudes
|
|
|
+//
|
|
|
+// Return value:
|
|
|
+// Number of multipolar expansion terms used for the calculations
|
|
|
+//********************************************************************************
|
|
|
+template <typename FloatType>
|
|
|
+void MultiLayerMie<FloatType>::calcScattCoeffs() {
|
|
|
+ isScaCoeffsCalc_ = false;
|
|
|
+ an_.clear();
|
|
|
+ bn_.clear();
|
|
|
+
|
|
|
+ const std::vector<FloatType>& x = size_param_;
|
|
|
+ const std::vector<std::complex<FloatType>>& m = refractive_index_;
|
|
|
+ const int& pl = PEC_layer_position_;
|
|
|
+ const int L = refractive_index_.size();
|
|
|
+
|
|
|
+ //************************************************************************//
|
|
|
+ // Calculate the index of the first layer. It can be either 0 (default) //
|
|
|
+ // or the index of the outermost PEC layer. In the latter case all layers //
|
|
|
+ // below the PEC are discarded. //
|
|
|
+ // ***********************************************************************//
|
|
|
+ int fl = (pl > 0) ? pl : 0;
|
|
|
+ if (nmax_preset_ <= 0)
|
|
|
+ nmax_ = calcNmax();
|
|
|
+ else
|
|
|
+ nmax_ = nmax_preset_;
|
|
|
+
|
|
|
+ std::complex<FloatType> 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 //
|
|
|
+ // are the arrays for representing D1, D3 and Q because they need a value //
|
|
|
+ // for the index 0 (zero), hence it is important to consider this shift //
|
|
|
+ // between different arrays. The change was done to optimize memory usage. //
|
|
|
+ //**************************************************************************//
|
|
|
+ // Allocate memory to the arrays
|
|
|
+ std::vector<std::complex<FloatType>> D1_mlxl(nmax_ + 1), D1_mlxlM1(nmax_ + 1),
|
|
|
+ D3_mlxl(nmax_ + 1), D3_mlxlM1(nmax_ + 1);
|
|
|
+
|
|
|
+ std::vector<std::vector<std::complex<FloatType>>> Q(L), Ha(L), Hb(L);
|
|
|
+
|
|
|
+ for (int l = 0; l < L; l++) {
|
|
|
+ Q[l].resize(nmax_ + 1, static_cast<FloatType>(0.0));
|
|
|
+ Ha[l].resize(nmax_, static_cast<FloatType>(0.0));
|
|
|
+ Hb[l].resize(nmax_, static_cast<FloatType>(0.0));
|
|
|
}
|
|
|
|
|
|
+ an_.resize(nmax_, static_cast<FloatType>(0.0));
|
|
|
+ bn_.resize(nmax_, static_cast<FloatType>(0.0));
|
|
|
|
|
|
- //**********************************************************************************//
|
|
|
- // This function calculates the logarithmic derivatives of the Riccati-Bessel //
|
|
|
- // functions (D1 and D3) for a complex argument (z). //
|
|
|
- // Equations (16a), (16b) and (18a) - (18d) //
|
|
|
- // //
|
|
|
- // Input parameters: //
|
|
|
- // z: Complex argument to evaluate D1 and D3 //
|
|
|
- // nmax_: Maximum number of terms to calculate D1 and D3 //
|
|
|
- // //
|
|
|
- // Output parameters: //
|
|
|
- // D1, D3: Logarithmic derivatives of the Riccati-Bessel functions //
|
|
|
- //**********************************************************************************//
|
|
|
- template <typename FloatType>
|
|
|
- void MultiLayerMie<FloatType>::calcD1D3(const std::complex<FloatType> z,
|
|
|
- std::vector<std::complex<FloatType> > &D1,
|
|
|
- std::vector<std::complex<FloatType> > &D3) {
|
|
|
- std::vector<std::complex<FloatType> > PsiZeta(nmax_+1);
|
|
|
- evalDownwardD1(z, D1);
|
|
|
- evalUpwardD3 (z, D1, D3, PsiZeta);
|
|
|
- }
|
|
|
-
|
|
|
+ std::vector<std::complex<FloatType>> PsiXL(nmax_ + 1), ZetaXL(nmax_ + 1);
|
|
|
|
|
|
- //**********************************************************************************//
|
|
|
- // This function calculates the Riccati-Bessel functions (Psi and Zeta) for a //
|
|
|
- // complex argument (z). //
|
|
|
- // Equations (20a) - (21b) //
|
|
|
- // //
|
|
|
- // Input parameters: //
|
|
|
- // z: Complex argument to evaluate Psi and Zeta //
|
|
|
- // nmax: Maximum number of terms to calculate Psi and Zeta //
|
|
|
- // //
|
|
|
- // Output parameters: //
|
|
|
- // Psi, Zeta: Riccati-Bessel functions //
|
|
|
- //**********************************************************************************//
|
|
|
- template <typename FloatType>
|
|
|
- void MultiLayerMie<FloatType>::calcPsiZeta(std::complex<FloatType> z,
|
|
|
- std::vector<std::complex<FloatType> > &Psi,
|
|
|
- std::vector<std::complex<FloatType> > &Zeta) {
|
|
|
- std::vector<std::complex<FloatType> > D1(nmax_ + 1), D3(nmax_ + 1),
|
|
|
- PsiZeta(nmax_+1);
|
|
|
- // First, calculate the logarithmic derivatives
|
|
|
- evalDownwardD1(z, D1);
|
|
|
- // Now, use the upward recurrence to calculate Psi equations (20ab)
|
|
|
- evalUpwardPsi(z, D1, Psi);
|
|
|
- // Now, use the upward recurrence to calculate Psi*Zeta equations (18ad)
|
|
|
- evalUpwardD3 (z, D1, D3, PsiZeta);
|
|
|
- for (unsigned int i = 0; i < Zeta.size(); i++) {
|
|
|
- Zeta[i] = PsiZeta[i]/Psi[i];
|
|
|
+ //*************************************************//
|
|
|
+ // 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<FloatType>(0.0, -1.0);
|
|
|
+ D3_mlxl[n] = std::complex<FloatType>(0.0, 1.0);
|
|
|
}
|
|
|
-// evalUpwardZeta(z, D3, Zeta);
|
|
|
+ } else { // Regular layer
|
|
|
+ z1 = x[fl] * m[fl];
|
|
|
+ // Calculate D1 and D3
|
|
|
+ calcD1D3(z1, D1_mlxl, D3_mlxl);
|
|
|
}
|
|
|
|
|
|
-
|
|
|
- template <typename FloatType>
|
|
|
- void MultiLayerMie<FloatType>::calcPiTauAllTheta(const double from_Theta, const double to_Theta,
|
|
|
- std::vector<std::vector<FloatType> > &Pi,
|
|
|
- std::vector<std::vector<FloatType> > &Tau) {
|
|
|
- const unsigned int perimeter_points = Pi.size();
|
|
|
- for (auto &val:Pi) val.resize(available_maximal_nmax_, static_cast<FloatType>(0.0));
|
|
|
- for (auto &val:Tau) val.resize(available_maximal_nmax_, static_cast<FloatType>(0.0));
|
|
|
- double delta_Theta = eval_delta<double>(perimeter_points, from_Theta, to_Theta);
|
|
|
- for (unsigned int i=0; i < perimeter_points; i++) {
|
|
|
- auto Theta = static_cast<FloatType>(from_Theta + i*delta_Theta);
|
|
|
- // Calculate angular functions Pi and Tau
|
|
|
- calcPiTau(nmm::cos(Theta), Pi[i], Tau[i]);
|
|
|
- }
|
|
|
+ //******************************************************************//
|
|
|
+ // Calculate Ha and Hb in the first layer - equations (7a) and (8a) //
|
|
|
+ //******************************************************************//
|
|
|
+ for (int n = 0; n < nmax_; n++) {
|
|
|
+ Ha[fl][n] = D1_mlxl[n + 1];
|
|
|
+ Hb[fl][n] = D1_mlxl[n + 1];
|
|
|
}
|
|
|
-
|
|
|
-
|
|
|
- //**********************************************************************************//
|
|
|
- // This function calculates Pi and Tau for a given value of cos(Theta). //
|
|
|
- // Equations (26a) - (26c) //
|
|
|
- // //
|
|
|
- // Input parameters: //
|
|
|
- // nmax_: Maximum number of terms to calculate Pi and Tau //
|
|
|
- // nTheta: Number of scattering angles //
|
|
|
- // Theta: Array containing all the scattering angles where the scattering //
|
|
|
- // amplitudes will be calculated //
|
|
|
- // //
|
|
|
- // Output parameters: //
|
|
|
- // Pi, Tau: Angular functions Pi and Tau, as defined in equations (26a) - (26c) //
|
|
|
- //**********************************************************************************//
|
|
|
- template <typename FloatType>
|
|
|
- void MultiLayerMie<FloatType>::calcPiTau(const FloatType &costheta,
|
|
|
- std::vector<FloatType> &Pi, std::vector<FloatType> &Tau) {
|
|
|
-
|
|
|
- int nmax = Pi.size();
|
|
|
- if (Pi.size() != Tau.size())
|
|
|
- throw std::invalid_argument("Error! Pi and Tau vectors should have the same size!");
|
|
|
-
|
|
|
- //****************************************************//
|
|
|
- // Equations (26a) - (26c) //
|
|
|
- //****************************************************//
|
|
|
- // Initialize Pi and Tau
|
|
|
- Pi[0] = 1.0; // n=1
|
|
|
- Tau[0] = costheta;
|
|
|
- // Calculate the actual values
|
|
|
- if (nmax > 1) {
|
|
|
- Pi[1] = 3*costheta*Pi[0]; //n=2
|
|
|
- Tau[1] = 2*costheta*Pi[1] - 3*Pi[0];
|
|
|
- for (int i = 2; i < nmax; i++) { //n=[3..nmax_]
|
|
|
- Pi[i] = ((i + i + 1)*costheta*Pi[i - 1] - (i + 1)*Pi[i - 2])/i;
|
|
|
- Tau[i] = (i + 1)*costheta*Pi[i] - (i + 2)*Pi[i - 1];
|
|
|
- }
|
|
|
- }
|
|
|
- } // end of MultiLayerMie::calcPiTau(...)
|
|
|
-
|
|
|
-
|
|
|
- //**********************************************************************************//
|
|
|
- // This function calculates vector spherical harmonics (eq. 4.50, p. 95 BH), //
|
|
|
- // required to calculate the near-field parameters. //
|
|
|
- // //
|
|
|
- // Input parameters: //
|
|
|
- // Rho: Radial distance //
|
|
|
- // Phi: Azimuthal angle //
|
|
|
- // Theta: Polar angle //
|
|
|
- // rn: Either the spherical Ricatti-Bessel function of first or third kind //
|
|
|
- // Dn: Logarithmic derivative of rn //
|
|
|
- // Pi, Tau: Angular functions Pi and Tau //
|
|
|
- // n: Order of vector spherical harmonics //
|
|
|
- // //
|
|
|
- // Output parameters: //
|
|
|
- // Mo1n, Me1n, No1n, Ne1n: Complex vector spherical harmonics //
|
|
|
- //**********************************************************************************//
|
|
|
- template <typename FloatType> template <typename evalType>
|
|
|
- void MultiLayerMie<FloatType>::calcSpherHarm(const std::complex<evalType> Rho, const evalType Theta, const evalType Phi,
|
|
|
- const std::complex<evalType> &rn, const std::complex<evalType> &Dn,
|
|
|
- const evalType &Pi, const evalType &Tau, const evalType &n,
|
|
|
- std::vector<std::complex<evalType> > &Mo1n, std::vector<std::complex<evalType> > &Me1n,
|
|
|
- std::vector<std::complex<evalType> > &No1n, std::vector<std::complex<evalType> > &Ne1n) {
|
|
|
-
|
|
|
- // using eq 4.50 in BH
|
|
|
- std::complex<evalType> c_zero(0.0, 0.0);
|
|
|
-
|
|
|
-// using nmm::sin;
|
|
|
-// using nmm::cos;
|
|
|
- auto sin_Phi = sin_t(Phi);
|
|
|
- auto cos_Phi = cos_t(Phi);
|
|
|
- auto sin_Theta = sin(Theta);
|
|
|
- Mo1n[0] = c_zero;
|
|
|
- Mo1n[1] = cos_Phi*Pi*rn/Rho;
|
|
|
- Mo1n[2] = -sin_Phi*Tau*rn/Rho;
|
|
|
-
|
|
|
- Me1n[0] = c_zero;
|
|
|
- Me1n[1] = -sin_Phi*Pi*rn/Rho;
|
|
|
- Me1n[2] = -cos_Phi*Tau*rn/Rho;
|
|
|
-
|
|
|
- No1n[0] = sin_Phi*(n*n + n)*sin_Theta*Pi*rn/Rho/Rho;
|
|
|
- No1n[1] = sin_Phi*Tau*Dn*rn/Rho;
|
|
|
- No1n[2] = cos_Phi*Pi*Dn*rn/Rho;
|
|
|
-
|
|
|
- Ne1n[0] = cos_Phi*(n*n + n)*sin_Theta*Pi*rn/Rho/Rho;
|
|
|
- Ne1n[1] = cos_Phi*Tau*Dn*rn/Rho;
|
|
|
- Ne1n[2] = -sin_Phi*Pi*Dn*rn/Rho;
|
|
|
- } // end of MultiLayerMie::calcSpherHarm(...)
|
|
|
-
|
|
|
-
|
|
|
- //**********************************************************************************//
|
|
|
- // This function calculates the scattering coefficients required to calculate //
|
|
|
- // both the near- and far-field parameters. //
|
|
|
- // //
|
|
|
- // Input parameters: //
|
|
|
- // L: Number of layers //
|
|
|
- // pl: Index of PEC layer. If there is none just send -1 //
|
|
|
- // x: Array containing the size parameters of the layers [0..L-1] //
|
|
|
- // m: Array containing the relative refractive indexes of the layers [0..L-1] //
|
|
|
- // nmax: Maximum number of multipolar expansion terms to be used for the //
|
|
|
- // calculations. Only use it if you know what you are doing, otherwise //
|
|
|
- // set this parameter to -1 and the function will calculate it. //
|
|
|
- // //
|
|
|
- // Output parameters: //
|
|
|
- // an, bn: Complex scattering amplitudes //
|
|
|
- // //
|
|
|
- // Return value: //
|
|
|
- // Number of multipolar expansion terms used for the calculations //
|
|
|
- //**********************************************************************************//
|
|
|
- template <typename FloatType>
|
|
|
- void MultiLayerMie<FloatType>::calcScattCoeffs() {
|
|
|
-
|
|
|
- isScaCoeffsCalc_ = false;
|
|
|
- an_.clear();
|
|
|
- bn_.clear();
|
|
|
-
|
|
|
- const std::vector<FloatType> &x = size_param_;
|
|
|
- const std::vector<std::complex<FloatType> > &m = refractive_index_;
|
|
|
- const int &pl = PEC_layer_position_;
|
|
|
- const int L = refractive_index_.size();
|
|
|
-
|
|
|
-
|
|
|
- //************************************************************************//
|
|
|
- // Calculate the index of the first layer. It can be either 0 (default) //
|
|
|
- // or the index of the outermost PEC layer. In the latter case all layers //
|
|
|
- // below the PEC are discarded. //
|
|
|
- // ***********************************************************************//
|
|
|
- int fl = (pl > 0) ? pl : 0;
|
|
|
- if (nmax_preset_ <= 0) nmax_ = calcNmax();
|
|
|
- else nmax_ = nmax_preset_;
|
|
|
-
|
|
|
- std::complex<FloatType> 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 //
|
|
|
- // are the arrays for representing D1, D3 and Q because they need a value //
|
|
|
- // for the index 0 (zero), hence it is important to consider this shift //
|
|
|
- // between different arrays. The change was done to optimize memory usage. //
|
|
|
- //**************************************************************************//
|
|
|
- // Allocate memory to the arrays
|
|
|
- std::vector<std::complex<FloatType> > D1_mlxl(nmax_ + 1), D1_mlxlM1(nmax_ + 1),
|
|
|
- D3_mlxl(nmax_ + 1), D3_mlxlM1(nmax_ + 1);
|
|
|
-
|
|
|
- std::vector<std::vector<std::complex<FloatType> > > Q(L), Ha(L), Hb(L);
|
|
|
-
|
|
|
- for (int l = 0; l < L; l++) {
|
|
|
- Q[l].resize(nmax_ + 1, static_cast<FloatType>(0.0));
|
|
|
- Ha[l].resize(nmax_, static_cast<FloatType>(0.0));
|
|
|
- Hb[l].resize(nmax_, static_cast<FloatType>(0.0));
|
|
|
- }
|
|
|
-
|
|
|
- an_.resize(nmax_, static_cast<FloatType>(0.0));
|
|
|
- bn_.resize(nmax_, static_cast<FloatType>(0.0));
|
|
|
-
|
|
|
- std::vector<std::complex<FloatType> > PsiXL(nmax_ + 1), ZetaXL(nmax_ + 1);
|
|
|
+ //*****************************************************//
|
|
|
+ // Iteration from the second layer to the last one (L) //
|
|
|
+ //*****************************************************//
|
|
|
+ std::complex<FloatType> Temp, Num, Denom;
|
|
|
+ std::complex<FloatType> G1, G2;
|
|
|
+ for (int l = fl + 1; l < L; l++) {
|
|
|
+ //************************************************************//
|
|
|
+ // Calculate D1 and D3 for z1 and z2 in the layers fl + 1..L //
|
|
|
+ //************************************************************//
|
|
|
+ z1 = x[l] * m[l];
|
|
|
+ z2 = x[l - 1] * m[l];
|
|
|
+ // Calculate D1 and D3 for z1
|
|
|
+ calcD1D3(z1, D1_mlxl, D3_mlxl);
|
|
|
+ // Calculate D1 and D3 for z2
|
|
|
+ calcD1D3(z2, D1_mlxlM1, D3_mlxlM1);
|
|
|
|
|
|
//*************************************************//
|
|
|
- // Calculate D1 and D3 for z1 in the first layer //
|
|
|
+ // Calculate Q, Ha and Hb in the layers fl + 1..L //
|
|
|
//*************************************************//
|
|
|
- if (fl == pl) { // PEC layer
|
|
|
- for (int n = 0; n <= nmax_; n++) {
|
|
|
- D1_mlxl[n] = std::complex<FloatType>(0.0, - 1.0);
|
|
|
- D3_mlxl[n] = std::complex<FloatType>(0.0, 1.0);
|
|
|
- }
|
|
|
- } else { // Regular layer
|
|
|
- z1 = x[fl]* m[fl];
|
|
|
- // Calculate D1 and D3
|
|
|
- calcD1D3(z1, D1_mlxl, D3_mlxl);
|
|
|
+ // Upward recurrence for Q - equations (19a) and (19b)
|
|
|
+ Num =
|
|
|
+ std::complex<FloatType>(nmm::exp(-2.0 * (z1.imag() - z2.imag())), 0.0) *
|
|
|
+ std::complex<FloatType>(
|
|
|
+ nmm::cos(-2.0 * z2.real()) - nmm::exp(-2.0 * z2.imag()),
|
|
|
+ nmm::sin(-2.0 * z2.real()));
|
|
|
+ Denom = std::complex<FloatType>(
|
|
|
+ nmm::cos(-2.0 * z1.real()) - nmm::exp(-2.0 * z1.imag()),
|
|
|
+ nmm::sin(-2.0 * z1.real()));
|
|
|
+ Q[l][0] = Num / Denom;
|
|
|
+
|
|
|
+ for (int n = 1; n <= nmax_; n++) {
|
|
|
+ Num = (z1 * D1_mlxl[n] + FloatType(n)) *
|
|
|
+ (FloatType(n) - z1 * D3_mlxl[n - 1]);
|
|
|
+ Denom = (z2 * D1_mlxlM1[n] + FloatType(n)) *
|
|
|
+ (FloatType(n) - z2 * D3_mlxlM1[n - 1]);
|
|
|
+ Q[l][n] = ((pow2(x[l - 1] / x[l]) * Q[l][n - 1]) * Num) / Denom;
|
|
|
+ }
|
|
|
+ // Upward recurrence for Ha and Hb - equations (7b), (8b) and (12) - (15)
|
|
|
+ for (int n = 1; n <= nmax_; n++) {
|
|
|
+ // Ha
|
|
|
+ if ((l - 1) == pl) { // The layer below the current one is a PEC layer
|
|
|
+ G1 = -D1_mlxlM1[n];
|
|
|
+ G2 = -D3_mlxlM1[n];
|
|
|
+ } else {
|
|
|
+ G1 = (m[l] * Ha[l - 1][n - 1]) - (m[l - 1] * D1_mlxlM1[n]);
|
|
|
+ G2 = (m[l] * Ha[l - 1][n - 1]) - (m[l - 1] * D3_mlxlM1[n]);
|
|
|
+ } // end of if PEC
|
|
|
+ Temp = Q[l][n] * G1;
|
|
|
+ Num = (G2 * D1_mlxl[n]) - (Temp * D3_mlxl[n]);
|
|
|
+ Denom = G2 - Temp;
|
|
|
+ Ha[l][n - 1] = Num / Denom;
|
|
|
+ // Hb
|
|
|
+ if ((l - 1) == pl) { // The layer below the current one is a PEC layer
|
|
|
+ G1 = Hb[l - 1][n - 1];
|
|
|
+ G2 = Hb[l - 1][n - 1];
|
|
|
+ } else {
|
|
|
+ G1 = (m[l - 1] * Hb[l - 1][n - 1]) - (m[l] * D1_mlxlM1[n]);
|
|
|
+ G2 = (m[l - 1] * Hb[l - 1][n - 1]) - (m[l] * D3_mlxlM1[n]);
|
|
|
+ } // end of if PEC
|
|
|
+
|
|
|
+ Temp = Q[l][n] * G1;
|
|
|
+ Num = (G2 * D1_mlxl[n]) - (Temp * D3_mlxl[n]);
|
|
|
+ Denom = (G2 - Temp);
|
|
|
+ Hb[l][n - 1] = (Num / Denom);
|
|
|
+ } // end of for Ha and Hb terms
|
|
|
+ } // end of for layers iteration
|
|
|
+
|
|
|
+ //**************************************//
|
|
|
+ // Calculate Psi and Zeta for XL //
|
|
|
+ //**************************************//
|
|
|
+ // Calculate PsiXL and ZetaXL
|
|
|
+ calcPsiZeta(std::complex<FloatType>(x[L - 1], 0.0), PsiXL, ZetaXL);
|
|
|
+
|
|
|
+ //*********************************************************************//
|
|
|
+ // Finally, we calculate the scattering coefficients (an and bn) and //
|
|
|
+ // the angular functions (Pi and Tau). Note that for these arrays the //
|
|
|
+ // first layer is 0 (zero), in future versions all arrays will follow //
|
|
|
+ // this convention to save memory. (13 Nov, 2014) //
|
|
|
+ //*********************************************************************//
|
|
|
+ FloatType a0 = 0, b0 = 0;
|
|
|
+ for (int n = 0; n < nmax_; n++) {
|
|
|
+ //********************************************************************//
|
|
|
+ // Expressions for calculating an and bn coefficients are not valid if //
|
|
|
+ // there is only one PEC layer (ie, for a simple PEC sphere). //
|
|
|
+ //********************************************************************//
|
|
|
+ if (pl < (L - 1)) {
|
|
|
+ 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<FloatType>(0.0, 0.0),
|
|
|
+ std::complex<FloatType>(1.0, 0.0), PsiXL[n + 1],
|
|
|
+ ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
|
|
|
+ bn_[n] = PsiXL[n + 1] / ZetaXL[n + 1];
|
|
|
+ }
|
|
|
+ if (n == 0) {
|
|
|
+ a0 = cabs(an_[0]);
|
|
|
+ b0 = cabs(bn_[0]);
|
|
|
+ }
|
|
|
+ if (n == nmax_ - 1 && nmax_preset_ <= 0 &&
|
|
|
+ (cabs(an_[n]) / a0 > convergence_threshold_ &&
|
|
|
+ cabs(bn_[n]) / b0 > convergence_threshold_)) {
|
|
|
+ std::cout << "Failed to converge in Mie series for nmax=" << nmax_
|
|
|
+ << std::endl;
|
|
|
+ std::cout << "convergence threshold: " << convergence_threshold_
|
|
|
+ << std::endl;
|
|
|
+ std::cout << "Mie series a[nmax]/a[1]:" << cabs(an_[n]) / a0
|
|
|
+ << " and b[nmax]/b[1]:" << cabs(bn_[n]) / b0 << std::endl;
|
|
|
}
|
|
|
|
|
|
-
|
|
|
- //******************************************************************//
|
|
|
- // Calculate Ha and Hb in the first layer - equations (7a) and (8a) //
|
|
|
- //******************************************************************//
|
|
|
- for (int n = 0; n < nmax_; n++) {
|
|
|
- Ha[fl][n] = D1_mlxl[n + 1];
|
|
|
- Hb[fl][n] = D1_mlxl[n + 1];
|
|
|
+ // TODO seems to provide not enough terms for near-field calclulation.
|
|
|
+ // if (cabs(an_[n]) / a0 < convergence_threshold_ &&
|
|
|
+ // cabs(bn_[n]) / b0 < convergence_threshold_) {
|
|
|
+ // if (nmax_preset_ <= 0) nmax_ = n;
|
|
|
+ // break;
|
|
|
+ // }
|
|
|
+
|
|
|
+ if (nmm::isnan(an_[n].real()) || nmm::isnan(an_[n].imag()) ||
|
|
|
+ nmm::isnan(bn_[n].real()) || nmm::isnan(bn_[n].imag())) {
|
|
|
+ std::cout
|
|
|
+ << "nmax value was changed due to unexpected error!!! New values is "
|
|
|
+ << n << " (was " << nmax_ << ")" << std::endl;
|
|
|
+ nmax_ = n;
|
|
|
+ break;
|
|
|
}
|
|
|
- //*****************************************************//
|
|
|
- // Iteration from the second layer to the last one (L) //
|
|
|
- //*****************************************************//
|
|
|
- std::complex<FloatType> Temp, Num, Denom;
|
|
|
- std::complex<FloatType> G1, G2;
|
|
|
- for (int l = fl + 1; l < L; l++) {
|
|
|
- //************************************************************//
|
|
|
- //Calculate D1 and D3 for z1 and z2 in the layers fl + 1..L //
|
|
|
- //************************************************************//
|
|
|
- z1 = x[l]*m[l];
|
|
|
- z2 = x[l - 1]*m[l];
|
|
|
- //Calculate D1 and D3 for z1
|
|
|
- calcD1D3(z1, D1_mlxl, D3_mlxl);
|
|
|
- //Calculate D1 and D3 for z2
|
|
|
- calcD1D3(z2, D1_mlxlM1, D3_mlxlM1);
|
|
|
-
|
|
|
- //*************************************************//
|
|
|
- //Calculate Q, Ha and Hb in the layers fl + 1..L //
|
|
|
- //*************************************************//
|
|
|
- // Upward recurrence for Q - equations (19a) and (19b)
|
|
|
- Num = std::complex<FloatType>(nmm::exp(-2.0*(z1.imag() - z2.imag())), 0.0)
|
|
|
- *std::complex<FloatType>(nmm::cos(-2.0*z2.real()) - nmm::exp(-2.0*z2.imag()), nmm::sin(-2.0*z2.real()));
|
|
|
- Denom = std::complex<FloatType>(nmm::cos(-2.0*z1.real()) - nmm::exp(-2.0*z1.imag()), nmm::sin(-2.0*z1.real()));
|
|
|
- Q[l][0] = Num/Denom;
|
|
|
-
|
|
|
- for (int n = 1; n <= nmax_; n++) {
|
|
|
- Num = (z1*D1_mlxl[n] + FloatType(n))*(FloatType(n) - z1*D3_mlxl[n - 1]);
|
|
|
- Denom = (z2*D1_mlxlM1[n] + FloatType(n))*(FloatType(n) - z2*D3_mlxlM1[n - 1]);
|
|
|
- Q[l][n] = ((pow2(x[l - 1]/x[l])* Q[l][n - 1])*Num)/Denom;
|
|
|
- }
|
|
|
- // Upward recurrence for Ha and Hb - equations (7b), (8b) and (12) - (15)
|
|
|
- for (int n = 1; n <= nmax_; n++) {
|
|
|
- //Ha
|
|
|
- if ((l - 1) == pl) { // The layer below the current one is a PEC layer
|
|
|
- G1 = -D1_mlxlM1[n];
|
|
|
- G2 = -D3_mlxlM1[n];
|
|
|
- } else {
|
|
|
- G1 = (m[l]*Ha[l - 1][n - 1]) - (m[l - 1]*D1_mlxlM1[n]);
|
|
|
- G2 = (m[l]*Ha[l - 1][n - 1]) - (m[l - 1]*D3_mlxlM1[n]);
|
|
|
- } // end of if PEC
|
|
|
- Temp = Q[l][n]*G1;
|
|
|
- Num = (G2*D1_mlxl[n]) - (Temp*D3_mlxl[n]);
|
|
|
- Denom = G2 - Temp;
|
|
|
- Ha[l][n - 1] = Num/Denom;
|
|
|
- //Hb
|
|
|
- if ((l - 1) == pl) { // The layer below the current one is a PEC layer
|
|
|
- G1 = Hb[l - 1][n - 1];
|
|
|
- G2 = Hb[l - 1][n - 1];
|
|
|
- } else {
|
|
|
- G1 = (m[l - 1]*Hb[l - 1][n - 1]) - (m[l]*D1_mlxlM1[n]);
|
|
|
- G2 = (m[l - 1]*Hb[l - 1][n - 1]) - (m[l]*D3_mlxlM1[n]);
|
|
|
- } // end of if PEC
|
|
|
-
|
|
|
- Temp = Q[l][n]*G1;
|
|
|
- Num = (G2*D1_mlxl[n]) - (Temp* D3_mlxl[n]);
|
|
|
- Denom = (G2- Temp);
|
|
|
- Hb[l][n - 1] = (Num/ Denom);
|
|
|
- } // end of for Ha and Hb terms
|
|
|
- } // end of for layers iteration
|
|
|
-
|
|
|
- //**************************************//
|
|
|
- //Calculate Psi and Zeta for XL //
|
|
|
- //**************************************//
|
|
|
- // Calculate PsiXL and ZetaXL
|
|
|
- calcPsiZeta(std::complex<FloatType>(x[L - 1],0.0), PsiXL, ZetaXL);
|
|
|
-
|
|
|
-
|
|
|
- //*********************************************************************//
|
|
|
- // Finally, we calculate the scattering coefficients (an and bn) and //
|
|
|
- // the angular functions (Pi and Tau). Note that for these arrays the //
|
|
|
- // first layer is 0 (zero), in future versions all arrays will follow //
|
|
|
- // this convention to save memory. (13 Nov, 2014) //
|
|
|
- //*********************************************************************//
|
|
|
- FloatType a0=0, b0=0;
|
|
|
- for (int n = 0; n < nmax_; n++) {
|
|
|
- //********************************************************************//
|
|
|
- //Expressions for calculating an and bn coefficients are not valid if //
|
|
|
- //there is only one PEC layer (ie, for a simple PEC sphere). //
|
|
|
- //********************************************************************//
|
|
|
- if (pl < (L - 1)) {
|
|
|
- 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<FloatType>(0.0, 0.0), std::complex<FloatType>(1.0, 0.0), PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
|
|
|
- bn_[n] = PsiXL[n + 1]/ZetaXL[n + 1];
|
|
|
- }
|
|
|
- if (n == 0) {a0 = cabs(an_[0]); b0 = cabs(bn_[0]);}
|
|
|
- if (n == nmax_ - 1 && nmax_preset_ <= 0
|
|
|
- && (cabs(an_[n]) / a0 > convergence_threshold_ &&
|
|
|
- cabs(bn_[n]) / b0 > convergence_threshold_)) {
|
|
|
- std::cout << "Failed to converge in Mie series for nmax="<<nmax_ << std::endl;
|
|
|
- std::cout << "convergence threshold: "<< convergence_threshold_ << std::endl;
|
|
|
- std::cout << "Mie series a[nmax]/a[1]:" << cabs(an_[n]) / a0 << " and b[nmax]/b[1]:" << cabs(bn_[n]) / b0 << std::endl;
|
|
|
|
|
|
- }
|
|
|
+ } // end of for an and bn terms
|
|
|
+ isScaCoeffsCalc_ = true;
|
|
|
+} // end of MultiLayerMie::calcScattCoeffs()
|
|
|
+
|
|
|
+//*******************************************************************************
|
|
|
+// This function calculates the actual scattering parameters and amplitudes
|
|
|
+//
|
|
|
+// Input parameters:
|
|
|
+// L: Number of layers
|
|
|
+// pl: Index of PEC layer. If there is none just send -1
|
|
|
+// x: Array containing the size parameters of the layers [0..L-1]
|
|
|
+// m: Array containing the relative refractive indexes of the layers [0..L-1]
|
|
|
+// nTheta: Number of scattering angles
|
|
|
+// Theta: Array containing all the scattering angles where the scattering
|
|
|
+// amplitudes will be calculated
|
|
|
+// nmax_: Maximum number of multipolar expansion terms to be used for the
|
|
|
+// calculations. Only use it if you know what you are doing, otherwise
|
|
|
+// set this parameter to -1 and the function will calculate it
|
|
|
+//
|
|
|
+// Output parameters:
|
|
|
+// Qext: Efficiency factor for extinction
|
|
|
+// Qsca: Efficiency factor for scattering
|
|
|
+// Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca)
|
|
|
+// Qbk: Efficiency factor for backscattering
|
|
|
+// Qpr: Efficiency factor for the radiation pressure
|
|
|
+// g: Asymmetry factor (g = (Qext-Qpr)/Qsca)
|
|
|
+// Albedo: Single scattering albedo (Albedo = Qsca/Qext)
|
|
|
+// S1, S2: Complex scattering amplitudes
|
|
|
+//
|
|
|
+// Return value:
|
|
|
+// Number of multipolar expansion terms used for the calculations
|
|
|
+//*******************************************************************************
|
|
|
+template <typename FloatType>
|
|
|
+void MultiLayerMie<FloatType>::RunMieCalculation() {
|
|
|
+ if (size_param_.size() != refractive_index_.size())
|
|
|
+ throw std::invalid_argument(
|
|
|
+ "Each size parameter should have only one index!");
|
|
|
+ if (size_param_.size() == 0)
|
|
|
+ throw std::invalid_argument("Initialize model first!");
|
|
|
+
|
|
|
+ const std::vector<FloatType>& x = size_param_;
|
|
|
+
|
|
|
+ // MarkUncalculated();
|
|
|
+
|
|
|
+ // Calculate scattering coefficients
|
|
|
+ if (!isScaCoeffsCalc_)
|
|
|
+ calcScattCoeffs();
|
|
|
+
|
|
|
+ // Initialize the scattering parameters
|
|
|
+ Qext_ = 0.0;
|
|
|
+ Qsca_ = 0.0;
|
|
|
+ Qabs_ = 0.0;
|
|
|
+ Qbk_ = 0.0;
|
|
|
+ Qpr_ = 0.0;
|
|
|
+
|
|
|
+ asymmetry_factor_ = 0.0;
|
|
|
+ albedo_ = 0.0;
|
|
|
+
|
|
|
+ // Initialize the scattering amplitudes
|
|
|
+ std::vector<std::complex<FloatType>> tmp1(theta_.size(),
|
|
|
+ std::complex<FloatType>(0.0, 0.0));
|
|
|
+ S1_.swap(tmp1);
|
|
|
+ S2_ = S1_;
|
|
|
+ // Precalculate cos(theta) - gives about 5% speed up.
|
|
|
+ std::vector<FloatType> costheta(theta_.size(), static_cast<FloatType>(0.0));
|
|
|
+ for (unsigned int t = 0; t < theta_.size(); t++) {
|
|
|
+ costheta[t] = nmm::cos(theta_[t]);
|
|
|
+ }
|
|
|
|
|
|
-// // TODO seems to provide not enough terms for near-field calclulation.
|
|
|
-// if (cabs(an_[n]) / a0 < convergence_threshold_ &&
|
|
|
-// cabs(bn_[n]) / b0 < convergence_threshold_) {
|
|
|
-// if (nmax_preset_ <= 0) nmax_ = n;
|
|
|
-// break;
|
|
|
-// }
|
|
|
-
|
|
|
- if (nmm::isnan(an_[n].real()) || nmm::isnan(an_[n].imag()) ||
|
|
|
- nmm::isnan(bn_[n].real()) || nmm::isnan(bn_[n].imag())
|
|
|
- ) {
|
|
|
- std::cout << "nmax value was changed due to unexpected error!!! New values is "<< n
|
|
|
- << " (was "<<nmax_<<")"<<std::endl;
|
|
|
- nmax_ = n;
|
|
|
- break;
|
|
|
+ std::vector<FloatType> Pi(nmax_), Tau(nmax_);
|
|
|
+
|
|
|
+ std::complex<FloatType> Qbktmp(0.0, 0.0);
|
|
|
+ std::vector<std::complex<FloatType>> 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
|
|
|
+ for (int n = nmax_ - 2; n >= 0; n--) {
|
|
|
+ // for (int n = 0; n < nmax_; n++) {
|
|
|
+ const int n1 = n + 1;
|
|
|
+ if (mode_n_ == Modes::kAll) {
|
|
|
+ // Equation (27)
|
|
|
+ Qext_ += (n1 + n1 + 1.0) * (an_[n].real() + bn_[n].real());
|
|
|
+ // Equation (28)
|
|
|
+ Qsca_ += (n1 + n1 + 1.0) *
|
|
|
+ (an_[n].real() * an_[n].real() + an_[n].imag() * an_[n].imag() +
|
|
|
+ bn_[n].real() * bn_[n].real() + bn_[n].imag() * bn_[n].imag());
|
|
|
+ // std::cout<<"n ="<< n1 << " ext:"<<Qext_ <<"
|
|
|
+ // sca:"<<Qsca_<<std::endl;
|
|
|
+ // Equation (29)
|
|
|
+ Qpr_ += ((n1 * (n1 + 2.0) / (n1 + 1.0)) *
|
|
|
+ ((an_[n] * std::conj(an_[n1]) + bn_[n] * std::conj(bn_[n1]))
|
|
|
+ .real()) +
|
|
|
+ ((n1 + n1 + 1.0) / (n1 * (n1 + 1.0))) *
|
|
|
+ (an_[n] * std::conj(bn_[n])).real());
|
|
|
+ // Equation (33)
|
|
|
+ Qbktmp += (FloatType)(n1 + n1 + 1.0) * (1.0 - 2.0 * (n1 % 2)) *
|
|
|
+ (an_[n] - bn_[n]);
|
|
|
+ // Calculate the scattering amplitudes (S1 and S2) Equations (25a) - (25b)
|
|
|
+ for (unsigned int t = 0; t < theta_.size(); t++) {
|
|
|
+ calcPiTau(costheta[t], Pi, Tau);
|
|
|
+ S1_[t] += calc_S1(n1, an_[n], bn_[n], Pi[n], Tau[n]);
|
|
|
+ S2_[t] += calc_S2(n1, an_[n], bn_[n], Pi[n], Tau[n]);
|
|
|
}
|
|
|
-
|
|
|
- } // end of for an and bn terms
|
|
|
- isScaCoeffsCalc_ = true;
|
|
|
- } // end of MultiLayerMie::calcScattCoeffs()
|
|
|
-
|
|
|
-
|
|
|
- //**********************************************************************************//
|
|
|
- // This function calculates the actual scattering parameters and amplitudes //
|
|
|
- // //
|
|
|
- // Input parameters: //
|
|
|
- // L: Number of layers //
|
|
|
- // pl: Index of PEC layer. If there is none just send -1 //
|
|
|
- // x: Array containing the size parameters of the layers [0..L-1] //
|
|
|
- // m: Array containing the relative refractive indexes of the layers [0..L-1] //
|
|
|
- // nTheta: Number of scattering angles //
|
|
|
- // Theta: Array containing all the scattering angles where the scattering //
|
|
|
- // amplitudes will be calculated //
|
|
|
- // nmax_: Maximum number of multipolar expansion terms to be used for the //
|
|
|
- // calculations. Only use it if you know what you are doing, otherwise //
|
|
|
- // set this parameter to -1 and the function will calculate it //
|
|
|
- // //
|
|
|
- // Output parameters: //
|
|
|
- // Qext: Efficiency factor for extinction //
|
|
|
- // Qsca: Efficiency factor for scattering //
|
|
|
- // Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca) //
|
|
|
- // Qbk: Efficiency factor for backscattering //
|
|
|
- // Qpr: Efficiency factor for the radiation pressure //
|
|
|
- // g: Asymmetry factor (g = (Qext-Qpr)/Qsca) //
|
|
|
- // Albedo: Single scattering albedo (Albedo = Qsca/Qext) //
|
|
|
- // S1, S2: Complex scattering amplitudes //
|
|
|
- // //
|
|
|
- // Return value: //
|
|
|
- // Number of multipolar expansion terms used for the calculations //
|
|
|
- //**********************************************************************************//
|
|
|
- template <typename FloatType>
|
|
|
- void MultiLayerMie<FloatType>::RunMieCalculation() {
|
|
|
- if (size_param_.size() != refractive_index_.size())
|
|
|
- throw std::invalid_argument("Each size parameter should have only one index!");
|
|
|
- if (size_param_.size() == 0)
|
|
|
- throw std::invalid_argument("Initialize model first!");
|
|
|
-
|
|
|
- const std::vector<FloatType> &x = size_param_;
|
|
|
-
|
|
|
- //MarkUncalculated();
|
|
|
-
|
|
|
- // Calculate scattering coefficients
|
|
|
- if (!isScaCoeffsCalc_) calcScattCoeffs();
|
|
|
-
|
|
|
- // Initialize the scattering parameters
|
|
|
- Qext_ = 0.0;
|
|
|
- Qsca_ = 0.0;
|
|
|
- Qabs_ = 0.0;
|
|
|
- Qbk_ = 0.0;
|
|
|
- Qpr_ = 0.0;
|
|
|
-
|
|
|
- asymmetry_factor_ = 0.0;
|
|
|
- albedo_ = 0.0;
|
|
|
-
|
|
|
- // Initialize the scattering amplitudes
|
|
|
- std::vector<std::complex<FloatType> > tmp1(theta_.size(),std::complex<FloatType>(0.0, 0.0));
|
|
|
- S1_.swap(tmp1);
|
|
|
- S2_ = S1_;
|
|
|
- // Precalculate cos(theta) - gives about 5% speed up.
|
|
|
- std::vector<FloatType> costheta(theta_.size(), static_cast<FloatType>(0.0));
|
|
|
- for (unsigned int t = 0; t < theta_.size(); t++) {
|
|
|
- costheta[t] = nmm::cos(theta_[t]);
|
|
|
+ continue;
|
|
|
}
|
|
|
-
|
|
|
- std::vector<FloatType> Pi(nmax_), Tau(nmax_);
|
|
|
-
|
|
|
- std::complex<FloatType> Qbktmp(0.0, 0.0);
|
|
|
- std::vector< std::complex<FloatType> > 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
|
|
|
- for (int n = nmax_ - 2; n >= 0; n--) {
|
|
|
-// for (int n = 0; n < nmax_; n++) {
|
|
|
- const int n1 = n + 1;
|
|
|
- if (mode_n_ == Modes::kAll) {
|
|
|
- // Equation (27)
|
|
|
- Qext_ += (n1 + n1 + 1.0) * (an_[n].real() + bn_[n].real());
|
|
|
- // Equation (28)
|
|
|
- Qsca_ += (n1 + n1 + 1.0) * (an_[n].real() * an_[n].real() + an_[n].imag() * an_[n].imag()
|
|
|
- + bn_[n].real() * bn_[n].real() + bn_[n].imag() * bn_[n].imag());
|
|
|
-// std::cout<<"n ="<< n1 << " ext:"<<Qext_ <<" sca:"<<Qsca_<<std::endl;
|
|
|
- // Equation (29)
|
|
|
- Qpr_ += ((n1 * (n1 + 2.0) / (n1 + 1.0)) * ((an_[n] * std::conj(an_[n1]) + bn_[n] * std::conj(bn_[n1])).real())
|
|
|
- + ((n1 + n1 + 1.0) / (n1 * (n1 + 1.0))) * (an_[n] * std::conj(bn_[n])).real());
|
|
|
- // Equation (33)
|
|
|
- Qbktmp += (FloatType) (n1 + n1 + 1.0) * (1.0 - 2.0 * (n1 % 2)) * (an_[n] - bn_[n]);
|
|
|
- // Calculate the scattering amplitudes (S1 and S2) Equations (25a) - (25b)
|
|
|
+ if (n1 == mode_n_) {
|
|
|
+ if (mode_type_ == Modes::kElectric || mode_type_ == Modes::kAll) {
|
|
|
+ Qext_ += (n1 + n1 + 1.0) * (an_[n].real());
|
|
|
+ Qsca_ += (n1 + n1 + 1.0) * (an_[n].real() * an_[n].real() +
|
|
|
+ an_[n].imag() * an_[n].imag());
|
|
|
+ Qpr_ += std::nan("");
|
|
|
+ Qbktmp +=
|
|
|
+ (FloatType)(n1 + n1 + 1.0) * (1.0 - 2.0 * (n1 % 2)) * (an_[n]);
|
|
|
for (unsigned int t = 0; t < theta_.size(); t++) {
|
|
|
calcPiTau(costheta[t], Pi, Tau);
|
|
|
- S1_[t] += calc_S1(n1, an_[n], bn_[n], Pi[n], Tau[n]);
|
|
|
- S2_[t] += calc_S2(n1, an_[n], bn_[n], Pi[n], Tau[n]);
|
|
|
+ S1_[t] += calc_S1(n1, an_[n], static_cast<std::complex<FloatType>>(0),
|
|
|
+ Pi[n], Tau[n]);
|
|
|
+ S2_[t] += calc_S2(n1, an_[n], static_cast<std::complex<FloatType>>(0),
|
|
|
+ Pi[n], Tau[n]);
|
|
|
}
|
|
|
- continue;
|
|
|
}
|
|
|
- if (n1 == mode_n_) {
|
|
|
- if (mode_type_ == Modes::kElectric || mode_type_ == Modes::kAll) {
|
|
|
- Qext_ += (n1 + n1 + 1.0) * (an_[n].real());
|
|
|
- Qsca_ += (n1 + n1 + 1.0) * (an_[n].real() * an_[n].real() + an_[n].imag() * an_[n].imag());
|
|
|
- Qpr_ += std::nan("");
|
|
|
- Qbktmp += (FloatType) (n1 + n1 + 1.0) * (1.0 - 2.0 * (n1 % 2)) * (an_[n]);
|
|
|
- for (unsigned int t = 0; t < theta_.size(); t++) {
|
|
|
- calcPiTau(costheta[t], Pi, Tau);
|
|
|
- S1_[t] += calc_S1(n1, an_[n], static_cast<std::complex<FloatType>>(0), Pi[n], Tau[n]);
|
|
|
- S2_[t] += calc_S2(n1, an_[n], static_cast<std::complex<FloatType>>(0), Pi[n], Tau[n]);
|
|
|
- }
|
|
|
- }
|
|
|
- if (mode_type_ == Modes::kMagnetic || mode_type_ == Modes::kAll) {
|
|
|
- Qext_ += (n1 + n1 + 1.0) * (bn_[n].real());
|
|
|
- Qsca_ += (n1 + n1 + 1.0) * (bn_[n].real() * bn_[n].real() + bn_[n].imag() * bn_[n].imag());
|
|
|
- Qpr_ += std::nan("");
|
|
|
- Qbktmp += (FloatType) (n1 + n1 + 1.0) * (1.0 - 2.0 * (n1 % 2)) * (bn_[n]);
|
|
|
- for (unsigned int t = 0; t < theta_.size(); t++) {
|
|
|
- calcPiTau(costheta[t], Pi, Tau);
|
|
|
- S1_[t] += calc_S1(n1, static_cast<std::complex<FloatType>>(0), bn_[n], Pi[n], Tau[n]);
|
|
|
- S2_[t] += calc_S2(n1, static_cast<std::complex<FloatType>>(0), bn_[n], Pi[n], Tau[n]);
|
|
|
- }
|
|
|
+ if (mode_type_ == Modes::kMagnetic || mode_type_ == Modes::kAll) {
|
|
|
+ Qext_ += (n1 + n1 + 1.0) * (bn_[n].real());
|
|
|
+ Qsca_ += (n1 + n1 + 1.0) * (bn_[n].real() * bn_[n].real() +
|
|
|
+ bn_[n].imag() * bn_[n].imag());
|
|
|
+ Qpr_ += std::nan("");
|
|
|
+ Qbktmp +=
|
|
|
+ (FloatType)(n1 + n1 + 1.0) * (1.0 - 2.0 * (n1 % 2)) * (bn_[n]);
|
|
|
+ for (unsigned int t = 0; t < theta_.size(); t++) {
|
|
|
+ calcPiTau(costheta[t], Pi, Tau);
|
|
|
+ S1_[t] += calc_S1(n1, static_cast<std::complex<FloatType>>(0), bn_[n],
|
|
|
+ Pi[n], Tau[n]);
|
|
|
+ S2_[t] += calc_S2(n1, static_cast<std::complex<FloatType>>(0), bn_[n],
|
|
|
+ Pi[n], Tau[n]);
|
|
|
}
|
|
|
}
|
|
|
}
|
|
|
- FloatType x2 = pow2(x.back());
|
|
|
- Qext_ = 2.0*(Qext_)/x2; // Equation (27)
|
|
|
- Qsca_ = 2.0*(Qsca_)/x2; // Equation (28)
|
|
|
- Qpr_ = Qext_ - 4.0*(Qpr_)/x2; // Equation (29)
|
|
|
- Qabs_ = Qext_ - Qsca_; // Equation (30)
|
|
|
- albedo_ = Qsca_/Qext_; // Equation (31)
|
|
|
- asymmetry_factor_ = (Qext_ - Qpr_)/Qsca_; // Equation (32)
|
|
|
- Qbk_ = (Qbktmp.real()*Qbktmp.real() + Qbktmp.imag()*Qbktmp.imag())/x2; // Equation (33)
|
|
|
-
|
|
|
- isMieCalculated_ = true;
|
|
|
}
|
|
|
-
|
|
|
+ FloatType x2 = pow2(x.back());
|
|
|
+ Qext_ = 2.0 * (Qext_) / x2; // Equation (27)
|
|
|
+ Qsca_ = 2.0 * (Qsca_) / x2; // Equation (28)
|
|
|
+ Qpr_ = Qext_ - 4.0 * (Qpr_) / x2; // Equation (29)
|
|
|
+ Qabs_ = Qext_ - Qsca_; // Equation (30)
|
|
|
+ albedo_ = Qsca_ / Qext_; // Equation (31)
|
|
|
+ asymmetry_factor_ = (Qext_ - Qpr_) / Qsca_; // Equation (32)
|
|
|
+ Qbk_ = (Qbktmp.real() * Qbktmp.real() + Qbktmp.imag() * Qbktmp.imag()) /
|
|
|
+ x2; // Equation (33)
|
|
|
+
|
|
|
+ isMieCalculated_ = true;
|
|
|
+}
|
|
|
|
|
|
} // end of namespace nmie
|
|
|
#endif // SRC_NMIE_BASIC_HPP_
|
