Cleaned code. Modified script to calculate flow lines; it still needs more work.

nmie.cc

@@ -1102,7 +1102,7 @@ namespace nmie {
       z = size_param_[l]*m[l];
       z1 = size_param_[l]*m1[l];
 
-      printf("\nz[%i] = %gr%+gi; z1[%i] = %gr%+gi\n", l, std::real(z), std::imag(z), l, std::real(z1), std::imag(z1));
+      //printf("\nz[%i] = %gr%+gi; z1[%i] = %gr%+gi\n", l, std::real(z), std::imag(z), l, std::real(z1), std::imag(z1));
 
       calcD1D3(z, D1z, D3z);
       calcD1D3(z1, D1z1, D3z1);
@@ -1112,7 +1112,7 @@ namespace nmie {
       for (int n = 0; n < nmax_; n++) {
         int n1 = n + 1;
 
-        printf("Psiz[%02i] = %11.4e,%11.4e\tZetaz[%02i] = %11.4e,%11.4e\tPsiz1[%02i] = %11.4e,%11.4e\tZetaz1[%02i] = %11.4e,%11.4e\n", n1, real(Psiz[n1]), imag(Psiz[n1]), n1, real(Zetaz[n1]), imag(Zetaz[n1]), n1, real(Psiz1[n1]), imag(Psiz1[n1]), n1, real(Zetaz1[n1]), imag(Zetaz1[n1]));
+        //printf("Psiz[%02i] = %11.4e,%11.4e\tZetaz[%02i] = %11.4e,%11.4e\tPsiz1[%02i] = %11.4e,%11.4e\tZetaz1[%02i] = %11.4e,%11.4e\n", n1, real(Psiz[n1]), imag(Psiz[n1]), n1, real(Zetaz[n1]), imag(Zetaz[n1]), n1, real(Psiz1[n1]), imag(Psiz1[n1]), n1, real(Zetaz1[n1]), imag(Zetaz1[n1]));
 
         denomZeta = Zetaz[n1]*(D1z[n1] - D3z[n1]);
         denomPsi  =  Psiz[n1]*(D1z[n1] - D3z[n1]);
@@ -1123,7 +1123,7 @@ namespace nmie {
         T3 = (dln_[l + 1][n]*D1z1[n1]*Psiz1[n1] - aln_[l + 1][n]*D3z1[n1]*Zetaz1[n1])*m[l]/m1[l];
         T4 =  cln_[l + 1][n]*D1z1[n1]*Psiz1[n1] - bln_[l + 1][n]*D3z1[n1]*Zetaz1[n1];
 
-        printf("T1[%02i] = %11.4e,%11.4e\tT2[%02i] = %11.4e,%11.4e\tT3[%02i] = %11.4e,%11.4e\tT4[%02i] = %11.4e,%11.4e\n", n, real(D1z[n1]*T1 + T3), imag(D1z[n1]*T1 + T3), n, real(D1z[n1]*T2 + T4), imag(D1z[n1]*T2 + T4), n, real(D3z[n1]*T2 + T4), imag(D3z[n1]*T2 + T4), n, real(D3z[n1]*T1 + T3), imag(D3z[n1]*T1 + T3));
+        //printf("T1[%02i] = %11.4e,%11.4e\tT2[%02i] = %11.4e,%11.4e\tT3[%02i] = %11.4e,%11.4e\tT4[%02i] = %11.4e,%11.4e\n", n, real(D1z[n1]*T1 + T3), imag(D1z[n1]*T1 + T3), n, real(D1z[n1]*T2 + T4), imag(D1z[n1]*T2 + T4), n, real(D3z[n1]*T2 + T4), imag(D3z[n1]*T2 + T4), n, real(D3z[n1]*T1 + T3), imag(D3z[n1]*T1 + T3));
 
         // aln
         aln_[l][n] = (D1z[n1]*T1 + T3)/denomZeta;
@@ -1138,22 +1138,20 @@ namespace nmie {
 
     // Check the result and change  aln_[0][n] and aln_[0][n] for exact zero
     for (int n = 0; n < nmax_; ++n) {
-      aln_[0][n] = 0.0;
-      bln_[0][n] = 0.0;
       //printf("n=%d, aln_=%g,%g,   bln_=%g,%g \n", n, real(aln_[0][n]), imag(aln_[0][n]),
       //   real(bln_[0][n]), imag(bln_[0][n]));
-      //if (std::abs(aln_[0][n]) < 1e-10) aln_[0][n] = 0.0;
-      //else throw std::invalid_argument("Unstable calculation of aln_[0][n]!");
-      //if (std::abs(bln_[0][n]) < 1e-10) bln_[0][n] = 0.0;
-      //else throw std::invalid_argument("Unstable calculation of bln_[0][n]!");
+      if (std::abs(aln_[0][n]) < 1e-10) aln_[0][n] = 0.0;
+      else throw std::invalid_argument("Unstable calculation of aln_[0][n]!");
+      if (std::abs(bln_[0][n]) < 1e-10) bln_[0][n] = 0.0;
+      else throw std::invalid_argument("Unstable calculation of bln_[0][n]!");
     }
 
-    for (int l = L; l >= 0; l--) {
-      for (int n = 0; n < nmax_; n++) {
-        printf("aln_[%02i, %02i] = %11.4e,%11.4e\tbln_[%02i, %02i] = %11.4e,%11.4e\tcln_[%02i, %02i] = %11.4e,%11.4e\tdln_[%02i, %02i] = %11.4e,%11.4e\n", l, n, real(aln_[l][n]), imag(aln_[l][n]), l, n, real(bln_[l][n]), imag(bln_[l][n]), l, n, real(cln_[l][n]), imag(cln_[l][n]), l, n, real(dln_[l][n]), imag(dln_[l][n]));
-      }
-      printf("\n");
-    }
+//    for (int l = L; l >= 0; l--) {
+//      for (int n = 0; n < nmax_; n++) {
+//        printf("aln_[%02i, %02i] = %11.4e,%11.4e\tbln_[%02i, %02i] = %11.4e,%11.4e\tcln_[%02i, %02i] = %11.4e,%11.4e\tdln_[%02i, %02i] = %11.4e,%11.4e\n", l, n, real(aln_[l][n]), imag(aln_[l][n]), l, n, real(bln_[l][n]), imag(bln_[l][n]), l, n, real(cln_[l][n]), imag(cln_[l][n]), l, n, real(dln_[l][n]), imag(dln_[l][n]));
+//      }
+//      printf("\n");
+//    }
 
     isExpCoeffsCalc_ = true;
   }  // end of   void MultiLayerMie::ExpanCoeffs()

tests/python/field-Ag-flow.py

@@ -62,6 +62,7 @@ def GetFlow(scale_x, scale_z, Ec, Hc, a, b, nmax):
         z_pos = flow_z[-1]
         x_idx = get_index(scale_x, x_pos)
         z_idx = get_index(scale_z, z_pos)
+        #print x_idx, z_idx
         S=np.cross(Ec[npts*z_idx+x_idx], np.conjugate(Hc[npts*z_idx+x_idx]) ).real
         #if (np.linalg.norm(S)> 1e-4):
         Snorm=S/np.linalg.norm(S)
@@ -69,22 +70,26 @@ def GetFlow(scale_x, scale_z, Ec, Hc, a, b, nmax):
         dpos = abs(scale_z[0]-scale_z[1])/2.0
         dx = dpos*Snorm[0]
         dz = dpos*Snorm[2]
-        x_pos = x_pos+dx
-        z_pos = z_pos+dz
+        if (dz < 0.0):
+            x_pos = x_pos-dx
+            z_pos = z_pos-dz
+        else:
+            x_pos = x_pos+dx
+            z_pos = z_pos+dz
         #3. Save result
         flow_x.append(x_pos)
         flow_z.append(z_pos)
     return flow_x, flow_z
 
 # # a)
-# WL=400 #nm
-# core_r = WL/20.0
-# epsilon_Ag = -2.0 + 10.0j
+#WL=400 #nm
+#core_r = WL/20.0
+#epsilon_Ag = -2.0 + 10.0j
 
 # # b)
-# WL=400 #nm
-# core_r = WL/20.0
-# epsilon_Ag = -2.0 + 1.0j
+#WL=400 #nm
+#core_r = WL/20.0
+#epsilon_Ag = -2.0 + 1.0j
 
 # c)
 WL=354 #nm
@@ -115,7 +120,7 @@ print "m =", m
 
 npts = 281
 
-factor=5
+factor=2
 scan = np.linspace(-factor*x[0, 0], factor*x[0, 0], npts)
 
 coordX, coordZ = np.meshgrid(scan, scan)
@@ -128,23 +133,11 @@ coord = np.vstack((coordX, coordY, coordZ)).transpose()
 
 terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(x, m)
 terms, E, H = fieldnlay(x, m, coord)
-Er = np.absolute(E)
-Hr = np.absolute(H)
-P=[]
-for n in range(0, len(E[0])):
-    P.append(np.linalg.norm( np.cross(E[0][n], np.conjugate(H[0][n]) ).real/2 ))
 
-print(min(P))
+P = np.array(map(lambda n: np.linalg.norm( np.cross(E[0][n], np.conjugate(H[0][n]) ).real ), range(0, len(E[0]))))
 
-# |E|/|Eo|
-Eabs = np.sqrt(Er[0, :, 0]**2 + Er[0, :, 1]**2 + Er[0, :, 2]**2)
 Ec = E[0, :, :]
 Hc = H[0, :, :]
-Eangle = np.angle(E[0, :, 0])/np.pi*180
-
-Habs= np.sqrt(Hr[0, :, 0]**2 + Hr[0, :, 1]**2 + Hr[0, :, 2]**2)
-Hangle = np.angle(H[0, :, 1])/np.pi*180
-
 
 
 try:
@@ -201,7 +194,7 @@ try:
 
     from matplotlib.path import Path
     #import matplotlib.patches as patches
-    flow_total = 31
+    flow_total = 41
     for flow in range(0,flow_total):
         flow_x, flow_z = GetFlow(scale_x, scale_z, Ec, Hc,
                                  min(scale_x)+flow*(scale_x[-1]-scale_x[0])/(flow_total-1),

tests/python/field.py

@@ -57,7 +57,7 @@ m[0, 1] = n2/nm
 print "x =", x
 print "m =", m
 
-npts = 281
+npts = 501
 
 scan = np.linspace(-4.0*x[0, 0], 4.0*x[0, 0], npts)