Logo Search packages:      
Sourcecode: fbasics version File versions  Download package

13B-HyperbolicDistribution.f

C
C
C======================================================================
C PACKAGE DENSITIES: DNIG
C======================================================================
C
      SUBROUTINE DNIG(density, x, n, alpha, beta, delta, mu)
C 
C     A function written by D. Wuertz
C     Purpose:
C       Return Hyperbolic Density Function PDF:
C       Needs: BESK1 Modified Bessel Function K1
C         |beta| <= alpha  - Shape Parameters
C         0 <= delta       - Scale Parameter
C         mu               - Location Parameter
C 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER i,n
      DIMENSION density(n), x(n) 
      DOUBLE PRECISION mu
      pi = 4.0d0*datan(1.0d0)
      do i=1,n
        xarg =alpha * dsqrt( delta**2 + (x(i)-mu)**2 )  
      density(i) = (alpha*delta/pi) *
     >    dexp(delta*dsqrt( alpha**2-beta**2 ) + beta*(x(i)-mu)) *
     >    BESK1(xarg) / dsqrt( delta**2 + (x(i)-mu)**2 )
C
      enddo
      END
C
C**********************************************************************
C
C     PACKAGE BESSEL FUNCTION K1
      SUBROUTINE CALCK1(ARG,RESULT,JINT)
C
C   This packet computes modified Bessel functions of the second kind
C   and order one,  K1(X)  and  EXP(X)*K1(X), for real arguments X.
C   It contains two function type subprograms, BESK1  and  BESEK1,
C   and one subroutine type subprogram, CALCK1.  The calling
C   statements for the primary entries are
C                   Y=BESK1(X)
C   and
C                   Y=BESEK1(X)
C   where the entry points correspond to the functions K1(X) and
C   EXP(X)*K1(X), respectively.  The routine CALCK1 is intended
C   for internal packet use only, all computations within the
C   packet being concentrated in this routine.  The function
C   subprograms invoke CALCK1 with the statement
C          CALL CALCK1(ARG,RESULT,JINT)
C   where the parameter usage is as follows
C
C      Function                      Parameters for CALCK1
C        Call             ARG                  RESULT          JINT
C     BESK1(ARG)  XLEAST .LT. ARG .LT. XMAX    K1(ARG)          1
C     BESEK1(ARG)     XLEAST .LT. ARG       EXP(ARG)*K1(ARG)    2
C
C   The main computation evaluates slightly modified forms of near 
C   minimax rational approximations generated by Russon and Blair, 
C   Chalk River (Atomic Energy of Canada Limited) Report AECL-3461, 
C   1969.  This transportable program is patterned after the 
C   machine-dependent FUNPACK packet NATSK1, but cannot match that
C   version for efficiency or accuracy.  This version uses rational
C   functions that theoretically approximate K-SUB-1(X) to at
C   least 18 significant decimal digits.  The accuracy achieved
C   depends on the arithmetic system, the compiler, the intrinsic
C   functions, and proper selection of the machine-dependent
C   constants.
C Explanation of machine-dependent constants
C   beta   = Radix for the floating-point system
C   minexp = Smallest representable power of beta
C   maxexp = Smallest power of beta that overflows
C   XLEAST = Smallest acceptable argument, i.e., smallest machine
C            number X such that 1/X is machine representable.
C   XSMALL = Argument below which BESK1(X) and BESEK1(X) may
C            each be represented by 1/X.  A safe value is the
C            largest X such that  1.0 + X = 1.0  to machine
C            precision.
C   XINF   = Largest positive machine number; approximately
C            beta**maxexp
C   XMAX   = Largest argument acceptable to BESK1;  Solution to
C            equation:  
C               W(X) * (1+3/8X-15/128X**2) = beta**minexp
C            where  W(X) = EXP(-X)*SQRT(PI/2X)
C
C     Approximate values for some important machines are:
C                           beta       minexp       maxexp
C  CRAY-1        (S.P.)       2        -8193         8191
C  Cyber 180/185 
C    under NOS   (S.P.)       2         -975         1070
C  IEEE (IBM/XT,
C    SUN, etc.)  (S.P.)       2         -126          128
C  IEEE (IBM/XT,
C    SUN, etc.)  (D.P.)       2        -1022         1024
C  IBM 3033      (D.P.)      16          -65           63
C  VAX D-Format  (D.P.)       2         -128          127
C  VAX G-Format  (D.P.)       2        -1024         1023
C
C                         XLEAST     XSMALL      XINF       XMAX 
C CRAY-1                1.84E-2466  3.55E-15  5.45E+2465  5674.858 
C Cyber 180/855
C   under NOS   (S.P.)  3.14E-294   1.77E-15  1.26E+322    672.789
C IEEE (IBM/XT,
C   SUN, etc.)  (S.P.)  1.18E-38    5.95E-8   3.40E+38      85.343
C IEEE (IBM/XT,
C   SUN, etc.)  (D.P.)  2.23D-308   1.11D-16  1.79D+308    705.343
C IBM 3033      (D.P.)  1.39D-76    1.11D-16  7.23D+75     177.855
C VAX D-Format  (D.P.)  5.88D-39    6.95D-18  1.70D+38      86.721
C VAX G-Format  (D.P.)  1.12D-308   5.55D-17  8.98D+307    706.728
C
C Error returns
C  The program returns the value XINF for ARG .LE. 0.0 and the
C   BESK1 entry returns the value 0.0 for ARG .GT. XMAX.
C  Intrinsic functions required are:
C     LOG, SQRT, EXP
C  Authors: W. J. Cody and Laura Stoltz
C           Mathematics and Computer Science Division
C           Argonne National Laboratory
C           Argonne, IL 60439
C  Latest modification: January 28, 1988
C
      INTEGER I,JINT
      DOUBLE PRECISION 
     1    ARG,F,G,ONE,P,PP,Q,QQ,RESULT,SUMF,SUMG,
     2    SUMP,SUMQ,X,XINF,XMAX,XLEAST,XSMALL,XX,ZERO
      DIMENSION P(5),Q(3),PP(11),QQ(9),F(5),G(3)
