/* struve.c * * Struve function * * * * SYNOPSIS: * * double v, x, y, struve(); * * y = struve( v, x ); * * * * DESCRIPTION: * * Computes the Struve function Hv(x) of order v, argument x. * Negative x is rejected unless v is an integer. * * This module also contains the hypergeometric functions 1F2 * and 3F0 and a routine for the Bessel function Yv(x) with * noninteger v. * * * * ACCURACY: * * Not accurately characterized, but spot checked against tables. * */ /* Cephes Math Library Release 2.81: June, 2000 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier */ #include <math.h> #define DEBUG 0 #ifdef ANSIPROT extern double gamma ( double ); extern double pow ( double, double ); extern double sqrt ( double ); extern double yn ( int, double ); extern double jv ( double, double ); extern double fabs ( double ); extern double floor ( double ); extern double sin ( double ); extern double cos ( double ); double yv ( double, double ); double onef2 (double, double, double, double, double * ); double threef0 (double, double, double, double, double * ); #else double gamma(), pow(), sqrt(), yn(), yv(), jv(), fabs(), floor(); double sin(), cos(); double onef2(), threef0(); #endif static double stop = 1.37e-17; extern double MACHEP; double onef2( a, b, c, x, err ) double a, b, c, x; double *err; { double n, a0, sum, t; double an, bn, cn, max, z; an = a; bn = b; cn = c; a0 = 1.0; sum = 1.0; n = 1.0; t = 1.0; max = 0.0; do { if( an == 0 ) goto done; if( bn == 0 ) goto error; if( cn == 0 ) goto error; if( (a0 > 1.0e34) || (n > 200) ) goto error; a0 *= (an * x) / (bn * cn * n); sum += a0; an += 1.0; bn += 1.0; cn += 1.0; n += 1.0; z = fabs( a0 ); if( z > max ) max = z; if( sum != 0 ) t = fabs( a0 / sum ); else t = z; } while( t > stop ); done: *err = fabs( MACHEP*max /sum ); #if DEBUG printf(" onef2 cancellation error %.5E\n", *err ); #endif goto xit; error: #if DEBUG printf("onef2 does not converge\n"); #endif *err = 1.0e38; xit: #if DEBUG printf("onef2( %.2E %.2E %.2E %.5E ) = %.3E %.6E\n", a, b, c, x, n, sum); #endif return(sum); } double threef0( a, b, c, x, err ) double a, b, c, x; double *err; { double n, a0, sum, t, conv, conv1; double an, bn, cn, max, z; an = a; bn = b; cn = c; a0 = 1.0; sum = 1.0; n = 1.0; t = 1.0; max = 0.0; conv = 1.0e38; conv1 = conv; do { if( an == 0.0 ) goto done; if( bn == 0.0 ) goto done; if( cn == 0.0 ) goto done; if( (a0 > 1.0e34) || (n > 200) ) goto error; a0 *= (an * bn * cn * x) / n; an += 1.0; bn += 1.0; cn += 1.0; n += 1.0; z = fabs( a0 ); if( z > max ) max = z; if( z >= conv ) { if( (z < max) && (z > conv1) ) goto done; } conv1 = conv; conv = z; sum += a0; if( sum != 0 ) t = fabs( a0 / sum ); else t = z; } while( t > stop ); done: t = fabs( MACHEP*max/sum ); #if DEBUG printf(" threef0 cancellation error %.5E\n", t ); #endif max = fabs( conv/sum ); if( max > t ) t = max; #if DEBUG printf(" threef0 convergence %.5E\n", max ); #endif goto xit; error: #if DEBUG printf("threef0 does not converge\n"); #endif t = 1.0e38; xit: #if DEBUG printf("threef0( %.2E %.2E %.2E %.5E ) = %.3E %.6E\n", a, b, c, x, n, sum); #endif *err = t; return(sum); } extern double PI; double struve( v, x ) double v, x; { double y, ya, f, g, h, t; double onef2err, threef0err; f = floor(v); if( (v < 0) && ( v-f == 0.5 ) ) { y = jv( -v, x ); f = 1.0 - f; g = 2.0 * floor(f/2.0); if( g != f ) y = -y; return(y); } t = 0.25*x*x; f = fabs(x); g = 1.5 * fabs(v); if( (f > 30.0) && (f > g) ) { onef2err = 1.0e38; y = 0.0; } else { y = onef2( 1.0, 1.5, 1.5+v, -t, &onef2err ); } if( (f < 18.0) || (x < 0.0) ) { threef0err = 1.0e38; ya = 0.0; } else { ya = threef0( 1.0, 0.5, 0.5-v, -1.0/t, &threef0err ); } f = sqrt( PI ); h = pow( 0.5*x, v-1.0 ); if( onef2err <= threef0err ) { g = gamma( v + 1.5 ); y = y * h * t / ( 0.5 * f * g ); return(y); } else { g = gamma( v + 0.5 ); ya = ya * h / ( f * g ); ya = ya + yv( v, x ); return(ya); } } /* Bessel function of noninteger order */ double yv( v, x ) double v, x; { double y, t; int n; y = floor( v ); if( y == v ) { n = v; y = yn( n, x ); return( y ); } t = PI * v; y = (cos(t) * jv( v, x ) - jv( -v, x ))/sin(t); return( y ); } /* Crossover points between ascending series and asymptotic series * for Struve function * * v x * * 0 19.2 * 1 18.95 * 2 19.15 * 3 19.3 * 5 19.7 * 10 21.35 * 20 26.35 * 30 32.31 * 40 40.0 */