Merge branch 'master' of github.com:ovidiopr/scattnlay

Conflicts:
	nmie.cc

compare.cc

@@ -237,16 +237,18 @@ int main(int argc, char *argv[]) {
     double from_coord = -3.0, to_coord = 3000.0;
     //double size=2.0*PI*1.0/6.0;
     double size=0.001;
+    double R = size/(2.0*PI);
     std::vector<double> range;
     // for (int i = 0; i < samples; ++i) {
       //range.push_back( from_coord + (to_coord-from_coord)/(static_cast<double>(samples)-1)*i );
     //range.push_back(size*0.01);
     //range.push_back(size*0.99999);
-    range.push_back(size/20.0);
+    range.push_back(R/2.0);
     //range.push_back(size*1.00001);
     //range.push_back(3);
     //printf("r=%g  ", range.back());
     //}
+    printf("r/2 = %g\n", R/2.0);
     int samples = range.size();
     std::vector<double> zero(samples, 0.0);
     std::vector<double> Xp, Yp, Zp;

nmie.cc

@@ -800,7 +800,13 @@ namespace nmie {
 
     // using eq 4.50 in BH
     std::complex<double> c_zero(0.0, 0.0);
+<<<<<<< HEAD
 
+=======
+    std::complex<double> deriv = Rho*dzn + zn;
+    //std::complex<double> deriv (6.6667e-04, 0.0);
+    //printf("n=%g, Phi=%g Pi=%g Tau=%g zn=%+er%+ei deriv=%+er%+ei\n", n, Phi, Pi, Tau, zn.real(), zn.imag(), deriv.real(), deriv.imag());
+>>>>>>> 645549af67fbee6f0583798534e030f27a500a44
     using std::sin;
     using std::cos;
     Mo1n[0] = c_zero;
@@ -1457,8 +1463,10 @@ namespace nmie {
       for (int i = size_param_.size() - 1; i >= 0 ; i--) {
         if (Rho <= size_param_[i]) {
           l = i;
+	  break;
         }
       }
+<<<<<<< HEAD
       ml = refractive_index_[l];
     }
 
@@ -1469,21 +1477,42 @@ namespace nmie {
     // Calculate spherical Bessel and Hankel functions and their derivatives
     //sbesjh(Rho*ml, jn, jnp, h1n, h1np);
 
+=======
+      ml = refractive_index_[l]; // TODO check the correct layer is used
+    }    
+    // Calculate spherical Bessel and Hankel functions and their derivatives
+    sbesjh(2.0*PI_*Rho*ml, jn, jnp, h1n, h1np);
+    printf("2.0*PI*Rho*ml = %10.5er%+10.5ei\n",std::real(2.0*PI_*Rho*ml), std::imag(2.0*PI_*Rho*ml));
+    
+>>>>>>> 645549af67fbee6f0583798534e030f27a500a44
     // Calculate angular functions Pi and Tau
     calcPiTau(std::cos(Theta), Pi, Tau);
-/*    printf("Thetd = %g, cos(Theta) = %g\n", Theta, std::cos(Theta));
+    printf("Thetd = %g, cos(Theta) = %g\n", Theta, std::cos(Theta));
+    printf("jn:\n");
+    for (auto p : jn) printf("%+11.4er%+11.4ei\n",p.real(), p.imag());
     printf("Pi:\n");
     for (auto p : Pi) printf("%11.4e\n",p);
     printf("Tau:\n");
     for (auto p : Tau) printf("%11.4e\n",p);
-*/
+
     for (int n = nmax_ - 2; n >= 0; n--) {
       int n1 = n + 1;
       double rn = static_cast<double>(n1);
 
       // using BH 4.12 and 4.50
+<<<<<<< HEAD
       calcSpherHarm(Rho, Theta, Phi, Psi[n1], D1n[n1], Pi[n], Tau[n], rn, M1o1n, M1e1n, N1o1n, N1e1n);
       calcSpherHarm(Rho, Theta, Phi, Zeta[n1], D3n[n1], Pi[n], Tau[n], rn, M3o1n, M3e1n, N3o1n, N3e1n);
+=======
+      calcSpherHarm(Rho, Theta, Phi, jn[n1], jnp[n1], Pi[n], Tau[n], rn, M1o1n, M1e1n, N1o1n, N1e1n);
+      calcSpherHarm(Rho, Theta, Phi, h1n[n1], h1np[n1], Pi[n], Tau[n], rn, M3o1n, M3e1n, N3o1n, N3e1n); //TODO uncomment line
+      auto deriv1 = -rn*jn[n1]+Rho*jn[n1-1];
+      auto deriv2 = Rho*jnp[n1] + jn[n1];
+      printf("n=%d   deriv1: %+11.4e   deriv2: %+11.4ei\n",n1, deriv1.real(), deriv2.real());
+      // printf("N1e1n[%d]: ", n1);
+      // for (auto p : N1e1n) printf("%+11.4er%+11.4ei\t",p.real(), p.imag());
+      // printf("\n");
+>>>>>>> 645549af67fbee6f0583798534e030f27a500a44
 
       // Total field in the lth layer: eqs. (1) and (2) in Yang, Appl. Opt., 42 (2003) 1710-1720
       std::complex<double> En = ipow[n1 % 4]*(rn + rn + 1.0)/(rn*rn + rn);