C  Mathematical constants
      DATA ONE/1.0D0/,ZERO/0.0D0/
C  Machine-dependent constants
      DATA XLEAST/2.23D-308/,XSMALL/1.11D-16/,XINF/1.79D+308/,
     1     XMAX/705.343D+0/
C  Coefficients for  XLEAST .LE.  ARG  .LE. 1.0
      DATA   P/ 4.8127070456878442310D-1, 9.9991373567429309922D+1,
     1          7.1885382604084798576D+3, 1.7733324035147015630D+5,
     2          7.1938920065420586101D+5/
      DATA   Q/-2.8143915754538725829D+2, 3.7264298672067697862D+4,
     1         -2.2149374878243304548D+6/
      DATA   F/-2.2795590826955002390D-1,-5.3103913335180275253D+1,
     1         -4.5051623763436087023D+3,-1.4758069205414222471D+5,
     2         -1.3531161492785421328D+6/
      DATA   G/-3.0507151578787595807D+2, 4.3117653211351080007D+4,
     2         -2.7062322985570842656D+6/
C  Coefficients for  1.0 .LT.  ARG
      DATA  PP/ 6.4257745859173138767D-2, 7.5584584631176030810D+0,
     1          1.3182609918569941308D+2, 8.1094256146537402173D+2,
     2          2.3123742209168871550D+3, 3.4540675585544584407D+3,
     3          2.8590657697910288226D+3, 1.3319486433183221990D+3,
     4          3.4122953486801312910D+2, 4.4137176114230414036D+1,
     5          2.2196792496874548962D+0/
      DATA  QQ/ 3.6001069306861518855D+1, 3.3031020088765390854D+2,
     1          1.2082692316002348638D+3, 2.1181000487171943810D+3,
     2          1.9448440788918006154D+3, 9.6929165726802648634D+2,
     3          2.5951223655579051357D+2, 3.4552228452758912848D+1,
     4          1.7710478032601086579D+0/
      X = ARG
      IF (X .LT. XLEAST) THEN
C  Error return for  ARG  .LT. XLEAST
            RESULT = XINF
         ELSE IF (X .LE. ONE) THEN
C  XLEAST .LE.  ARG  .LE. 1.0
            IF (X .LT. XSMALL) THEN
