C
C     ..................................................................
C
C        SUBROUTINE BESJ
C
C        PURPOSE
C           COMPUTE THE J BESSEL FUNCTION FOR A GIVEN ARGUMENT AND ORDER
C
C        USAGE
C           CALL BESJ(X,N,BJ,D,IER)
C
C        DESCRIPTION OF PARAMETERS
C           X  -THE ARGUMENT OF THE J BESSEL FUNCTION DESIRED
C           N  -THE ORDER OF THE J BESSEL FUNCTION DESIRED
C           BJ -THE RESULTANT J BESSEL FUNCTION
C           D  -REQUIRED ACCURACY
C           IER-RESULTANT ERROR CODE WHERE
C              IER=0  NO ERROR
C              IER=1  N IS NEGATIVE
C              IER=2  X IS NEGATIVE OR ZERO
C              IER=3  REQUIRED ACCURACY NOT OBTAINED
C              IER=4  RANGE OF N COMPARED TO X NOT CORRECT (SEE REMARKS)
C
C        REMARKS
C           N MUST BE GREATER THAN OR EQUAL TO ZERO, BUT IT MUST BE
C           LESS THAN
C              20+10*X-X** 2/3   FOR X LESS THAN OR EQUAL TO 15
C              90+X/2           FOR X GREATER THAN 15
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           RECURRENCE RELATION TECHNIQUE DESCRIBED BY H. GOLDSTEIN AND
C           R.M. THALER,'RECURRENCE TECHNIQUES FOR THE CALCULATION OF
C           BESSEL FUNCTIONS',M.T.A.C.,V.13,PP.102-108 AND I.A. STEGUN
C           AND M. ABRAMOWITZ,'GENERATION OF BESSEL FUNCTIONS ON HIGH
C           SPEED COMPUTERS',M.T.A.C.,V.11,1957,PP.255-257
C
C     ..................................................................
C
      SUBROUTINE BESJ(X,N,BJ,D,IER)
C
      BJ=.0
      IF(N)10,20,20
   10 IER=1
      RETURN
   20 IF(X)30,30,31
   30 IER=2
      RETURN
   31 IF(X-15.)32,32,34
   32 NTEST=20.+10.*X-X** 2/3
      GO TO 36
   34 NTEST=90.+X/2.
   36 IF(N-NTEST)40,38,38
   38 IER=4
      RETURN
   40 IER=0
      N1=N+1
      BPREV=.0
C
C     COMPUTE STARTING VALUE OF M
C
      IF(X-5.)50,60,60
   50 MA=X+6.
      GO TO 70
   60 MA=1.4*X+60./X
   70 MB=N+IFIX(X)/4+2
      MZERO=MAX0(MA,MB)
C
C     SET UPPER LIMIT OF M
C
      MMAX=NTEST
  100 DO 190 M=MZERO,MMAX,3
C
C     SET F(M),F(M-1)
C
      FM1=1.0E-28
      FM=.0
      ALPHA=.0
      IF(M-(M/2)*2)120,110,120
  110 JT=-1
      GO TO 130
  120 JT=1
  130 M2=M-2
      DO 160 K=1,M2
      MK=M-K
      BMK=2.*FLOAT(MK)*FM1/X-FM
      FM=FM1
      FM1=BMK
      IF(MK-N-1)150,140,150
  140 BJ=BMK
  150 JT=-JT
      S=1+JT
  160 ALPHA=ALPHA+BMK*S
      BMK=2.*FM1/X-FM
      IF(N)180,170,180
  170 BJ=BMK
  180 ALPHA=ALPHA+BMK
      BJ=BJ/ALPHA
      IF(ABS(BJ-BPREV)-ABS(D*BJ))200,200,190
  190 BPREV=BJ
      IER=3
  200 RETURN
      END