mesh on sphere

examples/go-cc-examples.sh

@@ -2,24 +2,32 @@
 path=$PWD
 PROGRAM='scattnlay-example.bin'
 
-file=example-minimal.cc
+file=example-eval-force.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 -lrt -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
-# Production for multiprecision
-#g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc read-spectra.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
 
-# Simplified for multiprecision
-#g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc ./read-spectra.cc -DMULTI_PRECISION=200 -lm -lrt -o $PROGRAM -march=native -mtune=native -msse4.2
+# file=example-get-Mie.cc
+# echo Compile $file with gcc
+# rm -f $PROGRAM
+# # Production for multiprecision
+# #g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc read-spectra.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
 
-# Simplified for double precision
-g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc ./read-spectra.cc -lm -lrt -o $PROGRAM -march=native -mtune=native -msse4.2
-echo Compilation done. Running...
-./$PROGRAM
+# # Simplified for multiprecision
+# #g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc ./read-spectra.cc -DMULTI_PRECISION=200 -lm -lrt -o $PROGRAM -march=native -mtune=native -msse4.2
+
+# # Simplified for double precision
+# g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc ./read-spectra.cc -lm -lrt -o $PROGRAM -march=native -mtune=native -msse4.2
+# echo Compilation done. Running...
+# ./$PROGRAM

src/nmie-applied-impl.hpp

