home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / ssp / spclfnct / besj.for next >
Encoding:
Text File  |  1985-11-29  |  2.6 KB  |  106 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE BESJ
  5. C
  6. C        PURPOSE
  7. C           COMPUTE THE J BESSEL FUNCTION FOR A GIVEN ARGUMENT AND ORDER
  8. C
  9. C        USAGE
  10. C           CALL BESJ(X,N,BJ,D,IER)
  11. C
  12. C        DESCRIPTION OF PARAMETERS
  13. C           X  -THE ARGUMENT OF THE J BESSEL FUNCTION DESIRED
  14. C           N  -THE ORDER OF THE J BESSEL FUNCTION DESIRED
  15. C           BJ -THE RESULTANT J BESSEL FUNCTION
  16. C           D  -REQUIRED ACCURACY
  17. C           IER-RESULTANT ERROR CODE WHERE
  18. C              IER=0  NO ERROR
  19. C              IER=1  N IS NEGATIVE
  20. C              IER=2  X IS NEGATIVE OR ZERO
  21. C              IER=3  REQUIRED ACCURACY NOT OBTAINED
  22. C              IER=4  RANGE OF N COMPARED TO X NOT CORRECT (SEE REMARKS)
  23. C
  24. C        REMARKS
  25. C           N MUST BE GREATER THAN OR EQUAL TO ZERO, BUT IT MUST BE
  26. C           LESS THAN
  27. C              20+10*X-X** 2/3   FOR X LESS THAN OR EQUAL TO 15
  28. C              90+X/2           FOR X GREATER THAN 15
  29. C
  30. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  31. C           NONE
  32. C
  33. C        METHOD
  34. C           RECURRENCE RELATION TECHNIQUE DESCRIBED BY H. GOLDSTEIN AND
  35. C           R.M. THALER,'RECURRENCE TECHNIQUES FOR THE CALCULATION OF
  36. C           BESSEL FUNCTIONS',M.T.A.C.,V.13,PP.102-108 AND I.A. STEGUN
  37. C           AND M. ABRAMOWITZ,'GENERATION OF BESSEL FUNCTIONS ON HIGH
  38. C           SPEED COMPUTERS',M.T.A.C.,V.11,1957,PP.255-257
  39. C
  40. C     ..................................................................
  41. C
  42.       SUBROUTINE BESJ(X,N,BJ,D,IER)
  43. C
  44.       BJ=.0
  45.       IF(N)10,20,20
  46.    10 IER=1
  47.       RETURN
  48.    20 IF(X)30,30,31
  49.    30 IER=2
  50.       RETURN
  51.    31 IF(X-15.)32,32,34
  52.    32 NTEST=20.+10.*X-X** 2/3
  53.       GO TO 36
  54.    34 NTEST=90.+X/2.
  55.    36 IF(N-NTEST)40,38,38
  56.    38 IER=4
  57.       RETURN
  58.    40 IER=0
  59.       N1=N+1
  60.       BPREV=.0
  61. C
  62. C     COMPUTE STARTING VALUE OF M
  63. C
  64.       IF(X-5.)50,60,60
  65.    50 MA=X+6.
  66.       GO TO 70
  67.    60 MA=1.4*X+60./X
  68.    70 MB=N+IFIX(X)/4+2
  69.       MZERO=MAX0(MA,MB)
  70. C
  71. C     SET UPPER LIMIT OF M
  72. C
  73.       MMAX=NTEST
  74.   100 DO 190 M=MZERO,MMAX,3
  75. C
  76. C     SET F(M),F(M-1)
  77. C
  78.       FM1=1.0E-28
  79.       FM=.0
  80.       ALPHA=.0
  81.       IF(M-(M/2)*2)120,110,120
  82.   110 JT=-1
  83.       GO TO 130
  84.   120 JT=1
  85.   130 M2=M-2
  86.       DO 160 K=1,M2
  87.       MK=M-K
  88.       BMK=2.*FLOAT(MK)*FM1/X-FM
  89.       FM=FM1
  90.       FM1=BMK
  91.       IF(MK-N-1)150,140,150
  92.   140 BJ=BMK
  93.   150 JT=-JT
  94.       S=1+JT
  95.   160 ALPHA=ALPHA+BMK*S
  96.       BMK=2.*FM1/X-FM
  97.       IF(N)180,170,180
  98.   170 BJ=BMK
  99.   180 ALPHA=ALPHA+BMK
  100.       BJ=BJ/ALPHA
  101.       IF(ABS(BJ-BPREV)-ABS(D*BJ))200,200,190
  102.   190 BPREV=BJ
  103.       IER=3
  104.   200 RETURN
  105.       END
  106.