subroutine aaaaaa ( ver ) !*****************************************************************************80 ! !! AAAAAA is the SLATEC Common Mathematical Library disclaimer and version. ! !***LIBRARY SLATEC !***CATEGORY Z !***TYPE ALL (AAAAAA-A) !***KEYWORDS DISCLAIMER, DOCUMENTATION, VERSION !***AUTHOR SLATEC Common Mathematical Library Committee !***DESCRIPTION ! ! The SLATEC Common Mathematical Library is issued by the following ! ! Air Force Weapons Laboratory, Albuquerque ! Lawrence Livermore National Laboratory, Livermore ! Los Alamos National Laboratory, Los Alamos ! National Institute of Standards and Technology, Washington ! National Energy Research Supercomputer Center, Livermore ! Oak Ridge National Laboratory, Oak Ridge ! Sandia National Laboratories, Albuquerque ! Sandia National Laboratories, Livermore ! ! All questions concerning the distribution of the library should be ! directed to the NATIONAL ENERGY SOFTWARE CENTER, 9700 Cass Ave., ! Argonne, Illinois 60439, and not to the authors of the subprograms. ! ! * * * * * Notice * * * * * ! ! This material was prepared as an account of work sponsored by the ! United States Government. Neither the United States, nor the ! Department of Energy, nor the Department of Defense, nor any of ! their employees, nor any of their contractors, subcontractors, or ! their employees, makes any warranty, expressed or implied, or ! assumes any legal liability or responsibility for the accuracy, ! completeness, or usefulness of any information, apparatus, product, ! or process disclosed, or represents that its use would not infringe ! upon privately owned rights. ! ! *Usage: ! ! CHARACTER * 16 VER ! ! call AAAAAA (VER) ! ! *Arguments: ! ! VER:OUT will contain the version number of the SLATEC CML. ! ! *Description: ! ! This routine contains the SLATEC Common Mathematical Library ! disclaimer and can be used to return the library version number. ! !***REFERENCES Kirby W. Fong, Thomas H. Jefferson, Tokihiko Suyehiro ! and Lee Walton, Guide to the SLATEC Common Mathema- ! tical Library, April 10, 1990. !***ROUTINES CALLED (NONE) !***REVISION HISTORY (YYMMDD) ! 800424 DATE WRITTEN ! 890414 REVISION DATE from Version 3.2 ! 890713 Routine modified to return version number. (WRB) ! 900330 Prologue converted to Version 4.0 format. (BAB) ! 920501 Reformatted the REFERENCES section. (WRB) ! 921215 Updated for Version 4.0. (WRB) ! 930701 Updated for Version 4.1. (WRB) !***END PROLOGUE AAAAAA CHARACTER * (*) VER !***FIRST EXECUTABLE STATEMENT AAAAAA VER = ' 4.1' return end function acosh ( x ) !*****************************************************************************80 ! !! ACOSH computes the arc hyperbolic cosine. ! !***LIBRARY SLATEC (FNLIB) !***CATEGORY C4C !***TYPE SINGLE PRECISION (ACOSH-S, DACOSH-D, CACOSH-C) !***KEYWORDS ACOSH, ARC HYPERBOLIC COSINE, ELEMENTARY FUNCTIONS, FNLIB, ! INVERSE HYPERBOLIC COSINE !***AUTHOR Fullerton, W., (LANL) !***DESCRIPTION ! ! ACOSH(X) computes the arc hyperbolic cosine of X. ! !***REFERENCES (NONE) !***ROUTINES CALLED R1MACH, XERMSG !***REVISION HISTORY (YYMMDD) ! 770401 DATE WRITTEN ! 890531 Changed all specific intrinsics to generic. (WRB) ! 890531 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) ! 900326 Removed duplicate information from DESCRIPTION section. ! (WRB) !***END PROLOGUE ACOSH real, parameter :: aln2 = 0.69314718055994530942E+00 real x real, save :: xmax = 0.0E+00 !***FIRST EXECUTABLE STATEMENT ACOSH if ( xmax == 0.0E+00 ) then xmax = 1.0E+00 / sqrt ( r1mach(3) ) end if if ( x < 1.0E+00 ) then call XERMSG ( 'SLATEC', 'ACOSH', 'X LESS THAN 1', 1, 2 ) end if if ( x < xmax ) then acosh = log ( x + sqrt ( x * x - 1.0E+00 ) ) end if if ( x >= xmax ) then acosh = aln2 + log ( x ) end if return end function ai ( x ) !*****************************************************************************80 ! !! AI evaluates the Airy function. ! !***LIBRARY SLATEC (FNLIB) !***CATEGORY C10D !***TYPE SINGLE PRECISION (AI-S, DAI-D) !***KEYWORDS AIRY FUNCTION, FNLIB, SPECIAL FUNCTIONS !***AUTHOR Fullerton, W., (LANL) !***DESCRIPTION ! ! AI(X) computes the Airy function Ai(X) ! Series for AIF on the interval -1.00000D+00 to 1.00000D+00 ! with weighted error 1.09E-19 ! log weighted error 18.96 ! significant figures required 17.76 ! decimal places required 19.44 ! ! Series for AIG on the interval -1.00000D+00 to 1.00000D+00 ! with weighted error 1.51E-17 ! log weighted error 16.82 ! significant figures required 15.19 ! decimal places required 17.27 ! !***REFERENCES (NONE) !***ROUTINES CALLED AIE, CSEVL, INITS, R1MACH, R9AIMP, XERMSG !***REVISION HISTORY (YYMMDD) ! 770701 DATE WRITTEN ! 890531 Changed all specific intrinsics to generic. (WRB) ! 890531 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) ! 900326 Removed duplicate information from DESCRIPTION section. ! (WRB) ! 920618 Removed space from variable names. (RWC, WRB) !***END PROLOGUE AI ! implicit none real ai real aie real, save, dimension ( 9 ) :: aifcs = (/ & -0.03797135849666999750E0, & 0.05919188853726363857E0, & 0.00098629280577279975E0, & 0.00000684884381907656E0, & 0.00000002594202596219E0, & 0.00000000006176612774E0, & 0.00000000000010092454E0, & 0.00000000000000012014E0, & 0.00000000000000000010E0 /) real, save, dimension ( 8 ) :: aigcs = (/ & 0.01815236558116127E0, & 0.02157256316601076E0, & 0.00025678356987483E0, & 0.00000142652141197E0, & 0.00000000457211492E0, & 0.00000000000952517E0, & 0.00000000000001392E0, & 0.00000000000000001E0 /) real csevl logical, save :: first = .true. integer inits integer, save :: naif integer, save :: naig real r1mach real theta real x real, save :: x3sml real xm real, save :: xmax real xmaxt real z !***FIRST EXECUTABLE STATEMENT AI if ( first ) then naif = inits ( aifcs, 9, 0.1 * r1mach(3) ) naig = inits ( aigcs, 8, 0.1 * r1mach(3) ) x3sml = r1mach(3)**0.3334 xmaxt = ( -1.5 * log ( r1mach(1) ) )**0.6667 xmax = xmaxt - xmaxt * log ( xmaxt ) & / ( 4.0 * sqrt ( xmaxt ) + 1.0 ) - 0.01 first = .false. end if if ( x < -1.0 ) then call r9aimp ( x, xm, theta ) ai = xm * cos ( theta ) else if ( x <= 1.0 ) then if ( abs ( x ) <= x3sml ) then z = 0.0 else z = x**3 end if ai = 0.375 + ( csevl ( z, aifcs, naif ) & - x * ( 0.25 + csevl ( z, aigcs, naig ) ) ) else if ( x <= xmax ) then ai = aie ( x ) * exp ( -2.0 * x * sqrt ( x ) / 3.0 ) else ai = 0.0 call xermsg ( 'SLATEC', 'AI', 'X so big AI underflows', 1, 1 ) end if return end function aie ( x ) !*****************************************************************************80 ! !! AIE calculates the Airy function for a negative argument... ! and an exponentially scaled Airy function for a non-negative argument. ! !***LIBRARY SLATEC (FNLIB) !***CATEGORY C10D !***TYPE SINGLE PRECISION (AIE-S, DAIE-D) !***KEYWORDS EXPONENTIALLY SCALED AIRY FUNCTION, FNLIB, ! SPECIAL FUNCTIONS !***AUTHOR Fullerton, W., (LANL) !***DESCRIPTION ! ! AIE(X) computes the exponentially scaled Airy function for ! non-negative X. It evaluates AI(X) for X <= 0.0 and ! EXP(ZETA)*AI(X) for X >= 0.0 where ZETA = (2.0/3.0)*(X**1.5). ! ! Series for AIF on the interval -1.00000D+00 to 1.00000D+00 ! with weighted error 1.09E-19 ! log weighted error 18.96 ! significant figures required 17.76 ! decimal places required 19.44 ! ! Series for AIG on the interval -1.00000D+00 to 1.00000D+00 ! with weighted error 1.51E-17 ! log weighted error 16.82 ! significant figures required 15.19 ! decimal places required 17.27 ! ! Series for AIP on the interval 0. to 1.00000D+00 ! with weighted error 5.10E-17 ! log weighted error 16.29 ! significant figures required 14.41 ! decimal places required 17.06 ! !***REFERENCES (NONE) !***ROUTINES CALLED CSEVL, INITS, R1MACH, R9AIMP !***REVISION HISTORY (YYMMDD) ! 770701 DATE WRITTEN ! 890206 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 920618 Removed space from variable names. (RWC, WRB) !***END PROLOGUE AIE DIMENSION AIFCS(9), AIGCS(8), AIPCS(34) LOGICAL FIRST SAVE AIFCS, AIGCS, AIPCS, NAIF, NAIG, & NAIP, X3SML, X32SML, XBIG, FIRST DATA AIFCS( 1) / -0.03797135849666999750E0 / DATA AIFCS( 2) / 0.05919188853726363857E0 / DATA AIFCS( 3) / 0.00098629280577279975E0 / DATA AIFCS( 4) / 0.00000684884381907656E0 / DATA AIFCS( 5) / 0.00000002594202596219E0 / DATA AIFCS( 6) / 0.00000000006176612774E0 / DATA AIFCS( 7) / 0.00000000000010092454E0 / DATA AIFCS( 8) / 0.00000000000000012014E0 / DATA AIFCS( 9) / 0.00000000000000000010E0 / DATA AIGCS( 1) / 0.01815236558116127E0 / DATA AIGCS( 2) / 0.02157256316601076E0 / DATA AIGCS( 3) / 0.00025678356987483E0 / DATA AIGCS( 4) / 0.00000142652141197E0 / DATA AIGCS( 5) / 0.00000000457211492E0 / DATA AIGCS( 6) / 0.00000000000952517E0 / DATA AIGCS( 7) / 0.00000000000001392E0 / DATA AIGCS( 8) / 0.00000000000000001E0 / DATA AIPCS( 1) / -0.0187519297793868E0 / DATA AIPCS( 2) / -0.0091443848250055E0 / DATA AIPCS( 3) / 0.0009010457337825E0 / DATA AIPCS( 4) / -0.0001394184127221E0 / DATA AIPCS( 5) / 0.0000273815815785E0 / DATA AIPCS( 6) / -0.0000062750421119E0 / DATA AIPCS( 7) / 0.0000016064844184E0 / DATA AIPCS( 8) / -0.0000004476392158E0 / DATA AIPCS( 9) / 0.0000001334635874E0 / DATA AIPCS(10) / -0.0000000420735334E0 / DATA AIPCS(11) / 0.0000000139021990E0 / DATA AIPCS(12) / -0.0000000047831848E0 / DATA AIPCS(13) / 0.0000000017047897E0 / DATA AIPCS(14) / -0.0000000006268389E0 / DATA AIPCS(15) / 0.0000000002369824E0 / DATA AIPCS(16) / -0.0000000000918641E0 / DATA AIPCS(17) / 0.0000000000364278E0 / DATA AIPCS(18) / -0.0000000000147475E0 / DATA AIPCS(19) / 0.0000000000060851E0 / DATA AIPCS(20) / -0.0000000000025552E0 / DATA AIPCS(21) / 0.0000000000010906E0 / DATA AIPCS(22) / -0.0000000000004725E0 / DATA AIPCS(23) / 0.0000000000002076E0 / DATA AIPCS(24) / -0.0000000000000924E0 / DATA AIPCS(25) / 0.0000000000000417E0 / DATA AIPCS(26) / -0.0000000000000190E0 / DATA AIPCS(27) / 0.0000000000000087E0 / DATA AIPCS(28) / -0.0000000000000040E0 / DATA AIPCS(29) / 0.0000000000000019E0 / DATA AIPCS(30) / -0.0000000000000009E0 / DATA AIPCS(31) / 0.0000000000000004E0 / DATA AIPCS(32) / -0.0000000000000002E0 / DATA AIPCS(33) / 0.0000000000000001E0 / DATA AIPCS(34) / -0.0000000000000000E0 / DATA FIRST /.TRUE./ !***FIRST EXECUTABLE STATEMENT AIE if ( first ) then eta = 0.1 * r1mach(3) naif = inits ( aifcs, 9, eta ) naig = inits ( aigcs, 8, eta ) naip = inits ( aipcs, 34, eta ) x3sml = eta**0.3333 x32sml = 1.3104 * x3sml**2 xbig = r1mach(2)**0.6666 first = .false. end if if ( x < -1.0 ) then call r9aimp ( x, xm, theta ) aie = xm * cos ( theta ) else if ( x <= 1.0 ) then if ( abs ( x ) <= x3sml ) then z = 0.0 else z = x**3 end if aie = 0.375 + ( csevl ( z, aifcs, naif ) & - x * ( 0.25 + csevl ( z, aigcs, naig ) ) ) if ( x32sml < x ) then aie = aie * exp ( 2.0 * x * sqrt ( x ) / 3.0 ) end if else sqrtx = sqrt ( x ) if ( x < xbig ) then z = 2.0 / ( x * sqrtx ) - 1.0 else z = -1.0 end if aie = ( 0.28125 + csevl ( z, aipcs, naip ) ) / sqrt ( sqrtx ) end if return end function ALBETA (A, B) ! !! ALBETA computes the natural logarithm of the complete Beta function. ! !***LIBRARY SLATEC (FNLIB) !***CATEGORY C7B !***TYPE SINGLE PRECISION (ALBETA-S, DLBETA-D, CLBETA-C) !***KEYWORDS FNLIB, LOGARITHM OF THE COMPLETE BETA FUNCTION, ! SPECIAL FUNCTIONS !***AUTHOR Fullerton, W., (LANL) !***DESCRIPTION ! ! ALBETA computes the natural log of the complete beta function. ! ! Input Parameters: ! A real and positive ! B real and positive ! !***REFERENCES (NONE) !***ROUTINES CALLED ALNGAM, ALNREL, GAMMA, R9LGMC, XERMSG !***REVISION HISTORY (YYMMDD) ! 770701 DATE WRITTEN ! 890531 Changed all specific intrinsics to generic. (WRB) ! 890531 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) ! 900326 Removed duplicate information from DESCRIPTION section. ! (WRB) ! 900727 Added EXTERNAL statement. (WRB) !***END PROLOGUE ALBETA EXTERNAL GAMMA SAVE SQ2PIL DATA SQ2PIL / 0.91893853320467274E0 / !***FIRST EXECUTABLE STATEMENT ALBETA P = MIN (A, B) Q = MAX (A, B) ! if (P <= 0.0) call XERMSG ('SLATEC', 'ALBETA', & 'BOTH ARGUMENTS MUST BE GT ZERO', 1, 2) if (P >= 10.0) go to 30 if (Q >= 10.0) go to 20 ! ! P AND Q ARE SMALL. ! ALBETA = LOG(GAMMA(P) * (GAMMA(Q)/GAMMA(P+Q)) ) return ! ! P IS SMALL, BUT Q IS BIG. ! 20 CORR = R9LGMC(Q) - R9LGMC(P+Q) ALBETA = ALNGAM(P) + CORR + P - P*LOG(P+Q) + & (Q-0.5)*ALNREL(-P/(P+Q)) return ! ! P AND Q ARE BIG. ! 30 CORR = R9LGMC(P) + R9LGMC(Q) - R9LGMC(P+Q) ALBETA = -0.5*LOG(Q) + SQ2PIL + CORR + (P-0.5)*LOG(P/(P+Q)) & + Q*ALNREL(-P/(P+Q)) return ! end subroutine ALGAMS (X, ALGAM, SGNGAM) ! !! ALGAMS computes the logarithm of the absolute value of the Gamma function. ! !***LIBRARY SLATEC (FNLIB) !***CATEGORY C7A !***TYPE SINGLE PRECISION (ALGAMS-S, DLGAMS-D) !***KEYWORDS ABSOLUTE VALUE OF THE LOGARITHM OF THE GAMMA FUNCTION, ! FNLIB, SPECIAL FUNCTIONS !***AUTHOR Fullerton, W., (LANL) !***DESCRIPTION ! ! Evaluates the logarithm of the absolute value of the gamma ! function. ! X - input argument ! ALGAM - result ! SGNGAM - is set to the sign of GAMMA(X) and will ! be returned at +1.0 or -1.0. ! !***REFERENCES (NONE) !***ROUTINES CALLED ALNGAM !***REVISION HISTORY (YYMMDD) ! 770701 DATE WRITTEN ! 890531 Changed all specific intrinsics to generic. (WRB) ! 890531 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) !***END PROLOGUE ALGAMS !***FIRST EXECUTABLE STATEMENT ALGAMS ALGAM = ALNGAM(X) SGNGAM = 1.0 if (X > 0.0) RETURN ! INT = MOD (-AINT(X), 2.0) + 0.1 if (INT == 0) SGNGAM = -1.0 ! return end function ALI (X) ! !! ALI computes the logarithmic integral. ! !***LIBRARY SLATEC (FNLIB) !***CATEGORY C5 !***TYPE SINGLE PRECISION (ALI-S, DLI-D) !***KEYWORDS FNLIB, LOGARITHMIC INTEGRAL, SPECIAL FUNCTIONS !***AUTHOR Fullerton, W., (LANL) !***DESCRIPTION ! ! ALI(X) computes the logarithmic integral; i.e., the ! integral from 0.0 to X of (1.0/ln(t))dt. ! !***REFERENCES (NONE) !***ROUTINES CALLED EI, XERMSG !***REVISION HISTORY (YYMMDD) ! 770601 DATE WRITTEN ! 890531 Changed all specific intrinsics to generic. (WRB) ! 890531 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) ! 900326 Removed duplicate information from DESCRIPTION section. ! (WRB) !***END PROLOGUE ALI !***FIRST EXECUTABLE STATEMENT ALI if (X <= 0.0) call XERMSG ('SLATEC', 'ALI', & 'LOG INTEGRAL UNDEFINED FOR X LE 0', 1, 2) if (X == 1.0) call XERMSG ('SLATEC', 'ALI', & 'LOG INTEGRAL UNDEFINED FOR X = 1', 2, 2) ! ALI = EI (LOG(X) ) ! return end function ALNGAM (X) ! !! ALNGAM computes the logarithm of the absolute value of the Gamma function. ! !***LIBRARY SLATEC (FNLIB) !***CATEGORY C7A !***TYPE SINGLE PRECISION (ALNGAM-S, DLNGAM-D, CLNGAM-C) !***KEYWORDS ABSOLUTE VALUE, COMPLETE GAMMA FUNCTION, FNLIB, LOGARITHM, ! SPECIAL FUNCTIONS !***AUTHOR Fullerton, W., (LANL) !***DESCRIPTION ! ! ALNGAM(X) computes the logarithm of the absolute value of the ! gamma function at X. ! !***REFERENCES (NONE) !***ROUTINES CALLED GAMMA, R1MACH, R9LGMC, XERMSG !***REVISION HISTORY (YYMMDD) ! 770601 DATE WRITTEN ! 890531 Changed all specific intrinsics to generic. (WRB) ! 890531 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) ! 900326 Removed duplicate information from DESCRIPTION section. ! (WRB) ! 900727 Added EXTERNAL statement. (WRB) !***END PROLOGUE ALNGAM LOGICAL FIRST EXTERNAL GAMMA SAVE SQ2PIL, SQPI2L, PI, XMAX, DXREL, FIRST DATA SQ2PIL / 0.91893853320467274E0/ DATA SQPI2L / 0.22579135264472743E0/ DATA PI / 3.14159265358979324E0/ DATA FIRST /.TRUE./ !***FIRST EXECUTABLE STATEMENT ALNGAM if (FIRST) THEN XMAX = R1MACH(2)/LOG(R1MACH(2)) DXREL = SQRT (R1MACH(4)) end if FIRST = .FALSE. ! Y = ABS(X) if (Y > 10.0) go to 20 ! ! LOG (ABS (GAMMA(X))) FOR ABS(X) <= 10.0 ! ALNGAM = LOG (ABS (GAMMA(X))) return ! ! LOG (ABS (GAMMA(X))) FOR ABS(X) > 10.0 ! 20 if (Y > XMAX) call XERMSG ('SLATEC', 'ALNGAM', & 'ABS(X) SO BIG ALNGAM OVERFLOWS', 2, 2) ! if (X > 0.) ALNGAM = SQ2PIL + (X-0.5)*LOG(X) - X + R9LGMC(Y) if (X > 0.) RETURN ! SINPIY = ABS (SIN(PI*Y)) if (SINPIY == 0.) call XERMSG ('SLATEC', 'ALNGAM', & 'X IS A NEGATIVE INTEGER', 3, 2) ! if (ABS((X-AINT(X-0.5))/X) < DXREL) call XERMSG ('SLATEC', & 'ALNGAM', 'ANSWER LT HALF PRECISION BECAUSE X TOO NEAR ' // & 'NEGATIVE INTEGER', 1, 1) ! ALNGAM = SQPI2L + (X-0.5)*LOG(Y) - X - LOG(SINPIY) - R9LGMC(Y) return ! end function ALNREL (X) ! !! ALNREL evaluates ln(1+X) accurate in the sense of relative error. ! !***LIBRARY SLATEC (FNLIB) !***CATEGORY C4B !***TYPE SINGLE PRECISION (ALNREL-S, DLNREL-D, CLNREL-C) !***KEYWORDS ELEMENTARY FUNCTIONS, FNLIB, LOGARITHM !***AUTHOR Fullerton, W., (LANL) !***DESCRIPTION ! ! ALNREL(X) evaluates ln(1+X) accurately in the sense of relative ! error when X is very small. This routine must be used to ! maintain relative error accuracy whenever X is small and ! accurately known. ! ! Series for ALNR on the interval -3.75000D-01 to 3.75000D-01 ! with weighted error 1.93E-17 ! log weighted error 16.72 ! significant figures required 16.44 ! decimal places required 17.40 ! !***REFERENCES (NONE) !***ROUTINES CALLED CSEVL, INITS, R1MACH, XERMSG !***REVISION HISTORY (YYMMDD) ! 770401 DATE WRITTEN ! 890531 Changed all specific intrinsics to generic. (WRB) ! 890531 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) ! 900326 Removed duplicate information from DESCRIPTION section. ! (WRB) !***END PROLOGUE ALNREL DIMENSION ALNRCS(23) LOGICAL FIRST SAVE ALNRCS, NLNREL, XMIN, FIRST DATA ALNRCS( 1) / 1.0378693562743770E0 / DATA ALNRCS( 2) / -.13364301504908918E0 / DATA ALNRCS( 3) / .019408249135520563E0 / DATA ALNRCS( 4) / -.003010755112753577E0 / DATA ALNRCS( 5) / .000486946147971548E0 / DATA ALNRCS( 6) / -.000081054881893175E0 / DATA ALNRCS( 7) / .000013778847799559E0 / DATA ALNRCS( 8) / -.000002380221089435E0 / DATA ALNRCS( 9) / .000000416404162138E0 / DATA ALNRCS(10) / -.000000073595828378E0 / DATA ALNRCS(11) / .000000013117611876E0 / DATA ALNRCS(12) / -.000000002354670931E0 / DATA ALNRCS(13) / .000000000425227732E0 / DATA ALNRCS(14) / -.000000000077190894E0 / DATA ALNRCS(15) / .000000000014075746E0 / DATA ALNRCS(16) / -.000000000002576907E0 / DATA ALNRCS(17) / .000000000000473424E0 / DATA ALNRCS(18) / -.000000000000087249E0 / DATA ALNRCS(19) / .000000000000016124E0 / DATA ALNRCS(20) / -.000000000000002987E0 / DATA ALNRCS(21) / .000000000000000554E0 / DATA ALNRCS(22) / -.000000000000000103E0 / DATA ALNRCS(23) / .000000000000000019E0 / DATA FIRST /.TRUE./ !***FIRST EXECUTABLE STATEMENT ALNREL if (FIRST) THEN NLNREL = INITS (ALNRCS, 23, 0.1*R1MACH(3)) XMIN = -1.0 + SQRT(R1MACH(4)) end if FIRST = .FALSE. ! if (X <= (-1.0)) call XERMSG ('SLATEC', 'ALNREL', 'X IS LE -1', & 2, 2) if (X < XMIN) call XERMSG ('SLATEC', 'ALNREL', & 'ANSWER LT HALF PRECISION BECAUSE X TOO NEAR -1', 1, 1) ! if (ABS(X) <= 0.375) ALNREL = X*(1. - & X*CSEVL (X/.375, ALNRCS, NLNREL)) if (ABS(X) > 0.375) ALNREL = LOG (1.0+X) ! return end function ASINH (X) ! !! ASINH computes the arc hyperbolic sine. ! !***LIBRARY SLATEC (FNLIB) !***CATEGORY C4C !***TYPE SINGLE PRECISION (ASINH-S, DASINH-D, CASINH-C) !***KEYWORDS ARC HYPERBOLIC SINE, ASINH, ELEMENTARY FUNCTIONS, FNLIB, ! INVERSE HYPERBOLIC SINE !***AUTHOR Fullerton, W., (LANL) !***DESCRIPTION ! ! ASINH(X) computes the arc hyperbolic sine of X. ! ! Series for ASNH on the interval 0. to 1.00000D+00 ! with weighted error 2.19E-17 ! log weighted error 16.66 ! significant figures required 15.60 ! decimal places required 17.31 ! !***REFERENCES (NONE) !***ROUTINES CALLED CSEVL, INITS, R1MACH !***REVISION HISTORY (YYMMDD) ! 770401 DATE WRITTEN ! 890531 Changed all specific intrinsics to generic. (WRB) ! 890531 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) !***END PROLOGUE ASINH DIMENSION ASNHCS(20) LOGICAL FIRST SAVE ALN2, ASNHCS, NTERMS, XMAX, SQEPS, FIRST DATA ALN2 /0.69314718055994530942E0/ DATA ASNHCS( 1) / -.12820039911738186E0 / DATA ASNHCS( 2) / -.058811761189951768E0 / DATA ASNHCS( 3) / .004727465432212481E0 / DATA ASNHCS( 4) / -.000493836316265361E0 / DATA ASNHCS( 5) / .000058506207058557E0 / DATA ASNHCS( 6) / -.000007466998328931E0 / DATA ASNHCS( 7) / .000001001169358355E0 / DATA ASNHCS( 8) / -.000000139035438587E0 / DATA ASNHCS( 9) / .000000019823169483E0 / DATA ASNHCS(10) / -.000000002884746841E0 / DATA ASNHCS(11) / .000000000426729654E0 / DATA ASNHCS(12) / -.000000000063976084E0 / DATA ASNHCS(13) / .000000000009699168E0 / DATA ASNHCS(14) / -.000000000001484427E0 / DATA ASNHCS(15) / .000000000000229037E0 / DATA ASNHCS(16) / -.000000000000035588E0 / DATA ASNHCS(17) / .000000000000005563E0 / DATA ASNHCS(18) / -.000000000000000874E0 / DATA ASNHCS(19) / .000000000000000138E0 / DATA ASNHCS(20) / -.000000000000000021E0 / DATA FIRST /.TRUE./ !***FIRST EXECUTABLE STATEMENT ASINH if (FIRST) THEN NTERMS = INITS (ASNHCS, 20, 0.1*R1MACH(3)) SQEPS = SQRT (R1MACH(3)) XMAX = 1.0/SQEPS end if FIRST = .FALSE. ! Y = ABS(X) if (Y > 1.0) go to 20 ! ASINH = X if (Y > SQEPS) ASINH = X*(1.0 + CSEVL (2.*X*X-1., ASNHCS,NTERMS)) return ! 20 if (Y < XMAX) ASINH = LOG (Y + SQRT(Y**2+1.)) if (Y >= XMAX) ASINH = ALN2 + LOG(Y) ASINH = SIGN (ASINH, X) ! return end subroutine ASYIK (X, FNU, KODE, FLGIK, RA, ARG, IN, Y) ! !! ASYIK is subsidiary to BESI and BESK. ! !***LIBRARY SLATEC !***TYPE SINGLE PRECISION (ASYIK-S, DASYIK-D) !***AUTHOR Amos, D. E., (SNLA) !***DESCRIPTION ! ! ASYIK computes Bessel functions I and K ! for arguments X > 0.0 and orders FNU >= 35 ! on FLGIK = 1 and FLGIK = -1 respectively. ! ! INPUT ! ! X - argument, X > 0.0E0 ! FNU - order of first Bessel function ! KODE - a parameter to indicate the scaling option ! KODE=1 returns Y(I)= I/SUB(FNU+I-1)/(X), I=1,IN ! or Y(I)= K/SUB(FNU+I-1)/(X), I=1,IN ! on FLGIK = 1.0E0 or FLGIK = -1.0E0 ! KODE=2 returns Y(I)=EXP(-X)*I/SUB(FNU+I-1)/(X), I=1,IN ! or Y(I)=EXP( X)*K/SUB(FNU+I-1)/(X), I=1,IN ! on FLGIK = 1.0E0 or FLGIK = -1.0E0 ! FLGIK - selection parameter for I or K function ! FLGIK = 1.0E0 gives the I function ! FLGIK = -1.0E0 gives the K function ! RA - SQRT(1.+Z*Z), Z=X/FNU ! ARG - argument of the leading exponential ! IN - number of functions desired, IN=1 or 2 ! ! OUTPUT ! ! Y - a vector whose first in components contain the sequence ! ! Abstract ! ASYIK implements the uniform asymptotic expansion of ! the I and K Bessel functions for FNU >= 35 and real ! X > 0.0E0. The forms are identical except for a change ! in sign of some of the terms. This change in sign is ! accomplished by means of the flag FLGIK = 1 or -1. ! !***SEE ALSO BESI, BESK !***ROUTINES CALLED R1MACH !***REVISION HISTORY (YYMMDD) ! 750101 DATE WRITTEN ! 890531 Changed all specific intrinsics to generic. (WRB) ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900328 Added TYPE section. (WRB) ! 910408 Updated the AUTHOR section. (WRB) !***END PROLOGUE ASYIK ! INTEGER IN, J, JN, K, KK, KODE, L REAL AK,AP,ARG,C, COEF,CON,ETX,FLGIK,FN, FNU,GLN,RA,S1,S2, & T, TOL, T2, X, Y, Z REAL R1MACH DIMENSION Y(*), C(65), CON(2) SAVE CON, C DATA CON(1), CON(2) / & 3.98942280401432678E-01, 1.25331413731550025E+00/ DATA C(1), C(2), C(3), C(4), C(5), C(6), C(7), C(8), C(9), C(10), & C(11), C(12), C(13), C(14), C(15), C(16), C(17), C(18), & C(19), C(20), C(21), C(22), C(23), C(24)/ & -2.08333333333333E-01, 1.25000000000000E-01, & 3.34201388888889E-01, -4.01041666666667E-01, & 7.03125000000000E-02, -1.02581259645062E+00, & 1.84646267361111E+00, -8.91210937500000E-01, & 7.32421875000000E-02, 4.66958442342625E+00, & -1.12070026162230E+01, 8.78912353515625E+00, & -2.36408691406250E+00, 1.12152099609375E-01, & -2.82120725582002E+01, 8.46362176746007E+01, & -9.18182415432400E+01, 4.25349987453885E+01, & -7.36879435947963E+00, 2.27108001708984E-01, & 2.12570130039217E+02, -7.65252468141182E+02, & 1.05999045252800E+03, -6.99579627376133E+02/ DATA C(25), C(26), C(27), C(28), C(29), C(30), C(31), C(32), & C(33), C(34), C(35), C(36), C(37), C(38), C(39), C(40), & C(41), C(42), C(43), C(44), C(45), C(46), C(47), C(48)/ & 2.18190511744212E+02, -2.64914304869516E+01, & 5.72501420974731E-01, -1.91945766231841E+03, & 8.06172218173731E+03, -1.35865500064341E+04, & 1.16553933368645E+04, -5.30564697861340E+03, & 1.20090291321635E+03, -1.08090919788395E+02, & 1.72772750258446E+00, 2.02042913309661E+04, & -9.69805983886375E+04, 1.92547001232532E+05, & -2.03400177280416E+05, 1.22200464983017E+05, & -4.11926549688976E+04, 7.10951430248936E+03, & -4.93915304773088E+02, 6.07404200127348E+00, & -2.42919187900551E+05, 1.31176361466298E+06, & -2.99801591853811E+06, 3.76327129765640E+06/ DATA C(49), C(50), C(51), C(52), C(53), C(54), C(55), C(56), & C(57), C(58), C(59), C(60), C(61), C(62), C(63), C(64), & C(65)/ & -2.81356322658653E+06, 1.26836527332162E+06, & -3.31645172484564E+05, 4.52187689813627E+04, & -2.49983048181121E+03, 2.43805296995561E+01, & 3.28446985307204E+06, -1.97068191184322E+07, & 5.09526024926646E+07, -7.41051482115327E+07, & 6.63445122747290E+07, -3.75671766607634E+07, & 1.32887671664218E+07, -2.78561812808645E+06, & 3.08186404612662E+05, -1.38860897537170E+04, & 1.10017140269247E+02/ !***FIRST EXECUTABLE STATEMENT ASYIK TOL = R1MACH(3) TOL = MAX(TOL,1.0E-15) FN = FNU Z = (3.0E0-FLGIK)/2.0E0 KK = INT(Z) DO 50 JN=1,IN if (JN == 1) go to 10 FN = FN - FLGIK Z = X/FN RA = SQRT(1.0E0+Z*Z) GLN = LOG((1.0E0+RA)/Z) ETX = KODE - 1 T = RA*(1.0E0-ETX) + ETX/(Z+RA) ARG = FN*(T-GLN)*FLGIK 10 COEF = EXP(ARG) T = 1.0E0/RA T2 = T*T T = T/FN T = SIGN(T,FLGIK) S2 = 1.0E0 AP = 1.0E0 L = 0 DO 30 K=2,11 L = L + 1 S1 = C(L) DO 20 J=2,K L = L + 1 S1 = S1*T2 + C(L) 20 CONTINUE AP = AP*T AK = AP*S1 S2 = S2 + AK if (MAX(ABS(AK),ABS(AP)) < TOL) go to 40 30 CONTINUE 40 CONTINUE T = ABS(T) Y(JN) = S2*COEF*SQRT(T)*CON(KK) 50 CONTINUE return end subroutine ASYJY (FUNJY, X, FNU, FLGJY, IN, Y, WK, IFLW) ! !! ASYJY is subsidiary to BESJ and BESY. ! !***LIBRARY SLATEC !***TYPE SINGLE PRECISION (ASYJY-S, DASYJY-D) !***AUTHOR Amos, D. E., (SNLA) !***DESCRIPTION ! ! ASYJY computes Bessel functions J and Y ! for arguments X > 0.0 and orders FNU >= 35.0 ! on FLGJY = 1 and FLGJY = -1 respectively ! ! INPUT ! ! FUNJY - external function JAIRY or YAIRY ! X - argument, X > 0.0E0 ! FNU - order of the first Bessel function ! FLGJY - selection flag ! FLGJY = 1.0E0 gives the J function ! FLGJY = -1.0E0 gives the Y function ! IN - number of functions desired, IN = 1 or 2 ! ! OUTPUT ! ! Y - a vector whose first in components contain the sequence ! IFLW - a flag indicating underflow or overflow ! return variables for BESJ only ! WK(1) = 1 - (X/FNU)**2 = W**2 ! WK(2) = SQRT(ABS(WK(1))) ! WK(3) = ABS(WK(2) - ATAN(WK(2))) or ! ABS(LN((1 + WK(2))/(X/FNU)) - WK(2)) ! = ABS((2/3)*ZETA**(3/2)) ! WK(4) = FNU*WK(3) ! WK(5) = (1.5*WK(3)*FNU)**(1/3) = SQRT(ZETA)*FNU**(1/3) ! WK(6) = SIGN(1.,W**2)*WK(5)**2 = SIGN(1.,W**2)*ZETA*FNU**(2/3) ! WK(7) = FNU**(1/3) ! ! Abstract ! ASYJY implements the uniform asymptotic expansion of ! the J and Y Bessel functions for FNU >= 35 and real ! X > 0.0E0. The forms are identical except for a change ! in sign of some of the terms. This change in sign is ! accomplished by means of the flag FLGJY = 1 or -1. On ! FLGJY = 1 the AIRY functions AI(X) and DAI(X) are ! supplied by the external function JAIRY, and on ! FLGJY = -1 the AIRY functions BI(X) and DBI(X) are ! supplied by the external function YAIRY. ! !***SEE ALSO BESJ, BESY !***ROUTINES CALLED I1MACH, R1MACH !***REVISION HISTORY (YYMMDD) ! 750101 DATE WRITTEN ! 890531 Changed all specific intrinsics to generic. (WRB) ! 891009 Removed unreferenced variable. (WRB) ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900328 Added TYPE section. (WRB) ! 910408 Updated the AUTHOR section. (WRB) !***END PROLOGUE ASYJY INTEGER I, IFLW, IN, J, JN,JR,JU,K, KB,KLAST,KMAX,KP1, KS, KSP1, & KSTEMP, L, LR, LRP1, ISETA, ISETB INTEGER I1MACH REAL ABW2, AKM, ALFA, ALFA1, ALFA2, AP, AR, ASUM, AZ, & BETA, BETA1, BETA2, BETA3, BR, BSUM, C, CON1, CON2, & CON548,CR,CRZ32, DFI,ELIM, DR,FI, FLGJY, FN, FNU, & FN2, GAMA, PHI, RCZ, RDEN, RELB, RFN2, RTZ, RZDEN, & SA, SB, SUMA, SUMB, S1, TA, TAU, TB, TFN, TOL, TOLS, T2, UPOL, & WK, X, XX, Y, Z, Z32 REAL R1MACH DIMENSION Y(*), WK(*), C(65) DIMENSION ALFA(26,4), BETA(26,5) DIMENSION ALFA1(26,2), ALFA2(26,2) DIMENSION BETA1(26,2), BETA2(26,2), BETA3(26,1) DIMENSION GAMA(26), KMAX(5), AR(8), BR(10), UPOL(10) DIMENSION CR(10), DR(10) EQUIVALENCE (ALFA(1,1),ALFA1(1,1)) EQUIVALENCE (ALFA(1,3),ALFA2(1,1)) EQUIVALENCE (BETA(1,1),BETA1(1,1)) EQUIVALENCE (BETA(1,3),BETA2(1,1)) EQUIVALENCE (BETA(1,5),BETA3(1,1)) SAVE TOLS, CON1, CON2, CON548, AR, BR, C, ALFA1, ALFA2, & BETA1, BETA2, BETA3, GAMA DATA TOLS /-6.90775527898214E+00/ DATA CON1,CON2,CON548/ & 6.66666666666667E-01, 3.33333333333333E-01, 1.04166666666667E-01/ DATA AR(1), AR(2), AR(3), AR(4), AR(5), AR(6), AR(7), & AR(8) / 8.35503472222222E-02, 1.28226574556327E-01, & 2.91849026464140E-01, 8.81627267443758E-01, 3.32140828186277E+00, & 1.49957629868626E+01, 7.89230130115865E+01, 4.74451538868264E+02/ DATA BR(1), BR(2), BR(3), BR(4), BR(5), BR(6), BR(7), BR(8), & BR(9), BR(10) /-1.45833333333333E-01,-9.87413194444444E-02, & -1.43312053915895E-01,-3.17227202678414E-01,-9.42429147957120E-01, & -3.51120304082635E+00,-1.57272636203680E+01,-8.22814390971859E+01, & -4.92355370523671E+02,-3.31621856854797E+03/ DATA C(1), C(2), C(3), C(4), C(5), C(6), C(7), C(8), C(9), C(10), & C(11), C(12), C(13), C(14), C(15), C(16), C(17), C(18), & C(19), C(20), C(21), C(22), C(23), C(24)/ & -2.08333333333333E-01, 1.25000000000000E-01, & 3.34201388888889E-01, -4.01041666666667E-01, & 7.03125000000000E-02, -1.02581259645062E+00, & 1.84646267361111E+00, -8.91210937500000E-01, & 7.32421875000000E-02, 4.66958442342625E+00, & -1.12070026162230E+01, 8.78912353515625E+00, & -2.36408691406250E+00, 1.12152099609375E-01, & -2.82120725582002E+01, 8.46362176746007E+01, & -9.18182415432400E+01, 4.25349987453885E+01, & -7.36879435947963E+00, 2.27108001708984E-01, & 2.12570130039217E+02, -7.65252468141182E+02, & 1.05999045252800E+03, -6.99579627376133E+02/ DATA C(25), C(26), C(27), C(28), C(29), C(30), C(31), C(32), & C(33), C(34), C(35), C(36), C(37), C(38), C(39), C(40), & C(41), C(42), C(43), C(44), C(45), C(46), C(47), C(48)/ & 2.18190511744212E+02, -2.64914304869516E+01, & 5.72501420974731E-01, -1.91945766231841E+03, & 8.06172218173731E+03, -1.35865500064341E+04, & 1.16553933368645E+04, -5.30564697861340E+03, & 1.20090291321635E+03, -1.08090919788395E+02, & 1.72772750258446E+00, 2.02042913309661E+04, & -9.69805983886375E+04, 1.92547001232532E+05, & -2.03400177280416E+05, 1.22200464983017E+05, & -4.11926549688976E+04, 7.10951430248936E+03, & -4.93915304773088E+02, 6.07404200127348E+00, & -2.42919187900551E+05, 1.31176361466298E+06, & -2.99801591853811E+06, 3.76327129765640E+06/ DATA C(49), C(50), C(51), C(52), C(53), C(54), C(55), C(56), & C(57), C(58), C(59), C(60), C(61), C(62), C(63), C(64), & C(65)/ & -2.81356322658653E+06, 1.26836527332162E+06, & -3.31645172484564E+05, 4.52187689813627E+04, & -2.49983048181121E+03, 2.43805296995561E+01, & 3.28446985307204E+06, -1.97068191184322E+07, & 5.09526024926646E+07, -7.41051482115327E+07, & 6.63445122747290E+07, -3.75671766607634E+07, & 1.32887671664218E+07, -2.78561812808645E+06, & 3.08186404612662E+05, -1.38860897537170E+04, & 1.10017140269247E+02/ DATA ALFA1(1,1), ALFA1(2,1), ALFA1(3,1), ALFA1(4,1), ALFA1(5,1), & ALFA1(6,1), ALFA1(7,1), ALFA1(8,1), ALFA1(9,1), ALFA1(10,1), & ALFA1(11,1),ALFA1(12,1),ALFA1(13,1),ALFA1(14,1),ALFA1(15,1), & ALFA1(16,1),ALFA1(17,1),ALFA1(18,1),ALFA1(19,1),ALFA1(20,1), & ALFA1(21,1),ALFA1(22,1),ALFA1(23,1),ALFA1(24,1),ALFA1(25,1), & ALFA1(26,1) /-4.44444444444444E-03,-9.22077922077922E-04, & -8.84892884892885E-05, 1.65927687832450E-04, 2.46691372741793E-04, & 2.65995589346255E-04, 2.61824297061501E-04, 2.48730437344656E-04, & 2.32721040083232E-04, 2.16362485712365E-04, 2.00738858762752E-04, & 1.86267636637545E-04, 1.73060775917876E-04, 1.61091705929016E-04, & 1.50274774160908E-04, 1.40503497391270E-04, 1.31668816545923E-04, & 1.23667445598253E-04, 1.16405271474738E-04, 1.09798298372713E-04, & 1.03772410422993E-04, 9.82626078369363E-05, 9.32120517249503E-05, & 8.85710852478712E-05, 8.42963105715700E-05, 8.03497548407791E-05/ DATA ALFA1(1,2), ALFA1(2,2), ALFA1(3,2), ALFA1(4,2), ALFA1(5,2), & ALFA1(6,2), ALFA1(7,2), ALFA1(8,2), ALFA1(9,2), ALFA1(10,2), & ALFA1(11,2),ALFA1(12,2),ALFA1(13,2),ALFA1(14,2),ALFA1(15,2), & ALFA1(16,2),ALFA1(17,2),ALFA1(18,2),ALFA1(19,2),ALFA1(20,2), & ALFA1(21,2),ALFA1(22,2),ALFA1(23,2),ALFA1(24,2),ALFA1(25,2), & ALFA1(26,2) / 6.93735541354589E-04, 2.32241745182922E-04, & -1.41986273556691E-05,-1.16444931672049E-04,-1.50803558053049E-04, & -1.55121924918096E-04,-1.46809756646466E-04,-1.33815503867491E-04, & -1.19744975684254E-04,-1.06184319207974E-04,-9.37699549891194E-05, & -8.26923045588193E-05,-7.29374348155221E-05,-6.44042357721016E-05, & -5.69611566009369E-05,-5.04731044303562E-05,-4.48134868008883E-05, & -3.98688727717599E-05,-3.55400532972042E-05,-3.17414256609022E-05, & -2.83996793904175E-05,-2.54522720634871E-05,-2.28459297164725E-05, & -2.05352753106481E-05,-1.84816217627666E-05,-1.66519330021394E-05/ DATA ALFA2(1,1), ALFA2(2,1), ALFA2(3,1), ALFA2(4,1), ALFA2(5,1), & ALFA2(6,1), ALFA2(7,1), ALFA2(8,1), ALFA2(9,1), ALFA2(10,1), & ALFA2(11,1),ALFA2(12,1),ALFA2(13,1),ALFA2(14,1),ALFA2(15,1), & ALFA2(16,1),ALFA2(17,1),ALFA2(18,1),ALFA2(19,1),ALFA2(20,1), & ALFA2(21,1),ALFA2(22,1),ALFA2(23,1),ALFA2(24,1),ALFA2(25,1), & ALFA2(26,1) /-3.54211971457744E-04,-1.56161263945159E-04, & 3.04465503594936E-05, 1.30198655773243E-04, 1.67471106699712E-04, & 1.70222587683593E-04, 1.56501427608595E-04, 1.36339170977445E-04, & 1.14886692029825E-04, 9.45869093034688E-05, 7.64498419250898E-05, & 6.07570334965197E-05, 4.74394299290509E-05, 3.62757512005344E-05, & 2.69939714979225E-05, 1.93210938247939E-05, 1.30056674793963E-05, & 7.82620866744497E-06, 3.59257485819352E-06, 1.44040049814252E-07, & -2.65396769697939E-06,-4.91346867098486E-06,-6.72739296091248E-06, & -8.17269379678658E-06,-9.31304715093561E-06,-1.02011418798016E-05/ DATA ALFA2(1,2), ALFA2(2,2), ALFA2(3,2), ALFA2(4,2), ALFA2(5,2), & ALFA2(6,2), ALFA2(7,2), ALFA2(8,2), ALFA2(9,2), ALFA2(10,2), & ALFA2(11,2),ALFA2(12,2),ALFA2(13,2),ALFA2(14,2),ALFA2(15,2), & ALFA2(16,2),ALFA2(17,2),ALFA2(18,2),ALFA2(19,2),ALFA2(20,2), & ALFA2(21,2),ALFA2(22,2),ALFA2(23,2),ALFA2(24,2),ALFA2(25,2), & ALFA2(26,2) / 3.78194199201773E-04, 2.02471952761816E-04, & -6.37938506318862E-05,-2.38598230603006E-04,-3.10916256027362E-04, & -3.13680115247576E-04,-2.78950273791323E-04,-2.28564082619141E-04, & -1.75245280340847E-04,-1.25544063060690E-04,-8.22982872820208E-05, & -4.62860730588116E-05,-1.72334302366962E-05, 5.60690482304602E-06, & 2.31395443148287E-05, 3.62642745856794E-05, 4.58006124490189E-05, & 5.24595294959114E-05, 5.68396208545815E-05, 5.94349820393104E-05, & 6.06478527578422E-05, 6.08023907788436E-05, 6.01577894539460E-05, & 5.89199657344698E-05, 5.72515823777593E-05, 5.52804375585853E-05/ DATA BETA1(1,1), BETA1(2,1), BETA1(3,1), BETA1(4,1), BETA1(5,1), & BETA1(6,1), BETA1(7,1), BETA1(8,1), BETA1(9,1), BETA1(10,1), & BETA1(11,1),BETA1(12,1),BETA1(13,1),BETA1(14,1),BETA1(15,1), & BETA1(16,1),BETA1(17,1),BETA1(18,1),BETA1(19,1),BETA1(20,1), & BETA1(21,1),BETA1(22,1),BETA1(23,1),BETA1(24,1),BETA1(25,1), & BETA1(26,1) / 1.79988721413553E-02, 5.59964911064388E-03, & 2.88501402231133E-03, 1.80096606761054E-03, 1.24753110589199E-03, & 9.22878876572938E-04, 7.14430421727287E-04, 5.71787281789705E-04, & 4.69431007606482E-04, 3.93232835462917E-04, 3.34818889318298E-04, & 2.88952148495752E-04, 2.52211615549573E-04, 2.22280580798883E-04, & 1.97541838033063E-04, 1.76836855019718E-04, 1.59316899661821E-04, & 1.44347930197334E-04, 1.31448068119965E-04, 1.20245444949303E-04, & 1.10449144504599E-04, 1.01828770740567E-04, 9.41998224204238E-05, & 8.74130545753834E-05, 8.13466262162801E-05, 7.59002269646219E-05/ DATA BETA1(1,2), BETA1(2,2), BETA1(3,2), BETA1(4,2), BETA1(5,2), & BETA1(6,2), BETA1(7,2), BETA1(8,2), BETA1(9,2), BETA1(10,2), & BETA1(11,2),BETA1(12,2),BETA1(13,2),BETA1(14,2),BETA1(15,2), & BETA1(16,2),BETA1(17,2),BETA1(18,2),BETA1(19,2),BETA1(20,2), & BETA1(21,2),BETA1(22,2),BETA1(23,2),BETA1(24,2),BETA1(25,2), & BETA1(26,2) /-1.49282953213429E-03,-8.78204709546389E-04, & -5.02916549572035E-04,-2.94822138512746E-04,-1.75463996970783E-04, & -1.04008550460816E-04,-5.96141953046458E-05,-3.12038929076098E-05, & -1.26089735980230E-05,-2.42892608575730E-07, 8.05996165414274E-06, & 1.36507009262147E-05, 1.73964125472926E-05, 1.98672978842134E-05, & 2.14463263790823E-05, 2.23954659232457E-05, 2.28967783814713E-05, & 2.30785389811178E-05, 2.30321976080909E-05, 2.28236073720349E-05, & 2.25005881105292E-05, 2.20981015361991E-05, 2.16418427448104E-05, & 2.11507649256221E-05, 2.06388749782171E-05, 2.01165241997082E-05/ DATA BETA2(1,1), BETA2(2,1), BETA2(3,1), BETA2(4,1), BETA2(5,1), & BETA2(6,1), BETA2(7,1), BETA2(8,1), BETA2(9,1), BETA2(10,1), & BETA2(11,1),BETA2(12,1),BETA2(13,1),BETA2(14,1),BETA2(15,1), & BETA2(16,1),BETA2(17,1),BETA2(18,1),BETA2(19,1),BETA2(20,1), & BETA2(21,1),BETA2(22,1),BETA2(23,1),BETA2(24,1),BETA2(25,1), & BETA2(26,1) / 5.52213076721293E-04, 4.47932581552385E-04, & 2.79520653992021E-04, 1.52468156198447E-04, 6.93271105657044E-05, & 1.76258683069991E-05,-1.35744996343269E-05,-3.17972413350427E-05, & -4.18861861696693E-05,-4.69004889379141E-05,-4.87665447413787E-05, & -4.87010031186735E-05,-4.74755620890087E-05,-4.55813058138628E-05, & -4.33309644511266E-05,-4.09230193157750E-05,-3.84822638603221E-05, & -3.60857167535411E-05,-3.37793306123367E-05,-3.15888560772110E-05, & -2.95269561750807E-05,-2.75978914828336E-05,-2.58006174666884E-05, & -2.41308356761280E-05,-2.25823509518346E-05,-2.11479656768913E-05/ DATA BETA2(1,2), BETA2(2,2), BETA2(3,2), BETA2(4,2), BETA2(5,2), & BETA2(6,2), BETA2(7,2), BETA2(8,2), BETA2(9,2), BETA2(10,2), & BETA2(11,2),BETA2(12,2),BETA2(13,2),BETA2(14,2),BETA2(15,2), & BETA2(16,2),BETA2(17,2),BETA2(18,2),BETA2(19,2),BETA2(20,2), & BETA2(21,2),BETA2(22,2),BETA2(23,2),BETA2(24,2),BETA2(25,2), & BETA2(26,2) /-4.74617796559960E-04,-4.77864567147321E-04, & -3.20390228067038E-04,-1.61105016119962E-04,-4.25778101285435E-05, & 3.44571294294968E-05, 7.97092684075675E-05, 1.03138236708272E-04, & 1.12466775262204E-04, 1.13103642108481E-04, 1.08651634848774E-04, & 1.01437951597662E-04, 9.29298396593364E-05, 8.40293133016090E-05, & 7.52727991349134E-05, 6.69632521975731E-05, 5.92564547323195E-05, & 5.22169308826976E-05, 4.58539485165361E-05, 4.01445513891487E-05, & 3.50481730031328E-05, 3.05157995034347E-05, 2.64956119950516E-05, & 2.29363633690998E-05, 1.97893056664022E-05, 1.70091984636413E-05/ DATA BETA3(1,1), BETA3(2,1), BETA3(3,1), BETA3(4,1), BETA3(5,1), & BETA3(6,1), BETA3(7,1), BETA3(8,1), BETA3(9,1), BETA3(10,1), & BETA3(11,1),BETA3(12,1),BETA3(13,1),BETA3(14,1),BETA3(15,1), & BETA3(16,1),BETA3(17,1),BETA3(18,1),BETA3(19,1),BETA3(20,1), & BETA3(21,1),BETA3(22,1),BETA3(23,1),BETA3(24,1),BETA3(25,1), & BETA3(26,1) / 7.36465810572578E-04, 8.72790805146194E-04, & 6.22614862573135E-04, 2.85998154194304E-04, 3.84737672879366E-06, & -1.87906003636972E-04,-2.97603646594555E-04,-3.45998126832656E-04, & -3.53382470916038E-04,-3.35715635775049E-04,-3.04321124789040E-04, & -2.66722723047613E-04,-2.27654214122820E-04,-1.89922611854562E-04, & -1.55058918599094E-04,-1.23778240761874E-04,-9.62926147717644E-05, & -7.25178327714425E-05,-5.22070028895634E-05,-3.50347750511901E-05, & -2.06489761035552E-05,-8.70106096849767E-06, 1.13698686675100E-06, & 9.16426474122779E-06, 1.56477785428873E-05, 2.08223629482467E-05/ DATA GAMA(1), GAMA(2), GAMA(3), GAMA(4), GAMA(5), & GAMA(6), GAMA(7), GAMA(8), GAMA(9), GAMA(10), & GAMA(11), GAMA(12), GAMA(13), GAMA(14), GAMA(15), & GAMA(16), GAMA(17), GAMA(18), GAMA(19), GAMA(20), & GAMA(21), GAMA(22), GAMA(23), GAMA(24), GAMA(25), & GAMA(26) / 6.29960524947437E-01, 2.51984209978975E-01, & 1.54790300415656E-01, 1.10713062416159E-01, 8.57309395527395E-02, & 6.97161316958684E-02, 5.86085671893714E-02, 5.04698873536311E-02, & 4.42600580689155E-02, 3.93720661543510E-02, 3.54283195924455E-02, & 3.21818857502098E-02, 2.94646240791158E-02, 2.71581677112934E-02, & 2.51768272973862E-02, 2.34570755306079E-02, 2.19508390134907E-02, & 2.06210828235646E-02, 1.94388240897881E-02, 1.83810633800683E-02, & 1.74293213231963E-02, 1.65685837786612E-02, 1.57865285987918E-02, & 1.50729501494096E-02, 1.44193250839955E-02, 1.38184805735342E-02/ !***FIRST EXECUTABLE STATEMENT ASYJY TA = R1MACH(3) TOL = MAX(TA,1.0E-15) TB = R1MACH(5) JU = I1MACH(12) if ( FLGJY == 1.0E0) go to 6 JR = I1MACH(11) ELIM = -2.303E0*TB*(JU+JR) go to 7 6 CONTINUE ELIM = -2.303E0*(TB*JU+3.0E0) 7 CONTINUE FN = FNU IFLW = 0 DO 170 JN=1,IN XX = X/FN WK(1) = 1.0E0 - XX*XX ABW2 = ABS(WK(1)) WK(2) = SQRT(ABW2) WK(7) = FN**CON2 if (ABW2 > 0.27750E0) go to 80 ! ! ASYMPTOTIC EXPANSION ! CASES NEAR X=FN, ABS(1.-(X/FN)**2) <= 0.2775 ! COEFFICIENTS OF ASYMPTOTIC EXPANSION BY SERIES ! ! ZETA AND TRUNCATION FOR A(ZETA) AND B(ZETA) SERIES ! ! KMAX IS TRUNCATION INDEX FOR A(ZETA) AND B(ZETA) SERIES=MAX(2,SA) ! SA = 0.0E0 if (ABW2 == 0.0E0) go to 10 SA = TOLS/LOG(ABW2) 10 SB = SA DO 20 I=1,5 AKM = MAX(SA,2.0E0) KMAX(I) = INT(AKM) SA = SA + SB 20 CONTINUE KB = KMAX(5) KLAST = KB - 1 SA = GAMA(KB) DO 30 K=1,KLAST KB = KB - 1 SA = SA*WK(1) + GAMA(KB) 30 CONTINUE Z = WK(1)*SA AZ = ABS(Z) RTZ = SQRT(AZ) WK(3) = CON1*AZ*RTZ WK(4) = WK(3)*FN WK(5) = RTZ*WK(7) WK(6) = -WK(5)*WK(5) if ( Z <= 0.0E0) go to 35 if ( WK(4) > ELIM) go to 75 WK(6) = -WK(6) 35 CONTINUE PHI = SQRT(SQRT(SA+SA+SA+SA)) ! ! B(ZETA) FOR S=0 ! KB = KMAX(5) KLAST = KB - 1 SB = BETA(KB,1) DO 40 K=1,KLAST KB = KB - 1 SB = SB*WK(1) + BETA(KB,1) 40 CONTINUE KSP1 = 1 FN2 = FN*FN RFN2 = 1.0E0/FN2 RDEN = 1.0E0 ASUM = 1.0E0 RELB = TOL*ABS(SB) BSUM = SB DO 60 KS=1,4 KSP1 = KSP1 + 1 RDEN = RDEN*RFN2 ! ! A(ZETA) AND B(ZETA) FOR S=1,2,3,4 ! KSTEMP = 5 - KS KB = KMAX(KSTEMP) KLAST = KB - 1 SA = ALFA(KB,KS) SB = BETA(KB,KSP1) DO 50 K=1,KLAST KB = KB - 1 SA = SA*WK(1) + ALFA(KB,KS) SB = SB*WK(1) + BETA(KB,KSP1) 50 CONTINUE TA = SA*RDEN TB = SB*RDEN ASUM = ASUM + TA BSUM = BSUM + TB if (ABS(TA) <= TOL .AND. ABS(TB) <= RELB) go to 70 60 CONTINUE 70 CONTINUE BSUM = BSUM/(FN*WK(7)) go to 160 ! 75 CONTINUE IFLW = 1 return ! 80 CONTINUE UPOL(1) = 1.0E0 TAU = 1.0E0/WK(2) T2 = 1.0E0/WK(1) if (WK(1) >= 0.0E0) go to 90 ! ! CASES FOR (X/FN) > SQRT(1.2775) ! WK(3) = ABS(WK(2)-ATAN(WK(2))) WK(4) = WK(3)*FN RCZ = -CON1/WK(4) Z32 = 1.5E0*WK(3) RTZ = Z32**CON2 WK(5) = RTZ*WK(7) WK(6) = -WK(5)*WK(5) go to 100 90 CONTINUE ! ! CASES FOR (X/FN) < SQRT(0.7225) ! WK(3) = ABS(LOG((1.0E0+WK(2))/XX)-WK(2)) WK(4) = WK(3)*FN RCZ = CON1/WK(4) if ( WK(4) > ELIM) go to 75 Z32 = 1.5E0*WK(3) RTZ = Z32**CON2 WK(7) = FN**CON2 WK(5) = RTZ*WK(7) WK(6) = WK(5)*WK(5) 100 CONTINUE PHI = SQRT((RTZ+RTZ)*TAU) TB = 1.0E0 ASUM = 1.0E0 TFN = TAU/FN RDEN=1.0E0/FN RFN2=RDEN*RDEN RDEN=1.0E0 UPOL(2) = (C(1)*T2+C(2))*TFN CRZ32 = CON548*RCZ BSUM = UPOL(2) + CRZ32 RELB = TOL*ABS(BSUM) AP = TFN KS = 0 KP1 = 2 RZDEN = RCZ L = 2 ISETA=0 ISETB=0 DO 140 LR=2,8,2 ! ! COMPUTE TWO U POLYNOMIALS FOR NEXT A(ZETA) AND B(ZETA) ! LRP1 = LR + 1 DO 120 K=LR,LRP1 KS = KS + 1 KP1 = KP1 + 1 L = L + 1 S1 = C(L) DO 110 J=2,KP1 L = L + 1 S1 = S1*T2 + C(L) 110 CONTINUE AP = AP*TFN UPOL(KP1) = AP*S1 CR(KS) = BR(KS)*RZDEN RZDEN = RZDEN*RCZ DR(KS) = AR(KS)*RZDEN 120 CONTINUE SUMA = UPOL(LRP1) SUMB = UPOL(LR+2) + UPOL(LRP1)*CRZ32 JU = LRP1 DO 130 JR=1,LR JU = JU - 1 SUMA = SUMA + CR(JR)*UPOL(JU) SUMB = SUMB + DR(JR)*UPOL(JU) 130 CONTINUE RDEN=RDEN*RFN2 TB = -TB if (WK(1) > 0.0E0) TB = ABS(TB) if (RDEN < TOL) go to 131 ASUM = ASUM + SUMA*TB BSUM = BSUM + SUMB*TB go to 140 131 if ( ISETA == 1) go to 132 if ( ABS(SUMA) < TOL) ISETA=1 ASUM=ASUM+SUMA*TB 132 if ( ISETB == 1) go to 133 if ( ABS(SUMB) < RELB) ISETB=1 BSUM=BSUM+SUMB*TB 133 if ( ISETA == 1 .AND. ISETB == 1) go to 150 140 CONTINUE 150 TB = WK(5) if (WK(1) > 0.0E0) TB = -TB BSUM = BSUM/TB ! 160 CONTINUE call FUNJY(WK(6), WK(5), WK(4), FI, DFI) TA=1.0E0/TOL TB=R1MACH(1)*TA*1.0E+3 if ( ABS(FI) > TB) go to 165 FI=FI*TA DFI=DFI*TA PHI=PHI*TOL 165 CONTINUE Y(JN) = FLGJY*PHI*(FI*ASUM+DFI*BSUM)/WK(7) FN = FN - FLGJY 170 CONTINUE return end function ATANH (X) ! !! ATANH computes the arc hyperbolic tangent. ! !***LIBRARY SLATEC (FNLIB) !***CATEGORY C4C !***TYPE SINGLE PRECISION (ATANH-S, DATANH-D, CATANH-C) !***KEYWORDS ARC HYPERBOLIC TANGENT, ATANH, ELEMENTARY FUNCTIONS, ! FNLIB, INVERSE HYPERBOLIC TANGENT !***AUTHOR Fullerton, W., (LANL) !***DESCRIPTION ! ! ATANH(X) computes the arc hyperbolic tangent of X. ! ! Series for ATNH on the interval 0. to 2.50000D-01 ! with weighted error 6.70E-18 ! log weighted error 17.17 ! significant figures required 16.01 ! decimal places required 17.76 ! !***REFERENCES (NONE) !***ROUTINES CALLED CSEVL, INITS, R1MACH, XERMSG !***REVISION HISTORY (YYMMDD) ! 770401 DATE WRITTEN ! 890531 Changed all specific intrinsics to generic. (WRB) ! 890531 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) ! 900326 Removed duplicate information from DESCRIPTION section. ! (WRB) !***END PROLOGUE ATANH DIMENSION ATNHCS(15) LOGICAL FIRST SAVE ATNHCS, NTERMS, DXREL, SQEPS, FIRST DATA ATNHCS( 1) / .094395102393195492E0 / DATA ATNHCS( 2) / .049198437055786159E0 / DATA ATNHCS( 3) / .002102593522455432E0 / DATA ATNHCS( 4) / .000107355444977611E0 / DATA ATNHCS( 5) / .000005978267249293E0 / DATA ATNHCS( 6) / .000000350506203088E0 / DATA ATNHCS( 7) / .000000021263743437E0 / DATA ATNHCS( 8) / .000000001321694535E0 / DATA ATNHCS( 9) / .000000000083658755E0 / DATA ATNHCS(10) / .000000000005370503E0 / DATA ATNHCS(11) / .000000000000348665E0 / DATA ATNHCS(12) / .000000000000022845E0 / DATA ATNHCS(13) / .000000000000001508E0 / DATA ATNHCS(14) / .000000000000000100E0 / DATA ATNHCS(15) / .000000000000000006E0 / DATA FIRST /.TRUE./ !***FIRST EXECUTABLE STATEMENT ATANH if (FIRST) THEN NTERMS = INITS (ATNHCS, 15, 0.1*R1MACH(3)) DXREL = SQRT (R1MACH(4)) SQEPS = SQRT (3.0*R1MACH(3)) end if FIRST = .FALSE. ! Y = ABS(X) if (Y >= 1.0) call XERMSG ('SLATEC', 'ATANH', 'ABS(X) GE 1', 2, & 2) ! if (1.0-Y < DXREL) call XERMSG ('SLATEC', 'ATANH', & 'ANSWER LT HALF PRECISION BECAUSE ABS(X) TOO NEAR 1', 1, 1) ! ATANH = X if (Y > SQEPS .AND. Y <= 0.5) ATANH = X*(1.0 + CSEVL (8.*X*X-1., & ATNHCS, NTERMS)) if (Y > 0.5) ATANH = 0.5*LOG((1.0+X)/(1.0-X)) ! return end subroutine AVINT (X, Y, N, XLO, XUP, ANS, IERR) ! !! AVINT integrates a function tabulated at arbitrarily spaced abscissas... ! using overlapping parabolas. ! !***LIBRARY SLATEC !***CATEGORY H2A1B2 !***TYPE SINGLE PRECISION (AVINT-S, DAVINT-D) !***KEYWORDS INTEGRATION, QUADRATURE, TABULATED DATA !***AUTHOR Jones, R. E., (SNLA) !***DESCRIPTION ! ! Abstract ! AVINT integrates a function tabulated at arbitrarily spaced ! abscissas. The limits of integration need not coincide ! with the tabulated abscissas. ! ! A method of overlapping parabolas fitted to the data is used ! provided that there are at least 3 abscissas between the ! limits of integration. AVINT also handles two special cases. ! If the limits of integration are equal, AVINT returns a result ! of zero regardless of the number of tabulated values. ! If there are only two function values, AVINT uses the ! trapezoid rule. ! ! Description of Parameters ! The user must dimension all arrays appearing in the call list ! X(N), Y(N). ! ! Input-- ! X - real array of abscissas, which must be in increasing ! order. ! Y - real array of functional values. i.e., Y(I)=FUNC(X(I)). ! N - the integer number of function values supplied. ! N >= 2 unless XLO = XUP. ! XLO - real lower limit of integration. ! XUP - real upper limit of integration. ! Must have XLO <= XUP. ! ! Output-- ! ANS - computed approximate value of integral ! IERR - a status code ! --normal code ! =1 means the requested integration was performed. ! --abnormal codes ! =2 means XUP was less than XLO. ! =3 means the number of X(I) between XLO and XUP ! (inclusive) was less than 3 and neither of the two ! special cases described in the Abstract occurred. ! No integration was performed. ! =4 means the restriction X(I+1) > X(I) was violated. ! =5 means the number N of function values was < 2. ! ANS is set to zero if IERR=2,3,4,or 5. ! ! AVINT is documented completely in SC-M-69-335 ! Original program from "Numerical Integration" by Davis & ! Rabinowitz. ! Adaptation and modifications for Sandia Mathematical Program ! Library by Rondall E. Jones. ! !***REFERENCES R. E. Jones, Approximate integrator of functions ! tabulated at arbitrarily spaced abscissas, ! Report SC-M-69-335, Sandia Laboratories, 1969. !***ROUTINES CALLED XERMSG !***REVISION HISTORY (YYMMDD) ! 690901 DATE WRITTEN ! 890831 Modified array declarations. (WRB) ! 890831 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) ! 900326 Removed duplicate information from DESCRIPTION section. ! (WRB) ! 920501 Reformatted the REFERENCES section. (WRB) !***END PROLOGUE AVINT ! DOUBLE PRECISION R3,RP5,SUM,SYL,SYL2,SYL3,SYU,SYU2,SYU3,X1,X2,X3 & ,X12,X13,X23,TERM1,TERM2,TERM3,A,B,C,CA,CB,CC DIMENSION X(*),Y(*) !***FIRST EXECUTABLE STATEMENT AVINT IERR=1 ANS =0.0 if (XLO-XUP) 3,100,200 3 if (N < 2) go to 215 DO 5 I=2,N if (X(I) <= X(I-1)) go to 210 if (X(I) > XUP) go to 6 5 CONTINUE 6 CONTINUE if (N >= 3) go to 9 ! ! SPECIAL N=2 CASE SLOPE = (Y(2)-Y(1))/(X(2)-X(1)) FL = Y(1) + SLOPE*(XLO-X(1)) FR = Y(2) + SLOPE*(XUP-X(2)) ANS = 0.5*(FL+FR)*(XUP-XLO) return 9 CONTINUE if (X(N-2) < XLO) go to 205 if (X(3) > XUP) go to 205 I = 1 10 if (X(I) >= XLO) go to 15 I = I+1 go to 10 15 INLFT = I I = N 20 if (X(I) <= XUP) go to 25 I = I-1 go to 20 25 INRT = I if ((INRT-INLFT) < 2) go to 205 ISTART = INLFT if (INLFT == 1) ISTART = 2 ISTOP = INRT if (INRT == N) ISTOP = N-1 ! R3 = 3.0D0 RP5= 0.5D0 SUM = 0.0 SYL = XLO SYL2= SYL*SYL SYL3= SYL2*SYL ! DO 50 I=ISTART,ISTOP X1 = X(I-1) X2 = X(I) X3 = X(I+1) X12 = X1-X2 X13 = X1-X3 X23 = X2-X3 TERM1 = DBLE(Y(I-1))/(X12*X13) TERM2 =-DBLE(Y(I)) /(X12*X23) TERM3 = DBLE(Y(I+1))/(X13*X23) A = TERM1+TERM2+TERM3 B = -(X2+X3)*TERM1 - (X1+X3)*TERM2 - (X1+X2)*TERM3 C = X2*X3*TERM1 + X1*X3*TERM2 + X1*X2*TERM3 if (I-ISTART) 30,30,35 30 CA = A CB = B CC = C go to 40 35 CA = 0.5*(A+CA) CB = 0.5*(B+CB) CC = 0.5*(C+CC) 40 SYU = X2 SYU2= SYU*SYU SYU3= SYU2*SYU SUM = SUM + CA*(SYU3-SYL3)/R3 + CB*RP5*(SYU2-SYL2) + CC*(SYU-SYL) CA = A CB = B CC = C SYL = SYU SYL2= SYU2 SYL3= SYU3 50 CONTINUE SYU = XUP ANS = SUM + CA*(SYU**3-SYL3)/R3 + CB*RP5*(SYU**2-SYL2) & + CC*(SYU-SYL) 100 RETURN 200 IERR=2 call XERMSG ('SLATEC', 'AVINT', & 'THE UPPER LIMIT OF INTEGRATION WAS NOT GREATER THAN THE ' // & 'LOWER LIMIT.', 4, 1) return 205 IERR=3 call XERMSG ('SLATEC', 'AVINT', & 'THERE WERE LESS THAN THREE FUNCTION VALUES BETWEEN THE ' // & 'LIMITS OF INTEGRATION.', 4, 1) return 210 IERR=4 call XERMSG ('SLATEC', 'AVINT', & 'THE ABSCISSAS WERE NOT STRICTLY INCREASING. MUST HAVE ' // & 'X(I-1) < X(I) FOR ALL I.', 4, 1) return 215 IERR=5 call XERMSG ('SLATEC', 'AVINT', & 'LESS THAN TWO FUNCTION VALUES WERE SUPPLIED.', 4, 1) return end subroutine BAKVEC (NM, N, T, E, M, Z, IERR) ! !! BAKVEC forms the eigenvectors of a certain real non-symmetric tridiagonal ! matrix from a symmetric tridiagonal matrix output from FIGI. ! !***LIBRARY SLATEC (EISPACK) !***CATEGORY D4C4 !***TYPE SINGLE PRECISION (BAKVEC-S) !***KEYWORDS EIGENVECTORS, EISPACK !***AUTHOR Smith, B. T., et al. !***DESCRIPTION ! ! This subroutine forms the eigenvectors of a NONSYMMETRIC ! TRIDIAGONAL matrix by back transforming those of the ! corresponding symmetric matrix determined by FIGI. ! ! On INPUT ! ! NM must be set to the row dimension of the two-dimensional ! array parameters, T and Z, as declared in the calling ! program dimension statement. NM is an INTEGER variable. ! ! N is the order of the matrix T. N is an INTEGER variable. ! N must be less than or equal to NM. ! ! T contains the nonsymmetric matrix. Its subdiagonal is ! stored in the last N-1 positions of the first column, ! its diagonal in the N positions of the second column, ! and its superdiagonal in the first N-1 positions of ! the third column. T(1,1) and T(N,3) are arbitrary. ! T is a two-dimensional REAL array, dimensioned T(NM,3). ! ! E contains the subdiagonal elements of the symmetric ! matrix in its last N-1 positions. E(1) is arbitrary. ! E is a one-dimensional REAL array, dimensioned E(N). ! ! M is the number of eigenvectors to be back transformed. ! M is an INTEGER variable. ! ! Z contains the eigenvectors to be back transformed ! in its first M columns. Z is a two-dimensional REAL ! array, dimensioned Z(NM,M). ! ! On OUTPUT ! ! T is unaltered. ! ! E is destroyed. ! ! Z contains the transformed eigenvectors in its first M columns. ! ! IERR is an INTEGER flag set to ! Zero for normal return, ! 2*N+I if E(I) is zero with T(I,1) or T(I-1,3) non-zero. ! In this case, the symmetric matrix is not similar ! to the original matrix, and the eigenvectors ! cannot be found by this program. ! ! Questions and comments should be directed to B. S. Garbow, ! APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY ! ------------------------------------------------------------------ ! !***REFERENCES B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow, ! Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen- ! system Routines - EISPACK Guide, Springer-Verlag, ! 1976. !***ROUTINES CALLED (NONE) !***REVISION HISTORY (YYMMDD) ! 760101 DATE WRITTEN ! 890831 Modified array declarations. (WRB) ! 890831 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 920501 Reformatted the REFERENCES section. (WRB) !***END PROLOGUE BAKVEC ! INTEGER I,J,M,N,NM,IERR REAL T(NM,3),E(*),Z(NM,*) ! !***FIRST EXECUTABLE STATEMENT BAKVEC IERR = 0 if (M == 0) go to 1001 E(1) = 1.0E0 if (N == 1) go to 1001 ! DO 100 I = 2, N if (E(I) /= 0.0E0) go to 80 if (T(I,1) /= 0.0E0 .OR. T(I-1,3) /= 0.0E0) go to 1000 E(I) = 1.0E0 go to 100 80 E(I) = E(I-1) * E(I) / T(I-1,3) 100 CONTINUE ! DO 120 J = 1, M ! DO 120 I = 2, N Z(I,J) = Z(I,J) * E(I) 120 CONTINUE ! go to 1001 ! .......... SET ERROR -- EIGENVECTORS CANNOT BE ! FOUND BY THIS PROGRAM .......... 1000 IERR = 2 * N + I 1001 RETURN end subroutine BALANC (NM, N, A, LOW, IGH, SCALE) ! !! BALANC balances a real general matrix and isolates eigenvalues when possible. ! !***LIBRARY SLATEC (EISPACK) !***CATEGORY D4C1A !***TYPE SINGLE PRECISION (BALANC-S, CBAL-C) !***KEYWORDS EIGENVECTORS, EISPACK !***AUTHOR Smith, B. T., et al. !***DESCRIPTION ! ! This subroutine is a translation of the ALGOL procedure BALANCE, ! NUM. MATH. 13, 293-304(1969) by Parlett and Reinsch. ! HANDBOOK FOR AUTO. COMP., Vol.II-LINEAR ALGEBRA, 315-326(1971). ! ! This subroutine balances a REAL matrix and isolates ! eigenvalues whenever possible. ! ! On INPUT ! ! NM must be set to the row dimension of the two-dimensional ! array parameter, A, as declared in the calling program ! dimension statement. NM is an INTEGER variable. ! ! N is the order of the matrix A. N is an INTEGER variable. ! N must be less than or equal to NM. ! ! A contains the input matrix to be balanced. A is a ! two-dimensional REAL array, dimensioned A(NM,N). ! ! On OUTPUT ! ! A contains the balanced matrix. ! ! LOW and IGH are two INTEGER variables such that A(I,J) ! is equal to zero if ! (1) I is greater than J and ! (2) J=1,...,LOW-1 or I=IGH+1,...,N. ! ! SCALE contains information determining the permutations and ! scaling factors used. SCALE is a one-dimensional REAL array, ! dimensioned SCALE(N). ! ! Suppose that the principal submatrix in rows LOW through IGH ! has been balanced, that P(J) denotes the index interchanged ! with J during the permutation step, and that the elements ! of the diagonal matrix used are denoted by D(I,J). Then ! SCALE(J) = P(J), for J = 1,...,LOW-1 ! = D(J,J), J = LOW,...,IGH ! = P(J) J = IGH+1,...,N. ! The order in which the interchanges are made is N to IGH+1, ! then 1 TO LOW-1. ! ! Note that 1 is returned for IGH if IGH is zero formally. ! ! The ALGOL procedure EXC contained in BALANCE appears in ! BALANC in line. (Note that the ALGOL roles of identifiers ! K,L have been reversed.) ! ! Questions and comments should be directed to B. S. Garbow, ! Applied Mathematics Division, ARGONNE NATIONAL LABORATORY ! ------------------------------------------------------------------ ! !***REFERENCES B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow, ! Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen- ! system Routines - EISPACK Guide, Springer-Verlag, ! 1976. !***ROUTINES CALLED (NONE) !***REVISION HISTORY (YYMMDD) ! 760101 DATE WRITTEN ! 890831 Modified array declarations. (WRB) ! 890831 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 920501 Reformatted the REFERENCES section. (WRB) !***END PROLOGUE BALANC ! INTEGER I,J,K,L,M,N,JJ,NM,IGH,LOW,IEXC REAL A(NM,*),SCALE(*) REAL C,F,G,R,S,B2,RADIX LOGICAL NOCONV ! !***FIRST EXECUTABLE STATEMENT BALANC RADIX = 16 ! B2 = RADIX * RADIX K = 1 L = N go to 100 ! .......... IN-LINE PROCEDURE FOR ROW AND ! COLUMN EXCHANGE .......... 20 SCALE(M) = J if (J == M) go to 50 ! DO 30 I = 1, L F = A(I,J) A(I,J) = A(I,M) A(I,M) = F 30 CONTINUE ! DO 40 I = K, N F = A(J,I) A(J,I) = A(M,I) A(M,I) = F 40 CONTINUE ! 50 go to (80,130), IEXC ! .......... SEARCH FOR ROWS ISOLATING AN EIGENVALUE ! AND PUSH THEM DOWN .......... 80 if (L == 1) go to 280 L = L - 1 ! .......... FOR J=L STEP -1 UNTIL 1 DO -- .......... 100 DO 120 JJ = 1, L J = L + 1 - JJ ! DO 110 I = 1, L if (I == J) go to 110 if (A(J,I) /= 0.0E0) go to 120 110 CONTINUE ! M = L IEXC = 1 go to 20 120 CONTINUE ! go to 140 ! .......... SEARCH FOR COLUMNS ISOLATING AN EIGENVALUE ! AND PUSH THEM LEFT .......... 130 K = K + 1 ! 140 DO 170 J = K, L ! DO 150 I = K, L if (I == J) go to 150 if (A(I,J) /= 0.0E0) go to 170 150 CONTINUE ! M = K IEXC = 2 go to 20 170 CONTINUE ! .......... NOW BALANCE THE SUBMATRIX IN ROWS K TO L .......... DO 180 I = K, L 180 SCALE(I) = 1.0E0 ! .......... ITERATIVE LOOP FOR NORM REDUCTION .......... 190 NOCONV = .FALSE. ! DO 270 I = K, L C = 0.0E0 R = 0.0E0 ! DO 200 J = K, L if (J == I) go to 200 C = C + ABS(A(J,I)) R = R + ABS(A(I,J)) 200 CONTINUE ! .......... GUARD AGAINST ZERO C OR R DUE TO UNDERFLOW .......... if (C == 0.0E0 .OR. R == 0.0E0) go to 270 G = R / RADIX F = 1.0E0 S = C + R 210 if (C >= G) go to 220 F = F * RADIX C = C * B2 go to 210 220 G = R * RADIX 230 if (C < G) go to 240 F = F / RADIX C = C / B2 go to 230 ! .......... NOW BALANCE .......... 240 if ((C + R) / F >= 0.95E0 * S) go to 270 G = 1.0E0 / F SCALE(I) = SCALE(I) * F NOCONV = .TRUE. ! DO 250 J = K, N 250 A(I,J) = A(I,J) * G ! DO 260 J = 1, L 260 A(J,I) = A(J,I) * F ! 270 CONTINUE ! if (NOCONV) go to 190 ! 280 LOW = K IGH = L return end subroutine BALBAK (NM, N, LOW, IGH, SCALE, M, Z) ! !! BALBAK forms the eigenvectors of a real general matrix from the ... ! eigenvectors of matrix output from BALANC. ! !***LIBRARY SLATEC (EISPACK) !***CATEGORY D4C4 !***TYPE SINGLE PRECISION (BALBAK-S, CBABK2-C) !***KEYWORDS EIGENVECTORS, EISPACK !***AUTHOR Smith, B. T., et al. !***DESCRIPTION ! ! This subroutine is a translation of the ALGOL procedure BALBAK, ! NUM. MATH. 13, 293-304(1969) by Parlett and Reinsch. ! HANDBOOK FOR AUTO. COMP., Vol.II-LINEAR ALGEBRA, 315-326(1971). ! ! This subroutine forms the eigenvectors of a REAL GENERAL ! matrix by back transforming those of the corresponding ! balanced matrix determined by BALANC. ! ! On INPUT ! ! NM must be set to the row dimension of the two-dimensional ! array parameter, Z, as declared in the calling program ! dimension statement. NM is an INTEGER variable. ! ! N is the number of components of the vectors in matrix Z. ! N is an INTEGER variable. N must be less than or equal ! to NM. ! ! LOW and IGH are INTEGER variables determined by BALANC. ! ! SCALE contains information determining the permutations and ! scaling factors used by BALANC. SCALE is a one-dimensional ! REAL array, dimensioned SCALE(N). ! ! M is the number of columns of Z to be back transformed. ! M is an INTEGER variable. ! ! Z contains the real and imaginary parts of the eigen- ! vectors to be back transformed in its first M columns. ! Z is a two-dimensional REAL array, dimensioned Z(NM,M). ! ! On OUTPUT ! ! Z contains the real and imaginary parts of the ! transformed eigenvectors in its first M columns. ! ! Questions and comments should be directed to B. S. Garbow, ! Applied Mathematics Division, ARGONNE NATIONAL LABORATORY ! ------------------------------------------------------------------ ! !***REFERENCES B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow, ! Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen- ! system Routines - EISPACK Guide, Springer-Verlag, ! 1976. !***ROUTINES CALLED (NONE) !***REVISION HISTORY (YYMMDD) ! 760101 DATE WRITTEN ! 890831 Modified array declarations. (WRB) ! 890831 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 920501 Reformatted the REFERENCES section. (WRB) !***END PROLOGUE BALBAK ! INTEGER I,J,K,M,N,II,NM,IGH,LOW REAL SCALE(*),Z(NM,*) REAL S ! !***FIRST EXECUTABLE STATEMENT BALBAK if (M == 0) go to 200 if (IGH == LOW) go to 120 ! DO 110 I = LOW, IGH S = SCALE(I) ! .......... LEFT HAND EIGENVECTORS ARE BACK TRANSFORMED ! if THE FOREGOING STATEMENT IS REPLACED BY ! S=1.0E0/SCALE(I). .......... DO 100 J = 1, M 100 Z(I,J) = Z(I,J) * S ! 110 CONTINUE ! ......... FOR I=LOW-1 STEP -1 UNTIL 1, ! IGH+1 STEP 1 UNTIL N DO -- .......... 120 DO 140 II = 1, N I = II if (I >= LOW .AND. I <= IGH) go to 140 if (I < LOW) I = LOW - II K = SCALE(I) if (K == I) go to 140 ! DO 130 J = 1, M S = Z(I,J) Z(I,J) = Z(K,J) Z(K,J) = S 130 CONTINUE ! 140 CONTINUE ! 200 RETURN end subroutine BANDR (NM, N, MB, A, D, E, E2, MATZ, Z) ! !! BANDR reduces a real symmetric band matrix to symmetric tridiagonal ... ! matrix and, optionally, accumulates orthogonal similarity transformations. ! !***LIBRARY SLATEC (EISPACK) !***CATEGORY D4C1B1 !***TYPE SINGLE PRECISION (BANDR-S) !***KEYWORDS EIGENVALUES, EIGENVECTORS, EISPACK !***AUTHOR Smith, B. T., et al. !***DESCRIPTION ! ! This subroutine is a translation of the ALGOL procedure BANDRD, ! NUM. MATH. 12, 231-241(1968) by Schwarz. ! HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 273-283(1971). ! ! This subroutine reduces a REAL SYMMETRIC BAND matrix ! to a symmetric tridiagonal matrix using and optionally ! accumulating orthogonal similarity transformations. ! ! On INPUT ! ! NM must be set to the row dimension of the two-dimensional ! array parameters, A and Z, as declared in the calling ! program dimension statement. NM is an INTEGER variable. ! ! N is the order of the matrix A. N is an INTEGER variable. ! N must be less than or equal to NM. ! ! MB is the (half) band width of the matrix, defined as the ! number of adjacent diagonals, including the principal ! diagonal, required to specify the non-zero portion of the ! lower triangle of the matrix. MB is less than or equal ! to N. MB is an INTEGER variable. ! ! A contains the lower triangle of the real symmetric band ! matrix. Its lowest subdiagonal is stored in the last ! N+1-MB positions of the first column, its next subdiagonal ! in the last N+2-MB positions of the second column, further ! subdiagonals similarly, and finally its principal diagonal ! in the N positions of the last column. Contents of storage ! locations not part of the matrix are arbitrary. A is a ! two-dimensional REAL array, dimensioned A(NM,MB). ! ! MATZ should be set to .TRUE. if the transformation matrix is ! to be accumulated, and to .FALSE. otherwise. MATZ is a ! LOGICAL variable. ! ! On OUTPUT ! ! A has been destroyed, except for its last two columns which ! contain a copy of the tridiagonal matrix. ! ! D contains the diagonal elements of the tridiagonal matrix. ! D is a one-dimensional REAL array, dimensioned D(N). ! ! E contains the subdiagonal elements of the tridiagonal ! matrix in its last N-1 positions. E(1) is set to zero. ! E is a one-dimensional REAL array, dimensioned E(N). ! ! E2 contains the squares of the corresponding elements of E. ! E2 may coincide with E if the squares are not needed. ! E2 is a one-dimensional REAL array, dimensioned E2(N). ! ! Z contains the orthogonal transformation matrix produced in ! the reduction if MATZ has been set to .TRUE. Otherwise, Z ! is not referenced. Z is a two-dimensional REAL array, ! dimensioned Z(NM,N). ! ! Questions and comments should be directed to B. S. Garbow, ! Applied Mathematics Division, ARGONNE NATIONAL LABORATORY ! ------------------------------------------------------------------ ! !***REFERENCES B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow, ! Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen- ! system Routines - EISPACK Guide, Springer-Verlag, ! 1976. !***ROUTINES CALLED (NONE) !***REVISION HISTORY (YYMMDD) ! 760101 DATE WRITTEN ! 890531 Changed all specific intrinsics to generic. (WRB) ! 890831 Modified array declarations. (WRB) ! 890831 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 920501 Reformatted the REFERENCES section. (WRB) !***END PROLOGUE BANDR ! INTEGER J,K,L,N,R,I1,I2,J1,J2,KR,MB,MR,M1,NM,N2,R1,UGL,MAXL,MAXR REAL A(NM,*),D(*),E(*),E2(*),Z(NM,*) REAL G,U,B1,B2,C2,F1,F2,S2,DMIN,DMINRT LOGICAL MATZ ! !***FIRST EXECUTABLE STATEMENT BANDR DMIN = 2.0E0**(-64) DMINRT = 2.0E0**(-32) ! .......... INITIALIZE DIAGONAL SCALING MATRIX .......... DO 30 J = 1, N 30 D(J) = 1.0E0 ! if (.NOT. MATZ) go to 60 ! DO 50 J = 1, N ! DO 40 K = 1, N 40 Z(J,K) = 0.0E0 ! Z(J,J) = 1.0E0 50 CONTINUE ! 60 M1 = MB - 1 if (M1 - 1) 900, 800, 70 70 N2 = N - 2 ! DO 700 K = 1, N2 MAXR = MIN(M1,N-K) ! .......... FOR R=MAXR STEP -1 UNTIL 2 DO -- .......... DO 600 R1 = 2, MAXR R = MAXR + 2 - R1 KR = K + R MR = MB - R G = A(KR,MR) A(KR-1,1) = A(KR-1,MR+1) UGL = K ! DO 500 J = KR, N, M1 J1 = J - 1 J2 = J1 - 1 if (G == 0.0E0) go to 600 B1 = A(J1,1) / G B2 = B1 * D(J1) / D(J) S2 = 1.0E0 / (1.0E0 + B1 * B2) if (S2 >= 0.5E0 ) go to 450 B1 = G / A(J1,1) B2 = B1 * D(J) / D(J1) C2 = 1.0E0 - S2 D(J1) = C2 * D(J1) D(J) = C2 * D(J) F1 = 2.0E0 * A(J,M1) F2 = B1 * A(J1,MB) A(J,M1) = -B2 * (B1 * A(J,M1) - A(J,MB)) - F2 + A(J,M1) A(J1,MB) = B2 * (B2 * A(J,MB) + F1) + A(J1,MB) A(J,MB) = B1 * (F2 - F1) + A(J,MB) ! DO 200 L = UGL, J2 I2 = MB - J + L U = A(J1,I2+1) + B2 * A(J,I2) A(J,I2) = -B1 * A(J1,I2+1) + A(J,I2) A(J1,I2+1) = U 200 CONTINUE ! UGL = J A(J1,1) = A(J1,1) + B2 * G if (J == N) go to 350 MAXL = MIN(M1,N-J1) ! DO 300 L = 2, MAXL I1 = J1 + L I2 = MB - L U = A(I1,I2) + B2 * A(I1,I2+1) A(I1,I2+1) = -B1 * A(I1,I2) + A(I1,I2+1) A(I1,I2) = U 300 CONTINUE ! I1 = J + M1 if (I1 > N) go to 350 G = B2 * A(I1,1) 350 if (.NOT. MATZ) go to 500 ! DO 400 L = 1, N U = Z(L,J1) + B2 * Z(L,J) Z(L,J) = -B1 * Z(L,J1) + Z(L,J) Z(L,J1) = U 400 CONTINUE ! go to 500 ! 450 U = D(J1) D(J1) = S2 * D(J) D(J) = S2 * U F1 = 2.0E0 * A(J,M1) F2 = B1 * A(J,MB) U = B1 * (F2 - F1) + A(J1,MB) A(J,M1) = B2 * (B1 * A(J,M1) - A(J1,MB)) + F2 - A(J,M1) A(J1,MB) = B2 * (B2 * A(J1,MB) + F1) + A(J,MB) A(J,MB) = U ! DO 460 L = UGL, J2 I2 = MB - J + L U = B2 * A(J1,I2+1) + A(J,I2) A(J,I2) = -A(J1,I2+1) + B1 * A(J,I2) A(J1,I2+1) = U 460 CONTINUE ! UGL = J A(J1,1) = B2 * A(J1,1) + G if (J == N) go to 480 MAXL = MIN(M1,N-J1) ! DO 470 L = 2, MAXL I1 = J1 + L I2 = MB - L U = B2 * A(I1,I2) + A(I1,I2+1) A(I1,I2+1) = -A(I1,I2) + B1 * A(I1,I2+1) A(I1,I2) = U 470 CONTINUE ! I1 = J + M1 if (I1 > N) go to 480 G = A(I1,1) A(I1,1) = B1 * A(I1,1) 480 if (.NOT. MATZ) go to 500 ! DO 490 L = 1, N U = B2 * Z(L,J1) + Z(L,J) Z(L,J) = -Z(L,J1) + B1 * Z(L,J) Z(L,J1) = U 490 CONTINUE ! 500 CONTINUE ! 600 CONTINUE ! if (MOD(K,64) /= 0) go to 700 ! .......... RESCALE TO AVOID UNDERFLOW OR OVERFLOW .......... DO 650 J = K, N if (D(J) >= DMIN) go to 650 MAXL = MAX(1,MB+1-J) ! DO 610 L = MAXL, M1 610 A(J,L) = DMINRT * A(J,L) ! if (J == N) go to 630 MAXL = MIN(M1,N-J) ! DO 620 L = 1, MAXL I1 = J + L I2 = MB - L A(I1,I2) = DMINRT * A(I1,I2) 620 CONTINUE ! 630 if (.NOT. MATZ) go to 645 ! DO 640 L = 1, N 640 Z(L,J) = DMINRT * Z(L,J) ! 645 A(J,MB) = DMIN * A(J,MB) D(J) = D(J) / DMIN 650 CONTINUE ! 700 CONTINUE ! .......... FORM SQUARE ROOT OF SCALING MATRIX .......... 800 DO 810 J = 2, N 810 E(J) = SQRT(D(J)) ! if (.NOT. MATZ) go to 840 ! DO 830 J = 1, N ! DO 820 K = 2, N 820 Z(J,K) = E(K) * Z(J,K) ! 830 CONTINUE ! 840 U = 1.0E0 ! DO 850 J = 2, N A(J,M1) = U * E(J) * A(J,M1) U = E(J) E2(J) = A(J,M1) ** 2 A(J,MB) = D(J) * A(J,MB) D(J) = A(J,MB) E(J) = A(J,M1) 850 CONTINUE ! D(1) = A(1,MB) E(1) = 0.0E0 E2(1) = 0.0E0 go to 1001 ! 900 DO 950 J = 1, N D(J) = A(J,MB) E(J) = 0.0E0 E2(J) = 0.0E0 950 CONTINUE ! 1001 RETURN end subroutine BANDV (NM, N, MBW, A, E21, M, W, Z, IERR, NV, RV, RV6) ! !! BANDV forms the eigenvectors of a real symmetric band matrix ... ! associated with a set of ordered approximate eigenvalues ! by inverse iteration. ! !***LIBRARY SLATEC (EISPACK) !***CATEGORY D4C3 !***TYPE SINGLE PRECISION (BANDV-S) !***KEYWORDS EIGENVECTORS, EISPACK !***AUTHOR Smith, B. T., et al. !***DESCRIPTION ! ! This subroutine finds those eigenvectors of a REAL SYMMETRIC ! BAND matrix corresponding to specified eigenvalues, using inverse ! iteration. The subroutine may also be used to solve systems ! of linear equations with a symmetric or non-symmetric band ! coefficient matrix. ! ! On INPUT ! ! NM must be set to the row dimension of the two-dimensional ! array parameters, A and Z, as declared in the calling ! program dimension statement. NM is an INTEGER variable. ! ! N is the order of the matrix A. N is an INTEGER variable. ! N must be less than or equal to NM. ! ! MBW is the number of columns of the array A used to store the ! band matrix. If the matrix is symmetric, MBW is its (half) ! band width, denoted MB and defined as the number of adjacent ! diagonals, including the principal diagonal, required to ! specify the non-zero portion of the lower triangle of the ! matrix. If the subroutine is being used to solve systems ! of linear equations and the coefficient matrix is not ! symmetric, it must however have the same number of adjacent ! diagonals above the main diagonal as below, and in this ! case, MBW=2*MB-1. MBW is an INTEGER variable. MB must not ! be greater than N. ! ! A contains the lower triangle of the symmetric band input ! matrix stored as an N by MB array. Its lowest subdiagonal ! is stored in the last N+1-MB positions of the first column, ! its next subdiagonal in the last N+2-MB positions of the ! second column, further subdiagonals similarly, and finally ! its principal diagonal in the N positions of column MB. ! If the subroutine is being used to solve systems of linear ! equations and the coefficient matrix is not symmetric, A is ! N by 2*MB-1 instead with lower triangle as above and with ! its first superdiagonal stored in the first N-1 positions of ! column MB+1, its second superdiagonal in the first N-2 ! positions of column MB+2, further superdiagonals similarly, ! and finally its highest superdiagonal in the first N+1-MB ! positions of the last column. Contents of storage locations ! not part of the matrix are arbitrary. A is a two-dimensional ! REAL array, dimensioned A(NM,MBW). ! ! E21 specifies the ordering of the eigenvalues and contains ! 0.0E0 if the eigenvalues are in ascending order, or ! 2.0E0 if the eigenvalues are in descending order. ! If the subroutine is being used to solve systems of linear ! equations, E21 should be set to 1.0E0 if the coefficient ! matrix is symmetric and to -1.0E0 if not. E21 is a REAL ! variable. ! ! M is the number of specified eigenvalues or the number of ! systems of linear equations. M is an INTEGER variable. ! ! W contains the M eigenvalues in ascending or descending order. ! If the subroutine is being used to solve systems of linear ! equations (A-W(J)*I)*X(J)=B(J), where I is the identity ! matrix, W(J) should be set accordingly, for J=1,2,...,M. ! W is a one-dimensional REAL array, dimensioned W(M). ! ! Z contains the constant matrix columns (B(J),J=1,2,...,M), if ! the subroutine is used to solve systems of linear equations. ! Z is a two-dimensional REAL array, dimensioned Z(NM,M). ! ! NV must be set to the dimension of the array parameter RV ! as declared in the calling program dimension statement. ! NV is an INTEGER variable. ! ! On OUTPUT ! ! A and W are unaltered. ! ! Z contains the associated set of orthogonal eigenvectors. ! Any vector which fails to converge is set to zero. If the ! subroutine is used to solve systems of linear equations, ! Z contains the solution matrix columns (X(J),J=1,2,...,M). ! ! IERR is an INTEGER flag set to ! Zero for normal return, ! -J if the eigenvector corresponding to the J-th ! eigenvalue fails to converge, or if the J-th ! system of linear equations is nearly singular. ! ! RV and RV6 are temporary storage arrays. If the subroutine ! is being used to solve systems of linear equations, the ! determinant (up to sign) of A-W(M)*I is available, upon ! return, as the product of the first N elements of RV. ! RV and RV6 are one-dimensional REAL arrays. Note that RV ! is dimensioned RV(NV), where NV must be at least N*(2*MB-1). ! RV6 is dimensioned RV6(N). ! ! Questions and comments should be directed to B. S. Garbow, ! Applied Mathematics Division, ARGONNE NATIONAL LABORATORY ! ------------------------------------------------------------------ ! !***REFERENCES B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow, ! Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen- ! system Routines - EISPACK Guide, Springer-Verlag, ! 1976. !***ROUTINES CALLED (NONE) !***REVISION HISTORY (YYMMDD) ! 760101 DATE WRITTEN ! 890531 Changed all specific intrinsics to generic. (WRB) ! 890831 Modified array declarations. (WRB) ! 890831 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 920501 Reformatted the REFERENCES section. (WRB) !***END PROLOGUE BANDV ! INTEGER I,J,K,M,N,R,II,IJ,JJ,KJ,MB,M1,NM,NV,IJ1,ITS,KJ1,MBW,M21 INTEGER IERR,MAXJ,MAXK,GROUP REAL A(NM,*),W(*),Z(NM,*),RV(*),RV6(*) REAL U,V,UK,XU,X0,X1,E21,EPS2,EPS3,EPS4,NORM,ORDER,S ! !***FIRST EXECUTABLE STATEMENT BANDV IERR = 0 if (M == 0) go to 1001 MB = MBW if (E21 < 0.0E0) MB = (MBW + 1) / 2 M1 = MB - 1 M21 = M1 + MB ORDER = 1.0E0 - ABS(E21) ! .......... FIND VECTORS BY INVERSE ITERATION .......... DO 920 R = 1, M ITS = 1 X1 = W(R) if (R /= 1) go to 100 ! .......... COMPUTE NORM OF MATRIX .......... NORM = 0.0E0 ! DO 60 J = 1, MB JJ = MB + 1 - J KJ = JJ + M1 IJ = 1 S = 0.0E0 ! DO 40 I = JJ, N S = S + ABS(A(I,J)) if (E21 >= 0.0E0) go to 40 S = S + ABS(A(IJ,KJ)) IJ = IJ + 1 40 CONTINUE ! NORM = MAX(NORM,S) 60 CONTINUE ! if (E21 < 0.0E0) NORM = 0.5E0 * NORM ! .......... EPS2 IS THE CRITERION FOR GROUPING, ! EPS3 REPLACES ZERO PIVOTS AND EQUAL ! ROOTS ARE MODIFIED BY EPS3, ! EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW .......... if (NORM == 0.0E0) NORM = 1.0E0 EPS2 = 1.0E-3 * NORM * ABS(ORDER) EPS3 = NORM 70 EPS3 = 0.5E0*EPS3 if (NORM + EPS3 > NORM) go to 70 UK = SQRT(REAL(N)) EPS3 = UK * EPS3 EPS4 = UK * EPS3 80 GROUP = 0 go to 120 ! .......... LOOK FOR CLOSE OR COINCIDENT ROOTS .......... 100 if (ABS(X1-X0) >= EPS2) go to 80 GROUP = GROUP + 1 if (ORDER * (X1 - X0) <= 0.0E0) X1 = X0 + ORDER * EPS3 ! .......... EXPAND MATRIX, SUBTRACT EIGENVALUE, ! AND INITIALIZE VECTOR .......... 120 DO 200 I = 1, N IJ = I + MIN(0,I-M1) * N KJ = IJ + MB * N IJ1 = KJ + M1 * N if (M1 == 0) go to 180 ! DO 150 J = 1, M1 if (IJ > M1) go to 125 if (IJ > 0) go to 130 RV(IJ1) = 0.0E0 IJ1 = IJ1 + N go to 130 125 RV(IJ) = A(I,J) 130 IJ = IJ + N II = I + J if (II > N) go to 150 JJ = MB - J if (E21 >= 0.0E0) go to 140 II = I JJ = MB + J 140 RV(KJ) = A(II,JJ) KJ = KJ + N 150 CONTINUE ! 180 RV(IJ) = A(I,MB) - X1 RV6(I) = EPS4 if (ORDER == 0.0E0) RV6(I) = Z(I,R) 200 CONTINUE ! if (M1 == 0) go to 600 ! .......... ELIMINATION WITH INTERCHANGES .......... DO 580 I = 1, N II = I + 1 MAXK = MIN(I+M1-1,N) MAXJ = MIN(N-I,M21-2) * N ! DO 360 K = I, MAXK KJ1 = K J = KJ1 + N JJ = J + MAXJ ! DO 340 KJ = J, JJ, N RV(KJ1) = RV(KJ) KJ1 = KJ 340 CONTINUE ! RV(KJ1) = 0.0E0 360 CONTINUE ! if (I == N) go to 580 U = 0.0E0 MAXK = MIN(I+M1,N) MAXJ = MIN(N-II,M21-2) * N ! DO 450 J = I, MAXK if (ABS(RV(J)) < ABS(U)) go to 450 U = RV(J) K = J 450 CONTINUE ! J = I + N JJ = J + MAXJ if (K == I) go to 520 KJ = K ! DO 500 IJ = I, JJ, N V = RV(IJ) RV(IJ) = RV(KJ) RV(KJ) = V KJ = KJ + N 500 CONTINUE ! if (ORDER /= 0.0E0) go to 520 V = RV6(I) RV6(I) = RV6(K) RV6(K) = V 520 if (U == 0.0E0) go to 580 ! DO 560 K = II, MAXK V = RV(K) / U KJ = K ! DO 540 IJ = J, JJ, N KJ = KJ + N RV(KJ) = RV(KJ) - V * RV(IJ) 540 CONTINUE ! if (ORDER == 0.0E0) RV6(K) = RV6(K) - V * RV6(I) 560 CONTINUE ! 580 CONTINUE ! .......... BACK SUBSTITUTION ! FOR I=N STEP -1 UNTIL 1 DO -- .......... 600 DO 630 II = 1, N I = N + 1 - II MAXJ = MIN(II,M21) if (MAXJ == 1) go to 620 IJ1 = I J = IJ1 + N JJ = J + (MAXJ - 2) * N ! DO 610 IJ = J, JJ, N IJ1 = IJ1 + 1 RV6(I) = RV6(I) - RV(IJ) * RV6(IJ1) 610 CONTINUE ! 620 V = RV(I) if (ABS(V) >= EPS3) go to 625 ! .......... SET ERROR -- NEARLY SINGULAR LINEAR SYSTEM .......... if (ORDER == 0.0E0) IERR = -R V = SIGN(EPS3,V) 625 RV6(I) = RV6(I) / V 630 CONTINUE ! XU = 1.0E0 if (ORDER == 0.0E0) go to 870 ! .......... ORTHOGONALIZE WITH RESPECT TO PREVIOUS ! MEMBERS OF GROUP .......... if (GROUP == 0) go to 700 ! DO 680 JJ = 1, GROUP J = R - GROUP - 1 + JJ XU = 0.0E0 ! DO 640 I = 1, N 640 XU = XU + RV6(I) * Z(I,J) ! DO 660 I = 1, N 660 RV6(I) = RV6(I) - XU * Z(I,J) ! 680 CONTINUE ! 700 NORM = 0.0E0 ! DO 720 I = 1, N 720 NORM = NORM + ABS(RV6(I)) ! if (NORM >= 0.1E0) go to 840 ! .......... IN-LINE PROCEDURE FOR CHOOSING ! A NEW STARTING VECTOR .......... if (ITS >= N) go to 830 ITS = ITS + 1 XU = EPS4 / (UK + 1.0E0) RV6(1) = EPS4 ! DO 760 I = 2, N 760 RV6(I) = XU ! RV6(ITS) = RV6(ITS) - EPS4 * UK go to 600 ! .......... SET ERROR -- NON-CONVERGED EIGENVECTOR .......... 830 IERR = -R XU = 0.0E0 go to 870 ! .......... NORMALIZE SO THAT SUM OF SQUARES IS ! 1 AND EXPAND TO FULL ORDER .......... 840 U = 0.0E0 ! DO 860 I = 1, N 860 U = U + RV6(I)**2 ! XU = 1.0E0 / SQRT(U) ! 870 DO 900 I = 1, N 900 Z(I,R) = RV6(I) * XU ! X0 = X1 920 CONTINUE ! 1001 RETURN end function BCRH (XLL, XRR, IZ, C, A, BH, F, SGN) ! !! BCRH is subsidiary to CBLKTR ! !***LIBRARY SLATEC !***TYPE SINGLE PRECISION (BCRH-S, BSRH-S) !***AUTHOR (UNKNOWN) !***SEE ALSO CBLKTR !***ROUTINES CALLED (NONE) !***COMMON BLOCKS CCBLK !***REVISION HISTORY (YYMMDD) ! 801001 DATE WRITTEN ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900402 Added TYPE section. (WRB) !***END PROLOGUE BCRH DIMENSION A(*) ,C(*) ,BH(*) COMMON /CCBLK/ NPP ,K ,EPS ,CNV , & NM ,NCMPLX ,IK !***FIRST EXECUTABLE STATEMENT BCRH XL = XLL XR = XRR DX = .5*ABS(XR-XL) 101 X = .5*(XL+XR) if (SGN*F(X,IZ,C,A,BH)) 103,105,102 102 XR = X go to 104 103 XL = X 104 DX = .5*DX if (DX-CNV) 105,105,101 105 BCRH = .5*(XL+XR) return end subroutine BDIFF (L, V) ! !! BDIFF is subsidiary to BSKIN ! !***LIBRARY SLATEC !***TYPE SINGLE PRECISION (BDIFF-S, DBDIFF-D) !***AUTHOR Amos, D. E., (SNLA) !***DESCRIPTION ! ! BDIFF computes the sum of B(L,K)*V(K)*(-1)**K where B(L,K) ! are the binomial coefficients. Truncated sums are computed by ! setting last part of the V vector to zero. On return, the binomial ! sum is in V(L). ! !***SEE ALSO BSKIN !***ROUTINES CALLED (NONE) !***REVISION HISTORY (YYMMDD) ! 820601 DATE WRITTEN ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900328 Added TYPE section. (WRB) !***END PROLOGUE BDIFF INTEGER I, J, K, L REAL V DIMENSION V(*) !***FIRST EXECUTABLE STATEMENT BDIFF if (L == 1) RETURN DO 20 J=2,L K = L DO 10 I=J,L V(K) = V(K-1) - V(K) K = K - 1 10 CONTINUE 20 CONTINUE return end subroutine BESI (X, ALPHA, KODE, N, Y, NZ) ! !! BESI computes an N member sequence of I Bessel functions ! I/SUB(ALPHA+K-1)/(X), K=1,...,N or scaled Bessel functions ! EXP(-X)*I/SUB(ALPHA+K-1)/(X), K=1,...,N for non-negative ! ALPHA and X. ! !***LIBRARY SLATEC !***CATEGORY C10B3 !***TYPE SINGLE PRECISION (BESI-S, DBESI-D) !***KEYWORDS I BESSEL FUNCTION, SPECIAL FUNCTIONS !***AUTHOR Amos, D. E., (SNLA) ! Daniel, S. L., (SNLA) !***DESCRIPTION ! ! Abstract ! BESI computes an N member sequence of I Bessel functions ! I/sub(ALPHA+K-1)/(X), K=1,...,N or scaled Bessel functions ! EXP(-X)*I/sub(ALPHA+K-1)/(X), K=1,...,N for non-negative ALPHA ! and X. A combination of the power series, the asymptotic ! expansion for X to infinity, and the uniform asymptotic ! expansion for NU to infinity are applied over subdivisions of ! the (NU,X) plane. For values not covered by one of these ! formulae, the order is incremented by an integer so that one ! of these formulae apply. Backward recursion is used to reduce ! orders by integer values. The asymptotic expansion for X to ! infinity is used only when the entire sequence (specifically ! the last member) lies within the region covered by the ! expansion. Leading terms of these expansions are used to test ! for over or underflow where appropriate. If a sequence is ! requested and the last member would underflow, the result is ! set to zero and the next lower order tried, etc., until a ! member comes on scale or all are set to zero. An overflow ! cannot occur with scaling. ! ! Description of Arguments ! ! Input ! X - X >= 0.0E0 ! ALPHA - order of first member of the sequence, ! ALPHA >= 0.0E0 ! KODE - a parameter to indicate the scaling option ! KODE=1 returns ! Y(K)= I/sub(ALPHA+K-1)/(X), ! K=1,...,N ! KODE=2 returns ! Y(K)=EXP(-X)*I/sub(ALPHA+K-1)/(X), ! K=1,...,N ! N - number of members in the sequence, N >= 1 ! ! Output ! Y - a vector whose first N components contain ! values for I/sub(ALPHA+K-1)/(X) or scaled ! values for EXP(-X)*I/sub(ALPHA+K-1)/(X), ! K=1,...,N depending on KODE ! NZ - number of components of Y set to zero due to ! underflow, ! NZ=0 , normal return, computation completed ! NZ /= 0, last NZ components of Y set to zero, ! Y(K)=0.0E0, K=N-NZ+1,...,N. ! ! Error Conditions ! Improper input arguments - a fatal error ! Overflow with KODE=1 - a fatal error ! Underflow - a non-fatal error (NZ /= 0) ! !***REFERENCES D. E. Amos, S. L. Daniel and M. K. Weston, CDC 6600 ! subroutines IBESS and JBESS for Bessel functions ! I(NU,X) and J(NU,X), X >= 0, NU >= 0, ACM ! Transactions on Mathematical Software 3, (1977), ! pp. 76-92. ! F. W. J. Olver, Tables of Bessel Functions of Moderate ! or Large Orders, NPL Mathematical Tables 6, Her ! Majesty's Stationery Office, London, 1962. !***ROUTINES CALLED ALNGAM, ASYIK, I1MACH, R1MACH, XERMSG !***REVISION HISTORY (YYMMDD) ! 750101 DATE WRITTEN ! 890531 Changed all specific intrinsics to generic. (WRB) ! 890531 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) ! 900326 Removed duplicate information from DESCRIPTION section. ! (WRB) ! 920501 Reformatted the REFERENCES section. (WRB) !***END PROLOGUE BESI ! INTEGER I, IALP, IN, INLIM, IS, I1, K, KK, KM, KODE, KT, & N, NN, NS, NZ INTEGER I1MACH REAL AIN, AK, AKM, ALPHA, ANS, AP, ARG, ATOL, TOLLN, DFN, & DTM, DX, EARG, ELIM, ETX, FLGIK,FN, FNF, FNI,FNP1,FNU,GLN,RA, & RTTPI, S, SX, SXO2, S1, S2, T, TA, TB, TEMP, TFN, TM, TOL, & TRX, T2, X, XO2, XO2L, Y, Z REAL R1MACH, ALNGAM DIMENSION Y(*), TEMP(3) SAVE RTTPI, INLIM DATA RTTPI / 3.98942280401433E-01/ DATA INLIM / 80 / !***FIRST EXECUTABLE STATEMENT BESI NZ = 0 KT = 1 ! I1MACH(15) REPLACES I1MACH(12) IN A DOUBLE PRECISION CODE ! I1MACH(14) REPLACES I1MACH(11) IN A DOUBLE PRECISION CODE RA = R1MACH(3) TOL = MAX(RA,1.0E-15) I1 = -I1MACH(12) GLN = R1MACH(5) ELIM = 2.303E0*(I1*GLN-3.0E0) ! TOLLN = -LN(TOL) I1 = I1MACH(11)+1 TOLLN = 2.303E0*GLN*I1 TOLLN = MIN(TOLLN,34.5388E0) if (N-1) 590, 10, 20 10 KT = 2 20 NN = N if (KODE < 1 .OR. KODE > 2) go to 570 if (X) 600, 30, 80 30 if (ALPHA) 580, 40, 50 40 Y(1) = 1.0E0 if (N == 1) RETURN I1 = 2 go to 60 50 I1 = 1 60 DO 70 I=I1,N Y(I) = 0.0E0 70 CONTINUE return 80 CONTINUE if (ALPHA < 0.0E0) go to 580 ! IALP = INT(ALPHA) FNI = IALP + N - 1 FNF = ALPHA - IALP DFN = FNI + FNF FNU = DFN IN = 0 XO2 = X*0.5E0 SXO2 = XO2*XO2 ETX = KODE - 1 SX = ETX*X ! ! DECISION TREE FOR REGION WHERE SERIES, ASYMPTOTIC EXPANSION FOR X ! TO INFINITY AND ASYMPTOTIC EXPANSION FOR NU TO INFINITY ARE ! APPLIED. ! if (SXO2 <= (FNU+1.0E0)) go to 90 if (X <= 12.0E0) go to 110 FN = 0.55E0*FNU*FNU FN = MAX(17.0E0,FN) if (X >= FN) go to 430 ANS = MAX(36.0E0-FNU,0.0E0) NS = INT(ANS) FNI = FNI + NS DFN = FNI + FNF FN = DFN IS = KT KM = N - 1 + NS if (KM > 0) IS = 3 go to 120 90 FN = FNU FNP1 = FN + 1.0E0 XO2L = LOG(XO2) IS = KT if (X <= 0.5E0) go to 230 NS = 0 100 FNI = FNI + NS DFN = FNI + FNF FN = DFN FNP1 = FN + 1.0E0 IS = KT if (N-1+NS > 0) IS = 3 go to 230 110 XO2L = LOG(XO2) NS = INT(SXO2-FNU) go to 100 120 CONTINUE ! ! OVERFLOW TEST ON UNIFORM ASYMPTOTIC EXPANSION ! if (KODE == 2) go to 130 if (ALPHA < 1.0E0) go to 150 Z = X/ALPHA RA = SQRT(1.0E0+Z*Z) GLN = LOG((1.0E0+RA)/Z) T = RA*(1.0E0-ETX) + ETX/(Z+RA) ARG = ALPHA*(T-GLN) if (ARG > ELIM) go to 610 if (KM == 0) go to 140 130 CONTINUE ! ! UNDERFLOW TEST ON UNIFORM ASYMPTOTIC EXPANSION ! Z = X/FN RA = SQRT(1.0E0+Z*Z) GLN = LOG((1.0E0+RA)/Z) T = RA*(1.0E0-ETX) + ETX/(Z+RA) ARG = FN*(T-GLN) 140 if (ARG < (-ELIM)) go to 280 go to 190 150 if (X > ELIM) go to 610 go to 130 ! ! UNIFORM ASYMPTOTIC EXPANSION FOR NU TO INFINITY ! 160 if (KM /= 0) go to 170 Y(1) = TEMP(3) return 170 TEMP(1) = TEMP(3) IN = NS KT = 1 I1 = 0 180 CONTINUE IS = 2 FNI = FNI - 1.0E0 DFN = FNI + FNF FN = DFN if ( I1 == 2) go to 350 Z = X/FN RA = SQRT(1.0E0+Z*Z) GLN = LOG((1.0E0+RA)/Z) T = RA*(1.0E0-ETX) + ETX/(Z+RA) ARG = FN*(T-GLN) 190 CONTINUE I1 = ABS(3-IS) I1 = MAX(I1,1) FLGIK = 1.0E0 call ASYIK(X,FN,KODE,FLGIK,RA,ARG,I1,TEMP(IS)) go to (180, 350, 510), IS ! ! SERIES FOR (X/2)**2 <= NU+1 ! 230 CONTINUE GLN = ALNGAM(FNP1) ARG = FN*XO2L - GLN - SX if (ARG < (-ELIM)) go to 300 EARG = EXP(ARG) 240 CONTINUE S = 1.0E0 if (X < TOL) go to 260 AK = 3.0E0 T2 = 1.0E0 T = 1.0E0 S1 = FN DO 250 K=1,17 S2 = T2 + S1 T = T*SXO2/S2 S = S + T if (ABS(T) < TOL) go to 260 T2 = T2 + AK AK = AK + 2.0E0 S1 = S1 + FN 250 CONTINUE 260 CONTINUE TEMP(IS) = S*EARG go to (270, 350, 500), IS 270 EARG = EARG*FN/XO2 FNI = FNI - 1.0E0 DFN = FNI + FNF FN = DFN IS = 2 go to 240 ! ! SET UNDERFLOW VALUE AND UPDATE PARAMETERS ! 280 Y(NN) = 0.0E0 NN = NN - 1 FNI = FNI - 1.0E0 DFN = FNI + FNF FN = DFN if (NN-1) 340, 290, 130 290 KT = 2 IS = 2 go to 130 300 Y(NN) = 0.0E0 NN = NN - 1 FNP1 = FN FNI = FNI - 1.0E0 DFN = FNI + FNF FN = DFN if (NN-1) 340, 310, 320 310 KT = 2 IS = 2 320 if (SXO2 <= FNP1) go to 330 go to 130 330 ARG = ARG - XO2L + LOG(FNP1) if (ARG < (-ELIM)) go to 300 go to 230 340 NZ = N - NN return ! ! BACKWARD RECURSION SECTION ! 350 CONTINUE NZ = N - NN 360 CONTINUE if ( KT == 2) go to 420 S1 = TEMP(1) S2 = TEMP(2) TRX = 2.0E0/X DTM = FNI TM = (DTM+FNF)*TRX if (IN == 0) go to 390 ! BACKWARD RECUR TO INDEX ALPHA+NN-1 DO 380 I=1,IN S = S2 S2 = TM*S2 + S1 S1 = S DTM = DTM - 1.0E0 TM = (DTM+FNF)*TRX 380 CONTINUE Y(NN) = S1 if (NN == 1) RETURN Y(NN-1) = S2 if (NN == 2) RETURN go to 400 390 CONTINUE ! BACKWARD RECUR FROM INDEX ALPHA+NN-1 TO ALPHA Y(NN) = S1 Y(NN-1) = S2 if (NN == 2) RETURN 400 K = NN + 1 DO 410 I=3,NN K = K - 1 Y(K-2) = TM*Y(K-1) + Y(K) DTM = DTM - 1.0E0 TM = (DTM+FNF)*TRX 410 CONTINUE return 420 Y(1) = TEMP(2) return ! ! ASYMPTOTIC EXPANSION FOR X TO INFINITY ! 430 CONTINUE EARG = RTTPI/SQRT(X) if (KODE == 2) go to 440 if (X > ELIM) go to 610 EARG = EARG*EXP(X) 440 ETX = 8.0E0*X IS = KT IN = 0 FN = FNU 450 DX = FNI + FNI TM = 0.0E0 if (FNI == 0.0E0 .AND. ABS(FNF) < TOL) go to 460 TM = 4.0E0*FNF*(FNI+FNI+FNF) 460 CONTINUE DTM = DX*DX S1 = ETX TRX = DTM - 1.0E0 DX = -(TRX+TM)/ETX T = DX S = 1.0E0 + DX ATOL = TOL*ABS(S) S2 = 1.0E0 AK = 8.0E0 DO 470 K=1,25 S1 = S1 + ETX S2 = S2 + AK DX = DTM - S2 AP = DX + TM T = -T*AP/S1 S = S + T if (ABS(T) <= ATOL) go to 480 AK = AK + 8.0E0 470 CONTINUE 480 TEMP(IS) = S*EARG if ( IS == 2) go to 360 IS = 2 FNI = FNI - 1.0E0 DFN = FNI + FNF FN = DFN go to 450 ! ! BACKWARD RECURSION WITH NORMALIZATION BY ! ASYMPTOTIC EXPANSION FOR NU TO INFINITY OR POWER SERIES. ! 500 CONTINUE ! COMPUTATION OF LAST ORDER FOR SERIES NORMALIZATION AKM = MAX(3.0E0-FN,0.0E0) KM = INT(AKM) TFN = FN + KM TA = (GLN+TFN-0.9189385332E0-0.0833333333E0/TFN)/(TFN+0.5E0) TA = XO2L - TA TB = -(1.0E0-1.0E0/TFN)/TFN AIN = TOLLN/(-TA+SQRT(TA*TA-TOLLN*TB)) + 1.5E0 IN = INT(AIN) IN = IN + KM go to 520 510 CONTINUE ! COMPUTATION OF LAST ORDER FOR ASYMPTOTIC EXPANSION NORMALIZATION T = 1.0E0/(FN*RA) AIN = TOLLN/(GLN+SQRT(GLN*GLN+T*TOLLN)) + 1.5E0 IN = INT(AIN) if (IN > INLIM) go to 160 520 CONTINUE TRX = 2.0E0/X DTM = FNI + IN TM = (DTM+FNF)*TRX TA = 0.0E0 TB = TOL KK = 1 530 CONTINUE ! ! BACKWARD RECUR UNINDEXED ! DO 540 I=1,IN S = TB TB = TM*TB + TA TA = S DTM = DTM - 1.0E0 TM = (DTM+FNF)*TRX 540 CONTINUE ! NORMALIZATION if (KK /= 1) go to 550 TA = (TA/TB)*TEMP(3) TB = TEMP(3) KK = 2 IN = NS if (NS /= 0) go to 530 550 Y(NN) = TB NZ = N - NN if (NN == 1) RETURN TB = TM*TB + TA K = NN - 1 Y(K) = TB if (NN == 2) RETURN DTM = DTM - 1.0E0 TM = (DTM+FNF)*TRX KM = K - 1 ! ! BACKWARD RECUR INDEXED ! DO 560 I=1,KM Y(K-1) = TM*Y(K) + Y(K+1) DTM = DTM - 1.0E0 TM = (DTM+FNF)*TRX K = K - 1 560 CONTINUE return ! ! ! 570 CONTINUE call XERMSG ('SLATEC', 'BESI', & 'SCALING OPTION, KODE, NOT 1 OR 2.', 2, 1) return 580 CONTINUE call XERMSG ('SLATEC', 'BESI', 'ORDER, ALPHA, LESS THAN ZERO.', & 2, 1) return 590 CONTINUE call XERMSG ('SLATEC', 'BESI', 'N LESS THAN ONE.', 2, 1) return 600 CONTINUE call XERMSG ('SLATEC', 'BESI', 'X LESS THAN ZERO.', 2, 1) return 610 CONTINUE call XERMSG ('SLATEC', 'BESI', & 'OVERFLOW, X TOO LARGE FOR KODE = 1.', 6, 1) return end function BESI0 (X) ! !! BESI0 computes the hyperbolic Bessel function of the first kind ... ! of order zero. ! !***LIBRARY SLATEC (FNLIB) !***CATEGORY C10B1 !***TYPE SINGLE PRECISION (BESI0-S, DBESI0-D) !***KEYWORDS FIRST KIND, FNLIB, HYPERBOLIC BESSEL FUNCTION, ! MODIFIED BESSEL FUNCTION, ORDER ZERO, SPECIAL FUNCTIONS !***AUTHOR Fullerton, W., (LANL) !***DESCRIPTION ! ! BESI0(X) computes the modified (hyperbolic) Bessel function ! of the first kind of order zero and real argument X. ! ! Series for BI0 on the interval 0. to 9.00000D+00 ! with weighted error 2.46E-18 ! log weighted error 17.61 ! significant figures required 17.90 ! decimal places required 18.15 ! !***REFERENCES (NONE) !***ROUTINES CALLED BESI0E, CSEVL, INITS, R1MACH, XERMSG !***REVISION HISTORY (YYMMDD) ! 770401 DATE WRITTEN ! 890531 Changed all specific intrinsics to generic. (WRB) ! 890531 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) ! 900326 Removed duplicate information from DESCRIPTION section. ! (WRB) !***END PROLOGUE BESI0 DIMENSION BI0CS(12) LOGICAL FIRST SAVE BI0CS, NTI0, XSML, XMAX, FIRST DATA BI0CS( 1) / -.07660547252839144951E0 / DATA BI0CS( 2) / 1.927337953993808270E0 / DATA BI0CS( 3) / .2282644586920301339E0 / DATA BI0CS( 4) / .01304891466707290428E0 / DATA BI0CS( 5) / .00043442709008164874E0 / DATA BI0CS( 6) / .00000942265768600193E0 / DATA BI0CS( 7) / .00000014340062895106E0 / DATA BI0CS( 8) / .00000000161384906966E0 / DATA BI0CS( 9) / .00000000001396650044E0 / DATA BI0CS(10) / .00000000000009579451E0 / DATA BI0CS(11) / .00000000000000053339E0 / DATA BI0CS(12) / .00000000000000000245E0 / DATA FIRST /.TRUE./ !***FIRST EXECUTABLE STATEMENT BESI0 if (FIRST) THEN NTI0 = INITS (BI0CS, 12, 0.1*R1MACH(3)) XSML = SQRT (4.5*R1MACH(3)) XMAX = LOG (R1MACH(2)) end if FIRST = .FALSE. ! Y = ABS(X) if (Y > 3.0) go to 20 ! BESI0 = 1.0 if (Y > XSML) BESI0 = 2.75 + CSEVL (Y*Y/4.5-1.0, BI0CS, NTI0) return ! 20 if (Y > XMAX) call XERMSG ('SLATEC', 'BESI0', & 'ABS(X) SO BIG I0 OVERFLOWS', 1, 2) ! BESI0 = EXP(Y) * BESI0E(X) ! return end function BESI0E (X) ! !! BESI0E computes the exponentially scaled modified (hyperbolic) ... ! Bessel function of the first kind of order zero. ! !***LIBRARY SLATEC (FNLIB) !***CATEGORY C10B1 !***TYPE SINGLE PRECISION (BESI0E-S, DBSI0E-D) !***KEYWORDS EXPONENTIALLY SCALED, FIRST KIND, FNLIB, ! HYPERBOLIC BESSEL FUNCTION, MODIFIED BESSEL FUNCTION, ! ORDER ZERO, SPECIAL FUNCTIONS !***AUTHOR Fullerton, W., (LANL) !***DESCRIPTION ! ! BESI0E(X) calculates the exponentially scaled modified (hyperbolic) ! Bessel function of the first kind of order zero for real argument X; ! i.e., EXP(-ABS(X))*I0(X). ! ! ! Series for BI0 on the interval 0. to 9.00000D+00 ! with weighted error 2.46E-18 ! log weighted error 17.61 ! significant figures required 17.90 ! decimal places required 18.15 ! ! ! Series for AI0 on the interval 1.25000D-01 to 3.33333D-01 ! with weighted error 7.87E-17 ! log weighted error 16.10 ! significant figures required 14.69 ! decimal places required 16.76 ! ! ! Series for AI02 on the interval 0. to 1.25000D-01 ! with weighted error 3.79E-17 ! log weighted error 16.42 ! significant figures required 14.86 ! decimal places required 17.09 ! !***REFERENCES (NONE) !***ROUTINES CALLED CSEVL, INITS, R1MACH !***REVISION HISTORY (YYMMDD) ! 770701 DATE WRITTEN ! 890313 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) !***END PROLOGUE BESI0E DIMENSION BI0CS(12), AI0CS(21), AI02CS(22) LOGICAL FIRST SAVE BI0CS, AI0CS, AI02CS, NTI0, NTAI0, NTAI02, XSML, FIRST DATA BI0CS( 1) / -.07660547252839144951E0 / DATA BI0CS( 2) / 1.927337953993808270E0 / DATA BI0CS( 3) / .2282644586920301339E0 / DATA BI0CS( 4) / .01304891466707290428E0 / DATA BI0CS( 5) / .00043442709008164874E0 / DATA BI0CS( 6) / .00000942265768600193E0 / DATA BI0CS( 7) / .00000014340062895106E0 / DATA BI0CS( 8) / .00000000161384906966E0 / DATA BI0CS( 9) / .00000000001396650044E0 / DATA BI0CS(10) / .00000000000009579451E0 / DATA BI0CS(11) / .00000000000000053339E0 / DATA BI0CS(12) / .00000000000000000245E0 / DATA AI0CS( 1) / .07575994494023796E0 / DATA AI0CS( 2) / .00759138081082334E0 / DATA AI0CS( 3) / .00041531313389237E0 / DATA AI0CS( 4) / .00001070076463439E0 / DATA AI0CS( 5) / -.00000790117997921E0 / DATA AI0CS( 6) / -.00000078261435014E0 / DATA AI0CS( 7) / .00000027838499429E0 / DATA AI0CS( 8) / .00000000825247260E0 / DATA AI0CS( 9) / -.00000001204463945E0 / DATA AI0CS(10) / .00000000155964859E0 / DATA AI0CS(11) / .00000000022925563E0 / DATA AI0CS(12) / -.00000000011916228E0 / DATA AI0CS(13) / .00000000001757854E0 / DATA AI0CS(14) / .00000000000112822E0 / DATA AI0CS(15) / -.00000000000114684E0 / DATA AI0CS(16) / .00000000000027155E0 / DATA AI0CS(17) / -.00000000000002415E0 / DATA AI0CS(18) / -.00000000000000608E0 / DATA AI0CS(19) / .00000000000000314E0 / DATA AI0CS(20) / -.00000000000000071E0 / DATA AI0CS(21) / .00000000000000007E0 / DATA AI02CS( 1) / .05449041101410882E0 / DATA AI02CS( 2) / .00336911647825569E0 / DATA AI02CS( 3) / .00006889758346918E0 / DATA AI02CS( 4) / .00000289137052082E0 / DATA AI02CS( 5) / .00000020489185893E0 / DATA AI02CS( 6) / .00000002266668991E0 / DATA AI02CS( 7) / .00000000339623203E0 / DATA AI02CS( 8) / .00000000049406022E0 / DATA AI02CS( 9) / .00000000001188914E0 / DATA AI02CS(10) / -.00000000003149915E0 / DATA AI02CS(11) / -.00000000001321580E0 / DATA AI02CS(12) / -.00000000000179419E0 / DATA AI02CS(13) / .00000000000071801E0 / DATA AI02CS(14) / .00000000000038529E0 / DATA AI02CS(15) / .00000000000001539E0 / DATA AI02CS(16) / -.00000000000004151E0 / DATA AI02CS(17) / -.00000000000000954E0 / DATA AI02CS(18) / .00000000000000382E0 / DATA AI02CS(19) / .00000000000000176E0 / DATA AI02CS(20) / -.00000000000000034E0 / DATA AI02CS(21) / -.00000000000000027E0 / DATA AI02CS(22) / .00000000000000003E0 / DATA FIRST /.TRUE./ !***FIRST EXECUTABLE STATEMENT BESI0E if (FIRST) THEN NTI0 = INITS (BI0CS, 12, 0.1*R1MACH(3)) NTAI0 = INITS (AI0CS, 21, 0.1*R1MACH(3)) NTAI02 = INITS (AI02CS, 22, 0.1*R1MACH(3)) XSML = SQRT (4.5*R1MACH(3)) end if FIRST = .FALSE. ! Y = ABS(X) if (Y > 3.0) go to 20 ! BESI0E = 1.0 - X if (Y > XSML) BESI0E = EXP(-Y) * ( 2.75 + & CSEVL (Y*Y/4.5-1.0, BI0CS, NTI0) ) return ! 20 if (Y <= 8.) BESI0E = (.375 + CSEVL ((48./Y-11.)/5., AI0CS, NTAI0) & ) / SQRT(Y) if (Y > 8.) BESI0E = (.375 + CSEVL (16./Y-1., AI02CS, NTAI02)) & / SQRT(Y) ! return end function BESI1 (X) ! !! BESI1 computes the modified (hyperbolic) Bessel function of the ... ! first kind of order one. ! !***LIBRARY SLATEC (FNLIB) !***CATEGORY C10B1 !***TYPE SINGLE PRECISION (BESI1-S, DBESI1-D) !***KEYWORDS FIRST KIND, FNLIB, HYPERBOLIC BESSEL FUNCTION, ! MODIFIED BESSEL FUNCTION, ORDER ONE, SPECIAL FUNCTIONS !***AUTHOR Fullerton, W., (LANL) !***DESCRIPTION ! ! BESI1(X) calculates the modified (hyperbolic) Bessel function ! of the first kind of order one for real argument X. ! ! Series for BI1 on the interval 0. to 9.00000D+00 ! with weighted error 2.40E-17 ! log weighted error 16.62 ! significant figures required 16.23 ! decimal places required 17.14 ! !***REFERENCES (NONE) !***ROUTINES CALLED BESI1E, CSEVL, INITS, R1MACH, XERMSG !***REVISION HISTORY (YYMMDD) ! 770401 DATE WRITTEN ! 890531 Changed all specific intrinsics to generic. (WRB) ! 890531 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) ! 900326 Removed duplicate information from DESCRIPTION section. ! (WRB) !***END PROLOGUE BESI1 DIMENSION BI1CS(11) LOGICAL FIRST SAVE BI1CS, NTI1, XMIN, XSML, XMAX, FIRST DATA BI1CS( 1) / -.001971713261099859E0 / DATA BI1CS( 2) / .40734887667546481E0 / DATA BI1CS( 3) / .034838994299959456E0 / DATA BI1CS( 4) / .001545394556300123E0 / DATA BI1CS( 5) / .000041888521098377E0 / DATA BI1CS( 6) / .000000764902676483E0 / DATA BI1CS( 7) / .000000010042493924E0 / DATA BI1CS( 8) / .000000000099322077E0 / DATA BI1CS( 9) / .000000000000766380E0 / DATA BI1CS(10) / .000000000000004741E0 / DATA BI1CS(11) / .000000000000000024E0 / DATA FIRST /.TRUE./ !***FIRST EXECUTABLE STATEMENT BESI1 if (FIRST) THEN NTI1 = INITS (BI1CS, 11, 0.1*R1MACH(3)) XMIN = 2.0*R1MACH(1) XSML = SQRT (4.5*R1MACH(3)) XMAX = LOG (R1MACH(2)) end if FIRST = .FALSE. ! Y = ABS(X) if (Y > 3.0) go to 20 ! BESI1 = 0.0 if (Y == 0.0) return ! if (Y <= XMIN) call XERMSG ('SLATEC', 'BESI1', & 'ABS(X) SO SMALL I1 UNDERFLOWS', 1, 1) if (Y > XMIN) BESI1 = 0.5*X if (Y > XSML) BESI1 = X * (.875 + CSEVL(Y*Y/4.5-1., BI1CS, NTI1)) return ! 20 if (Y > XMAX) call XERMSG ('SLATEC', 'BESI1', & 'ABS(X) SO BIG I1 OVERFLOWS', 2, 2) ! BESI1 = EXP(Y) * BESI1E(X) ! return end function BESI1E (X) ! !! BESI1E computes the exponentially scaled modified (hyperbolic) ... ! Bessel function of the first kind of order one. ! !***LIBRARY SLATEC (FNLIB) !***CATEGORY C10B1 !***TYPE SINGLE PRECISION (BESI1E-S, DBSI1E-D) !***KEYWORDS EXPONENTIALLY SCALED, FIRST KIND, FNLIB, ! HYPERBOLIC BESSEL FUNCTION, MODIFIED BESSEL FUNCTION, ! ORDER ONE, SPECIAL FUNCTIONS !***AUTHOR Fullerton, W., (LANL) !***DESCRIPTION ! ! BESI1E(X) calculates the exponentially scaled modified (hyperbolic) ! Bessel function of the first kind of order one for real argument X; ! i.e., EXP(-ABS(X))*I1(X). ! ! Series for BI1 on the interval 0. to 9.00000D+00 ! with weighted error 2.40E-17 ! log weighted error 16.62 ! significant figures required 16.23 ! decimal places required 17.14 ! ! Series for AI1 on the interval 1.25000D-01 to 3.33333D-01 ! with weighted error 6.98E-17 ! log weighted error 16.16 ! significant figures required 14.53 ! decimal places required 16.82 ! ! Series for AI12 on the interval 0. to 1.25000D-01 ! with weighted error 3.55E-17 ! log weighted error 16.45 ! significant figures required 14.69 ! decimal places required 17.12 ! !***REFERENCES (NONE) !***ROUTINES CALLED CSEVL, INITS, R1MACH, XERMSG !***REVISION HISTORY (YYMMDD) ! 770401 DATE WRITTEN ! 890210 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) ! 900326 Removed duplicate information from DESCRIPTION section. ! (WRB) ! 920618 Removed space from variable names. (RWC, WRB) !***END PROLOGUE BESI1E DIMENSION BI1CS(11), AI1CS(21), AI12CS(22) LOGICAL FIRST SAVE BI1CS, AI1CS, AI12CS, NTI1, NTAI1, NTAI12, XMIN, XSML, FIRST DATA BI1CS( 1) / -.001971713261099859E0 / DATA BI1CS( 2) / .40734887667546481E0 / DATA BI1CS( 3) / .034838994299959456E0 / DATA BI1CS( 4) / .001545394556300123E0 / DATA BI1CS( 5) / .000041888521098377E0 / DATA BI1CS( 6) / .000000764902676483E0 / DATA BI1CS( 7) / .000000010042493924E0 / DATA BI1CS( 8) / .000000000099322077E0 / DATA BI1CS( 9) / .000000000000766380E0 / DATA BI1CS(10) / .000000000000004741E0 / DATA BI1CS(11) / .000000000000000024E0 / DATA AI1CS( 1) / -.02846744181881479E0 / DATA AI1CS( 2) / -.01922953231443221E0 / DATA AI1CS( 3) / -.00061151858579437E0 / DATA AI1CS( 4) / -.00002069971253350E0 / DATA AI1CS( 5) / .00000858561914581E0 / DATA AI1CS( 6) / .00000104949824671E0 / DATA AI1CS( 7) / -.00000029183389184E0 / DATA AI1CS( 8) / -.00000001559378146E0 / DATA AI1CS( 9) / .00000001318012367E0 / DATA AI1CS(10) / -.00000000144842341E0 / DATA AI1CS(11) / -.00000000029085122E0 / DATA AI1CS(12) / .00000000012663889E0 / DATA AI1CS(13) / -.00000000001664947E0 / DATA AI1CS(14) / -.00000000000166665E0 / DATA AI1CS(15) / .00000000000124260E0 / DATA AI1CS(16) / -.00000000000027315E0 / DATA AI1CS(17) / .00000000000002023E0 / DATA AI1CS(18) / .00000000000000730E0 / DATA AI1CS(19) / -.00000000000000333E0 / DATA AI1CS(20) / .00000000000000071E0 / DATA AI1CS(21) / -.00000000000000006E0 / DATA AI12CS( 1) / .02857623501828014E0 / DATA AI12CS( 2) / -.00976109749136147E0 / DATA AI12CS( 3) / -.00011058893876263E0 / DATA AI12CS( 4) / -.00000388256480887E0 / DATA AI12CS( 5) / -.00000025122362377E0 / DATA AI12CS( 6) / -.00000002631468847E0 / DATA AI12CS( 7) / -.00000000383538039E0 / DATA AI12CS( 8) / -.00000000055897433E0 / DATA AI12CS( 9) / -.00000000001897495E0 / DATA AI12CS(10) / .00000000003252602E0 / DATA AI12CS(11) / .00000000001412580E0 / DATA AI12CS(12) / .00000000000203564E0 / DATA AI12CS(13) / -.00000000000071985E0 / DATA AI12CS(14) / -.00000000000040836E0 / DATA AI12CS(15) / -.00000000000002101E0 / DATA AI12CS(16) / .00000000000004273E0 / DATA AI12CS(17) / .00000000000001041E0 / DATA AI12CS(18) / -.00000000000000382E0 / DATA AI12CS(19) / -.00000000000000186E0 / DATA AI12CS(20) / .00000000000000033E0 / DATA AI12CS(21) / .00000000000000028E0 / DATA AI12CS(22) / -.00000000000000003E0 / DATA FIRST /.TRUE./ !***FIRST EXECUTABLE STATEMENT BESI1E if (FIRST) THEN NTI1 = INITS (BI1CS, 11, 0.1*R1MACH(3)) NTAI1 = INITS (AI1CS, 21, 0.1*R1MACH(3)) NTAI12 = INITS (AI12CS, 22, 0.1*R1MACH(3)) ! XMIN = 2.0*R1MACH(1) XSML = SQRT (4.5*R1MACH(3)) end if FIRST = .FALSE. ! Y = ABS(X) if (Y > 3.0) go to 20 ! BESI1E = 0.0 if (Y == 0.0) return ! if (Y <= XMIN) call XERMSG ('SLATEC', 'BESI1E', & 'ABS(X) SO SMALL I1 UNDERFLOWS', 1, 1) if (Y > XMIN) BESI1E = 0.5*X if (Y > XSML) BESI1E = X * (.875 + CSEVL(Y*Y/4.5-1., BI1CS,NTI1)) BESI1E = EXP(-Y) * BESI1E return ! 20 if (Y <= 8.) BESI1E = (.375 + CSEVL ((48./Y-11.)/5., AI1CS, NTAI1) & ) / SQRT(Y) if (Y > 8.) BESI1E = (.375 + CSEVL (16./Y-1.0, AI12CS, NTAI12)) & / SQRT(Y) BESI1E = SIGN (BESI1E, X) ! return end subroutine BESJ (X, ALPHA, N, Y, NZ) ! !! BESJ computes an N member sequence of J Bessel functions ... ! J/SUB(ALPHA+K-1)/(X), K=1,...,N for non-negative ALPHA and X. ! !***LIBRARY SLATEC !***CATEGORY C10A3 !***TYPE SINGLE PRECISION (BESJ-S, DBESJ-D) !***KEYWORDS J BESSEL FUNCTION, SPECIAL FUNCTIONS !***AUTHOR Amos, D. E., (SNLA) ! Daniel, S. L., (SNLA) ! Weston, M. K., (SNLA) !***DESCRIPTION ! ! Abstract ! BESJ computes an N member sequence of J Bessel functions ! J/sub(ALPHA+K-1)/(X), K=1,...,N for non-negative ALPHA and X. ! A combination of the power series, the asymptotic expansion ! for X to infinity and the uniform asymptotic expansion for ! NU to infinity are applied over subdivisions of the (NU,X) ! plane. For values of (NU,X) not covered by one of these ! formulae, the order is incremented or decremented by integer ! values into a region where one of the formulae apply. Backward ! recursion is applied to reduce orders by integer values except ! where the entire sequence lies in the oscillatory region. In ! this case forward recursion is stable and values from the ! asymptotic expansion for X to infinity start the recursion ! when it is efficient to do so. Leading terms of the series ! and uniform expansion are tested for underflow. If a sequence ! is requested and the last member would underflow, the result ! is set to zero and the next lower order tried, etc., until a ! member comes on scale or all members are set to zero. ! Overflow cannot occur. ! ! Description of Arguments ! ! Input ! X - X >= 0.0E0 ! ALPHA - order of first member of the sequence, ! ALPHA >= 0.0E0 ! N - number of members in the sequence, N >= 1 ! ! Output ! Y - a vector whose first N components contain ! values for J/sub(ALPHA+K-1)/(X), K=1,...,N ! NZ - number of components of Y set to zero due to ! underflow, ! NZ=0 , normal return, computation completed ! NZ /= 0, last NZ components of Y set to zero, ! Y(K)=0.0E0, K=N-NZ+1,...,N. ! ! Error Conditions ! Improper input arguments - a fatal error ! Underflow - a non-fatal error (NZ /= 0) ! !***REFERENCES D. E. Amos, S. L. Daniel and M. K. Weston, CDC 6600 ! subroutines IBESS and JBESS for Bessel functions ! I(NU,X) and J(NU,X), X >= 0, NU >= 0, ACM ! Transactions on Mathematical Software 3, (1977), ! pp. 76-92. ! F. W. J. Olver, Tables of Bessel Functions of Moderate ! or Large Orders, NPL Mathematical Tables 6, Her ! Majesty's Stationery Office, London, 1962. !***ROUTINES CALLED ALNGAM, ASYJY, I1MACH, JAIRY, R1MACH, XERMSG !***REVISION HISTORY (YYMMDD) ! 750101 DATE WRITTEN ! 890531 Changed all specific intrinsics to generic. (WRB) ! 890531 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) ! 900326 Removed duplicate information from DESCRIPTION section. ! (WRB) ! 920501 Reformatted the REFERENCES section. (WRB) !***END PROLOGUE BESJ EXTERNAL JAIRY INTEGER I,IALP,IDALP,IFLW,IN,INLIM,IS,I1,I2,K,KK,KM,KT,N,NN, & NS,NZ INTEGER I1MACH REAL AK,AKM,ALPHA,ANS,AP,ARG,COEF,DALPHA,DFN,DTM,EARG, & ELIM1,ETX,FIDAL,FLGJY,FN,FNF,FNI,FNP1,FNU,FNULIM, & GLN,PDF,PIDT,PP,RDEN,RELB,RTTP,RTWO,RTX,RZDEN, & S,SA,SB,SXO2,S1,S2,T,TA,TAU,TB,TEMP,TFN,TM,TOL, & TOLLN,TRX,TX,T1,T2,WK,X,XO2,XO2L,Y,RTOL,SLIM SAVE RTWO, PDF, RTTP, PIDT, PP, INLIM, FNULIM REAL R1MACH, ALNGAM DIMENSION Y(*), TEMP(3), FNULIM(2), PP(4), WK(7) DATA RTWO,PDF,RTTP,PIDT / 1.34839972492648E+00, & 7.85398163397448E-01, 7.97884560802865E-01, 1.57079632679490E+00/ DATA PP(1), PP(2), PP(3), PP(4) / 8.72909153935547E+00, & 2.65693932265030E-01, 1.24578576865586E-01, 7.70133747430388E-04/ DATA INLIM / 150 / DATA FNULIM(1), FNULIM(2) / 100.0E0, 60.0E0 / !***FIRST EXECUTABLE STATEMENT BESJ NZ = 0 KT = 1 NS=0 ! I1MACH(14) REPLACES I1MACH(11) IN A DOUBLE PRECISION CODE ! I1MACH(15) REPLACES I1MACH(12) IN A DOUBLE PRECISION CODE TA = R1MACH(3) TOL = MAX(TA,1.0E-15) I1 = I1MACH(11) + 1 I2 = I1MACH(12) TB = R1MACH(5) ELIM1 = -2.303E0*(I2*TB+3.0E0) RTOL=1.0E0/TOL SLIM=R1MACH(1)*1.0E+3*RTOL ! TOLLN = -LN(TOL) TOLLN = 2.303E0*TB*I1 TOLLN = MIN(TOLLN,34.5388E0) if (N-1) 720, 10, 20 10 KT = 2 20 NN = N if (X) 730, 30, 80 30 if (ALPHA) 710, 40, 50 40 Y(1) = 1.0E0 if (N == 1) RETURN I1 = 2 go to 60 50 I1 = 1 60 DO 70 I=I1,N Y(I) = 0.0E0 70 CONTINUE return 80 CONTINUE if (ALPHA < 0.0E0) go to 710 ! IALP = INT(ALPHA) FNI = IALP + N - 1 FNF = ALPHA - IALP DFN = FNI + FNF FNU = DFN XO2 = X*0.5E0 SXO2 = XO2*XO2 ! ! DECISION TREE FOR REGION WHERE SERIES, ASYMPTOTIC EXPANSION FOR X ! TO INFINITY AND ASYMPTOTIC EXPANSION FOR NU TO INFINITY ARE ! APPLIED. ! if (SXO2 <= (FNU+1.0E0)) go to 90 TA = MAX(20.0E0,FNU) if (X > TA) go to 120 if (X > 12.0E0) go to 110 XO2L = LOG(XO2) NS = INT(SXO2-FNU) + 1 go to 100 90 FN = FNU FNP1 = FN + 1.0E0 XO2L = LOG(XO2) IS = KT if (X <= 0.50E0) go to 330 NS = 0 100 FNI = FNI + NS DFN = FNI + FNF FN = DFN FNP1 = FN + 1.0E0 IS = KT if (N-1+NS > 0) IS = 3 go to 330 110 ANS = MAX(36.0E0-FNU,0.0E0) NS = INT(ANS) FNI = FNI + NS DFN = FNI + FNF FN = DFN IS = KT if (N-1+NS > 0) IS = 3 go to 130 120 CONTINUE RTX = SQRT(X) TAU = RTWO*RTX TA = TAU + FNULIM(KT) if (FNU <= TA) go to 480 FN = FNU IS = KT ! ! UNIFORM ASYMPTOTIC EXPANSION FOR NU TO INFINITY ! 130 CONTINUE I1 = ABS(3-IS) I1 = MAX(I1,1) FLGJY = 1.0E0 call ASYJY(JAIRY,X,FN,FLGJY,I1,TEMP(IS),WK,IFLW) if ( IFLW /= 0) go to 380 go to (320, 450, 620), IS 310 TEMP(1) = TEMP(3) KT = 1 320 IS = 2 FNI = FNI - 1.0E0 DFN = FNI + FNF FN = DFN if ( I1 == 2) go to 450 go to 130 ! ! SERIES FOR (X/2)**2 <= NU+1 ! 330 CONTINUE GLN = ALNGAM(FNP1) ARG = FN*XO2L - GLN if (ARG < (-ELIM1)) go to 400 EARG = EXP(ARG) 340 CONTINUE S = 1.0E0 if (X < TOL) go to 360 AK = 3.0E0 T2 = 1.0E0 T = 1.0E0 S1 = FN DO 350 K=1,17 S2 = T2 + S1 T = -T*SXO2/S2 S = S + T if (ABS(T) < TOL) go to 360 T2 = T2 + AK AK = AK + 2.0E0 S1 = S1 + FN 350 CONTINUE 360 CONTINUE TEMP(IS) = S*EARG go to (370, 450, 610), IS 370 EARG = EARG*FN/XO2 FNI = FNI - 1.0E0 DFN = FNI + FNF FN = DFN IS = 2 go to 340 ! ! SET UNDERFLOW VALUE AND UPDATE PARAMETERS ! UNDERFLOW CAN ONLY OCCUR FOR NS=0 SINCE THE ORDER MUST BE ! LARGER THAN 36. THEREFORE, NS NEED NOT BE CONSIDERED. ! 380 Y(NN) = 0.0E0 NN = NN - 1 FNI = FNI - 1.0E0 DFN = FNI + FNF FN = DFN if (NN-1) 440, 390, 130 390 KT = 2 IS = 2 go to 130 400 Y(NN) = 0.0E0 NN = NN - 1 FNP1 = FN FNI = FNI - 1.0E0 DFN = FNI + FNF FN = DFN if (NN-1) 440, 410, 420 410 KT = 2 IS = 2 420 if (SXO2 <= FNP1) go to 430 go to 130 430 ARG = ARG - XO2L + LOG(FNP1) if (ARG < (-ELIM1)) go to 400 go to 330 440 NZ = N - NN return ! ! BACKWARD RECURSION SECTION ! 450 CONTINUE if ( NS /= 0) go to 451 NZ = N - NN if (KT == 2) go to 470 ! BACKWARD RECUR FROM INDEX ALPHA+NN-1 TO ALPHA Y(NN) = TEMP(1) Y(NN-1) = TEMP(2) if (NN == 2) RETURN 451 CONTINUE TRX = 2.0E0/X DTM = FNI TM = (DTM+FNF)*TRX AK=1.0E0 TA=TEMP(1) TB=TEMP(2) if ( ABS(TA) > SLIM) go to 455 TA=TA*RTOL TB=TB*RTOL AK=TOL 455 CONTINUE KK=2 IN=NS-1 if ( IN == 0) go to 690 if ( NS /= 0) go to 670 K=NN-2 DO 460 I=3,NN S=TB TB=TM*TB-TA TA=S Y(K)=TB*AK K=K-1 DTM = DTM - 1.0E0 TM = (DTM+FNF)*TRX 460 CONTINUE return 470 Y(1) = TEMP(2) return ! ! ASYMPTOTIC EXPANSION FOR X TO INFINITY WITH FORWARD RECURSION IN ! OSCILLATORY REGION X > MAX(20, NU), PROVIDED THE LAST MEMBER ! OF THE SEQUENCE IS ALSO IN THE REGION. ! 480 CONTINUE IN = INT(ALPHA-TAU+2.0E0) if (IN <= 0) go to 490 IDALP = IALP - IN - 1 KT = 1 go to 500 490 CONTINUE IDALP = IALP IN = 0 500 IS = KT FIDAL = IDALP DALPHA = FIDAL + FNF ARG = X - PIDT*DALPHA - PDF SA = SIN(ARG) SB = COS(ARG) COEF = RTTP/RTX ETX = 8.0E0*X 510 CONTINUE DTM = FIDAL + FIDAL DTM = DTM*DTM TM = 0.0E0 if (FIDAL == 0.0E0 .AND. ABS(FNF) < TOL) go to 520 TM = 4.0E0*FNF*(FIDAL+FIDAL+FNF) 520 CONTINUE TRX = DTM - 1.0E0 T2 = (TRX+TM)/ETX S2 = T2 RELB = TOL*ABS(T2) T1 = ETX S1 = 1.0E0 FN = 1.0E0 AK = 8.0E0 DO 530 K=1,13 T1 = T1 + ETX FN = FN + AK TRX = DTM -