C ALGORITHM 639, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 12, NO. 1, PP. 24-25. C Remark: VOL. 28, NO. 3. PP. 285-300 #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Doc/ # Doc/MakefileF90.Dp # Fortran90/ # Fortran90/Dp/ # Fortran90/Dp/Drivers/ # Fortran90/Dp/Drivers/data1 # Fortran90/Dp/Drivers/data10 # Fortran90/Dp/Drivers/data11 # Fortran90/Dp/Drivers/data12 # Fortran90/Dp/Drivers/data13 # Fortran90/Dp/Drivers/data14 # Fortran90/Dp/Drivers/data15 # Fortran90/Dp/Drivers/data16 # Fortran90/Dp/Drivers/data17 # Fortran90/Dp/Drivers/data18 # Fortran90/Dp/Drivers/data19 # Fortran90/Dp/Drivers/data2 # Fortran90/Dp/Drivers/data20 # Fortran90/Dp/Drivers/data3 # Fortran90/Dp/Drivers/data4 # Fortran90/Dp/Drivers/data5 # Fortran90/Dp/Drivers/data6 # Fortran90/Dp/Drivers/data7 # Fortran90/Dp/Drivers/data8 # Fortran90/Dp/Drivers/data9 # Fortran90/Dp/Drivers/driver.f90 # Fortran90/Dp/Drivers/res1 # Fortran90/Dp/Drivers/res10 # Fortran90/Dp/Drivers/res11 # Fortran90/Dp/Drivers/res12 # Fortran90/Dp/Drivers/res13 # Fortran90/Dp/Drivers/res14 # Fortran90/Dp/Drivers/res15 # Fortran90/Dp/Drivers/res16 # Fortran90/Dp/Drivers/res17 # Fortran90/Dp/Drivers/res18 # Fortran90/Dp/Drivers/res19 # Fortran90/Dp/Drivers/res2 # Fortran90/Dp/Drivers/res20 # Fortran90/Dp/Drivers/res3 # Fortran90/Dp/Drivers/res4 # Fortran90/Dp/Drivers/res5 # Fortran90/Dp/Drivers/res6 # Fortran90/Dp/Drivers/res7 # Fortran90/Dp/Drivers/res8 # Fortran90/Dp/Drivers/res9 # Fortran90/Dp/Drivers/src_715.f # Fortran90/Dp/Src/ # Fortran90/Dp/Src/oscsub.f90 # Fortran90/Dp/Src/prec.f90 # Fortran90/Dp/Src/src.f90 # Original/ # Original/639 # This archive created: Fri May 10 10:58:08 2002 export PATH; PATH=/bin:$PATH if test ! -d 'Doc' then mkdir 'Doc' fi cd 'Doc' if test -f 'MakefileF90.Dp' then echo shar: will not over-write existing file "'MakefileF90.Dp'" else cat << "SHAR_EOF" > 'MakefileF90.Dp' # Fortran 90 compiler and linker F90 = epcf90 F90OPTS = -sloppy -C -d5 -g -temp=/tmp -u F90 = f90 F90OPTS = -xprofile=tcov F90 = f90 F90OPTS = -C F90LINK = $(F90) F90LINKOPTS = $(F90OPTS) all: res Objs1= prec.o oscsub.o src.o src_715.o driver.o driver1: $(Objs1) $(F90LINK) $(F90LINKOPTS) -o driver1 $(Objs1) res: driver1 data1 data2 data3 data4 data5 data6 data7 data8 data9 data10 data11 data12 data13 data14 data15 data16 data17 data18 data19 data20 driver1 < data1 >res1 driver1 < data2 >res2 driver1 < data3 >res3 driver1 < data4 >res4 driver1 < data5 >res5 driver1 < data6 >res6 driver1 < data7 >res7 driver1 < data8 >res8 driver1 < data9 >res9 driver1 < data10 >res10 driver1 < data11 >res11 driver1 < data12 >res12 driver1 < data13 >res13 driver1 < data14 >res14 driver1 < data15 >res15 driver1 < data16 >res16 driver1 < data17 >res17 driver1 < data18 >res18 driver1 < data19 >res19 driver1 < data20 >res20 port.o: port.f $(F90) $(F90OPTS) -c port.f src.o: src.f90 prec.o $(F90) $(F90OPTS) -c src.f90 driver.o: driver.f90 src.o oscsub.o prec.o $(F90) $(F90OPTS) -c driver.f90 oscsub.o: oscsub.f90 prec.o $(F90) $(F90OPTS) -c oscsub.f90 src_715.o: src_715.f $(F90) $(F90OPTS) -c src_715.f prec.o: prec.f90 $(F90) $(F90OPTS) -c prec.f90 SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Fortran90' then mkdir 'Fortran90' fi cd 'Fortran90' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test -f 'data1' then echo shar: will not over-write existing file "'data1'" else cat << "SHAR_EOF" > 'data1' 1 1.0D0 SIN(X)/X Test 1: Successful execution of function 1. 0.0 0.0 1.0e-13 6.2831853071795864769 100 14 -9 0 0 0 F SHAR_EOF fi # end of overwriting check if test -f 'data10' then echo shar: will not over-write existing file "'data10'" else cat << "SHAR_EOF" > 'data10' 4 1.0D0 CONST(1 - EXP(X)) * J0(X) / X . (QUADPACK, PAGE 118) Test 10: use NDIM2 = 11 to illustrate variant behaviour 0.0 0.0 1.0e-13 6.2831853071795864769 100 11 -9 0 40 0 T The previously chosen element (16,14) is unreachable. The equally good, but marginally more expensive element (25,11) is chosen instead. SHAR_EOF fi # end of overwriting check if test -f 'data11' then echo shar: will not over-write existing file "'data11'" else cat << "SHAR_EOF" > 'data11' 1 1.0D0 SIN(X)/X Test 11: Illegal period provided. It should be 2*pi. 0.0 0.0 1.0e-13 -3.1415926535897932384 100 14 -9 0 0 0 F SHAR_EOF fi # end of overwriting check if test -f 'data12' then echo shar: will not over-write existing file "'data12'" else cat << "SHAR_EOF" > 'data12' 2 COS(X) 1.0D0/(1+X**2) Test 12: Incorrect period defined. It should be 2*pi. 0.0 0.0 1.0e-13 3.1415926535897932384 100 14 64 0 40 0 F SHAR_EOF fi # end of overwriting check if test -f 'data13' then echo shar: will not over-write existing file "'data13'" else cat << "SHAR_EOF" > 'data13' 3 1.0D0 EXP(-X*SINH(PI/2)*J1(X) Test 13: Incorrect period provided -- not detected. It should be 2*pi. 0.0 0.0 1.0e-10 3.1415926535897932384 100 14 -9 0 40 0 T This easy calculation ended during the sign change grace period. (You can't lose them all.) SHAR_EOF fi # end of overwriting check if test -f 'data14' then echo shar: will not over-write existing file "'data14'" else cat << "SHAR_EOF" > 'data14' 4 1.0D0 CONST*(1 - EXP(X)) * J0(X) / X . (QUADPACK, page 118) Test 14: Incorrect period provided. It should be 2*pi. 0.0 0.0 1.0e-13 3.1415926535897932384 100 14 -9 0 40 0 T Since GPER = 1, wrong period was not noticed, but sign change pattern ultimately led to termination after the grace period. SHAR_EOF fi # end of overwriting check if test -f 'data15' then echo shar: will not over-write existing file "'data15'" else cat << "SHAR_EOF" > 'data15' 1 1.0D0 SIN(X)/X Test 15: Illegal input value (NQUAD=0) 0.0 0.0 1.0e-13 6.2831853071795864769 100 14 0 0 0 0 F SHAR_EOF fi # end of overwriting check if test -f 'data16' then echo shar: will not over-write existing file "'data16'" else cat << "SHAR_EOF" > 'data16' 1 1.0D0 SIN(X)/X Test 16: Illegal input value (NQUAD=1) 0.0 0.0 1.0e-13 6.2831853071795864769 100 14 1 0 0 0 F SHAR_EOF fi # end of overwriting check if test -f 'data17' then echo shar: will not over-write existing file "'data17'" else cat << "SHAR_EOF" > 'data17' 1 1.0D0 SIN(X)/X Test 17: Illegal input value (NDIM1=0) 0.0 0.0 1.0e-13 6.2831853071795864769 0 14 -9 0 0 0 F SHAR_EOF fi # end of overwriting check if test -f 'data18' then echo shar: will not over-write existing file "'data18'" else cat << "SHAR_EOF" > 'data18' 1 1.0D0 SIN(X)/X Test 18: Illegal input value (NDIM2=0) 0.0 0.0 1.0e-13 6.2831853071795864769 100 0 -9 0 0 0 F SHAR_EOF fi # end of overwriting check if test -f 'data19' then echo shar: will not over-write existing file "'data19'" else cat << "SHAR_EOF" > 'data19' 1 1.0D0 SIN(X)/X Test 19: Failure return from user supplied GAUSS routine 0.0 0.0 1.0e-13 6.2831853071795864769 100 14 -7 0 0 0 F SHAR_EOF fi # end of overwriting check if test -f 'data2' then echo shar: will not over-write existing file "'data2'" else cat << "SHAR_EOF" > 'data2' 2 COS(X) 1.0D0/(1+X**2) Test 2: Successful execution of function 2. 0.0 0.0 1.0e-13 6.2831853071795864769 100 14 64 0 0 0 F SHAR_EOF fi # end of overwriting check if test -f 'data20' then echo shar: will not over-write existing file "'data20'" else cat << "SHAR_EOF" > 'data20' 1 1.0D0 SIN(X)/X Test 20: Use value of RFIRST other than zero 0.0 3.1415926535897932384 1.0e-13 6.2831853071795864769 100 14 -9 0 0 0 F SHAR_EOF fi # end of overwriting check if test -f 'data3' then echo shar: will not over-write existing file "'data3'" else cat << "SHAR_EOF" > 'data3' 3 1.0D0 EXP(-X*SINH(PI/2)*J1(X) Test 3: Successful execution of function 3. 0.0 0.0 1.0e-13 6.2831853071795864769 100 14 -9 0 0 0 F SHAR_EOF fi # end of overwriting check if test -f 'data4' then echo shar: will not over-write existing file "'data4'" else cat << "SHAR_EOF" > 'data4' 4 1.0D0 CONST(1 - EXP(X)) * J0(X) / X . (QUADPACK, PAGE 118) Test 4: Successful execution of function 4. 0.0 0.0 1.0e-13 6.2831853071795864769 100 14 -9 0 0 0 F SHAR_EOF fi # end of overwriting check if test -f 'data5' then echo shar: will not over-write existing file "'data5'" else cat << "SHAR_EOF" > 'data5' 5 1.0D0 (1 + 100*EXP(-X**2/PI**2))*J1(X) Test 5: Successful execution of function 5. 0.0 0.0 1.0e-13 6.2831853071795864769 100 14 -9 0 0 0 F SHAR_EOF fi # end of overwriting check if test -f 'data6' then echo shar: will not over-write existing file "'data6'" else cat << "SHAR_EOF" > 'data6' 6 (COS**4(X) - 3/8) 1.0D0/(1 + X**2) Test 6: Successful execution of function 6. 0.0 0.0 1.0e-11 3.1415926535897932384 200 14 -9 0 0 0 F SHAR_EOF fi # end of overwriting check if test -f 'data7' then echo shar: will not over-write existing file "'data7'" else cat << "SHAR_EOF" > 'data7' 1 1.0D0 SIN(X)/X Test 7: Successful execution of function 1 (extra output). 0.0 0.0 1.0e-13 6.2831853071795864769 100 14 -9 12 40 12 F SHAR_EOF fi # end of overwriting check if test -f 'data8' then echo shar: will not over-write existing file "'data8'" else cat << "SHAR_EOF" > 'data8' 2 COS(X) 1.0D0/(1+X**2) Test 8: Function 2 with NQUAD = -9. Generates incorrect result. 0.0 0.0 1.0e-13 6.2831853071795864769 100 14 -9 0 40 0 T This quadrature rule is inadequate. Note the steady convergence of the table elements below to a result which is in fact wrong. Routine cannot know this. SHAR_EOF fi # end of overwriting check if test -f 'data9' then echo shar: will not over-write existing file "'data9'" else cat << "SHAR_EOF" > 'data9' 3 1.0D0 EXP(-X*SINH(PI/2)*J1(X) Test 9: Set NDIM = 5 to illustrate non convergence. 0.0 0.0 1.0e-13 6.2831853071795864769 5 14 -9 0 0 0 F SHAR_EOF fi # end of overwriting check if test -f 'driver.f90' then echo shar: will not over-write existing file "'driver.f90'" else cat << "SHAR_EOF" > 'driver.f90' MODULE constants_639 USE prec, ONLY : wp REAL (wp), PARAMETER :: one = 1.0_wp, two = 2.0_wp, zero = 0.0_wp, & four = 4.0_wp, sixteen = 16.0_wp, half = 0.5_wp, & three_eights = 0.375_wp, six = 6.0_wp, one_twenty = 120.0_wp, & ten_to_minus_4 = 10.0E-4_wp, ten_to_minus_8 = 10.0E-8_wp, & quarter = 0.25_wp, hundred = 100.0_wp END MODULE constants_639 PROGRAM driver_639 ! .. Use Statements .. USE calgo_639, ONLY : oscint USE prec, ONLY : wp USE oscsub, ONLY : alpha, const6, ebase, halfpi, nb, nix, numfun, pi USE constants_639, ONLY : four, half, one, two, zero ! .. ! .. Local Scalars .. REAL (wp) :: azero, cin, eps, period, result, rfirst INTEGER :: alloc, endv, i, jp, jpmk, jpr, kpr, ndim1, ndim2, nquad LOGICAL :: comments CHARACTER (50) :: form CHARACTER (80) :: gname CHARACTER (80) :: hname, notes, error_mess ! .. ! .. Local Arrays .. REAL (wp), ALLOCATABLE :: work(:,:) INTEGER :: istate(6) CHARACTER (80) :: com(3) ! .. ! .. External Functions .. REAL (wp), EXTERNAL :: exact, gper, hfun ! .. ! .. External Subroutines .. EXTERNAL g5and9 ! .. ! .. Intrinsic Functions .. INTRINSIC atan, exp, log, sqrt, min ! CONSTANTS, ETC.USED IN SUBROUTINES. ebase = exp(one) pi = four*atan(one) halfpi = pi*half cin = one + sqrt(two) const6 = one/log(cin) alpha = zero nb = 2 10 READ (*,*,end=20) numfun READ (*,'(a)') gname READ (*,'(a)') hname READ (*,'(a)') notes READ (*,*) azero, rfirst, eps, period READ (*,*) ndim1, ndim2, nquad READ (*,*) nix(1:3) READ (*,*) comments IF (comments) THEN READ (*,'(a)') com(1:3) END IF ALLOCATE (work(ndim1,ndim2),STAT=alloc) IF (alloc/=0) THEN WRITE (*,'("Allocation fails in driver_639")') STOP END IF WRITE (*,90000) numfun, gname, hname WRITE (*,90040) notes WRITE (*,90050) azero, period, rfirst, eps, nquad WRITE (*,90060) ndim1, ndim2 WRITE (*,90070) IF (nix(1)>0 .OR. nix(3)>0) WRITE (*,90080) CALL oscint(azero,period,rfirst,eps,nquad,work,result,istate,g5and9, & hfun,gper) WRITE (*,90090) IF (istate(1)>=0) THEN WRITE (*,90110) result, exact(), exact() - result END IF WRITE (*,90120) (istate(i),i=1,5) WRITE (*,90140) istate(6) WRITE (*,90020) IF (istate(1)>=0) THEN WRITE (*,90100) ELSE WRITE (*,'(1x,''ERROR: '',a)') error_mess(istate(1)) END IF IF (comments) THEN WRITE (*,90130) com(1), com(2), com(3) END IF WRITE (*,90030) IF (nix(2)>0) THEN WRITE (*,90010) jp = istate(3) jpmk = jp - istate(4) DO jpr = 1, min(istate(3),nix(2)) endv = jp - jpr IF (jpr>jpmk) THEN endv = endv + 1 END IF endv = min(endv,ndim2) WRITE (form,fmt='("(1x,i3,",i3,"(7d10.2,/4X))")') endv WRITE (*,fmt=form) jpr, (work(jpr,kpr),kpr=1,endv) END DO END IF DEALLOCATE (work) GO TO 10 90000 FORMAT (/' NUMFUN = ',1X,I1,10X,'F(X) = G(X).H(X).'/22X,'G(X) = ',A/, & 22X,'H(X) = ',A) 90010 FORMAT (/25X,'Finite Average Table'/25X,'--------------------') 90020 FORMAT ( & ' ************************** Exit status **************************'/) 90030 FORMAT (/ & ' *****************************************************************') 90040 FORMAT (/' INPUT PARAMETERS-'//1X,A80/) 90050 FORMAT (4X,'AZERO',12X,'PERIOD',12X,'RFIRST',12X,'EPS',12X,'NQUAD', & 12X/D13.6,5X,D12.6,6X,D12.6,6X,D12.6,6X,I3/) 90060 FORMAT (/4X,'NDIM1',10X,'NDIM2'/5X,I3,11X,I3/) 90070 FORMAT (' ----------------------') 90080 FORMAT (/' Code check output from HFUN, GPER, QRULE'/) 90090 FORMAT (/' PRINCIPAL OUTPUT PARAMETERS-'//) 90100 FORMAT (1X,' Successful exit from OSCINT'/) 90110 FORMAT (5X,'RESULT',20X,'Exact',21X,'Difference'/D17.8,10X,D16.8,10X, & D16.8,10X//) 90120 FORMAT (4X,'ISTATE(1)',5X,'ISTATE(2)',5X,'ISTATE(3)',5X,'ISTATE(4)',8X, & 'ISTATE(5)'/5X,'(IFAIL)',7X,'(LSIGCH)',8X,'(N)',7X,'(Column No.)',6X, & '(Row No.)'/I9,11X,I3,11X,I3,11X,I3,13X,I4/) 90130 FORMAT (3(/1X,A)) 90140 FORMAT (4X,' No. of function values (ISTATE(6)) = ',I6/) 20 END PROGRAM driver_639 CHARACTER (len=80) FUNCTION error_mess(ifail) USE calgo_639, ONLY : err_gauss_routine, err_max_intervals, & err_sign_change, err_input_period, err_gper_function, err_input_params INTEGER, INTENT (IN) :: ifail SELECT CASE (ifail) CASE (err_gauss_routine) error_mess = 'Error reported from user supplied function GAUSS' CASE (err_max_intervals) error_mess = 'Maximum number of intervals exceeded' CASE (err_sign_change) error_mess = & 'Normal sign change pattern is violated after the grace period' CASE (err_input_period) error_mess = 'Supplied value of PERIOD <= 10**(-5)' CASE (err_gper_function) error_mess = & 'Routine has detected that GPER is not periodic with given period' CASE (err_input_params) error_mess = 'Input parameter error (NQUAD = 0,1 or WORK & &array not at least (1,1))' CASE DEFAULT error_mess = 'ILLEGAL ERROR RETURN -- PLEASE REPORT' END SELECT END FUNCTION error_mess FUNCTION exact() USE prec, ONLY : wp USE oscsub, ONLY : alpha, ebase, halfpi, numfun, pi USE constants_639, ONLY : four, hundred, one, quarter, sixteen, two, & zero ! EXACT IS THE VALUE OF THE INTEGRAL OF F(X) OVER (0,INFINITY). ! HERE, OF COURSE, F(X) = H(X).G(X) GIVEN BY HFUN AND GPER BELOW. ! .. Function Return Value .. REAL (wp) :: exact ! .. ! .. Intrinsic Functions .. INTRINSIC cosh, exp SELECT CASE (numfun) CASE (1) exact = halfpi CASE (2) exact = pi/(two*ebase) CASE (3) exact = exp(-(one+alpha)*halfpi)/cosh(halfpi) CASE (4) exact = one CASE (5) exact = one + hundred*(one-exp(-pi*pi*quarter)) CASE (6) exact = pi/sixteen*(ebase**(-4)+four*ebase**(-2)) CASE DEFAULT exact = zero END SELECT END FUNCTION exact FUNCTION hfun(x) USE prec, ONLY : wp USE oscsub, ONLY : alpha, const6, ebase, halfpi, nb, nix, numfun, pi USE constants_639, ONLY : hundred, one, one_twenty, six, ten_to_minus_4, & ten_to_minus_8, two, zero ! .. Function Return Value .. REAL (wp) :: hfun ! .. ! .. Scalar Arguments .. REAL (wp) :: x ! .. ! .. Local Scalars .. INTEGER :: ncalc ! .. ! .. Local Arrays .. REAL (wp) :: b(nb) ! .. ! .. External Subroutines .. EXTERNAL rjbesl ! .. ! .. Intrinsic Functions .. INTRINSIC abs, exp, sin ! ALPHA - FRACTIONAL PART OF ORDER FOR WHICH BESSEL FUNCTION IS TO ! BE CALCULATED, 0 <= ALPHA < 1 ! NB NUMBER OF BESSEL FUNCTIONS IN SEQUENCE. ! IF (NUMFUN .EQ. 1) F(X) = SIN(X)/X ! IF (NUMFUN .EQ. 2) F(X) = COS(X)/(1+X**2) ! IF (NUMFUN .EQ. 3) F(X) = EXP(-X*SINH(PI/2)*J1(X) ! IF (NUMFUN .EQ. 4) F(X) = QUADPACK 5.2.3. PAGE118 ! IF (NUMFUN .EQ. 5) F(X) = (1 + 100*EXP(-X**2/PI**2))*J1(X) ! IF (NUMFUN .EQ. 6) F(X) = (COS**4(X) - 3/8)/(1 + X**2) ! WHEN NUMFUN IS 1,3,4 OR 5, H(X) = F(X). ! WHEN NUMFUN IS 2 OR 6, H(X) = 1/(1 + X**2) SELECT CASE (numfun) CASE (1) IF (abs(x)>ten_to_minus_4) THEN hfun = sin(x)/x ELSE hfun = one - x**2/six + x**4/one_twenty END IF CASE (2) hfun = one/(one+(x**two)) CASE (3) CALL rjbesl(x,alpha,nb,b,ncalc) hfun = ebase**(-x*sinh(halfpi))*b(2) CASE (4) CALL rjbesl(x,alpha,nb,b,ncalc) hfun = const6*b(1) IF (x>=ten_to_minus_8) THEN hfun = ((one-exp(-x))/x)*hfun END IF CASE (5) CALL rjbesl(x,alpha,nb,b,ncalc) hfun = (one+hundred*exp(-(x**two)/(pi**two)))*b(2) CASE (6) hfun = one/(one+(x**two)) CASE DEFAULT hfun = zero END SELECT IF (nix(3)>0) THEN nix(3) = nix(3) - 1 WRITE (*,'('' X = '',D16.8,'' HFUN = '',D16.8)') x, hfun END IF END FUNCTION hfun FUNCTION gper(x) USE prec, ONLY : wp USE oscsub, ONLY : nix, numfun USE constants_639, ONLY : four, one, three_eights, zero ! .. Function Return Value .. REAL (wp) :: gper ! .. ! .. Scalar Arguments .. REAL (wp) :: x SELECT CASE (numfun) CASE (1,3,4,5) gper = one CASE (2) gper = cos(x) CASE (6) gper = (cos(x))**four - three_eights CASE DEFAULT gper = zero END SELECT IF (nix(3)/=0) THEN WRITE (*,'('' X = '',D16.8,'' GPER = '',D16.8)') x, gper END IF RETURN END FUNCTION gper SUBROUTINE g5and9(nquad,weight,abscis,ierr) USE prec, ONLY : wp ! THIS IS AN EXTRACT FROM A QUADRATURE ROUTINE CONSTRUCTED ONLY FOR ! USE IN A DRIVER WHICH ILLUSTRATES OSCINT. A NAG LIBRARY ! SUBSCRIBER MAY REPLACE THIS BY D01BCF FOR GENERAL USE. ! THE FOLLOWING DATA ARE WEIGHTS AND ABSCISSAS OF THE FIVE POINT ! AND THE NINE POINT GAUSS-LEGENDRE QUADRATURE RULES RESPECTIVELY. ! NORMALISED TO THE INTERVAL (-1,1). ! .. Scalar Arguments .. INTEGER :: ierr, nquad ! .. ! .. Array Arguments .. REAL (wp) :: abscis(nquad), weight(nquad) ! .. ! .. Local Arrays .. REAL (wp) :: absc5(5), absc9(9), wt5(5), wt9(9) ! .. ! .. Data Statements .. DATA wt5/.2369268850561891_wp, .4786286704993665_wp, & .5688888888888889_wp, .4786286704993665_wp, .2369268850561891_wp/ DATA absc5/ -.9061798459386640_wp, -.5384693101056831_wp, 0.0_wp, & .5384693101056830_wp, .9061798459386640_wp/ DATA wt9/.0812743883615745_wp, .1806481606948574_wp, & .2606106964029355_wp, .3123470770400029_wp, .3302393550012598_wp, & .3123470770400029_wp, .2606106964029356_wp, .1806481606948575_wp, & .8127438836157467D-01/ DATA absc9/ -.9681602395076261_wp, -.8360311073266358_wp, & -.6133714327005904_wp, -.3242534234038090_wp, 0.0_wp, & .3242534234038087_wp, .6133714327005902_wp, .8360311073266357_wp, & .9681602395076260_wp/ ! .. SELECT CASE (nquad) CASE (5) abscis(1:5) = absc5 weight(1:5) = wt5 ierr = 0 CASE (9) abscis(1:9) = absc9 weight(1:9) = wt9 ierr = 0 CASE DEFAULT ierr = 59 END SELECT END SUBROUTINE g5and9 SHAR_EOF fi # end of overwriting check if test -f 'res1' then echo shar: will not over-write existing file "'res1'" else cat << "SHAR_EOF" > 'res1' NUMFUN = 1 F(X) = G(X).H(X). G(X) = 1.0D0 H(X) = SIN(X)/X INPUT PARAMETERS- Test 1: Successful execution of function 1. AZERO PERIOD RFIRST EPS NQUAD 0.000000D+00 0.628319D+01 0.000000D+00 0.100000D-12 -9 NDIM1 NDIM2 100 14 ---------------------- PRINCIPAL OUTPUT PARAMETERS- RESULT Exact Difference 0.15707963D+01 0.15707963D+01 0.99920072D-14 ISTATE(1) ISTATE(2) ISTATE(3) ISTATE(4) ISTATE(5) (IFAIL) (LSIGCH) (N) (Column No.) (Row No.) 0 0 31 14 18 No. of function values (ISTATE(6)) = 279 ************************** Exit status ************************** Successful exit from OSCINT ***************************************************************** SHAR_EOF fi # end of overwriting check if test -f 'res10' then echo shar: will not over-write existing file "'res10'" else cat << "SHAR_EOF" > 'res10' NUMFUN = 4 F(X) = G(X).H(X). G(X) = 1.0D0 H(X) = CONST(1 - EXP(X)) * J0(X) / X . (QUADPACK, PAGE 118) INPUT PARAMETERS- Test 10: use NDIM2 = 11 to illustrate variant behaviour AZERO PERIOD RFIRST EPS NQUAD 0.000000D+00 0.628319D+01 0.000000D+00 0.100000D-12 -9 NDIM1 NDIM2 100 11 ---------------------- PRINCIPAL OUTPUT PARAMETERS- RESULT Exact Difference 0.10000000D+01 0.10000000D+01 -0.18207658D-13 ISTATE(1) ISTATE(2) ISTATE(3) ISTATE(4) ISTATE(5) (IFAIL) (LSIGCH) (N) (Column No.) (Row No.) 0 0 35 11 25 No. of function values (ISTATE(6)) = 315 ************************** Exit status ************************** Successful exit from OSCINT The previously chosen element (16,14) is unreachable. The equally good, but marginally more expensive element (25,11) is chosen instead. ***************************************************************** Finite Average Table -------------------- 1 0.11D+01 0.48D+00 0.21D+00 0.98D-01 0.46D-01 0.22D-01 0.10D-01 0.49D-02 0.23D-02 0.11D-02 0.54D-03 2 -0.17D+00 -0.50D-01 -0.17D-01 -0.66D-02 -0.26D-02 -0.11D-02 -0.46D-03 -0.20D-03 -0.87D-04 -0.38D-04 -0.17D-04 3 0.71D-01 0.15D-01 0.42D-02 0.13D-02 0.46D-03 0.17D-03 0.63D-04 0.25D-04 0.10D-04 0.41D-05 0.17D-05 4 -0.41D-01 -0.68D-02 -0.15D-02 -0.42D-03 -0.12D-03 -0.40D-04 -0.14D-04 -0.49D-05 -0.18D-05 -0.68D-06 -0.26D-06 5 0.27D-01 0.37D-02 0.71D-03 0.17D-03 0.44D-04 0.13D-04 0.40D-05 0.13D-05 0.44D-06 0.15D-06 0.56D-07 6 -0.20D-01 -0.23D-02 -0.38D-03 -0.78D-04 -0.19D-04 -0.49D-05 -0.14D-05 -0.41D-06 -0.13D-06 -0.43D-07 -0.14D-07 7 0.15D-01 0.15D-02 0.22D-03 0.41D-04 0.89D-05 0.21D-05 0.55D-06 0.15D-06 0.45D-07 0.14D-07 0.44D-08 8 -0.12D-01 -0.11D-02 -0.14D-03 -0.24D-04 -0.46D-05 -0.10D-05 -0.25D-06 -0.63D-07 -0.17D-07 -0.50D-08 -0.15D-08 9 0.99D-02 0.78D-03 0.94D-04 0.14D-04 0.26D-05 0.53D-06 0.12D-06 0.29D-07 0.74D-08 0.20D-08 0.57D-09 10 -0.83D-02 -0.60D-03 -0.65D-04 -0.91D-05 -0.15D-05 -0.29D-06 -0.61D-07 -0.14D-07 -0.34D-08 -0.87D-09 -0.23D-09 11 0.71D-02 0.47D-03 0.47D-04 0.61D-05 0.95D-06 0.17D-06 0.33D-07 0.71D-08 0.16D-08 0.40D-09 0.10D-09 12 -0.62D-02 -0.37D-03 -0.35D-04 -0.42D-05 -0.61D-06 -0.10D-06 -0.19D-07 -0.39D-08 -0.84D-09 -0.19D-09 -0.48D-10 13 0.54D-02 0.30D-03 0.26D-04 0.30D-05 0.41D-06 0.64D-07 0.11D-07 0.22D-08 0.45D-09 0.10D-09 0.23D-10 14 -0.48D-02 -0.25D-03 -0.20D-04 -0.22D-05 -0.28D-06 -0.41D-07 -0.69D-08 -0.13D-08 -0.25D-09 -0.53D-10 -0.12D-10 15 0.43D-02 0.21D-03 0.16D-04 0.16D-05 0.20D-06 0.28D-07 0.44D-08 0.77D-09 0.15D-09 0.29D-10 0.63D-11 16 -0.39D-02 -0.18D-03 -0.13D-04 -0.12D-05 -0.14D-06 -0.19D-07 -0.29D-08 -0.48D-09 -0.86D-10 -0.17D-10 -0.35D-11 17 0.35D-02 0.15D-03 0.10D-04 0.93D-06 0.10D-06 0.13D-07 0.19D-08 0.30D-09 0.53D-10 0.98D-11 0.20D-11 18 -0.32D-02 -0.13D-03 -0.85D-05 -0.73D-06 -0.76D-07 -0.93D-08 -0.13D-08 -0.20D-09 -0.33D-10 -0.59D-11 -0.11D-11 19 0.30D-02 0.11D-03 0.70D-05 0.57D-06 0.57D-07 0.67D-08 0.90D-09 0.13D-09 0.21D-10 0.37D-11 0.68D-12 20 -0.27D-02 -0.10D-03 -0.59D-05 -0.46D-06 -0.44D-07 -0.49D-08 -0.63D-09 -0.90D-10 -0.14D-10 -0.23D-11 -0.41D-12 21 0.25D-02 0.89D-04 0.50D-05 0.37D-06 0.34D-07 0.37D-08 0.45D-09 0.62D-10 0.92D-11 0.15D-11 0.26D-12 22 -0.24D-02 -0.79D-04 -0.42D-05 -0.30D-06 -0.27D-07 -0.28D-08 -0.33D-09 -0.43D-10 -0.63D-11 -0.98D-12 -0.16D-12 23 0.22D-02 0.71D-04 0.36D-05 0.25D-06 0.21D-07 0.21D-08 0.24D-09 0.31D-10 0.43D-11 0.65D-12 0.11D-12 24 -0.21D-02 -0.63D-04 -0.31D-05 -0.21D-06 -0.17D-07 -0.16D-08 -0.18D-09 -0.22D-10 -0.30D-11 -0.44D-12 -0.69D-13 25 0.19D-02 0.57D-04 0.27D-05 0.17D-06 0.14D-07 0.13D-08 0.14D-09 0.16D-10 0.21D-11 0.30D-12 0.46D-13 26 -0.18D-02 -0.52D-04 -0.24D-05 -0.15D-06 -0.11D-07 -0.10D-08 -0.10D-09 -0.12D-10 -0.15D-11 -0.21D-12 27 0.17D-02 0.47D-04 0.21D-05 0.12D-06 0.91D-08 0.80D-09 0.80D-10 0.89D-11 0.11D-11 28 -0.16D-02 -0.43D-04 -0.18D-05 -0.10D-06 -0.75D-08 -0.64D-09 -0.62D-10 -0.67D-11 29 0.15D-02 0.39D-04 0.16D-05 0.90D-07 0.62D-08 0.51D-09 0.48D-10 30 -0.15D-02 -0.36D-04 -0.14D-05 -0.77D-07 -0.52D-08 -0.42D-09 31 0.14D-02 0.33D-04 0.13D-05 0.67D-07 0.44D-08 32 -0.13D-02 -0.31D-04 -0.11D-05 -0.58D-07 33 0.13D-02 0.28D-04 0.10D-05 34 -0.12D-02 -0.26D-04 35 0.12D-02 SHAR_EOF fi # end of overwriting check if test -f 'res11' then echo shar: will not over-write existing file "'res11'" else cat << "SHAR_EOF" > 'res11' NUMFUN = 1 F(X) = G(X).H(X). G(X) = 1.0D0 H(X) = SIN(X)/X INPUT PARAMETERS- Test 11: Illegal period provided. It should be 2*pi. AZERO PERIOD RFIRST EPS NQUAD 0.000000D+00 -.314159D+01 0.000000D+00 0.100000D-12 -9 NDIM1 NDIM2 100 14 ---------------------- PRINCIPAL OUTPUT PARAMETERS- ISTATE(1) ISTATE(2) ISTATE(3) ISTATE(4) ISTATE(5) (IFAIL) (LSIGCH) (N) (Column No.) (Row No.) -3000 0 0 0 0 No. of function values (ISTATE(6)) = 0 ************************** Exit status ************************** ERROR: Supplied value of PERIOD <= 10**(-5) ***************************************************************** SHAR_EOF fi # end of overwriting check if test -f 'res12' then echo shar: will not over-write existing file "'res12'" else cat << "SHAR_EOF" > 'res12' NUMFUN = 2 F(X) = G(X).H(X). G(X) = COS(X) H(X) = 1.0D0/(1+X**2) INPUT PARAMETERS- Test 12: Incorrect period defined. It should be 2*pi. AZERO PERIOD RFIRST EPS NQUAD 0.000000D+00 0.314159D+01 0.000000D+00 0.100000D-12 64 NDIM1 NDIM2 100 14 ---------------------- PRINCIPAL OUTPUT PARAMETERS- ISTATE(1) ISTATE(2) ISTATE(3) ISTATE(4) ISTATE(5) (IFAIL) (LSIGCH) (N) (Column No.) (Row No.) -5000 3 3 3 1 No. of function values (ISTATE(6)) = 253 ************************** Exit status ************************** ERROR: Routine has detected that GPER is not periodic with given period ***************************************************************** Finite Average Table -------------------- 1 0.75D+00 0.31D+00 0.10D+00 2 -0.14D+00 -0.10D+00 3 -0.69D-01 SHAR_EOF fi # end of overwriting check if test -f 'res13' then echo shar: will not over-write existing file "'res13'" else cat << "SHAR_EOF" > 'res13' NUMFUN = 3 F(X) = G(X).H(X). G(X) = 1.0D0 H(X) = EXP(-X*SINH(PI/2)*J1(X) INPUT PARAMETERS- Test 13: Incorrect period provided -- not detected. It should be 2*pi. AZERO PERIOD RFIRST EPS NQUAD 0.000000D+00 0.314159D+01 0.000000D+00 0.100000D-09 -9 NDIM1 NDIM2 100 14 ---------------------- PRINCIPAL OUTPUT PARAMETERS- RESULT Exact Difference 0.82847664D-01 0.82847664D-01 -0.38996584D-14 ISTATE(1) ISTATE(2) ISTATE(3) ISTATE(4) ISTATE(5) (IFAIL) (LSIGCH) (N) (Column No.) (Row No.) 3 8 9 2 8 No. of function values (ISTATE(6)) = 81 ************************** Exit status ************************** Successful exit from OSCINT This easy calculation ended during the sign change grace period. (You can't lose them all.) ***************************************************************** Finite Average Table -------------------- 1 0.76D-01 0.41D-01 0.22D-01 0.12D-01 0.64D-02 0.34D-02 0.18D-02 0.95D-03 2 0.63D-02 0.32D-02 0.16D-02 0.81D-03 0.40D-03 0.20D-03 0.10D-03 3 0.39D-04 0.18D-04 0.84D-05 0.39D-05 0.18D-05 0.80D-06 4 -0.26D-05 -0.13D-05 -0.67D-06 -0.34D-06 -0.17D-06 5 -0.23D-07 -0.11D-07 -0.49D-08 -0.23D-08 6 0.15D-08 0.74D-09 0.37D-09 7 0.14D-10 0.66D-11 8 -0.90D-12 -0.46D-12 9 -0.90D-14 SHAR_EOF fi # end of overwriting check if test -f 'res14' then echo shar: will not over-write existing file "'res14'" else cat << "SHAR_EOF" > 'res14' NUMFUN = 4 F(X) = G(X).H(X). G(X) = 1.0D0 H(X) = CONST*(1 - EXP(X)) * J0(X) / X . (QUADPACK, page 118) INPUT PARAMETERS- Test 14: Incorrect period provided. It should be 2*pi. AZERO PERIOD RFIRST EPS NQUAD 0.000000D+00 0.314159D+01 0.000000D+00 0.100000D-12 -9 NDIM1 NDIM2 100 14 ---------------------- PRINCIPAL OUTPUT PARAMETERS- ISTATE(1) ISTATE(2) ISTATE(3) ISTATE(4) ISTATE(5) (IFAIL) (LSIGCH) (N) (Column No.) (Row No.) -200 11 11 6 4 No. of function values (ISTATE(6)) = 108 ************************** Exit status ************************** ERROR: Normal sign change pattern is violated after the grace period Since GPER = 1, wrong period was not noticed, but sign change pattern ultimately led to termination after the grace period. ***************************************************************** Finite Average Table -------------------- 1 0.11D+01 0.56D+00 0.25D+00 0.92D-01 0.22D-01 -0.67D-03 -0.32D-02 -0.62D-03 0.13D-02 0.17D-02 2 0.55D-01 -0.54D-01 -0.70D-01 -0.49D-01 -0.23D-01 -0.57D-02 0.20D-02 0.33D-02 0.20D-02 3 -0.16D+00 -0.86D-01 -0.28D-01 0.26D-02 0.12D-01 0.96D-02 0.46D-02 0.71D-03 4 -0.80D-02 0.30D-01 0.33D-01 0.21D-01 0.77D-02 -0.41D-03 -0.32D-02 5 0.68D-01 0.35D-01 0.86D-02 -0.54D-02 -0.85D-02 -0.60D-02 6 0.26D-02 -0.18D-01 -0.19D-01 -0.12D-01 -0.35D-02 0.13D-02 7 -0.39D-01 -0.20D-01 -0.39D-02 0.46D-02 0.61D-02 8 -0.12D-02 0.13D-01 0.13D-01 0.76D-02 9 0.26D-01 0.13D-01 0.21D-02 10 0.67D-03 -0.92D-02 11 -0.19D-01 SHAR_EOF fi # end of overwriting check if test -f 'res15' then echo shar: will not over-write existing file "'res15'" else cat << "SHAR_EOF" > 'res15' NUMFUN = 1 F(X) = G(X).H(X). G(X) = 1.0D0 H(X) = SIN(X)/X INPUT PARAMETERS- Test 15: Illegal input value (NQUAD=0) AZERO PERIOD RFIRST EPS NQUAD 0.000000D+00 0.628319D+01 0.000000D+00 0.100000D-12 0 NDIM1 NDIM2 100 14 ---------------------- PRINCIPAL OUTPUT PARAMETERS- ISTATE(1) ISTATE(2) ISTATE(3) ISTATE(4) ISTATE(5) (IFAIL) (LSIGCH) (N) (Column No.) (Row No.) -6000 0 0 0 0 No. of function values (ISTATE(6)) = 0 ************************** Exit status ************************** ERROR: Input parameter error (NQUAD = 0,1 or WORK array not at least (1,1)) ***************************************************************** SHAR_EOF fi # end of overwriting check if test -f 'res16' then echo shar: will not over-write existing file "'res16'" else cat << "SHAR_EOF" > 'res16' NUMFUN = 1 F(X) = G(X).H(X). G(X) = 1.0D0 H(X) = SIN(X)/X INPUT PARAMETERS- Test 16: Illegal input value (NQUAD=1) AZERO PERIOD RFIRST EPS NQUAD 0.000000D+00 0.628319D+01 0.000000D+00 0.100000D-12 1 NDIM1 NDIM2 100 14 ---------------------- PRINCIPAL OUTPUT PARAMETERS- ISTATE(1) ISTATE(2) ISTATE(3) ISTATE(4) ISTATE(5) (IFAIL) (LSIGCH) (N) (Column No.) (Row No.) -6000 0 0 0 0 No. of function values (ISTATE(6)) = 0 ************************** Exit status ************************** ERROR: Input parameter error (NQUAD = 0,1 or WORK array not at least (1,1)) ***************************************************************** SHAR_EOF fi # end of overwriting check if test -f 'res17' then echo shar: will not over-write existing file "'res17'" else cat << "SHAR_EOF" > 'res17' NUMFUN = 1 F(X) = G(X).H(X). G(X) = 1.0D0 H(X) = SIN(X)/X INPUT PARAMETERS- Test 17: Illegal input value (NDIM1=0) AZERO PERIOD RFIRST EPS NQUAD 0.000000D+00 0.628319D+01 0.000000D+00 0.100000D-12 -9 NDIM1 NDIM2 0 14 ---------------------- PRINCIPAL OUTPUT PARAMETERS- ISTATE(1) ISTATE(2) ISTATE(3) ISTATE(4) ISTATE(5) (IFAIL) (LSIGCH) (N) (Column No.) (Row No.) -6000 0 0 0 0 No. of function values (ISTATE(6)) = 0 ************************** Exit status ************************** ERROR: Input parameter error (NQUAD = 0,1 or WORK array not at least (1,1)) ***************************************************************** SHAR_EOF fi # end of overwriting check if test -f 'res18' then echo shar: will not over-write existing file "'res18'" else cat << "SHAR_EOF" > 'res18' NUMFUN = 1 F(X) = G(X).H(X). G(X) = 1.0D0 H(X) = SIN(X)/X INPUT PARAMETERS- Test 18: Illegal input value (NDIM2=0) AZERO PERIOD RFIRST EPS NQUAD 0.000000D+00 0.628319D+01 0.000000D+00 0.100000D-12 -9 NDIM1 NDIM2 100 0 ---------------------- PRINCIPAL OUTPUT PARAMETERS- ISTATE(1) ISTATE(2) ISTATE(3) ISTATE(4) ISTATE(5) (IFAIL) (LSIGCH) (N) (Column No.) (Row No.) -6000 0 0 0 0 No. of function values (ISTATE(6)) = 0 ************************** Exit status ************************** ERROR: Input parameter error (NQUAD = 0,1 or WORK array not at least (1,1)) ***************************************************************** SHAR_EOF fi # end of overwriting check if test -f 'res19' then echo shar: will not over-write existing file "'res19'" else cat << "SHAR_EOF" > 'res19' NUMFUN = 1 F(X) = G(X).H(X). G(X) = 1.0D0 H(X) = SIN(X)/X INPUT PARAMETERS- Test 19: Failure return from user supplied GAUSS routine AZERO PERIOD RFIRST EPS NQUAD 0.000000D+00 0.628319D+01 0.000000D+00 0.100000D-12 -7 NDIM1 NDIM2 100 14 ---------------------- PRINCIPAL OUTPUT PARAMETERS- ISTATE(1) ISTATE(2) ISTATE(3) ISTATE(4) ISTATE(5) (IFAIL) (LSIGCH) (N) (Column No.) (Row No.) -4000 0 0 0 0 No. of function values (ISTATE(6)) = 0 ************************** Exit status ************************** ERROR: Error reported from user supplied function GAUSS ***************************************************************** SHAR_EOF fi # end of overwriting check if test -f 'res2' then echo shar: will not over-write existing file "'res2'" else cat << "SHAR_EOF" > 'res2' NUMFUN = 2 F(X) = G(X).H(X). G(X) = COS(X) H(X) = 1.0D0/(1+X**2) INPUT PARAMETERS- Test 2: Successful execution of function 2. AZERO PERIOD RFIRST EPS NQUAD 0.000000D+00 0.628319D+01 0.000000D+00 0.100000D-12 64 NDIM1 NDIM2 100 14 ---------------------- PRINCIPAL OUTPUT PARAMETERS- RESULT Exact Difference 0.57786367D+00 0.57786367D+00 -0.12878587D-13 ISTATE(1) ISTATE(2) ISTATE(3) ISTATE(4) ISTATE(5) (IFAIL) (LSIGCH) (N) (Column No.) (Row No.) 0 0 27 13 15 No. of function values (ISTATE(6)) = 1702 ************************** Exit status ************************** Successful exit from OSCINT ***************************************************************** SHAR_EOF fi # end of overwriting check if test -f 'res20' then echo shar: will not over-write existing file "'res20'" else cat << "SHAR_EOF" > 'res20' NUMFUN = 1 F(X) = G(X).H(X). G(X) = 1.0D0 H(X) = SIN(X)/X INPUT PARAMETERS- Test 20: Use value of RFIRST other than zero AZERO PERIOD RFIRST EPS NQUAD 0.000000D+00 0.628319D+01 0.314159D+01 0.100000D-12 -9 NDIM1 NDIM2 100 14 ---------------------- PRINCIPAL OUTPUT PARAMETERS- RESULT Exact Difference 0.15707963D+01 0.15707963D+01 0.99920072D-14 ISTATE(1) ISTATE(2) ISTATE(3) ISTATE(4) ISTATE(5) (IFAIL) (LSIGCH) (N) (Column No.) (Row No.) 0 0 31 14 18 No. of function values (ISTATE(6)) = 279 ************************** Exit status ************************** Successful exit from OSCINT ***************************************************************** SHAR_EOF fi # end of overwriting check if test -f 'res3' then echo shar: will not over-write existing file "'res3'" else cat << "SHAR_EOF" > 'res3' NUMFUN = 3 F(X) = G(X).H(X). G(X) = 1.0D0 H(X) = EXP(-X*SINH(PI/2)*J1(X) INPUT PARAMETERS- Test 3: Successful execution of function 3. AZERO PERIOD RFIRST EPS NQUAD 0.000000D+00 0.628319D+01 0.000000D+00 0.100000D-12 -9 NDIM1 NDIM2 100 14 ---------------------- PRINCIPAL OUTPUT PARAMETERS- RESULT Exact Difference 0.82847664D-01 0.82847664D-01 0.10607210D-11 ISTATE(1) ISTATE(2) ISTATE(3) ISTATE(4) ISTATE(5) (IFAIL) (LSIGCH) (N) (Column No.) (Row No.) 0 1 7 2 6 No. of function values (ISTATE(6)) = 63 ************************** Exit status ************************** Successful exit from OSCINT ***************************************************************** SHAR_EOF fi # end of overwriting check if test -f 'res4' then echo shar: will not over-write existing file "'res4'" else cat << "SHAR_EOF" > 'res4' NUMFUN = 4 F(X) = G(X).H(X). G(X) = 1.0D0 H(X) = CONST(1 - EXP(X)) * J0(X) / X . (QUADPACK, PAGE 118) INPUT PARAMETERS- Test 4: Successful execution of function 4. AZERO PERIOD RFIRST EPS NQUAD 0.000000D+00 0.628319D+01 0.000000D+00 0.100000D-12 -9 NDIM1 NDIM2 100 14 ---------------------- PRINCIPAL OUTPUT PARAMETERS- RESULT Exact Difference 0.10000000D+01 0.10000000D+01 0.14099832D-13 ISTATE(1) ISTATE(2) ISTATE(3) ISTATE(4) ISTATE(5) (IFAIL) (LSIGCH) (N) (Column No.) (Row No.) 0 0 29 14 16 No. of function values (ISTATE(6)) = 261 ************************** Exit status ************************** Successful exit from OSCINT ***************************************************************** SHAR_EOF fi # end of overwriting check if test -f 'res5' then echo shar: will not over-write existing file "'res5'" else cat << "SHAR_EOF" > 'res5' NUMFUN = 5 F(X) = G(X).H(X). G(X) = 1.0D0 H(X) = (1 + 100*EXP(-X**2/PI**2))*J1(X) INPUT PARAMETERS- Test 5: Successful execution of function 5. AZERO PERIOD RFIRST EPS NQUAD 0.000000D+00 0.628319D+01 0.000000D+00 0.100000D-12 -9 NDIM1 NDIM2 100 14 ---------------------- PRINCIPAL OUTPUT PARAMETERS- RESULT Exact Difference 0.92519503D+02 0.92519503D+02 0.76170181D-11 ISTATE(1) ISTATE(2) ISTATE(3) ISTATE(4) ISTATE(5) (IFAIL) (LSIGCH) (N) (Column No.) (Row No.) 0 0 30 14 17 No. of function values (ISTATE(6)) = 270 ************************** Exit status ************************** Successful exit from OSCINT ***************************************************************** SHAR_EOF fi # end of overwriting check if test -f 'res6' then echo shar: will not over-write existing file "'res6'" else cat << "SHAR_EOF" > 'res6' NUMFUN = 6 F(X) = G(X).H(X). G(X) = (COS**4(X) - 3/8) H(X) = 1.0D0/(1 + X**2) INPUT PARAMETERS- Test 6: Successful execution of function 6. AZERO PERIOD RFIRST EPS NQUAD 0.000000D+00 0.314159D+01 0.000000D+00 0.100000D-10 -9 NDIM1 NDIM2 200 14 ---------------------- PRINCIPAL OUTPUT PARAMETERS- RESULT Exact Difference 0.10988835D+00 0.10988835D+00 0.43723067D-09 ISTATE(1) ISTATE(2) ISTATE(3) ISTATE(4) ISTATE(5) (IFAIL) (LSIGCH) (N) (Column No.) (Row No.) 0 0 190 4 187 No. of function values (ISTATE(6)) = 1710 ************************** Exit status ************************** Successful exit from OSCINT ***************************************************************** SHAR_EOF fi # end of overwriting check if test -f 'res7' then echo shar: will not over-write existing file "'res7'" else cat << "SHAR_EOF" > 'res7' NUMFUN = 1 F(X) = G(X).H(X). G(X) = 1.0D0 H(X) = SIN(X)/X INPUT PARAMETERS- Test 7: Successful execution of function 1 (extra output). AZERO PERIOD RFIRST EPS NQUAD 0.000000D+00 0.628319D+01 0.000000D+00 0.100000D-12 -9 NDIM1 NDIM2 100 14 ---------------------- Code check output from HFUN, GPER, QRULE X = 0.50013779D-01 HFUN = 0.99958316D+00 X = 0.50013779D-01 GPER = 0.10000000D+01 X = 0.25756173D+00 HFUN = 0.98898027D+00 X = 0.25756173D+00 GPER = 0.10000000D+01 X = 0.60731473D+00 HFUN = 0.93965187D+00 X = 0.60731473D+00 GPER = 0.10000000D+01 X = 0.10614602D+01 HFUN = 0.82251637D+00 X = 0.10614602D+01 GPER = 0.10000000D+01 X = 0.15707963D+01 HFUN = 0.63661977D+00 X = 0.15707963D+01 GPER = 0.10000000D+01 X = 0.20801324D+01 HFUN = 0.41971771D+00 X = 0.20801324D+01 GPER = 0.10000000D+01 X = 0.25342779D+01 HFUN = 0.22517831D+00 X = 0.25342779D+01 GPER = 0.10000000D+01 X = 0.28840309D+01 HFUN = 0.88322033D-01 X = 0.28840309D+01 GPER = 0.10000000D+01 X = 0.30915789D+01 HFUN = 0.16170679D-01 X = 0.30915789D+01 GPER = 0.10000000D+01 X = 0.31916064D+01 GPER = 0.10000000D+01 X = 0.31916064D+01 HFUN = -0.15663877D-01 X = 0.33991544D+01 GPER = 0.10000000D+01 X = 0.33991544D+01 HFUN = -0.74937307D-01 X = 0.37489074D+01 GPER = 0.10000000D+01 X = 0.37489074D+01 HFUN = -0.15222153D+00 PRINCIPAL OUTPUT PARAMETERS- RESULT Exact Difference 0.15707963D+01 0.15707963D+01 0.99920072D-14 ISTATE(1) ISTATE(2) ISTATE(3) ISTATE(4) ISTATE(5) (IFAIL) (LSIGCH) (N) (Column No.) (Row No.) 0 0 31 14 18 No. of function values (ISTATE(6)) = 279 ************************** Exit status ************************** Successful exit from OSCINT ***************************************************************** Finite Average Table -------------------- 1 0.19D+01 0.71D+00 0.31D+00 0.14D+00 0.67D-01 0.32D-01 0.15D-01 0.74D-02 0.36D-02 0.17D-02 0.85D-03 0.42D-03 0.20D-03 0.10D-03 2 -0.43D+00 -0.89D-01 -0.26D-01 -0.87D-02 -0.32D-02 -0.13D-02 -0.51D-03 -0.21D-03 -0.91D-04 -0.39D-04 -0.17D-04 -0.77D-05 -0.34D-05 -0.16D-05 3 0.26D+00 0.37D-01 0.83D-02 0.23D-02 0.71D-03 0.24D-03 0.85D-04 0.32D-04 0.12D-04 0.48D-05 0.19D-05 0.79D-06 0.33D-06 0.14D-06 4 -0.18D+00 -0.20D-01 -0.37D-02 -0.87D-03 -0.23D-03 -0.69D-04 -0.22D-04 -0.73D-05 -0.26D-05 -0.93D-06 -0.35D-06 -0.13D-06 -0.51D-07 -0.20D-07 5 0.14D+00 0.13D-01 0.20D-02 0.40D-03 0.95D-04 0.25D-04 0.72D-05 0.22D-05 0.71D-06 0.24D-06 0.82D-07 0.29D-07 0.11D-07 0.40D-08 6 -0.12D+00 -0.89D-02 -0.12D-02 -0.21D-03 -0.45D-04 -0.11D-04 -0.28D-05 -0.79D-06 -0.23D-06 -0.73D-07 -0.24D-07 -0.79D-08 -0.27D-08 -0.96D-09 7 0.98D-01 0.65D-02 0.77D-03 0.12D-03 0.23D-04 0.51D-05 0.12D-05 0.32D-06 0.88D-07 0.26D-07 0.78D-08 0.25D-08 0.80D-09 0.27D-09 8 -0.85D-01 -0.50D-02 -0.53D-03 -0.76D-04 -0.13D-04 -0.26D-05 -0.59D-06 -0.14D-06 -0.37D-07 -0.10D-07 -0.29D-08 -0.86D-09 -0.27D-09 -0.84D-10 9 0.75D-01 0.39D-02 0.38D-03 0.49D-04 0.79D-05 0.15D-05 0.30D-06 0.69D-07 0.17D-07 0.43D-08 0.12D-08 0.33D-09 0.97D-10 0.29D-10 10 -0.67D-01 -0.32D-02 -0.28D-03 -0.33D-04 -0.50D-05 -0.86D-06 -0.17D-06 -0.35D-07 -0.81D-08 -0.20D-08 -0.51D-09 -0.14D-09 -0.38D-10 -0.11D-10 11 0.61D-01 0.26D-02 0.21D-03 0.24D-04 0.32D-05 0.52D-06 0.96D-07 0.19D-07 0.41D-08 0.96D-09 0.23D-09 0.60D-10 0.16D-10 0.44D-11 12 -0.55D-01 -0.22D-02 -0.16D-03 -0.17D-04 -0.22D-05 -0.33D-06 -0.57D-07 -0.11D-07 -0.22D-08 -0.49D-09 -0.11D-09 -0.28D-10 -0.71D-11 -0.19D-11 13 0.51D-01 0.19D-02 0.13D-03 0.13D-04 0.15D-05 0.22D-06 0.36D-07 0.64D-08 0.12D-08 0.26D-09 0.58D-10 0.14D-10 0.33D-11 0.85D-12 14 -0.47D-01 -0.16D-02 -0.11D-03 -0.96D-05 -0.11D-05 -0.15D-06 -0.23D-07 -0.39D-08 -0.73D-09 -0.15D-09 -0.31D-10 -0.70D-11 -0.16D-11 -0.40D-12 15 0.44D-01 0.14D-02 0.86D-04 0.74D-05 0.80D-06 0.10D-06 0.15D-07 0.24D-08 0.43D-09 0.83D-10 0.17D-10 0.37D-11 0.83D-12 0.20D-12 16 -0.41D-01 -0.12D-02 -0.71D-04 -0.58D-05 -0.59D-06 -0.72D-07 -0.10D-07 -0.16D-08 -0.27D-09 -0.49D-10 -0.97D-11 -0.20D-11 -0.44D-12 -0.10D-12 17 0.39D-01 0.11D-02 0.60D-04 0.46D-05 0.45D-06 0.52D-07 0.70D-08 0.10D-08 0.17D-09 0.30D-10 0.57D-11 0.11D-11 0.24D-12 0.53D-13 18 -0.36D-01 -0.98D-03 -0.50D-04 -0.37D-05 -0.34D-06 -0.38D-07 -0.49D-08 -0.70D-09 -0.11D-09 -0.19D-10 -0.34D-11 -0.65D-12 -0.13D-12 -0.28D-13 19 0.34D-01 0.88D-03 0.43D-04 0.30D-05 0.27D-06 0.28D-07 0.35D-08 0.48D-09 0.72D-10 0.12D-10 0.21D-11 0.39D-12 0.76D-13 20 -0.33D-01 -0.80D-03 -0.37D-04 -0.25D-05 -0.21D-06 -0.21D-07 -0.25D-08 -0.33D-09 -0.49D-10 -0.77D-11 -0.13D-11 -0.24D-12 21 0.31D-01 0.72D-03 0.32D-04 0.21D-05 0.17D-06 0.16D-07 0.19D-08 0.24D-09 0.33D-10 0.51D-11 0.83D-12 22 -0.30D-01 -0.66D-03 -0.28D-04 -0.17D-05 -0.13D-06 -0.13D-07 -0.14D-08 -0.17D-09 -0.23D-10 -0.34D-11 23 0.28D-01 0.60D-03 0.25D-04 0.14D-05 0.11D-06 0.99D-08 0.10D-08 0.12D-09 0.16D-10 24 -0.27D-01 -0.55D-03 -0.22D-04 -0.12D-05 -0.89D-07 -0.78D-08 -0.80D-09 -0.92D-10 25 0.26D-01 0.51D-03 0.19D-04 0.10D-05 0.74D-07 0.62D-08 0.61D-09 26 -0.25D-01 -0.47D-03 -0.17D-04 -0.90D-06 -0.61D-07 -0.50D-08 27 0.24D-01 0.44D-03 0.15D-04 0.78D-06 0.51D-07 28 -0.23D-01 -0.41D-03 -0.14D-04 -0.68D-06 29 0.22D-01 0.38D-03 0.12D-04 30 -0.22D-01 -0.35D-03 31 0.21D-01 SHAR_EOF fi # end of overwriting check if test -f 'res8' then echo shar: will not over-write existing file "'res8'" else cat << "SHAR_EOF" > 'res8' NUMFUN = 2 F(X) = G(X).H(X). G(X) = COS(X) H(X) = 1.0D0/(1+X**2) INPUT PARAMETERS- Test 8: Function 2 with NQUAD = -9. Generates incorrect result. AZERO PERIOD RFIRST EPS NQUAD 0.000000D+00 0.628319D+01 0.000000D+00 0.100000D-12 -9 NDIM1 NDIM2 100 14 ---------------------- PRINCIPAL OUTPUT PARAMETERS- RESULT Exact Difference 0.57786243D+00 0.57786367D+00 0.12492351D-05 ISTATE(1) ISTATE(2) ISTATE(3) ISTATE(4) ISTATE(5) (IFAIL) (LSIGCH) (N) (Column No.) (Row No.) 0 0 27 13 15 No. of function values (ISTATE(6)) = 243 ************************** Exit status ************************** Successful exit from OSCINT This quadrature rule is inadequate. Note the steady convergence of the table elements below to a result which is in fact wrong. Routine cannot know this. ***************************************************************** Finite Average Table -------------------- 1 0.61D+00 0.29D+00 0.14D+00 0.64D-01 0.31D-01 0.15D-01 0.71D-02 0.35D-02 0.17D-02 0.81D-03 0.39D-03 0.19D-03 0.93D-04 0.46D-04 2 -0.39D-01 -0.15D-01 -0.64D-02 -0.27D-02 -0.12D-02 -0.53D-03 -0.24D-03 -0.11D-03 -0.50D-04 -0.23D-04 -0.11D-04 -0.50D-05 -0.23D-05 -0.11D-05 3 0.84D-02 0.27D-02 0.93D-03 0.34D-03 0.13D-03 0.53D-04 0.22D-04 0.90D-05 0.38D-05 0.16D-05 0.71D-06 0.31D-06 0.14D-06 0.61D-07 4 -0.30D-02 -0.80D-03 -0.24D-03 -0.78D-04 -0.27D-04 -0.97D-05 -0.36D-05 -0.14D-05 -0.54D-06 -0.22D-06 -0.88D-07 -0.36D-07 -0.15D-07 -0.64D-08 5 0.14D-02 0.32D-03 0.84D-04 0.24D-04 0.75D-05 0.25D-05 0.85D-06 0.30D-06 0.11D-06 0.41D-07 0.16D-07 0.60D-08 0.24D-08 0.95D-09 6 -0.78D-03 -0.15D-03 -0.36D-04 -0.92D-05 -0.26D-05 -0.78D-06 -0.25D-06 -0.81D-07 -0.28D-07 -0.97D-08 -0.35D-08 -0.13D-08 -0.48D-09 -0.18D-09 7 0.47D-03 0.82D-04 0.17D-04 0.40D-05 0.10D-05 0.29D-06 0.84D-07 0.26D-07 0.83D-08 0.27D-08 0.92D-09 0.32D-09 0.11D-09 0.41D-10 8 -0.31D-03 -0.48D-04 -0.91D-05 -0.20D-05 -0.46D-06 -0.12D-06 -0.33D-07 -0.94D-08 -0.28D-08 -0.88D-09 -0.28D-09 -0.93D-10 -0.31D-10 -0.11D-10 9 0.21D-03 0.30D-04 0.52D-05 0.10D-05 0.23D-06 0.54D-07 0.14D-07 0.37D-08 0.11D-08 0.31D-09 0.96D-10 0.30D-10 0.97D-11 0.32D-11 10 -0.15D-03 -0.20D-04 -0.31D-05 -0.57D-06 -0.12D-06 -0.26D-07 -0.64D-08 -0.16D-08 -0.44D-09 -0.12D-09 -0.35D-10 -0.11D-10 -0.33D-11 -0.10D-11 11 0.11D-03 0.13D-04 0.20D-05 0.34D-06 0.65D-07 0.14D-07 0.31D-08 0.75D-09 0.19D-09 0.51D-10 0.14D-10 0.41D-11 0.12D-11 0.37D-12 12 -0.85D-04 -0.94D-05 -0.13D-05 -0.21D-06 -0.37D-07 -0.75D-08 -0.16D-08 -0.37D-09 -0.89D-10 -0.23D-10 -0.60D-11 -0.17D-11 -0.47D-12 -0.14D-12 13 0.66D-04 0.68D-05 0.88D-06 0.13D-06 0.23D-07 0.43D-08 0.87D-09 0.19D-09 0.44D-10 0.11D-10 0.27D-11 0.72D-12 0.20D-12 0.55D-13 14 -0.52D-04 -0.51D-05 -0.61D-06 -0.87D-07 -0.14D-07 -0.25D-08 -0.49D-09 -0.10D-09 -0.23D-10 -0.53D-11 -0.13D-11 -0.32D-12 -0.85D-13 15 0.42D-04 0.38D-05 0.44D-06 0.59D-07 0.90D-08 0.15D-08 0.28D-09 0.57D-10 0.12D-10 0.27D-11 0.63D-12 0.15D-12 0.39D-13 16 -0.35D-04 -0.30D-05 -0.32D-06 -0.41D-07 -0.60D-08 -0.97D-09 -0.17D-09 -0.33D-10 -0.67D-11 -0.14D-11 -0.32D-12 -0.76D-13 17 0.29D-04 0.23D-05 0.24D-06 0.29D-07 0.40D-08 0.62D-09 0.11D-09 0.19D-10 0.38D-11 0.78D-12 0.17D-12 18 -0.24D-04 -0.18D-05 -0.18D-06 -0.21D-07 -0.28D-08 -0.41D-09 -0.67D-10 -0.12D-10 -0.22D-11 -0.44D-12 19 0.20D-04 0.15D-05 0.14D-06 0.15D-07 0.20D-08 0.28D-09 0.44D-10 0.74D-11 0.13D-11 20 -0.17D-04 -0.12D-05 -0.11D-06 -0.11D-07 -0.14D-08 -0.19D-09 -0.29D-10 -0.47D-11 21 0.15D-04 0.10D-05 0.85D-07 0.86D-08 0.10D-08 0.13D-09 0.19D-10 22 -0.13D-04 -0.83D-06 -0.67D-07 -0.66D-08 -0.75D-09 -0.95D-10 23 0.11D-04 0.69D-06 0.54D-07 0.51D-08 0.56D-09 24 -0.99D-05 -0.58D-06 -0.44D-07 -0.40D-08 25 0.88D-05 0.50D-06 0.36D-07 26 -0.78D-05 -0.42D-06 27 0.69D-05 SHAR_EOF fi # end of overwriting check if test -f 'res9' then echo shar: will not over-write existing file "'res9'" else cat << "SHAR_EOF" > 'res9' NUMFUN = 3 F(X) = G(X).H(X). G(X) = 1.0D0 H(X) = EXP(-X*SINH(PI/2)*J1(X) INPUT PARAMETERS- Test 9: Set NDIM = 5 to illustrate non convergence. AZERO PERIOD RFIRST EPS NQUAD 0.000000D+00 0.628319D+01 0.000000D+00 0.100000D-12 -9 NDIM1 NDIM2 5 14 ---------------------- PRINCIPAL OUTPUT PARAMETERS- ISTATE(1) ISTATE(2) ISTATE(3) ISTATE(4) ISTATE(5) (IFAIL) (LSIGCH) (N) (Column No.) (Row No.) -100 1 5 2 4 No. of function values (ISTATE(6)) = 45 ************************** Exit status ************************** ERROR: Maximum number of intervals exceeded ***************************************************************** SHAR_EOF fi # end of overwriting check if test -f 'src_715.f' then echo shar: will not over-write existing file "'src_715.f'" else cat << "SHAR_EOF" > 'src_715.f' C--**--CH3394--715--P:Gen--30:3:2000 C--**--CH2209--715--C:Gen--12:9:1999 CS REAL FUNCTION ANORM(ARG) DOUBLE PRECISION FUNCTION ANORM(ARG) C------------------------------------------------------------------ C C This function evaluates the normal distribution function: C C / x C 1 | -t*t/2 C P(x) = ----------- | e dt C sqrt(2 pi) | C /-oo C C The main computation evaluates near-minimax approximations C derived from those in "Rational Chebyshev approximations for C the error function" by W. J. Cody, Math. Comp., 1969, 631-637. C This transportable program uses rational functions that C theoretically approximate the normal distribution function to C at least 18 significant decimal digits. The accuracy achieved C depends on the arithmetic system, the compiler, the intrinsic C functions, and proper selection of the machine-dependent C constants. C C******************************************************************* C******************************************************************* C C Explanation of machine-dependent constants. Let C C XMIN = the smallest positive floating-point number. C C Then the following machine-dependent constants must be declared C in DATA statements. IEEE values are provided as a default. C C EPS = argument below which anorm(x) may be represented by C 0.5 and above which x*x will not underflow. C A conservative value is the largest machine number X C such that 1.0 + X = 1.0 to machine precision. C XLOW = the most negative argument for which ANORM does not C vanish. This is the negative of the solution to C W(x) * (1-1/x**2) = XMIN, C where W(x) = exp(-x*x/2)/[x*sqrt(2*pi)]. C XUPPR = positive argument beyond which anorm = 1.0. A C conservative value is the solution to the equation C exp(-x*x/2) = EPS, C i.e., XUPPR = sqrt[-2 ln(eps)]. C C Approximate values for some important machines are: C C XMIN EPS XLOW XUPPR C C CDC 7600 (S.P.) 3.13E-294 7.11E-15 -36.641 8.072 C CRAY-1 (S.P.) 4.58E-246 7.11E-157 -106.521 26.816 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 1.18E-38 5.96E-8 -12.949 5.768 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 2.23D-308 1.11D-16 -37.519 8.572 C IBM 195 (D.P.) 5.40D-79 1.39D-17 -18.781 8.811 C VAX D-Format (D.P.) 2.94D-39 1.39D-17 -13.055 8.811 C VAX G-Format (D.P.) 5.56D-309 1.11D-16 -37.556 8.572 C C******************************************************************* C******************************************************************* C C Error returns C C The program returns ANORM = 0 for ARG .LE. XLOW. C C C Intrinsic functions required are: C C ABS, AINT, EXP C C C Author: W. J. Cody C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C C Latest modification: March 15, 1992 C C------------------------------------------------------------------ INTEGER I CS REAL DOUBLE PRECISION 1 A,ARG,B,C,D,DEL,EPS,HALF,P,ONE,Q,RESULT,SIXTEN, 2 SQRPI,THRSH,ROOT32,X,XLOW,XDEN,XNUM,Y,XSQ,XUPPR,ZERO DIMENSION A(5),B(4),C(9),D(8),P(6),Q(5) C------------------------------------------------------------------ C Mathematical constants C C SQRPI = 1 / sqrt(2*pi), ROOT32 = sqrt(32), and C THRSH is the argument for which anorm = 0.75. C------------------------------------------------------------------ CS DATA ONE,HALF,ZERO,SIXTEN/1.0E0,0.5E0,0.0E0,1.60E1/, CS 1 SQRPI/3.9894228040143267794E-1/,THRSH/0.66291E0/, CS 2 ROOT32/5.656854248E0/ DATA ONE,HALF,ZERO,SIXTEN/1.0D0,0.5D0,0.0D0,1.60D1/, 1 SQRPI/3.9894228040143267794D-1/,THRSH/0.66291D0/, 2 ROOT32/5.656854248D0/ C------------------------------------------------------------------ C Machine-dependent constants C------------------------------------------------------------------ CS DATA EPS/5.96E-8/,XLOW/-12.949E0/,XUPPR/5.768E0/ DATA EPS/1.11D-16/,XLOW/-37.519D0/,XUPPR/8.572D0/ C------------------------------------------------------------------ C Coefficients for approximation in first interval C------------------------------------------------------------------ CS DATA A/2.2352520354606839287E00,1.6102823106855587881E02, CS 1 1.0676894854603709582E03,1.8154981253343561249E04, CS 2 6.5682337918207449113E-2/ CS DATA B/4.7202581904688241870E01,9.7609855173777669322E02, CS 1 1.0260932208618978205E04,4.5507789335026729956E04/ DATA A/2.2352520354606839287D00,1.6102823106855587881D02, 1 1.0676894854603709582D03,1.8154981253343561249D04, 2 6.5682337918207449113D-2/ DATA B/4.7202581904688241870D01,9.7609855173777669322D02, 1 1.0260932208618978205D04,4.5507789335026729956D04/ C------------------------------------------------------------------ C Coefficients for approximation in second interval C------------------------------------------------------------------ CS DATA C/3.9894151208813466764E-1,8.8831497943883759412E00, CS 1 9.3506656132177855979E01,5.9727027639480026226E02, CS 2 2.4945375852903726711E03,6.8481904505362823326E03, CS 3 1.1602651437647350124E04,9.8427148383839780218E03, CS 4 1.0765576773720192317E-8/ CS DATA D/2.2266688044328115691E01,2.3538790178262499861E02, CS 1 1.5193775994075548050E03,6.4855582982667607550E03, CS 2 1.8615571640885098091E04,3.4900952721145977266E04, CS 3 3.8912003286093271411E04,1.9685429676859990727E04/ DATA C/3.9894151208813466764D-1,8.8831497943883759412D00, 1 9.3506656132177855979D01,5.9727027639480026226D02, 2 2.4945375852903726711D03,6.8481904505362823326D03, 3 1.1602651437647350124D04,9.8427148383839780218D03, 4 1.0765576773720192317D-8/ DATA D/2.2266688044328115691D01,2.3538790178262499861D02, 1 1.5193775994075548050D03,6.4855582982667607550D03, 2 1.8615571640885098091D04,3.4900952721145977266D04, 3 3.8912003286093271411D04,1.9685429676859990727D04/ C------------------------------------------------------------------ C Coefficients for approximation in third interval C------------------------------------------------------------------ CS DATA P/2.1589853405795699E-1,1.274011611602473639E-1, CS 1 2.2235277870649807E-2,1.421619193227893466E-3, CS 2 2.9112874951168792E-5,2.307344176494017303E-2/ CS DATA Q/1.28426009614491121E00,4.68238212480865118E-1, CS 1 6.59881378689285515E-2,3.78239633202758244E-3, CS 2 7.29751555083966205E-5/ DATA P/2.1589853405795699D-1,1.274011611602473639D-1, 1 2.2235277870649807D-2,1.421619193227893466D-3, 2 2.9112874951168792D-5,2.307344176494017303D-2/ DATA Q/1.28426009614491121D00,4.68238212480865118D-1, 1 6.59881378689285515D-2,3.78239633202758244D-3, 2 7.29751555083966205D-5/ C------------------------------------------------------------------ X = ARG Y = ABS(X) IF (Y .LE. THRSH) THEN C------------------------------------------------------------------ C Evaluate anorm for |X| <= 0.66291 C------------------------------------------------------------------ XSQ = ZERO IF (Y .GT. EPS) XSQ = X * X XNUM = A(5)*XSQ XDEN = XSQ DO 20 I = 1, 3 XNUM = (XNUM + A(I)) * XSQ XDEN = (XDEN + B(I)) * XSQ 20 CONTINUE RESULT = X * (XNUM + A(4)) / (XDEN + B(4)) RESULT = HALF + RESULT C------------------------------------------------------------------ C Evaluate anorm for 0.66291 <= |X| <= sqrt(32) C------------------------------------------------------------------ ELSE IF (Y .LE. ROOT32) THEN XNUM = C(9)*Y XDEN = Y DO 120 I = 1, 7 XNUM = (XNUM + C(I)) * Y XDEN = (XDEN + D(I)) * Y 120 CONTINUE RESULT = (XNUM + C(8)) / (XDEN + D(8)) XSQ = AINT(Y*SIXTEN)/SIXTEN DEL = (Y-XSQ)*(Y+XSQ) RESULT = EXP(-XSQ*XSQ*HALF)*EXP(-DEL*HALF)*RESULT IF (X .GT. ZERO) RESULT = ONE - RESULT C------------------------------------------------------------------ C Evaluate anorm for |X| > sqrt(32) C------------------------------------------------------------------ ELSE RESULT = ZERO IF ((X .GE. XLOW) .AND. (X .LT. XUPPR)) THEN XSQ = ONE / (X * X) XNUM = P(6)*XSQ XDEN = XSQ DO 240 I = 1, 4 XNUM = (XNUM + P(I)) * XSQ XDEN = (XDEN + Q(I)) * XSQ 240 CONTINUE RESULT = XSQ *(XNUM + P(5)) / (XDEN + Q(5)) RESULT = (SQRPI - RESULT) / Y XSQ = AINT(X*SIXTEN)/SIXTEN DEL = (X-XSQ)*(X+XSQ) RESULT = EXP(-XSQ*XSQ*HALF)*EXP(-DEL*HALF)*RESULT END IF IF (X .GT. ZERO) RESULT = ONE - RESULT END IF C------------------------------------------------------------------ C Fix up for negative argument, erf, etc. C------------------------------------------------------------------ ANORM = RESULT C---------- Last card of ANORM ---------- END SUBROUTINE CALCI0(ARG,RESULT,JINT) C-------------------------------------------------------------------- C C This packet computes modified Bessel functions of the first kind C and order zero, I0(X) and EXP(-ABS(X))*I0(X), for real C arguments X. It contains two function type subprograms, BESI0 C and BESEI0, and one subroutine type subprogram, CALCI0. C The calling statements for the primary entries are C C Y=BESI0(X) C and C Y=BESEI0(X) C C where the entry points correspond to the functions I0(X) and C EXP(-ABS(X))*I0(X), respectively. The routine CALCI0 is C intended for internal packet use only, all computations within C the packet being concentrated in this routine. The function C subprograms invoke CALCI0 with the statement C CALL CALCI0(ARG,RESULT,JINT) C where the parameter usage is as follows C C Function Parameters for CALCI0 C Call ARG RESULT JINT C C BESI0(ARG) ABS(ARG) .LE. XMAX I0(ARG) 1 C BESEI0(ARG) any real ARG EXP(-ABS(ARG))*I0(ARG) 2 C C The main computation evaluates slightly modified forms of C minimax approximations generated by Blair and Edwards, Chalk C River (Atomic Energy of Canada Limited) Report AECL-4928, C October, 1974. This transportable program is patterned after C the machine-dependent FUNPACK packet NATSI0, but cannot match C that version for efficiency or accuracy. This version uses C rational functions that theoretically approximate I-SUB-0(X) C to at least 18 significant decimal digits. The accuracy C achieved depends on the arithmetic system, the compiler, the C intrinsic functions, and proper selection of the machine- C dependent constants. C C******************************************************************* C******************************************************************* C C Explanation of machine-dependent constants. Let C C beta = Radix for the floating-point system C maxexp = Smallest power of beta that overflows C C Then the following machine-dependent constants must be declared C in DATA statements. IEEE values are provided as a default. C C XSMALL = Positive argument such that 1.0 - X = 1.0 to C machine precision for all ABS(X) .LE. XSMALL. C XINF = Largest positive machine number; approximately C beta**maxexp C XMAX = Largest argument acceptable to BESI0; Solution to C equation: C W(X) * (1+1/(8*X)+9/(128*X**2) = beta**maxexp C where W(X) = EXP(X)/SQRT(2*PI*X) C C C Approximate values for some important machines are: C C beta maxexp XSMALL C C CRAY-1 (S.P.) 2 8191 3.55E-15 C Cyber 180/855 C under NOS (S.P.) 2 1070 3.55E-15 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 2 128 2.98E-8 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 2 1024 5.55D-17 C IBM 3033 (D.P.) 16 63 6.95D-18 C VAX (S.P.) 2 127 2.98E-8 C VAX D-Format (D.P.) 2 127 6.95D-18 C VAX G-Format (D.P.) 2 1023 5.55D-17 C C C XINF XMAX C C CRAY-1 (S.P.) 5.45E+2465 5682.810 C Cyber 180/855 C under NOS (S.P.) 1.26E+322 745.893 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 3.40E+38 91.900 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 1.79D+308 713.986 C IBM 3033 (D.P.) 7.23D+75 178.182 C VAX (S.P.) 1.70D+38 91.203 C VAX D-Format (D.P.) 1.70D+38 91.203 C VAX G-Format (D.P.) 8.98D+307 713.293 C C******************************************************************* C******************************************************************* C C Error returns C C The program returns XINF for BESI0 for ABS(ARG) .GT. XMAX. C C C Intrinsic functions required are: C C ABS, SQRT, EXP C C C Authors: W. J. Cody and L. Stoltz C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C C Latest modification: March 12, 1992 C C-------------------------------------------------------------------- INTEGER I,JINT CS REAL DOUBLE PRECISION 1 A,ARG,B,EXP40,FORTY,ONE,ONE5,P,PP,Q,QQ,RESULT, 2 REC15,SUMP,SUMQ,TWO25,X,XINF,XMAX,XSMALL,XX DIMENSION P(15),PP(8),Q(5),QQ(7) C-------------------------------------------------------------------- C Mathematical constants C-------------------------------------------------------------------- CS DATA ONE/1.0E0/,ONE5/15.0E0/,EXP40/2.353852668370199854E17/, CS 1 FORTY/40.0E0/,REC15/6.6666666666666666666E-2/, CS 2 TWO25/225.0E0/ DATA ONE/1.0D0/,ONE5/15.0D0/,EXP40/2.353852668370199854D17/, 1 FORTY/40.0D0/,REC15/6.6666666666666666666D-2/, 2 TWO25/225.0D0/ C-------------------------------------------------------------------- C Machine-dependent constants C-------------------------------------------------------------------- CS DATA XSMALL/2.98E-8/,XINF/3.40E38/,XMAX/91.9E0/ DATA XSMALL/5.55D-17/,XINF/1.79D308/,XMAX/713.986D0/ C-------------------------------------------------------------------- C Coefficients for XSMALL .LE. ABS(ARG) .LT. 15.0 C-------------------------------------------------------------------- CS DATA P/-5.2487866627945699800E-18,-1.5982226675653184646E-14, CS 1 -2.6843448573468483278E-11,-3.0517226450451067446E-08, CS 2 -2.5172644670688975051E-05,-1.5453977791786851041E-02, CS 3 -7.0935347449210549190E+00,-2.4125195876041896775E+03, CS 4 -5.9545626019847898221E+05,-1.0313066708737980747E+08, CS 5 -1.1912746104985237192E+10,-8.4925101247114157499E+11, CS 6 -3.2940087627407749166E+13,-5.5050369673018427753E+14, CS 7 -2.2335582639474375249E+15/ CS DATA Q/-3.7277560179962773046E+03, 6.5158506418655165707E+06, CS 1 -6.5626560740833869295E+09, 3.7604188704092954661E+12, CS 2 -9.7087946179594019126E+14/ DATA P/-5.2487866627945699800D-18,-1.5982226675653184646D-14, 1 -2.6843448573468483278D-11,-3.0517226450451067446D-08, 2 -2.5172644670688975051D-05,-1.5453977791786851041D-02, 3 -7.0935347449210549190D+00,-2.4125195876041896775D+03, 4 -5.9545626019847898221D+05,-1.0313066708737980747D+08, 5 -1.1912746104985237192D+10,-8.4925101247114157499D+11, 6 -3.2940087627407749166D+13,-5.5050369673018427753D+14, 7 -2.2335582639474375249D+15/ DATA Q/-3.7277560179962773046D+03, 6.5158506418655165707D+06, 1 -6.5626560740833869295D+09, 3.7604188704092954661D+12, 2 -9.7087946179594019126D+14/ C-------------------------------------------------------------------- C Coefficients for 15.0 .LE. ABS(ARG) C-------------------------------------------------------------------- CS DATA PP/-3.9843750000000000000E-01, 2.9205384596336793945E+00, CS 1 -2.4708469169133954315E+00, 4.7914889422856814203E-01, CS 2 -3.7384991926068969150E-03,-2.6801520353328635310E-03, CS 3 9.9168777670983678974E-05,-2.1877128189032726730E-06/ CS DATA QQ/-3.1446690275135491500E+01, 8.5539563258012929600E+01, CS 1 -6.0228002066743340583E+01, 1.3982595353892851542E+01, CS 2 -1.1151759188741312645E+00, 3.2547697594819615062E-02, CS 3 -5.5194330231005480228E-04/ DATA PP/-3.9843750000000000000D-01, 2.9205384596336793945D+00, 1 -2.4708469169133954315D+00, 4.7914889422856814203D-01, 2 -3.7384991926068969150D-03,-2.6801520353328635310D-03, 3 9.9168777670983678974D-05,-2.1877128189032726730D-06/ DATA QQ/-3.1446690275135491500D+01, 8.5539563258012929600D+01, 1 -6.0228002066743340583D+01, 1.3982595353892851542D+01, 2 -1.1151759188741312645D+00, 3.2547697594819615062D-02, 3 -5.5194330231005480228D-04/ C-------------------------------------------------------------------- X = ABS(ARG) IF (X .LT. XSMALL) THEN RESULT = ONE ELSE IF (X .LT. ONE5) THEN C-------------------------------------------------------------------- C XSMALL .LE. ABS(ARG) .LT. 15.0 C-------------------------------------------------------------------- XX = X * X SUMP = P(1) DO 50 I = 2, 15 SUMP = SUMP * XX + P(I) 50 CONTINUE XX = XX - TWO25 SUMQ = ((((XX+Q(1))*XX+Q(2))*XX+Q(3))*XX+Q(4))*XX+Q(5) RESULT = SUMP / SUMQ IF (JINT .EQ. 2) RESULT = RESULT * EXP(-X) ELSE IF (X .GE. ONE5) THEN IF ((JINT .EQ. 1) .AND. (X .GT. XMAX)) THEN RESULT = XINF ELSE C-------------------------------------------------------------------- C 15.0 .LE. ABS(ARG) C-------------------------------------------------------------------- XX = ONE / X - REC15 SUMP = ((((((PP(1)*XX+PP(2))*XX+PP(3))*XX+PP(4))*XX+ 1 PP(5))*XX+PP(6))*XX+PP(7))*XX+PP(8) SUMQ = ((((((XX+QQ(1))*XX+QQ(2))*XX+QQ(3))*XX+ 1 QQ(4))*XX+QQ(5))*XX+QQ(6))*XX+QQ(7) RESULT = SUMP / SUMQ IF (JINT .EQ. 2) THEN RESULT = (RESULT - PP(1)) / SQRT(X) ELSE C-------------------------------------------------------------------- C Calculation reformulated to avoid premature overflow C-------------------------------------------------------------------- IF (X .LE.(XMAX-ONE5)) THEN A = EXP(X) B = ONE ELSE A = EXP(X-FORTY) B = EXP40 END IF RESULT = ((RESULT*A-PP(1)*A)/SQRT(X))*B END IF END IF END IF C-------------------------------------------------------------------- C Return for ABS(ARG) .LT. XSMALL C-------------------------------------------------------------------- RETURN C----------- Last line of CALCI0 ----------- END CS REAL FUNCTION BESI0(X) DOUBLE PRECISION FUNCTION BESI0(X) C-------------------------------------------------------------------- C C This long precision subprogram computes approximate values for C modified Bessel functions of the first kind of order zero for C arguments ABS(ARG) .LE. XMAX (see comments heading CALCI0). C C-------------------------------------------------------------------- INTEGER JINT CS REAL DOUBLE PRECISION 1 X, RESULT C-------------------------------------------------------------------- JINT=1 CALL CALCI0(X,RESULT,JINT) BESI0=RESULT RETURN C---------- Last line of BESI0 ---------- END CS REAL FUNCTION BESEI0(X) DOUBLE PRECISION FUNCTION BESEI0(X) C-------------------------------------------------------------------- C C This function program computes approximate values for the C modified Bessel function of the first kind of order zero C multiplied by EXP(-ABS(X)), where EXP is the C exponential function, ABS is the absolute value, and X C is any argument. C C-------------------------------------------------------------------- INTEGER JINT CS REAL DOUBLE PRECISION 1 X, RESULT C-------------------------------------------------------------------- JINT=2 CALL CALCI0(X,RESULT,JINT) BESEI0=RESULT RETURN C---------- Last line of BESEI0 ---------- END SUBROUTINE CALCI1(ARG,RESULT,JINT) C-------------------------------------------------------------------- C C This packet computes modified Bessel functioons of the first kind C and order one, I1(X) and EXP(-ABS(X))*I1(X), for real C arguments X. It contains two function type subprograms, BESI1 C and BESEI1, and one subroutine type subprogram, CALCI1. C The calling statements for the primary entries are C C Y=BESI1(X) C and C Y=BESEI1(X) C C where the entry points correspond to the functions I1(X) and C EXP(-ABS(X))*I1(X), respectively. The routine CALCI1 is C intended for internal packet use only, all computations within C the packet being concentrated in this routine. The function C subprograms invoke CALCI1 with the statement C CALL CALCI1(ARG,RESULT,JINT) C where the parameter usage is as follows C C Function Parameters for CALCI1 C Call ARG RESULT JINT C C BESI1(ARG) ABS(ARG) .LE. XMAX I1(ARG) 1 C BESEI1(ARG) any real ARG EXP(-ABS(ARG))*I1(ARG) 2 C C The main computation evaluates slightly modified forms of C minimax approximations generated by Blair and Edwards, Chalk C River (Atomic Energy of Canada Limited) Report AECL-4928, C October, 1974. This transportable program is patterned after C the machine-dependent FUNPACK packet NATSI1, but cannot match C that version for efficiency or accuracy. This version uses C rational functions that theoretically approximate I-SUB-1(X) C to at least 18 significant decimal digits. The accuracy C achieved depends on the arithmetic system, the compiler, the C intrinsic functions, and proper selection of the machine- C dependent constants. C C******************************************************************* C******************************************************************* C C Explanation of machine-dependent constants. Let C C beta = Radix for the floating-point system C maxexp = Smallest power of beta that overflows C C Then the following machine-dependent constants must be declared C in DATA statements. IEEE values are provided as a default. C C XSMALL = Positive argument such that 1.0 - X = 1.0 to C machine precision for all ABS(X) .LE. XSMALL. C XINF = Largest positive machine number; approximately C beta**maxexp C XMAX = Largest argument acceptable to BESI1; Solution to C equation: C EXP(X) * (1-3/(8*X)) / SQRT(2*PI*X) = beta**maxexp C C C Approximate values for some important machines are: C C beta maxexp XSMALL C C CRAY-1 (S.P.) 2 8191 3.55E-15 C Cyber 180/855 C under NOS (S.P.) 2 1070 3.55E-15 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 2 128 2.98E-8 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 2 1024 5.55D-17 C IBM 3033 (D.P.) 16 63 6.95D-18 C VAX (S.P.) 2 127 2.98E-8 C VAX D-Format (D.P.) 2 127 6.95D-18 C VAX G-Format (D.P.) 2 1023 5.55D-17 C C C XINF XMAX C C CRAY-1 (S.P.) 5.45E+2465 5682.810 C Cyber 180/855 C under NOS (S.P.) 1.26E+322 745.894 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 3.40E+38 91.906 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 1.79D+308 713.987 C IBM 3033 (D.P.) 7.23D+75 178.185 C VAX (S.P.) 1.70D+38 91.209 C VAX D-Format (D.P.) 1.70D+38 91.209 C VAX G-Format (D.P.) 8.98D+307 713.293 C C******************************************************************* C******************************************************************* C C Error returns C C The program returns the value XINF for ABS(ARG) .GT. XMAX. C C C Intrinsic functions required are: C C ABS, SQRT, EXP C C C Authors: W. J. Cody and L. Stoltz C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C C Latest modification: March 13, 1992 C C-------------------------------------------------------------------- INTEGER J,JINT CS REAL DOUBLE PRECISION 1 A,ARG,B,EXP40,FORTY,HALF,ONE,ONE5,P,PBAR,PP,Q,QQ,REC15, 2 RESULT,SUMP,SUMQ,TWO25,X,XINF,XMAX,XSMALL,XX,ZERO DIMENSION P(15),PP(8),Q(5),QQ(6) C-------------------------------------------------------------------- C Mathematical constants C-------------------------------------------------------------------- CS DATA ONE/1.0E0/,ONE5/15.0E0/,EXP40/2.353852668370199854E17/, CS 1 FORTY/40.0E0/,REC15/6.6666666666666666666E-2/, CS 2 TWO25/225.0E0/,HALF/0.5E0/,ZERO/0.0E0/ DATA ONE/1.0D0/,ONE5/15.0D0/,EXP40/2.353852668370199854D17/, 1 FORTY/40.0D0/,REC15/6.6666666666666666666D-2/, 2 TWO25/225.0D0/,HALF/0.5D0/,ZERO/0.0D0/ C-------------------------------------------------------------------- C Machine-dependent constants C-------------------------------------------------------------------- CS DATA XSMALL/2.98E-8/,XINF/3.4E38/,XMAX/91.906E0/ DATA XSMALL/5.55D-17/,XINF/1.79D308/,XMAX/713.987D0/ C-------------------------------------------------------------------- C Coefficients for XSMALL .LE. ABS(ARG) .LT. 15.0 C-------------------------------------------------------------------- CS DATA P/-1.9705291802535139930E-19,-6.5245515583151902910E-16, CS 1 -1.1928788903603238754E-12,-1.4831904935994647675E-09, CS 2 -1.3466829827635152875E-06,-9.1746443287817501309E-04, CS 3 -4.7207090827310162436E-01,-1.8225946631657315931E+02, CS 4 -5.1894091982308017540E+04,-1.0588550724769347106E+07, CS 5 -1.4828267606612366099E+09,-1.3357437682275493024E+11, CS 6 -6.9876779648010090070E+12,-1.7732037840791591320E+14, CS 7 -1.4577180278143463643E+15/ DATA P/-1.9705291802535139930D-19,-6.5245515583151902910D-16, 1 -1.1928788903603238754D-12,-1.4831904935994647675D-09, 2 -1.3466829827635152875D-06,-9.1746443287817501309D-04, 3 -4.7207090827310162436D-01,-1.8225946631657315931D+02, 4 -5.1894091982308017540D+04,-1.0588550724769347106D+07, 5 -1.4828267606612366099D+09,-1.3357437682275493024D+11, 6 -6.9876779648010090070D+12,-1.7732037840791591320D+14, 7 -1.4577180278143463643D+15/ CS DATA Q/-4.0076864679904189921E+03, 7.4810580356655069138E+06, CS 1 -8.0059518998619764991E+09, 4.8544714258273622913E+12, CS 2 -1.3218168307321442305E+15/ DATA Q/-4.0076864679904189921D+03, 7.4810580356655069138D+06, 1 -8.0059518998619764991D+09, 4.8544714258273622913D+12, 2 -1.3218168307321442305D+15/ C-------------------------------------------------------------------- C Coefficients for 15.0 .LE. ABS(ARG) C-------------------------------------------------------------------- CS DATA PP/-6.0437159056137600000E-02, 4.5748122901933459000E-01, CS 1 -4.2843766903304806403E-01, 9.7356000150886612134E-02, CS 2 -3.2457723974465568321E-03,-3.6395264712121795296E-04, CS 3 1.6258661867440836395E-05,-3.6347578404608223492E-07/ DATA PP/-6.0437159056137600000D-02, 4.5748122901933459000D-01, 1 -4.2843766903304806403D-01, 9.7356000150886612134D-02, 2 -3.2457723974465568321D-03,-3.6395264712121795296D-04, 3 1.6258661867440836395D-05,-3.6347578404608223492D-07/ CS DATA QQ/-3.8806586721556593450E+00, 3.2593714889036996297E+00, CS 1 -8.5017476463217924408E-01, 7.4212010813186530069E-02, CS 2 -2.2835624489492512649E-03, 3.7510433111922824643E-05/ DATA QQ/-3.8806586721556593450D+00, 3.2593714889036996297D+00, 1 -8.5017476463217924408D-01, 7.4212010813186530069D-02, 2 -2.2835624489492512649D-03, 3.7510433111922824643D-05/ CS DATA PBAR/3.98437500E-01/ DATA PBAR/3.98437500D-01/ C-------------------------------------------------------------------- X = ABS(ARG) IF (X .LT. XSMALL) THEN C-------------------------------------------------------------------- C Return for ABS(ARG) .LT. XSMALL C-------------------------------------------------------------------- RESULT = HALF * X ELSE IF (X .LT. ONE5) THEN C-------------------------------------------------------------------- C XSMALL .LE. ABS(ARG) .LT. 15.0 C-------------------------------------------------------------------- XX = X * X SUMP = P(1) DO 50 J = 2, 15 SUMP = SUMP * XX + P(J) 50 CONTINUE XX = XX - TWO25 SUMQ = ((((XX+Q(1))*XX+Q(2))*XX+Q(3))*XX+Q(4)) 1 *XX+Q(5) RESULT = (SUMP / SUMQ) * X IF (JINT .EQ. 2) RESULT = RESULT * EXP(-X) ELSE IF ((JINT .EQ. 1) .AND. (X .GT. XMAX)) THEN RESULT = XINF ELSE C-------------------------------------------------------------------- C 15.0 .LE. ABS(ARG) C-------------------------------------------------------------------- XX = ONE / X - REC15 SUMP = ((((((PP(1)*XX+PP(2))*XX+PP(3))*XX+ 1 PP(4))*XX+PP(5))*XX+PP(6))*XX+PP(7))*XX+PP(8) SUMQ = (((((XX+QQ(1))*XX+QQ(2))*XX+QQ(3))*XX+ 1 QQ(4))*XX+QQ(5))*XX+QQ(6) RESULT = SUMP / SUMQ IF (JINT .NE. 1) THEN RESULT = (RESULT + PBAR) / SQRT(X) ELSE C-------------------------------------------------------------------- C Calculation reformulated to avoid premature overflow C-------------------------------------------------------------------- IF (X .GT. XMAX-ONE5) THEN A = EXP(X-FORTY) B = EXP40 ELSE A = EXP(X) B = ONE END IF RESULT = ((RESULT * A + PBAR * A) / 1 SQRT(X)) * B C-------------------------------------------------------------------- C Error return for ABS(ARG) .GT. XMAX C-------------------------------------------------------------------- END IF END IF IF (ARG .LT. ZERO) RESULT = -RESULT RETURN C----------- Last line of CALCI1 ----------- END CS REAL FUNCTION BESI1(X) DOUBLE PRECISION FUNCTION BESI1(X) C-------------------------------------------------------------------- C C This long precision subprogram computes approximate values for C modified Bessel functions of the first kind of order one for C arguments ABS(ARG) .LE. XMAX (see comments heading CALCI1). C C-------------------------------------------------------------------- INTEGER JINT CS REAL DOUBLE PRECISION 1 X, RESULT C-------------------------------------------------------------------- JINT=1 CALL CALCI1(X,RESULT,JINT) BESI1=RESULT RETURN C---------- Last line of BESI1 ---------- END CS REAL FUNCTION BESEI1(X) DOUBLE PRECISION FUNCTION BESEI1(X) C-------------------------------------------------------------------- C C This function program computes approximate values for the C modified Bessel function of the first kind of order one C multiplied by EXP(-ABS(X)), where EXP is the C exponential function, ABS is the absolute value, and X C is any argument. C C-------------------------------------------------------------------- INTEGER JINT CS REAL DOUBLE PRECISION 1 X, RESULT C-------------------------------------------------------------------- JINT=2 CALL CALCI1(X,RESULT,JINT) BESEI1=RESULT RETURN C---------- Last line of BESEI1 ---------- END SUBROUTINE CALJY0(ARG,RESULT,JINT) C--------------------------------------------------------------------- C C This packet computes zero-order Bessel functions of the first and C second kind (J0 and Y0), for real arguments X, where 0 < X <= XMAX C for Y0, and |X| <= XMAX for J0. It contains two function-type C subprograms, BESJ0 and BESY0, and one subroutine-type C subprogram, CALJY0. The calling statements for the primary C entries are: C C Y = BESJ0(X) C and C Y = BESY0(X), C C where the entry points correspond to the functions J0(X) and Y0(X), C respectively. The routine CALJY0 is intended for internal packet C use only, all computations within the packet being concentrated in C this one routine. The function subprograms invoke CALJY0 with C the statement C CALL CALJY0(ARG,RESULT,JINT), C where the parameter usage is as follows: C C Function Parameters for CALJY0 C call ARG RESULT JINT C C BESJ0(ARG) |ARG| .LE. XMAX J0(ARG) 0 C BESY0(ARG) 0 .LT. ARG .LE. XMAX Y0(ARG) 1 C C The main computation uses unpublished minimax rational C approximations for X .LE. 8.0, and an approximation from the C book Computer Approximations by Hart, et. al., Wiley and Sons, C New York, 1968, for arguments larger than 8.0 Part of this C transportable packet is patterned after the machine-dependent C FUNPACK program BESJ0(X), but cannot match that version for C efficiency or accuracy. This version uses rational functions C that are theoretically accurate to at least 18 significant decimal C digits for X <= 8, and at least 18 decimal places for X > 8. The C accuracy achieved depends on the arithmetic system, the compiler, C the intrinsic functions, and proper selection of the machine- C dependent constants. C C******************************************************************* C C The following machine-dependent constants must be declared in C DATA statements. IEEE values are provided as a default. C C XINF = largest positive machine number C XMAX = largest acceptable argument. The functions AINT, SIN C and COS must perform properly for ABS(X) .LE. XMAX. C We recommend that XMAX be a small integer multiple of C sqrt(1/eps), where eps is the smallest positive number C such that 1+eps > 1. C XSMALL = positive argument such that 1.0-(X/2)**2 = 1.0 C to machine precision for all ABS(X) .LE. XSMALL. C We recommend that XSMALL < sqrt(eps)/beta, where beta C is the floating-point radix (usually 2 or 16). C C Approximate values for some important machines are C C eps XMAX XSMALL XINF C C CDC 7600 (S.P.) 7.11E-15 1.34E+08 2.98E-08 1.26E+322 C CRAY-1 (S.P.) 7.11E-15 1.34E+08 2.98E-08 5.45E+2465 C IBM PC (8087) (S.P.) 5.96E-08 8.19E+03 1.22E-04 3.40E+38 C IBM PC (8087) (D.P.) 1.11D-16 2.68D+08 3.72D-09 1.79D+308 C IBM 195 (D.P.) 2.22D-16 6.87D+09 9.09D-13 7.23D+75 C UNIVAC 1108 (D.P.) 1.73D-18 4.30D+09 2.33D-10 8.98D+307 C VAX 11/780 (D.P.) 1.39D-17 1.07D+09 9.31D-10 1.70D+38 C C******************************************************************* C******************************************************************* C C Error Returns C C The program returns the value zero for X .GT. XMAX, and returns C -XINF when BESLY0 is called with a negative or zero argument. C C C Intrinsic functions required are: C C ABS, AINT, COS, LOG, SIN, SQRT C C C Latest modification: March 13, 1992 C C Author: W. J. Cody C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C C-------------------------------------------------------------------- INTEGER I,JINT CS REAL DOUBLE PRECISION 1 ARG,AX,CONS,DOWN,EIGHT,FIVE5,FOUR,ONE,ONEOV8,PI2,PJ0, 2 PJ1,PLG,PROD,PY0,PY1,PY2,P0,P1,P17,QJ0,QJ1,QLG,QY0,QY1, 3 QY2,Q0,Q1,RESJ,RESULT,R0,R1,SIXTY4,THREE,TWOPI,TWOPI1, 4 TWOPI2,TWO56,UP,W,WSQ,XDEN,XINF,XMAX,XNUM,XSMALL,XJ0, 5 XJ1,XJ01,XJ02,XJ11,XJ12,XY,XY0,XY01,XY02,XY1,XY11,XY12, 6 XY2,XY21,XY22,Z,ZERO,ZSQ DIMENSION PJ0(7),PJ1(8),PLG(4),PY0(6),PY1(7),PY2(8),P0(6),P1(6), 1 QJ0(5),QJ1(7),QLG(4),QY0(5),QY1(6),QY2(7),Q0(5),Q1(5) C------------------------------------------------------------------- C Mathematical constants C CONS = ln(.5) + Euler's gamma C------------------------------------------------------------------- CS DATA ZERO,ONE,THREE,FOUR,EIGHT/0.0E0,1.0E0,3.0E0,4.0E0,8.0E0/, CS 1 FIVE5,SIXTY4,ONEOV8,P17/5.5E0,64.0E0,0.125E0,1.716E-1/, CS 2 TWO56,CONS/256.0E0,-1.1593151565841244881E-1/, CS 3 PI2,TWOPI/6.3661977236758134308E-1,6.2831853071795864769E0/, CS 4 TWOPI1,TWOPI2/6.28125E0,1.9353071795864769253E-3/ DATA ZERO,ONE,THREE,FOUR,EIGHT/0.0D0,1.0D0,3.0D0,4.0D0,8.0D0/, 1 FIVE5,SIXTY4,ONEOV8,P17/5.5D0,64.0D0,0.125D0,1.716D-1/, 2 TWO56,CONS/256.0D0,-1.1593151565841244881D-1/, 3 PI2,TWOPI/6.3661977236758134308D-1,6.2831853071795864769D0/, 4 TWOPI1,TWOPI2/6.28125D0,1.9353071795864769253D-3/ C------------------------------------------------------------------- C Machine-dependent constants C------------------------------------------------------------------- CS DATA XMAX/8.19E+03/,XSMALL/1.22E-04/,XINF/3.40E+38/ DATA XMAX/2.68D+08/,XSMALL/3.72D-09/,XINF/1.79D+308/ C------------------------------------------------------------------- C Zeroes of Bessel functions C------------------------------------------------------------------- CS DATA XJ0/2.4048255576957727686E+0/,XJ1/5.5200781102863106496E+0/, CS 1 XY0/8.9357696627916752158E-1/,XY1/3.9576784193148578684E+0/, CS 2 XY2/7.0860510603017726976E+0/, CS 3 XJ01/ 616.0E+0/, XJ02/-1.4244423042272313784E-03/, CS 4 XJ11/1413.0E+0/, XJ12/ 5.4686028631064959660E-04/, CS 5 XY01/ 228.0E+0/, XY02/ 2.9519662791675215849E-03/, CS 6 XY11/1013.0E+0/, XY12/ 6.4716931485786837568E-04/, CS 7 XY21/1814.0E+0/, XY22/ 1.1356030177269762362E-04/ DATA XJ0/2.4048255576957727686D+0/,XJ1/5.5200781102863106496D+0/, 1 XY0/8.9357696627916752158D-1/,XY1/3.9576784193148578684D+0/, 2 XY2/7.0860510603017726976D+0/, 3 XJ01/ 616.0D+0/, XJ02/-1.4244423042272313784D-03/, 4 XJ11/1413.0D+0/, XJ12/ 5.4686028631064959660D-04/, 5 XY01/ 228.0D+0/, XY02/ 2.9519662791675215849D-03/, 6 XY11/1013.0D+0/, XY12/ 6.4716931485786837568D-04/, 7 XY21/1814.0D+0/, XY22/ 1.1356030177269762362D-04/ C------------------------------------------------------------------- C Coefficients for rational approximation to ln(x/a) C-------------------------------------------------------------------- CS DATA PLG/-2.4562334077563243311E+01,2.3642701335621505212E+02, CS 1 -5.4989956895857911039E+02,3.5687548468071500413E+02/ CS DATA QLG/-3.5553900764052419184E+01,1.9400230218539473193E+02, CS 1 -3.3442903192607538956E+02,1.7843774234035750207E+02/ DATA PLG/-2.4562334077563243311D+01,2.3642701335621505212D+02, 1 -5.4989956895857911039D+02,3.5687548468071500413D+02/ DATA QLG/-3.5553900764052419184D+01,1.9400230218539473193D+02, 1 -3.3442903192607538956D+02,1.7843774234035750207D+02/ C------------------------------------------------------------------- C Coefficients for rational approximation of C J0(X) / (X**2 - XJ0**2), XSMALL < |X| <= 4.0 C-------------------------------------------------------------------- CS DATA PJ0/6.6302997904833794242E+06,-6.2140700423540120665E+08, CS 1 2.7282507878605942706E+10,-4.1298668500990866786E+11, CS 2 -1.2117036164593528341E-01, 1.0344222815443188943E+02, CS 3 -3.6629814655107086448E+04/ CS DATA QJ0/4.5612696224219938200E+05, 1.3985097372263433271E+08, CS 1 2.6328198300859648632E+10, 2.3883787996332290397E+12, CS 2 9.3614022392337710626E+02/ DATA PJ0/6.6302997904833794242D+06,-6.2140700423540120665D+08, 1 2.7282507878605942706D+10,-4.1298668500990866786D+11, 2 -1.2117036164593528341D-01, 1.0344222815443188943D+02, 3 -3.6629814655107086448D+04/ DATA QJ0/4.5612696224219938200D+05, 1.3985097372263433271D+08, 1 2.6328198300859648632D+10, 2.3883787996332290397D+12, 2 9.3614022392337710626D+02/ C------------------------------------------------------------------- C Coefficients for rational approximation of C J0(X) / (X**2 - XJ1**2), 4.0 < |X| <= 8.0 C------------------------------------------------------------------- CS DATA PJ1/4.4176707025325087628E+03, 1.1725046279757103576E+04, CS 1 1.0341910641583726701E+04,-7.2879702464464618998E+03, CS 2 -1.2254078161378989535E+04,-1.8319397969392084011E+03, CS 3 4.8591703355916499363E+01, 7.4321196680624245801E+02/ CS DATA QJ1/3.3307310774649071172E+02,-2.9458766545509337327E+03, CS 1 1.8680990008359188352E+04,-8.4055062591169562211E+04, CS 2 2.4599102262586308984E+05,-3.5783478026152301072E+05, CS 3 -2.5258076240801555057E+01/ DATA PJ1/4.4176707025325087628D+03, 1.1725046279757103576D+04, 1 1.0341910641583726701D+04,-7.2879702464464618998D+03, 2 -1.2254078161378989535D+04,-1.8319397969392084011D+03, 3 4.8591703355916499363D+01, 7.4321196680624245801D+02/ DATA QJ1/3.3307310774649071172D+02,-2.9458766545509337327D+03, 1 1.8680990008359188352D+04,-8.4055062591169562211D+04, 2 2.4599102262586308984D+05,-3.5783478026152301072D+05, 3 -2.5258076240801555057D+01/ C------------------------------------------------------------------- C Coefficients for rational approximation of C (Y0(X) - 2 LN(X/XY0) J0(X)) / (X**2 - XY0**2), C XSMALL < |X| <= 3.0 C-------------------------------------------------------------------- CS DATA PY0/1.0102532948020907590E+04,-2.1287548474401797963E+06, CS 1 2.0422274357376619816E+08,-8.3716255451260504098E+09, CS 2 1.0723538782003176831E+11,-1.8402381979244993524E+01/ CS DATA QY0/6.6475986689240190091E+02, 2.3889393209447253406E+05, CS 1 5.5662956624278251596E+07, 8.1617187777290363573E+09, CS 2 5.8873865738997033405E+11/ DATA PY0/1.0102532948020907590D+04,-2.1287548474401797963D+06, 1 2.0422274357376619816D+08,-8.3716255451260504098D+09, 2 1.0723538782003176831D+11,-1.8402381979244993524D+01/ DATA QY0/6.6475986689240190091D+02, 2.3889393209447253406D+05, 1 5.5662956624278251596D+07, 8.1617187777290363573D+09, 2 5.8873865738997033405D+11/ C------------------------------------------------------------------- C Coefficients for rational approximation of C (Y0(X) - 2 LN(X/XY1) J0(X)) / (X**2 - XY1**2), C 3.0 < |X| <= 5.5 C-------------------------------------------------------------------- CS DATA PY1/-1.4566865832663635920E+04, 4.6905288611678631510E+06, CS 1 -6.9590439394619619534E+08, 4.3600098638603061642E+10, CS 2 -5.5107435206722644429E+11,-2.2213976967566192242E+13, CS 3 1.7427031242901594547E+01/ CS DATA QY1/ 8.3030857612070288823E+02, 4.0669982352539552018E+05, CS 1 1.3960202770986831075E+08, 3.4015103849971240096E+10, CS 2 5.4266824419412347550E+12, 4.3386146580707264428E+14/ DATA PY1/-1.4566865832663635920D+04, 4.6905288611678631510D+06, 1 -6.9590439394619619534D+08, 4.3600098638603061642D+10, 2 -5.5107435206722644429D+11,-2.2213976967566192242D+13, 3 1.7427031242901594547D+01/ DATA QY1/ 8.3030857612070288823D+02, 4.0669982352539552018D+05, 1 1.3960202770986831075D+08, 3.4015103849971240096D+10, 2 5.4266824419412347550D+12, 4.3386146580707264428D+14/ C------------------------------------------------------------------- C Coefficients for rational approximation of C (Y0(X) - 2 LN(X/XY2) J0(X)) / (X**2 - XY2**2), C 5.5 < |X| <= 8.0 C-------------------------------------------------------------------- CS DATA PY2/ 2.1363534169313901632E+04,-1.0085539923498211426E+07, CS 1 2.1958827170518100757E+09,-1.9363051266772083678E+11, CS 2 -1.2829912364088687306E+11, 6.7016641869173237784E+14, CS 3 -8.0728726905150210443E+15,-1.7439661319197499338E+01/ CS DATA QY2/ 8.7903362168128450017E+02, 5.3924739209768057030E+05, CS 1 2.4727219475672302327E+08, 8.6926121104209825246E+10, CS 2 2.2598377924042897629E+13, 3.9272425569640309819E+15, CS 3 3.4563724628846457519E+17/ DATA PY2/ 2.1363534169313901632D+04,-1.0085539923498211426D+07, 1 2.1958827170518100757D+09,-1.9363051266772083678D+11, 2 -1.2829912364088687306D+11, 6.7016641869173237784D+14, 3 -8.0728726905150210443D+15,-1.7439661319197499338D+01/ DATA QY2/ 8.7903362168128450017D+02, 5.3924739209768057030D+05, 1 2.4727219475672302327D+08, 8.6926121104209825246D+10, 2 2.2598377924042897629D+13, 3.9272425569640309819D+15, 3 3.4563724628846457519D+17/ C------------------------------------------------------------------- C Coefficients for Hart,s approximation, |X| > 8.0 C------------------------------------------------------------------- CS DATA P0/3.4806486443249270347E+03, 2.1170523380864944322E+04, CS 1 4.1345386639580765797E+04, 2.2779090197304684302E+04, CS 2 8.8961548424210455236E-01, 1.5376201909008354296E+02/ CS DATA Q0/3.5028735138235608207E+03, 2.1215350561880115730E+04, CS 1 4.1370412495510416640E+04, 2.2779090197304684318E+04, CS 2 1.5711159858080893649E+02/ CS DATA P1/-2.2300261666214198472E+01,-1.1183429920482737611E+02, CS 1 -1.8591953644342993800E+02,-8.9226600200800094098E+01, CS 2 -8.8033303048680751817E-03,-1.2441026745835638459E+00/ CS DATA Q1/1.4887231232283756582E+03, 7.2642780169211018836E+03, CS 1 1.1951131543434613647E+04, 5.7105024128512061905E+03, CS 2 9.0593769594993125859E+01/ DATA P0/3.4806486443249270347D+03, 2.1170523380864944322D+04, 1 4.1345386639580765797D+04, 2.2779090197304684302D+04, 2 8.8961548424210455236D-01, 1.5376201909008354296D+02/ DATA Q0/3.5028735138235608207D+03, 2.1215350561880115730D+04, 1 4.1370412495510416640D+04, 2.2779090197304684318D+04, 2 1.5711159858080893649D+02/ DATA P1/-2.2300261666214198472D+01,-1.1183429920482737611D+02, 1 -1.8591953644342993800D+02,-8.9226600200800094098D+01, 2 -8.8033303048680751817D-03,-1.2441026745835638459D+00/ DATA Q1/1.4887231232283756582D+03, 7.2642780169211018836D+03, 1 1.1951131543434613647D+04, 5.7105024128512061905D+03, 2 9.0593769594993125859D+01/ C------------------------------------------------------------------- C Check for error conditions C------------------------------------------------------------------- AX = ABS(ARG) IF ((JINT .EQ. 1) .AND. (ARG .LE. ZERO)) THEN RESULT = -XINF GO TO 2000 ELSE IF (AX .GT. XMAX) THEN RESULT = ZERO GO TO 2000 END IF IF (AX .GT. EIGHT) GO TO 800 IF (AX .LE. XSMALL) THEN IF (JINT .EQ. 0) THEN RESULT = ONE ELSE RESULT = PI2 * (LOG(AX) + CONS) END IF GO TO 2000 END IF C------------------------------------------------------------------- C Calculate J0 for appropriate interval, preserving C accuracy near the zero of J0 C------------------------------------------------------------------- ZSQ = AX * AX IF (AX .LE. FOUR) THEN XNUM = (PJ0(5) * ZSQ + PJ0(6)) * ZSQ + PJ0(7) XDEN = ZSQ + QJ0(5) DO 50 I = 1, 4 XNUM = XNUM * ZSQ + PJ0(I) XDEN = XDEN * ZSQ + QJ0(I) 50 CONTINUE PROD = ((AX - XJ01/TWO56) - XJ02) * (AX + XJ0) ELSE WSQ = ONE - ZSQ / SIXTY4 XNUM = PJ1(7) * WSQ + PJ1(8) XDEN = WSQ + QJ1(7) DO 220 I = 1, 6 XNUM = XNUM * WSQ + PJ1(I) XDEN = XDEN * WSQ + QJ1(I) 220 CONTINUE PROD = (AX + XJ1) * ((AX - XJ11/TWO56) - XJ12) END IF RESULT = PROD * XNUM / XDEN IF (JINT .EQ. 0) GO TO 2000 C------------------------------------------------------------------- C Calculate Y0. First find RESJ = pi/2 ln(x/xn) J0(x), C where xn is a zero of Y0 C------------------------------------------------------------------- IF (AX .LE. THREE) THEN UP = (AX-XY01/TWO56)-XY02 XY = XY0 ELSE IF (AX .LE. FIVE5) THEN UP = (AX-XY11/TWO56)-XY12 XY = XY1 ELSE UP = (AX-XY21/TWO56)-XY22 XY = XY2 END IF DOWN = AX + XY IF (ABS(UP) .LT. P17*DOWN) THEN W = UP/DOWN WSQ = W*W XNUM = PLG(1) XDEN = WSQ + QLG(1) DO 320 I = 2, 4 XNUM = XNUM*WSQ + PLG(I) XDEN = XDEN*WSQ + QLG(I) 320 CONTINUE RESJ = PI2 * RESULT * W * XNUM/XDEN ELSE RESJ = PI2 * RESULT * LOG(AX/XY) END IF C------------------------------------------------------------------- C Now calculate Y0 for appropriate interval, preserving C accuracy near the zero of Y0 C------------------------------------------------------------------- IF (AX .LE. THREE) THEN XNUM = PY0(6) * ZSQ + PY0(1) XDEN = ZSQ + QY0(1) DO 340 I = 2, 5 XNUM = XNUM * ZSQ + PY0(I) XDEN = XDEN * ZSQ + QY0(I) 340 CONTINUE ELSE IF (AX .LE. FIVE5) THEN XNUM = PY1(7) * ZSQ + PY1(1) XDEN = ZSQ + QY1(1) DO 360 I = 2, 6 XNUM = XNUM * ZSQ + PY1(I) XDEN = XDEN * ZSQ + QY1(I) 360 CONTINUE ELSE XNUM = PY2(8) * ZSQ + PY2(1) XDEN = ZSQ + QY2(1) DO 380 I = 2, 7 XNUM = XNUM * ZSQ + PY2(I) XDEN = XDEN * ZSQ + QY2(I) 380 CONTINUE END IF RESULT = RESJ + UP * DOWN * XNUM / XDEN GO TO 2000 C------------------------------------------------------------------- C Calculate J0 or Y0 for |ARG| > 8.0 C------------------------------------------------------------------- 800 Z = EIGHT / AX W = AX / TWOPI W = AINT(W) + ONEOV8 W = (AX - W * TWOPI1) - W * TWOPI2 ZSQ = Z * Z XNUM = P0(5) * ZSQ + P0(6) XDEN = ZSQ + Q0(5) UP = P1(5) * ZSQ + P1(6) DOWN = ZSQ + Q1(5) DO 850 I = 1, 4 XNUM = XNUM * ZSQ + P0(I) XDEN = XDEN * ZSQ + Q0(I) UP = UP * ZSQ + P1(I) DOWN = DOWN * ZSQ + Q1(I) 850 CONTINUE R0 = XNUM / XDEN R1 = UP / DOWN IF (JINT .EQ. 0) THEN RESULT = SQRT(PI2/AX) * (R0*COS(W) - Z*R1*SIN(W)) ELSE RESULT = SQRT(PI2/AX) * (R0*SIN(W) + Z*R1*COS(W)) END IF 2000 RETURN C---------- Last line of CALJY0 ---------- END CS REAL FUNCTION BESJ0(X) DOUBLE PRECISION FUNCTION BESJ0(X) C-------------------------------------------------------------------- C C This subprogram computes approximate values for Bessel functions C of the first kind of order zero for arguments |X| <= XMAX C (see comments heading CALJY0). C C-------------------------------------------------------------------- INTEGER JINT CS REAL X, RESULT DOUBLE PRECISION X, RESULT C-------------------------------------------------------------------- JINT=0 CALL CALJY0(X,RESULT,JINT) BESJ0 = RESULT RETURN C---------- Last line of BESJ0 ---------- END CS REAL FUNCTION BESY0(X) DOUBLE PRECISION FUNCTION BESY0(X) C-------------------------------------------------------------------- C C This subprogram computes approximate values for Bessel functions C of the second kind of order zero for arguments 0 < X <= XMAX C (see comments heading CALJY0). C C-------------------------------------------------------------------- INTEGER JINT CS REAL X, RESULT DOUBLE PRECISION X, RESULT C-------------------------------------------------------------------- JINT=1 CALL CALJY0(X,RESULT,JINT) BESY0 = RESULT RETURN C---------- Last line of BESY0 ---------- END SUBROUTINE CALJY1(ARG,RESULT,JINT) C--------------------------------------------------------------------- C C This packet computes first-order Bessel functions of the first and C second kind (J1 and Y1), for real arguments X, where 0 < X <= XMAX C for Y1, and |X| <= XMAX for J1. It contains two function-type C subprograms, BESJ1 and BESY1, and one subroutine-type C subprogram, CALJY1. The calling statements for the primary C entries are: C C Y = BESJ1(X) C and C Y = BESY1(X), C C where the entry points correspond to the functions J1(X) and Y1(X), C respectively. The routine CALJY1 is intended for internal packet C use only, all computations within the packet being concentrated in C this one routine. The function subprograms invoke CALJY1 with C the statement C CALL CALJY1(ARG,RESULT,JINT), C where the parameter usage is as follows: C C Function Parameters for CALJY1 C call ARG RESULT JINT C C BESJ1(ARG) |ARG| .LE. XMAX J1(ARG) 0 C BESY1(ARG) 0 .LT. ARG .LE. XMAX Y1(ARG) 1 C C The main computation uses unpublished minimax rational C approximations for X .LE. 8.0, and an approximation from the C book Computer Approximations by Hart, et. al., Wiley and Sons, C New York, 1968, for arguments larger than 8.0 Part of this C transportable packet is patterned after the machine-dependent C FUNPACK program BESJ1(X), but cannot match that version for C efficiency or accuracy. This version uses rational functions C that are theoretically accurate to at least 18 significant decimal C digits for X <= 8, and at least 18 decimal places for X > 8. The C accuracy achieved depends on the arithmetic system, the compiler, C the intrinsic functions, and proper selection of the machine- C dependent constants. C C******************************************************************* C C The following machine-dependent constants must be declared in C DATA statements. IEEE values are provided as a default. C C XINF = largest positive machine number C XMAX = largest acceptable argument. The functions AINT, SIN C and COS must perform properly for ABS(X) .LE. XMAX. C We recommend that XMAX be a small integer multiple of C sqrt(1/eps), where eps is the smallest positive number C such that 1+eps > 1. C XSMALL = positive argument such that 1.0-(1/2)(X/2)**2 = 1.0 C to machine precision for all ABS(X) .LE. XSMALL. C We recommend that XSMALL < sqrt(eps)/beta, where beta C is the floating-point radix (usually 2 or 16). C C Approximate values for some important machines are C C eps XMAX XSMALL XINF C C CDC 7600 (S.P.) 7.11E-15 1.34E+08 2.98E-08 1.26E+322 C CRAY-1 (S.P.) 7.11E-15 1.34E+08 2.98E-08 5.45E+2465 C IBM PC (8087) (S.P.) 5.96E-08 8.19E+03 1.22E-04 3.40E+38 C IBM PC (8087) (D.P.) 1.11D-16 2.68D+08 3.72D-09 1.79D+308 C IBM 195 (D.P.) 2.22D-16 6.87D+09 9.09D-13 7.23D+75 C UNIVAC 1108 (D.P.) 1.73D-18 4.30D+09 2.33D-10 8.98D+307 C VAX 11/780 (D.P.) 1.39D-17 1.07D+09 9.31D-10 1.70D+38 C C******************************************************************* C******************************************************************* C C Error Returns C C The program returns the value zero for X .GT. XMAX, and returns C -XINF when BESLY1 is called with a negative or zero argument. C C C Intrinsic functions required are: C C ABS, AINT, COS, LOG, SIN, SQRT C C C Author: W. J. Cody C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C C Latest modification: March 13, 1992 C C-------------------------------------------------------------------- INTEGER I,JINT DIMENSION PJ0(7),PJ1(8),PLG(4),PY0(7),PY1(9),P0(6),P1(6), 1 QJ0(5),QJ1(7),QLG(4),QY0(6),QY1(8),Q0(6),Q1(6) CS REAL DOUBLE PRECISION 1 ARG,AX,DOWN,EIGHT,FOUR,HALF,PI2,PJ0,PJ1,PLG,PROD,PY0, 2 PY1,P0,P1,P17,QJ0,QJ1,QLG,QY0,QY1,Q0,Q1,RESJ,RESULT, 3 RTPI2,R0,R1,THROV8,TWOPI,TWOPI1,TWOPI2,TWO56,UP,W,WSQ, 4 XDEN,XINF,XMAX,XNUM,XSMALL,XJ0,XJ1,XJ01,XJ02,XJ11,XJ12, 5 XY,XY0,XY01,XY02,XY1,XY11,XY12,Z,ZERO,ZSQ C------------------------------------------------------------------- C Mathematical constants C------------------------------------------------------------------- CS DATA EIGHT/8.0E0/, CS 1 FOUR/4.0E0/,HALF/0.5E0/,THROV8/0.375E0/, CS 2 PI2/6.3661977236758134308E-1/,P17/1.716E-1/ CS 3 TWOPI/6.2831853071795864769E+0/,ZERO/0.0E0/, CS 4 TWOPI1/6.28125E0/,TWOPI2/1.9353071795864769253E-03/ CS 5 TWO56/256.0E+0/,RTPI2/7.9788456080286535588E-1/ DATA EIGHT/8.0D0/, 1 FOUR/4.0D0/,HALF/0.5D0/,THROV8/0.375D0/, 2 PI2/6.3661977236758134308D-1/,P17/1.716D-1/ 3 TWOPI/6.2831853071795864769D+0/,ZERO/0.0D0/, 4 TWOPI1/6.28125D0/,TWOPI2/1.9353071795864769253D-03/ 5 TWO56/256.0D+0/,RTPI2/7.9788456080286535588D-1/ C------------------------------------------------------------------- C Machine-dependent constants C------------------------------------------------------------------- CS DATA XMAX/8.19E+03/,XSMALL/1.22E-04/,XINF/3.40E+38/ DATA XMAX/2.68D+08/,XSMALL/3.72D-09/,XINF/1.79D+308/ C------------------------------------------------------------------- C Zeroes of Bessel functions C------------------------------------------------------------------- CS DATA XJ0/3.8317059702075123156E+0/,XJ1/7.0155866698156187535E+0/, CS 1 XY0/2.1971413260310170351E+0/,XY1/5.4296810407941351328E+0/, CS 2 XJ01/ 981.0E+0/, XJ02/-3.2527979248768438556E-04/, CS 3 XJ11/1796.0E+0/, XJ12/-3.8330184381246462950E-05/, CS 4 XY01/ 562.0E+0/, XY02/ 1.8288260310170351490E-03/, CS 5 XY11/1390.0E+0/, XY12/-6.4592058648672279948E-06/ DATA XJ0/3.8317059702075123156D+0/,XJ1/7.0155866698156187535D+0/, 1 XY0/2.1971413260310170351D+0/,XY1/5.4296810407941351328D+0/, 2 XJ01/ 981.0D+0/, XJ02/-3.2527979248768438556D-04/, 3 XJ11/1796.0D+0/, XJ12/-3.8330184381246462950D-05/, 4 XY01/ 562.0D+0/, XY02/ 1.8288260310170351490D-03/, 5 XY11/1390.0D+0/, XY12/-6.4592058648672279948D-06/ C------------------------------------------------------------------- C Coefficients for rational approximation to ln(x/a) C-------------------------------------------------------------------- CS DATA PLG/-2.4562334077563243311E+01,2.3642701335621505212E+02, CS 1 -5.4989956895857911039E+02,3.5687548468071500413E+02/ CS DATA QLG/-3.5553900764052419184E+01,1.9400230218539473193E+02, CS 1 -3.3442903192607538956E+02,1.7843774234035750207E+02/ DATA PLG/-2.4562334077563243311D+01,2.3642701335621505212D+02, 1 -5.4989956895857911039D+02,3.5687548468071500413D+02/ DATA QLG/-3.5553900764052419184D+01,1.9400230218539473193D+02, 1 -3.3442903192607538956D+02,1.7843774234035750207D+02/ C------------------------------------------------------------------- C Coefficients for rational approximation of C J1(X) / (X * (X**2 - XJ0**2)), XSMALL < |X| <= 4.0 C-------------------------------------------------------------------- CS DATA PJ0/9.8062904098958257677E+05,-1.1548696764841276794E+08, CS 1 6.6781041261492395835E+09,-1.4258509801366645672E+11, CS 2 -4.4615792982775076130E+03, 1.0650724020080236441E+01, CS 3 -1.0767857011487300348E-02/ CS DATA QJ0/5.9117614494174794095E+05, 2.0228375140097033958E+08, CS 1 4.2091902282580133541E+10, 4.1868604460820175290E+12, CS 2 1.0742272239517380498E+03/ DATA PJ0/9.8062904098958257677D+05,-1.1548696764841276794D+08, 1 6.6781041261492395835D+09,-1.4258509801366645672D+11, 2 -4.4615792982775076130D+03, 1.0650724020080236441D+01, 3 -1.0767857011487300348D-02/ DATA QJ0/5.9117614494174794095D+05, 2.0228375140097033958D+08, 1 4.2091902282580133541D+10, 4.1868604460820175290D+12, 2 1.0742272239517380498D+03/ C------------------------------------------------------------------- C Coefficients for rational approximation of C J1(X) / (X * (X**2 - XJ1**2)), 4.0 < |X| <= 8.0 C------------------------------------------------------------------- CS DATA PJ1/4.6179191852758252280E+00,-7.1329006872560947377E+03, CS 1 4.5039658105749078904E+06,-1.4437717718363239107E+09, CS 2 2.3569285397217157313E+11,-1.6324168293282543629E+13, CS 3 1.1357022719979468624E+14, 1.0051899717115285432E+15/ CS DATA QJ1/1.1267125065029138050E+06, 6.4872502899596389593E+08, CS 1 2.7622777286244082666E+11, 8.4899346165481429307E+13, CS 2 1.7128800897135812012E+16, 1.7253905888447681194E+18, CS 3 1.3886978985861357615E+03/ DATA PJ1/4.6179191852758252280D+00,-7.1329006872560947377D+03, 1 4.5039658105749078904D+06,-1.4437717718363239107D+09, 2 2.3569285397217157313D+11,-1.6324168293282543629D+13, 3 1.1357022719979468624D+14, 1.0051899717115285432D+15/ DATA QJ1/1.1267125065029138050D+06, 6.4872502899596389593D+08, 1 2.7622777286244082666D+11, 8.4899346165481429307D+13, 2 1.7128800897135812012D+16, 1.7253905888447681194D+18, 3 1.3886978985861357615D+03/ C------------------------------------------------------------------- C Coefficients for rational approximation of C (Y1(X) - 2 LN(X/XY0) J1(X)) / (X**2 - XY0**2), C XSMALL < |X| <= 4.0 C-------------------------------------------------------------------- CS DATA PY0/2.2157953222280260820E+05,-5.9157479997408395984E+07, CS 1 7.2144548214502560419E+09,-3.7595974497819597599E+11, CS 2 5.4708611716525426053E+12, 4.0535726612579544093E+13, CS 3 -3.1714424660046133456E+02/ CS DATA QY0/8.2079908168393867438E+02, 3.8136470753052572164E+05, CS 1 1.2250435122182963220E+08, 2.7800352738690585613E+10, CS 2 4.1272286200406461981E+12, 3.0737873921079286084E+14/ DATA PY0/2.2157953222280260820D+05,-5.9157479997408395984D+07, 1 7.2144548214502560419D+09,-3.7595974497819597599D+11, 2 5.4708611716525426053D+12, 4.0535726612579544093D+13, 3 -3.1714424660046133456D+02/ DATA QY0/8.2079908168393867438D+02, 3.8136470753052572164D+05, 1 1.2250435122182963220D+08, 2.7800352738690585613D+10, 2 4.1272286200406461981D+12, 3.0737873921079286084D+14/ C-------------------------------------------------------------------- C Coefficients for rational approximation of C (Y1(X) - 2 LN(X/XY1) J1(X)) / (X**2 - XY1**2), C 4.0 < |X| <= 8.0 C-------------------------------------------------------------------- CS DATA PY1/ 1.9153806858264202986E+06,-1.1957961912070617006E+09, CS 1 3.7453673962438488783E+11,-5.9530713129741981618E+13, CS 2 4.0686275289804744814E+15,-2.3638408497043134724E+16, CS 3 -5.6808094574724204577E+18, 1.1514276357909013326E+19, CS 4 -1.2337180442012953128E+03/ CS DATA QY1/ 1.2855164849321609336E+03, 1.0453748201934079734E+06, CS 1 6.3550318087088919566E+08, 3.0221766852960403645E+11, CS 2 1.1187010065856971027E+14, 3.0837179548112881950E+16, CS 3 5.6968198822857178911E+18, 5.3321844313316185697E+20/ DATA PY1/ 1.9153806858264202986D+06,-1.1957961912070617006D+09, 1 3.7453673962438488783D+11,-5.9530713129741981618D+13, 2 4.0686275289804744814D+15,-2.3638408497043134724D+16, 3 -5.6808094574724204577D+18, 1.1514276357909013326D+19, 4 -1.2337180442012953128D+03/ DATA QY1/ 1.2855164849321609336D+03, 1.0453748201934079734D+06, 1 6.3550318087088919566D+08, 3.0221766852960403645D+11, 2 1.1187010065856971027D+14, 3.0837179548112881950D+16, 3 5.6968198822857178911D+18, 5.3321844313316185697D+20/ C------------------------------------------------------------------- C Coefficients for Hart,s approximation, |X| > 8.0 C------------------------------------------------------------------- CS DATA P0/-1.0982405543459346727E+05,-1.5235293511811373833E+06, CS 1 -6.6033732483649391093E+06,-9.9422465050776411957E+06, CS 2 -4.4357578167941278571E+06,-1.6116166443246101165E+03/ CS DATA Q0/-1.0726385991103820119E+05,-1.5118095066341608816E+06, CS 1 -6.5853394797230870728E+06,-9.9341243899345856590E+06, CS 2 -4.4357578167941278568E+06,-1.4550094401904961825E+03/ CS DATA P1/ 1.7063754290207680021E+03, 1.8494262873223866797E+04, CS 1 6.6178836581270835179E+04, 8.5145160675335701966E+04, CS 2 3.3220913409857223519E+04, 3.5265133846636032186E+01/ CS DATA Q1/ 3.7890229745772202641E+04, 4.0029443582266975117E+05, CS 1 1.4194606696037208929E+06, 1.8194580422439972989E+06, CS 2 7.0871281941028743574E+05, 8.6383677696049909675E+02/ DATA P0/-1.0982405543459346727D+05,-1.5235293511811373833D+06, 1 -6.6033732483649391093D+06,-9.9422465050776411957D+06, 2 -4.4357578167941278571D+06,-1.6116166443246101165D+03/ DATA Q0/-1.0726385991103820119D+05,-1.5118095066341608816D+06, 1 -6.5853394797230870728D+06,-9.9341243899345856590D+06, 2 -4.4357578167941278568D+06,-1.4550094401904961825D+03/ DATA P1/ 1.7063754290207680021D+03, 1.8494262873223866797D+04, 1 6.6178836581270835179D+04, 8.5145160675335701966D+04, 2 3.3220913409857223519D+04, 3.5265133846636032186D+01/ DATA Q1/ 3.7890229745772202641D+04, 4.0029443582266975117D+05, 1 1.4194606696037208929D+06, 1.8194580422439972989D+06, 2 7.0871281941028743574D+05, 8.6383677696049909675D+02/ C------------------------------------------------------------------- C Check for error conditions C------------------------------------------------------------------- AX = ABS(ARG) IF ((JINT .EQ. 1) .AND. ((ARG .LE. ZERO) .OR. 1 ((ARG .LT. HALF) .AND. (AX*XINF .LT. PI2)))) THEN RESULT = -XINF GO TO 2000 ELSE IF (AX .GT. XMAX) THEN RESULT = Z