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