1000'  NAME--BESSEL
1010'
1020'  DESCRIPTION--DEFINES A FUNCTION FOR CALCULATING BESSEL
1030'  FUNCTIONS OF THE FIRST KIND (J) OF ANY ORDER FOR ANY REAL
1040'  ARGUMENT. 
1050'
1060'  SOURCE--INTEGRATION ROUTINE BASED ON SIMPSON'S RULE AND
1070'  INTEGRATES THE FUNCTION GIVEN IN 'HANDBOOK OF MATHEMATICAL
1080'  FUNCTIONS',N.B.S. APPLIED MATH SERIES #55,SECTION 9.1.22.
1090'
1100'  INSTRUCTIONS--MAKES USE OF THE TWO INTERNAL FUNCTIONS FNF AND
1110'  FNS. THE FIRST FUNCTION DEFINES THE INTEGRAND, THE SECOND THE
1120'  HYPERBOLIC SIN. THE ROUTINE MAKES USE OF THE
1130'  FOLLOWING VARIABLES:
1140'  E,H7,H8,H9,I9,L9,F0,F4,F9,T8,T9,X8, AND X9.
1150'  TO USE THE FUNCTION THE CALLING PROGRAM SHOULD HAVE THE 
1160'  ORDER IN N AND THE ARGUMENT IN Z, THEN CALL
1170'  FNB WITH AN ARGUMENT EQUAL TO THE ACCEPTABLE ERROR, FOR EXAMPLE:
1180'   
1190'  80 LET N=2
1200'  90 LET Z=2.5
1210'  100 LET P=FNB(.000001)
1220'
1230'  THIS WOULD ASSIGN TO P THE BESSEL FUNCTION OF ORDER 2 WITH
1240'  ARGUMENT 2.5 WITH AN ERROR OF .000001 AT MOST.
1250'
1260'
1270'  *  *  *  *  *  *   MAIN PROGRAM   *  *  *  *  *  *  *  *  *  
1280'
1290 DEF FNB(E)
1300 LET I9=0
1310 LET X8=X9=0
1320 LET F0=FNF(X8)
1330 LET H9=(3.14159265-X9)/4
1340 REM H9 IS THE STEP SIZE, 1/4 RANGE INITIALLY
1350 LET L9=0
1360 LET X8=X9+H9
1370 LET F9=F0+4*FNF(X8)
1380 LET X8=X8+H9
1390 LET T9=FNF(X8)
1400 LET F9=F9+2*T9
1410 LET X8=X8+H9
1420 LET F9=F9+4*FNF(X8)
1430 LET X8=X8+H9
1440 LET F4=FNF(X8)
1450 LET F9=F9+F4
1460 REM F9 =FO+4*F1+2*F2+4*F3+F4
1470 LET T9=8*T9+2*(F0+F4)-F9
1480 REM T9 IS THE TRUNCATION ERROR=F0-4*F1+6*F2-4*F3+F4
1490 IF ABS(T9)<E THEN 1540
1500 LET H9=H9/2
1510 LET L9=1
1520 REM L9=0 ONLY WHEN THE UPPER LIMIT OF THE INTEGRAL =3.14159265
1530 GOTO 1360
1540 LET X9=X8
1550 LET  F0=F4
1560 LET I9=H9*(F9-T9/15)+I9
1570 REM I9 IS ACCUMALATING THE SUM OF THE INDIVIDUAL INTEGRALS
1580 IF L9 = 0 THEN 1650
1590 IF ABS(T9)>E/32 THEN 1620
1600 REM TEST TO SEE IF STEP LENGTH CAN BE DOUBLED
1610 LET H9=H9*2
1620 IF SGN(3.14159265-X9-4*H9)<>SGN(H9) THEN 1330
1630 REM TEST FOR UPPER LIMIT
1640 GOTO 1360
1650 LET FNB=I9/3
1660 FNEND
1670 DEF FNF(X8)
1680 LET H8=COS(Z*SIN(X8)-N*X8)/3.14159265
1690 IF N=INT(N) THEN 1780
1700 IF Z<>0 THEN 1730
1710 LET H8=0
1720 GOTO 1780
1730 IF X8=0 THEN 1780
1740 LET H7=1/X8-0.318309886
1750 REM THE SECOND INTEGRAL HAS BEEN TRANSFORMED SO THAT THE LIMITS
1760 REM ARE ZERO TO PI. THE TRANSFORMATION IS  X8=1/(T+1/PI)
1770 LET H8=H8-SIN(N*3.14159265)/(X8*X8*3.14159265)*EXP(-Z*FNS(H7)-N*H7)
1780 LET FNF=H8
1790 FNEND
1800 DEF FNS(X8)
1810 LET T8=-EXP(-X8)
1820 LET FNS=(T8-1/T8)/2
1830 FNEND
1840 END