fixed passing by pointer

Makefile

@@ -28,10 +28,10 @@ source:
 cython: scattnlay.pyx
     # create c++ code for double precision module
 	$(CYTHON) --cplus scattnlay.pyx -o $(SRCDIR)/scattnlay.cpp
-	# create c++ code for MP module
-	ln -s scattnlay.pyx scattnlay_mp.pyx
-	$(CYTHON) --cplus scattnlay_mp.pyx -o $(SRCDIR)/scattnlay_mp.cpp
-	rm scattnlay_mp.pyx
+	# # create c++ code for MP module
+	# ln -s scattnlay.pyx scattnlay_mp.pyx
+	# $(CYTHON) --cplus scattnlay_mp.pyx -o $(SRCDIR)/scattnlay_mp.cpp
+	# rm scattnlay_mp.pyx
 
 python_ext: $(SRCDIR)/nmie.cc $(SRCDIR)/py_nmie.cc $(SRCDIR)/scattnlay.cpp $(SRCDIR)/scattnlay_mp.cpp
 	$(PYTHON) setup.py build_ext --inplace

examples/calc-expansion-coeffs-spectra.py

@@ -0,0 +1,77 @@
+#!/usr/bin/env python
+# -*- 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/>.
+
+from scattnlay import fieldnlay, scattnlay, expansioncoeffs, scattcoeffs
+
+import numpy as np
+import cmath
+
+WL=550 #nm
+core_r = 102
+#core_r = 120
+
+index_NP = 4.0
+x = np.ones((1), dtype = np.float64)
+m = np.ones((1), dtype = np.complex128)
+
+import matplotlib.pyplot as plt
+
+# dipole - 1 -- 2
+# quad   - 2 -- 4
+# octo   - 3 -- 8
+# hex    - 4 -- 16
+# 32     - 5 -- 32
+npts = 151
+ext = ".png"
+# npts = 351
+# ext = ".pdf"
+x[0] = 2.0*np.pi*core_r/WL#/4.0*3.0
+m[0] = index_NP
+#    for mode_type, field_to_plot, WL, mode_n,  crossplane, isStream in plot_params :
+comment='bulk-NP-R'+str(core_r)+'nm-WL'+str(WL)
+# scattnlay function process several NP designs in one call
+x = np.array([x])
+m = np.array([m])
+#print(x.shape)
+terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(x, m)
+nmax = terms
+#print(nmax)
+terms, an, bn = scattcoeffs(x, m, nmax)
+#print(terms, an)
+terms, aln, bln, cln, dln = expansioncoeffs(x, m, nmax)
+print(terms, dln)
+# fig, axs = plt.subplots(1,1)#, sharey=True, sharex=True)
+# fig.tight_layout()
+# fig.subplots_adjust(hspace=0.3, wspace=-0.1)
+
+# # plt.savefig("Egor3/"+"%02d"%(i)+
+# #             comment+"-R"+str(int(round(x[-1]*WL/2.0/np.pi)))+"-"+crossplane+"-"
+# #                 +field_to_plot+"-mode"+mode+mt+st+ext,pad_inches=0.02, bbox_inches='tight')
+# plt.clf()
+# plt.close()
+
+#print("end")

examples/fieldplot.py

@@ -286,8 +286,8 @@ def fieldplot(fig, ax, x, m, WL, comment='', WL_units=' ', crossplane='XZ',
         #     cbar.ax.set_yticklabels(['%3.0f' % (a) for a in scale_ticks])
         # else:
         #     cbar.ax.set_yticklabels(['%g' % (a) for a in scale_ticks])
-        # pos = list(cbar.ax.get_position().bounds)
-        #fig.text(pos[0] - 0.02, 0.925, '|E|/|E$_0$|', fontsize = 14)
+        # # pos = list(cbar.ax.get_position().bounds)
+        # #fig.text(pos[0] - 0.02, 0.925, '|E|/|E$_0$|', fontsize = 14)
         lp2 = -10.0
         lp1 = -1.0
         if crossplane == 'XZ':

scattnlay.pyx

@@ -43,17 +43,17 @@ cdef inline double *npy2c(np.ndarray a):
 
 cdef extern from "py_nmie.h":
     cdef int ScattCoeffs(int L, int pl, vector[double] x, vector[complex] m, int nmax, double anr[], double ani[], double bnr[], double bni[])
-    cdef int ExpansionCoeffs( int L,  int pl, vector[double] x,  vector[complex] m,
-                     int nmax,
-                    vector[vector[double] ]  alnr,
-                    vector[vector[double] ]  alni,
-                    vector[vector[double] ]  blnr,
-                    vector[vector[double] ]  blni,
-                    vector[vector[double] ]  clnr,
-                    vector[vector[double] ]  clni,
-                    vector[vector[double] ]  dlnr,
-                    vector[vector[double] ]  dlni)
-    # cdef int ExpansionCoeffs(int L, int pl, vector[double] x, vector[complex] m, int nmax, double alnr[], double alni[], double blnr[], double blni[], double clnr[], double clni[], double dlnr[], double dlni[])
+    # cdef int ExpansionCoeffs( int L,  int pl, vector[double] x,  vector[complex] m,
+    #                  int nmax,
+    #                 vector[vector[double] ]&  alnr,
+    #                 vector[vector[double] ]&  alni,
+    #                 vector[vector[double] ]&  blnr,
+    #                 vector[vector[double] ]&  blni,
+    #                 vector[vector[double] ]&  clnr,
+    #                 vector[vector[double] ]&  clni,
+    #                 vector[vector[double] ]&  dlnr,
+    #                 vector[vector[double] ]&  dlni)
+    cdef int ExpansionCoeffs(int L, int pl, vector[double] x, vector[complex] m, int nmax, double alnr[], double alni[], double blnr[], double blni[], double clnr[], double clni[], double dlnr[], double dlni[])
     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 mode_n, int mode_type, 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[])
 
