From 6e1421428bc1d7c0d0209c7271f8a81607aaf78c Mon Sep 17 00:00:00 2001 From: Konstantin Ladutenko Date: Wed, 10 Jan 2024 19:37:30 +0300 Subject: [PATCH] split web interface from nmie-applied (#74) * split web interface from nmie-applied * add windows * add windows to actions * more build tests * test asser * fix assert --- .github/workflows/cmake.yml | 18 +- examples/example-get-Mie.cc | 289 ++++++++------- examples/example-minimal.cc | 72 ++-- examples/go-cc-examples.sh | 22 +- src/nmie-applied-impl.hpp | 702 ++++++++++++++++++------------------ src/nmie-applied.cc | 1 - src/nmie-applied.hpp | 429 +++++++++++++--------- src/nmie-js-wrapper.cc | 217 ++++++----- src/nmie-precision.hpp | 183 +++++----- src/nmie-web-impl.hpp | 79 ++++ src/nmie-web.hpp | 56 +++ tests/c++/go-speed-test.sh | 6 +- tests/c++/speed-test.cc | 282 ++++++++------- tests/test_py.py | 7 +- 14 files changed, 1343 insertions(+), 1020 deletions(-) create mode 100644 src/nmie-web-impl.hpp create mode 100644 src/nmie-web.hpp diff --git a/.github/workflows/cmake.yml b/.github/workflows/cmake.yml index 8a749791..4c9c0100 100644 --- a/.github/workflows/cmake.yml +++ b/.github/workflows/cmake.yml @@ -11,7 +11,7 @@ env: BUILD_TYPE: Release jobs: - macosPython: + macos_Python: runs-on: macOS-latest steps: - uses: actions/checkout@v3 @@ -23,6 +23,18 @@ jobs: working-directory: ${{github.workspace}} run: tox run + # windows_Python: + # runs-on: windows-latest + # steps: + # - uses: actions/checkout@v3 + + # - name: Install tox + # run: pip3 install tox + + # - name: Python initial test + # working-directory: ${{github.workspace}} + # run: tox run + ubuntu_Python_wo_Boost: runs-on: ubuntu-latest steps: @@ -68,6 +80,10 @@ jobs: working-directory: ${{github.workspace}}/build run: ctest -C ${{env.BUILD_TYPE}} --output-on-failure + - name: isBuilding without cmake + working-directory: ${{github.workspace}}/examples + run: ./go-cc-examples.sh + ctest_wo_Boost: runs-on: ubuntu-latest steps: diff --git a/examples/example-get-Mie.cc b/examples/example-get-Mie.cc index 8d50a845..a5be1633 100644 --- a/examples/example-get-Mie.cc +++ b/examples/example-get-Mie.cc @@ -1,49 +1,48 @@ //**********************************************************************************// -// Copyright (C) 2009-2015 Ovidio Pena // -// Copyright (C) 2013-2015 Konstantin Ladutenko // +// Copyright (C) 2009-2015 Ovidio Pena // Copyright +// (C) 2013-2015 Konstantin Ladutenko // // // -// This file is part of scattnlay // +// 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 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. // +// This program is distributed in the hope that it will be useful, // but +// WITHOUT ANY WARRANTY; without even the implied warranty of // +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU +// General Public License for more details. // // // -// The only additional remark is that we expect that all publications // -// describing work using this software, or all commercial products // -// using it, cite the following reference: // -// [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by // -// a multilayered sphere," Computer Physics Communications, // -// vol. 180, Nov. 2009, pp. 2348-2354. // +// The only additional remark is that we expect that all publications // +// describing work using this software, or all commercial products // using +// it, cite the following reference: // +// [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by // +// a multilayered sphere," Computer Physics Communications, // vol. 180, +// Nov. 2009, pp. 2348-2354. // // // -// You should have received a copy of the GNU General Public License // -// along with this program. If not, see . // +// You should have received a copy of the GNU General Public License // along +// with this program. If not, see . // //**********************************************************************************// // This program returns expansion coefficents of Mie series #include #include #include -#include "../src/nmie-applied.hpp" #include "../src/nmie-applied-impl.hpp" #include "../src/nmie-precision.hpp" #include "./read-spectra.h" // template inline T pow2(const T value) {return value*value;} -int main(int argc, char *argv[]) { - using namespace nmie ; +int main(int argc, char* argv[]) { + using namespace nmie; try { read_spectra::ReadSpectra Si_index, Ag_index; read_spectra::ReadSpectra plot_core_index_, plot_TiN_; std::string core_filename("Si-int.txt"); - //std::string core_filename("Ag.txt"); - //std::string TiN_filename("TiN.txt"); + // std::string core_filename("Ag.txt"); + // std::string TiN_filename("TiN.txt"); std::string TiN_filename("Ag-int.txt"); - //std::string TiN_filename("Si.txt"); + // std::string TiN_filename("Si.txt"); std::string shell_filename(core_filename); nmie::MultiLayerMieApplied multi_layer_mie; @@ -51,105 +50,113 @@ int main(int argc, char *argv[]) { const std::complex epsilon_Ag(-8.5014154589, 0.7585845411); const std::complex index_Si = std::sqrt(epsilon_Si); const std::complex index_Ag = std::sqrt(epsilon_Ag); - double WL=500; //nm - double core_width = 5.27; //nm Si - double inner_width = 8.22; //nm Ag - double outer_width = 67.91; //nm Si + double WL = 500; // nm + double core_width = 5.27; // nm Si + double inner_width = 8.22; // nm Ag + double outer_width = 67.91; // nm Si bool isSiAgSi = true; double delta_width = 25.0; - //bool isSiAgSi = false; + // bool isSiAgSi = false; if (isSiAgSi) { - core_width = 5.27; //nm Si - inner_width = 8.22; //nm Ag - outer_width = 67.91; //nm Si + core_width = 5.27; // nm Si + inner_width = 8.22; // nm Ag + outer_width = 67.91; // nm Si multi_layer_mie.AddTargetLayer(core_width, index_Si); multi_layer_mie.AddTargetLayer(inner_width, index_Ag); - multi_layer_mie.AddTargetLayer(outer_width+delta_width, index_Si); + multi_layer_mie.AddTargetLayer(outer_width + delta_width, index_Si); } else { - inner_width = 31.93; //nm Ag - outer_width = 4.06; //nm Si + inner_width = 31.93; // nm Ag + outer_width = 4.06; // nm Si multi_layer_mie.AddTargetLayer(inner_width, index_Ag); - multi_layer_mie.AddTargetLayer(outer_width+delta_width, index_Si); + multi_layer_mie.AddTargetLayer(outer_width + delta_width, index_Si); } - for (int dd = 0; dd<50; ++dd) { + for (int dd = 0; dd < 50; ++dd) { delta_width = dd; - FILE *fp; - std::string fname = "absorb-layered-spectra-d"+std::to_string(dd)+".dat"; - fp = fopen(fname.c_str(), "w"); + FILE* fp; + std::string fname = + "absorb-layered-spectra-d" + std::to_string(dd) + ".dat"; + fp = fopen(fname.c_str(), "w"); - multi_layer_mie.SetWavelength(WL); - multi_layer_mie.RunMieCalculation(); - - double Qabs = static_cast(multi_layer_mie.GetQabs()); - printf("Qabs = %g\n", Qabs); - std::vector< std::vector > > aln, bln, cln, dln; - multi_layer_mie.GetExpanCoeffs(aln, bln, cln, dln); - std::vector< std::vector > > d_aln = - nmie::ConvertComplexVectorVector(aln); - std::string str = std::string("#WL "); - for (int l = 0; l &Si = Si_data[i].second; - const std::complex &Ag = Ag_data[i].second; - str+=std::to_string(WL); - multi_layer_mie.ClearTarget(); - if (isSiAgSi) { - multi_layer_mie.AddTargetLayer(core_width, Si); - multi_layer_mie.AddTargetLayer(inner_width, Ag); - multi_layer_mie.AddTargetLayer(outer_width+delta_width, Si); - } else { - inner_width = 31.93; //nm Ag - outer_width = 4.06; //nm Si - multi_layer_mie.AddTargetLayer(inner_width, Ag); - multi_layer_mie.AddTargetLayer(outer_width+delta_width, Si); - } multi_layer_mie.SetWavelength(WL); multi_layer_mie.RunMieCalculation(); - multi_layer_mie.GetQabs(); - multi_layer_mie.GetExpanCoeffs(aln, bln, cln, dln); - for (int l = 0; l(pow2(std::abs(aln[l][n]))+ - pow2(std::abs(dln[l][n])))) - + " " - + std::to_string(static_cast(pow2(std::abs(bln[l][n])) - + pow2(std::abs(cln[l][n])) )); - // str+=" "+std::to_string(aln[l][n].real() - pow2(std::abs(aln[l][n])) - // +dln[l][n].real() - pow2(std::abs(dln[l][n]))) - // + " " - // + std::to_string(bln[l][n].real() - pow2(std::abs(bln[l][n])) - // +cln[l][n].real() - pow2(std::abs(cln[l][n])) ); - } + double Qabs = static_cast(multi_layer_mie.GetQabs()); + printf("Qabs = %g\n", Qabs); + std::vector > > aln, bln, cln, + dln; + multi_layer_mie.GetExpanCoeffs(aln, bln, cln, dln); + std::vector > > d_aln = + nmie::ConvertComplexVectorVector(aln); + std::string str = std::string("#WL "); + for (int l = 0; l < d_aln.size(); ++l) { + for (int n = 0; n < 3; ++n) { + str += "|a|^2+|d|^2_ln" + std::to_string(l) + std::to_string(n) + + " " + "|b|^2+|c|^2_ln" + std::to_string(l) + + std::to_string(n) + " "; + } } - str+="\n"; + str += "\n"; fprintf(fp, "%s", str.c_str()); + fprintf(fp, "# |a|+|d|"); str.clear(); - } - fclose(fp); + double from_WL = 400; + double to_WL = 600; + int total_points = 401; + double delta_WL = std::abs(to_WL - from_WL) / (total_points - 1.0); + Si_index.ReadFromFile(core_filename) + .ResizeToComplex(from_WL, to_WL, total_points) + .ToIndex(); + Ag_index.ReadFromFile(TiN_filename) + .ResizeToComplex(from_WL, to_WL, total_points) + .ToIndex(); + auto Si_data = Si_index.GetIndex(); + auto Ag_data = Ag_index.GetIndex(); + for (int i = 0; i < Si_data.size(); ++i) { + const double& WL = Si_data[i].first; + const std::complex& Si = Si_data[i].second; + const std::complex& Ag = Ag_data[i].second; + str += std::to_string(WL); + multi_layer_mie.ClearTarget(); + if (isSiAgSi) { + multi_layer_mie.AddTargetLayer(core_width, Si); + multi_layer_mie.AddTargetLayer(inner_width, Ag); + multi_layer_mie.AddTargetLayer(outer_width + delta_width, Si); + } else { + inner_width = 31.93; // nm Ag + outer_width = 4.06; // nm Si + multi_layer_mie.AddTargetLayer(inner_width, Ag); + multi_layer_mie.AddTargetLayer(outer_width + delta_width, Si); + } + multi_layer_mie.SetWavelength(WL); + multi_layer_mie.RunMieCalculation(); + multi_layer_mie.GetQabs(); + multi_layer_mie.GetExpanCoeffs(aln, bln, cln, dln); + for (int l = 0; l < aln.size(); ++l) { + for (int n = 0; n < 3; ++n) { + str += " " + + std::to_string(static_cast( + pow2(std::abs(aln[l][n])) + pow2(std::abs(dln[l][n])))) + + " " + + std::to_string(static_cast( + pow2(std::abs(bln[l][n])) + pow2(std::abs(cln[l][n])))); + + // str+=" "+std::to_string(aln[l][n].real() - + // pow2(std::abs(aln[l][n])) +dln[l][n].real() - + // pow2(std::abs(dln[l][n]))) + // + " " + // + std::to_string(bln[l][n].real() - pow2(std::abs(bln[l][n])) + // +cln[l][n].real() - pow2(std::abs(cln[l][n])) + // ); + } + } + str += "\n"; + fprintf(fp, "%s", str.c_str()); + str.clear(); + } + + fclose(fp); } // WL = 500; @@ -161,49 +168,57 @@ int main(int argc, char *argv[]) { // printf("\n Scattering"); // for (int l = 0; l // -// Copyright (C) 2013-2015 Konstantin Ladutenko // +// Copyright (C) 2009-2015 Ovidio Pena // Copyright +// (C) 2013-2015 Konstantin Ladutenko // // // -// This file is part of scattnlay // +// 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 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. // +// This program is distributed in the hope that it will be useful, // but +// WITHOUT ANY WARRANTY; without even the implied warranty of // +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU +// General Public License for more details. // // // -// The only additional remark is that we expect that all publications // -// describing work using this software, or all commercial products // -// using it, cite the following reference: // -// [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by // -// a multilayered sphere," Computer Physics Communications, // -// vol. 180, Nov. 2009, pp. 2348-2354. // +// The only additional remark is that we expect that all publications // +// describing work using this software, or all commercial products // using +// it, cite the following reference: // +// [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by // +// a multilayered sphere," Computer Physics Communications, // vol. 180, +// Nov. 2009, pp. 2348-2354. // // // -// You should have received a copy of the GNU General Public License // -// along with this program. If not, see . // +// You should have received a copy of the GNU General Public License // along +// with this program. If not, see . // //**********************************************************************************// // This program evaluates absorption of a triple layered nanoparticle +#include #include #include +#include +#include #include -#include "../src/nmie-applied.hpp" + #include "../src/nmie-applied-impl.hpp" -int main(int , char **) { +int main(int, char**) { try { nmie::MultiLayerMieApplied multi_layer_mie; const std::complex epsilon_Si(18.4631066585, 0.6259727805); const std::complex epsilon_Ag(-8.5014154589, 0.7585845411); const std::complex index_Si = std::sqrt(epsilon_Si); const std::complex index_Ag = std::sqrt(epsilon_Ag); - double WL=500; //nm - double core_width = 5.27; //nm Si - double inner_width = 8.22; //nm Ag - double outer_width = 67.91; //nm Si - core_width = 5.27; //nm Si - inner_width = 8.22; //nm Ag - outer_width = 67.91; //nm Si + double WL = 500; // nm + double core_width = 5.27; // nm Si + double inner_width = 8.22; // nm Ag + double outer_width = 67.91; // nm Si + core_width = 5.27; // nm Si + inner_width = 8.22; // nm Ag + outer_width = 67.91; // nm Si multi_layer_mie.AddTargetLayer(core_width, index_Si); multi_layer_mie.AddTargetLayer(inner_width, index_Ag); multi_layer_mie.AddTargetLayer(outer_width, index_Si); multi_layer_mie.SetWavelength(WL); multi_layer_mie.RunMieCalculation(); double Qabs = multi_layer_mie.GetQabs(); - printf("Qabs = %g\n", Qabs); - } catch( const std::invalid_argument &ia ) { + std::stringstream stream; + stream << std::fixed << std::setprecision(10) << Qabs; + auto Qabs_str = stream.str(); + printf("Qabs = %s\n", Qabs_str.c_str()); + assert(Qabs_str == "3.1415556911"); + + } catch (const std::invalid_argument& ia) { // Will catch if multi_layer_mie fails or other errors. std::cerr << "Invalid argument: " << ia.what() << std::endl; return -1; } - return 0; + return 0; } - - diff --git a/examples/go-cc-examples.sh b/examples/go-cc-examples.sh index 549a0693..2a225f2a 100755 --- a/examples/go-cc-examples.sh +++ b/examples/go-cc-examples.sh @@ -10,23 +10,23 @@ PROGRAM='scattnlay-example.bin' # echo Compilation done. Running... # time ./$PROGRAM -file=test-surf-integral.cc -echo Compile $file with gcc -rm -f $PROGRAM -g++ -Ofast -std=c++11 $file ../src/shell-generator.cc -lm -lrt -o $PROGRAM -march=native -mtune=native -msse4.2 -# g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc ../src/shell-generator.cc -lm -lrt -o $PROGRAM -march=native -mtune=native -msse4.2 -echo Compilation done. Running... -./$PROGRAM - - -# file=example-minimal.cc +# file=test-surf-integral.cc # echo Compile $file with gcc # rm -f $PROGRAM -# g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc -lm -lrt -o $PROGRAM -march=native -mtune=native -msse4.2 +# g++ -Ofast -std=c++11 $file ../src/shell-generator.cc -lm -lrt -o $PROGRAM -march=native -mtune=native -msse4.2 +# # g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc ../src/shell-generator.cc -lm -lrt -o $PROGRAM -march=native -mtune=native -msse4.2 # echo Compilation done. Running... # ./$PROGRAM +file=example-minimal.cc +echo Compile $file with gcc +rm -f $PROGRAM +g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc -lm -o $PROGRAM -march=native -mtune=native -msse4.2 +echo Compilation done. Running... +./$PROGRAM + + # file=example-get-Mie.cc # echo Compile $file with gcc # rm -f $PROGRAM diff --git a/src/nmie-applied-impl.hpp b/src/nmie-applied-impl.hpp index 37303d74..c0dbbc39 100644 --- a/src/nmie-applied-impl.hpp +++ b/src/nmie-applied-impl.hpp @@ -1,369 +1,389 @@ #ifndef SRC_NMIE_APPLIED_IMPL_HPP_ #define SRC_NMIE_APPLIED_IMPL_HPP_ -//**********************************************************************************// -// Copyright (C) 2009-2018 Ovidio Pena // -// Copyright (C) 2013-2018 Konstantin Ladutenko // -// // -// This file is part of scattnlay // -// // -// This program is free software: you can redistribute it and/or modify // -// it under the terms of the GNU General Public License as published by // -// the Free Software Foundation, either version 3 of the License, or // -// (at your option) any later version. // -// // -// This program is distributed in the hope that it will be useful, // -// but WITHOUT ANY WARRANTY; without even the implied warranty of // -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // -// GNU General Public License for more details. // -// // -// The only additional remark is that we expect that all publications // -// describing work using this software, or all commercial products // -// using it, cite 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 . // -// // -// @brief Wrapper class around nMie function for ease of use // -// // -//**********************************************************************************// -#include +//********************************************************************************** +// Copyright (C) 2009-2018 Ovidio Pena +// Copyright (C) 2013-2018 Konstantin Ladutenko +// +// This file is part of scattnlay +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// The only additional remark is that we expect that all publications +// describing work using this software, or all commercial products +// using it, cite 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 . +// +// @brief Wrapper class around nMie function for ease of use +// +//********************************************************************************** #include +#include #include #include #include #include -namespace nmie { - // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // - template - void MultiLayerMieApplied::GetFailed() { - FloatType faild_x = 9.42477796076938; - //FloatType faild_x = 9.42477796076937; - std::complex z(faild_x, 0.0); - std::vector nmax_local_array = {20, 100, 500, 2500}; - for (auto nmax_local : nmax_local_array) { - std::vector > D1_failed(nmax_local + 1); - // Downward recurrence for D1 - equations (16a) and (16b) - D1_failed[nmax_local] = std::complex(0.0, 0.0); - const std::complex zinv = std::complex(1.0, 0.0)/z; - for (int n = nmax_local; n > 0; n--) { - D1_failed[n - 1] = FloatType(n)*zinv - 1.0/(D1_failed[n] + FloatType(n)*zinv); - } - printf("Faild D1[0] from reccurence (z = %16.14f, nmax = %d): %g\n", - faild_x, nmax_local, D1_failed[0].real()); - } - printf("Faild D1[0] from continued fraction (z = %16.14f): %g\n", faild_x, - calcD1confra(0,z).real()); - //D1[nmax_] = calcD1confra(nmax_, z); - +#include "nmie-applied.hpp" - } - // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // - template - void MultiLayerMieApplied::AddTargetLayer(FloatType width, std::complex layer_index) { - this->MarkUncalculated(); - if (width <= 0) - throw std::invalid_argument("Layer width should be positive!"); - target_width_.push_back(width); - target_index_.push_back(layer_index); - } // end of void MultiLayerMieApplied::AddTargetLayer(...) - // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // - template - void MultiLayerMieApplied::AddTargetLayerReIm(FloatType width, - FloatType re_layer_index, FloatType im_layer_index) { - this->MarkUncalculated(); - if (width <= 0) - throw std::invalid_argument("Layer width should be positive!"); - target_width_.push_back(width); - target_index_.push_back(std::complex(re_layer_index,im_layer_index)); - } // end of void MultiLayerMieApplied::AddTargetLayer(...) - // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // - template - void MultiLayerMieApplied::SetTargetPEC(FloatType radius) { - this->MarkUncalculated(); - if (target_width_.size() != 0 || target_index_.size() != 0) - throw std::invalid_argument("Error! Define PEC target radius before any other layers!"); - // Add layer of any index... - AddTargetLayer(radius, std::complex(0.0, 0.0)); - // ... and mark it as PEC - this->SetPECLayer(0); - } - // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // - template - void MultiLayerMieApplied::SetCoatingIndex(std::vector > index) { - this->MarkUncalculated(); - coating_index_.clear(); - for (auto value : index) coating_index_.push_back(value); - } // end of void MultiLayerMieApplied::SetCoatingIndex(std::vector index); - // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // - template - void MultiLayerMieApplied::SetCoatingWidth(std::vector width) { - this->MarkUncalculated(); - coating_width_.clear(); - for (auto w : width) - if (w <= 0) - throw std::invalid_argument("Coating width should be positive!"); - else coating_width_.push_back(w); - } - // end of void MultiLayerMieApplied::SetCoatingWidth(...); - // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // - template - void MultiLayerMieApplied::SetWidthSP(const std::vector& size_parameter) { - this->MarkUncalculated(); - this->size_param_.clear(); - FloatType prev_size_parameter = 0.0; - for (auto layer_size_parameter : size_parameter) { - if (layer_size_parameter <= 0.0) - throw std::invalid_argument("Size parameter should be positive!"); - if (prev_size_parameter > layer_size_parameter) - throw std::invalid_argument - ("Size parameter for next layer should be larger than the previous one!"); - prev_size_parameter = layer_size_parameter; - this->size_param_.push_back(layer_size_parameter); - } - } - // end of void MultiLayerMieApplied::SetWidthSP(...); - // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // - template - void MultiLayerMieApplied::SetIndexSP(const std::vector< std::complex >& index) { - this->MarkUncalculated(); - //refractive_index_.clear(); - this->refractive_index_ = index; - // for (auto value : index) refractive_index_.push_back(value); - } // end of void MultiLayerMieApplied::SetIndexSP(...); - // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // - template - void MultiLayerMieApplied::SetFieldPointsSP(const std::vector< std::vector >& coords_sp) { - if (coords_sp.size() != 3) - throw std::invalid_argument("Error! Wrong dimension of field monitor points!"); - if (coords_sp[0].size() != coords_sp[1].size() || coords_sp[0].size() != coords_sp[2].size()) - throw std::invalid_argument("Error! Missing coordinates for field monitor points!"); - this->coords_ = coords_sp; - // for (int i = 0; i < coords_sp_[0].size(); ++i) { - // printf("%g, %g, %g\n", coords_sp_[0][i], coords_sp_[1][i], coords_sp_[2][i]); - // } - } // end of void MultiLayerMieApplied::SetFieldPointsSP(...) - // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // - template - void MultiLayerMieApplied::GenerateSizeParameter() { - this->MarkUncalculated(); - this->size_param_.clear(); - FloatType radius = 0.0; - for (auto width : target_width_) { - radius += width; - this->size_param_.push_back(2*nmie::PI_*radius/wavelength_); - } - for (auto width : coating_width_) { - radius += width; - this->size_param_.push_back(2*nmie::PI_*radius/wavelength_); +namespace nmie { +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // +template +void MultiLayerMieApplied::GetFailed() { + FloatType faild_x = 9.42477796076938; + // FloatType faild_x = 9.42477796076937; + std::complex z(faild_x, 0.0); + std::vector nmax_local_array = {20, 100, 500, 2500}; + for (auto nmax_local : nmax_local_array) { + std::vector > D1_failed(nmax_local + 1); + // Downward recurrence for D1 - equations (16a) and (16b) + D1_failed[nmax_local] = std::complex(0.0, 0.0); + const std::complex zinv = std::complex(1.0, 0.0) / z; + for (int n = nmax_local; n > 0; n--) { + D1_failed[n - 1] = + FloatType(n) * zinv - 1.0 / (D1_failed[n] + FloatType(n) * zinv); } - this->total_radius_ = radius; - } // end of void MultiLayerMieApplied::GenerateSizeParameter(); - // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // - template - void MultiLayerMieApplied::GenerateIndex() { - this->MarkUncalculated(); - this->refractive_index_.clear(); - for (auto index : this->target_index_) - this->refractive_index_.push_back(index); - for (auto index : this->coating_index_) - this->refractive_index_.push_back(index); - } // end of void MultiLayerMieApplied::GenerateIndex(); - // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // - template - FloatType MultiLayerMieApplied::GetTotalRadius() { - if (!this->isMieCalculated()) GenerateSizeParameter(); - return this->total_radius_; - } // end of FloatType MultiLayerMieApplied::GetTotalRadius(); - // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // - template - std::vector< std::vector > - MultiLayerMieApplied::GetSpectra(FloatType from_WL, FloatType to_WL, int samples) { - if (!this->isMieCalculated()) - throw std::invalid_argument("You should run calculations before result request!"); - std::vector< std::vector > spectra; - FloatType step_WL = (to_WL - from_WL)/static_cast(samples); - FloatType wavelength_backup = wavelength_; - long fails = 0; - for (FloatType WL = from_WL; WL < to_WL; WL += step_WL) { - wavelength_ = WL; - try { - RunMieCalculation(); - } catch(const std::invalid_argument& ia) { - fails++; - continue; - } - //printf("%3.1f ",WL); - spectra.push_back(std::vector({wavelength_, this->GetQext(), - this->GetQsca(), this->GetQabs(), this->GetQbk()})); - } // end of for each WL in spectra - printf("Spectrum has %li fails\n",fails); - wavelength_ = wavelength_backup; - return spectra; - } - // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // - template - void MultiLayerMieApplied::ClearTarget() { - this->MarkUncalculated(); - this->target_width_.clear(); - this->target_index_.clear(); + printf("Faild D1[0] from reccurence (z = %16.14f, nmax = %d): %g\n", + faild_x, nmax_local, D1_failed[0].real()); } - // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // - template - void MultiLayerMieApplied::ClearCoating() { - this->MarkUncalculated(); - this->coating_width_.clear(); - this->coating_index_.clear(); + printf("Faild D1[0] from continued fraction (z = %16.14f): %g\n", faild_x, + calcD1confra(0, z).real()); + // D1[nmax_] = calcD1confra(nmax_, z); +} +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // +template +void MultiLayerMieApplied::AddTargetLayer( + FloatType width, + std::complex layer_index) { + this->MarkUncalculated(); + if (width <= 0) + throw std::invalid_argument("Layer width should be positive!"); + target_width_.push_back(width); + target_index_.push_back(layer_index); +} // end of void MultiLayerMieApplied::AddTargetLayer(...) +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // +template +void MultiLayerMieApplied::AddTargetLayerReIm( + FloatType width, + FloatType re_layer_index, + FloatType im_layer_index) { + this->MarkUncalculated(); + if (width <= 0) + throw std::invalid_argument("Layer width should be positive!"); + target_width_.push_back(width); + target_index_.push_back( + std::complex(re_layer_index, im_layer_index)); +} // end of void MultiLayerMieApplied::AddTargetLayer(...) +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // +template +void MultiLayerMieApplied::SetTargetPEC(FloatType radius) { + this->MarkUncalculated(); + if (target_width_.size() != 0 || target_index_.size() != 0) + throw std::invalid_argument( + "Error! Define PEC target radius before any other layers!"); + // Add layer of any index... + AddTargetLayer(radius, std::complex(0.0, 0.0)); + // ... and mark it as PEC + this->SetPECLayer(0); +} +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // +template +void MultiLayerMieApplied::SetCoatingIndex( + std::vector > index) { + this->MarkUncalculated(); + coating_index_.clear(); + for (auto value : index) + coating_index_.push_back(value); +} // end of void + // MultiLayerMieApplied::SetCoatingIndex(std::vector + // index); +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // +template +void MultiLayerMieApplied::SetCoatingWidth( + std::vector width) { + this->MarkUncalculated(); + coating_width_.clear(); + for (auto w : width) + if (w <= 0) + throw std::invalid_argument("Coating width should be positive!"); + else + coating_width_.push_back(w); +} +// end of void MultiLayerMieApplied::SetCoatingWidth(...); +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // +template +void MultiLayerMieApplied::SetWidthSP( + const std::vector& size_parameter) { + this->MarkUncalculated(); + this->size_param_.clear(); + FloatType prev_size_parameter = 0.0; + for (auto layer_size_parameter : size_parameter) { + if (layer_size_parameter <= 0.0) + throw std::invalid_argument("Size parameter should be positive!"); + if (prev_size_parameter > layer_size_parameter) + throw std::invalid_argument( + "Size parameter for next layer should be larger than the previous " + "one!"); + prev_size_parameter = layer_size_parameter; + this->size_param_.push_back(layer_size_parameter); } - // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // - template - void MultiLayerMieApplied::ClearLayers() { - this->MarkUncalculated(); - this->ClearTarget(); - this->ClearCoating(); +} +// end of void MultiLayerMieApplied::SetWidthSP(...); +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // +template +void MultiLayerMieApplied::SetIndexSP( + const std::vector >& index) { + this->MarkUncalculated(); + // refractive_index_.clear(); + this->refractive_index_ = index; + // for (auto value : index) refractive_index_.push_back(value); +} // end of void MultiLayerMieApplied::SetIndexSP(...); +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // +template +void MultiLayerMieApplied::SetFieldPointsSP( + const std::vector >& coords_sp) { + if (coords_sp.size() != 3) + throw std::invalid_argument( + "Error! Wrong dimension of field monitor points!"); + if (coords_sp[0].size() != coords_sp[1].size() || + coords_sp[0].size() != coords_sp[2].size()) + throw std::invalid_argument( + "Error! Missing coordinates for field monitor points!"); + this->coords_ = coords_sp; + // for (int i = 0; i < coords_sp_[0].size(); ++i) { + // printf("%g, %g, %g\n", coords_sp_[0][i], coords_sp_[1][i], + // coords_sp_[2][i]); + // } +} // end of void MultiLayerMieApplied::SetFieldPointsSP(...) +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // +template +void MultiLayerMieApplied::GenerateSizeParameter() { + this->MarkUncalculated(); + this->size_param_.clear(); + FloatType radius = 0.0; + for (auto width : target_width_) { + radius += width; + this->size_param_.push_back(2 * nmie::PI_ * radius / wavelength_); } - // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // - template - void MultiLayerMieApplied::ClearAllDesign() { - this->MarkUncalculated(); - this->ClearLayers(); - this->size_param_.clear(); - this->refractive_index_.clear(); + for (auto width : coating_width_) { + radius += width; + this->size_param_.push_back(2 * nmie::PI_ * radius / wavelength_); } - // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // - // Computational core - // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // - template - void MultiLayerMieApplied::ConvertToSP() { - this->MarkUncalculated(); - if (target_width_.size() + coating_width_.size() == 0) - return; // Nothing to convert, we suppose that SP was set directly + this->total_radius_ = radius; +} // end of void MultiLayerMieApplied::GenerateSizeParameter(); +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // +template +void MultiLayerMieApplied::GenerateIndex() { + this->MarkUncalculated(); + this->refractive_index_.clear(); + for (auto index : this->target_index_) + this->refractive_index_.push_back(index); + for (auto index : this->coating_index_) + this->refractive_index_.push_back(index); +} // end of void MultiLayerMieApplied::GenerateIndex(); +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // +template +FloatType MultiLayerMieApplied::GetTotalRadius() { + if (!this->isMieCalculated()) GenerateSizeParameter(); - GenerateIndex(); - if (this->size_param_.size() != this->refractive_index_.size()) - throw std::invalid_argument("Ivalid conversion of width to size parameter units!/n"); - } - // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // - template - void MultiLayerMieApplied::RunMieCalculation() { - ConvertToSP(); - this->MultiLayerMie::RunMieCalculation(); - } - - template - void MultiLayerMieApplied::RunFieldCalculationPolar(const int outer_arc_points, - const int radius_points, - const double from_Rho, const double to_Rho, - const double from_Theta, const double to_Theta, - const double from_Phi, const double to_Phi, - const int isIgnoreAvailableNmax) { - ConvertToSP(); // Converts to size parameter units only the particle design, - // so we need to convert input parameters too... - const FloatType a = 2*nmie::PI_/wavelength_; - this->MultiLayerMie::RunFieldCalculationPolar(outer_arc_points, radius_points, a*from_Rho, a*to_Rho, - from_Theta, to_Theta, from_Phi, to_Phi, - isIgnoreAvailableNmax == 0 ? false : true); - } - + return this->total_radius_; +} // end of FloatType MultiLayerMieApplied::GetTotalRadius(); +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // +template +std::vector > +MultiLayerMieApplied::GetSpectra(FloatType from_WL, + FloatType to_WL, + int samples) { + if (!this->isMieCalculated()) + throw std::invalid_argument( + "You should run calculations before result request!"); + std::vector > spectra; + FloatType step_WL = (to_WL - from_WL) / static_cast(samples); + FloatType wavelength_backup = wavelength_; + long fails = 0; + for (FloatType WL = from_WL; WL < to_WL; WL += step_WL) { + wavelength_ = WL; + try { + RunMieCalculation(); + } catch (const std::invalid_argument& ia) { + fails++; + continue; + } + // printf("%3.1f ",WL); + spectra.push_back( + std::vector({wavelength_, this->GetQext(), this->GetQsca(), + this->GetQabs(), this->GetQbk()})); + } // end of for each WL in spectra + printf("Spectrum has %li fails\n", fails); + wavelength_ = wavelength_backup; + return spectra; +} +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // template -void MultiLayerMieApplied::RunFieldCalculationCartesian(const int first_side_points, - const int second_side_points, - const double relative_side_length, - const int plane_selected, - const double at_x, const double at_y, - const double at_z, - const int isIgnoreAvailableNmax) { -// std::cout<<'test'<::ClearTarget() { + this->MarkUncalculated(); + this->target_width_.clear(); + this->target_index_.clear(); +} +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // +template +void MultiLayerMieApplied::ClearCoating() { + this->MarkUncalculated(); + this->coating_width_.clear(); + this->coating_index_.clear(); +} +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // +template +void MultiLayerMieApplied::ClearLayers() { + this->MarkUncalculated(); + this->ClearTarget(); + this->ClearCoating(); +} +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // +template +void MultiLayerMieApplied::ClearAllDesign() { + this->MarkUncalculated(); + this->ClearLayers(); + this->size_param_.clear(); + this->refractive_index_.clear(); +} +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // +// Computational core +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // +template +void MultiLayerMieApplied::ConvertToSP() { + this->MarkUncalculated(); + if (target_width_.size() + coating_width_.size() == 0) + return; // Nothing to convert, we suppose that SP was set directly + GenerateSizeParameter(); + GenerateIndex(); + if (this->size_param_.size() != this->refractive_index_.size()) + throw std::invalid_argument( + "Ivalid conversion of width to size parameter units!/n"); +} +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // +template +void MultiLayerMieApplied::RunMieCalculation() { ConvertToSP(); - this->MultiLayerMie::RunFieldCalculationCartesian( first_side_points, second_side_points, - relative_side_length, plane_selected, - at_x, at_y, at_z, - isIgnoreAvailableNmax == 0 ? false : true); - + this->MultiLayerMie::RunMieCalculation(); } -//from https://toughengineer.github.io/demo/dsp/fft-perf/ -template -emscripten::val toJSFloat64Array(const std::vector &v) { - emscripten::val view{ emscripten::typed_memory_view(v.size(), v.data()) }; // make a view of local object - - auto result = emscripten::val::global("Float64Array").new_(v.size()); // make a JS typed array - result.call("set", view); // set typed array values "on the JS side" using the memory view - - return result; +template +void MultiLayerMieApplied::RunFieldCalculationPolar( + const int outer_arc_points, + const int radius_points, + const double from_Rho, + const double to_Rho, + const double from_Theta, + const double to_Theta, + const double from_Phi, + const double to_Phi, + const int isIgnoreAvailableNmax) { + ConvertToSP(); // Converts to size parameter units only the particle design, + // so we need to convert input parameters too... + const FloatType a = 2 * nmie::PI_ / wavelength_; + this->MultiLayerMie::RunFieldCalculationPolar( + outer_arc_points, radius_points, a * from_Rho, a * to_Rho, from_Theta, + to_Theta, from_Phi, to_Phi, isIgnoreAvailableNmax == 0 ? false : true); } template -emscripten::val MultiLayerMieApplied::GetFieldEabs() { - auto Eabs = this->MultiLayerMie::GetFieldEabs(); - return toJSFloat64Array(Eabs); +void MultiLayerMieApplied::RunFieldCalculationCartesian( + const int first_side_points, + const int second_side_points, + const double relative_side_length, + const int plane_selected, + const double at_x, + const double at_y, + const double at_z, + const int isIgnoreAvailableNmax) { + // std::cout<<'test'<MultiLayerMie::RunFieldCalculationCartesian( + first_side_points, second_side_points, relative_side_length, + plane_selected, at_x, at_y, at_z, + isIgnoreAvailableNmax == 0 ? false : true); } + +// ********************************************************************** // // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // - template - void MultiLayerMieApplied::GetExpanCoeffs( std::vector< std::vector > >& aln, std::vector< std::vector > >& bln, std::vector< std::vector > >& cln, std::vector< std::vector > >& dln) { - ConvertToSP(); // Case of call before running full Mie calculation. - // Calculate scattering coefficients an_ and bn_ - this->calcScattCoeffs(); - // Calculate expansion coefficients aln_, bln_, cln_, and dln_ - this->calcExpanCoeffs(); - aln = this->aln_; - bln = this->bln_; - cln = this->cln_; - dln = this->dln_; +// ********************************************************************** // +template +void MultiLayerMieApplied::GetExpanCoeffs( + std::vector > >& aln, + std::vector > >& bln, + std::vector > >& cln, + std::vector > >& dln) { + ConvertToSP(); // Case of call before running full Mie calculation. + // Calculate scattering coefficients an_ and bn_ + this->calcScattCoeffs(); + // Calculate expansion coefficients aln_, bln_, cln_, and dln_ + this->calcExpanCoeffs(); + aln = this->aln_; + bln = this->bln_; + cln = this->cln_; + dln = this->dln_; - } // end of void MultiLayerMieApplied::GetExpanCoeffs( ...) - // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // +} // end of void MultiLayerMieApplied::GetExpanCoeffs( ...) +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // } // end of namespace nmie #endif // SRC_NMIE_APPLIED_IMPL_HPP_ diff --git a/src/nmie-applied.cc b/src/nmie-applied.cc index 2859ecf6..f8c0204b 100644 --- a/src/nmie-applied.cc +++ b/src/nmie-applied.cc @@ -34,7 +34,6 @@ #include #include -#include "nmie-applied.hpp" #include "nmie-applied-impl.hpp" namespace nmie { diff --git a/src/nmie-applied.hpp b/src/nmie-applied.hpp index 0fd2b8de..40ffc4cf 100644 --- a/src/nmie-applied.hpp +++ b/src/nmie-applied.hpp @@ -1,35 +1,35 @@ #ifndef SRC_NMIE_APPLIED_HPP_ #define SRC_NMIE_APPLIED_HPP_ -//**********************************************************************************// -// Copyright (C) 2009-2018 Ovidio Pena // -// Copyright (C) 2013-2018 Konstantin Ladutenko // -// // -// This file is part of scattnlay // -// // -// This program is free software: you can redistribute it and/or modify // -// it under the terms of the GNU General Public License as published by // -// the Free Software Foundation, either version 3 of the License, or // -// (at your option) any later version. // -// // -// This program is distributed in the hope that it will be useful, // -// but WITHOUT ANY WARRANTY; without even the implied warranty of // -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // -// GNU General Public License for more details. // -// // -// The only additional remark is that we expect that all publications // -// describing work using this software, or all commercial products // -// using it, cite 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 . // -//**********************************************************************************// +//********************************************************************************** +// Copyright (C) 2009-2018 Ovidio Pena +// Copyright (C) 2013-2018 Konstantin Ladutenko +// +// This file is part of scattnlay +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// The only additional remark is that we expect that all publications +// describing work using this software, or all commercial products +// using it, cite 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 . +//********************************************************************************** #include #include @@ -37,161 +37,232 @@ #include #include -#include "nmie.hpp" #include "nmie-basic.hpp" #include "nmie-nearfield.hpp" +#include "nmie.hpp" -#include -#include +namespace nmie { +int nMieApplied(const unsigned int L, + const int pl, + std::vector& x, + std::vector >& m, + const unsigned int nTheta, + std::vector& Theta, + const int nmax, + double* Qext, + double* Qsca, + double* Qabs, + double* Qbk, + double* Qpr, + double* g, + double* Albedo, + std::vector >& S1, + std::vector >& S2); +int nMieApplied(const unsigned int L, + std::vector& x, + std::vector >& m, + const unsigned int nTheta, + std::vector& Theta, + double* Qext, + double* Qsca, + double* Qabs, + double* Qbk, + double* Qpr, + double* g, + double* Albedo, + std::vector >& S1, + std::vector >& S2); +int nMieApplied(const unsigned int L, + const int pl, + std::vector& x, + std::vector >& m, + const unsigned int nTheta, + std::vector& Theta, + double* Qext, + double* Qsca, + double* Qabs, + double* Qbk, + double* Qpr, + double* g, + double* Albedo, + std::vector >& S1, + std::vector >& S2); +int nMieApplied(const unsigned int L, + std::vector& x, + std::vector >& m, + const unsigned int nTheta, + std::vector& Theta, + const int nmax, + double* Qext, + double* Qsca, + double* Qabs, + double* Qbk, + double* Qpr, + double* g, + double* Albedo, + std::vector >& S1, + std::vector >& S2); -namespace nmie { +template +class MultiLayerMieApplied : public MultiLayerMie { + // Will throw for any error! + public: + void RunMieCalculation(); + void RunFieldCalculationPolar(const int outer_arc_points, + const int radius_points, + const double from_Rho, + const double to_Rho, + const double from_Theta, + const double to_Theta, + const double from_Phi, + const double to_Phi, + const int isIgnoreAvailableNmax); + void RunFieldCalculationCartesian(const int first_side_points, + const int second_side_points, + const double relative_side_length, + const int plane_selected, + const double at_x, + const double at_y, + const double at_z, + const int isIgnoreAvailableNmax); + // emscripten::val GetFieldEabs(); + + void GetFailed(); + long iformat = 0; + bool output = true; + void prn(FloatType 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(FloatType wavelength) { + this->MultiLayerMie::MarkUncalculated(); + 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(FloatType layer_width, + std::complex layer_index); + void AddTargetLayerReIm(FloatType layer_width, + FloatType re_layer_index, + FloatType im_layer_index); + void AddCoatingLayer(FloatType layer_width, + std::complex layer_index); + // For all layers + void SetTargetWidth(std::vector width); + void SetTargetIndex(std::vector > index); + void SetTargetPEC(FloatType radius); + void SetCoatingWidth(std::vector width); + void SetCoatingIndex(std::vector > index); + void SetFieldPoints(std::vector > coords); + + // Set parameters in size parameter units + void SetWidthSP(const std::vector& width); + void SetIndexSP(const std::vector >& index); + void SetFieldPointsSP(const std::vector >& coords_sp); + + // Set common parameters + void SetModeNmaxAndType(int mode_n, int mode_type) { + this->MultiLayerMie::SetModeNmaxAndType(mode_n, mode_type); + }; + void SetAnglesForPattern(FloatType from_angle, + FloatType to_angle, + int samples); + std::vector GetAngles(); + + void ClearTarget(); + void ClearCoating(); + void ClearLayers(); + void ClearAllDesign(); // Layers + SP + index_ + + // Applied units requests + FloatType GetTotalRadius(); + FloatType GetTargetRadius(); + FloatType GetCoatingWidth(); + std::vector GetTargetLayersWidth(); + std::vector > GetTargetLayersIndex(); + std::vector GetCoatingLayersWidth(); + std::vector > GetCoatingLayersIndex(); + std::vector > GetFieldPoints(); + // ext, sca, abs, bk + std::vector > GetSpectra(FloatType from_WL, + FloatType to_WL, + int samples); + FloatType GetRCSext(); + FloatType GetRCSsca(); + FloatType GetRCSabs(); + FloatType GetRCSbk(); + std::vector GetPatternEk(); + std::vector GetPatternHk(); + std::vector GetPatternUnpolarized(); + + // Dimensionless + FloatType GetQsca() { return this->MultiLayerMie::GetQsca(); }; + FloatType GetQabs() { return this->MultiLayerMie::GetQabs(); }; + FloatType GetQext() { return this->MultiLayerMie::GetQext(); }; + // Size parameter units + std::vector GetLayerWidthSP(); + // Same as to get target and coating index + std::vector > GetLayerIndex(); + std::vector > GetFieldPointsSP(); + // Do we need normalize field to size parameter? + /* std::vector > > GetFieldESP(); */ + /* std::vector > > GetFieldHSP(); */ + std::vector > GetSpectraSP( + FloatType from_SP, + FloatType to_SP, + int samples); // WL,ext, sca, abs, bk + + std::vector GetPatternEkSP(); + std::vector GetPatternHkSP(); + std::vector GetPatternUnpolarizedSP(); + + void GetExpanCoeffs(std::vector > >& aln, + std::vector > >& bln, + std::vector > >& cln, + std::vector > >& dln); + + // 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 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); + + FloatType wavelength_ = 1.0; + FloatType total_radius_ = 0.0; + /// Width and index for each layer of the structure + std::vector target_width_, coating_width_; + std::vector > target_index_, coating_index_; + + std::vector > coords_sp_; - int nMieApplied(const unsigned int L, const int pl, std::vector &x, std::vector >& m, const unsigned int nTheta, std::vector& Theta, const int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector >& S1, std::vector >& S2); - int nMieApplied(const unsigned int L, std::vector& x, std::vector >& m, const unsigned int nTheta, std::vector& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector >& S1, std::vector >& S2); - int nMieApplied(const unsigned int L, const int pl, std::vector& x, std::vector >& m, const unsigned int nTheta, std::vector& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector >& S1, std::vector >& S2); - int nMieApplied(const unsigned int L, std::vector& x, std::vector >& m, const unsigned int nTheta, std::vector& Theta, const int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector >& S1, std::vector >& S2); - - - template - class MultiLayerMieApplied : public MultiLayerMie { - // Will throw for any error! - public: - void RunMieCalculation(); - void RunFieldCalculationPolar(const int outer_arc_points, - const int radius_points, - const double from_Rho, const double to_Rho, - const double from_Theta, const double to_Theta, - const double from_Phi, const double to_Phi, - const int isIgnoreAvailableNmax); - void RunFieldCalculationCartesian(const int first_side_points, - const int second_side_points, - const double relative_side_length, - const int plane_selected, - const double at_x, const double at_y, - const double at_z, - const int isIgnoreAvailableNmax); - emscripten::val GetFieldEabs(); - - void GetFailed(); - long iformat = 0; - bool output = true; - void prn(FloatType 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(FloatType wavelength) { - this->MultiLayerMie::MarkUncalculated(); - 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(FloatType layer_width, std::complex layer_index); - void AddTargetLayerReIm(FloatType layer_width, - FloatType re_layer_index, FloatType im_layer_index); - void AddCoatingLayer(FloatType layer_width, std::complex layer_index); - // For all layers - void SetTargetWidth(std::vector width); - void SetTargetIndex(std::vector< std::complex > index); - void SetTargetPEC(FloatType 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 SetModeNmaxAndType(int mode_n, int mode_type) { - this->MultiLayerMie::SetModeNmaxAndType(mode_n, mode_type);}; - void SetAnglesForPattern(FloatType from_angle, FloatType to_angle, int samples); - std::vector GetAngles(); - - void ClearTarget(); - void ClearCoating(); - void ClearLayers(); - void ClearAllDesign(); //Layers + SP + index_ - - // Applied units requests - FloatType GetTotalRadius(); - FloatType GetTargetRadius(); - FloatType GetCoatingWidth(); - std::vector GetTargetLayersWidth(); - std::vector< std::complex > GetTargetLayersIndex(); - std::vector GetCoatingLayersWidth(); - std::vector< std::complex > GetCoatingLayersIndex(); - std::vector< std::vector > GetFieldPoints(); - std::vector< std::vector > GetSpectra(FloatType from_WL, FloatType to_WL, int samples); // ext, sca, abs, bk - FloatType GetRCSext(); - FloatType GetRCSsca(); - FloatType GetRCSabs(); - FloatType GetRCSbk(); - std::vector GetPatternEk(); - std::vector GetPatternHk(); - std::vector GetPatternUnpolarized(); - - // Dimensionless - FloatType GetQsca(){ return this->MultiLayerMie::GetQsca();}; - FloatType GetQabs(){ return this->MultiLayerMie::GetQabs();}; - FloatType GetQext(){ return this->MultiLayerMie::GetQext();}; - // 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(FloatType from_SP, FloatType to_SP, int samples); // WL,ext, sca, abs, bk - - - std::vector GetPatternEkSP(); - std::vector GetPatternHkSP(); - std::vector GetPatternUnpolarizedSP(); - - void GetExpanCoeffs - (std::vector< std::vector > >& aln, - std::vector< std::vector > >& bln, - std::vector< std::vector > >& cln, - std::vector< std::vector > >& dln); - - - // 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 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); - - FloatType wavelength_ = 1.0; - FloatType 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_; - - std::vector< std::vector > coords_sp_; - - }; // end of class MultiLayerMie +}; // end of class MultiLayerMie } // end of namespace nmie #endif // SRC_NMIE_APPLIED_HPP diff --git a/src/nmie-js-wrapper.cc b/src/nmie-js-wrapper.cc index 792999f8..fc0b54d3 100644 --- a/src/nmie-js-wrapper.cc +++ b/src/nmie-js-wrapper.cc @@ -1,100 +1,111 @@ -//**********************************************************************************// -// Copyright (C) 2009-2019 Ovidio Pena // -// Copyright (C) 2013-2019 Konstantin Ladutenko // -// // -// This file is part of scattnlay // -// // -// This program is free software: you can redistribute it and/or modify // -// it under the terms of the GNU General Public License as published by // -// the Free Software Foundation, either version 3 of the License, or // -// (at your option) any later version. // -// // -// This program is distributed in the hope that it will be useful, // -// but WITHOUT ANY WARRANTY; without even the implied warranty of // -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // -// GNU General Public License for more details. // -// // -// The only additional remark is that we expect that all publications // -// describing work using this software, or all commercial products // -// using it, cite 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 . // -// // +//********************************************************************************** +// Copyright (C) 2009-2019 Ovidio Pena +// Copyright (C) 2013-2019 Konstantin Ladutenko +// +// This file is part of scattnlay +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// The only additional remark is that we expect that all publications +// describing work using this software, or all commercial products +// using it, cite 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 . +// // @brief Wrapper to JS -// // -//**********************************************************************************// -#include "nmie-applied.hpp" -#include "nmie-applied-impl.hpp" +// +//********************************************************************************** #include "nmie-precision.hpp" +#include "nmie-web-impl.hpp" using namespace emscripten; -nmie::MultiLayerMieApplied ml_mie; +nmie::MultiLayerMieWeb ml_mie; -EMSCRIPTEN_BINDINGS (c) { - class_>("nmie") - .constructor<>() - .function("SetWavelength", &nmie::MultiLayerMieApplied::SetWavelength) - .function("AddTargetLayerReIm",&nmie::MultiLayerMieApplied::AddTargetLayerReIm) - .function("SetModeNmaxAndType",&nmie::MultiLayerMieApplied::SetModeNmaxAndType) - .function("ClearTarget",&nmie::MultiLayerMieApplied::ClearTarget) - .function("RunMieCalculation",&nmie::MultiLayerMieApplied::RunMieCalculation) - .function("RunFieldCalculationPolar",&nmie::MultiLayerMieApplied::RunFieldCalculationPolar) - .function("RunFieldCalculationCartesian",&nmie::MultiLayerMieApplied::RunFieldCalculationCartesian) - .function("GetFieldEabs",&nmie::MultiLayerMieApplied::GetFieldEabs) - .function("GetQsca",&nmie::MultiLayerMieApplied::GetQsca) - .function("GetQext",&nmie::MultiLayerMieApplied::GetQext) - .function("GetQabs",&nmie::MultiLayerMieApplied::GetQabs) -// .function("bf",&nmie::MultiLayerMieApplied::bf) - ; +EMSCRIPTEN_BINDINGS(c) { + class_>("nmie") + .constructor<>() + .function("SetWavelength", &nmie::MultiLayerMieWeb::SetWavelength) + .function("AddTargetLayerReIm", + &nmie::MultiLayerMieWeb::AddTargetLayerReIm) + .function("SetModeNmaxAndType", + &nmie::MultiLayerMieWeb::SetModeNmaxAndType) + .function("ClearTarget", &nmie::MultiLayerMieWeb::ClearTarget) + .function("RunMieCalculation", + &nmie::MultiLayerMieWeb::RunMieCalculation) + .function("RunFieldCalculationPolar", + &nmie::MultiLayerMieWeb::RunFieldCalculationPolar) + .function("RunFieldCalculationCartesian", + &nmie::MultiLayerMieWeb::RunFieldCalculationCartesian) + .function("GetFieldEabs", &nmie::MultiLayerMieWeb::GetFieldEabs) + .function("GetQsca", &nmie::MultiLayerMieWeb::GetQsca) + .function("GetQext", &nmie::MultiLayerMieWeb::GetQext) + .function("GetQabs", &nmie::MultiLayerMieWeb::GetQabs) + // .function("bf",&nmie::MultiLayerMieWeb::bf) + ; } -//namespace nmie { - //**********************************************************************************// - // This function emulates a C call to calculate 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 // - //**********************************************************************************// -// int nMieApplied(const unsigned int L, const int pl, std::vector &x, std::vector > &m, const unsigned int nTheta, std::vector &Theta, const int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector >& S1, std::vector >& S2) { +// namespace nmie { +//********************************************************************************** +// This function emulates a C call to calculate 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 +//********************************************************************************** +// int nMieWeb(const unsigned int L, const int pl, std::vector &x, +// std::vector > &m, const unsigned int nTheta, +// std::vector &Theta, const int nmax, double *Qext, double *Qsca, +// double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, +// std::vector >& S1, std::vector >& +// S2) { // // if (x.size() != L || m.size() != L) -// throw std::invalid_argument("Declared number of layers do not fit x and m!"); +// throw std::invalid_argument("Declared number of layers do not fit x +// and m!"); // if (Theta.size() != nTheta) -// throw std::invalid_argument("Declared number of sample for Theta is not correct!"); +// throw std::invalid_argument("Declared number of sample for Theta is +// not correct!"); // try { -// MultiLayerMieApplied ml_mie; +// MultiLayerMieWeb ml_mie; // ml_mie.SetLayersSize(ConvertVector(x)); // ml_mie.SetLayersIndex(ConvertComplexVector(m)); // ml_mie.SetAngles(ConvertVector(Theta)); @@ -122,17 +133,35 @@ EMSCRIPTEN_BINDINGS (c) { // } // return 0; // } -// int nMieApplied(const unsigned int L, std::vector& x, std::vector >& m, const unsigned int nTheta, std::vector& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector >& S1, std::vector >& S2) { -// return nmie::nMieApplied(L, -1, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2); +// int nMieWeb(const unsigned int L, std::vector& x, +// std::vector >& m, const unsigned int nTheta, +// std::vector& Theta, double *Qext, double *Qsca, double *Qabs, double +// *Qbk, double *Qpr, double *g, double *Albedo, +// std::vector >& S1, std::vector >& +// S2) { +// return nmie::nMieWeb(L, -1, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, +// Qbk, Qpr, g, Albedo, S1, S2); // } -// int nMieApplied(const unsigned int L, const int pl, std::vector& x, std::vector >& m, const unsigned int nTheta, std::vector& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector >& S1, std::vector >& S2) { -// return nmie::nMieApplied(L, pl, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2); +// int nMieWeb(const unsigned int L, const int pl, std::vector& x, +// std::vector >& m, const unsigned int nTheta, +// std::vector& Theta, double *Qext, double *Qsca, double *Qabs, double +// *Qbk, double *Qpr, double *g, double *Albedo, +// std::vector >& S1, std::vector >& +// S2) { +// return nmie::nMieWeb(L, pl, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, +// Qbk, Qpr, g, Albedo, S1, S2); // } -// int nMieApplied(const unsigned int L, std::vector& x, std::vector >& m, const unsigned int nTheta, std::vector& Theta, const int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector >& S1, std::vector >& S2) { -// return nmie::nMieApplied(L, -1, x, m, nTheta, Theta, nmax, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2); +// int nMieWeb(const unsigned int L, std::vector& x, +// std::vector >& m, const unsigned int nTheta, +// std::vector& Theta, const int nmax, double *Qext, double *Qsca, +// double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, +// std::vector >& S1, std::vector >& +// S2) { +// return nmie::nMieWeb(L, -1, x, m, nTheta, Theta, nmax, Qext, Qsca, Qabs, +// Qbk, Qpr, g, Albedo, S1, S2); // } - // ********************************************************************** // - // ********************************************************************** // - // ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // //} // end of namespace nmie diff --git a/src/nmie-precision.hpp b/src/nmie-precision.hpp index f59d2019..1851b856 100644 --- a/src/nmie-precision.hpp +++ b/src/nmie-precision.hpp @@ -1,117 +1,130 @@ #ifndef SRC_NMIE_PRECISION_H_ #define SRC_NMIE_PRECISION_H_ //**********************************************************************************// -// Copyright (C) 2009-2018 Ovidio Pena // -// Copyright (C) 2013-2018 Konstantin Ladutenko // +// Copyright (C) 2009-2018 Ovidio Pena // Copyright +// (C) 2013-2018 Konstantin Ladutenko // // // -// This file is part of scattnlay // +// 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 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. // +// 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. // +// 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 . // +// You should have received a copy of the GNU General Public License // along +// with this program. If not, see . // //**********************************************************************************// //**********************************************************************************// -// 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. // +// 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. // +// 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. // // // -// Hereinafter all equations numbers refer to [2] // +// Hereinafter all equations numbers refer to [2] // //**********************************************************************************// +#include +#include #ifdef MULTI_PRECISION -#include #include +#include #endif // MULTI_PRECISION namespace nmie { - #ifdef MULTI_PRECISION - namespace nmm = boost::multiprecision; - typedef nmm::number > FloatType; - #else - namespace nmm = std; - typedef double FloatType; - //typedef float FloatType; - #endif // MULTI_PRECISION - - template T sin_t(T v) { - if (std::is_same::value) return static_cast(std::sin(static_cast(v))); - return static_cast(nmm::sin(static_cast(v))); - } - template T cos_t(T v) { - if (std::is_same::value) return static_cast(std::cos(static_cast(v))); - return static_cast(nmm::cos(static_cast(v))); - } - template T sqrt_t(T v) { - if (std::is_same::value) return static_cast(std::sqrt(static_cast(v))); - return static_cast(nmm::sqrt(static_cast(v))); - } - +#ifdef MULTI_PRECISION +namespace nmm = boost::multiprecision; +typedef nmm::number > FloatType; +#else +namespace nmm = std; +typedef double FloatType; +// typedef float FloatType; +#endif // MULTI_PRECISION +template +T sin_t(T v) { + if (std::is_same::value) + return static_cast(std::sin(static_cast(v))); + return static_cast(nmm::sin(static_cast(v))); +} +template +T cos_t(T v) { + if (std::is_same::value) + return static_cast(std::cos(static_cast(v))); + return static_cast(nmm::cos(static_cast(v))); +} +template +T sqrt_t(T v) { + if (std::is_same::value) + return static_cast(std::sqrt(static_cast(v))); + return static_cast(nmm::sqrt(static_cast(v))); +} template - std::vector ConvertVector(const std::vector x) { - std::vector new_x; - for (auto element : x) { - new_x.push_back(static_cast(element)); - } - return new_x; +std::vector ConvertVector(const std::vector x) { + std::vector new_x; + for (auto element : x) { + new_x.push_back(static_cast(element)); } + return new_x; +} - template - std::complex ConvertComplex(std::complex z) { - return std::complex(static_cast(z.real()), - static_cast(z.imag())); - } +template +std::complex ConvertComplex(std::complex z) { + return std::complex(static_cast(z.real()), + static_cast(z.imag())); +} - template - std::vector > ConvertComplexVector(std::vector > x) { - std::vector > new_x; - for (auto element : x) { - new_x.push_back(std::complex(static_cast(element.real()), - static_cast(element.imag()) ) ); - } - return new_x; +template +std::vector > ConvertComplexVector( + std::vector > x) { + std::vector > new_x; + for (auto element : x) { + new_x.push_back( + std::complex(static_cast(element.real()), + static_cast(element.imag()))); } + return new_x; +} - template - std::vector > > ConvertComplexVectorVector(std::vector > > x) { - std::vector > > new_x; - std::vector > new_y; - for (auto y : x) { - new_y.clear(); - for (auto element : y) { - new_y.push_back(std::complex(static_cast(element.real()), - static_cast(element.imag()) ) ); - } - new_x.push_back(new_y); +template +std::vector > > +ConvertComplexVectorVector( + std::vector > > x) { + std::vector > > new_x; + std::vector > new_y; + for (auto y : x) { + new_y.clear(); + for (auto element : y) { + new_y.push_back( + std::complex(static_cast(element.real()), + static_cast(element.imag()))); } - return new_x; + new_x.push_back(new_y); } + return new_x; +} } // end of namespace nmie #endif // SRC_NMIE_PRECISION_H_ diff --git a/src/nmie-web-impl.hpp b/src/nmie-web-impl.hpp new file mode 100644 index 00000000..8edb1aeb --- /dev/null +++ b/src/nmie-web-impl.hpp @@ -0,0 +1,79 @@ +#ifndef SRC_NMIE_WEB_IMPL_HPP_ +#define SRC_NMIE_WEB_IMPL_HPP_ +//********************************************************************************** +// Copyright (C) 2009-2024 Ovidio Pena +// Copyright (C) 2013-2024 Konstantin Ladutenko +// +// This file is part of scattnlay +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// The only additional remark is that we expect that all publications +// describing work using this software, or all commercial products +// using it, cite 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 . +// +// @brief Wrapper class around nMie function for ease of use +// +//********************************************************************************** +#include +#include +#include +#include +#include +#include + +#include +#include +#include "nmie-web.hpp" + +namespace nmie { +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // + +// from https://toughengineer.github.io/demo/dsp/fft-perf/ +template +emscripten::val toJSFloat64Array(const std::vector& v) { + emscripten::val view{emscripten::typed_memory_view( + v.size(), v.data())}; // make a view of local object + + auto result = emscripten::val::global("Float64Array") + .new_(v.size()); // make a JS typed array + result.call( + "set", + view); // set typed array values "on the JS side" using the memory view + + return result; +} + +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // +template +emscripten::val MultiLayerMieWeb::GetFieldEabs() { + auto Eabs = this->MultiLayerMie::GetFieldEabs(); + return toJSFloat64Array(Eabs); +} +// ********************************************************************** // +// ********************************************************************** // +// ********************************************************************** // +} // end of namespace nmie +#endif // SRC_NMIE_WEB_IMPL_HPP_ diff --git a/src/nmie-web.hpp b/src/nmie-web.hpp new file mode 100644 index 00000000..8255a832 --- /dev/null +++ b/src/nmie-web.hpp @@ -0,0 +1,56 @@ +#ifndef SRC_NMIE_WEB_HPP_ +#define SRC_NMIE_WEB_HPP_ +//********************************************************************************** +// Copyright (C) 2009-2024 Ovidio Pena +// Copyright (C) 2013-2024 Konstantin Ladutenko +// +// This file is part of scattnlay +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// The only additional remark is that we expect that all publications +// describing work using this software, or all commercial products +// using it, cite 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 . +//********************************************************************************** + +#include +#include +#include +#include +#include + +#include "nmie-applied-impl.hpp" + +#include +#include + +namespace nmie { + +template +class MultiLayerMieWeb : public MultiLayerMieApplied { + // Will throw for any error! + public: + emscripten::val GetFieldEabs(); + +}; // end of class MultiLayerMieWeb + +} // end of namespace nmie +#endif // SRC_NMIE_WEB_HPP diff --git a/tests/c++/go-speed-test.sh b/tests/c++/go-speed-test.sh index b256cd21..b7d6e775 100755 --- a/tests/c++/go-speed-test.sh +++ b/tests/c++/go-speed-test.sh @@ -15,7 +15,11 @@ rm -f $PROGRAM file=speed-test-applied.cc # g++ -Ofast -std=c++11 $file ../../src/nmie.cc ../../src/nmie-applied.cc -DMULTI_PRECISION=200 -lm -lrt -o $PROGRAM /usr/lib/libtcmalloc_minimal.so.4 -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2 -g++ -Ofast -std=c++11 $file ../../src/nmie.cc ../../src/nmie-applied.cc -lm -lrt -o $PROGRAM /usr/lib/libtcmalloc_minimal.so.4 -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2 + +# with libtcmalloc +# g++ -Ofast -std=c++11 $file ../../src/nmie.cc ../../src/nmie-applied.cc -lm -lrt -o $PROGRAM /usr/lib/libtcmalloc_minimal.so.4 -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2 + +g++ -Ofast -std=c++11 $file ../../src/nmie.cc ../../src/nmie-applied.cc -lm -o $PROGRAM -march=native -mtune=native echo Should be: echo test01, +1.41154e+00, +4.17695e-01, +9.93844e-01, +1.59427e-01, +1.25809e+00, +3.67376e-01, +2.95915e-01 diff --git a/tests/c++/speed-test.cc b/tests/c++/speed-test.cc index f91d7376..d4b516ac 100644 --- a/tests/c++/speed-test.cc +++ b/tests/c++/speed-test.cc @@ -1,30 +1,34 @@ -//**********************************************************************************// -// Copyright (C) 2009-2016 Ovidio Pena // -// Copyright (C) 2013-2016 Konstantin Ladutenko // -// // -// This file is part of scattnlay // -// // -// This program is free software: you can redistribute it and/or modify // -// it under the terms of the GNU General Public License as published by // -// the Free Software Foundation, either version 3 of the License, or // -// (at your option) any later version. // -// // -// This program is distributed in the hope that it will be useful, // -// but WITHOUT ANY WARRANTY; without even the implied warranty of // -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // -// GNU General Public License for more details. // -// // -// The only additional remark is that we expect that all publications // -// describing work using this software, or all commercial products // -// using it, cite the following reference: // -// [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by // -// a multilayered sphere," Computer Physics Communications, // -// vol. 180, Nov. 2009, pp. 2348-2354. // -// // -// You should have received a copy of the GNU General Public License // -// along with this program. If not, see . // -//**********************************************************************************// +//********************************************************************************** +// Copyright (C) 2009-2016 Ovidio Pena +// Copyright (C) 2013-2016 Konstantin Ladutenko +// +// This file is part of scattnlay +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// The only additional remark is that we expect that all publications +// describing work using this software, or all commercial products +// using it, cite the following reference: +// [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by +// a multilayered sphere," Computer Physics Communications, +// vol. 180, Nov. 2009, pp. 2348-2354. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . +//********************************************************************************** +#include +#include +#include +#include #include #include #include @@ -32,47 +36,54 @@ #include #include #include -#include -#include -#include -#include -//sudo aptitude install libgoogle-perftools-dev -//#include -#include "../../src/nmie.hpp" -//#include "../../src/nmie-precision.hpp" -//#include "../../src/nmie-impl.hpp" +// sudo aptitude install libgoogle-perftools-dev +// #include +#include "../../src/nmie-applied-impl.hpp" timespec diff(timespec start, timespec end); -const double PI=3.14159265358979323846; -template inline T pow2(const T value) {return value*value;} - +const double PI = 3.14159265358979323846; +template +inline T pow2(const T value) { + return value * value; +} -//***********************************************************************************// -// This is the main function of 'scattnlay', here we read the parameters as // -// arguments passed to the program which should be executed with the following // -// syntaxis: // -// ./scattnlay -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] [-t ti tf nt] [-c comment] // -// // -// When all the parameters were correctly passed we setup the integer L (the // -// number of layers) and the arrays x and m, containing the size parameters and // -// refractive indexes of the layers, respectively and call the function nMie. // -// If the calculation is successful the results are printed with the following // -// format: // -// // -// * If no comment was passed: // -// 'Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo' // -// // -// * If a comment was passed: // -// 'comment, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo' // -//***********************************************************************************// -int main(int argc, char *argv[]) { +//*********************************************************************************** +// This is the main function of 'scattnlay', here we read the parameters as +// arguments passed to the program which should be executed with the following +// syntaxis: +// ./scattnlay -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] [-t ti tf nt] [-c +// comment] +// +// When all the parameters were correctly passed we setup the integer L (the +// number of layers) and the arrays x and m, containing the size parameters and +// refractive indexes of the layers, respectively and call the function nMie. +// If the calculation is successful the results are printed with the following +// format: +// +// * If no comment was passed: +// 'Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo' +// +// * If a comment was passed: +// 'comment, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo' +//*********************************************************************************** +int main(int argc, char* argv[]) { try { std::vector args; args.assign(argv, argv + argc); - std::string error_msg(std::string("Insufficient parameters.\nUsage: ") + args[0] - + " -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] " - + "[-t ti tf nt] [-c comment]\n"); - enum mode_states {read_L, read_x, read_mr, read_mi, read_ti, read_tf, read_nt, read_comment}; + std::string error_msg(std::string("Insufficient parameters.\nUsage: ") + + args[0] + + " -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] " + + "[-t ti tf nt] [-c comment]\n"); + enum mode_states { + read_L, + read_x, + read_mr, + read_mi, + read_ti, + read_tf, + read_nt, + read_comment + }; // for (auto arg : args) std::cout<< arg <( tmp_mr,std::stod(arg) )); - mode = read_x; - continue; + m.push_back(std::complex(tmp_mr, std::stod(arg))); + mode = read_x; + continue; } if (mode == read_ti) { - ti = std::stod(arg); - mode = read_tf; - continue; + ti = std::stod(arg); + mode = read_tf; + continue; } if (mode == read_tf) { - tf = std::stod(arg); - mode = read_nt; - continue; + tf = std::stod(arg); + mode = read_nt; + continue; } if (mode == read_nt) { - nt = std::stoi(arg); + nt = std::stoi(arg); Theta.resize(nt); S1.resize(nt); S2.resize(nt); S1w.resize(nt); S2w.resize(nt); - continue; + continue; } - if (mode == read_comment) { - comment = arg; + if (mode == read_comment) { + comment = arg; has_comment = 1; - continue; + continue; } } - if ( (x.size() != m.size()) || (L != x.size()) ) - throw std::invalid_argument(std::string("Broken structure!\n") - +error_msg); - if ( (0 == m.size()) || ( 0 == x.size()) ) - throw std::invalid_argument(std::string("Empty structure!\n") - +error_msg); + if ((x.size() != m.size()) || (L != x.size())) + throw std::invalid_argument(std::string("Broken structure!\n") + + error_msg); + if ((0 == m.size()) || (0 == x.size())) + throw std::invalid_argument(std::string("Empty structure!\n") + + error_msg); if (nt < 0) { printf("Error reading Theta.\n"); return -1; } else if (nt == 1) { - Theta[0] = ti*PI/180.0; + Theta[0] = ti * PI / 180.0; } else { for (i = 0; i < nt; i++) { - Theta[i] = (ti + (double)i*(tf - ti)/(nt - 1))*PI/180.0; + Theta[i] = (ti + (double)i * (tf - ti) / (nt - 1)) * PI / 180.0; } } timespec time1, time2; @@ -184,59 +197,60 @@ int main(int argc, char *argv[]) { long ctime_nsec, best_c; long cpptime_sec, ctime_sec; long repeats = 150; - //HeapProfilerStart("heapprof"); + // HeapProfilerStart("heapprof"); printf("Start\n"); do { clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time1); - for (int i = 0; i 0) { - // printf(" Theta, S1.r, S1.i, S2.r, S2.i\n"); + // printf(" Theta, S1.r, S1.i, S2.r, S2.i\n"); // for (i = 0; i < nt; i++) { - // printf("%6.2f, %+.5e, %+.5e, %+.5e, %+.5e \n", Theta[i]*180.0/PI, S1w[i].real(), S1w[i].imag(), S2w[i].real(), S2w[i].imag()); + // printf("%6.2f, %+.5e, %+.5e, %+.5e, %+.5e \n", Theta[i]*180.0/PI, + // S1w[i].real(), S1w[i].imag(), S2w[i].real(), S2w[i].imag()); // } // } - } catch( const std::invalid_argument &ia ) { + } catch (const std::invalid_argument& ia) { // Will catch if multi_layer_mie fails or other errors. std::cerr << "Invalid argument: " << ia.what() << std::endl; return -1; } - return 0; + return 0; } - - -timespec diff(timespec start, timespec end) -{ - timespec temp; - if ((end.tv_nsec-start.tv_nsec)<0) { - temp.tv_sec = end.tv_sec-start.tv_sec-1; - temp.tv_nsec = 1000000000+end.tv_nsec-start.tv_nsec; - } else { - temp.tv_sec = end.tv_sec-start.tv_sec; - temp.tv_nsec = end.tv_nsec-start.tv_nsec; - } - return temp; +timespec diff(timespec start, timespec end) { + timespec temp; + if ((end.tv_nsec - start.tv_nsec) < 0) { + temp.tv_sec = end.tv_sec - start.tv_sec - 1; + temp.tv_nsec = 1000000000 + end.tv_nsec - start.tv_nsec; + } else { + temp.tv_sec = end.tv_sec - start.tv_sec; + temp.tv_nsec = end.tv_nsec - start.tv_nsec; + } + return temp; } diff --git a/tests/test_py.py b/tests/test_py.py index 83cf1865..941aa633 100644 --- a/tests/test_py.py +++ b/tests/test_py.py @@ -56,14 +56,15 @@ def test_bulk_multilayer(self): continue print("Using solver: ", solver) for case in test_cases: - print("test case:", case) solver.SetLayersSize(case[0]) solver.SetLayersIndex(case[1]) solver.RunMieCalculation() Qext = solver.GetQext() Qsca = solver.GetQsca() - print("ext tol:", (case[2] - Qext) / Qext) - print("sca tol:", (case[3] - Qsca) / Qsca) + if case == test_cases[0]: + print("test case:", case) + print("ext tol:", (case[2] - Qext) / Qext) + print("sca tol:", (case[3] - Qsca) / Qsca) self.assertTrue((case[2] - Qext) / Qext < tol) self.assertTrue((case[3] - Qsca) / Qsca < tol)