update doc

README.md

@@ -78,24 +78,33 @@ fieldnlay -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] -p xi xf nx yi yf ny zi zf n
   
 3. C++ library
 
+Scattnlay "Hello world!" example:
+
 ```C++
     try {
-      MultiLayerMie multi_layer_mie;
-      multi_layer_mie.SetLayersSize(x);
-      multi_layer_mie.SetLayersIndex(m);
-
+      nmie::MultiLayerMieApplied<double> multi_layer_mie;  
+      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();
-
-      *Qsca = multi_layer_mie.GetQsca();
-      *Qabs = multi_layer_mie.GetQabs();
-    } catch(const std::invalid_argument& ia) {
+      double Qabs = multi_layer_mie.GetQabs();
+      printf("Qabs = %g\n", Qabs);
+    } catch( const std::invalid_argument& ia ) {
       // Will catch if  multi_layer_mie fails or other errors.
       std::cerr << "Invalid argument: " << ia.what() << std::endl;
-      throw std::invalid_argument(ia);
       return -1;
     }
 ```
 
+The complete `example-minimal.cc` and a bit more complicated
+`example-get-Mie.cc` can be found in example directory along with
+`go-cc-examples.sh` script with build commands.
+
+`example-get-Mie.cc` can be compiled using double precision or
+multiple precision (just include `-DMULTI_PRECISION=200` to use 200
+digits for calculations). 
+
 Papers
 ------
 

examples/example-minimal.cc

@@ -0,0 +1,63 @@
+//**********************************************************************************//
+//    Copyright (C) 2009-2015  Ovidio Pena <ovidio@bytesfall.com>                   //
+//    Copyright (C) 2013-2015  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 absorption of a triple layered nanoparticle
+#include <complex>
+#include <cstdio>
+#include <string>
+#include "../src/nmie-applied.hpp"
+#include "../src/nmie-applied-impl.hpp"
+
+int main(int argc, char *argv[]) {
+  try {
+    nmie::MultiLayerMieApplied<double> multi_layer_mie;  
+    const std::complex<double> epsilon_Si(18.4631066585, 0.6259727805);
+    const std::complex<double> epsilon_Ag(-8.5014154589, 0.7585845411);
+    const std::complex<double> index_Si = std::sqrt(epsilon_Si);
+    const std::complex<double> 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
+    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 ) {
+    // Will catch if  multi_layer_mie fails or other errors.
+    std::cerr << "Invalid argument: " << ia.what() << std::endl;
+    return -1;
+  }  
+    return 0;
+}
+
+

examples/go-cc-examples.sh

@@ -2,15 +2,24 @@
 path=$PWD
 PROGRAM='scattnlay-example.bin'
 
-echo Compile with gcc
+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
-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
-#g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc read-spectra.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
+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
 
-# g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc ./read-spectra.cc -lm -lrt -o $PROGRAM -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2
-echo Compilation done.
+# 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
-# #result

examples/read-spectra.cc

@@ -0,0 +1,146 @@
+/**
+ * @file   read-spectra.cc
+ * @author Konstantin Ladutenko <kostyfisik at gmail (.) com>
+ * @date   Wed Mar 11 11:51:26 2015
+ * 
+ * @copyright 2015 Konstantin Ladutenko
+ *
+ * @brief  Read complex spectra from file in format 'WL real imag'
+ * 
+ * read-spectra 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.
+ *
+ * read-spectra 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 read-spectra.  If not, see <http://www.gnu.org/licenses/>.
+ * 
+ */
+
+#include <algorithm>
+#include <complex>
+#include <cstdio>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <stdexcept>
+#include <iostream>
+#include <vector>
+#include "read-spectra.h"
+namespace read_spectra {
+  template<class T> inline T pow2(const T value) {return value*value;}
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  ReadSpectra& ReadSpectra::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 ReadSpectra::ReadFromFile(std::string fname)
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  /// Cut the spectra to the range and convert it to std::complex<double>
+  ReadSpectra& ReadSpectra::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
+  ReadSpectra& ReadSpectra::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 ReadSpectra::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
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+
+}  // end of namespace read_spectra
+