C  Return for small ARG
                  RESULT = ONE / X
               ELSE
                  XX = X * X
                  SUMP = ((((P(1)*XX + P(2))*XX + P(3))*XX + P(4))*XX
     1                   + P(5))*XX + Q(3)
                  SUMQ = ((XX + Q(1))*XX + Q(2))*XX + Q(3)
                  SUMF = (((F(1)*XX + F(2))*XX + F(3))*XX + F(4))*XX
     1                   + F(5)
                  SUMG = ((XX + G(1))*XX + G(2))*XX + G(3)
                  RESULT = (XX * LOG(X) * SUMF/SUMG + SUMP/SUMQ) / X
                  IF (JINT .EQ. 2) RESULT = RESULT * EXP(X)
            END IF
         ELSE IF ((JINT .EQ. 1) .AND. (X .GT. XMAX)) THEN
C  Error return for  ARG  .GT. XMAX
            RESULT = ZERO
         ELSE
C  1.0 .LT.  ARG
            XX = ONE / X
            SUMP = PP(1)
            DO 120 I = 2, 11
               SUMP = SUMP * XX + PP(I)
  120       CONTINUE
            SUMQ = XX
            DO 140 I = 1, 8
               SUMQ = (SUMQ + QQ(I)) * XX
  140       CONTINUE
            SUMQ = SUMQ + QQ(9)
            RESULT = SUMP / SUMQ / SQRT(X)
            IF (JINT .EQ. 1) RESULT = RESULT * EXP(-X)
      END IF
      RETURN
      END
C
C**********************************************************************
C
      DOUBLE PRECISION FUNCTION BESK1(X)
C
C     This function program computes approximate values for the
C     modified Bessel function of the second kind of order one
C     for arguments  XLEAST .LE. ARG .LE. XMAX.
C
      INTEGER JINT
      DOUBLE PRECISION X, RESULT
      JINT = 1
      CALL CALCK1(X,RESULT,JINT)
      BESK1 = RESULT
      RETURN
      END
C
C**********************************************************************
C
      DOUBLE PRECISION FUNCTION BESEK1(X)
C
C     This function program computes approximate values for the
C     modified Bessel function of the second kind of order one
C     multiplied by the exponential function, for arguments
C     XLEAST .LE. ARG .LE. XMAX.
C
      INTEGER JINT
      DOUBLE PRECISION X,RESULT
      JINT = 2
      CALL CALCK1(X,RESULT,JINT)
      BESEK1 = RESULT
      RETURN
      END
C
C
C
C======================================================================
C PACKAGE PROBABILITIES: PNIG
C======================================================================
C
C
C
C     SMI INFINITE INTEGRATION      
      SUBROUTINE DEHINT (FUNC, A, EPS, V)
C
C     Integrate FUNC over (A,infinity) by         
C     the double exponential formula (DE formula)   
C               in double precision                 
C               By Masatake Mori                   
C     Before the first call of this routine         
C     HIAB must be called                           
C     ==== input ====                               
C     FUNC...name of function subprogram for        
C            integrand                              
C     A...lower bound of integration                
C     EPS...absolute error tolerance                
C     ==== output ====                              
C     V...integral of FUNC                          
C
      REAL*8 FUNC,A,EPS,V
      REAL*8 ONE,HALF
      REAL*8 H
      REAL*8 VOLD,VNEW,WM,WP
      REAL*8 EPS0,EPSQ,EPSM,EPSP
      REAL*8 AM,A0,AP,BM,B0,BP
C
      COMMON /COMDEF/ AM(640,3),A0(3),AP(640,3),
     $                BM(640,3),B0(3),BP(640,3) 
      COMMON /COMDEN/ NPOW,NEND
      COMMON /COMCNT/ NEVAL,L
C
      PARAMETER (ONE = 1)
      PARAMETER (HALF = ONE / 2)
C
      DATA EPS0 /1.0D-32/
      DATA EPSM /0/, EPSP /0/
C-added
      wm = 0.0d0
      wp = 0.0d0
C
      NEVAL = 0
      IF (ABS(EPS) .GE. EPS0) THEN
        EPSV = ABS(EPS)
      ELSE
        EPSV = EPS0 
      END IF
      EPSQ = 0.2 * SQRT(EPSV) 
      H = HALF
      IS = 2**NPOW
      IH = IS
      L = 1
  101 CONTINUE
      KM = 0
      KP = 0
      NM = 0
      NP = 0
      VNEW = 0