@@ -0,0 +1,416 @@
+#ifndef SRC_NMIE_APPLIED_IMPL_HPP_
+#define SRC_NMIE_APPLIED_IMPL_HPP_
+///
+/// @file   nmie-applied-impl.hpp
+/// @author Ladutenko Konstantin <kostyfisik at gmail (.) com>
+/// @date   Tue Sep  3 00:38:27 2013
+/// @copyright 2013-2016 Ladutenko Konstantin
+///
+/// nmie is free software: you can redistribute it and/or modify
+/// it under the terms of the GNU General Public License as published by
+/// the Free Software Foundation, either version 3 of the License, or
+/// (at your option) any later version.
+///
+/// nmie-wrapper is distributed in the hope that it will be useful,
+/// but WITHOUT ANY WARRANTY; without even the implied warranty of
+/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+/// GNU General Public License for more details.
+///
+/// You should have received a copy of the GNU General Public License
+/// along with nmie-wrapper.  If not, see <http://www.gnu.org/licenses/>.
+///
+/// nmie uses nmie.c from scattnlay by Ovidio Pena
+/// <ovidio@bytesfall.com> . He has an additional condition to 
+/// his library:
+//    The only additional condition is that we expect that all publications         //
+//    describing  work using this software , or all commercial products             //
+//    using it, cite the following reference:                                       //
+//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by           //
+//        a multilayered sphere," Computer Physics Communications,                  //
+//        vol. 180, Nov. 2009, pp. 2348-2354.                                       //
+///
+/// @brief  Wrapper class around nMie function for ease of use
+///
+#include "nmie-applied.hpp"
+#include "nmie-precision.hpp"
+#include <array>
+#include <algorithm>
+#include <cstdio>
+#include <cstdlib>
+#include <stdexcept>
+#include <vector>
+
+namespace nmie {  
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  template <typename FloatType>
+  void MultiLayerMieApplied<FloatType>::GetFailed() {
+    FloatType faild_x = 9.42477796076938;
+    //FloatType faild_x = 9.42477796076937;
+    std::complex<FloatType> z(faild_x, 0.0);
+    std::vector<int> nmax_local_array = {20, 100, 500, 2500};
+    for (auto nmax_local : nmax_local_array) {
+      std::vector<std::complex<FloatType> > D1_failed(nmax_local + 1);
+      // Downward recurrence for D1 - equations (16a) and (16b)
+      D1_failed[nmax_local] = std::complex<FloatType>(0.0, 0.0);
+      const std::complex<FloatType> zinv = std::complex<FloatType>(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);
+  
+    
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  template <typename FloatType>
+  void MultiLayerMieApplied<FloatType>::AddTargetLayer(FloatType width, std::complex<FloatType> 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<FloatType>::AddTargetLayer(...)  
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  template <typename FloatType>
+  void MultiLayerMieApplied<FloatType>::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<FloatType>(0.0, 0.0));
+    // ... and mark it as PEC
+    this->SetPECLayer(0);
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  template <typename FloatType>
+  void MultiLayerMieApplied<FloatType>::SetCoatingIndex(std::vector<std::complex<FloatType> > index) {
+    this->MarkUncalculated();
+    coating_index_.clear();
+    for (auto value : index) coating_index_.push_back(value);
+  }  // end of void MultiLayerMieApplied<FloatType>::SetCoatingIndex(std::vector<complex> index);  
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  template <typename FloatType>
+  void MultiLayerMieApplied<FloatType>::SetCoatingWidth(std::vector<FloatType> 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<FloatType>::SetCoatingWidth(...);
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  template <typename FloatType>
+  void MultiLayerMieApplied<FloatType>::SetWidthSP(const std::vector<FloatType>& 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<FloatType>::SetWidthSP(...);
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  template <typename FloatType>
+  void MultiLayerMieApplied<FloatType>::SetIndexSP(const std::vector< std::complex<FloatType> >& index) {
+    this->MarkUncalculated();
+    //refractive_index_.clear();
+    this->refractive_index_ = index;
+    // for (auto value : index) refractive_index_.push_back(value);
+  }  // end of void MultiLayerMieApplied<FloatType>::SetIndexSP(...);  
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  template <typename FloatType>
+  void MultiLayerMieApplied<FloatType>::SetFieldPointsSP(const std::vector< std::vector<FloatType> >& 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!");
+    coords_sp_ = 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<FloatType>::SetFieldPointsSP(...)
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  template <typename FloatType>
+  void MultiLayerMieApplied<FloatType>::GenerateSizeParameter() {
+    this->MarkUncalculated();
+    this->size_param_.clear();
+    FloatType radius = 0.0;
+    for (auto width : target_width_) {
+      radius += width;
+      this->size_param_.push_back(2*this->PI_*radius/wavelength_);
+    }
+    for (auto width : coating_width_) {
+      radius += width;
+      this->size_param_.push_back(2*this->PI_*radius/wavelength_);
+    }
+    this->total_radius_ = radius;
+  }  // end of void MultiLayerMieApplied<FloatType>::GenerateSizeParameter();
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  template <typename FloatType>
+  void MultiLayerMieApplied<FloatType>::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<FloatType>::GenerateIndex();
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  template <typename FloatType>
+  FloatType MultiLayerMieApplied<FloatType>::GetTotalRadius() {
+    if (!this->isMieCalculated())  GenerateSizeParameter();
+    return this->total_radius_;      
+  }  // end of FloatType MultiLayerMieApplied<FloatType>::GetTotalRadius();
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  template <typename FloatType>
+  std::vector< std::vector<FloatType> >
+  MultiLayerMieApplied<FloatType>::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<FloatType> > spectra;
+    FloatType step_WL = (to_WL - from_WL)/static_cast<FloatType>(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<FloatType>({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 <typename FloatType>
+  void MultiLayerMieApplied<FloatType>::ClearTarget() {
+    this->MarkUncalculated();
+    this->target_width_.clear();
+    this->target_index_.clear();
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  template <typename FloatType>
+  void MultiLayerMieApplied<FloatType>::ClearCoating() {
+    this->MarkUncalculated();
+    this->coating_width_.clear();
+    this->coating_index_.clear();
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  template <typename FloatType>
+  void MultiLayerMieApplied<FloatType>::ClearLayers() {
+    this->MarkUncalculated();
+    this->ClearTarget();
+    this->ClearCoating();
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  template <typename FloatType>
+  void MultiLayerMieApplied<FloatType>::ClearAllDesign() {
+    this->MarkUncalculated();
+    this->ClearLayers();
+    this->size_param_.clear();
+    this->refractive_index_.clear();
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  //                         Computational core
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  //**********************************************************************************//
+  // Function CONFRA ported from MIEV0.f (Wiscombe,1979)
+  // Ref. to NCAR Technical Notes, Wiscombe, 1979
+  /*
+c         Compute Bessel function ratio A-sub-N from its
+c         continued fraction using Lentz method
+
+c         ZINV = Reciprocal of argument of A
+
+
+c    I N T E R N A L    V A R I A B L E S
+c    ------------------------------------
+
+c    CAK      Term in continued fraction expansion of A (Eq. R25)
+c     a_k
+
+c    CAPT     Factor used in Lentz iteration for A (Eq. R27)
+c     T_k
+
+c    CNUMER   Numerator   in capT  (Eq. R28A)
+c     N_k
+c    CDENOM   Denominator in capT  (Eq. R28B)
+c     D_k
+
+c    CDTD     Product of two successive denominators of capT factors
+c                 (Eq. R34C)
+c     xi_1
+
+c    CNTN     Product of two successive numerators of capT factors
+c                 (Eq. R34B)
+c     xi_2
+
+c    EPS1     Ill-conditioning criterion
+c    EPS2     Convergence criterion
+
+c    KK       Subscript k of cAk  (Eq. R25B)
+c     k
+
+c    KOUNT    Iteration counter (used to prevent infinite looping)
+
+c    MAXIT    Max. allowed no. of iterations
+
+c    MM + 1  and - 1, alternately
+*/
+  template <typename FloatType>
+  std::complex<FloatType> MultiLayerMieApplied<FloatType>::calcD1confra(const int N, const std::complex<FloatType> z) {
+  // NTMR -> nmax_ - 1  \\TODO nmax_ ?
+    //int N = nmax_ - 1;
+    int KK, KOUNT, MAXIT = 10000, MM;
+    //    FloatType EPS1=1.0e-2;
+    FloatType EPS2=1.0e-8;
+    std::complex<FloatType> CAK, CAPT, CDENOM, CDTD, CNTN, CNUMER;
+    std::complex<FloatType> one = std::complex<FloatType>(1.0,0.0);
+    std::complex<FloatType> ZINV = one/z;
+// c                                 ** Eq. R25a
+    std::complex<FloatType> CONFRA = static_cast<std::complex<FloatType> >(N + 1)*ZINV;   //debug ZINV
+    MM = - 1; 
+    KK = 2*N +3; //debug 3
+// c                                 ** Eq. R25b, k=2
+    CAK    = static_cast<std::complex<FloatType> >(MM*KK)*ZINV; //debug -3 ZINV
+    CDENOM = CAK;
+    CNUMER = CDENOM + one/CONFRA; //-3zinv+z
+    KOUNT  = 1;
+    //10 CONTINUE
+    do {      ++KOUNT;
+      if (KOUNT > MAXIT) {
+        printf("re(%g):im(%g)\t\n", CONFRA.real(), CONFRA.imag());
+        throw std::invalid_argument("ConFra--Iteration failed to converge!\n");
+      }
+      MM *= - 1;      KK += 2;  //debug  mm=1 kk=5
+      CAK = static_cast<std::complex<FloatType> >(MM*KK)*ZINV; //    ** Eq. R25b //debug 5zinv
+     //  //c ** Eq. R32    Ill-conditioned case -- stride two terms instead of one
+     //  if (std::abs(CNUMER/CAK) >= EPS1 ||  std::abs(CDENOM/CAK) >= EPS1) {
+     //         //c                       ** Eq. R34
+     //         CNTN   = CAK*CNUMER + 1.0;
+     //         CDTD   = CAK*CDENOM + 1.0;
+     //         CONFRA = (CNTN/CDTD)*CONFRA; // ** Eq. R33
+     //         MM  *= - 1;        KK  += 2;
+     //         CAK = static_cast<std::complex<FloatType> >(MM*KK)*ZINV; // ** Eq. R25b
+     //         //c                        ** Eq. R35
+     //         CNUMER = CAK + CNUMER/CNTN;
+     //         CDENOM = CAK + CDENOM/CDTD;
+     //         ++KOUNT;
+     //         //GO TO  10
+     //         continue;
+     // } else { //c                           *** Well-conditioned case
+      {
+        CAPT   = CNUMER/CDENOM; // ** Eq. R27 //debug (-3zinv + z)/(-3zinv)
+        // printf("re(%g):im(%g)**\t", CAPT.real(), CAPT.imag());
+       CONFRA = CAPT*CONFRA; // ** Eq. R26
+       //if (N == 0) {output=true;printf(" re:");prn(CONFRA.real());printf(" im:"); prn(CONFRA.imag());output=false;};
+       //c                                  ** Check for convergence; Eq. R31
+       if (std::abs(CAPT.real() - 1.0) >= EPS2 ||  std::abs(CAPT.imag()) >= EPS2) {
+//c                                        ** Eq. R30
+         CNUMER = CAK + one/CNUMER;
+         CDENOM = CAK + one/CDENOM;
+         continue;
+         //GO TO  10
+       }  // end of if < eps2
+      }
+      break;
+    } while(1);    
+    //if (N == 0)  printf(" return confra for z=(%g,%g)\n", ZINV.real(), ZINV.imag());
+    return CONFRA;
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  template <typename FloatType>
+  void MultiLayerMieApplied<FloatType>::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 <typename FloatType>
+  void MultiLayerMieApplied<FloatType>::RunMieCalculation() {
+    ConvertToSP();
+    this->MultiLayerMie<FloatType>::RunMieCalculation(); 
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  template <typename FloatType>
+  void MultiLayerMieApplied<FloatType>::GetExpanCoeffs( std::vector< std::vector<std::complex<FloatType> > >& aln, std::vector< std::vector<std::complex<FloatType> > >& bln, std::vector< std::vector<std::complex<FloatType> > >& cln, std::vector< std::vector<std::complex<FloatType> > >& 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<FloatType>::GetExpanCoeffs( ...)
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+
+}  // end of namespace nmie
+#endif  // SRC_NMIE_APPLIED_IMPL_HPP_

src/shell-generator.cc

@@ -0,0 +1,383 @@
+//**********************************************************************************//
+//    Copyright (C) 2009-2016  Ovidio Pena <ovidio@bytesfall.com>                   //
+//    Copyright (C) 2013-2016  Konstantin Ladutenko <kostyfisik@gmail.com>          //
+//                                                                                  //
+//    This file is part of scattnlay                                                //
+//                                                                                  //
+//    This program is free software: you can redistribute it and/or modify          //
+//    it under the terms of the GNU General Public License as published by          //
+//    the Free Software Foundation, either version 3 of the License, or             //
+//    (at your option) any later version.                                           //
+//                                                                                  //
+//    This program is distributed in the hope that it will be useful,               //
+//    but WITHOUT ANY WARRANTY; without even the implied warranty of                //
+//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                 //
+//    GNU General Public License for more details.                                  //
+//                                                                                  //
+//    The only additional remark is that we expect that all publications            //
+//    describing work using this software, or all commercial products               //
+//    using it, cite 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 <http://www.gnu.org/licenses/>.         //
+//**********************************************************************************//
+//   @brief  Generates points for integration on sphere surface
+
+#include <algorithm>
+#include <cmath>
+#include <complex>
+#include <cstdio>
+#include <fstream>
+#include <iostream>
+#include <set>
+#include <sstream>
+#include <stdexcept>
+#include <string>
+#include <vector>
+#include "shell-generator.hpp"
+namespace shell_generator {
+  template<class T> inline T pow2(const T value) {return value*value;}
+
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  double ShellGenerator::IntegrateGauss(double charge, double shift) {
+    double integral = 0.0;
+    auto E = E_[0];
+    auto H = H_[0];
+
+    // return field at point p from the charge, located at (shift, 0,0)
+    auto field = [](double charge, double shift, std::vector<double> p){
+      double r = std::sqrt(pow2(p[0]-shift) + pow2(p[1]) + pow2(p[2]) );
+      double ampl = charge/pow2(r);
+      std::vector<double> field = {ampl*(p[0]-shift)/r, ampl*(p[1])/r, ampl*(p[2])/r};
+      return field;
+    };
+    
+    for(const auto face : faces_){
+      std::vector<double> mean_vector = {0.0, 0.0, 0.0}, mean_point = {0.0, 0.0, 0.0};
+      //Get mean
+      for (int vert = 0; vert<3; ++vert) { //vertice
+        auto point = vertices_[face[vert]];
+        auto E0 = field(charge, shift, point);
+        for (int i=0; i<3; ++i) {
+          mean_vector[i] += E0[i]/3.0;
+          mean_point[i] += point[i]/3.0;
+        }
+      }
+      // Vector to unit product
+      double r = norm(mean_point);
+      std::vector<double> unit = { mean_point[0]/r, mean_point[1]/r, mean_point[2]/r};
+      std::cout << norm(unit) << std::endl;
+      for (int i =0; i < 3; ++i) 
+        integral += face_area_*unit[i]*mean_vector[i];
+    }
+    return integral;
+  }
+
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  double ShellGenerator::Integrate() {
+    double integral = 0.0;
+    auto E = E_[0];
+    auto H = H_[0];
+    for(const auto face : faces_){
+      for (int i = 0; i<3; ++i) {
+        auto mean = (E_[face[0]][i]+E_[face[1]][i]+E_[face[2]][i])/3.0;
+        E[i] = mean;
+        mean = (H_[face[0]][i]+H_[face[1]][i]+H_[face[2]][i])/3.0;
+        H[i] = mean;
+      }
+      integral += // std::abs(E[0]*H[0])*
+        face_area_;
+    }
+    return integral;
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  std::vector< std::vector<double> > ShellGenerator::GetVerticesT() {
+    std::vector< std::vector<double> > vertices_t;
+    vertices_t.resize(3); 
+    for(const auto vert : vertices_){
+      vertices_t[0].push_back(vert[0]);
+      vertices_t[1].push_back(vert[1]);
+      vertices_t[2].push_back(vert[2]);
+    }
+    
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  double ShellGenerator::norm(std::vector<double> a){
+    double norm_value = 0;
+    for (auto coord:a)
+      norm_value += pow2(coord);
+    return std::sqrt(norm_value);
+   }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+   void ShellGenerator::Rescale(double scale) {
+     for(auto& vert : vertices_){
+       double factor = scale/norm(vert);
+       //std::cout<< factor <<std::endl;
+       for (auto &coord:vert) {
+         coord*=factor;
+       }
+     }
+     const double pi = 3.1415926535897932384626433832795;
+     face_area_ = 4.0*pi*pow2(scale)/faces_.size();
+   }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+   void ShellGenerator::PrintVerts() {
+   std::cout << "Verts coords:" << std::endl;
+    for(auto vert : vertices_){
+      std::cout <<"(";
+      for (auto coord:vert) std::cout<<coord<<",";
+      std::cout <<"),";      
+    }
+    std::cout << std::endl;
+   }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  void ShellGenerator::Refine() {
+    for (auto edge : edges_) {
+      auto p0 = vertices_[edge[0]];
+      auto p1 = vertices_[edge[1]];
+      std::vector<double> new_point = {
+        (p0[0]+p1[0])/2.0,
+        (p0[1]+p1[1])/2.0,
+        (p0[2]+p1[2])/2.0};
+      vertices_.push_back(new_point);      
+    }
+    GenerateEdges();
+    GenerateFaces();
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  void ShellGenerator::GenerateFaces() {
+    faces_.clear();
+    for (int i = 0; i < edges_.size(); ++i) {
+      const auto ie = edges_[i];
+      for(int j = i + 1; j < edges_.size(); ++j) {
+        const auto je = edges_[j];
+        for(int k = j + 1; k < edges_.size(); ++k) {
+          const auto ke = edges_[k];
+          std::set<long int> edges = {ie[0],ie[1],
+                                 je[0],je[1],
+                                 ke[0],ke[1]};
+          if (edges.size() != 3) continue;
+          std::vector<long int> face(edges.begin(), edges.end());
+          // std::cout << ie[0]<<"-"<<ie[1] << ":"
+          //             << je[0]<<"-"<<je[1] << ":"
+          //             << ke[0]<<"-"<<ke[1] 
+          //             << std::endl;
+          // std::cout << face[0]<<"-"<<face[1] << "-"<<face[2]<<std::endl;
+          faces_.push_back(face);
+        }
+      }
+    }
+    std::cout << "Faces: "<< faces_.size() <<std::endl;
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  void ShellGenerator::GenerateEdges() {
+    std::cout << "Vertices: "<< vertices_.size() <<std::endl;
+    edges_.clear();
+    EvalEdgeLength();
+    for (int i = 0; i < vertices_.size(); ++i)
+      for(int j = i + 1; j < vertices_.size(); ++j) {
+        if (dist(vertices_[i],vertices_[j]) > 1.001*edge_length_) continue;
+        std::vector<long int> edge = {i,j};
+        // std::cout << i<< " " << j<<std::endl;
+        edges_.push_back(edge);
+      }
+    std::cout << "Edges: "<< edges_.size() <<std::endl;
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  void ShellGenerator::EvalEdgeLength() {
+    auto p0 = vertices_[0];
+    std::vector<double> zero(p0.size(), 0.0);
+    double min_dist = 42.0*dist(zero, p0);
+    for (auto point : vertices_) {
+      if (point == p0) continue;
+      double new_dist = dist(p0, point);
+      if (new_dist < min_dist) min_dist = new_dist;
+    }
+    std::cout << "Edge length = " << min_dist << std::endl;
+    edge_length_ = min_dist;
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  double ShellGenerator::dist(std::vector<double> a, std::vector<double> b) {
+    unsigned int len = a.size();
+    if (b.size() != len)
+      throw std::invalid_argument("Error! Vector need to be the same size!\n");
+    double distance = 0.0;
+    for (int i = 0; i<len; ++i) {
+      distance += pow2(a[i]-b[i]);
+    }
+    return std::sqrt(distance);
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // @brief set up regular icosahedron
+  void ShellGenerator::SetInitialVertices() {
+    double a = 0.0;
+    double b = 1.0;
+    double c = (1+std::sqrt(5.0))/2.0;
+    std::vector< std::vector<double> > points = {
+      {a, b, c},
+      {a, b,-c},
+      {a,-b, c},
+      {a,-b,-c},
+      
+      { b, c,a},
+      { b,-c,a},
+      {-b, c,a},
+      {-b,-c,a},
+
+      { c,a, b},
+      {-c,a, b},
+      { c,a,-b},
+      {-c,a,-b}
+    };
+    vertices_ = std::move(points);
+    // for (auto v : vertices_) {
+    //   for (auto p : v)
+    //     std::cout<< p<<std::endl;
+    // }
+  }
+
+
+
+  // ShellGenerator& ShellGenerator::ReadFromFile(std::string fname) {
+  //   //std::cout<<"Reading file: "<< fname << std::endl;
+  //   std::ifstream infile(fname.c_str());
+  //   data_.clear();
+  //   std::string line;
+  //   while (std::getline(infile, line))
+  //     {
+  //       if (line.front() == '#') continue; //do not read comments
+  //       if (line.find('#') != std::string::npos) 
+  //         throw std::invalid_argument("Error! Comments should be marked with # in the begining of the line!\n");
+  //       std::istringstream iss(line);	
+  //       double wl, re, im;
+  //       if (!(iss >> wl >> re >> im)) throw std::invalid_argument("Error! Unexpected format of the line!\n");
+  //       data_.push_back(std::vector<double>({wl,re,im}));
+  //       //std::cout<<wl<<' '<<re<<' '<<im<<std::endl;
+  //     }  // end of wile reading file 
+  //   std::sort(data_.begin(), data_.end(),
+  //             [](const std::vector<double>& a, const std::vector<double>& b) {
+  //       	return a.front() < b.front();
+  //             });
+  //   return *this;
+  // }  // end of void ShellGenerator::ReadFromFile(std::string fname)
+  // // ********************************************************************** //
+  // // ********************************************************************** //
+  // // ********************************************************************** //
+  // /// Cut the spectra to the range and convert it to std::complex<double>
+  // ShellGenerator& ShellGenerator::ResizeToComplex(double from_wl, double to_wl, int samples) {
+  //   if (data_.size() < 2) throw std::invalid_argument("Nothing to resize!/n");
+  //   if (data_.front()[0] > from_wl || data_.front()[0] > to_wl ||
+  //       data_.back()[0] < from_wl || data_.back()[0] < to_wl ||
+  //       from_wl > to_wl)
+  //     throw std::invalid_argument("Invalid range to resize spectra!/n");
+  //   if (samples < 1) throw std::invalid_argument("Not enough samples!/n");
+  //   std::vector<double> wl_sampled(samples, 0.0);
+  //   if (samples == 1) {
+  //     wl_sampled[0] = (from_wl + to_wl)/2.0;
+  //   } else {
+  //     for (int i =0; i<samples; ++i)
+  //       wl_sampled[i] = from_wl
+  //         + (to_wl-from_wl)*static_cast<double>(i)/static_cast<double>(samples-1);
+  //   }  // end of setting wl_sampled
+  //   data_complex_.clear();
+  //   int j = 0;
+  //   for (int i = 0; i < data_.size(); ++i) {
+  //     const double& wl_i = data_[i][0];
+  //     const double& wl_s = wl_sampled[j];
+  //     if (wl_i < wl_s) continue;
+  //     else {
+  //       const double& wl_prev = data_[i-1][0];	
+  //       const double& re_prev = data_[i-1][1];
+  //       const double& im_prev = data_[i-1][2];
+  //       const double& re_i = data_[i][1];
+  //       const double& im_i = data_[i][2];
+  //       // Linear approximation
+  //       double re_s = re_i - (re_i-re_prev)*(wl_i-wl_s)/(wl_i-wl_prev);
+  //       double im_s = im_i - (im_i-im_prev)*(wl_i-wl_s)/(wl_i-wl_prev);
+
+  //       auto tmp = std::make_pair(wl_s, std::complex<double>(re_s,im_s));
+  //       data_complex_.push_back(tmp);	
+	       
+  //       ++j;
+  //       --i; // Next sampled point(j) can be in the same i .. i-1 region
+  //       // All sampled wavelengths has a value
+  //       if (j >= wl_sampled.size()) break;  
+  //     }
+  //   }
+  //   if (data_complex_.size() == 0)
+  //     throw std::invalid_argument("No points in spectra for the defined range!/n");
+  //   if (data_complex_.size() != samples)
+  //     throw std::invalid_argument("Was not able to get all samples!/n");
+  //   return *this;
+  // }
+  // // ********************************************************************** //
+  // // ********************************************************************** //
+  // // ********************************************************************** //
+  // /// from relative permittivity to refractive index
+  // ShellGenerator& ShellGenerator::ToIndex() {
+  //   data_complex_index_.clear();
+  //   for (auto row : data_complex_) {
+  //     const double wl = row.first;
+  //     const double e1 = row.second.real();
+  //     const double e2 = row.second.imag();
+  //     const double n = std::sqrt( (std::sqrt(pow2(e1)+pow2(e2)) + e1) /2.0 );
+  //     const double k = std::sqrt( (std::sqrt(pow2(e1)+pow2(e2)) - e1) /2.0 );
+  //     auto tmp = std::make_pair(wl, std::complex<double>(n,k));
+  //     data_complex_index_.push_back(tmp);	
+  //   }
+  //   return *this;
+  // }
+
+  // // ********************************************************************** //
+  // // ********************************************************************** //
+  // // ********************************************************************** //
+  // void ShellGenerator::PrintData() {
+  //   if (data_complex_.size() == 0) 
+  //     throw std::invalid_argument("Nothing to print!");
+  //   for (auto row : data_complex_) {
+  //     printf("wl:%g\tre:%g\tim:%g\n", row.first, row.second.real(),
+  //            row.second.imag());
+  //   }  // end of for each row
+  // }
+  // // ********************************************************************** //
+  // // ********************************************************************** //
+  // // ********************************************************************** //
+
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  void ShellGenerator::Init() {
+    SetInitialVertices();
+    GenerateEdges();
+    GenerateFaces();
+  }
+
+}  // end of namespace read_spectra
+

src/shell-generator.hpp

@@ -0,0 +1,70 @@
+#ifndef SRC_SHELL_GENERATOR_H_
+#define SRC_SHELL_GENERATOR_H_
+//**********************************************************************************//
+//    Copyright (C) 2009-2016  Ovidio Pena <ovidio@bytesfall.com>                   //
+//    Copyright (C) 2013-2016  Konstantin Ladutenko <kostyfisik@gmail.com>          //
+//                                                                                  //
+//    This file is part of scattnlay                                                //
+//                                                                                  //
+//    This program is free software: you can redistribute it and/or modify          //
+//    it under the terms of the GNU General Public License as published by          //
+//    the Free Software Foundation, either version 3 of the License, or             //
+//    (at your option) any later version.                                           //
+//                                                                                  //
+//    This program is distributed in the hope that it will be useful,               //
+//    but WITHOUT ANY WARRANTY; without even the implied warranty of                //
+//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                 //
+//    GNU General Public License for more details.                                  //
+//                                                                                  //
+//    The only additional remark is that we expect that all publications            //
+//    describing work using this software, or all commercial products               //
+//    using it, cite 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 <http://www.gnu.org/licenses/>.         //
+//**********************************************************************************//
+//   @brief  Generates points for integration on sphere surface
+#include <complex>
+#include <string>
+#include <utility>
+#include <vector>
+namespace shell_generator {
+  class ShellGenerator {  // will throw for any error
+   public:
+    /* ShellGenerator& ReadFromFile(std::string filename); */
+    /* ShellGenerator& ResizeToComplex(double from_wl, double to_wl, int samples); */
+    /* ShellGenerator& ToIndex(); */
+    double dist(std::vector<double> a, std::vector<double> b);
+    double norm(std::vector<double> a);
+    std::vector< std::vector<double> > GetVertices(){return vertices_;};
+    std::vector< std::vector<double> > GetVerticesT();
+    std::vector< std::vector<long int> > GetEdges(){return edges_;};
+    std::vector< std::vector<long int> > GetFaces(){return faces_;};
+    void EvalEdgeLength();
+    void GenerateEdges();
+    void GenerateFaces();
+    void Init();
+    double Integrate();
+    double IntegrateGauss(double charge, double dist);
+    void PrintVerts();
+    void Refine();
+    void Rescale(double scale);
+    void SetField(std::vector<std::vector< std::complex<double> > > &E,
+                  std::vector<std::vector< std::complex<double> > > &H) {E_ = E; H_=H;};
+    void SetInitialVertices();
+  private:
+    std::vector<std::vector< std::complex<double> > > E_, H_;
+    double edge_length_ = 0.0;
+    double face_area_ = 0.0;
+    std::vector< std::vector<double> > vertices_;
+    std::vector< std::vector<long int> > edges_;
+    std::vector< std::vector<long int> > faces_;
+    // std::vector< std::pair< double, std::complex<double> > > data_complex_;
+    // std::vector< std::pair< double, std::complex<double> > > data_complex_index_;
+    // void PermittivityToIndex();
+  };  // end of class ShellGenerator
+}  // end of namespase read_spectra
+#endif  // SRC_SHELL_GENERATOR_H_