//// jv.cpp//// $Id: jv.cpp,v 1.4 2006/01/22 13:40:37 ediap Exp $// #include <cmath>#include <itpp/base/itassert.h>#include <itpp/base/scalfunc.h>#include <itpp/base/bessel/bessel_internal.h>using namespace itpp;// This is slightly modified routine from the Cephes library, see http://www.netlib.org/cephes///  // According to licence agreement this software can be used freely.///* *	Bessel function of noninteger order * * double v, x, y, jv(); * * y = jv( v, x ); * * DESCRIPTION: * * Returns Bessel function of order v of the argument, * where v is real.  Negative x is allowed if v is an integer. * * Several expansions are included: the ascending power * series, the Hankel expansion, and two transitional * expansions for large v.  If v is not too large, it * is reduced by recurrence to a region of best accuracy. * The transitional expansions give 12D accuracy for v > 500. * * ACCURACY: * Results for integer v are indicated by *, where x and v * both vary from -125 to +125.  Otherwise, * x ranges from 0 to 125, v ranges as indicated by "domain." * Error criterion is absolute, except relative when |jv()| > 1. * * arithmetic  v domain  x domain    # trials      peak       rms *    IEEE      0,125     0,125      100000      4.6e-15    2.2e-16 *    IEEE   -125,0       0,125       40000      5.4e-11    3.7e-13 *    IEEE      0,500     0,500       20000      4.4e-15    4.0e-16 * Integer v: *    IEEE   -125,125   -125,125      50000      3.5e-15*   1.9e-16* *//*Cephes Math Library Release 2.8:  June, 2000Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier*/#ifdef _MSC_VER#pragma warning(disable:4996)#endif#define MAXGAM 171.624376956302725static double recur(double *, double, double *, int);static double jvs(double, double);static double hankel(double, double);static double jnx(double, double);static double jnt(double, double);#define MAXNUM 1.79769313486231570815E308    /* 2**1024*(1-MACHEP) */#define MACHEP 1.11022302462515654042E-16   /* 2**-53 */#define MAXLOG 7.08396418532264106224E2     /* log 2**1022 */#define MINLOG -7.08396418532264106224E2    /* log 2**-1022 */#define PI 3.14159265358979323846       /* pi */#define BIG  1.44115188075855872E+17// ---------------------------- jv() -------------------------------------------------------double jv(double n, double x){  double k, q, t, y, an;  int i, sign, nint;  nint = 0;	/* Flag for integer n */  sign = 1;	/* Flag for sign inversion */  an = fabs( n );  y = floor( an );  if( y == an )    {      nint = 1;      i = int(an - 16384.0 * floor( an/16384.0 ));      if( n < 0.0 )	{	  if( i & 1 )	    sign = -sign;	  n = an;	}      if( x < 0.0 )	{	  if( i & 1 )	    sign = -sign;	  x = -x;	}      if( n == 0.0 ) // use 0th order bessel function	return( j0(x) );      if( n == 1.0 ) // use 1th order bessel function	return( sign * j1(x) );    }  it_error_if( (x < 0.0) && (y != an), "besselj:: negative values only allowed for integer orders.");  y = fabs(x);  if( y < MACHEP )    goto underf;  k = 3.6 * sqrt(y);  t = 3.6 * sqrt(an);  if( (y < t) && (an > 21.0) )    return( sign * jvs(n,x) );  if( (an < k) && (y > 21.0) )    return( sign * hankel(n,x) );  if( an < 500.0 )    {      /* Note: if x is too large, the continued       * fraction will fail; but then the       * Hankel expansion can be used.       */      if( nint != 0 )	{	  k = 0.0;	  q = recur( &n, x, &k, 1 );	  if( k == 0.0 )	    {	      y = j0(x)/q;	      goto done;	    }	  if( k == 1.0 )	    {	      y = j1(x)/q;	      goto done;	    }	}      if( an > 2.0 * y )	goto rlarger;      if( (n >= 0.0) && (n < 20.0)	  && (y > 6.0) && (y < 20.0) )	{	  /* Recur backwards from a larger value of n	   */	rlarger:	  k = n;	  y = y + an + 1.0;	  if( y < 30.0 )	    y = 30.0;	  y = n + floor(y-n);	  q = recur( &y, x, &k, 0 );	  y = jvs(y,x) * q;	  goto done;	}      if( k <= 30.0 )	{	  k = 2.0;	}      else if( k < 90.0 )	{	  k = (3*k)/4;	}      if( an > (k + 3.0) )	{	  if( n < 0.0 )	    k = -k;	  q = n - floor(n);	  k = floor(k) + q;	  if( n > 0.0 )	    q = recur( &n, x, &k, 1 );	  else	    {	      t = k;	      k = n;	      q = recur( &t, x, &k, 1 );	      k = t;	    }	  if( q == 0.0 )	    {	    underf:	      y = 0.0;	      goto done;	    }	}      else	{	  k = n;	  q = 1.0;	}      /* boundary between convergence of       * power series and Hankel expansion       */      y = fabs(k);      if( y < 26.0 )	t = (0.0083*y + 0.09)*y + 12.9;      else	t = 0.9 * y;      if( x > t )	y = hankel(k,x);      else	y = jvs(k,x);      if( n > 0.0 )	y /= q;      else	y *= q;    }  else    {      /* For large n, use the uniform expansion       * or the transitional expansion.       * But if x is of the order of n**2,       * these may blow up, whereas the       * Hankel expansion will then work.       */      if( n < 0.0 )	{	  it_warning("besselj:: partial loss of precision");	  y = 0.0;	  goto done;	}      t = x/n;      t /= n;      if( t > 0.3 )	y = hankel(n,x);      else	y = jnx(n,x);    } done:	return( sign * y);}/* Reduce the order by backward recurrence. * AMS55 #9.1.27 and 9.1.73. */static double recur(double *n, double x, double *newn, int cancel){  double pkm2, pkm1, pk, qkm2, qkm1;  /* double pkp1; */  double k, ans, qk, xk, yk, r, t, kf;  static double big = BIG;  int nflag, ctr;  /* continued fraction for Jn(x)/Jn-1(x)  */  if( *n < 0.0 )    nflag = 1;  else    nflag = 0; fstart:  pkm2 = 0.0;  qkm2 = 1.0;  pkm1 = x;  qkm1 = *n + *n;  xk = -x * x;  yk = qkm1;  ans = 1.0;  ctr = 0;  do    {      yk += 2.0;      pk = pkm1 * yk +  pkm2 * xk;      qk = qkm1 * yk +  qkm2 * xk;      pkm2 = pkm1;      pkm1 = pk;      qkm2 = qkm1;      qkm1 = qk;      if( qk != 0 )	r = pk/qk;      else	r = 0.0;      if( r != 0 )	{	  t = fabs( (ans - r)/r );	  ans = r;	}      else	t = 1.0;      if( ++ctr > 1000 )	{	  it_warning("besselj:: Underflow");	  //mtherr( "jv", UNDERFLOW );	  goto done;	}      if( t < MACHEP )	goto done;      if( fabs(pk) > big )	{	  pkm2 /= big;	  pkm1 /= big;	  qkm2 /= big;	  qkm1 /= big;	}    }  while( t > MACHEP ); done:  /* Change n to n-1 if n < 0 and the continued fraction is small   */  if( nflag > 0 )    {      if( fabs(ans) < 0.125 )	{	  nflag = -1;	  *n = *n - 1.0;	  goto fstart;	}    }  kf = *newn;  /* backward recurrence   *              2k   *  J   (x)  =  --- J (x)  -  J   (x)   *   k-1         x   k         k+1   */  pk = 1.0;  pkm1 = 1.0/ans;  k = *n - 1.0;  r = 2 * k;  do    {      pkm2 = (pkm1 * r  -  pk * x) / x;      /*	pkp1 = pk; */      pk = pkm1;      pkm1 = pkm2;      r -= 2.0;      /*	t = fabs(pkp1) + fabs(pk);	if( (k > (kf + 2.5)) && (fabs(pkm1) < 0.25*t) )	{	k -= 1.0;	t = x*x;	pkm2 = ( (r*(r+2.0)-t)*pk - r*x*pkp1 )/t;	pkp1 = pk;	pk = pkm1;	pkm1 = pkm2;	r -= 2.0;	}      */      k -= 1.0;    }  while( k > (kf + 0.5) );  /* Take the larger of the last two iterates   * on the theory that it may have less cancellation error.   */  if( cancel )    {      if( (kf >= 0.0) && (fabs(pk) > fabs(pkm1)) )	{	  k += 1.0;	  pkm2 = pk;	}    }  *newn = k;  return( pkm2 );}/* Ascending power series for Jv(x). * AMS55 #9.1.10. */static double jvs(double n, double x){  double t, u, y, z, k;  int ex;  z = -x * x / 4.0;  u = 1.0;  y = u;  k = 1.0;  t = 1.0;  while( t > MACHEP )    {      u *= z / (k * (n+k));      y += u;      k += 1.0;      if( y != 0 )	t = fabs( u/y );    }  t = frexp( 0.5*x, &ex );