C
C     ---- initial step ---- 
C          integrate with mesh size = 0.5
C          and check decay of integrand
C
      DO 110 I = IS, NEND, IH 
        IF (KP .LE. 1) THEN
          WP = FUNC(AP(I,L) + A) * BP(I,L)
          NEVAL = NEVAL + 1
          VNEW = VNEW + WP
          IF (ABS(WP) .LE. EPSV) THEN
            KP = KP + 1
            IF (KP .GE. 2) THEN
              NP = I - IH
              GO TO 111
            END IF
          ELSE
            KP = 0
          END IF
        END IF
  110 CONTINUE
  111 CONTINUE
      IF (L .LE. 2) THEN
        IF (NP .EQ. 0) THEN
          L = L + 1 
          GO TO 101 
        END IF
      END IF
      DO 120 I = IS, NEND, IH 
        IF (KM .LE. 1) THEN
          WM = FUNC(AM(I,L) + A) * BM(I,L)
          NEVAL = NEVAL + 1
          VNEW = VNEW + WM
          IF (ABS(WM) .LE. EPSV) THEN
            KM = KM + 1
            IF (KM .GE. 2) THEN
              NM = I - IH
              GO TO 121
            END IF
          ELSE
            KM = 0
          END IF
        END IF
  120 CONTINUE
  121 CONTINUE
      VNEW = VNEW + FUNC(A0(L) + A) * B0(L)
      NEVAL = NEVAL + 1
      IF (NM .EQ. 0) THEN
        NM = NEND
        EPSM = 0.2 * SQRT(ABS(WM))
CW       WRITE (*,2001)
CW 2001  FORMAT (' (DEHINT) slow decay on negative side') 
      END IF
      IF (NP .EQ. 0) THEN
        NP = NEND
        EPSP = 0.2 * SQRT(ABS(WP))
CW      WRITE (*,2002)
CW 2002 FORMAT (' (DEHINT) slow decay on positive side') 
      END IF
      EPSQ = MAX(EPSQ, EPSM, EPSP)
C
C     ---- general step ---- 
C
      VOLD = H * VNEW
      DO 20 MSTEP = 1, NPOW
        VNEW = 0.0
        IH = IS
        IS = IS / 2 
        DO 540 I = IS, NM, IH 
          VNEW = VNEW
     $         + FUNC(AM(I,L) + A) * BM(I,L)
          NEVAL = NEVAL + 1
  540   CONTINUE
        DO 550 I = IS, NP, IH 
          VNEW = VNEW
     $         + FUNC(AP(I,L) + A) * BP(I,L)
          NEVAL = NEVAL + 1
  550   CONTINUE
        VNEW = (VOLD + H * VNEW) * HALF 
        IF (ABS(VNEW - VOLD) .LT. EPSQ) THEN
C
C        ---- converged and return ----
C
          V = VNEW
          RETURN
        END IF
        H = H * HALF
        VOLD = VNEW 
   20 CONTINUE
CW      WRITE (*,2003)
CW 2003 FORMAT (' (DEHINT) insufficient mesh',
CW   $        ' refinement')
      V = VNEW
      RETURN
      END 
C
C**********************************************************************
C
      SUBROUTINE HIAB()
C
C     Generate points and weights of the DE formula   
C     over half infinite interval (A,infinity)        
C     Before the first call of DEHINT this routine    
C     must be called                                  
C      L for DE half infinite transformation          
C         L = 1  x = exp(0.5*t - exp(-t))             
C         L = 2  x = exp(t - exp(-t))                 
C         L = 3  x = exp(2 * sinh t)                  
C
      REAL*8 ONE,HALF
      REAL*8 H,EH,EN,ENI
      REAL*8 SH,CH
      REAL*8 AM,A0,AP,BM,B0,BP
      COMMON /COMDEF/ AM(640,3),A0(3),AP(640,3),
     >       BM(640,3),B0(3),BP(640,3) 
      COMMON /COMDEN/ NPOW,NEND
      PARAMETER (ONE = 1)
      PARAMETER (HALF = ONE / 2)
      PARAMETER (N6 = 6)
      NPOW = N6
      NEND = 5 * 2**(NPOW+1)
      H = ONE / 2**(NPOW+1)
      EH = EXP(H)
