Python version compiling too.

nmie.cc

@@ -229,7 +229,7 @@ int sbesjh(std::complex<double> z, int nmax, std::vector<std::complex<double> >&
 //   bj, by: Spherical Bessel functions                                             //
 //   bd: Logarithmic derivative                                                     //
 //**********************************************************************************//
-void sphericalBessel(std::complex<double> z, int nmax, std::vector<std::complex<double>>& bj, std::vector<std::complex<double>>& by, std::vector<std::complex<double>>& bd) {
+void sphericalBessel(std::complex<double> z, int nmax, std::vector<std::complex<double> >& bj, std::vector<std::complex<double> >& by, std::vector<std::complex<double> >& bd) {
 
     std::vector<std::complex<double> > jn, jnp, h1n, h1np;
     jn.resize(nmax);
@@ -250,7 +250,7 @@ void sphericalBessel(std::complex<double> z, int nmax, std::vector<std::complex<
 // external scattering field = incident + scattered
 // BH p.92 (4.37), 94 (4.45), 95 (4.50)
 // assume: medium is non-absorbing; refim = 0; Uabs = 0
-int fieldExt(int nmax, double Rho, double Phi, double Theta, std::vector<double> Pi, std::vector<double> Tau,
+void fieldExt(int nmax, double Rho, double Phi, double Theta, std::vector<double> Pi, std::vector<double> Tau,
              std::vector<std::complex<double> > an, std::vector<std::complex<double> > bn,
 		     std::vector<std::complex<double> >& E, std::vector<std::complex<double> >& H)  {
 

py_nmie.cc

@@ -60,20 +60,33 @@ int nMie(int L, int pl, std::vector<double> x, std::vector<std::complex<double>
 // std::vector<std::complex<double> >
 int nField(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m, int nmax,
            int nCoords, std::vector<double> Xp, std::vector<double> Yp, std::vector<double> Zp,
-           double Er[], double Ei[], double Hr[], double Hi[]) {
+           double Erx[], double Ery[], double Erz[], double Eix[], double Eiy[], double Eiz[],
+           double Hrx[], double Hry[], double Hrz[], double Hix[], double Hiy[], double Hiz[]) {
 
   int i, result;
-  std::vector<std::complex<double> > E, H;
+  std::vector<std::vector<std::complex<double> > > E, H;
   E.resize(nCoords);
   H.resize(nCoords);
+  for (i = 0; i < nCoords; i++) {
+    E[i].resize(3);
+    H[i].resize(3);
+  }
 
   result = nField(L, pl, x, m, nmax, nCoords, Xp, Yp, Zp, E, H);
 
   for (i = 0; i < nCoords; i++) {
-    Er[i] = E[i].real();
-    Ei[i] = E[i].imag();
-    Hr[i] = H[i].real();
-    Hi[i] = H[i].imag();
+    Erx[i] = E[i][0].real();
+    Ery[i] = E[i][1].real();
+    Erz[i] = E[i][2].real();
+    Eix[i] = E[i][0].imag();
+    Eiy[i] = E[i][1].imag();
+    Eiz[i] = E[i][2].imag();
+    Hrx[i] = H[i][0].real();
+    Hry[i] = H[i][1].real();
+    Hrz[i] = H[i][2].real();
+    Hix[i] = H[i][0].imag();
+    Hiy[i] = H[i][1].imag();
+    Hiz[i] = H[i][2].imag();
   }
 
   return result;

py_nmie.h

@@ -34,5 +34,6 @@ int nMie(int L, int pl, std::vector<double> x, std::vector<std::complex<double>
 
 int nField(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m, int nmax,
            int nCoords, std::vector<double> Xp, std::vector<double> Yp, std::vector<double> Zp,
-           double Er[], double Ei[], double Hr[], double Hi[]);
+           double Erx[], double Ery[], double Erz[], double Eix[], double Eiy[], double Eiz[],
+           double Hrx[], double Hry[], double Hrz[], double Hix[], double Hiy[], double Hiz[]);
 

scattnlay.pyx

@@ -51,7 +51,7 @@ cdef inline double *npy2c(np.ndarray a):
 
 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 Er[], double Ei[], double Hr[], double Hi[])
+    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 Py_ssize_t i
@@ -92,24 +92,44 @@ def fieldnlay(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] E = np.zeros((x.shape[0], coords.shape[0]), dtype = np.complex128)
-    cdef np.ndarray[np.complex128_t, ndim = 2] H = np.zeros((x.shape[0], coords.shape[0]), dtype = np.complex128)
-
-    cdef np.ndarray[np.float64_t, ndim = 1] Er
-    cdef np.ndarray[np.float64_t, ndim = 1] Ei
-    cdef np.ndarray[np.float64_t, ndim = 1] Hr
-    cdef np.ndarray[np.float64_t, ndim = 1] Hi
+    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]):
-        Er = np.zeros(coords.shape[0], dtype = np.float64)
-        Ei = np.zeros(coords.shape[0], dtype = np.float64)
-        Hr = np.zeros(coords.shape[0], dtype = np.float64)
-        Hi = 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(Er), npy2c(Ei), npy2c(Hr), npy2c(Hi))
-
-        E[i] = Er.copy('C') + 1.0j*Ei.copy('C')
-        H[i] = Hr.copy('C') + 1.0j*Hi.copy('C')
+        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][0] = Erx.copy('C') + 1.0j*Eix.copy('C')
+        E[i][1] = Ery.copy('C') + 1.0j*Eiy.copy('C')
+        E[i][2] = Erz.copy('C') + 1.0j*Eiz.copy('C')
+        H[i][0] = Hrx.copy('C') + 1.0j*Hix.copy('C')
+        H[i][1] = Hry.copy('C') + 1.0j*Hiy.copy('C')
+        H[i][2] = Hrz.copy('C') + 1.0j*Hiz.copy('C')
 
     return terms, E, H