1 subroutine SBESJN(X,ALPHA,NUM,BJ) 2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. 3c ALL RIGHTS RESERVED. 4c Based on Government Sponsored Research NAS7-03001. 5c>> 1996-03-30 SBESJN Krogh Removed INT from type statement. 6C>> 1995-11-13 SBESJN Krogh Converted SFTRAN to Fortran 7C>> 1995-10-24 SBESJN Krogh Removed blanks in numbers for C conversion. 8C>> 1994-10-19 SBESJN Krogh Changes to use M77CON 9C>> 1994-04-19 SBESJN CLL Edited to make DP & SP files similar. 10C>> 1992-03-13 SBESJN FTK Removed implicit statements. 11C>> 1989-08-09 SBESJN CLL More accurate HALFPI for Cray 12C>> 1986-03-18 SBESJN Lawson Initial code. 13c--S replaces "?": ?BESJN, ?GAMMA, ?BESPQ, ?ERV1 14C 15C COMPUTE THE J-BESSEL FUNCTIONS OF X FOR ORDER ALPHA THROUGH 16C ALPHA + NUM - 1 in steps of one. 17C STORE THE RESULT INTO BJ(I),I=1,...,NUM. 18C Require X .ge. 0, ALPHA .ge. 0, NUM .ge. 1 19C 20c Original code due to E. W. Ng and W. V. Snyder, JPL, 1973. 21c Modified by Ng and S. Singletary, 1974. 22C C.Lawson & S.Chan, JPL, 1984 Apr 05: 23c Changed calling sequence. 24c Adapted to SFTRAN3 and Fortran 77. 25c Added code to avoid overflow during the recurrsion. 26c Added code to use Taylor series for small x. 27c March 23. Improved accuracy of backward recursion. 28c Changed to use P and Q instead of R and THETA 29c in the region of large X. 30C ------------------------------------------------------------------ 31C 32 save EPS, SMALL, XPQ, BIG, HICUT 33c ---------- 34 external SGAMMA, SBESPQ, R1MACH, ERMSG, SERV1, IERV1 35c ---------- 36 integer I,I1,II,J,K,M,MU,NMAX,NMIN,NUM 37 real R1MACH, SGAMMA 38 real ALPHA,BESV,BESVM1,BESVP1,BIG,BJ(NUM),BJNU,BOT 39 real C11293,C16,C1P5,C59,CHI,CP6 40 real DR,EM,EMU,EPS,ETA,FAC,FAC2,FK,FKM1,FKP1,FV,FVM1 41 real FVP1,G,GNU,HALF,HALFPI,HALFX,HICUT 42 real ONE,P,Q,SCALE,SMALL,SUM,T1,T2,TENTH,TERM 43 real TERM1,TEST,TOP,TWO,TWODX,V,VPMU,X,XPQ,ZERO 44 parameter(ZERO=0.E0, ONE=1.E0, TWO=2.E0, C16=16.E0) 45 parameter( TENTH = 0.1E0, HALF = 0.5E0) 46 parameter(HALFPI = 1.57079632679489661923132169163975144E0) 47 parameter(C11293 = 1.1293E0, CP6 = 0.60206E0, C59 = 0.59E0) 48 parameter(C1P5 = 1.5E0) 49 data EPS / ZERO / 50C ------------------------------------------------------------------ 51C 52c Set environment parameters. This should happen only 53c the first time this subroutine is called. 54c 55 if( EPS .EQ. ZERO ) then 56 EPS = R1MACH(3) 57 HICUT = ONE /(EPS*C16) 58 XPQ = C11293 * (CP6 - log10(EPS)) - C59 59 SMALL = C16 * R1MACH(1) / EPS 60 BIG = R1MACH(2)/TWO 61 end if 62C 63c Test validity of input values. 64c 65 if (X .lt. ZERO .or. ALPHA .lt. ZERO .or. NUM .lt. 1) then 66c 67c Error 1. 68 call ERMSG('SBESJN',1,0, 69 * 'REQUIRE X.GE.0, ALPHA.GE.0, NUM.GE.1',',') 70 else 71c Begin computation. 72 NMIN = int(ALPHA) 73 V = ALPHA - real(NMIN) 74 NMAX = NMIN + NUM - 1 75C 76 if( X .le. TENTH ) then 77c 78c ********************* Code for the small X case. ********************* 79c 80 if( X .eq. ZERO ) then 81c Special for X .eq. 0. 82 do 80 I = 1, NUM 83 BJ(I) = ZERO 84 80 continue 85 if (ALPHA .eq. 0) BJ(1) = ONE 86 else 87c Here use Taylor series for small x. 88C 89 GNU = ALPHA 90 HALFX = HALF*X 91 FAC2 = -HALFX*HALFX 92 TERM1 = (HALFX**GNU) / SGAMMA(GNU + ONE) 93 do 150 I = 1, NUM 94c Sum the series for the Bessel fcn J sub GNU of X given 95c in Eq 9.1.10, page 360, of AMS 55. 96c 1984 March 9, JPL, C. L. Lawson. 97c 98 if(TERM1 .eq. ZERO) then 99 BJNU = ZERO 100 else 101c 102 SUM = ZERO 103 TOP = FAC2 104 BOT = GNU + ONE 105 T1 = ONE 106 T2 = BOT 107 TERM = TOP/BOT 108c 109 100 if ( abs(TERM) .GT. EPS ) then 110 SUM = SUM + TERM 111 TOP = TOP * FAC2 112 T1 = T1 + ONE 113 T2 = T2 + ONE 114 BOT = BOT * T1 * T2 115 TERM = TOP/BOT 116 go to 100 117 end if 118 BJNU = TERM1 + TERM1 * SUM 119 end if 120c 121 BJ(I) = BJNU 122 if( BJNU .eq. ZERO ) then 123C 124c Here current result has underflowed to zero, so we will 125c set the rest of the results to zero also. 126c 127 do 120 J = I+1, NUM 128 BJ(J) = ZERO 129 120 continue 130 return 131 end if 132 GNU = GNU + ONE 133 TERM1 = TERM1 * (HALFX / GNU) 134 150 continue 135 end if 136 return 137 else if( X .lt. MAX(real(NMAX+1), XPQ) ) then 138c 139c ********************* Code for the middle X case. ******************** 140c 141C 142C J-TYPE BESSEL FUNCTIONS FOLLOW THE RECURRENCE RELATION 143C F(V-1,X)=(2*V/X)*F(V,X)-F(V+1,X). 144C 145 TWODX = TWO / X 146 MU = MAX( int(X)+1, NMAX) 147 DR = TWODX * (V+real(MU)) 148 FKP1 = ONE 149 FK = ZERO 150C 151C RECUR FORWARD UNTIL FKP1 IS GREATER THAN PRECISION OF ARITHMETIC. 152C 153 200 if ( EPS * abs(FKP1) .le. ONE ) then 154 MU = MU + 1 155 DR = DR + TWODX 156 FKM1 = FK 157 FK = FKP1 158 FKP1 = DR * FK - FKM1 159 go to 200 160 end if 161C 162C WE ARE NOW ASSURED THAT BACKWARD RECURRENCE FROM MU WILL YIELD 163C ACCURATE RESULTS. 164C 165C GUARANTEE EVEN MU 166 if (MOD(MU,2) .ne. 0) MU = MU + 1 167 FVM1 = SMALL 168 FV = ZERO 169 ETA = ONE 170 SUM = FVM1 171 M = MU / 2 172 EM = real(M) 173 EMU = real(MU) 174 FAC = (V + EMU) * TWODX 175C 176c Set TEST = largest value that can be multiplied by 177c FAC without risking overflow. The present value of 178c FAC is the largest that will occur during the recursion. 179c TEST will be used to protect against overflow during 180c the recursion. 181c 182 TEST = BIG / MAX(ONE, FAC) 183C 184C Loop while MU .GT. ZERO 185 220 continue 186 FVP1 = FV 187 FV = FVM1 188 if( abs(FV) .GT. TEST )then 189 FV = FV / SUM 190 FVP1 = FVP1 / SUM 191 I1 = MAX( 1, MU - NMIN + 1 ) 192 do 230 II = I1, NUM 193 BJ(II) = BJ(II) / SUM 194 230 continue 195 SUM = ONE 196 end if 197 FVM1 = FAC * FV - FVP1 198 MU = MU -1 199 EMU = EMU - ONE 200 FAC = (V + EMU) * TWODX 201 if (MU .ge. NMIN .AND. MU .le. NMAX) BJ(MU-NMIN+1) = FVM1 202 if (MOD(MU,2) .eq. 0) then 203 if (V .eq. ZERO) then 204 SUM = SUM + FVM1 205 if (MU .eq. 0) then 206 SCALE = ONE / SUM 207 go to 250 208 end if 209 SUM = SUM + FVM1 210 else 211 if (MU .ne. 0) then 212 VPMU = V + EMU 213 ETA = ETA * (EM/(V+(EM-ONE)))*(VPMU/(VPMU+TWO)) 214 SUM = SUM + FVM1 * ETA 215 EM = EM - ONE 216 else 217c 218c Here MU = 0 and EM = 0NE. Thus the expression for 219c updating ETA reduces to the following simpler 220c expression. 221c 222 ETA = ETA / (V + TWO) 223 SUM = SUM + FVM1 * ETA 224 SCALE = SGAMMA(V+ONE) / ETA * SUM * TWODX ** V 225 SCALE = ONE / SCALE 226 go to 250 227 end if 228 end if 229 end if 230 go to 220 231 250 continue 232C 233C NORMALIZE BJ() TO GET VALUES OF J-BESSEL FUNCTION. 234C 235 do 260 I = 1, NUM 236 BJ(I) = BJ(I) * SCALE 237 260 continue 238 return 239 else if( X .lt. HICUT ) then 240c 241c ********************* Code for the large X case. ********************* 242c 243c Here we have X .ge. XPQ, and V in [0.,1.). 244c The asymptotic series for the auxiliary functions P and Q can be 245c used. From these we will compute J(V,X) and J(V+1,X) and 246c then recur forward. Reference: NBS AMS 55 Eqs 9.2.5 & 9.2.6 247c 248 call SBESPQ (X,V, P,Q) 249 CHI = X - (V + HALF) * HALFPI 250 BESV = sqrt(ONE / (HALFPI*X)) * (P*COS(CHI) - Q*SIN(CHI)) 251C 252 if( NMAX .GT. 0 ) then 253 call SBESPQ (X,V+ONE, P,Q) 254 CHI = X - (V + C1P5) * HALFPI 255 BESVP1 = sqrt(ONE / (HALFPI*X))*(P*COS(CHI)-Q*SIN(CHI)) 256 end if 257c 258 TWODX = TWO / X 259c 260c Given BESV = J(V,X), BESVP1 = J(V+1,X), TWODX = 2/X, 261c NMIN, NUM, NMAX = NMIN + NUM -1, X, ALPHA, and BIG. 262c Recur forward and store J(NMIN+V) thru J(NMAX+V) in 263c BJ(1) thru BJ(NUM). 264c There should be no overflow posibility in this forward 265c recursion since NMAX .le. X - 1, and in this region the 266c magnitude of the J function is less than one. 267c 268 if( NMIN .eq. 0 ) then 269 BJ(1) = BESV 270 if( NMAX .GT. 0 ) then 271 BJ(2) = BESVP1 272 end if 273 else if( NMIN .eq. 1 ) then 274 BJ(1) = BESVP1 275 end if 276c 277 if( NMAX .GT. 1 ) then 278 G = V * TWODX 279c 280c Note: In the following statement, 3-NMIN can be nonpositive. 281c 282 do 300 K = 3-NMIN, NUM 283 BESVM1 = BESV 284 BESV = BESVP1 285 G = G + TWODX 286 BESVP1 = G * BESV - BESVM1 287 if( K .ge. 1) BJ(K) = BESVP1 288 300 continue 289 end if 290 return 291 else 292c Error 2. 293 call ERMSG('SBESJN', 2, 0, 294 * 'Cannot obtain any accuracy when X exceeds HICUT.', ',') 295 call SERV1('HICUT', HICUT, ',') 296 end if 297 end if 298 call SERV1('X',X,',') 299 call SERV1('ALPHA',ALPHA,',') 300 call IERV1('NUM',NUM,'.') 301 return 302 end 303