C
C LE=1 ---- DE transformation x = exp(0.5*t-exp(-t)) ---- 
C
      A0(1) = EXP(-ONE)
      B0(1) = 1.5D0 * A0(1)
      EN = 1.0
      DO 10 I = 1, NEND
        EN = EN * EH
        ENI = 1 / EN
        SH = HALF * H * I
        AP(I,1) = EXP(SH - ENI)
        BP(I,1) = (HALF + ENI) * AP(I,1)
        AM(I,1) = EXP(-SH - EN)
        BM(I,1) = (HALF + EN) * AM(I,1) 
   10 CONTINUE
C
C LE=2 ---- DE transformation x = exp(t-exp(-t)) ----
C
      A0(2) = EXP(-ONE)
      B0(2) = 2.0 * A0(2)
      EN = 1.0
      DO 20 I = 1, NEND
        EN = EN * EH
        ENI = 1 / EN
        SH = H * I
        AP(I,2) = EXP(SH - ENI)
        BP(I,2) = (1 + ENI) * AP(I,2)
        AM(I,2) = EXP(-SH - EN)
        BM(I,2) = (1 + EN) * AM(I,2)
   20 CONTINUE
C
C---- LE=3 / DE transformation x = exp(2*sinh t) ----
C
      A0(3) = 1.0
      B0(3) = 2.0
      EN = 1.0
      DO 30 I = 1, NEND
        EN = EH * EN
        ENI = 1 / EN
        SH = EN - ENI
        CH = EN + ENI
        AP(I,3) = EXP(SH)
        BP(I,3) = CH * AP(I,3)
        AM(I,3) = 1 / AP(I,3) 
        BM(I,3) = CH * AM(I,3)
   30 CONTINUE
C
      RETURN
      END 
C
C**********************************************************************
C
      SUBROUTINE PNIG(c, x, n, alpha, beta, delta, mu) 
C
      DOUBLE PRECISION fdnig, c(n), x(n)
      DOUBLE PRECISION alpha, beta, delta, mu, a, v, eps
      DOUBLE PRECISION salpha, sbeta, sdelta, smu
      INTEGER i, n, neval, l
      EXTERNAL FDNIG
      COMMON /COMCNT/ NEVAL,L
      COMMON /S/ salpha, sbeta, sdelta, smu
     C
      EPS = 1.0D-12 
      salpha = alpha
      sbeta = beta
      sdelta = delta
      smu = mu
      do i=1, n
         CALL HIAB()
         A = x(i)
         CALL DEHINT (FDNIG, A, EPS, V)
         c(i) = 1.0d0 - V
      enddo
      RETURN
      END 
C
C**********************************************************************
C
      DOUBLE PRECISION FUNCTION FDNIG(X)
C
C     Function to integrate:
C
      DOUBLE PRECISION salpha, sbeta, sdelta, smu
      DOUBLE PRECISION x, c(1), xx(1)
      COMMON /S/ salpha, sbeta, sdelta, smu
C
      xx(1) = x
      CALL DNIG(c,xx,1,salpha,sbeta,sdelta,smu)
      FDNIG=c(1)
C
      RETURN
      END
C
C
C**********************************************************************
C
      SUBROUTINE IKNA(N,X,NM,BI,DI,BK,DK)
C
C       Purpose: Compute modified Bessel functions In(x) and
C                Kn(x), and their derivatives
C       Input:   x --- Argument of In(x) and Kn(x)
C                n --- Order of In(x) and Kn(x)
C       Output:  BI(n) --- In(x)
C                DI(n) --- In

'(x)C                BK(n) --- Kn(x)C                DK(n) --- Kn'(x)
C                NM --- Highest order computed
C       Routines called:
C            (1) IK01A for computing I0(x),I1(x),K0(x) & K1(x)
C            (2) MSTA1 and MSTA2 for computing the starting 
C                point for backward recurrence
C
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DIMENSION BI(0:N),DI(0:N),BK(0:N),DK(0:N)
C-dded
      f=0.0d0
C
        NM=N
        IF (X.LE.1.0D-100) THEN
           DO 10 K=0,N
              BI(K)=0.0D0
              DI(K)=0.0D0
              BK(K)=1.0D+300
