try forces

examples/example-eval-force.cc

@@ -77,13 +77,14 @@ int main(int argc, char *argv[]) {
         // auto Hs = nmie::ConvertComplexVectorVector<double>(multi_layer_mie.GetFieldHs());
         shell.SetField(E,H);
         //shell.SetFieldSph(Es,Hs);
-        // auto F = shell.Integrate();
+        //auto F = shell.Integrate();
         //auto F = shell.IntegrateByFaces();
-        //auto F = shell.IntegrateByComp();
-        //        auto F = shell.IntegrateByCompReal();
-        //std::cout << "integrate_R:\t" << scale*WL/(2.0*pi);
-        //std::cout<<"\tforce:\t" <<F[0]<<"\t"<< F[1] <<"\t"<<F[2] << std::endl;
+        // auto F = shell.IntegrateByComp();
+        // //        auto F = shell.IntegrateByCompReal();
+        // std::cout << "integrate_R:\t" << scale*WL/(2.0*pi);
+        // std::cout<<"\tforce:\t" <<F[0]<<"\t"<< F[1] <<"\t"<<F[2] << std::endl;
 
+        // std::cout<<"clean"<<std::endl;
         auto F = shell.IntegrateGauss(2.54,12.03);
         //std::cout<<"\tcharge:\t" <<F<< std::endl;
         

examples/force-quad.py

@@ -0,0 +1,101 @@
+#!/usr/bin/env python3
+# -*- coding: UTF-8 -*-
+#
+#    Copyright (C) 2018  Konstantin Ladutenko <kostyfisik@gmail.com>
+#
+#    This file is part of python-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/>.
+
+# WIP: try to evaluate forces.
+
+from scattnlay import fieldnlay
+from scattnlay import scattnlay
+import cmath
+import matplotlib.pyplot as plt
+import numpy as np
+import scattnlay
+import quadpy
+
+
+index_Ag = 4.0
+
+nm = 1.0
+WL_units='nm'
+
+x = np.ones((1), dtype = np.float64)
+m = np.ones((1), dtype = np.complex128)
+
+core_r = 120
+WL = 550 
+
+x[0] = 2.0*np.pi*core_r/WL#/4.0*3.0
+m[0] = index_Ag/nm
+
+R = x[0]*3.01
+
+comment='bulk-NP-WL'+str(WL)+WL_units
+
+quad_ord = 19
+#quad_ord = 131
+
+coord = quadpy.sphere.Lebedev(3).points
+# coord = np.vstack((coordX, coordY, coordZ)).transpose()
+def force(in_coord):
+    coord = in_coord.T
+    terms, Eall, Hall = fieldnlay(np.array([x]), np.array([m]), coord)
+    E_all = Eall[0, :, :]
+    H_all = Hall[0, :, :]
+      # std::vector<double> P = (1/(2.0))
+      #   *real(
+      #         dot(unit,E)*vconj(E) +
+      #         dot(unit,H)*vconj(H) +
+      #         (-1.0/2.0)*(dot(E,vconj(E))
+      #                     +dot(H,vconj(H))
+      #                     )*unit
+      #         );
+    unit_all = coord/ R
+    P = np.array([
+        ( (1/(2.0))
+        *np.real(
+            np.dot(unit,E)*np.conj(E) +
+            np.dot(unit,H)*np.conj(H) +
+              (-1.0/2.0)*(np.dot(E,np.conj(E))
+                          +np.dot(H,np.conj(H))
+                          )*unit
+              )
+        )
+        for unit,E,H in zip(unit_all,E_all, H_all)
+        ])
+    return P.T
+    
+# P = np.array(map(lambda n: np.linalg.norm(np.cross(Ec[n], Hc[n])).real,
+#                      range(0, len(E[0]))))
+
+val = quadpy.sphere.integrate(
+    force,
+    [0.0, 0.0, 0.0], R,
+    quadpy.sphere.Lebedev(quad_ord)
+    )
+
+
+
+print(val)

examples/go-cc-examples.sh

@@ -2,13 +2,13 @@
 path=$PWD
 PROGRAM='scattnlay-example.bin'
 
-# 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 -Wall
-# # 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...
-# time ./$PROGRAM
+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 -Wall
+# 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...
+time ./$PROGRAM
 
 # file=test-surf-integral.cc
 # echo Compile $file with gcc
@@ -27,16 +27,16 @@ PROGRAM='scattnlay-example.bin'
 # ./$PROGRAM
 
 
-file=example-get-Mie.cc
-echo Compile $file with gcc
-rm -f $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
 
-# Simplified for double precision
-clang++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc ./read-spectra.cc -lm -o $PROGRAM -march=native -mtune=native -msse4.2
-echo Compilation done. Running...
-./$PROGRAM
+# # Simplified for double precision
+# clang++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc ./read-spectra.cc -lm -o $PROGRAM -march=native -mtune=native -msse4.2
+# echo Compilation done. Running...
+# ./$PROGRAM