c++ python nMie

go.sh

@@ -94,7 +94,10 @@ PROGRAM='../../../scattnlay'
 # ./$file.bin
 
 rm scattnlay.so
+export CFLAGS='-std=c++11'
 python setup.py build_ext --inplace
 cp scattnlay.so tests/python/
 cd tests/python/
 ./field.py
+./test01.py
+./test01-wrapper.py

nmie-wrapper.cc

@@ -48,7 +48,7 @@ namespace nmie {
   // ********************************************************************** //
   // ********************************************************************** //
   //emulate C call.
-  int nMie_wrapper(int L, const std::vector<double>& x, const std::vector<std::complex<double> >& m, int nTheta, const std::vector<double>& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
+  int nMie_wrapper(int L,  std::vector<double>& x,  std::vector<std::complex<double> >& m, int nTheta, std::vector<double>& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
     
     if (x.size() != L || m.size() != L)
         throw std::invalid_argument("Declared number of layers do not fit x and m!");
@@ -1267,8 +1267,8 @@ c    MM + 1  and - 1, alternately
 
     isMieCalculated_ = true;
     nmax_used_ = nmax_;
-    printf("Run Mie result: Qext = %g, Qsca = %g, Qabs = %g, Qbk = %g \n",
-               GetQext(), GetQsca(), GetQabs(), GetQbk());
+    // printf("Run Mie result: Qext = %g, Qsca = %g, Qabs = %g, Qbk = %g \n",
+    //            GetQext(), GetQsca(), GetQabs(), GetQbk());
     //return nmax;
   }
   

nmie-wrapper.h

@@ -54,7 +54,7 @@
 
 namespace nmie {
 
-  int nMie_wrapper(int L, const std::vector<double>& x, const std::vector<std::complex<double> >& m, int nTheta, const std::vector<double>& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2);
+  int nMie_wrapper(int L,  std::vector<double>& x,  std::vector<std::complex<double> >& m, int nTheta,  std::vector<double>& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2);
   int nField(const int L, const int pl, const std::vector<double>& x, const std::vector<std::complex<double> >& m, const int nmax, const int ncoord, const std::vector<double>& Xp, const std::vector<double>& Yp, const std::vector<double>& Zp,  std::vector<std::vector<std::complex<double> > >& E, std::vector<std::vector<std::complex<double> > >& H);
 
 

scattnlay.pyx

@@ -21,7 +21,8 @@
 #
 #    You should have received a copy of the GNU General Public License
 #    along with this program.  If not, see <http://www.gnu.org/licenses/>.
-
+# distutils: language = c++
+# distutils: sources = nmie-wrapper.cc
 from __future__ import division
 import numpy as np
 cimport numpy as np
@@ -50,12 +51,112 @@ cdef inline double *npy2c(np.ndarray a):
 
     # Return data pointer
     return <double *>(a.data)
-
-cdef extern from "py_nmie.h":
-    cdef int nMie(int L, int pl, vector[double] x, vector[complex] m, int nTheta, vector[double] Theta, int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, double S1r[], double S1i[], double S2r[], double S2i[])
-    cdef int nField(int L, int pl, vector[double] x, vector[complex] m, int nmax, int nCoords, vector[double] Xp, vector[double] Yp, vector[double] Zp, double Erx[], double Ery[], double Erz[], double Eix[], double Eiy[], double Eiz[], double Hrx[], double Hry[], double Hrz[], double Hix[], double Hiy[], double Hiz[])
-
-def scattnlay(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t, ndim = 2] m, np.ndarray[np.float64_t, ndim = 1] theta = np.zeros(0, dtype = np.float64), np.int_t pl = -1, np.int_t nmax = -1):
+##############################################################################
+##############################################################################
+##############################################################################
+##############################################################################
+
+# cdef extern from "py_nmie.h":
+#     cdef int nMie(int L, int pl, vector[double] x, vector[double complex] m, int nTheta, vector[double] Theta, int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, double S1r[], double S1i[], double S2r[], double S2i[])
+#     cdef int nField(int L, int pl, vector[double] x, vector[double complex] m, int nmax, int nCoords, vector[double] Xp, vector[double] Yp, vector[double] Zp, double Erx[], double Ery[], double Erz[], double Eix[], double Eiy[], double Eiz[], double Hrx[], double Hry[], double Hrz[], double Hix[], double Hiy[], double Hiz[])
+##############################################################################
+##############################################################################
+##############################################################################
+##############################################################################
+
+cdef extern from "nmie-wrapper.h"  namespace "nmie":
+    cdef int nMie_wrapper(int L, const vector[double] x, const vector[double complex] m , int nTheta, const vector[double] Theta, double *qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, vector[double complex] S1, vector[double complex] S2);
+
+
+##############################################################################
+##############################################################################
+##############################################################################
+##############################################################################
+# def scattnlay(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t, ndim = 2] m, np.ndarray[np.float64_t, ndim = 1] theta = np.zeros(0, dtype = np.float64), np.int_t pl = -1, np.int_t nmax = -1):
+#     cdef Py_ssize_t i
+
+#     cdef np.ndarray[np.int_t, ndim = 1] terms = np.zeros(x.shape[0], dtype = np.int)
+
+#     cdef np.ndarray[np.float64_t, ndim = 1] Qext = np.zeros(x.shape[0], dtype = np.float64)
+#     cdef np.ndarray[np.float64_t, ndim = 1] Qabs = np.zeros(x.shape[0], dtype = np.float64)
+#     cdef np.ndarray[np.float64_t, ndim = 1] Qsca = np.zeros(x.shape[0], dtype = np.float64)
+#     cdef np.ndarray[np.float64_t, ndim = 1] Qbk = np.zeros(x.shape[0], dtype = np.float64)
+#     cdef np.ndarray[np.float64_t, ndim = 1] Qpr = np.zeros(x.shape[0], dtype = np.float64)
+#     cdef np.ndarray[np.float64_t, ndim = 1] g = np.zeros(x.shape[0], dtype = np.float64)
+#     cdef np.ndarray[np.float64_t, ndim = 1] Albedo = np.zeros(x.shape[0], dtype = np.float64)
+
+#     cdef np.ndarray[np.complex128_t, ndim = 2] S1 = np.zeros((x.shape[0], theta.shape[0]), dtype = np.complex128)
+#     cdef np.ndarray[np.complex128_t, ndim = 2] S2 = np.zeros((x.shape[0], theta.shape[0]), dtype = np.complex128)
+
+#     cdef np.ndarray[np.float64_t, ndim = 1] S1r
+#     cdef np.ndarray[np.float64_t, ndim = 1] S1i
+#     cdef np.ndarray[np.float64_t, ndim = 1] S2r
+#     cdef np.ndarray[np.float64_t, ndim = 1] S2i
+
+#     for i in range(x.shape[0]):
+#         S1r = np.zeros(theta.shape[0], dtype = np.float64)
+#         S1i = np.zeros(theta.shape[0], dtype = np.float64)
+#         S2r = np.zeros(theta.shape[0], dtype = np.float64)
+#         S2i = np.zeros(theta.shape[0], dtype = np.float64)
+
+#         terms[i] = nMie(x.shape[1], pl, x[i].copy('C'), m[i].copy('C'), theta.shape[0], theta.copy('C'), nmax, &Qext[i], &Qsca[i], &Qabs[i], &Qbk[i], &Qpr[i], &g[i], &Albedo[i], npy2c(S1r), npy2c(S1i), npy2c(S2r), npy2c(S2i))
+
+#         S1[i] = S1r.copy('C') + 1.0j*S1i.copy('C')
+#         S2[i] = S2r.copy('C') + 1.0j*S2i.copy('C')
+
+#     return terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2
+# ##############################################################################
+# ##############################################################################
+# ##############################################################################
+# ##############################################################################
+
+# #def fieldnlay(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t, ndim = 2] m, np.ndarray[np.float64_t, ndim = 2] coords = np.zeros((0, 3), dtype = np.float64), np.int_t pl = 0, np.int_t nmax = 0):
+# def fieldnlay(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t, ndim = 2] m, np.ndarray[np.float64_t, ndim = 2] coords, np.int_t pl = 0, np.int_t nmax = 0):
+#     cdef Py_ssize_t i
+
+#     cdef np.ndarray[np.int_t, ndim = 1] terms = np.zeros(x.shape[0], dtype = np.int)
+
+#     cdef np.ndarray[np.complex128_t, ndim = 3] E = np.zeros((x.shape[0], coords.shape[0], 3), dtype = np.complex128)
+#     cdef np.ndarray[np.complex128_t, ndim = 3] H = np.zeros((x.shape[0], coords.shape[0], 3), dtype = np.complex128)
+
+#     cdef np.ndarray[np.float64_t, ndim = 1] Erx
+#     cdef np.ndarray[np.float64_t, ndim = 1] Ery
+#     cdef np.ndarray[np.float64_t, ndim = 1] Erz
+#     cdef np.ndarray[np.float64_t, ndim = 1] Eix
+#     cdef np.ndarray[np.float64_t, ndim = 1] Eiy
+#     cdef np.ndarray[np.float64_t, ndim = 1] Eiz
+#     cdef np.ndarray[np.float64_t, ndim = 1] Hrx
+#     cdef np.ndarray[np.float64_t, ndim = 1] Hry
+#     cdef np.ndarray[np.float64_t, ndim = 1] Hrz
+#     cdef np.ndarray[np.float64_t, ndim = 1] Hix
+#     cdef np.ndarray[np.float64_t, ndim = 1] Hiy
+#     cdef np.ndarray[np.float64_t, ndim = 1] Hiz
+
+#     for i in range(x.shape[0]):
+#         Erx = np.zeros(coords.shape[0], dtype = np.float64)
+#         Ery = np.zeros(coords.shape[0], dtype = np.float64)
+#         Erz = np.zeros(coords.shape[0], dtype = np.float64)
+#         Eix = np.zeros(coords.shape[0], dtype = np.float64)
+#         Eiy = np.zeros(coords.shape[0], dtype = np.float64)
+#         Eiz = np.zeros(coords.shape[0], dtype = np.float64)
+#         Hrx = np.zeros(coords.shape[0], dtype = np.float64)
+#         Hry = np.zeros(coords.shape[0], dtype = np.float64)
+#         Hrz = np.zeros(coords.shape[0], dtype = np.float64)
+#         Hix = np.zeros(coords.shape[0], dtype = np.float64)
+#         Hiy = np.zeros(coords.shape[0], dtype = np.float64)
+#         Hiz = np.zeros(coords.shape[0], dtype = np.float64)
+
+#         terms[i] = nField(x.shape[1], pl, x[i].copy('C'), m[i].copy('C'), nmax, coords.shape[0], coords[:, 0].copy('C'), coords[:, 1].copy('C'), coords[:, 2].copy('C'), npy2c(Erx), npy2c(Ery), npy2c(Erz), npy2c(Eix), npy2c(Eiy), npy2c(Eiz), npy2c(Hrx), npy2c(Hry), npy2c(Hrz), npy2c(Hix), npy2c(Hiy), npy2c(Hiz))
+
+#         E[i] = np.vstack((Erx.copy('C') + 1.0j*Eix.copy('C'), Ery.copy('C') + 1.0j*Eiy.copy('C'), Erz.copy('C') + 1.0j*Eiz.copy('C'))).transpose()
+#         H[i] = np.vstack((Hrx.copy('C') + 1.0j*Hix.copy('C'), Hry.copy('C') + 1.0j*Hiy.copy('C'), Hrz.copy('C') + 1.0j*Hiz.copy('C'))).transpose()
+
+#     return terms, E, H
+##############################################################################
+##############################################################################
+##############################################################################
+##############################################################################
+def scattnlay_wrapper(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t, ndim = 2] m, np.ndarray[np.float64_t, ndim = 1] theta = np.zeros(0, dtype = np.float64), np.int_t pl = -1, np.int_t nmax = -1):
     cdef Py_ssize_t i
 
     cdef np.ndarray[np.int_t, ndim = 1] terms = np.zeros(x.shape[0], dtype = np.int)
@@ -82,53 +183,15 @@ def scattnlay(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t,
         S2r = np.zeros(theta.shape[0], dtype = np.float64)
         S2i = np.zeros(theta.shape[0], dtype = np.float64)
 
-        terms[i] = nMie(x.shape[1], pl, x[i].copy('C'), m[i].copy('C'), theta.shape[0], theta.copy('C'), nmax, &Qext[i], &Qsca[i], &Qabs[i], &Qbk[i], &Qpr[i], &g[i], &Albedo[i], npy2c(S1r), npy2c(S1i), npy2c(S2r), npy2c(S2i))
+        terms[i] = nMie_wrapper(x.shape[1], x[i].copy('C'),
+                                m[i].copy('C'), theta.shape[0],
+                                theta.copy('C'), &Qext[i], &Qsca[i],
+                                &Qabs[i], &Qbk[i], &Qpr[i], &g[i],
+                                &Albedo[i],
+                                S1r, S2r)
 
-        S1[i] = S1r.copy('C') + 1.0j*S1i.copy('C')
-        S2[i] = S2r.copy('C') + 1.0j*S2i.copy('C')
+        S1[i] = S1r.copy('C')
+        S2[i] = S2r.copy('C')
 
     return terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2
 
-#def fieldnlay(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t, ndim = 2] m, np.ndarray[np.float64_t, ndim = 2] coords = np.zeros((0, 3), dtype = np.float64), np.int_t pl = 0, np.int_t nmax = 0):
-def fieldnlay(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t, ndim = 2] m, np.ndarray[np.float64_t, ndim = 2] coords, np.int_t pl = 0, np.int_t nmax = 0):
-    cdef Py_ssize_t i
-
-    cdef np.ndarray[np.int_t, ndim = 1] terms = np.zeros(x.shape[0], dtype = np.int)
-
-    cdef np.ndarray[np.complex128_t, ndim = 3] E = np.zeros((x.shape[0], coords.shape[0], 3), dtype = np.complex128)
-    cdef np.ndarray[np.complex128_t, ndim = 3] H = np.zeros((x.shape[0], coords.shape[0], 3), dtype = np.complex128)
-
-    cdef np.ndarray[np.float64_t, ndim = 1] Erx
-    cdef np.ndarray[np.float64_t, ndim = 1] Ery
-    cdef np.ndarray[np.float64_t, ndim = 1] Erz
-    cdef np.ndarray[np.float64_t, ndim = 1] Eix
-    cdef np.ndarray[np.float64_t, ndim = 1] Eiy
-    cdef np.ndarray[np.float64_t, ndim = 1] Eiz
-    cdef np.ndarray[np.float64_t, ndim = 1] Hrx
-    cdef np.ndarray[np.float64_t, ndim = 1] Hry
-    cdef np.ndarray[np.float64_t, ndim = 1] Hrz
-    cdef np.ndarray[np.float64_t, ndim = 1] Hix
-    cdef np.ndarray[np.float64_t, ndim = 1] Hiy
-    cdef np.ndarray[np.float64_t, ndim = 1] Hiz
-
-    for i in range(x.shape[0]):
-        Erx = np.zeros(coords.shape[0], dtype = np.float64)
-        Ery = np.zeros(coords.shape[0], dtype = np.float64)
-        Erz = np.zeros(coords.shape[0], dtype = np.float64)
-        Eix = np.zeros(coords.shape[0], dtype = np.float64)
-        Eiy = np.zeros(coords.shape[0], dtype = np.float64)
-        Eiz = np.zeros(coords.shape[0], dtype = np.float64)
-        Hrx = np.zeros(coords.shape[0], dtype = np.float64)
-        Hry = np.zeros(coords.shape[0], dtype = np.float64)
-        Hrz = np.zeros(coords.shape[0], dtype = np.float64)
-        Hix = np.zeros(coords.shape[0], dtype = np.float64)
-        Hiy = np.zeros(coords.shape[0], dtype = np.float64)
-        Hiz = np.zeros(coords.shape[0], dtype = np.float64)
-
-        terms[i] = nField(x.shape[1], pl, x[i].copy('C'), m[i].copy('C'), nmax, coords.shape[0], coords[:, 0].copy('C'), coords[:, 1].copy('C'), coords[:, 2].copy('C'), npy2c(Erx), npy2c(Ery), npy2c(Erz), npy2c(Eix), npy2c(Eiy), npy2c(Eiz), npy2c(Hrx), npy2c(Hry), npy2c(Hrz), npy2c(Hix), npy2c(Hiy), npy2c(Hiz))
-
-        E[i] = np.vstack((Erx.copy('C') + 1.0j*Eix.copy('C'), Ery.copy('C') + 1.0j*Eiy.copy('C'), Erz.copy('C') + 1.0j*Eiz.copy('C'))).transpose()
-        H[i] = np.vstack((Hrx.copy('C') + 1.0j*Hix.copy('C'), Hry.copy('C') + 1.0j*Hiy.copy('C'), Hrz.copy('C') + 1.0j*Hiz.copy('C'))).transpose()
-
-    return terms, E, H
-

setup.py

@@ -36,23 +36,31 @@ from distutils.core import setup
 from distutils.extension import Extension
 import numpy as np
 
-setup(name = __mod__,
-      version = __version__,
-      description = __title__,
-      long_description="""The Python version of scattnlay, a computer implementation of the algorithm for the calculation of electromagnetic \
-radiation scattering by a multilayered sphere developed by Yang. It has been shown that the program is effective, \
-resulting in very accurate values of scattering efficiencies for a wide range of size parameters, which is a \
-considerable improvement over previous implementations of similar algorithms. For details see: \
-O. Pena, U. Pal, Comput. Phys. Commun. 180 (2009) 2348-2354.""",
-      author = __author__,
-      author_email = __email__,
-      maintainer = __author__,
-      maintainer_email = __email__,
-      keywords = ['Mie scattering', 'Multilayered sphere', 'Efficiency factors', 'Cross-sections'],
-      url = __url__,
-      license = 'GPL',
-      platforms = 'any',
-      ext_modules = [Extension("scattnlay", ["nmie.cc", "py_nmie.cc", "scattnlay.cc"], language = "c++", include_dirs = [np.get_include()])], 
-      extra_compile_args=['-std=c++11']
-)
+# setup(name = __mod__,
+#       version = __version__,
+#       description = __title__,
+#       long_description="""The Python version of scattnlay, a computer implementation of the algorithm for the calculation of electromagnetic \
+# radiation scattering by a multilayered sphere developed by Yang. It has been shown that the program is effective, \
+# resulting in very accurate values of scattering efficiencies for a wide range of size parameters, which is a \
+# considerable improvement over previous implementations of similar algorithms. For details see: \
+# O. Pena, U. Pal, Comput. Phys. Commun. 180 (2009) 2348-2354.""",
+#       author = __author__,
+#       author_email = __email__,
+#       maintainer = __author__,
+#       maintainer_email = __email__,
+#       keywords = ['Mie scattering', 'Multilayered sphere', 'Efficiency factors', 'Cross-sections'],
+#       url = __url__,
+#       license = 'GPL',
+#       platforms = 'any',
+#       ext_modules = [Extension("scattnlay", ["nmie.cc", "nmie-wrapper.cc", "py_nmie.cc", "scattnlay.cc"], language = "c++", include_dirs = [np.get_include()])], 
+#       extra_compile_args=['-std=c++11']
+# )
 
+from distutils.core import setup
+from Cython.Build import cythonize
+setup(ext_modules = cythonize(
+           "scattnlay.pyx",                 # our Cython source
+           sources=["nmie-wrapper.cc"],  # additional source file(s)
+           language="c++",             # generate C++ code
+           extra_compile_args=['-std=c++11'],
+      ))

tests/python/field.py

@@ -33,6 +33,11 @@
 # 400 nm. Maximum of the surface plasmon resonance (and,
 # hence, of electric field) is expected under those
 # conditions.
+import scattnlay
+
+import os
+path = os.path.dirname(scattnlay.__file__)
+print(scattnlay.__file__)
 
 from scattnlay import fieldnlay
 import numpy as np

tests/python/test01.py

@@ -47,9 +47,17 @@
 # v2/Vt=0.26
 # v3/Vt=0.044
 # v4/Vt=0.3666
+import scattnlay
+
+import os
+path = os.path.dirname(scattnlay.__file__)
+print(scattnlay.__file__)
 
 from scattnlay import scattnlay
 import numpy as np
+# import os
+# import inspect
+# inspect.getfile(scattnlay)
 
 x = np.ones((400, 5), dtype = np.float64)
 x[:, 4] = np.arange(0.25, 100.25, 0.25)
@@ -86,8 +94,8 @@ try:
     plt.ylabel('Albedo')
 
     plt.xlabel('X')
-
-    plt.show()
+    
+    #plt.show()
 finally:
     np.savetxt("test01.txt", result, fmt = "%.5f")
     print result