10            DK(K)=-1.0D+300
           BI(0)=1.0D0
           DI(1)=0.5D0
           RETURN
        ENDIF
        CALL IK01A(X,BI0,DI0,BI1,DI1,BK0,DK0,BK1,DK1)
        BI(0)=BI0
        BI(1)=BI1
        BK(0)=BK0
        BK(1)=BK1
        DI(0)=DI0
        DI(1)=DI1
        DK(0)=DK0
        DK(1)=DK1
        IF (N.LE.1) RETURN
        IF (X.GT.40.0.AND.N.LT.INT(0.25*X)) THEN
           H0=BI0
           H1=BI1
           DO 15 K=2,N
           H=-2.0D0*(K-1.0D0)/X*H1+H0
           BI(K)=H
           H0=H1
15         H1=H
        ELSE
           M=MSTA1(X,200)
           IF (M.LT.N) THEN
              NM=M
           ELSE
              M=MSTA2(X,N,15)
           ENDIF
           F0=0.0D0
           F1=1.0D-100
           DO 20 K=M,0,-1
              F=2.0D0*(K+1.0D0)*F1/X+F0
              IF (K.LE.NM) BI(K)=F
              F0=F1
20            F1=F
           S0=BI0/F
           DO 25 K=0,NM
25            BI(K)=S0*BI(K)
        ENDIF
        G0=BK0
        G1=BK1
        DO 30 K=2,NM
           G=2.0D0*(K-1.0D0)/X*G1+G0
           BK(K)=G
           G0=G1
30         G1=G
        DO 40 K=2,NM
           DI(K)=BI(K-1)-K/X*BI(K)
40         DK(K)=-BK(K-1)-K/X*BK(K)
        RETURN
        END
C
C**********************************************************************
C
      SUBROUTINE IK01A(X,BI0,DI0,BI1,DI1,BK0,DK0,BK1,DK1)
C
C       Purpose: Compute modified Bessel functions I0(x), I1(1),
C                K0(x) and K1(x), and their derivatives
C       Input :  x   --- Argument
C       Output:  BI0 --- I0(x)
C                DI0 --- I0

'(x)C                BI1 --- I1(x)C                DI1 --- I1'(x)
C                BK0 --- K0(x)
C                DK0 --- K0

'(x)C                BK1 --- K1(x)C                DK1 --- K1'(x)
C
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DIMENSION A(12),B(12),A1(8)
        PI=3.141592653589793D0
        EL=0.5772156649015329D0
        X2=X*X
        IF (X.EQ.0.0D0) THEN
           BI0=1.0D0
           BI1=0.0D0
           BK0=1.0D+300
           BK1=1.0D+300
           DI0=0.0D0
           DI1=0.5D0
           DK0=-1.0D+300
           DK1=-1.0D+300
           RETURN
        ELSE IF (X.LE.18.0) THEN
           BI0=1.0D0
           R=1.0D0
           DO 15 K=1,50
              R=0.25D0*R*X2/(K*K)
              BI0=BI0+R
              IF (DABS(R/BI0).LT.1.0D-15) GO TO 20
15         CONTINUE
20         BI1=1.0D0
           R=1.0D0
           DO 25 K=1,50
              R=0.25D0*R*X2/(K*(K+1))
              BI1=BI1+R
              IF (DABS(R/BI1).LT.1.0D-15) GO TO 30
25         CONTINUE
30         BI1=0.5D0*X*BI1
        ELSE
           DATA A/0.125D0,7.03125D-2,
     &            7.32421875D-2,1.1215209960938D-1,
     &            2.2710800170898D-1,5.7250142097473D-1,
     &            1.7277275025845D0,6.0740420012735D0,
     &            2.4380529699556D01,1.1001714026925D02,
     &            5.5133589612202D02,3.0380905109224D03/
           DATA B/-0.375D0,-1.171875D-1,
     &            -1.025390625D-1,-1.4419555664063D-1,
     &            -2.7757644653320D-1,-6.7659258842468D-1,
     &            -1.9935317337513D0,-6.8839142681099D0,
     &            -2.7248827311269D01,-1.2159789187654D02,
     &            -6.0384407670507D02,-3.3022722944809D03/
           K0=12
           IF (X.GE.35.0) K0=9
           IF (X.GE.50.0) K0=7
           CA=DEXP(X)/DSQRT(2.0D0*PI*X)
           BI0=1.0D0
           XR=1.0D0/X
           DO 35 K=1,K0