@@ -62,8 +62,10 @@ def scattcoeffs(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t
 
     cdef np.ndarray[np.int_t, ndim = 1] terms = np.zeros(x.shape[0], dtype = np.int)
 
-    cdef np.ndarray[np.complex128_t, ndim = 2] an = np.zeros((x.shape[0], nmax), dtype = np.complex128)
-    cdef np.ndarray[np.complex128_t, ndim = 2] bn = np.zeros((x.shape[0], nmax), dtype = np.complex128)
+    cdef np.ndarray[np.complex128_t, ndim = 2] an = np.zeros(
+        (x.shape[0], nmax), dtype = np.complex128)
+    cdef np.ndarray[np.complex128_t, ndim = 2] bn = np.zeros(
+        (x.shape[0], nmax), dtype = np.complex128)
 
     cdef np.ndarray[np.float64_t, ndim = 1] anr
     cdef np.ndarray[np.float64_t, ndim = 1] ani
@@ -88,38 +90,37 @@ def expansioncoeffs(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex1
     cdef Py_ssize_t l
 
     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] aln = np.zeros((x.shape[0], x.shape[1]+1, nmax), dtype = np.complex128)
     cdef np.ndarray[np.complex128_t, ndim = 3] bln = np.zeros((x.shape[0], x.shape[1]+1, nmax), dtype = np.complex128)
     cdef np.ndarray[np.complex128_t, ndim = 3] cln = np.zeros((x.shape[0], x.shape[1]+1, nmax), dtype = np.complex128)
     cdef np.ndarray[np.complex128_t, ndim = 3] dln = np.zeros((x.shape[0], x.shape[1]+1, nmax), dtype = np.complex128)
 
-    cdef np.ndarray[np.float64_t, ndim = 2] alnr
-    cdef np.ndarray[np.float64_t, ndim = 2] alni
-    cdef np.ndarray[np.float64_t, ndim = 2] blnr
-    cdef np.ndarray[np.float64_t, ndim = 2] blni
-    cdef np.ndarray[np.float64_t, ndim = 2] clnr
-    cdef np.ndarray[np.float64_t, ndim = 2] clni
-    cdef np.ndarray[np.float64_t, ndim = 2] dlnr
-    cdef np.ndarray[np.float64_t, ndim = 2] dlni
+    cdef np.ndarray[np.float64_t, ndim = 1] alnr
+    cdef np.ndarray[np.float64_t, ndim = 1] alni
+    cdef np.ndarray[np.float64_t, ndim = 1] blnr
+    cdef np.ndarray[np.float64_t, ndim = 1] blni
+    cdef np.ndarray[np.float64_t, ndim = 1] clnr
+    cdef np.ndarray[np.float64_t, ndim = 1] clni
+    cdef np.ndarray[np.float64_t, ndim = 1] dlnr
+    cdef np.ndarray[np.float64_t, ndim = 1] dlni
 
     for i in range(x.shape[0]):
-        alnr = np.zeros((x.shape[1]+1,nmax), dtype = np.float64)
-        alni = np.zeros((x.shape[1]+1,nmax), dtype = np.float64)
-        blnr = np.zeros((x.shape[1]+1,nmax), dtype = np.float64)
-        blni = np.zeros((x.shape[1]+1,nmax), dtype = np.float64)
-        clnr = np.zeros((x.shape[1]+1,nmax), dtype = np.float64)
-        clni = np.zeros((x.shape[1]+1,nmax), dtype = np.float64)
-        dlnr = np.zeros((x.shape[1]+1,nmax), dtype = np.float64)
-        dlni = np.zeros((x.shape[1]+1,nmax), dtype = np.float64)
-
-        terms[i] = ExpansionCoeffs(x.shape[1], pl, x[i].copy('C'), m[i].copy('C'), nmax, alnr, alni, blnr, blni, clnr, clni, dlnr, dlni)
-        # terms[i] = ExpansionCoeffs(x.shape[1], pl, x[i].copy('C'), m[i].copy('C'), nmax, npy2c(alnr), npy2c(alni), npy2c(blnr), npy2c(blni), npy2c(clnr), npy2c(clni), npy2c(dlnr), npy2c(dlni))
-
-        for l in range(x.shape[1]):
-            aln[l][i] = alnr[l].copy('C') + 1.0j*alni[l].copy('C')
-            bln[l][i] = blnr[l].copy('C') + 1.0j*blni[l].copy('C')
-
+        alnr = np.zeros((x.shape[1]+1)*nmax, dtype = np.float64)
+        alni = np.zeros((x.shape[1]+1)*nmax, dtype = np.float64)
+        blnr = np.zeros((x.shape[1]+1)*nmax, dtype = np.float64)
+        blni = np.zeros((x.shape[1]+1)*nmax, dtype = np.float64)
+        clnr = np.zeros((x.shape[1]+1)*nmax, dtype = np.float64)
+        clni = np.zeros((x.shape[1]+1)*nmax, dtype = np.float64)
+        dlnr = np.zeros((x.shape[1]+1)*nmax, dtype = np.float64)
+        dlni = np.zeros((x.shape[1]+1)*nmax, dtype = np.float64)
+
+        #terms[i] = ExpansionCoeffs(x.shape[1], pl, x[i].copy('C'), m[i].copy('C'), nmax, alnr, alni, blnr, blni, clnr, clni, dlnr, dlni)
+        terms[i] = ExpansionCoeffs(x.shape[1], pl, x[i].copy('C'), m[i].copy('C'), nmax, npy2c(alnr), npy2c(alni), npy2c(blnr), npy2c(blni), npy2c(clnr), npy2c(clni), npy2c(dlnr), npy2c(dlni))
+        for l in range(x.shape[1]+1):
+            aln[i][l] = alnr[l*nmax:(l+1)*nmax].copy('C') + 1.0j*alni[l*nmax:(l+1)*nmax].copy('C')
+            bln[i][l] = blnr[l*nmax:(l+1)*nmax].copy('C') + 1.0j*blni[l*nmax:(l+1)*nmax].copy('C')
+            cln[i][l] = clnr[l*nmax:(l+1)*nmax].copy('C') + 1.0j*clni[l*nmax:(l+1)*nmax].copy('C')
+            dln[i][l] = dlnr[l*nmax:(l+1)*nmax].copy('C') + 1.0j*dlni[l*nmax:(l+1)*nmax].copy('C')
     return terms, aln, bln, cln, dln
 
 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 nmax = -1, np.int_t pl = -1):

