/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:32:03 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "sbesyn.h" #include #include /* PARAMETER translations */ #define C11293 1.1293e0 #define C16 16.e0 #define C1P5 1.5e0 #define C2BYPI 0.63661977236758134307553505349005744e0 #define C59 0.59e0 #define CP6 0.60206e0 #define CUTLOW 0.012e0 #define EC 0.57721566490153286060651209008240243e0 #define FOUR 4.0e0 #define HALF 0.5e0 #define HALFPI 1.57079632679489661923132169163975144e0 #define LDIMA 95 #define LOGPI 1.14472988584940017414342735135305869e0 #define LOGTWO 0.69314718055994530941723212145817657e0 #define ND2 10 #define ONE 1.e0 #define PI 3.14159265358979323846264338327950288e0 #define THREE 3.0e0 #define THSJ 0.12e0 #define TWO 2.e0 #define ZERO 0.e0 /* end of PARAMETER translations */ void /*FUNCTION*/ sbesyn( float x, float alpha, long num, float by[]) { LOGICAL32 flag, j1smal; long int _l0, i, ii, k, m, mu, musave, nmax, nmin; float aj[LDIMA], arg, bgam, bgamsq, chi, d1, d2, d2val, d3, dr, em, em1, emu, en1, en2, en3, eta, fac, fk, fkm1, fkp1, fv, fvm1, fvp1, g, g0, g1, gmax, gnu, p, piv2, psi1, psiz, q, q2dxpv, scale, sum, temp, test, thvdx, twodx, v, v2, vpmu, xlog, ylog, yv, yvm1, yvp1, z; static float big, biglog, hicut, small, xpq; static float eps = ZERO; static float d2ser[ND2-(0)+1]={.3674669051966159615185e00,-.1782903980807269842231e01, .9411644685122855908427e00,-.1958865250248747878077e01,.1557306621108283294475e01, -.2521051413546812096437e01,.2184502685635110914145e01,-.3140067153452674402872e01, .2817401038921461361158e01,-.3772198775799670081858e01,.3452780604492584575060e01}; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Aj = &aj[0] - 1; float *const By = &by[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. *>> 2009-10-28 SBESYN Krogh/Snyder Moved BGAMSQ = BGAM**2 inside else. *>> 2000-01-05 SBESYN Krogh Changed c$ comments to c. (Cray f90 prob.) *>> 1995-11-13 SBESYN Krogh Converted from SFTRAN to Fortran *>> 1994-10-19 SBESYN Krogh Changes to use M77CON *>> 1994-04-19 SBESYN CLL Edited to make DP & SP files similar. *>> 1991-01-14 SBESYN CLL Removed duplicate data statement. *>> 1989-08-09 SBESYN CLL More accurate constants for Cray *>> 1986-03-18 SBESYN Lawson Initial code. *--S replaces "?": ?BESYN, ?BESPQ, ?ERV1, ?GAMMA, ?LGAMA * * This subr computes the Y Bessel functions of X for * NUM orders from ALPHA through ALPHA + NUM - 1. * The results will be stored in BY(I), I = 1,...,NUM. * * Require X .gt. 0., ALPHA .ge. 0., NUM .ge. 1. * * The original subroutines SBYNU and BESJ/BESY were * designed and programmed by E. W. Ng and W. V. Snyder, JPL, * in 1973. Modified by Ng and S. Singletary, JPL, 1974. * * 1984 April 16, JPL, C. Lawson and S. Chan. Original subrs * combined into single subroutine. Converted to SFTRAN3 * and Fortran 77 for the JPL MATH77 library. * * ------------------------------------------------------------------ * * As the machine precision, EPS, increases so does XPQ, * and thus so does the requirement for the storage dimension, * LDIMA. Here are some values of -LOG10(EPS), XPQ, and LDIMA. * * -LOG10(EPS)= 5 10 15 20 25 30 * XPQ = 5.74 11.38 17.03 22.68 28.32 33.97 * LDIMA = 17 33 49 63 79 95 * * Since LDIMA cannot be adjusted at run-time, we are * setting it to 95 to handle the case of CDC double precision * which is probably the largest precision system likely * to be encountered. Note that the relative precision of results * from this subr cannot exceed about 16 or 17 decimal places * because of the limited accuracy of the polynomial coeffs used to * compute G0, and also the limited precision of the * subprograms referenced for gamma and loggamma. * ------------------------------------------------------------------ * Subprograms used: TANH, LOG10, LOG, EXP, COS, SIN, SQRT */ /* ---------- */ /* LOGPI = ln(pi) */ /* C2BYPI = 2/pi */ /* EC = Euler's constant. */ /* LOGTWO = Ln(2) */ /* ------------------------------------------------------------------ * * Set environmental parameters. * */ if (eps == ZERO) { eps = FLT_EPSILON/FLT_RADIX; hicut = ONE/(eps*C16); small = C16*FLT_MIN/eps; xpq = C11293*(CP6 - log10f( eps )) - C59; big = FLT_MAX/TWO; biglog = logf( big ); } /* ------------------------------------------------------------------ * * Compute V, NMIN, and NMAX. * */ nmin = (long)( alpha ); v = alpha - (float)( nmin ); nmax = nmin + num - 1; /* ------------------------------------------------------------------ * * Test validity of given argument values. * */ if ((x < ZERO || alpha < ZERO) || num < 1) { /* Error 1. */ ermsg( "SBESYN", 1, 0, "Require X .gt. 0, ALPHA .ge. 0, NUM .ge. 1" , ',' ); goto L_400; /* ------------------------------------------------------------------ * * Branch on size of X. * */ } else if (x == ZERO) { /* Error 6. */ for (i = 1; i <= num; i++) { By[i] = -big; } ermsg( "SBESYN", 6, 0, "When X = 0., function value is -INFINITY." , ',' ); goto L_400; } else if (x < eps) { /* ********************* Code for very small X case. ******************** * * Use a single term expression for Y, valid for X very close to zero. * Ref NBS AMS 55 Eqs 9.1.8 & 9.1.9. * For GNU = 0, Y = (2/pi) * (EC + Ln(X/2)), {EC = Euler's const.} * For GNU .gt. 0 Y = -(1/pi) * Gamma(GNU) * (X/2)**(-GNU) * */ xlog = logf( x ); gnu = alpha; for (i = 1; i <= num; i++) { if (gnu == ZERO) { By[i] = C2BYPI*(EC + xlog - LOGTWO); } else { ylog = slgama( gnu ) - gnu*(xlog - LOGTWO) - LOGPI; if (ylog < biglog) { By[i] = -expf( ylog ); } else { /* Error 5. */ for (ii = i; ii <= num; ii++) { By[ii] = -big; } ermsg( "SBESYN", 5, 0, "Results exceed overflow limit from BY(I) on." , ',' ); ierv1( "I", i, ',' ); goto L_400; } } gnu += ONE; } return; } else { twodx = TWO/x; if (x <= 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). * */ mu = (long)( x ) + 1; dr = twodx*(v + (float)( mu )); fkp1 = ONE; fk = ZERO; /* RECUR FORWARD UNTIL FKP1 IS GREATER THAN PRECISION OF ARITHMETIC. * */ L_210: if (eps*fabsf( fkp1 ) <= ONE) { mu += 1; dr += twodx; fkm1 = fk; fk = fkp1; fkp1 = dr*fk - fkm1; goto L_210; } /* WE ARE NOW ASSURED THAT BACKWARD RECURRENCE FROM MU WILL YIELD * ACCURATE RESULTS. * * GUARANTEE EVEN MU */ if ((mu%2) != 0) mu += 1; musave = mu; /* Test for Error 3. * This error should never happen. Large MU would be due to * large X. But X is not larger than XPQ here. * See explanation at the beginning of this subroutine * of the relation of XPQ and LDIMA to the machine EPS. * */ if (mu + 1 > LDIMA) { ermsg( "SBESYN", 3, 0, "Need larger dimension, LDIMA, to process given X." , ',' ); ermor( "Require LDIMA .ge. MU + 1", ',' ); ierv1( "MU", mu, ',' ); ierv1( "LDIMA", LDIMA, ',' ); goto L_400; } fvm1 = small; Aj[mu + 1] = fvm1; 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_230: ; fvp1 = fv; fv = fvm1; if (fabsf( fv ) > test) { /* Rescale */ fv /= sum; fvp1 /= sum; for (ii = mu + 1; ii <= musave; ii++) { Aj[ii] /= sum; } sum = ONE; } fvm1 = fac*fv - fvp1; mu -= 1; emu -= ONE; fac = (v + emu)*twodx; Aj[mu + 1] = fvm1; if ((mu%2) == 0) { if (v == ZERO) { sum += fvm1; if (mu == 0) { scale = ONE/sum; goto L_260; } 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; bgam = sgamma( v + ONE ); q2dxpv = powf(twodx,v); scale = (bgam/eta)*sum*q2dxpv; scale = ONE/scale; goto L_260; } } } goto L_230; L_260: ; /* NORMALIZE AJ() TO GET VALUES OF J-BESSEL FUNCTION. * */ for (i = 1; i <= (musave + 1); i++) { Aj[i] *= scale; } mu = musave; /* Compute Y Bessel functions, making use of previously computed J * Bessel functions and other previously computed values, MU, BGAM, * Q2DXPV, TWODX, V. * * Here V is in the range [0.,1.). The quantities G0 and G1 depend * on X and V and are unbounded as V approaches 1. Therefore we * make the the change of variables * V2 = V if V .le. 0.5 and * V2 = 1 - V if V .gt. 0.5 * Then G0 and G1 are computed as functions of X and V2 with V2 in * the range (-0.5, 0.5]. * * Compute G0 and G1. * */ v2 = v; if (v == ZERO) { z = EC - logf( twodx ); g0 = z/HALFPI; g1 = TWO/HALFPI; bgamsq = ONE; q2dxpv = ONE; } else { bgamsq = SQ(bgam); if (v > HALF) { /* Use the transformed variable, V2 = V - 1. * Make corresponding transformation of Q2DXPV & BGAMSQ. * */ v2 = (v - HALF) - HALF; q2dxpv /= twodx; bgamsq = bgamsq/SQ(v); } piv2 = PI*v2; /* Here V2 is in [-.5, .5]. Now test against CUTLOW = 0.012 * */ if (fabsf( v2 ) < CUTLOW) { /* Here we compute * G0 = (ONE / TAN(PIV2)) - Q2DXPV**2 * BGAMSQ / PIV2 * by a formulation that retains accuracy for * V2 close to zero. * > The no. of coeffs from D2SER() used to compute D2VAL * could be fewer on lower precision computers, however this * computation is only done about 2.4% of the time so the * potential time saving would probably not be noticeable. * * This method was derived by C. Lawson and W. V. Snyder, * JPL, 1984 Apr 15. * * First compute EM1 = (2/X)**(2*V2) - 1 * = exp(2 * V2 * log(2/X)) - 1 * */ arg = TWO*v2*logf( twodx ); if (fabsf( arg ) < LOGTWO) { temp = tanhf( HALF*arg ); em1 = TWO*temp/(ONE - temp); } else { em1 = expf( arg ) - ONE; } /* Evaluate taylor series for * D2VAL = (PIV2 * cotan(PIV2) - BGAMSQ) / PIV2 * */ d2val = d2ser[ND2]; for (i = ND2 - 1; i >= 0; i--) { d2val = d2ser[i] + v2*d2val; } g0 = d2val - bgamsq*(em1/piv2); g1 = (SQ(q2dxpv)/HALFPI)*bgamsq*(TWO + v2)/(ONE - v2); } else { g0 = (ONE/tanf( piv2 )) - SQ(q2dxpv)*bgamsq/piv2; g1 = (SQ(q2dxpv)/HALFPI)*bgamsq*(TWO + v2)/(ONE - v2); } } /* ---------------------------------- * * COMPUTE YO FROM SUM(J'S) FORM * */ en3 = v2 + ONE; en2 = v2 + en3; en1 = v2 + FOUR; d1 = TWO; d2 = d1 - v2; d3 = d1 + v2; flag = FALSE; /* THSJ = 0.12 */ j1smal = fabsf( Aj[1] ) < THSJ; if (j1smal || v2 < ZERO) { flag = TRUE; /* Y(V2+1,X) MUST ALSO BE COMPUTED BY A SUM * */ thvdx = THREE*v2/x; psiz = -bgamsq*SQ(q2dxpv)/(HALFPI*x); psi1 = g0 - HALF*g1; } if (v2 >= ZERO) { m = 3; yv = g0*Aj[1]; if (j1smal) { yvp1 = psiz*Aj[1] + psi1*Aj[2]; } } else { z = twodx*v*Aj[1] - Aj[2]; yv = g0*z; m = 2; yvp1 = psiz*z + psi1*Aj[1]; } for (i = m; i <= mu; i += 2) { yv += g1*Aj[i]; g = g1; g1 *= -(en1/d1)*(en2/d2)*(en3/d3); en1 += TWO; en2 += ONE; en3 += ONE; d1 += ONE; d2 += ONE; d3 += TWO; if (flag) { yvp1 += thvdx*g*Aj[i] + HALF*(g - g1)*Aj[i + 1]; } } if (v2 < ZERO) { z = yvp1; yvp1 = v*z*twodx - yv; yv = z; } else if (!j1smal) { /* NOW COMPUTE Y(V+1) * WRONSKIAN PROVIDED NOT NEAR A ZERO OF J * */ yvp1 = (yv*Aj[2] - ONE/(x*HALFPI))/Aj[1]; } goto L_350; } 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 Y(V,X) and Y(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; yv = sqrtf( ONE/(HALFPI*x) )*(p*sinf( chi ) + q*cosf( chi )); if (nmax > 0) { sbespq( x, v + ONE, &p, &q ); chi = x - (v + C1P5)*HALFPI; yvp1 = sqrtf( ONE/(HALFPI*x) )*(p*sinf( chi ) + q*cosf( chi )); } goto L_350; } else { /* Error 2. */ ermsg( "SBESYN", 2, 0, "Cannot obtain any accuracy when X exceeds HICUT." , ',' ); serv1( "HICUT", hicut, ',' ); goto L_400; } } L_350: ; /* Do forward recursion * Given YV = Y(V,X), YVP1 = Y(V+1,X), TWODX = 2/X, NMIN, NUM, NMAX = * NMIN + NUM -1, X, ALPHA, and BIG. Recur forward and store * Y(NMIN+V) thru Y(NMAX+V) in BY(1) thru BY(NUM). * */ if (nmin == 0) { By[1] = yv; if (nmax > 0) { By[2] = yvp1; } } else if (nmin == 1) { By[1] = yvp1; } if (nmax > 1) { g = v*twodx; gmax = g + twodx*(float)( nmax - 1 ); test = big/fmaxf( ONE, gmax ); /* Note: In the following statement, 3-NMIN can be nonpositive. * */ for (k = 3 - nmin; k <= num; k++) { yvm1 = yv; yv = yvp1; if (fabsf( yv ) > test) { /* The recursion has reached the overflow limit. * Set remaining elts of BY() to a large negative value * and issue error message. * */ for (ii = max( k, 1 ); ii <= num; ii++) { By[ii] = -big; } /* Error 4. */ ermsg( "SBESYN", 4, 0, "Results exceed overflow limit from BY(I) on." , ',' ); ierv1( "I", max( k, 1 ), ',' ); goto L_400; } g += twodx; yvp1 = g*yv - yvm1; if (k >= 1) By[k] = yvp1; } } return; /* Error return */ L_400: ; serv1( "X", x, ',' ); serv1( "ALPHA", alpha, ',' ); ierv1( "NUM", num, '.' ); return; } /* end of function */