8988 REM REVISED APRIL ,1972 ANDY ROTH C.C.G.
8999 ON ERROR GO TO 9500
9000 REM ***** BESSEL ***** MATHEMATICS PROGRAM *****
9001 REM ***** VERSION 1 ****: 7/31/69 *****
9002 REM BESSEL FUNCTION OF THE FIRST KIND
9003 INPUT"WHAT IS THE ORDER";N
9005 INPUT"WHAT IS THE ARGUEMENT";Z
9007 INPUT "WHAT IS THE ACCEPTABLE ERROR";E
9009   I9=0
9010   X8=0
9011   X9=0
9012 GOSUB 9052
9013   F0=H8
9014 H9=(PI-X9)/4
9015 REM H9 IS THE STEP SIZE, 1/4 RANGE INITIALLT
9016   L9=0
9017   X8=X9+H9
9018 GOSUB 9052
9019   F9=F0+4*H8
9020   X8=X8+H9
9021 GOSUB 9052
9022 T9=H8
9023   F9=F9+2*T9
9024   X8=X8+H9
9025 GOSUB 9052
9026   F9=F9+4*H8
9027   X8=X8+H9
9028 GOSUB 9052
9029   F4=H8
9030   F9=F9+F4
9031 REM F9=F0+4*F1+2*F2+4*F3+F4
9032   T9=8*T9+2*(F0+F4)-F9
9033 REM T9 IS THE TRUNCATION ERROR=F0-4*F1+6*F2-4*F3+F4
9034 IF ABS(T9)<E THEN 9039
9035   H9=H9/2
9036   L9=1
9037 REM L9=0 ONLY WHEN THE UPPER LIMIT OF THE INTEGRAL =3.14159265
9038 GOTO 9017
9039   X9=X8
9040   F0=F4
9041   I9=H9*(F9-T9/15)+I9
9042 REM I9 IS ACCUMATLATING THE SUM OF THE INDIVIDUAL INTEGRALS
9043 IF L9=0 THEN 9050
9044 IF ABS(T9)>E/32 THEN 9047
9045 REM TEST TO SEE IF STEP LENGTH CAN BE DOUBLED
9046   H9=H9*2
9047 IF SGN(3.14159-X9-4*H9) <> SGN(H9) THEN 9014
9048 REM TEST FOR UPPER LIMIT
9049 GOTO 9017
9050   J=I9/3
9051 GOTO 9066
9052 REM SUBROUTINE BESSEL INTEGRAND
9053   H8=COS(Z*SIN(X8)-N*X8)/3.14159
9054 IF N=INT(N) THEN 9065
9055 IF Z <> 0 THEN 9058
9056   H8=0
9057 GOTO 9065
9058 IF X8=0 THEN 9065
9059   H7=1/X8-.31831
9060 REM THE SECOND INTEGRAL HAS BEEN TRANSFORMED SO THAT THE LIMITS
9061 REM ARE ZERO TO PI. THE TRANSFORMATION IS X8=1/(T+1/PI)
9062   T8=EXP(X8)
9063 DEF FNS(X8)=(T8-1/T8)/2
9064   H8=H8-SIN(N*3.14159)/(X8*X8*3.14159)*EXP(-Z*FNS(H7)-N*H7)
9065 RETURN
9066 PRINT "N="N;"="Z;"J="J
9499 GO TO 9800
9500 PRINT "YOU BLEW IT ! TRY AGAIN":RESUME 9003
9800 INPUT"AGAIN";R$
9810 IF R$="YES" THEN 9003
9820 IF R$<>"NO" THEN PRINT"YES OR NO PLEASE":GO TO 9800
9999 END