Bugfix/fix asymptotic bessel

Makes the following corrections to multiprecision bessel function implementation:

• Include second series scaling with exp(-z) in asymptotic calculation of I
• Correct normalization of besinorm
• Make precision-dependent condition for selecting Taylor series or asymptotic approximation for I and K

The first change is the only one to affect the rest of the IW2D-scripts in their current form, as the normalized K and I are not used, and precision is hard coded. The first change is particularly important for small Re(z). All changes are based on the work of Amos: AMOS, D. E, Computation of Bessel functions of complex argument. SAND83-0086, Sandia National Laboratories, Albuquerque, N.M. (May 1983).

ImpedanceWake2D/multiprecision.cc

@@ -124,7 +124,7 @@ template<unsigned int Precision> void cmatrixgemm(int m, int n, int k,
 
     // use Taylor series and asymptotic expansions from Abramowitz-Stegun. 
     
-    amp::campf<Precision> power,sumi,sumk,zover2,z2,sumf,zover2m;
+    amp::campf<Precision> power,exponential,imagj,sumi,sumk,zover2,z2,sumf,zover2m;
     unsigned int condition,factm=1;
     amp::ampf<Precision> fact,abspower,k,expoabs,absz2,psi,m1powm,m1powk,mu,mMP;
     //timeval c1,c2; // clock ticks
@@ -140,7 +140,11 @@ template<unsigned int Precision> void cmatrixgemm(int m, int n, int k,
     z2=amp::csqr(zover2); // z^2/4
     absz2=amp::abscomplex(zover2); // |z/2|
     
-    if (absz2<30) {
+    // Decide whether to use taylor expansion or asymptotic formula for the calculation.
+    // Note: Asymptotic formula will not converge for abs(z) lower than what this condition accepts.
+    // Limit given in:
+    // Amos, Donald E. Computation of Bessel functions of complex argument. No. SAND-83-0086. Sandia National Lab.(SNL-NM), Albuquerque, NM (United States), 1983.
+    if ((absz2 < (2.4+0.36*Precision)/2) && (absz2 < mMP*mMP/4)) {
       // use Taylor series for both I and K
 
       expoabs=amp::exp(amp::abscomplex(z2)); // exp(|z^2/4|)
@@ -192,8 +196,9 @@ template<unsigned int Precision> void cmatrixgemm(int m, int n, int k,
       besk = sumf + m1powm*(-clog(zover2)*besi + zover2m*sumk/2);
       // normalized ones
       power=cexp(z);
-      besinorm = besi/power;
       besknorm = besk*power;
+      power=amp::exp(amp::abs(z.x));
+      besinorm = besi/power;
       
     } else {
       // use asymptotic formulae for large |z| (note: accuracy control is less good)
@@ -219,14 +224,16 @@ template<unsigned int Precision> void cmatrixgemm(int m, int n, int k,
       if (amp::atan2(z.y,z.x)>=amp::halfpi<Precision>()) std::cout << "Pb in bessel: argument of z=" <<
       		amp::atan2(z.y,z.x).toDec().c_str() << "\n";
       
-      // final results (normalized)
+      // final results
       power=1/csqrt(amp::twopi<Precision>()*z);
-      besinorm = sumi*power;
+      exponential=cexp(z);
+      imagj.x = 0; imagj.y = 1;
+
       besknorm = sumk*power*amp::pi<Precision>();
-      // unormalized ones
-      power=cexp(z);
-      besi = besinorm*power;
-      besk = besknorm/power;
+      besi = power*(exponential*sumi + imagj*cexp(imagj*amp::pi<Precision>()*mMP)/exponential*sumk);
+      besk = besknorm/exponential;
+      exponential = amp::exp( -amp::abs(z.x));
+      besinorm = besi*exponential;
       
     }