Browse Source

Wrong argument for bessel function

Konstantin Ladutenko 10 years ago
parent
commit
daebcf25d4
2 changed files with 10 additions and 8 deletions
  1. 3 1
      compare.cc
  2. 7 7
      nmie.cc

+ 3 - 1
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/2.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;

+ 7 - 7
nmie.cc

@@ -1461,23 +1461,23 @@ namespace nmie {
       for (int i = size_param_.size() - 1; i >= 0 ; i--) {
         if (Rho <= size_param_[i]) {
           l = i;
+	  break;
         }
       }
-      ml = refractive_index_[l];
-    }
-
+      ml = refractive_index_[l]; // TODO check the correct layer is used
+    }    
     // Calculate spherical Bessel and Hankel functions and their derivatives
-    sbesjh(Rho*ml, jn, jnp, h1n, h1np);
-    printf("Rho*ml = %10.5er%+10.5ei\n",std::real(Rho*ml), std::imag(Rho*ml));
+    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));
     
     // 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("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);