setup.py

@@ -57,11 +57,12 @@ O. Pena, U. Pal, Comput. Phys. Commun. 180 (2009) 2348-2354.""",
                                ["src/nmie.cc", "src/py_nmie.cc", "src/scattnlay.cpp"],
                                language = "c++",
                                include_dirs = [np.get_include()], 
-                               extra_compile_args=['-std=c++11']),
-                     Extension("scattnlay_mp",
-                               ["src/nmie.cc", "src/py_nmie.cc", "src/scattnlay_mp.cpp"],
-                               language = "c++",
-                               include_dirs = [np.get_include()], 
-                               extra_compile_args=['-std=c++11', '-DMULTI_PRECISION=100'])]
+                               extra_compile_args=['-std=c++11'])
+                     # ,Extension("scattnlay_mp",
+                     #           ["src/nmie.cc", "src/py_nmie.cc", "src/scattnlay_mp.cpp"],
+                     #           language = "c++",
+                     #           include_dirs = [np.get_include()], 
+                     #           extra_compile_args=['-std=c++11', '-DMULTI_PRECISION=100'])
+      ]
 )
 

src/py_nmie.cc

@@ -59,17 +59,22 @@ int ScattCoeffs(const unsigned int L, const int pl,
 // return the results (useful for python).  This is a workaround
 // because I have not been able to return the results using
 // std::vector<std::vector<std::complex<double> > >
+// int ExpansionCoeffs(const unsigned int L, const int pl,
+//                     std::vector<double>& x, std::vector<std::complex<double> >& m,
+//                     const int nmax,
+//                     std::vector<std::vector<double> >&  alnr,
+//                     std::vector<std::vector<double> >&  alni,
+//                     std::vector<std::vector<double> >&  blnr,
+//                     std::vector<std::vector<double> >&  blni,
+//                     std::vector<std::vector<double> >&  clnr,
+//                     std::vector<std::vector<double> >&  clni,
+//                     std::vector<std::vector<double> >&  dlnr,
+//                     std::vector<std::vector<double> >&  dlni)
 int ExpansionCoeffs(const unsigned int L, const int pl,
                     std::vector<double>& x, std::vector<std::complex<double> >& m,
                     const int nmax,
-                    std::vector<std::vector<double> >&  alnr,
-                    std::vector<std::vector<double> >&  alni,
-                    std::vector<std::vector<double> >&  blnr,
-                    std::vector<std::vector<double> >&  blni,
-                    std::vector<std::vector<double> >&  clnr,
-                    std::vector<std::vector<double> >&  clni,
-                    std::vector<std::vector<double> >&  dlnr,
-                    std::vector<std::vector<double> >&  dlni) {
+                    double alnr[], double alni[], double blnr[], double blni[],
+                    double clnr[], double clni[], double dlnr[], double dlni[]) {
 
   int result;
   std::vector< std::vector<std::complex<double> > > aln, bln, cln, dln;
@@ -86,16 +91,18 @@ int ExpansionCoeffs(const unsigned int L, const int pl,
 
   result = nmie::ExpansionCoeffs(L, pl, x, m, nmax, aln, bln, cln, dln);
 
+  
   for (unsigned int l = 0; l <= L; l++) {
     for (int i = 0; i < result; i++) {
-      alnr[l][i] = aln[l][i].real();
-      alni[l][i] = aln[l][i].imag();
-      blnr[l][i] = bln[l][i].real();
-      blni[l][i] = bln[l][i].imag();
-      clnr[l][i] = cln[l][i].real();
-      clni[l][i] = cln[l][i].imag();
-      dlnr[l][i] = dln[l][i].real();
-      dlni[l][i] = dln[l][i].imag();
+      alnr[l*(result)+i] = aln[l][i].real();
+      alni[l*(result)+i] = aln[l][i].imag();
+      blnr[l*(result)+i] = bln[l][i].real();
+      blni[l*(result)+i] = bln[l][i].imag();
+      clnr[l*(result)+i] = cln[l][i].real();
+      clni[l*(result)+i] = cln[l][i].imag();
+      dlnr[l*(result)+i] = dln[l][i].real();
+      //std::cout<<l*(result)+i<<": "<<dlnr[l*(result)+i]<<std::endl;
+      dlni[l*(result)+i] = dln[l][i].imag();
     }
   }
   return result;

src/py_nmie.h

@@ -33,22 +33,23 @@ int ScattCoeffs(const unsigned int L, const int pl,
                 const int nmax,
                 double anr[], double ani[], double bnr[], double bni[]);
 
-int ExpansionCoeffs(const unsigned int L, const int pl,
-                    std::vector<double>& x, std::vector<std::complex<double> >& m,
-                    const int nmax,
-                    std::vector<std::vector<double> >&  alnr,
-                    std::vector<std::vector<double> >&  alni,
-                    std::vector<std::vector<double> >&  blnr,
-                    std::vector<std::vector<double> >&  blni,
-                    std::vector<std::vector<double> >&  clnr,
-                    std::vector<std::vector<double> >&  clni,
-                    std::vector<std::vector<double> >&  dlnr,
-                    std::vector<std::vector<double> >&  dlni);
+//std::vector<std::vector<std::vector<double> > >
 /* int ExpansionCoeffs(const unsigned int L, const int pl, */
 /*                     std::vector<double>& x, std::vector<std::complex<double> >& m, */
 /*                     const int nmax, */
-/*                     double alnr[], double alni[], double blnr[], double blni[], */
-/*                     double clnr[], double clni[], double dlnr[], double dlni[]); */
+/*                     std::vector<std::vector<double> >&  alnr, */
+/*                     std::vector<std::vector<double> >&  alni, */
+/*                     std::vector<std::vector<double> >&  blnr, */
+/*                     std::vector<std::vector<double> >&  blni, */
+/*                     std::vector<std::vector<double> >&  clnr, */
+/*                     std::vector<std::vector<double> >&  clni, */
+/*                     std::vector<std::vector<double> >&  dlnr, */
+/*                     std::vector<std::vector<double> >&  dlni); */
+int ExpansionCoeffs(const unsigned int L, const int pl,
+                    std::vector<double>& x, std::vector<std::complex<double> >& m,
+                    const int nmax,
+                    double alnr[], double alni[], double blnr[], double blni[],
+                    double clnr[], double clni[], double dlnr[], double dlni[]);
 
 int nMie(const int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m,
          const int nTheta, std::vector<double>& Theta, const int nmax,