35            BI0=BI0+A(K)*XR**K
           BI0=CA*BI0
           BI1=1.0D0
           DO 40 K=1,K0
40            BI1=BI1+B(K)*XR**K
           BI1=CA*BI1
        ENDIF
        IF (X.LE.9.0D0) THEN
           CT=-(DLOG(X/2.0D0)+EL)
           BK0=0.0D0
C-added
           ww = BK0
C
           W0=0.0D0
           R=1.0D0
           DO 65 K=1,50
              W0=W0+1.0D0/K
              R=0.25D0*R/(K*K)*X2
              BK0=BK0+R*(W0+CT)
              IF (DABS((BK0-WW)/BK0).LT.1.0D-15) GO TO 70
65            WW=BK0
70         BK0=BK0+CT
        ELSE
           DATA A1/0.125D0,0.2109375D0,
     &             1.0986328125D0,1.1775970458984D01,
     &             2.1461706161499D02,5.9511522710323D03,
     &             2.3347645606175D05,1.2312234987631D07/
           CB=0.5D0/X
           XR2=1.0D0/X2
           BK0=1.0D0
           DO 75 K=1,8
75            BK0=BK0+A1(K)*XR2**K
           BK0=CB*BK0/BI0
        ENDIF
        BK1=(1.0D0/X-BI1*BK0)/BI0
        DI0=BI1
        DI1=BI0-BI1/X
        DK0=-BK1
        DK1=-BK0-BK1/X
        RETURN
        END
C
C**********************************************************************
C
      INTEGER FUNCTION MSTA1(X,MP)
C
C       Purpose: Determine the starting point for backward  
C                recurrence such that the magnitude of    
C                Jn(x) at that point is about 10^(-MP)
C       Input :  x     --- Argument of Jn(x)
C                MP    --- Value of magnitude
C       Output:  MSTA1 --- Starting point   
C
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        A0=DABS(X)
        N0=INT(1.1*A0)+1
        F0=ENVJ(N0,A0)-MP
        N1=N0+5
        F1=ENVJ(N1,A0)-MP
        DO 10 IT=1,20             
           NN=N1-(N1-N0)/(1.0D0-F0/F1)                  
           F=ENVJ(NN,A0)-MP
           IF(ABS(NN-N1).LT.1) GO TO 20
           N0=N1
           F0=F1
           N1=NN
 10        F1=F
 20     MSTA1=NN
        RETURN
        END
C
C**********************************************************************
C
      INTEGER FUNCTION MSTA2(X,N,MP)
C
C       Purpose: Determine the starting point for backward
C                recurrence such that all Jn(x) has MP
C                significant digits
C       Input :  x  --- Argument of Jn(x)
C                n  --- Order of Jn(x)
C                MP --- Significant digit
C       Output:  MSTA2 --- Starting point
C
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        A0=DABS(X)
        HMP=0.5D0*MP
        EJN=ENVJ(N,A0)
        IF (EJN.LE.HMP) THEN
           OBJ=MP
           N0=INT(1.1*A0)
        ELSE
           OBJ=HMP+EJN
           N0=N
        ENDIF
        F0=ENVJ(N0,A0)-OBJ
        N1=N0+5
        F1=ENVJ(N1,A0)-OBJ
        DO 10 IT=1,20
           NN=N1-(N1-N0)/(1.0D0-F0/F1)
           F=ENVJ(NN,A0)-OBJ
           IF (ABS(NN-N1).LT.1) GO TO 20
           N0=N1
           F0=F1
           N1=NN
10         F1=F
20      MSTA2=NN+10
        RETURN
        END
C
C*******************************************************************************
C
      REAL*8 FUNCTION ENVJ(N,X)
C
      DOUBLE PRECISION X
      ENVJ=0.5D0*DLOG10(6.28D0*N)-N*DLOG10(1.36D0*X/N)
C
      RETURN
      END
C
C*******************************************************************************
C


Generated by  Doxygen 1.6.0   Back to index