/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:57 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "sbesjn.h" #include #include /* PARAMETER translations */ #define C11293 1.1293e0 #define C16 16.e0 #define C1P5 1.5e0 #define C59 0.59e0 #define CP6 0.60206e0 #define HALF 0.5e0 #define HALFPI 1.57079632679489661923132169163975144e0 #define ONE 1.e0 #define TENTH 0.1e0 #define TWO 2.e0 #define ZERO 0.e0 /* end of PARAMETER translations */ void /*FUNCTION*/ sbesjn( float x, float alpha, long num, float bj[]) { long int _l0, i, i1, ii, j, k, m, mu, nmax, nmin; float besv, besvm1, besvp1, bjnu, bot, chi, dr, em, emu, eta, fac, fac2, fk, fkm1, fkp1, fv, fvm1, fvp1, g, gnu, halfx, p, q, scale, sum, t1, t2, term, term1, test, top, twodx, v, vpmu; static float big, hicut, small, xpq; static float eps = ZERO; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Bj = &bj[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 1996-03-30 SBESJN Krogh Removed INT from type statement. *>> 1995-11-13 SBESJN Krogh Converted SFTRAN to Fortran *>> 1995-10-24 SBESJN Krogh Removed blanks in numbers for C conversion. *>> 1994-10-19 SBESJN Krogh Changes to use M77CON *>> 1994-04-19 SBESJN CLL Edited to make DP & SP files similar. *>> 1992-03-13 SBESJN FTK Removed implicit statements. *>> 1989-08-09 SBESJN CLL More accurate HALFPI for Cray *>> 1986-03-18 SBESJN Lawson Initial code. *--S replaces "?": ?BESJN, ?GAMMA, ?BESPQ, ?ERV1 * * COMPUTE THE J-BESSEL FUNCTIONS OF X FOR ORDER ALPHA THROUGH * ALPHA + NUM - 1 in steps of one. * STORE THE RESULT INTO BJ(I),I=1,...,NUM. * Require X .ge. 0, ALPHA .ge. 0, NUM .ge. 1 * * Original code due to E. W. Ng and W. V. Snyder, JPL, 1973. * Modified by Ng and S. Singletary, 1974. * C.Lawson & S.Chan, JPL, 1984 Apr 05: * Changed calling sequence. * Adapted to SFTRAN3 and Fortran 77. * Added code to avoid overflow during the recurrsion. * Added code to use Taylor series for small x. * March 23. Improved accuracy of backward recursion. * Changed to use P and Q instead of R and THETA * in the region of large X. * ------------------------------------------------------------------ * */ /* ---------- */ /* ---------- */ /* ------------------------------------------------------------------ * * Set environment parameters. This should happen only * the first time this subroutine is called. * */ if (eps == ZERO) { eps = FLT_EPSILON/FLT_RADIX; hicut = ONE/(eps*C16); xpq = C11293*(CP6 - log10f( eps )) - C59; small = C16*FLT_MIN/eps; big = FLT_MAX/TWO; } /* Test validity of input values. * */ if ((x < ZERO || alpha < ZERO) || num < 1) { /* Error 1. */ ermsg( "SBESJN", 1, 0, "REQUIRE X.GE.0, ALPHA.GE.0, NUM.GE.1" , ',' ); } else { /* Begin computation. */ nmin = (long)( alpha ); v = alpha - (float)( nmin ); nmax = nmin + num - 1; if (x <= TENTH) { /* ********************* Code for the small X case. ********************* * */ if (x == ZERO) { /* Special for X .eq. 0. */ for (i = 1; i <= num; i++) { Bj[i] = ZERO; } if (alpha == 0) Bj[1] = ONE; } else { /* Here use Taylor series for small x. * */ gnu = alpha; halfx = HALF*x; fac2 = -halfx*halfx; term1 = (powf(halfx,gnu))/sgamma( gnu + ONE ); for (i = 1; i <= num; i++) { /* Sum the series for the Bessel fcn J sub GNU of X given * in Eq 9.1.10, page 360, of AMS 55. * 1984 March 9, JPL, C. L. Lawson. * */ if (term1 == ZERO) { bjnu = ZERO; } else { sum = ZERO; top = fac2; bot = gnu + ONE; t1 = ONE; t2 = bot; term = top/bot; L_100: if (fabsf( term ) > eps) { sum += term; top *= fac2; t1 += ONE; t2 += ONE; bot *= t1*t2; term = top/bot; goto L_100; } bjnu = term1 + term1*sum; } Bj[i] = bjnu; if (bjnu == ZERO) { /* Here current result has underflowed to zero, so we will * set the rest of the results to zero also. * */ for (j = i + 1; j <= num; j++) { Bj[j] = ZERO; } return; } gnu += ONE; term1 *= halfx/gnu; } } return; } else if (x < fmaxf( (float)( nmax + 1 ), xpq )) { /* ********************* Code for the middle X case. ******************** * * * J-TYPE BESSEL FUNCTIONS FOLLOW THE RECURRENCE RELATION * F(V-1,X)=(2*V/X)*F(V,X)-F(V+1,X). * */ twodx = TWO/x; mu = max( (long)( x ) + 1, nmax ); dr = twodx*(v + (float)( mu )); fkp1 = ONE; fk = ZERO; /* RECUR FORWARD UNTIL FKP1 IS GREATER THAN PRECISION OF ARITHMETIC. * */ L_200: if (eps*fabsf( fkp1 ) <= ONE) { mu += 1; dr += twodx; fkm1 = fk; fk = fkp1; fkp1 = dr*fk - fkm1; goto L_200; } /* WE ARE NOW ASSURED THAT BACKWARD RECURRENCE FROM MU WILL YIELD * ACCURATE RESULTS. * * GUARANTEE EVEN MU */ if ((mu%2) != 0) mu += 1; fvm1 = small; fv = ZERO; eta = ONE; sum = fvm1; m = mu/2; em = (float)( m ); emu = (float)( mu ); fac = (v + emu)*twodx; /* Set TEST = largest value that can be multiplied by * FAC without risking overflow. The present value of * FAC is the largest that will occur during the recursion. * TEST will be used to protect against overflow during * the recursion. * */ test = big/fmaxf( ONE, fac ); /* Loop while MU .GT. ZERO */ L_220: ; fvp1 = fv; fv = fvm1; if (fabsf( fv ) > test) { fv /= sum; fvp1 /= sum; i1 = max( 1, mu - nmin + 1 ); for (ii = i1; ii <= num; ii++) { Bj[ii] /= sum; } sum = ONE; } fvm1 = fac*fv - fvp1; mu -= 1; emu -= ONE; fac = (v + emu)*twodx; if (mu >= nmin && mu <= nmax) Bj[mu - nmin + 1] = fvm1; if ((mu%2) == 0) { if (v == ZERO) { sum += fvm1; if (mu == 0) { scale = ONE/sum; goto L_250; } sum += fvm1; } else { if (mu != 0) { vpmu = v + emu; eta *= (em/(v + (em - ONE)))*(vpmu/(vpmu + TWO)); sum += fvm1*eta; em -= ONE; } else { /* Here MU = 0 and EM = 0NE. Thus the expression for * updating ETA reduces to the following simpler * expression. * */ eta /= v + TWO; sum += fvm1*eta; scale = sgamma( v + ONE )/eta*sum*powf(twodx,v); scale = ONE/scale; goto L_250; } } } goto L_220; L_250: ; /* NORMALIZE BJ() TO GET VALUES OF J-BESSEL FUNCTION. * */ for (i = 1; i <= num; i++) { Bj[i] *= scale; } return; } else if (x < hicut) { /* ********************* Code for the large X case. ********************* * * Here we have X .ge. XPQ, and V in [0.,1.). * The asymptotic series for the auxiliary functions P and Q can be * used. From these we will compute J(V,X) and J(V+1,X) and * then recur forward. Reference: NBS AMS 55 Eqs 9.2.5 & 9.2.6 * */ sbespq( x, v, &p, &q ); chi = x - (v + HALF)*HALFPI; besv = sqrtf( ONE/(HALFPI*x) )*(p*cosf( chi ) - q*sinf( chi )); if (nmax > 0) { sbespq( x, v + ONE, &p, &q ); chi = x - (v + C1P5)*HALFPI; besvp1 = sqrtf( ONE/(HALFPI*x) )*(p*cosf( chi ) - q* sinf( chi )); } twodx = TWO/x; /* Given BESV = J(V,X), BESVP1 = J(V+1,X), TWODX = 2/X, * NMIN, NUM, NMAX = NMIN + NUM -1, X, ALPHA, and BIG. * Recur forward and store J(NMIN+V) thru J(NMAX+V) in * BJ(1) thru BJ(NUM). * There should be no overflow posibility in this forward * recursion since NMAX .le. X - 1, and in this region the * magnitude of the J function is less than one. * */ if (nmin == 0) { Bj[1] = besv; if (nmax > 0) { Bj[2] = besvp1; } } else if (nmin == 1) { Bj[1] = besvp1; } if (nmax > 1) { g = v*twodx; /* Note: In the following statement, 3-NMIN can be nonpositive. * */ for (k = 3 - nmin; k <= num; k++) { besvm1 = besv; besv = besvp1; g += twodx; besvp1 = g*besv - besvm1; if (k >= 1) Bj[k] = besvp1; } } return; } else { /* Error 2. */ ermsg( "SBESJN", 2, 0, "Cannot obtain any accuracy when X exceeds HICUT." , ',' ); serv1( "HICUT", hicut, ',' ); } } serv1( "X", x, ',' ); serv1( "ALPHA", alpha, ',' ); ierv1( "NUM", num, '.' ); return; } /* end of function */