Gauss integration verified vs Gauss law

.gitignore

@@ -41,3 +41,5 @@ cov-int
 build*
 .ipynb_checkpoints
 tests/python/field.txt
+examples/*.pdf
+examples/*.png

examples/go-cc-examples.sh

@@ -2,7 +2,15 @@
 path=$PWD
 PROGRAM='scattnlay-example.bin'
 
-file=example-eval-force.cc
+# file=example-eval-force.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=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
@@ -10,6 +18,7 @@ g++ -Ofast -std=c++11 $file  ../src/shell-generator.cc -lm -lrt -o $PROGRAM -mar
 echo Compilation done. Running...
 ./$PROGRAM
 
+
 # file=example-minimal.cc
 # echo Compile $file with gcc
 # rm -f $PROGRAM

examples/test-surf-integral.cc

@@ -0,0 +1,56 @@
+//**********************************************************************************//
+//    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/>.         //
+//**********************************************************************************//
+//   This program evaluates forces acting on the nanoparticle under irradiaton.
+#include <complex>
+#include <cstdio>
+#include <string>
+#include <iostream>
+#include "../src/shell-generator.hpp"
+int main(int argc, char *argv[]) {
+  try {
+    shell_generator::ShellGenerator shell;
+    shell.Init();
+    shell.Refine();
+    //shell.Refine();
+    //shell.Refine();
+    double scale = 1.9;  //Integration sphere radius.
+    double shift = 1.1;  //Shift of the charge relative to the sphere center.
+    double charge = 11.0; //Coulomb charge.
+    shell.Rescale(scale);
+    shell.PrintVerts();
+    auto points = shell.GetVertices();
+    double charge_s = shell.IntegrateGaussSimple(charge, shift);
+    std::cout << "Accuracy ( 1==ideal ): " << charge_s/charge << std::endl; 
+  } 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;
+}
+
+

src/shell-generator.cc

@@ -40,19 +40,43 @@
 #include "shell-generator.hpp"
 namespace shell_generator {
   template<class T> inline T pow2(const T value) {return value*value;}
-
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  double ShellGenerator::IntegrateGaussSimple(double charge, double shift) {
+    double integral = 0.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]) );
+      const double pi = 3.1415926535897932384626433832795;
+      double ampl = charge/(4.0*pi*pow2(r));      
+      std::vector<double> field = {ampl*(p[0]-shift)/r, ampl*(p[1])/r, ampl*(p[2])/r};
+      return field;
+    };
+    //simple 
+    for (auto vert :vertices_) {
+      auto E0 = field(charge, shift, vert);
+      // Vector to unit product
+      double r = norm(vert);
+      std::vector<double> unit = { vert[0]/r, vert[1]/r, vert[2]/r};
+      // std::cout << norm(unit) << std::endl;
+      for (int i =0; i < 3; ++i) 
+        integral += per_vertice_area_*unit[i]*E0[i];
+    }
+    return integral;
+  }
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
   double ShellGenerator::IntegrateGauss(double charge, double shift) {
+    if (faces_.size() == 0)
+      throw std::invalid_argument("Error! Faces were not initialized!\n");
     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);
+      const double pi = 3.1415926535897932384626433832795;
+      double ampl = charge/(4.0*pi*pow2(r));      
       std::vector<double> field = {ampl*(p[0]-shift)/r, ampl*(p[1])/r, ampl*(p[2])/r};
       return field;
     };
@@ -71,9 +95,10 @@ namespace shell_generator {
       // 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;
+      // std::cout << norm(unit) << std::endl;
       for (int i =0; i < 3; ++i) 
-        integral += face_area_*unit[i]*mean_vector[i];
+        integral += face_area_*
+          unit[i]*mean_vector[i];
     }
     return integral;
   }
@@ -131,7 +156,9 @@ namespace shell_generator {
        }
      }
      const double pi = 3.1415926535897932384626433832795;
-     face_area_ = 4.0*pi*pow2(scale)/faces_.size();
+     double area = 4.0*pi*pow2(scale); 
+     face_area_ = area/faces_.size();
+     per_vertice_area_ = area/vertices_.size(); 
    }
   // ********************************************************************** //
   // ********************************************************************** //
@@ -159,7 +186,7 @@ namespace shell_generator {
       vertices_.push_back(new_point);      
     }
     GenerateEdges();
-    GenerateFaces();
+    // GenerateFaces();
   }
   // ********************************************************************** //
   // ********************************************************************** //
@@ -262,121 +289,13 @@ namespace shell_generator {
     //     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();
+    //GenerateFaces();
   }
 
 }  // end of namespace read_spectra

src/shell-generator.hpp

@@ -49,6 +49,7 @@ namespace shell_generator {
     void Init();
     double Integrate();
     double IntegrateGauss(double charge, double dist);
+    double IntegrateGaussSimple(double charge, double dist);
     void PrintVerts();
     void Refine();
     void Rescale(double scale);
@@ -59,6 +60,7 @@ namespace shell_generator {
     std::vector<std::vector< std::complex<double> > > E_, H_;
     double edge_length_ = 0.0;
     double face_area_ = 0.0;
+    double per_vertice_area_ = 0.0;
     std::vector< std::vector<double> > vertices_;
     std::vector< std::vector<long int> > edges_;
     std::vector< std::vector<long int> > faces_;