SUBROUTINE RIBESL(X, ALPHA, NB, IZE, B, NCALC) c*********************************************************************72 c cc RIBESL computes Bessel I functions with real (non integer) order. c C THIS ROUTINE CALCULATES BESSEL FUNCTIONS I SUB(N+ALPHA) (X) C FOR NON-NEGATIVE ARGUMENT X, AND NON-NEGATIVE ORDER N+ALPHA, C WITH OR WITHOUT EXPONENTIAL SCALING. C C C EXPLANATION OF VARIABLES IN THE CALLING SEQUENCE C C X - WORKING PRECISION NON-NEGATIVE REAL ARGUMENT FOR WHICH C I'S OR EXPONENTIALLY SCALED I'S (I*EXP(-X)) C ARE TO BE CALCULATED. IF I'S ARE TO BE CALCULATED, C X MUST BE LESS THAN EXPARG (SEE BELOW). C ALPHA - WORKING PRECISION FRACTIONAL PART OF ORDER FOR WHICH C I'S OR EXPONENTIALLY SCALED I'S (I*EXP(-X)) ARE C TO BE CALCULATED. 0 .LE. ALPHA .LT. 1.0. C NB - INTEGER NUMBER OF FUNCTIONS TO BE CALCULATED, NB .GT. 0. C THE FIRST FUNCTION CALCULATED IS OF ORDER ALPHA, AND THE C LAST IS OF ORDER (NB - 1 + ALPHA). C IZE - INTEGER TYPE. IZE = 1 IF UNSCALED I'S ARE TO CALCULATED, C AND 2 IF EXPONENTIALLY SCALED I'S ARE TO BE CALCULATED. C B - WORKING PRECISION OUTPUT VECTOR OF LENGTH NB. IF THE ROUTINE C TERMINATES NORMALLY (NCALC=NB), THE VECTOR B CONTAINS THE C FUNCTIONS I/ALPHA/(X) THROUGH I/NB-1+ALPHA/(X), OR THE C CORRESPONDING EXPONENTIALLY SCALED FUNCTIONS. C NCALC - INTEGER OUTPUT VARIABLE INDICATING POSSIBLE ERRORS. C BEFORE USING THE VECTOR B, THE USER SHOULD CHECK THAT C NCALC=NB, I.E., ALL ORDERS HAVE BEEN CALCULATED TO C THE DESIRED ACCURACY. SEE ERROR RETURNS BELOW. C C C******************************************************************* C******************************************************************* C C EXPLANATION OF MACHINE-DEPENDENT CONSTANTS C C NSIG - DECIMAL SIGNIFICANCE DESIRED. SHOULD BE SET TO C IFIX(ALOG10(2)*NBIT+1), WHERE NBIT IS THE NUMBER OF C BITS IN THE MANTISSA OF A WORKING PRECISION VARIABLE. C SETTING NSIG LOWER WILL RESULT IN DECREASED ACCURACY C WHILE SETTING NSIG HIGHER WILL INCREASE CPU TIME C WITHOUT INCREASING ACCURACY. THE TRUNCATION ERROR C IS LIMITED TO A RELATIVE ERROR OF T=.5*10**(-NSIG). C ENTEN - 10.0 ** K, WHERE K IS THE LARGEST INTEGER SUCH THAT C ENTEN IS MACHINE-REPRESENTABLE IN WORKING PRECISION. C ENSIG - 10.0 ** NSIG. C RTNSIG - 10.0 ** (-K) FOR THE SMALLEST INTEGER K SUCH THAT C K .GE. NSIG/4. C ENMTEN - THE SMALLEST ABS(X) SUCH THAT X/4 DOES NOT UNDERFLOW. C XLARGE - UPPER LIMIT ON THE MAGNITUDE OF X WHEN IZE=2. BEAR C IN MIND THAT IF ABS(X)=N, THEN AT LEAST N ITERATIONS C OF THE BACKWARD RECURSION WILL BE EXECUTED. C EXPARG - LARGEST WORKING PRECISION ARGUMENT THAT THE LIBRARY C EXP ROUTINE CAN HANDLE AND UPPER LIMIT ON THE C MAGNITUDE OF X WHEN IZE=1. C C APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE: C C IBM/195 CDC/7600 UNIVAC/1108 VAX 11/780 (UNIX) C (D.P.) (S.P.,RNDG) (D.P.) (S.P.) (D.P.) C C NSIG 16 14 18 8 17 C ENTEN 1.0D75 1.0E322 1.0D307 1.0E38 1.0D38 C ENSIG 1.0D16 1.0E14 1.0D18 1.0E8 1.0D17 C RTNSIG 1.0D-4 1.0E-4 1.0D-5 1.0E-2 1.0D-4 C ENMTEN 2.2D-78 1.0E-290 1.2D-308 1.2E-37 1.2D-37 C XLARGE 1.0D4 1.0E4 1.0D4 1.0E4 1.0D4 C EXPARG 174.0D0 740.0E0 709.0D0 88.0E0 88.0D0 C C******************************************************************* C******************************************************************* C C C ERROR RETURNS C C IN CASE OF AN ERROR, NCALC .NE. NB, AND NOT ALL I'S ARE C CALCULATED TO THE DESIRED ACCURACY. C C NCALC .LT. 0: AN ARGUMENT IS OUT OF RANGE. FOR EXAMPLE, C NB .LE. 0, IZE IS NOT 1 OR 2, OR IZE=1 AND ABS(X) .GE. EXPARG. C IN THIS CASE, THE B-VECTOR IS NOT CALCULATED, AND NCALC IS C SET TO MIN0(NB,0)-1 SO THAT NCALC .NE. NB. C C NB .GT. NCALC .GT. 0: NOT ALL REQUESTED FUNCTION VALUES COULD C BE CALCULATED ACCURATELY. THIS USUALLY OCCURS BECAUSE NB IS C MUCH LARGER THAN ABS(X). IN THIS CASE, B(N) IS CALCULATED C TO THE DESIRED ACCURACY FOR N .LE. NCALC, BUT PRECISION C IS LOST FOR NCALC .LT. N .LE. NB. IF B(N) DOES NOT VANISH C FOR N .GT. NCALC (BECAUSE IT IS TOO SMALL TO BE REPRESENTED), C AND B(N)/B(NCALC) = 10**(-K), THEN ONLY THE FIRST NSIG-K C SIGNIFICANT FIGURES OF B(N) CAN BE TRUSTED. C C C OTHER SUBPROGRAMS REQUIRED (SINGLE PRECISION VERSION) C C EXP,GAMMA,AMAX1,SQRT,FLOAT,IFIX,MIN0 C C OTHER SUBPROGRAMS REQUIRED (DOUBLE PRECISION VERSION) C C DBLE,DEXP,DGAMMA,DMAX1,DSQRT,FLOAT,IFIX,MIN0,SNGL C C C ACKNOWLEDGEMENT C C THIS PROGRAM IS BASED ON A PROGRAM WRITTEN BY DAVID J. C SOOKNE (2) THAT COMPUTES VALUES OF THE BESSEL FUNCTIONS J OR C I OF REAL ARGUMENT AND INTEGER ORDER. MODIFICATIONS INCLUDE C THE RESTRICTION OF THE COMPUTATION TO THE I BESSEL FUNCTION C OF NON-NEGATIVE REAL ARGUMENT, THE EXTENSION OF THE COMPUTATION C TO ARBITRARY POSITIVE ORDER, THE INCLUSION OF OPTIONAL C EXPONENTIAL SCALING, AND THE ELIMINATION OF MOST UNDERFLOW. C C REFERENCES C C (1) OLVER, F. W. J., AND SOOKNE, D. J., "A NOTE ON BACKWARD C RECURRENCE ALGORITHMS," MATH. COMP. 26, 1972, PP 941-947. C C (2) SOOKNE, D. J., "BESSEL FUNCTIONS OF REAL ARGUMENT AND C INTEGER ORDER," NBS JOUR. OF RES. B. 77B, 1973, PP 125-132. C C C MODIFIED BY: W. J. CODY C APPLIED MATHEMATICS DIVISION C ARGONNE NATIONAL LABORATORY C ARGONNE, IL 60439 C C LATEST MODIFICATION: MAY 18, 1982 C C------------------------------------------------------------------- INTEGER IZE, L, MAGX, N, NB, NBMX, NCALC, NEND, NSIG, NSTART CS REAL ALPHA,B,EM,EMPAL,EMP2AL,EN,ENMTEN,ENSIG, CS 2 ENTEN,EXPARG,GAMMA,HALF,HALFX,ONE,P,PLAST,POLD,PSAVE,PSAVEL, CS 3 RTNSIG,SUM,TEMPA,TEMPB,TEMPC,TEST,TOVER,TWO,X,XLARGE,ZERO DOUBLE PRECISION ALPHA, B, EM, EMPAL, EMP2AL, EN, ENMTEN, ENSIG, * ENTEN, EXPARG, DGAMMA, HALF, HALFX, ONE, P, PLAST, POLD, PSAVE, * PSAVEL, RTNSIG, SUM, TEMPA, TEMPB, TEMPC, TEST, TOVER, TWO, X, * XLARGE, ZERO DIMENSION B(NB) C------------------------------------------------------------------- C MATHEMATICAL CONSTANTS C------------------------------------------------------------------- CS DATA ONE,TWO,ZERO,HALF/1.0E0,2.0E0,0.0E0,0.5E0/ DATA ONE, TWO, ZERO, HALF /1.0D0,2.0D0,0.0D0,0.5D0/ C------------------------------------------------------------------- C MACHINE DEPENDENT PARAMETERS C------------------------------------------------------------------- CS DATA NSIG,XLARGE,EXPARG / 7,1.0E4,88.0E0/ CS DATA ENTEN,ENSIG,RTNSIG/1.0E38,1.0E7,1.0E-2/ CS DATA ENMTEN/1.2E-37/ DATA NSIG, XLARGE, EXPARG /17,1.0D4,88.0D0/ DATA ENTEN, ENSIG, RTNSIG /1.0D38,1.0D17,1.0D-4/ DATA ENMTEN /1.2D-37/ C------------------------------------------------------------------- CS MAGX = IFIX(X) MAGX = IFIX(SNGL(X)) IF ((NB.GT.0) .AND. (X.GE.ZERO) .AND. (ALPHA.GE.ZERO) .AND. * (ALPHA.LT.ONE) .AND. (((IZE.EQ.1) .AND. (X.LE.EXPARG)) .OR. * ((IZE.EQ.2) .AND. (X.LE.XLARGE)))) GO TO 10 C------------------------------------------------------------------- C ERROR RETURN -- X,NB,OR IZE IS OUT OF RANGE C------------------------------------------------------------------- NCALC = MIN0(NB,0) - 1 RETURN C------------------------------------------------------------------- C USE 2-TERM ASCENDING SERIES FOR SMALL X C------------------------------------------------------------------- 10 NCALC = NB IF (X.LT.RTNSIG) GO TO 210 C------------------------------------------------------------------- C INITIALIZE THE FORWARD SWEEP, THE P-SEQUENCE OF OLVER C------------------------------------------------------------------- NBMX = NB - MAGX N = MAGX + 1 CS EN = FLOAT(N+N) + (ALPHA+ALPHA) EN = DBLE(FLOAT(N+N)) + (ALPHA+ALPHA) PLAST = ONE P = EN/X C------------------------------------------------------------------- C CALCULATE GENERAL SIGNIFICANCE TEST C------------------------------------------------------------------- TEST = ENSIG + ENSIG CS IF (2*MAGX .GT. 5*NSIG) TEST = SQRT(TEST*P) IF (2*MAGX.GT.5*NSIG) TEST = DSQRT(TEST*P) CS IF (2*MAGX .LE. 5*NSIG) TEST = TEST / 1.585E0**MAGX IF (2*MAGX.LE.5*NSIG) TEST = TEST/1.585D0**MAGX IF (NBMX.LT.3) GO TO 30 C------------------------------------------------------------------- C CALCULATE P-SEQUENCE UNTIL N = NB-1. CHECK FOR POSSIBLE OVERFLOW. C------------------------------------------------------------------- TOVER = ENTEN/ENSIG NSTART = MAGX + 2 NEND = NB - 1 DO 20 N=NSTART,NEND EN = EN + TWO POLD = PLAST PLAST = P P = EN*PLAST/X + POLD IF (P.GT.TOVER) GO TO 40 20 CONTINUE N = NEND CS EN = FLOAT(N+N) + (ALPHA+ALPHA) EN = DBLE(FLOAT(N+N)) + (ALPHA+ALPHA) C------------------------------------------------------------------- C CALCULATE SPECIAL SIGNIFICANCE TEST FOR NBMX.GT.2. C------------------------------------------------------------------- CS TEST = AMAX1(TEST,SQRT(PLAST*ENSIG)*SQRT(P+P)) TEST = DMAX1(TEST,DSQRT(PLAST*ENSIG)*DSQRT(P+P)) C------------------------------------------------------------------- C CALCULATE P-SEQUENCE UNTIL SIGNIFICANCE TEST PASSES C------------------------------------------------------------------- 30 N = N + 1 EN = EN + TWO POLD = PLAST PLAST = P P = EN*PLAST/X + POLD IF (P.LT.TEST) GO TO 30 GO TO 80 C------------------------------------------------------------------- C TO AVOID OVERFLOW, DIVIDE P-SEQUENCE BY TOVER. CALCULATE C P-SEQUENCE UNTIL ABS(P).GT.1. C------------------------------------------------------------------- 40 TOVER = ENTEN P = P/TOVER PLAST = PLAST/TOVER PSAVE = P PSAVEL = PLAST NSTART = N + 1 50 N = N + 1 EN = EN + TWO POLD = PLAST PLAST = P P = EN*PLAST/X + POLD IF (P.LE.ONE) GO TO 50 TEMPB = EN/X C------------------------------------------------------------------- C CALCULATE BACKWARD TEST, AND FIND NCALC, THE HIGHEST N C SUCH THAT THE TEST IS PASSED. C------------------------------------------------------------------- TEST = POLD*PLAST*(HALF-HALF/(TEMPB*TEMPB))/ENSIG P = PLAST*TOVER N = N - 1 EN = EN - TWO NEND = MIN0(NB,N) DO 60 L=NSTART,NEND NCALC = L POLD = PSAVEL PSAVEL = PSAVE PSAVE = EN*PSAVEL/X + POLD IF (PSAVE*PSAVEL.GT.TEST) GO TO 70 60 CONTINUE NCALC = NEND + 1 70 NCALC = NCALC - 1 C------------------------------------------------------------------- C INITIALIZE THE BACKWARD RECURSION AND THE NORMALIZATION SUM C------------------------------------------------------------------- 80 N = N + 1 EN = EN + TWO TEMPB = ZERO TEMPA = ONE/P CS EM = FLOAT(N) - ONE EM = DBLE(FLOAT(N)) - ONE EMPAL = EM + ALPHA EMP2AL = (EM-ONE) + (ALPHA+ALPHA) SUM = TEMPA*EMPAL*EMP2AL/EM NEND = N - NB IF (NEND) 130, 110, 90 C------------------------------------------------------------------- C RECUR BACKWARD VIA DIFFERENCE EQUATION, CALCULATING (BUT C NOT STORING) B(N), UNTIL N = NB. C------------------------------------------------------------------- 90 DO 100 L=1,NEND N = N - 1 EN = EN - TWO TEMPC = TEMPB TEMPB = TEMPA TEMPA = (EN*TEMPB)/X + TEMPC EM = EM - ONE EMP2AL = EMP2AL - ONE IF (N.EQ.1) GO TO 110 IF (N.EQ.2) EMP2AL = ONE EMPAL = EMPAL - ONE SUM = (SUM+TEMPA*EMPAL)*EMP2AL/EM 100 CONTINUE C------------------------------------------------------------------- C STORE B(NB) C------------------------------------------------------------------- 110 B(N) = TEMPA IF (NB.GT.1) GO TO 120 SUM = (SUM+SUM) + TEMPA GO TO 190 C------------------------------------------------------------------- C CALCULATE AND STORE B(NB-1) C------------------------------------------------------------------- 120 N = N - 1 EN = EN - TWO B(N) = (EN*TEMPA)/X + TEMPB IF (N.EQ.1) GO TO 180 EM = EM - ONE EMP2AL = EMP2AL - ONE IF (N.EQ.2) EMP2AL = ONE EMPAL = EMPAL - ONE SUM = (SUM+B(N)*EMPAL)*EMP2AL/EM GO TO 150 C------------------------------------------------------------------- C N.LT.NB, SO STORE B(N) AND SET HIGHER ORDERS TO ZERO C------------------------------------------------------------------- 130 B(N) = TEMPA NEND = -NEND DO 140 L=1,NEND B(N+L) = ZERO 140 CONTINUE 150 NEND = N - 2 IF (NEND.EQ.0) GO TO 170 C------------------------------------------------------------------- C CALCULATE VIA DIFFERENCE EQUATION AND STORE B(N), UNTIL N = 2 C------------------------------------------------------------------- DO 160 L=1,NEND N = N - 1 EN = EN - TWO B(N) = (EN*B(N+1))/X + B(N+2) EM = EM - ONE EMP2AL = EMP2AL - ONE IF (N.EQ.2) EMP2AL = ONE EMPAL = EMPAL - ONE SUM = (SUM+B(N)*EMPAL)*EMP2AL/EM 160 CONTINUE C------------------------------------------------------------------- C CALCULATE B(1) C------------------------------------------------------------------- 170 B(1) = TWO*EMPAL*B(2)/X + B(3) 180 SUM = (SUM+SUM) + B(1) C------------------------------------------------------------------- C NORMALIZE. DIVIDE ALL B(N) BY SUM. C------------------------------------------------------------------- 190 CONTINUE CS IF (ALPHA.NE.ZERO)SUM=SUM*GAMMA(ONE+ALPHA)*(X*HALF)**(-ALPHA) IF (ALPHA.NE.ZERO) SUM = SUM*DGAMMA(ONE+ALPHA)*(X*HALF)**(-ALPHA) CS IF (IZE .EQ. 1) SUM = SUM * EXP(-X) IF (IZE.EQ.1) SUM = SUM*DEXP(-X) TEMPA = ENMTEN IF (SUM.GT.ONE) TEMPA = TEMPA*SUM DO 200 N=1,NB IF (B(N).LT.TEMPA) B(N) = ZERO B(N) = B(N)/SUM 200 CONTINUE RETURN C------------------------------------------------------------------- C TWO-TERM ASCENDING SERIES FOR SMALL X C------------------------------------------------------------------- 210 TEMPA = ONE EMPAL = ONE + ALPHA HALFX = ZERO IF (X.GT.ENMTEN) HALFX = HALF*X CS IF (ALPHA .NE. ZERO) TEMPA = HALFX ** ALPHA / GAMMA(EMPAL) IF (ALPHA.NE.ZERO) TEMPA = HALFX**ALPHA/DGAMMA(EMPAL) CS IF (IZE .EQ. 2) TEMPA = TEMPA * EXP(-X) IF (IZE.EQ.2) TEMPA = TEMPA*DEXP(-X) TEMPB = ZERO IF ((X+ONE).GT.ONE) TEMPB = HALFX*HALFX B(1) = TEMPA + TEMPA*TEMPB/EMPAL IF ((X.NE.ZERO) .AND. (B(1).EQ.ZERO)) NCALC = 0 IF (NB.EQ.1) GO TO 250 IF (X.GT.ZERO) GO TO 230 DO 220 N=2,NB B(N) = ZERO 220 CONTINUE GO TO 250 C------------------------------------------------------------------- C CALCULATE HIGHER ORDER FUNCTIONS C------------------------------------------------------------------- 230 TEMPC = HALFX TOVER = (ENMTEN+ENMTEN)/X IF (TEMPB.NE.ZERO) TOVER = ENMTEN/TEMPB DO 240 N=2,NB TEMPA = TEMPA/EMPAL EMPAL = EMPAL + ONE TEMPA = TEMPA*TEMPC IF (TEMPA.LE.TOVER*EMPAL) TEMPA = ZERO B(N) = TEMPA + TEMPA*TEMPB/EMPAL IF ((B(N).EQ.ZERO) .AND. (NCALC.GT.N)) NCALC = N - 1 240 CONTINUE 250 RETURN C ---------- LAST CARD OF RIBESL ---------- END DOUBLE PRECISION FUNCTION DGAMMA(X) c*********************************************************************72 C cc DGAMMA computes the gamma function. c C THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A REAL ARGUMENT X. C COMPUTATION IS BASED ON AN ALGORITHM OUTLINED IN W. J. CODY, C 'AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL FUNCTIONS', C LECTURE NOTES IN MATHEMATICS, 506, NUMERICAL ANALYSIS DUNDEE, C 1975, G. A. WATSON (ED.), SPRINGER VERLAG, BERLIN, 1976. THE C PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA C FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS. COEFFICIENTS C FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED. C THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM HART, ET. AL., C COMPUTER APPROXIMATIONS, WILEY AND SONS, NEW YORK, 1968. LOWER C ORDER APPROXIMATIONS CAN BE SUBSTITUTED FOR THESE ON MACHINES C WITH LESS PRECISE ARITHMETIC. C C C******************************************************************* C******************************************************************* C C EXPLANATION OF MACHINE-DEPENDENT CONSTANTS C C EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT C 1.0 + EPS .GT. 1.0 C XBIG - THE LARGEST ARGUMENT FOR WHICH GAMMA(X) IS REPRESENTABLE C IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION C GAMMA(XBIG) = XINF. C XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT C 1/XMININ IS MACHINE REPRESENTABLE. C XINF - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER. C C APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE: C C IBM/195 CDC/7600 UNIVAC/1108 VAX 11/780 (UNIX) C (D.P.) (S.P.,RNDG) (D.P.) (S.P.) (D.P.) C C EPS 2.221D-16 3.553E-15 1.735D-18 5.961E-08 1.388D-16 C XBIG 57.574 177.802 171.489 34.844 34.844 C XMININ 1.382D-76 3.132E-294 1.113D-308 5.883E-39 5.883D-39 C XINF 7.23D+75 1.26E+322 8.98D+307 1.70E+38 1.70D+38 C C******************************************************************* C******************************************************************* C C C ERROR RETURNS C C THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR C WHEN OVERFLOW WOULD OCCUR. THE COMPUTATION IS BELIEVED C TO BE FREE OF UNDERFLOW AND OVERFLOW. C C C C OTHER SUBPROGRAMS REQUIRED (SINGLE PRECISION VERSION) C C ALOG,EXP,FLOAT,IFIX,SIN C C OTHER SUBPROGRAMS REQUIRED (DOUBLE PRECISION VERSION) C C DBLE,DEXP,DLOG,DSIN,FLOAT,IFIX,SNGL C C C C AUTHOR: W. J. CODY C APPLIED MATHEMATICS DIVISION C ARGONNE NATIONAL LABORATORY C ARGONNE, IL 60439 C C LATEST MODIFICATION: MAY 18, 1982 C C---------------------------------------------------------------------- CS REAL C,EPS,FACT,HALF,ONE,P,PI,Q,RES,SQRTPI, CS 1 SUM,TWELVE,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO DOUBLE PRECISION C, EPS, FACT, HALF, ONE, P, PI, Q, RES, SQRTPI, * SUM, TWELVE, X, XBIG, XDEN, XINF, XMININ, XNUM, Y, Y1, YSQ, Z, * ZERO INTEGER I, J, N LOGICAL PARITY DIMENSION C(7), P(8), Q(8) C---------------------------------------------------------------------- C MATHEMATICAL CONSTANTS C---------------------------------------------------------------------- CS DATA ONE,HALF,TWELVE,ZERO/1.0E0,0.5E0,12.0E0,0.0E0/ CS DATA SQRTPI/0.9189385332046727417803297E0/ CS DATA PI/3.1415926535897932384626434E0/ DATA ONE, HALF, TWELVE, ZERO /1.0D0,0.5D0,12.0D0,0.0D0/ DATA SQRTPI /0.9189385332046727417803297D0/ DATA PI /3.1415926535897932384626434D0/ C---------------------------------------------------------------------- C MACHINE DEPENDENT PARAMETERS C---------------------------------------------------------------------- CS DATA XBIG,XMININ,EPS/34.844E0,5.883E-39,5.961E-08/ CS DATA XINF/1.7014E38/ DATA XBIG, XMININ, EPS /34.844D0,5.883D-39,1.388D-16/ DATA XINF /1.7014D38/ C---------------------------------------------------------------------- C NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX C APPROXIMATION OVER (1,2). C---------------------------------------------------------------------- CS DATA P/-1.71618513886549492533811E+0,2.47656508055759199108314E+1, CS 1 -3.79804256470945635097577E+2,6.29331155312818442661052E+2, CS 2 8.66966202790413211295064E+2,-3.14512729688483675254357E+4, CS 3 -3.61444134186911729807069E+4,6.64561438202405440627855E+4/ CS DATA Q/-3.08402300119738975254353E+1,3.15350626979604161529144E+2, CS 1 -1.01515636749021914166146E+3,-3.10777167157231109440444E+3, CS 2 2.25381184209801510330112E+4,4.75584627752788110767815E+3, CS 3 -1.34659959864969306392456E+5,-1.15132259675553483497211E+5/ DATA P /-1.71618513886549492533811D+0,2.47656508055759199108314D+1 * ,-3.79804256470945635097577D+2,6.29331155312818442661052D+2, * 8.66966202790413211295064D+2,-3.14512729688483675254357D+4, * -3.61444134186911729807069D+4,6.64561438202405440627855D+4/ DATA Q /-3.08402300119738975254353D+1,3.15350626979604161529144D+2 * ,-1.01515636749021914166146D+3,-3.10777167157231109440444D+3, * 2.25381184209801510330112D+4,4.75584627752788110767815D+3, * -1.34659959864969306392456D+5,-1.15132259675553483497211D+5/ C---------------------------------------------------------------------- C COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF). C---------------------------------------------------------------------- CS DATA C/-1.910444077728E-03,8.4171387781295E-04, CS 1 -5.952379913043012E-04,7.93650793500350248E-04, CS 2 -2.777777777777681622553E-03,8.333333333333333331554247E-02, CS 3 5.7083835261E-03/ DATA C /-1.910444077728D-03,8.4171387781295D-04, * -5.952379913043012D-04,7.93650793500350248D-04, * -2.777777777777681622553D-03,8.333333333333333331554247D-02, * 5.7083835261D-03/ C---------------------------------------------------------------------- PARITY = .FALSE. FACT = ONE N = 0 Y = X IF (Y.GT.ZERO) GO TO 10 C---------------------------------------------------------------------- C ARGUMENT IS NEGATIVE C---------------------------------------------------------------------- Y = -X CS J = IFIX(Y) J = IFIX(SNGL(Y)) CS RES = Y - FLOAT(J) RES = Y - DBLE(FLOAT(J)) IF (RES.EQ.ZERO) GO TO 100 IF (J.NE.(J/2)*2) PARITY = .TRUE. CS FACT = -PI / SIN(PI*RES) FACT = -PI/DSIN(PI*RES) Y = Y + ONE C---------------------------------------------------------------------- C ARGUMENT IS POSITIVE C---------------------------------------------------------------------- 10 IF (Y.LT.EPS) GO TO 90 IF (Y.GE.TWELVE) GO TO 70 Y1 = Y IF (Y.GE.ONE) GO TO 20 C---------------------------------------------------------------------- C 0.0 .LT. ARGUMENT .LT. 1.0 C---------------------------------------------------------------------- Z = Y Y = Y + ONE GO TO 30 C---------------------------------------------------------------------- C 1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY C---------------------------------------------------------------------- CS210 N = IFIX(Y) - 1 20 N = IFIX(SNGL(Y)) - 1 CS Y = Y - FLOAT(N) Y = Y - DBLE(FLOAT(N)) Z = Y - ONE C---------------------------------------------------------------------- C EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0 C---------------------------------------------------------------------- 30 XNUM = ZERO XDEN = ONE DO 40 I=1,8 XNUM = (XNUM+P(I))*Z XDEN = XDEN*Z + Q(I) 40 CONTINUE RES = XNUM/XDEN + ONE IF (Y.EQ.Y1) GO TO 110 IF (Y1.GT.Y) GO TO 50 C---------------------------------------------------------------------- C ADJUST RESULT FOR CASE 0.0 .LT. ARGUMENT .LT. 1.0 C---------------------------------------------------------------------- RES = RES/Y1 GO TO 110 C---------------------------------------------------------------------- C ADJUST RESULT FOR CASE 2.0 .LT. ARGUMENT .LT. 12.0 C---------------------------------------------------------------------- 50 DO 60 I=1,N RES = RES*Y Y = Y + ONE 60 CONTINUE GO TO 110 C---------------------------------------------------------------------- C EVALUATE FOR ARGUMENT .GE. 12.0, C---------------------------------------------------------------------- 70 IF (Y.GT.XBIG) GO TO 100 YSQ = Y*Y SUM = C(7) DO 80 I=1,6 SUM = SUM/YSQ + C(I) 80 CONTINUE SUM = SUM/Y - Y + SQRTPI CS SUM = SUM + (Y-HALF)*ALOG(Y) SUM = SUM + (Y-HALF)*DLOG(Y) CS RES = EXP(SUM) RES = DEXP(SUM) GO TO 110 C---------------------------------------------------------------------- C ARGUMENT .LT. EPS C---------------------------------------------------------------------- 90 IF (Y.LT.XMININ) GO TO 100 RES = ONE/Y GO TO 110 C---------------------------------------------------------------------- C RETURN FOR SINGULARITIES, EXTREME ARGUMENTS, ETC. C---------------------------------------------------------------------- CS700 GAMMA = XINF 100 DGAMMA = XINF GO TO 120 C---------------------------------------------------------------------- C FINAL ADJUSTMENTS AND RETURN C---------------------------------------------------------------------- 110 IF (PARITY) RES = -RES IF (FACT.NE.ONE) RES = FACT/RES CS GAMMA = RES DGAMMA = RES 120 RETURN C ---------- LAST CARD OF GAMMA ---------- END subroutine timestamp ( ) c*********************************************************************72 c cc TIMESTAMP prints out the current YMDHMS date as a timestamp. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 12 June 2014 c c Author: c c John Burkardt c c Parameters: c c None c implicit none character * ( 8 ) ampm integer d character * ( 8 ) date integer h integer m integer mm character * ( 9 ) month(12) integer n integer s character * ( 10 ) time integer y save month data month / & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' / call date_and_time ( date, time ) read ( date, '(i4,i2,i2)' ) y, m, d read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm if ( h .lt. 12 ) then ampm = 'AM' else if ( h .eq. 12 ) then if ( n .eq. 0 .and. s .eq. 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h .lt. 12 ) then ampm = 'PM' else if ( h .eq. 12 ) then if ( n .eq. 0 .and. s .eq. 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, & '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, & trim ( ampm ) return end