tpc/generator/src.f

0001        SUBROUTINE DKIA(IFAC,X,A,DKI,DKID,IERRO)
0002 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0003 CCCC THE PURPOUSE OF THE ROUTINE IS THE CALCULATION OF THE
0004 CCCC MODIFIED BESSEL FUNCTIONS K_IA(X) AND K'_IA(X),
0005 CCCC WHERE I IS THE IMAGINARY UNIT AND X IS THE ARGUMENT OF THE
0006 CCCC FUNCTIONS. WE WILL REFER TO A AS THE ORDER OF THE FUNCTIONS.
0007 CCCC
0008 CCCC THE ROUTINE HAS THE OPTION OF COMPUTING SCALED FUNCTIONS.
0009 CCCC THIS SCALING CAN BE USED TO ENLARGE THE RANGE OF
0010 CCCC COMPUTATION.
0011 CCCC THE SCALED FUNCTIONS ARE DEFINED AS FOLLOWS
0012 CCCC     (S STANDS FOR SCALED  AND
0013 CCCC      L=SQRT{X**2-A**2} + A*ARCSIN(A/X)):
0014 CCCC
0015 CCCC                 EXP(L)*K_IA(X),            IF X >=ABS(A)
0016 CCCC    SK_IA(X)  =
0017 CCCC                 EXP(ABS(A)*PI/2)*K_IA(X),  IF X < ABS(A)
0018 CCCC
0019 CCCC
0020 CCCC                 EXP(L)*K'_IA(X),           IF X >=ABS(A)
0021 CCCC    SK'_IA(X) =
0022 CCCC                 EXP(ABS(A)*PI/2)*K'_IA(X), IF X < ABS(A)
0023 CCCC
0024 CCCC THE RANGE OF THE PARAMETERS (X,A) FOR THE COMPUTATION OF
0025 CCCC SCALED FUNCTIONS IS:
0026 CCCC        0 < X <= 1500,  -1500 <= A <= 1500.
0027 CCCC FOR FUNCTIONS WITHOUT SCALING, THE RANGE IS USUALLY LARGER
0028 CCCC THAN
0029 CCCC        0 < X <= 500,  -400 <= A <= 400,
0030 CCCC DEPENDING ON THE MACHINE OVERFLOW/UNDERFLOW PARAMETERS, WHICH
0031 CCCC ARE SET UP BY THE ROUTINE D1MACH.
0032 CCCC
0033 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0034 CCCC              METHODS OF COMPUTATION:
0035 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0036 CCCC    1) SERIES,  IF  ABS(A) > 0.044*ABS(X-3.1)**1.9+(X-3.1)
0037 CCCC    2) CONTINUED FRACTION,
0038 CCCC        IF X>3 AND ABS(A) < 380*(ABS((X-3)/2300))**0.572
0039 CCCC    3) AIRY-TYPE ASYMPTOTIC EXPANSION,
0040 CCCC        IF ABS(A) > 0.4*X+7.5 AND ABS(A) < 3.7*X-10
0041 CCCC    4) QUADRATURES,
0042 CCCC        IN THE REST OF CASES
0043 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0044 CCCC        DESCRIPTION OF INPUT/OUTPUT VARIABLES:
0045 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0046 CCCC  INPUTS:
0047 CCCC    X      : ARGUMENT OF THE FUNCTIONS
0048 CCC     A      : ORDER OF THE FUNCTIONS
0049 CCCC    IFAC   : INTEGER FLAG FOR THE SCALING
0050 CCCC      * IF IFAC=1, THE CODE COMPUTES KIA(X), KIA'(X)
0051 CCCC      * OTHERWISE, THE CODE COMPUTES SCALED KIA(X), KIA'(X)
0052 CCCC  OUTPUTS
0053 CCCC    DKI    :  KIA(X) FUNCTION
0054 CCCC    DKID   :  DERIVATIVE OF THE KIA(X) FUNCTION
0055 CCCC    IERRO: ERROR FLAG
0056 CCCC     * IF IERRO=0, COMPUTATION SUCCESSFUL.
0057 CCCC     * IF IERRO=1, COMPUTATION OUT OF RANGE.
0058 CCCC     * IF IERRO=2, VARIABLES X AND/OR A, OUT OF RANGE.
0059 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0060 CCCC        ACCURACY:
0061 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0062 CCCC THE RELATIVE ACCURACY IS:
0063 CCCC  BETTER THAN 10**(-13)   FOR (X,A) IN (0,200]X[-200,200];
0064 CCCC  BETTER THAN 5.10**(-13) FOR (X,A) IN (0,500]X[-500,500];
0065 CCCC  CLOSE TO    10**(-12)   FOR (X,A) IN (0,1500]X[-1500,1500].
0066 CCCC NEAR THE ZEROS OF THE FUNCTIONS (THERE ARE INFINITELY
0067 CCCC MANY OF THEM IN THE OSCILLATORY REGION) RELATIVE PRECISION
0068 CCCC LOOSES MEANING AND ONLY ABSOLUTE PRECISION MAKES SENSE.
0069 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0070 C     AUTHORS:
0071 C        AMPARO GIL    (U. CANTABRIA, SANTANDER, SPAIN).
0072 C                      E-MAIL: AMPARO.GIL@UNICAN.ES
0073 C        JAVIER SEGURA (U. CANTABRIA, SANTANDER, SPAIN).
0074 C                      E-MAIL: SEGURAJJ@UNICAN.ES
0075 C        NICO M. TEMME (CWI, AMSTERDAM, THE NETHERLANDS).
0076 C                      E-MAIL: NICO.TEMME@CWI.NL
0077 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0078 C    REFERENCES:
0079 C     THIS IS THE COMPANION SOFTWARE OF THE ARTICLES
0080 C      1)'COMPUTING SOLUTIONS OF THE MODIFIED BESSEL DIFFERENTIAL
0081 C         EQUATION FOR IMAGINARY ORDERS AND POSITIVE ARGUMENTS',
0082 C         A. GIL, J. SEGURA, N.M. TEMME
0083 C         ACM TRANS. MATH. SOFT. (2004)
0084 C      2)'MODIFIED BESSEL FUNCTIONS OF IMAGINARY ORDER AND
0085 C         POSITIVE ARGUMENT',
0086 C         A. GIL, J. SEGURA, N.M. TEMME
0087 C         ACM TRANS. MATH. SOFT. (2004)
0088 C    ADDITIONAL REFERENCES:
0089 C     - 'COMPUTATION OF THE MODIFIED BESSEL FUNCTION OF THE
0090 C         THIRD KIND FOR IMAGINARY ORDERS'
0091 C         A. GIL, J. SEGURA, N.M. TEMME
0092 C         J. COMPUT. PHYS. 175 (2002) 398-411
0093 C     - 'COMPUTATION OF THE MODIFIED BESSEL FUNCTIONS OF THE
0094 C         THIRD KIND OF IMAGINARY ORDERS:
0095 C         UNIFORM AIRY-TYPE ASYMPTOTIC EXPANSION'
0096 C         A. GIL, J. SEGURA, N.M. TEMME, PROCEEDINGS OPSFA 2001
0097 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0098 CCC        LOCAL VARIABLES:
0099 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0100 C      D:    X-3.1
0101 C      DF1:  0.044*ABS(D)**1.9+(X-3.1)
0102 C      DF2:  380*(ABS((X-3)/2300))**0.572
0103 C      DF3:  0.4*X+7.5
0104 C      DF4:  3.7*X-10
0105 C      PNU:  ORDER OF THE FUNCTION
0106 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0107        DOUBLE PRECISION A,D,DF1,DF2,DF3,DF4,DKI,
0108      + DKID,PNU,X
0109        INTEGER IERRO,IFAC
0110        IERRO=0
0111        PNU=A
0112        IF ((X.GT.1500.D0).OR.(X.LE.0.D0)) THEN
0113          IERRO=2
0114          DKI=0.D0
0115          DKID=0.D0
0116        ENDIF
0117        IF (ABS(PNU).GT.1500.D0) THEN
0118          IERRO=2
0119          DKI=0.D0
0120          DKID=0.D0
0121        ENDIF
0122        IF (IERRO.EQ.0) THEN
0123 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0124 CCC THESE FUNCTIONS ARE EVEN FUNCTIONS IN THE  C
0125 CCC PARAMETER A (=PNU)                         C
0126 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0127          IF (PNU.LT.0.D0) PNU=-PNU
0128          D=X-3.1D0
0129          DF1=0.044D0*ABS(D)**1.9D0+D
0130          DF2=380.D0*(ABS((X-3.D0)/2300.D0))**0.572D0
0131          DF3=0.4D0*X+7.5D0
0132          DF4=3.7D0*X-10.D0
0133          IF (PNU.GT.DF1) THEN
0134            CALL SERKIA(IFAC,X,PNU,DKI,DKID,IERRO)
0135          ELSEIF ((X.GT.3.D0).AND.(PNU.LT.DF2)) THEN
0136            CALL FRAKIA(IFAC,X,PNU,DKI,DKID,IERRO)
0137          ELSEIF ((PNU.GT.DF3).AND.(PNU.LT.DF4)) THEN
0138            CALL AIEXKI(IFAC,X,PNU,DKI,DKID,IERRO)
0139          ELSE
0140            CALL DKINT(IFAC,X,PNU,DKI,DKID,IERRO)
0141          ENDIF
0142        ENDIF
0143        RETURN
0144        END
0145
0146        SUBROUTINE DKINT(IFAC,XX,PNUA,DKINF,DKIND,IERRO)
0147 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0148 CCC  CALCULATION OF THE FUNCTIONS K,K' USING NON-OSCILLATING   C
0149 CCC  INTEGRAL REPRESENTATIONS                                  C
0150 CCC                                                            C
0151 CCC  INPUTS:                                                   C
0152 CCC      XX:   ARGUMENT OF THE FUNCTIONS                       C
0153 CCC      PNUA: ORDER OF THE FUNCTIONS                          C
0154 CCC      IFAC:                                                 C
0155 CCC               =1,  NON SCALED FUNCTIONS                    C
0156 CCC               OTYHERWISE,  SCALED FUNCTIONS                C
0157 CCC  OUTPUTS:                                                  C
0158 CCC      DKINF: K FUNCTION                                     C
0159 CCC      DKIND: DERIVATIVE OF THE K FUNCTION                   C
0160 CCC      IERRO: ERROR FLAG                                     C
0161 CCC      * IF IERRO=0, COMPUTATION SUCCESSFUL.                 C
0162 CCC      * IF IERRO=1, COMPUTATION OUT OF RANGE.               C
0163 CCC      * IF IERRO=2, ARGUMENT AND/OR ORDER,  OUT OF RANGE.   C
0164 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0165 CCC  METHOD OF COMPUTATION:
0166 CCC  * IF XX>=PNUA, THE NON-OSCILLATING INTEGRAL REPRESENTATIONS
0167 CCC                 FOR THE MONOTONIC REGION ARE USED
0168 CCC  * IF XX<PNUA,  THE NON-OSCILLATING INTEGRAL REPRESENTATIONS
0169 CCC                 FOR THE OSCILLATORY REGION ARE USED
0170 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0171 CCC    LOCAL VARIABLES:
0172 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0173 C CHI      : X*SINH(MU)-PNU*MU
0174 C CONTR1   : CONTRIBUTION OF THE SEMI-INFINITE INTEGRAL
0175 C            IN THE OSCILLATORY CASE (INCLUDING ADDITIONAL
0176 C            FACTORS: CONTR1=FOS1*FACTORS).
0177 C COSCHI   : COS(CHI)
0178 C COSH2M   : COSH(2*MU)
0179 C COSHM    : COSH(MU)
0180 C COSTH    : COS(THET), THET=ASIN(PNU/X)
0181 C DF1      : FACTOR (DEPENDING ON IFAC)
0182 C DMAIN    : PI*PNU*0.5
0183 C DMU      : SOLUTION MU OF COSH(MU)=PNU/X
0184 C DMU2     : 2*MU
0185 C DMU3     : 3*MU
0186 C DMU5     : 5*MU
0187 C DMU7     : 7*MU
0188 C DMU9     : 9*MU
0189 C DMUFAC   : MU*COSH(MU)-SINH(MU)
0190 C DMUTAN   : MU-TANH(MU)
0191 C FDOMIN   : X*(COS(THET)+THET*SIN(THET))
0192 C FOS1     : CONTRIBUTION OF THE SEMI-INFINITE INTEGRAL
0193 C            IN THE OSCILLATORY CASE (KIA(X)).
0194 C FOSD1    : CONTRIBUTION OF THE SEMI-INFINITE INTEGRAL
0195 C            IN THE OSCILLATORY CASE (KIA'(X)).
0196 C HIR      : MONOTONIC CASE, CONTRIBUTION OF THE INTEGRAL
0197 C            (KIA(X)).
0198 C HIRD     : MONOTONIC CASE, CONTRIBUTION OF THE INTEGRAL
0199 C            (KIA'(X)).
0200 C PI       : 3.1415...
0201 C PINU     : PI*PNU
0202 C PNU      : ORDER OF THE FUNCTION
0203 C SINCHI   : SIN(CHI)
0204 C SINH2M   : SINH(2*MU)
0205 C SINHM    : SINH(MU)
0206 C SINTH    : SIN(THET)
0207 C THET     : ASIN(PNU/X)
0208 C UNDER    : UNDERFLOW NUMBER
0209 C X        : ARGUMENT OF THE FUNCTION
0210 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0211        DOUBLE PRECISION CHI,CONTR1,COSCHI,COSH2M,COSHM,COSTH,
0212      + D1MACH,DF1,DKIND,DKINF,DMAIN,DMU,DMU2,DMU3,DMU5,DMU7,
0213      + DMU9,DMUFAC,DMUTAN,FDOMIN,FOS1,FOSD1,HIR,HIRD,PI,
0214      + PINU,PNU,PNUA,SINCHI,SINH2M,SINHM,SINTH,THET,UNDER,
0215      + X,XX
0216        INTEGER IERRO,IFAC
0217        COMMON/ARGU/X,PNU
0218        COMMON/PARMON/THET,SINTH,COSTH,FDOMIN
0219        COMMON/PAROS1/DMU,COSHM,SINHM,DMUFAC,COSCHI,SINCHI
0220        COMMON/PAROS2/COSH2M,SINH2M,DMAIN
0221        COMMON/PAROS3/DMUTAN
0222 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
0223        UNDER=D1MACH(1)*1.D+8
0224        X=XX
0225        PNU=PNUA
0226        IERRO=0
0227        IF (X.GE.PNU) THEN
0228 CCCCCCCCCCCCCCCCCCCCCCCCCC
0229 CCC  MONOTONIC REGION   CC
0230 CCCCCCCCCCCCCCCCCCCCCCCCCC
0231          SINTH=PNU/X
0232          THET=ASIN(SINTH)
0233          COSTH=COS(THET)
0234          FDOMIN=X*(COSTH+THET*SINTH)
0235          IF (IFAC.EQ.1) THEN
0236            IF (-FDOMIN.LE.LOG(UNDER)) IERRO=1
0237          ENDIF
0238          IF (IERRO.EQ.0) THEN
0239 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0240 CCC CALCULATION OF KIA:
0241 CCC    CALL TO THE TRAPEZOIDAL ROUTINE  C
0242 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0243            CALL TRAPRE(1,HIR)
0244            IF (IFAC.EQ.1) THEN
0245              DKINF=0.5D0*HIR*EXP(-FDOMIN)
0246            ELSE
0247              DKINF=0.5D0*HIR
0248            ENDIF
0249 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0250 CCC CALCULATION OF KIA':
0251 CCC    CALL TO THE TRAPEZOIDAL ROUTINE  C
0252 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0253            CALL TRAPRE(10,HIRD)
0254            IF (IFAC.EQ.1) THEN
0255              DKIND=-HIRD*EXP(-FDOMIN)
0256            ELSE
0257              DKIND=-HIRD
0258            ENDIF
0259          ELSE
0260            DKINF=0.D0
0261            DKIND=0.D0
0262          ENDIF
0263        ELSE
0264 CCCCCCCCCCCCCCCCCCCCCCCCCCCC
0265 CCC  OSCILLATORY REGION   CC
0266 CCCCCCCCCCCCCCCCCCCCCCCCCCCC
0267          PI=ACOS(-1.D0)
0268          IF (IFAC.EQ.1) THEN
0269            IF ((-PI*PNU*0.5D0).LE.LOG(UNDER)) IERRO=1
0270          ENDIF
0271          IF (IERRO.EQ.0) THEN
0272            COSHM=PNU/X
0273            DMU=LOG(COSHM+SQRT((COSHM-1.D0)*(COSHM+1.D0)))
0274            DMU2=2.D0*DMU
0275            COSH2M=COSH(DMU2)
0276            SINHM=SINH(DMU)
0277            SINH2M=SINH(DMU2)
0278            IF (DMU.GT.0.1D0) THEN
0279              CHI=X*SINHM-PNU*DMU
0280              DMUFAC=DMU*COSHM-SINHM
0281            ELSE
0282              DMU2=DMU*DMU
0283              DMU3=DMU2*DMU
0284              DMU5=DMU3*DMU2
0285              DMU7=DMU5*DMU2
0286              DMU9=DMU7*DMU2
0287              CHI=-2.D0*X*(1.D0/6.D0*DMU3+1.D0/60.D0*DMU5+
0288      +        3.D0/5040.D0*DMU7+4.D0/362880.D0*DMU9)
0289              DMUFAC=DMU3/3.D0+DMU5/30.D0+DMU7/840.D0+DMU9/45360.D0
0290            ENDIF
0291            DMUTAN=DMU-TANH(DMU)
0292            COSCHI=COS(CHI)
0293            SINCHI=SIN(CHI)
0294            PINU=PI*PNU
0295            DMAIN=PINU*0.5D0
0296 CCCCCCCCCCCCCCCCCC
0297 CCC  INTEGRALS  CC
0298 CCCCCCCCCCCCCCCCCC
0299            CALL TRAPRE(2,FOS1)
0300 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0301 CCC THEN, THE KIA FUNCTION IS GIVEN BY ... CC
0302 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0303            IF (IFAC.EQ.1) THEN
0304              DF1=EXP(-DMAIN)
0305              CONTR1=DF1*FOS1
0306            ELSE
0307              DF1=1.D0
0308              CONTR1=DF1*FOS1
0309            ENDIF
0310            DKINF=CONTR1
0311 CCCCCCCCCCCCCCCCCCCCCCCCCC
0312 CCC CALCULATION OF KIA' CC
0313 CCCCCCCCCCCCCCCCCCCCCCCCCC
0314 CCCCCCCCCCCCCCCCCC
0315 CCC  INTEGRALS  CC
0316 CCCCCCCCCCCCCCCCCC
0317            CALL TRAPRE(20,FOSD1)
0318 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0319 CCC THEN, KIA' IS GIVEN BY ...
0320 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0321            IF (IFAC.EQ.1) THEN
0322              DF1=EXP(-DMAIN)
0323              CONTR1=DF1*FOSD1
0324            ELSE
0325              DF1=1.D0
0326              CONTR1=DF1*FOSD1
0327            ENDIF
0328            DKIND=CONTR1
0329          ELSE
0330            DKINF=0.D0
0331            DKIND=0.D0
0332          ENDIF
0333        ENDIF
0334        RETURN
0335        END
0336
0337        DOUBLE PRECISION FUNCTION FA(U)
0338 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0339 CCC  COMPUTATION OF THE INTEGRAND FOR THE
0340 CCC  INTEGRAL REPRESENTATION OF THE K FUNCTION
0341 CCC  IN THE MONOTONIC REGION.
0342 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0343 CCC        INPUT VARIABLE:
0344 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0345 C    U       : ARGUMENT OF THE FUNCTION
0346 C             (VARIABLE OF INTEGRATION).
0347 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0348 CCC        LOCAL VARIABLES:
0349 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0350 C    PHIR    : X*COSH(U)*COS(V(U))+PNU*V(U)
0351 C              -FDOMIN, WHERE
0352 C              V(U)=ASIN(U/SINH(U)*SIN(THET))
0353 C              AND FDOMIN=X*(COS(THET)+
0354 C                         THET*SIN(THET))
0355 C    UNDER   : UNDERFLOW NUMBER
0356 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0357        DOUBLE PRECISION D1MACH,PHIR,U,UNDER
0358 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
0359        UNDER=D1MACH(1)*1.D+8
0360        IF ((-PHIR(U)).LE.LOG(UNDER)) THEN
0361          FA=0.D0
0362        ELSE
0363          FA=EXP(-PHIR(U))
0364        ENDIF
0365        RETURN
0366        END
0367
0368        DOUBLE PRECISION FUNCTION FAD(T)
0369 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0370 CCC  COMPUTATION OF THE INTEGRAND FOR THE
0371 CCC  INTEGRAL REPRESENTATION OF THE DERIVATIVE OF
0372 CCC  THE K FUNCTION IN THE MONOTONIC REGION.
0373 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0374 CCC        INPUT VARIABLE:
0375 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0376 C    T       : ARGUMENT OF THE FUNCTION
0377 C             (VARIABLE OF INTEGRATION).
0378 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0379 CCC        LOCAL VARIABLES:
0380 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0381 C ANG      :  THET-V(U)
0382 C ANGH     :  ANG/2
0383 C COSTH    :  COS(THET)
0384 C COSVU    :  COS(V(U))
0385 C DFAC     :  COS(THET)+(FU1+2.D0*SIN(ANG/2)
0386 C             *SIN(ANG/2))/COS(VU)
0387 C DJACO    :  COSH(T)*EXP(S)/(1.D0+EXP(S))
0388 C EXPO     :  EXP(-PHIR(U))
0389 C FDOMIN   :  X*(COS(THET)+THET*SIN(THET))
0390 C FU1      :  2*SINH(U/2)**2
0391 C FUAC     :  U**3/6+U**5/120+U**7/5040
0392 C PHIR(U)  :  X*COSH(U)*COS(V(U))+PNU*V(U)-FDOMIN
0393 C S        :  SINH(T)
0394 C SINANH   :  SIN(ANG/2)
0395 C SINHU    :  SINH(U)
0396 C SINHUH   :  SINH(U/2)
0397 C SINTH    :  SIN(THET)
0398 C THET     :  ASIN(PNU/X)
0399 C U        :  LOG(1+EXP(S))
0400 C U2       :  U**2
0401 C U3       :  U**3
0402 C U5       :  U**5
0403 C U7       :  U**7
0404 C UH       :  U/2
0405 C UNDER    :  UNDERFLOW NUMBER
0406 C V(U)     :  ASIN(U/SINH(U)*SIN(THET))
0407 C VU       :  V(U)
0408 C Y        :  EXP(S)
0409 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0410        DOUBLE PRECISION ANG,ANGH,COSTH,COSVU,D1MACH,
0411      + DFAC,DJACO,EXPO,FDOMIN,FU1,FUAC,PHIR,S,SINANH,
0412      + SINHU,SINHUH,SINTH,T,THET,U,U2,U3,U5,U7,UH,UNDER,
0413      + V,VU,Y
0414        COMMON/PARMON/THET,SINTH,COSTH,FDOMIN
0415 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
0416        UNDER=D1MACH(1)*1.D+8
0417        S=SINH(T)
0418        Y=EXP(S)
0419        U=LOG(1.D0+Y)
0420        DJACO=COSH(T)*Y/(1.D0+Y)
0421        IF (ABS(U).LT.1.D-1) THEN
0422          IF (ABS(U).LT.UNDER) THEN
0423            DFAC=COSTH
0424          ELSE
0425            UH=U*0.5D0
0426            SINHUH=SINH(UH)
0427            FU1=2.D0*SINHUH*SINHUH
0428            U2=U*U
0429            U3=U2*U
0430            U5=U2*U3
0431            U7=U5*U2
0432            FUAC=U3/6.D0+U5/120.D0+U7/5040.D0
0433            SINHU=SINH(U)
0434            VU=V(U)
0435            COSVU=COS(VU)
0436            ANG=ASIN(-SINTH/(COSTH*U/SINHU+COSVU)*
0437      +        (SINHU+U)*FUAC/(SINHU*SINHU))
0438            ANGH=0.5D0*ANG
0439            SINANH=SIN(ANGH)
0440            DFAC=COSTH+(FU1+2.D0*SINANH*SINANH)/COSVU
0441          ENDIF
0442        ELSE
0443          VU=V(U)
0444          ANG=THET-VU
0445          FU1=COSH(U)-1.D0
0446          ANGH=0.5D0*ANG
0447          SINANH=SIN(ANGH)
0448          DFAC=COSTH+(FU1+2.D0*SINANH*SINANH)/COS(VU)
0449        ENDIF
0450        IF ((-PHIR(U)).LE.LOG(UNDER)) THEN
0451          EXPO=0.D0
0452        ELSE
0453          EXPO=EXP(-PHIR(U))
0454        ENDIF
0455        FAD=EXPO*DFAC*DJACO
0456        RETURN
0457        END
0458
0459        DOUBLE PRECISION FUNCTION FDTAU2(T)
0460 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0461 CCC  COMPUTATION OF THE INTEGRAND FOR THE
0462 CCC  INTEGRAL FOR THE DERIVATIVE OF THE K FUNCTION
0463 CCC  IN EQ.(37) OF THE REFERENCE:
0464 CCC  'COMPUTING SOLUTIONS OF THE MODIFIED BESSEL
0465 CCC   DIFFERENTIAL EQUATION FOR IMAGINARY ORDERS
0466 CCC   AND POSITIVE ARGUMENTS',
0467 CCC   A. GIL, J. SEGURA, N.M. TEMME
0468 CCC   ACM TRANS. MATH. SOFT. (2004)
0469 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0470 CCC        INPUT VARIABLE:
0471 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0472 C    T    : ARGUMENT OF THE FUNCTION
0473 C           (VARIABLE OF INTEGRATION).
0474 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0475 CCC         LOCAL VARIABLES:
0476 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0477 C ARGU    : (COSH(MU)*U-DMUFAC)/SINH(U)
0478 C ARGU2   : ARGU**2
0479 C COSCHI  : COSH(CHI), CHI=X*SINH(MU)-PNU*MU
0480 C COSH2M  : COSH(2*MU)
0481 C COSHM   : COSH(MU)
0482 C COSHU   : COSH(U)
0483 C COSS    : COS(SIGMA(U))
0484 C D1      : COSH(MU)-DMUFAC
0485 C DELTA   : -SIN(CHI)*COS(SIGMA(U))*COSH(U)
0486 C           +COS(CHI)*SIN(SIGMA(U))*SINH(U)
0487 C DERI    : (COSH(MU)/SINH(U)-D1*COSH(U)/SINH(U)**2)
0488 C           /SQRT(1-ARGU2)
0489 C DJACO   : COSH(T)*EXP(S)/SQRT(1+EXP(S)**2)
0490 C DMAIN   : PI*PNU*0.5
0491 C DMU     : SOLUTION MU OF COSH(MU)=PNU/X
0492 C DMU2    : 2*MU
0493 C DMUFAC  : MU*COSH(MU)-SINH(MU)
0494 C EXPON   : EXP(-(PHIB(U)-PI*PNU*0.5))
0495 C F1      : SINH(U-MU)/(U-MU)
0496 C G1      : Z/6+Z3/120+Z5/5040+Z7/362880, Z=2*Y
0497 C GAMMA   : SIN(CHI)*SIN(SIGMA(U))*SINH(U)+
0498 C           COS(CHI)*COS(SIGMA(U))*COSH(U)
0499 C PHIB(U) : X*COSH(U)*COS(SIGMA(U))+PNU*SIGMA(U),
0500 C           WHERE X=ARGUMENT OF THE FUNCTION,
0501 C           SIGMA(U)=ARCSIN((COSH(MU)*U-DMUFAC)/SINH(U)
0502 C RESTO   : -SIN(CHI)*COS(SIGMA(U))*COSH(U)
0503 C           +COS(CHI)*SIN(SIGMA(U))*SINH(U)+
0504 C           (SIN(CHI)*SIN(SIGMA(U))*SINH(U)+
0505 C            COS(CHI)*COS(SIGMA(U))*COSH(U))*DERI
0506 C S       : SINH(T)
0507 C SIGMA(U): ASIN((COSH(MU)*U-DMUFAC)/SINH(U))
0508 C SIGMAU  : SIGMA(U)
0509 C SINCHI  : SIN(CHI)
0510 C SINH2M  : SINH(2*MU)
0511 C SINHM   : SINH(MU)
0512 C SINHU   : SINH(U)
0513 C SINHU2  : 2*SINH(U)
0514 C SINS    : SIN(SIGMA(U))
0515 C U       : MU+LOG(X+SQRT(X**2+1))
0516 C UNDER   : UNDERFLOW NUMBER
0517 C X       : EXP(S)
0518 C Y       : U-MU
0519 C Z       : 2*Y
0520 C Z2      : Z**2
0521 C Z3      : Z**3
0522 C Z5      : Z**5
0523 C Z7      : Z**7
0524 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0525        DOUBLE PRECISION ARGU,ARGU2,COSCHI,COSH2M,COSHM,
0526      + COSHU,COSS,D1,D1MACH,DELTA,DERI,DJACO,DMAIN,DMU,
0527      + DMU2,DMUFAC,DMUTAN,EXPON,F1,G1,GAMMA,PHIB,RESTO,
0528      + S,SIGMA,SIGMAU,SINCHI,SINH2M,SINHM,SINHU,SINHU2,
0529      + SINS,T,U,UNDER,X,Y,Z,Z2,Z3,Z5,Z7
0530        COMMON/PAROS1/DMU,COSHM,SINHM,DMUFAC,COSCHI,SINCHI
0531        COMMON/PAROS2/COSH2M,SINH2M,DMAIN
0532        COMMON/PAROS3/DMUTAN
0533 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
0534        UNDER=D1MACH(1)*1.D+8
0535        S=SINH(T)
0536        X=EXP(S)
0537        U=DMUTAN+LOG(X+SQRT(X*X+1.D0))
0538        Y=U-DMU
0539        COSHU=COSH(U)
0540        SINHU=SINH(U)
0541        IF (ABS(Y).LE.1.D-1) THEN
0542          IF (ABS(Y).GT.UNDER) THEN
0543            F1=SINH(Y)/Y
0544          ELSE
0545            F1=1.D0
0546          ENDIF
0547          Z=Y*2.D0
0548          Z2=Z*Z
0549          Z3=Z2*Z
0550          Z5=Z3*Z2
0551          Z7=Z5*Z2
0552          G1=Z/6.D0+Z3/120.D0+Z5/5040.D0+Z7/362880.D0
0553          DMU2=2.D0*DMU
0554          DERI=SINHU/(F1-COSHM*COSHU)*SQRT(COSH(DMU2)*F1*F1+
0555      +   2.D0*SINH(DMU2)*G1-COSHM*COSHM)
0556          DERI=1.D0/DERI
0557        ELSE
0558          D1=COSHM*U-DMUFAC
0559          ARGU=D1/SINHU
0560          ARGU2=ARGU*ARGU
0561          SINHU2=SINHU*SINHU
0562          DERI=1.D0/SQRT(1.D0-ARGU2)*(COSHM/SINHU-
0563      +   D1*COSHU/SINHU2)
0564          IF (U.LT.DMU) DERI=-DERI
0565        ENDIF
0566        DJACO=COSH(T)*X/SQRT(1.D0+X*X)
0567        IF ((-(PHIB(U)-DMAIN)).LE.LOG(UNDER)) THEN
0568          EXPON=0.D0
0569          FDTAU2=0.D0
0570        ELSE
0571          EXPON=EXP(-(PHIB(U)-DMAIN))
0572          SIGMAU=SIGMA(U)
0573          SINS=SIN(SIGMAU)
0574          COSS=COS(SIGMAU)
0575          GAMMA=COSCHI*SINS*SINHU-SINCHI*COSS*COSHU
0576          DELTA=-COSCHI*COSS*COSHU-SINCHI*SINS*SINHU
0577          RESTO=DELTA+GAMMA*DERI
0578          FDTAU2=EXPON*RESTO*DJACO
0579        ENDIF
0580        RETURN
0581        END
0582
0583        DOUBLE PRECISION FUNCTION FSTAU2(T)
0584 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0585 CCC  COMPUTATION OF THE INTEGRAND FOR THE
0586 CCC  INTEGRAL FOR THE K FUNCTION IN EQ.(37) OF
0587 CCC  THE REFERENCE:
0588 CCC  'COMPUTING SOLUTIONS OF THE MODIFIED BESSEL
0589 CCC   DIFFERENTIAL EQUATION FOR IMAGINARY ORDERS
0590 CCC   AND POSITIVE ARGUMENTS',
0591 CCC   A. GIL, J. SEGURA, N.M. TEMME
0592 CCC   ACM TRANS. MATH. SOFT. (2004)
0593 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0594 CCC        INPUT VARIABLE:
0595 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0596 C    T    : ARGUMENT OF THE FUNCTION
0597 C           (VARIABLE OF INTEGRATION).
0598 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0599 CCC        LOCAL VARIABLES:
0600 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0601 C ARGU    : (COSH(MU)*U-DMUFAC)/SINH(U)
0602 C ARGU2   : ARGU**2
0603 C COSCHI  : COSH(CHI), CHI=X*SINH(MU)-PNU*MU
0604 C COSH2M  : COSH(2*MU)
0605 C COSHM   : COSH(MU)
0606 C D1      : COSH(MU)-DMUFAC
0607 C DERI    : (COSH(MU)/SINH(U)-D1*COSH(U)/SINH(U)**2)
0608 C           /SQRT(1-ARGU2)
0609 C DJACO   : COSH(T)*EXP(S)/SQRT(1+EXP(S)**2)
0610 C DMAIN   : PI*PNU*0.5
0611 C DMU     : SOLUTION MU OF COSH(MU)=PNU/X
0612 C DMU2    : 2*MU
0613 C DMUFAC  : MU*COSH(MU)-SINH(MU)
0614 C EXPON   : EXP(-(PHIB(U)-PI*PNU*0.5))
0615 C F1      : SINH(U-MU)/(U-MU)
0616 C G1      : Z/6+Z3/120+Z5/5040+Z7/362880, Z=2*Y
0617 C PHIB(U) : X*COSH(U)*COS(SIGMA(U))+PNU*SIGMA(U),
0618 C           WHERE X=ARGUMENT OF THE FUNCTION,
0619 C           SIGMA(U)=ARCSIN((COSH(MU)*U-DMUFAC)/SINH(U)
0620 C RESTO   : COS(CHI)+SIN(CHI)*DERI
0621 C S       : SINH(T)
0622 C SINCHI  : SIN(CHI)
0623 C SINH2M  : SINH(2*MU)
0624 C SINHM   : SINH(MU)
0625 C SINHU   : SINH(U)
0626 C SINHU2  : 2*SINH(U)
0627 C U       : MU+LOG(X+SQRT(X**2+1))
0628 C UNDER   : UNDERFLOW NUMBER
0629 C X       : EXP(S)
0630 C Y       : U-MU
0631 C Z       : 2*Y
0632 C Z2      : Z**2
0633 C Z3      : Z**3
0634 C Z5      : Z**5
0635 C Z7      : Z**7
0636 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0637        DOUBLE PRECISION ARGU,ARGU2,COSCHI,COSH2M,
0638      + COSHM,D1,D1MACH,DERI,DJACO,DMAIN,DMU,DMU2,
0639      + DMUFAC,DMUTAN,EXPON,F1,G1,PHIB,RESTO,S,SINCHI,
0640      + SINH2M,SINHM,SINHU,SINHU2,T,U,UNDER,X,Y,Z,Z2,
0641      + Z3,Z5,Z7
0642        COMMON/PAROS1/DMU,COSHM,SINHM,DMUFAC,COSCHI,SINCHI
0643        COMMON/PAROS2/COSH2M,SINH2M,DMAIN
0644        COMMON/PAROS3/DMUTAN
0645 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
0646        UNDER=D1MACH(1)*1.D+8
0647        S=SINH(T)
0648        X=EXP(S)
0649        U=DMUTAN+LOG(X+SQRT(X*X+1.D0))
0650        Y=U-DMU
0651        IF (ABS(Y).GT.UNDER) THEN
0652          F1=SINH(Y)/Y
0653        ELSE
0654          F1=1.D0
0655        ENDIF
0656        Z=Y*2.D0
0657        Z2=Z*Z
0658        IF (ABS(Y).LE.0.1D0) THEN
0659          Z3=Z2*Z
0660          Z5=Z3*Z2
0661          Z7=Z5*Z2
0662          G1=Z/6.D0+Z3/120.D0+Z5/5040.D0+Z7/362880.D0
0663          DMU2=2.D0*DMU
0664          DERI=SINH(U)/(F1-COSHM*COSH(U))*
0665      +      SQRT(COSH(DMU2)*F1*F1+2.D0*SINH(DMU2)*G1-
0666      +      COSHM*COSHM)
0667          DERI=1.D0/DERI
0668        ELSE
0669          D1=COSHM*U-DMUFAC
0670          SINHU=SINH(U)
0671          ARGU=D1/SINHU
0672          ARGU2=ARGU*ARGU
0673          SINHU2=SINHU*SINHU
0674          DERI=1.D0/SQRT(1.D0-ARGU2)*
0675      +    (COSHM/SINHU-D1*COSH(U)/SINHU2)
0676          IF (U.LT.DMU) DERI=-DERI
0677        ENDIF
0678        DJACO=COSH(T)*X/SQRT(1.D0+X*X)
0679        IF ((-(PHIB(U)-DMAIN)).LE.LOG(UNDER)) THEN
0680          EXPON=0.D0
0681          FSTAU2=0.D0
0682        ELSE
0683          EXPON=EXP(-(PHIB(U)-DMAIN))
0684          RESTO=(COSCHI+SINCHI*DERI)
0685          FSTAU2=EXPON*RESTO*DJACO
0686        ENDIF
0687        RETURN
0688        END
0689
0690        SUBROUTINE TRAPRE(IC,TI)
0691 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0692 CCC IMPLEMENTATION OF AN ADAPTIVE TRAPEZOIDAL RULE
0693 CCC FOR COMPUTING THE INTEGRAL REPRESENTATIONS OF
0694 CCC THE FUNCTIONS, FOR THE DIFFERENT REGIONS
0695 CCC (MONOTONIC OR OSCILLATORY).
0696 CCC
0697 CCC   INPUT VARIABLE:
0698 CCC      IC: DEPENDING ON THE VALUES OF IC,
0699 CCC          DIFFERENT INTEGRALS ARE COMPUTED:
0700 CCC          *IC=1, K FUNCTION, MONOTONIC REGION
0701 CCC          *IC=2, K FUNCTION, OSCILLATORY REGION
0702 CCC          *IC=10, DERIVATIVE OF THE K FUNCTION,
0703 CCC                  MONOTONIC REGION.
0704 CCC          *IC=20, DERIVATIVE OF THE K FUNCTION,
0705 CCC                  OSCILLATORY REGION.
0706 CCC   OUTPUT VARIABLE:
0707 CCC          TI,  RESULT OF THE INTEGRAL
0708 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0709 CCC            LOCAL VARIABLES
0710 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0711 C A        :  LOWER INTEGRATION LIMIT
0712 C B        :  UPPER INTEGRATION LIMIT
0713 C DELTA    :  CALCULATES THE RELATIVE PRECISION
0714 C EPS      :  RELATIVE PRECISION PARAMETER USED IN
0715 C             THE CALCULATION
0716 C H        :  INTEGRATION STEP
0717 C PNU      :  ORDER OF THE FUNCTION
0718 C SUM      :  ACCUMULATES THE ELEMENTARY
0719 C             CONTRIBUTIONS
0720 C TIN      :  EVALUATED INTEGRAL
0721 C X        :  ARGUMENT
0722 C XAC      :  INTEGRATION ABCISSA
0723 C Z        :  X/PNU
0724 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0725        DOUBLE PRECISION A,B,D1MACH,DELTA,EPS,FINTE,H,
0726      + PNU,SUM,TI,TIN,X,XAC,Z
0727        INTEGER I,IC,IFI,N
0728        COMMON/ARGU/X,PNU
0729        EPS=D1MACH(3)
0730        IF (EPS.LT.1.D-14) EPS=1.D-14
0731        N=0
0732 CCCCC   INTEGRATION LIMITS: A,B
0733        Z=X/PNU
0734        IF ((Z.GT.0.999D0).AND.(Z.LT.1.001D0)) THEN
0735          IF ((IC.EQ.1).OR.(IC.EQ.10)) THEN
0736            A=-3.5D0
0737          ELSE
0738            A=-4.5D0
0739          ENDIF
0740        ELSE
0741          IF ((IC.EQ.2).OR.(IC.EQ.20)) THEN
0742            A=-2.5D0
0743          ELSE
0744            A=-4.5D0
0745          ENDIF
0746        ENDIF
0747        B=-A
0748        H=B-A
0749        TI=0.5D0*H*(FINTE(IC,A)+FINTE(IC,B))
0750        DELTA=1.D0+EPS
0751 11     N=N+1
0752        H=0.5D0*H
0753        IF (N.EQ.1) THEN
0754          IFI=1
0755        ELSE
0756          IFI=2*IFI
0757        ENDIF
0758        SUM=0.D0
0759        DO 3 I=1,IFI
0760          XAC=A+DBLE(2*I-1)*H
0761          SUM=SUM+FINTE(IC,XAC)
0762  3     CONTINUE
0763        TIN=0.5D0*TI+H*SUM
0764        IF ((TIN.NE.0.D0).AND.(N.GT.4)) THEN
0765          DELTA=ABS(1.D0-TI/TIN)
0766        ENDIF
0767        TI=TIN
0768        IF ((DELTA.GT.EPS).AND.(N.LT.9)) GOTO 11
0769        RETURN
0770        END
0771
0772        DOUBLE PRECISION FUNCTION FINTE(IC,T)
0773 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0774 CCC  AUXILIARY FUNCTION FOR THE SUBROUTINE TRAPRE.
0775 CCC  THIS FUNCTION CALLS THE DIFFERENT FUNCTIONS
0776 CCC  CONTAINING THE INTEGRANDS WHICH ARE INTEGRATED
0777 CCC  BY TRAPRE.
0778 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0779 CCC   INPUT:                                      C
0780 CCC      IC: INTEGER PARAMETER. ITS ADMISSIBLE    C
0781 CCC          VALUES ARE THE SAME AS IN THE        C
0782 CCC          SUBROUTINE TRAPRE.                   C
0783 CCC      T:  INTEGRATION ABSCISSA                 C
0784 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0785 CCC       LOCAL VARIABLES:
0786 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0787 C  FA      : MONOTONIC PART CONTRIBUTION
0788 C            (KIA(X) FUNCTION).
0789 C  FAD     : MONOTONIC PART CONTRIBUTION
0790 C            (KIA'(X) FUNCTION).
0791 C  FDTAU2  : OSCILLATORY PART: TAU CONTRIBUTION ROUTINE
0792 C            (KIA'(X) FUNCTION). SEMI-INFINITE INTEGRAL.
0793 C                USED FOR LARGE PNU
0794 C  FSTAU2  : OSCILLATORY PART: TAU CONTRIBUTION ROUTINE
0795 C            (KIA(X) FUNCTION). SEMI-INFINITE INTEGRAL.
0796 C                  USED FOR LARGE PNU
0797 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0798        DOUBLE PRECISION FA,FAD,FDTAU2,FSTAU2,T
0799        INTEGER IC
0800        IF (IC.EQ.1) THEN
0801          FINTE=FA(T)
0802        ENDIF
0803        IF (IC.EQ.2) THEN
0804          FINTE=FSTAU2(T)
0805        ENDIF
0806        IF (IC.EQ.10) THEN
0807          FINTE=FAD(T)
0808        ENDIF
0809        IF (IC.EQ.20) THEN
0810          FINTE=FDTAU2(T)
0811        ENDIF
0812        END
0813
0814        SUBROUTINE FRAKIA(IFAC,X,PNU,PSER,PSERD,IERRO)
0815 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0816 CCC   IMPLEMENTATION OF THE CONTINUED FRACTION     C
0817 CCC   METHOD FOR THE CALCULATION OF K AND K'       C
0818 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0819 CCC   INPUT VARIABLES:                             C
0820 CCC      X:   ARGUMENT OF THE FUNCTIONS            C
0821 CCC      PNU: ORDER OF THE FUNCTIONS               C
0822 CCC      IFAC:                                     C
0823 CCC               =1,  NON SCALED FUNCTIONS        C
0824 CCC               OTHERWISE,  SCALED FUNCTIONS     C
0825 CCC   OUTPUT VARIABLES:                            C
0826 CCC      PSER: K FUNCTION                          C
0827 CCC      PSERD: DERIVATIVE OF THE K FUNCTION       C
0828 CCC      IERRO: ERROR FLAG                         C
0829 CCC      * IF IERRO=0, COMPUTATION SUCCESSFUL.     C
0830 CCC      * IF IERRO=1, COMPUTATION OUT OF RANGE.   C
0831 CCC      * IF IERRO=2, ARGUMENT AND/OR ORDER,      C
0832 CCC                    OUT OF RANGE.               C
0833 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0834 CCC           LOCAL VARIABLES:                     C
0835 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0836 C A      : (PNU), ORDER OF THE FUNCTIONS
0837 C A2     : A**2
0838 C AA     : AUXILIARY VARIABLE IN THE IMPLEMENTATION
0839 C          OF THE LENZ-THOMPSON ALGORITHM
0840 C AB     : AUXILIARY VARIABLE IN THE IMPLEMENTATION
0841 C          OF THE LENZ-THOMPSON ALGORITHM
0842 C B      : AUXILIARY VARIABLE IN THE IMPLEMENTATION
0843 C          OF THE LENZ-THOMPSON ALGORITHM
0844 C CC     : AUXILIARY VARIABLE IN THE IMPLEMENTATION
0845 C          OF THE LENZ-THOMPSON ALGORITHM
0846 C COSTH  : COS(THET)
0847 C D      : AUXILIARY VARIABLE IN THE IMPLEMENTATION
0848 C          OF THE LENZ-THOMPSON ALGORITHM
0849 C DELS   : AUXILIARY VARIABLE IN THE IMPLEMENTATION
0850 C          OF THE LENZ-THOMPSON ALGORITHM
0851 C DELTA  : AUXILIARY VARIABLE IN THE IMPLEMENTATION
0852 C          OF THE LENZ-THOMPSON ALGORITHM
0853 C FC     : CONTINUED FRACTION FOR K(PNU+1)/K(PNU)
0854 C FDOMIN : X*(COS(THET)+THET*SIN(THET))
0855 C K      : KIA FUNCTION
0856 C KP     : KIA' FUNCTION
0857 C PI     : 3.1415..
0858 C PIA    : PI*PNU
0859 C PISQ   : SQRT(PI)
0860 C PRECI  : RELATIVE PRECISION USED IN THE CALCULATION
0861 C Q0B    : AUXILIARY VARIABLE IN THE IMPLEMENTATION
0862 C          OF THE LENZ-THOMPSON ALGORITHM
0863 C Q1B    : AUXILIARY VARIABLE IN THE IMPLEMENTATION
0864 C          OF THE LENZ-THOMPSON ALGORITHM
0865 C QB     : AUXILIARY VARIABLE IN THE IMPLEMENTATION
0866 C          OF THE LENZ-THOMPSON ALGORITHM
0867 C QQB    : AUXILIARY VARIABLE IN THE IMPLEMENTATION
0868 C          OF THE LENZ-THOMPSON ALGORITHM
0869 C S      : AUXILIARY VARIABLE IN THE IMPLEMENTATION
0870 C          OF THE LENZ-THOMPSON ALGORITHM
0871 C SINTH  : SIN(THET)
0872 C THET   : ASIN(PNU/X)
0873 C UNDER  : UNDERFLOW NUMBER
0874 C Z0     : 1/(1+S)
0875 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0876        DOUBLE PRECISION A,A2,AA,AB,B,CC,COSTH,D,D1MACH,
0877      + DELS,DELTA,FC,FDOMIN,K,KP,PI,PIA,PISQ,PNU,PRECI,
0878      + PSER,PSERD,Q0B,Q1B,QB,QQB,S,SINTH,THET,UNDER,X,Z0
0879        INTEGER IERRO,IFAC,MM
0880 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
0881        UNDER=D1MACH(1)*1.D+8
0882        PI=ACOS(-1.D0)
0883        IERRO=0
0884        IF (IFAC.EQ.1) THEN
0885          IF (X.GE.PNU) THEN
0886            SINTH=PNU/X
0887            THET=ASIN(SINTH)
0888            COSTH=COS(THET)
0889            FDOMIN=X*(COSTH+THET*SINTH)
0890            IF (-FDOMIN.LE.LOG(UNDER)) IERRO=1
0891          ELSE
0892            IF ((-PI*PNU*0.5D0).LE.LOG(UNDER)) IERRO=1
0893          ENDIF
0894        ENDIF
0895        IF (IERRO.EQ.0) THEN
0896          PRECI=D1MACH(3)*10
0897          A=PNU
0898          A2=A*A
0899          PISQ=PI**0.5D0
0900 CCC CF FOR K(NU+1)/K(NU)
0901          MM=2
0902          AB=-(0.25D0+A2)
0903          AA=AB
0904          CC=-AB
0905          Q0B=0.D0
0906          Q1B=1.D0
0907          QQB=CC
0908          B=2.D0*(1.D0+X)
0909          D=1.D0/B
0910          DELTA=D
0911          FC=DELTA
0912          S=QQB*DELTA
0913          AB=-2.D0+AB
0914          B=B+2.D0
0915  91      D=1.D0/(B+AB*D)
0916          DELTA=(B*D-1.D0)*DELTA
0917          FC=FC+DELTA
0918          CC=-AB*CC/MM
0919          QB=(Q0B-(B-2.D0)*Q1B)/AB
0920          Q0B=Q1B
0921          Q1B=QB
0922          QQB=QQB+CC*Q1B
0923          DELS=QQB*DELTA
0924          S=S+DELS
0925          B=B+2.D0
0926          AB=-2.D0*MM+AB
0927          MM=MM+1
0928          IF (MM.LT.10000) THEN
0929            IF (ABS(DELS/S).GT.PRECI) GOTO 91
0930          ENDIF
0931          Z0=1.D0/(1.D0+S)
0932          IF (IFAC.EQ.1) THEN
0933            K=PISQ*(2.D0*X)**(-0.5D0)*EXP(-X)*Z0
0934            KP=-K/X*(0.5D0+X-(A2+0.25D0)*FC)
0935          ELSE
0936            IF (X.LT.A) THEN
0937              PIA=PI*A
0938              K=PISQ*(2.D0*X)**(-0.5D0)*EXP(-X+PIA*0.5D0)*Z0
0939              KP=-K/X*(0.5D0+X-(A2+0.25D0)*FC)
0940            ELSE
0941              SINTH=A/X
0942              THET=ASIN(SINTH)
0943              COSTH=COS(THET)
0944              FDOMIN=X*(COSTH+THET*SINTH)
0945              K=PISQ*(2.D0*X)**(-0.5D0)*EXP(-X+FDOMIN)*Z0
0946              KP=-K/X*(0.5D0+X-(A2+0.25D0)*FC)
0947            ENDIF
0948          ENDIF
0949          PSER=K
0950          PSERD=KP
0951        ELSE
0952          PSER=0.D0
0953          PSERD=0.D0
0954        ENDIF
0955        RETURN
0956        END
0957
0958        SUBROUTINE SERKIA(IFAC,X,PNU,PSER,PSERD,IERRO)
0959 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0960 CCC   CALCULATION OF POWER SERIES FOR K, K'
0961 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0962 CCC   INPUT VARIABLES:                              C
0963 CCC      X:   ARGUMENT OF THE FUNCTIONS             C
0964 CCC      PNU: ORDER OF THE FUNCTIONS                C
0965 CCC      IFAC:                                      C
0966 CCC               =1,  NON SCALED FUNCTIONS         C
0967 CCC               OTHERWISE,  SCALED FUNCTIONS      C
0968 CCC   OUTPUT VARIABLES:                             C
0969 CCC      PSER: K FUNCTION                           C
0970 CCC      PSERD: DERIVATIVE OF THE K FUNCTION        C
0971 CCC      IERRO: ERROR FLAG                          C
0972 CCC      * IF IERRO=0, COMPUTATION SUCCESSFUL.      C
0973 CCC      * IF IERRO=1, COMPUTATION OUT OF RANGE.    C
0974 CCC      * IF IERRO=2, ARGUMENT AND/OR ORDER,       C
0975 CCC                    OUT OF RANGE.                C
0976 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0977 CCC         LOCAL VARIABLES:                        C
0978 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0979 C A       : (PNU) ORDER OF THE FUNCTIONS
0980 C A2      : A**2
0981 C A2H     : A**2/2
0982 C A2N     : A**(2*N)
0983 C ACCP    : ACCUMULATES THE P COEFFICIENTS
0984 C ACCQ    : ACCUMULATES THE Q COEFFICIENTS
0985 C ARGU    : SIG(0)-A*LOG(X/2)
0986 C C       : X**(2*K)/K!
0987 C COCI    : 1/(K**2+A**2)
0988 C COSTH   : COS(THET)
0989 C DELTAK  : ACCUMULATES THE CONTRIBUTION FOR THE
0990 C           KIA(X) FUNCTION
0991 C DELTKP  : ACCUMULATES THE CONTRIBUTION FOR THE
0992 C           KIA'(X) FUNCTION
0993 C DF1     : FACTOR (DEPENDING ON IFAC)
0994 C ETA0    : PARAMETER FOR THE CALCULATION OF THE
0995 C           COULOMB PHASE SHIFT
0996 C ETA02   : ETA0**2
0997 C F(K)    : SIN(PHI(A,K)-A*LOG(X/2))
0998 C           /(A**2*(1+A**2)...(K**2+A**2))**1/2,
0999 C           WHERE PHI(A,K)=PHASE(GAMMA(1+K+IA))
1000 C FDOMIN  : X*(COS(THET)+THET*SIN(THET))
1001 C K       : CONTRIBUTION TO THE KIA(X) FUNCTION
1002 C KP      : CONTRIBUTION TO THE KIA'(X) FUNCTION
1003 C OVER    : OVERFLOW NUMBER
1004 C P0      : PARAMETERS FOR THE CALCULATION OF THE COULOMB
1005 C           PHASE SHIFT
1006 C P1      : PARAMETERS FOR THE CALCULATION OF THE COULOMB
1007 C           PHASE SHIFT
1008 C P2      : PARAMETERS FOR THE CALCULATION OF THE COULOMB
1009 C           PHASE SHIFT
1010 C PI      : 3.1415...
1011 C PIA     : PI*A
1012 C PIA2    : 2*PI*A
1013 C PRECI   : RELATIVE PRECISION USED IN THE CALCULATION
1014 C Q0      : PARAMETERS FOR THE CALCULATION OF THE COULOMB
1015 C           PHASE SHIFT
1016 C Q1      : PARAMETERS FOR THE CALCULATION OF THE COULOMB
1017 C           PHASE SHIFT
1018 C Q2      : PARAMETERS FOR THE CALCULATION OF THE COULOMB
1019 C           PHASE SHIFT
1020 C R(K)    : F(K)*A/TAN(PHI(A,K)-A*LOG(X/2))
1021 C SIG0    : COULOMB PHASE SHIFT
1022 C SINTH   : SIN(THET)
1023 C THET    : ASIN(A/X)
1024 C UNDER   : UNDERFLOW NUMBER
1025 C X2      : X*X
1026 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1027        DOUBLE PRECISION A,A2,A2N,ACCP,ACCQ,ARGU,C,COCI,COSTH,
1028      + D1MACH,DDS,DEE,DELTAK,DELTKP,DF1,DSMALL,ETA0,ETA02,
1029      + EULER,F(0:500),FDOMIN,K,KP,OVER,P0(0:9),P1(0:9),P2(0:6),
1030      + PI,PIA,PIA2,PNU,PRECI,PSER,PSERD,Q0(0:9),Q1(0:9),Q2(0:6),
1031      + R(0:500),SIG0,SINTH,THET,UNDER,X,X2
1032        INTEGER IERRO,IFAC,L,M,N
1033        SAVE P0,Q0,P1,Q1,P2,Q2
1034        DATA P0/1.08871504904797411683D5,3.64707573081160914640D5,
1035      +        4.88801471582878013158D5,3.36275736298197324009D5,
1036      +        1.26899226277838479804D5,
1037      +        2.60795543527084582682D4,2.73352480554497990544D3,
1038      +        1.26447543569902963184D2,
1039      +        1.85446022125533909390D0,1.90716219990037648146D-3/
1040        DATA Q0/6.14884786346071135090D5,2.29801588515708014282D6,
1041      +        3.50310844128424021934D6,
1042      +        2.81194990286041080264D6,1.28236441994358406742D6,
1043      +        3.35209348711803753154D5,
1044      +        4.84319580247948701171D4,3.54877039006873206531D3,
1045      +        1.11207201299804390166D2,1.D0/
1046        DATA P1/-1.044100987526487618670D10,-1.508574107180079913696D10,
1047      +        -5.582652833355901160542D9,4.052529174369477275446D8,
1048      +         5.461712273118594275192D8,
1049      +         9.510404403068169395714D7,6.281126609997342119416D6,
1050      +         1.651178048950518520416D5,
1051      +         1.498824421329341285521D3,2.974686506595477984776D0/
1052        DATA Q1/1.808868161493543887787D10,3.869142051704700267785D10,
1053      +        3.003264575147162634046D10,1.075554651494601843525D10,
1054      +        1.901298501823290694245D9,
1055      +        1.665999832151229472632D8,6.952188089169487375936D6,
1056      +        1.253235080625688652718D5,7.904420414560291396996D2,1.D0/
1057        DATA P2/7.08638611024520906826D-3,-6.54026368947801591128D-2,
1058      +         2.92684143106158043933D-1,4.66821392319665609167D0,
1059      +        -3.43943790382690949054D0,
1060      +        -7.72786486869252994370D0,-9.88841771200290647461D-01/
1061        DATA Q2/-7.08638611024520908189D-3,6.59931690706339630254D-2,
1062      +        -2.98754421632058618922D-1,-4.63752355513412248006D0,
1063      +        3.79700454098863541593D0,7.06184065426336718524D0,1.D0/
1064        PI=ACOS(-1.D0)
1065        ETA0=1.8055470716051069198764D0
1066        EULER=0.577215664901532860606512D0
1067        PRECI=D1MACH(3)*10
1068 CCC MACHINE DEPENDENT CONSTANT (OVERFLOW NUMBER)
1069        OVER=D1MACH(2)*1.D-8
1070 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
1071        UNDER=D1MACH(1)*1.D+8
1072        IERRO=0
1073        IF (IFAC.EQ.1) THEN
1074          IF (X.GE.PNU) THEN
1075            SINTH=PNU/X
1076            THET=ASIN(SINTH)
1077            COSTH=COS(THET)
1078            FDOMIN=X*(COSTH+THET*SINTH)
1079            IF (-FDOMIN.LE.LOG(UNDER)) IERRO=1
1080          ELSE
1081            IF ((-PI*PNU*0.5D0).LE.LOG(UNDER)) IERRO=1
1082          ENDIF
1083        ENDIF
1084 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1085 CCC COEFFICIENTS FOR THE CALCULATION OF THE COULOMB PHASE SHIFT
1086 CCC FROM CODY & HILLSTROM, MATH. COMPUT. 24(111) 1970
1087 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1088        IF (IERRO.EQ.0) THEN
1089          A=PNU
1090          ETA02=ETA0*ETA0
1091          A2=A*A
1092          N=0
1093          ACCP=0.D0
1094          ACCQ=0.D0
1095          IF (A.LE.2.D0) THEN
1096  33        A2N=A2**N
1097            ACCP=ACCP+P0(N)*A2N
1098            ACCQ=ACCQ+Q0(N)*A2N
1099            N=N+1
1100            IF (N.LE.9) GOTO 33
1101            SIG0=A*(A2-ETA02)*ACCP/ACCQ
1102          ELSE
1103            IF ((A.GT.2.D0).AND.(A.LE.4.D0)) THEN
1104  44          A2N=A2**N
1105              ACCP=ACCP+P1(N)*A2N
1106              ACCQ=ACCQ+Q1(N)*A2N
1107              N=N+1
1108              IF (N.LE.9) GOTO 44
1109              SIG0=A*ACCP/ACCQ
1110            ELSE
1111  55          A2N=A2**N
1112              ACCP=ACCP+P2(N)/A2N
1113              ACCQ=ACCQ+Q2(N)/A2N
1114              N=N+1
1115              IF (N.LE.6) GOTO 55
1116              SIG0=ATAN(A)*0.5D0+A*(LOG(1.D0+A2)*0.5D0+ACCP/ACCQ)
1117            ENDIF
1118          ENDIF
1119 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1120 CCC EVALUATION OF F(0), R(0), R(1)
1121 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1122          PIA=PI*A
1123          ARGU=SIG0-A*LOG(X*0.5D0)
1124          R(0)=COS(ARGU)
1125          R(1)=1.D0/(1.D0+A2)*(R(0)-A*SIN(ARGU))
1126          IF (A.LT.UNDER) THEN
1127            F(0)=-(EULER+LOG(X*0.5D0))
1128          ELSE
1129            F(0)=1.D0/A*SIN(ARGU)
1130          ENDIF
1131          C=1.D0
1132          K=F(0)
1133          KP=-0.5D0*R(0)
1134 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1135 CCC RECURSION FOR F(K), R(K), C(K)
1136 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1137          X2=0.25D0*X*X
1138          M=1
1139          COCI=1.D0/(M*M+A2)
1140          DELTAK=K*10
1141          DELTKP=KP*10
1142  66      F(M)=(M*F(M-1)+R(M-1))*COCI
1143          C=X2*C/M
1144          DELTAK=F(M)*C
1145          K=K+DELTAK
1146          DELTKP=(M*F(M)-0.5D0*R(M))*C
1147          KP=KP+DELTKP
1148          M=M+1
1149          COCI=1.D0/(M*M+A2)
1150          R(M)=((2.D0*M-1.D0)*R(M-1)-R(M-2))*COCI
1151          IF (K.GT.OVER) M=500
1152          IF (KP.GT.OVER) M=500
1153          IF (M.LT.500) THEN
1154            IF ((ABS(DELTAK/K).GT.PRECI).OR.
1155      +       (ABS(DELTKP/KP).GT.PRECI)) GOTO 66
1156          ENDIF
1157          PIA2=2.D0*PIA
1158          IF (IFAC.EQ.1) THEN
1159            IF (-PIA2.LE.LOG(UNDER)) THEN
1160              DEE=0.D0
1161            ELSE
1162              DEE=EXP(-PIA2)
1163            ENDIF
1164            IF (A.LT.1.D-1) THEN
1165              IF (A.LT.UNDER) THEN
1166                DF1=1.D0
1167              ELSE
1168                L=0
1169                DDS=1.D0
1170                DSMALL=1.D0
1171  47            L=L+1
1172                DDS=-DDS*PIA2/(L+1.D0)
1173                DSMALL=DSMALL+DDS
1174                IF (ABS(DDS/DSMALL).GT.PRECI) GOTO 47
1175                DF1=EXP(PIA*0.5D0)*SQRT(DSMALL)
1176              ENDIF
1177            ELSE
1178              DF1=EXP(PIA*0.5D0)*SQRT((1.D0-DEE)/PIA2)
1179            ENDIF
1180            PSER=K/DF1
1181            PSERD=KP*2.D0/X/DF1
1182          ELSE
1183            IF (X.LT.A) THEN
1184              IF (-PIA2.LE.LOG(UNDER)) THEN
1185                DEE=0.D0
1186              ELSE
1187                DEE=EXP(-PIA2)
1188              ENDIF
1189              IF (A.LT.1.D-1) THEN
1190                IF (A.LT.UNDER) THEN
1191                  DF1=1.D0
1192                ELSE
1193                  L=0
1194                  DDS=1.D0
1195                  DSMALL=1.D0
1196  48              L=L+1
1197                  DDS=-DDS*PIA2/(L+1.D0)
1198                  DSMALL=DSMALL+DDS
1199                  IF (ABS(DDS/DSMALL).GT.PRECI) GOTO 48
1200                  DF1=SQRT(DSMALL)
1201                ENDIF
1202              ELSE
1203                DF1=SQRT((1.D0-DEE)/PIA2)
1204              ENDIF
1205            ELSE
1206              SINTH=A/X
1207              THET=ASIN(SINTH)
1208              COSTH=COS(THET)
1209              FDOMIN=X*(COSTH+THET*SINTH)
1210              IF (-PIA2.LE.LOG(UNDER)) THEN
1211                DEE=0.D0
1212              ELSE
1213                DEE=EXP(-PIA2)
1214              ENDIF
1215              IF (A.LT.1.D-1) THEN
1216                IF (A.LT.UNDER) THEN
1217                  DF1=1.D0
1218                ELSE
1219                  L=0
1220                  DDS=1.D0
1221                  DSMALL=1.D0
1222  49              L=L+1
1223                  DDS=-DDS*PIA2/(L+1.D0)
1224                  DSMALL=DSMALL+DDS
1225                  IF (ABS(DDS/DSMALL).GT.PRECI) GOTO 49
1226                  DF1=EXP(PIA*0.5D0-FDOMIN)*SQRT(DSMALL)
1227                ENDIF
1228              ELSE
1229                DF1=EXP(PIA*0.5D0-FDOMIN)*SQRT((1.D0-DEE)/PIA2)
1230              ENDIF
1231            ENDIF
1232            PSER=K/DF1
1233            PSERD=KP*2.D0/X/DF1
1234          ENDIF
1235        ELSE
1236          PSER=0.D0
1237          PSERD=0.D0
1238        ENDIF
1239        RETURN
1240        END
1241
1242        SUBROUTINE AIEXKI(IFAC,X,A,DKAI,DKAID,IERROK)
1243 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1244 CCC AIRY-TYPE ASYMPTOTIC EXPANSION FOR THE K AND K'      C
1245 CCC FUNCTIONS                                            C
1246 CCC                                                      C
1247 CCC   INPUT VARIABLES:                                   C
1248 CCC      X:   ARGUMENT OF THE FUNCTIONS                  C
1249 CCC      A:   ORDER OF THE FUNCTIONS                     C
1250 CCC      IFAC:                                           C
1251 CCC               =1,  NON SCALED FUNCTIONS              C
1252 CCC               OTHERWISE,  SCALED FUNCTIONS           C
1253 CCC   OUTPUT VARIABLES:                                  C
1254 CCC      DKAI: K FUNCTION                                C
1255 CCC      DKAID: DERIVATIVE OF THE K FUNCTION             C
1256 CCC      IERROK: ERROR FLAG                              C
1257 CCC      * IF IERROK=0, COMPUTATION SUCCESSFUL.          C
1258 CCC      * IF IERROK=1, COMPUTATION OUT OF RANGE.        C
1259 CCC      * IF IERROK=2, ARGUMENT AND/OR ORDER,           C
1260 CCC                     OUT OF RANGE.                    C
1261 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1262 CCC           LOCAL VARIABLES:                           C
1263 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1264 C A0EX    : EXACT VALUE OF THE COEFFICIENT A_0
1265 C A13     : A**(1/3)
1266 C A2      : A**2
1267 C A23     : A**(2/3)
1268 C A2K     : A**(2*K)
1269 C AC      : COEFFICIENTS A_S(PSI) USED IN THE EXPANSIONS
1270 C AII     : IMAGINARY PART OF THE AIRY FUNCTION AI(Z)
1271 C AIR     : REAL PART OF THE AIRY FUNCTION AI(Z)
1272 C APIHAL  : A*PI/2
1273 C ARG     : -A**(2/3)*PSI
1274 C AS      : ACCUMULATES THE CONTRIBUTION OF THE AC COEFFICIENTS
1275 C           (FOR THE KIA(X) FUNCTION)
1276 C ASP     : ACCUMULATES THE CONTRIBUTION OF THE AC COEFFICIENTS
1277 C           (FOR THE KIA'(X) FUNCTION)
1278 C B0EX    : EXACT VALUE OF THE COEFFICIENT B_0
1279 C B0PEX   : EXACT VALUE OF THE COEFFICIENT B'_0
1280 C BC      : COEFFICIENTS B_S(PSI) USED IN THE EXPANSIONS
1281 C BS      : ACCUMULATES THE CONTRIBUTION OF THE BC COEFFICIENTS
1282 C           (FOR THE KIA(X) FUNCTION)
1283 C BSO     : ACCUMULATES THE OLD VALUES OF BS
1284 C BSP     : ACCUMULATES THE CONTRIBUTION OF THE BC COEFFICIENTS
1285 C           (FOR THE KIA'(X) FUNCTION)
1286 C BSPO    : ACCUMULATES THE OLD VALUES OF BSP
1287 C C0EX    : EXACT VALUE OF THE COEFFICIENT C_0
1288 C CHI     : ACCUMULATES THE CONTRIBUTION OF THE CHIN COEFFICIENTS
1289 C CHIEX   : EXACT VALUE OF THE FUNCTION CHI(PSI)
1290 C CHIN    : COEFFICIENTS FOR THE EXPANSION OF THE FUNCTION
1291 C           CHI(PSI)
1292 C COSTH   : COS(THET)
1293 C D0EX    : EXACT VALUE OF THE COEFFICIENT D_0
1294 C DAII    : IMAGINARY PART OF THE AIRY FUNCTION AI'(Z)
1295 C DAIR    : REAL PART OF THE AIRY FUNCTION AI'(Z)
1296 C DF      : 2**(1/3)
1297 C DZZ     : (ABS(1-Z**2))**(1/2)
1298 C ETA     : PSI/2**(1/3)
1299 C ETAJ    : ETA**J
1300 C ETAK    : ETA**K
1301 C ETAL    : ETA**L
1302 C EXPAM   : EXP(A*PI/2-FDOMIN)
1303 C EXPAPI  : EXP(A*PI/2)
1304 C F2      : (-)**K/A**(2*K)
1305 C F4      : (A**(1/3))**4
1306 C FAC     : FACTOR FOR THE KIA(X) FUNCTION
1307 C FACD    : FACTOR FOR THE KIA'(X) FUNCTION
1308 C FDOMIN  : X*(COS(THET)+THET*SIN(THET))
1309 C IERRO   : ERROR FLAG FOR THE AIRY FUNCTIONS
1310 C IFACA   :
1311 C           *IF IFACA=1, COMPUTATION OF UNSCALED AIRY AI(Z),AI'(Z)
1312 C           *IF IFACA=2, COMPUTATION OF SCALED AIRY AI(Z),AI'(Z)
1313 C IFUN    :
1314 C           *IF IFUN=1, COMPUTATION OF AIRY AI(Z)
1315 C           *IF IFUN=2, COMPUTATION OF AIRY AI'(Z)
1316 C PHI     : COEFFICIENTS FOR THE EXPANSION OF THE FUNCTION
1317 C           PHI(PSI)
1318 C PHIEX   : EXACT VALUE OF THE FUNCTION PHI(PSI)
1319 C PHIEX2  : PHI(PSI)**2
1320 C PHIS    : VALUE OF THE FUNCTION PHI(PSI)
1321 C PI      : 3.1415...
1322 C PIHALF  : PI/2
1323 C PSI     :
1324 C          *IF 0<=Z<=1, 2/3*PSI**3/2=LOG((1+SQRT(1-Z**2))/Z)
1325 C                       -SQRT(1-Z**2)
1326 C          *IF Z>=1, 2/3*(-PSI)**3/2=SQRT(Z**2-1)-ARCCOS(1/Z)
1327 C PSI12   : SQRT(PSI)
1328 C PSI2    : PSI*PSI
1329 C PSI3    : PSI**3
1330 C SAS     : ACCUMULATES THE CONTRIBUTION OF A_S(PSI)
1331 C SBS     : ACCUMULATES THE CONTRIBUTION OF B_S(PSI)
1332 C SCS     : ACCUMULATES THE CONTRIBUTION OF C_S(PSI)
1333 C SDS     : ACCUMULATES THE CONTRIBUTION OF D_S(PSI)
1334 C SIG     : (-)**K
1335 C SINTH   : SIN(THET)
1336 C THET    : ASIN(A/X)
1337 C UNDER   : UNDERFLOW NUMBER
1338 C Y       : Z-1
1339 C Z       : X/A
1340 C Z2      : Z**2
1341 C Z21M    : 1-Z**2
1342 C ZDS     : SQRT(1-Z**2)
1343 C ZMASF   : EXPANSION IN TERMS OF (Z-1) FOR THE
1344 C           CALCULATION OF PSI(Z)
1345 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1346        DOUBLE PRECISION A,A0EX,A13,A2,A23,A2K,AC(0:5,0:20),
1347      + AII,AIR,APIHAL,ARG,AS,ASP,B0EX,B0PEX,BC(0:5,0:20),
1348      + BS,BSO,BSP,BSPO,C0EX,CHI,CHIEX,CHIN(0:26),COSTH,
1349      + COZMAS(0:12),D0EX,D1MACH,DAII,DAIR,DF,DKAI,DKAID
1350        DOUBLE PRECISION DZZ,ETA,ETAJ,ETAK,ETAL,EXPAM,
1351      + EXPAPI,F2,F4,FAC,FACD,FDOMIN,PHI(0:35),PHIEX,PHIEX2,
1352      + PHIS,PI,PIHALF,PSI,PSI12,PSI2,PSI3,SAS,SBS,SCS,SDS,
1353      + SIG,SINTH,THET,UNDER,X,Y,Z,Z2,Z21M,ZDS,ZMASF
1354        INTEGER IERRO,IERROK,IFAC,IFACA,IFUN,INDA(0:5),INDB(0:5),
1355      + J,K,L
1356        SAVE PHI,CHIN,AC,BC
1357 CCCC VALUES OF THE COEFFICIENTS
1358        DATA PHI/1.D0,0.2D0,.25714285714285714286D-1,
1359      + -.56507936507936507937D-2,-.39367965367965367965D-2,
1360      + -.5209362066504924D-3,.3708541264731741D-3,
1361      + .2123827840293627D-3,.2150629788753145D-4,
1362      + -.2636904062981799D-4,-.1405469826493129D-4,
1363      + -.1149328899029441D-5,.1972641193938624D-5,
1364      + .1014324305593329D-5,.7034331100192200D-7,
1365      + -.1525044777392676D-6,-.7677866256900572D-7,
1366      + -.4669842638693018D-8,.1206673645965329D-7,
1367      + .5990877668092344D-8,.3269102150077715D-9,
1368      + -.97138350244428335095D-9,-.47745934295232233834D-9,
1369      + -.23750318642839155779D-10,.79244598109106655567D-10,
1370      + .38653584230817865528D-10,.17734105846426807873D-11,
1371      + -.65332110030864751956D-11,-.31673512094772575686D-11,
1372      + -.13524195177121030660D-12,.54324103217903338951D-12,
1373      + .26204918647967626464D-12,.10488134973664584898D-13,
1374      + -.45490420121539001435D-13,-.21851238232690420878D-13,
1375      + -.82456080260379042800D-15/
1376        DATA CHIN/.2D0,.1142857142857142857142857D-1,
1377      + -.2438095238095238095238095D-1,-.1003471449185734900020614D-1,
1378      + .8811404468547325690182833D-3,.2318339655809043564145605D-2,
1379      + .7794494413564441575646057D-3,-.1373952504077228744949558D-3,
1380      + -.2162301322540308393022684D-3,-.6585004634375583559042795D-4,
1381      + .1502851367097217578058824D-4,.1991904617871647675455262D-4,
1382      + .5731972706719818525615427D-5,-.1496427612747891044606885D-5,
1383      + -.1821231830428939670133992D-5,-.5054254875930882534538217D-6,
1384      + .1433283859497625931203415D-6,.1657681319078678321113361D-6,
1385      + .4485592642702540575627044D-7,-.1346138242826094098161839D-7,
1386      + -.1504601862773535117708677D-7,-.3995660198654805921651406D-8,
1387      + .1250124931952495738882300D-8,.1363187174221864073749532D-8,
1388      + .3567608459777506132029204D-9,-.1152749290139859167119863D-9,
1389      + -.1233547289799408517691696D-9/
1390        DATA COZMAS/.9428090415820634D0,-.4242640687119285D0,
1391      + .2904188565587606D0,-.2234261009999160D0,
1392      + .1821754153944745D0,-.1539625154198624D0,
1393      + .1333737583085222D0,-.1176617834148007D0,
1394      + .1052687194772381D0,-.9524025714638822D-1,
1395      + .8695738197073783D-1,-.8000034897653656D-1,
1396      + .7407421797273086D-1/
1397        DATA AC(0,0)/1.D0/
1398        DATA AC(1,0)/-.44444444444444444445D-2/
1399        DATA AC(1,1)/-.18441558441558441558D-2/
1400        DATA AC(1,2)/.11213675213675213675D-2/
1401        DATA AC(1,3)/.13457752124418791086D-2/
1402        DATA AC(1,4)/.0003880626562979504D0/
1403        DATA AC(1,5)/-.0001830686723781799D0/
1404        DATA AC(1,6)/-.0001995460887806733D0/
1405        DATA AC(1,7)/-.00005256191234041587D0/
1406        DATA AC(1,8)/.00002460619652459158D0/
1407        DATA AC(1,9)/.00002519246780924541D0/
1408        DATA AC(1,10)/.6333157376533242D-5/
1409        DATA AC(1,11)/-.2957485733830202D-5/
1410        DATA AC(1,12)/-.2925255920564838D-5/
1411        DATA AC(1,13)/-.7159702610502009D-6/
1412        DATA AC(1,14)/.3331510720390949D-6/
1413        DATA AC(1,15)/.3227670475692310D-6/
1414        DATA AC(1,16)/.7767729381664199D-7/
1415        DATA AC(1,17)/-.3600954237921120D-7/
1416        DATA AC(1,18)/-.3441724449034226D-7/
1417        DATA AC(1,19)/-.8188194356398772D-8/
1418        DATA AC(1,20)/.3783148485152038D-8/
1419        DATA AC(2,0)/.69373554135458897365D-3/
1420        DATA AC(2,1)/.46448349036584330703D-3/
1421        DATA AC(2,2)/-.42838130171535112460D-3/
1422        DATA AC(2,3)/-.0007026702868771135D0/
1423        DATA AC(2,4)/-.0002632580046778811D0/
1424        DATA AC(2,5)/.0001663853666288703D0/
1425        DATA AC(2,6)/.0002212087687818584D0/
1426        DATA AC(2,7)/.00007020345615329662D0/
1427        DATA AC(2,8)/-.00004000421782540614D0/
1428        DATA AC(2,9)/-.00004786324966453962D0/
1429        DATA AC(2,10)/-.00001394600741473631D0/
1430        DATA AC(2,11)/.7536186591273727D-5/
1431        DATA AC(2,12)/.8478502161067667D-5/
1432        DATA AC(2,13)/.2345355228453912D-5/
1433        DATA AC(2,14)/-.1225943294710883D-5/
1434        DATA AC(2,15)/-.1325082343401027D-5/
1435        DATA AC(2,16)/-.3539954776569997D-6/
1436        DATA AC(2,17)/.1808291719376674D-6/
1437        DATA AC(2,18)/.1900383515233655D-6/
1438        DATA AC(3,0)/-.35421197145774384076D-3/
1439        DATA AC(3,1)/-.31232252789031883276D-3/
1440        DATA AC(3,2)/.3716442237502298D-3/
1441        DATA AC(3,3)/.0007539269155977733D0/
1442        DATA AC(3,4)/.0003408300059444739D0/
1443        DATA AC(3,5)/-.0002634968172069594D0/
1444        DATA AC(3,6)/-.0004089275726648432D0/
1445        DATA AC(3,7)/-.0001501108759563460D0/
1446        DATA AC(3,8)/.00009964015205538056D0/
1447        DATA AC(3,9)/.0001352492955751283D0/
1448        DATA AC(3,10)/.00004443117087272903D0/
1449        DATA AC(3,11)/-.00002713205071914117D0/
1450        DATA AC(3,12)/-.00003396796969771860D0/
1451        DATA AC(3,13)/-.00001040708865273043D0/
1452        DATA AC(3,14)/.6024639065414414D-5/
1453        DATA AC(3,15)/.7143919607846883D-5/
1454        DATA AC(4,0)/.378194199201773D-3/
1455        DATA AC(4,1)/.000404943905523634D0/
1456        DATA AC(4,2)/-.000579130526946453D0/
1457        DATA AC(4,3)/-.00138017901171011D0/
1458        DATA AC(4,4)/-.000722520056780091D0/
1459        DATA AC(4,5)/.000651265924036825D0/
1460        DATA AC(4,6)/.00114674563328389D0/
1461        DATA AC(4,7)/.000474423189340405D0/
1462        DATA AC(4,8)/-.000356495172735468D0/
1463        DATA AC(4,9)/-.000538157791035111D0/
1464        DATA AC(4,10)/-.000195687390661225D0/
1465        DATA AC(4,11)/.000132563525210293D0/
1466        DATA AC(4,12)/.000181949256267291D0/
1467        DATA AC(5,0)/-.69114139728829416760D-3/
1468        DATA AC(5,1)/-.00085995326611774383285D0/
1469        DATA AC(5,2)/.0014202335568143511489D0/
1470        DATA AC(5,3)/.0038535426995603052443D0/
1471        DATA AC(5,4)/.0022752811642901374595D0/
1472        DATA AC(5,5)/-.0023219572034556988366D0/
1473        DATA AC(5,6)/-.0045478643394434635622D0/
1474        DATA AC(5,7)/-.0020824431758272449919D0/
1475        DATA AC(5,8)/.0017370443573195808719D0/
1476        DATA BC(0,0)/.14285714285714285714D-1/
1477        DATA BC(0,1)/.88888888888888888889D-2/
1478        DATA BC(0,2)/.20482374768089053803D-2/
1479        DATA BC(0,3)/-.57826617826617826618D-3/
1480        DATA BC(0,4)/-.60412089799844901886D-3/
1481        DATA BC(0,5)/-.0001472685745626922D0/
1482        DATA BC(0,6)/.00005324102148009784D0/
1483        DATA BC(0,7)/.00005206561006583416D0/
1484        DATA BC(0,8)/.00001233115050894939D0/
1485        DATA BC(0,9)/-.4905932728531366D-5/
1486        DATA BC(0,10)/-.4632230987136350D-5/
1487        DATA BC(0,11)/-.1077174523455235D-5/
1488        DATA BC(0,12)/.4475963978932822D-6/
1489        DATA BC(0,13)/.4152586188464624D-6/
1490        DATA BC(0,14)/.9555819293589234D-7/
1491        DATA BC(0,15)/-.4060599208403059D-7/
1492        DATA BC(0,16)/-.3731367187988482D-7/
1493        DATA BC(0,17)/-.8532670645553778D-8/
1494        DATA BC(0,18)/.3673017245573624D-8/
1495        DATA BC(0,19)/.3355960460784536D-8/
1496        DATA BC(0,20)/.7643107095110475D-9/
1497        DATA BC(1,0)/-.11848595848595848596D-2/
1498        DATA BC(1,1)/-.13940630797773654917D-2/
1499        DATA BC(1,2)/-.48141005586383737645D-3/
1500        DATA BC(1,3)/.26841705366016142958D-3/
1501        DATA BC(1,4)/.0003419706982709903D0/
1502        DATA BC(1,5)/.0001034548234902078D0/
1503        DATA BC(1,6)/-.00005418191982095504D0/
1504        DATA BC(1,7)/-.00006202184830690167D0/
1505        DATA BC(1,8)/-.00001724885886056087D0/
1506        DATA BC(1,9)/.8744675992887053D-5/
1507        DATA BC(1,10)/.9420684216180929D-5/
1508        DATA BC(1,11)/.2494922112085850D-5/
1509        DATA BC(1,12)/-.1238458608836357D-5/
1510        DATA BC(1,13)/-.1285461713809769D-5/
1511        DATA BC(1,14)/-.3299710862537507D-6/
1512        DATA BC(1,15)/.1613441105788315D-6/
1513        DATA BC(1,16)/.1633623194402374D-6/
1514        DATA BC(1,17)/.4104252949605779D-7/
1515        DATA BC(1,18)/-.1984317042326989D-7/
1516        DATA BC(1,19)/-.1973948142769706D-7/
1517        DATA BC(1,20)/-.4882194808588752D-8/
1518        DATA BC(2,0)/.43829180944898810994D-3/
1519        DATA BC(2,1)/.71104865116708668943D-3/
1520        DATA BC(2,2)/.31858383945387580576D-3/
1521        DATA BC(2,3)/-.0002404809426804458D0/
1522        DATA BC(2,4)/-.0003722966038621536D0/
1523        DATA BC(2,5)/-.0001352752059595618D0/
1524        DATA BC(2,6)/.00008691694372704142D0/
1525        DATA BC(2,7)/.0001158750753591377D0/
1526        DATA BC(2,8)/.00003724965927846447D0/
1527        DATA BC(2,9)/-.00002198334949606935D0/
1528        DATA BC(2,10)/-.00002686449633870452D0/
1529        DATA BC(2,11)/-.8023061612032524D-5/
1530        DATA BC(2,12)/.4494756592180126D-5/
1531        DATA BC(2,13)/.5193504763856015D-5/
1532        DATA BC(2,14)/.1477156191529617D-5/
1533        DATA BC(2,15)/-.7988793826096784D-6/
1534        DATA BC(3,0)/-.37670439477105454219D-3/
1535        DATA BC(3,1)/-.75856271658798642365D-3/
1536        DATA BC(3,2)/-.0004103253968775048D0/
1537        DATA BC(3,3)/.0003791263310429010D0/
1538        DATA BC(3,4)/.0006850981673903450D0/
1539        DATA BC(3,5)/.0002878310571932216D0/
1540        DATA BC(3,6)/-.0002157010636115705D0/
1541        DATA BC(3,7)/-.0003260863991373500D0/
1542        DATA BC(3,8)/-.0001181317008748678D0/
1543        DATA BC(3,9)/.00007887526841582582D0/
1544        DATA BC(3,10)/.0001072081833420685D0/
1545        DATA BC(3,11)/.00003544595251288735D0/
1546        DATA BC(3,12)/-.00002201447920733824D0/
1547        DATA BC(3,13)/-.00002789336359620813D0/
1548        DATA BC(4,0)/.00058453330122076187255D0/
1549        DATA BC(4,1)/.0013854690422372401251D0/
1550        DATA BC(4,2)/.00086830374184946900245D0/
1551        DATA BC(4,3)/-.00093502904801345951693D0/
1552        DATA BC(4,4)/-.0019175486005525492059D0/
1553        DATA BC(4,5)/-.00090795047113308137941D0/
1554        DATA BC(4,6)/.00077050429806392235104D0/
1555        DATA BC(4,7)/.0012953100128255983488D0/
1556        DATA BC(4,8)/.00051933869471899550762D0/
1557        DATA BC(4,9)/-.00038482631948759834653D0/
1558        DATA BC(4,10)/-.00057335393099012476502D0/
1559        DATA BC(5,0)/-.0014301070053470410656D0/
1560        DATA BC(5,1)/-.0038637811942002539408D0/
1561        DATA BC(5,2)/-.0027322816261168328245D0/
1562        DATA BC(5,3)/.0033294980346743452748D0/
1563        DATA BC(5,4)/.0075972237660887795911D0/
1564        DATA BC(5,5)/.0039816655673062060620D0/
1565        DATA BC(5,6)/-.0037510180460986006595D0/
1566        DATA INDA/0,20,18,15,12,8/
1567        DATA INDB/20,20,15,13,10,6/
1568 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
1569        UNDER=D1MACH(1)*1.D+8
1570 CCCC CONSTANTS CCCCCCCCCCCCCC
1571        PI=ACOS(-1.D0)
1572        PIHALF=PI*0.5D0
1573        IERROK=0
1574        IF (IFAC.EQ.1) THEN
1575          IF (X.GE.A) THEN
1576            SINTH=A/X
1577            THET=ASIN(SINTH)
1578            COSTH=COS(THET)
1579            FDOMIN=X*(COSTH+THET*SINTH)
1580            IF (-FDOMIN.LE.LOG(UNDER)) IERROK=1
1581          ELSE
1582            IF ((-PI*A*0.5D0).LE.LOG(UNDER)) IERROK=1
1583          ENDIF
1584        ENDIF
1585        IF (IERROK.EQ.0) THEN
1586          DF=2.D0**(1.D0/3.D0)
1587 CCCC VARIABLES CCCCCCCCCCCCCCC
1588          Z=X/A
1589          A2=A*A
1590          A13=A**(1.D0/3.D0)
1591          A23=A13*A13
1592          IF (IFAC.EQ.1) THEN
1593            APIHAL=A*PIHALF
1594            EXPAPI=EXP(-APIHAL)
1595            FAC=PI*EXPAPI/A13
1596            FACD=2.D0*PI*EXPAPI/Z/A23
1597          ELSE
1598            IF (X.LT.A) THEN
1599              FAC=PI/A13
1600              FACD=2.D0*PI/Z/A23
1601            ELSE
1602              SINTH=A/X
1603              THET=ASIN(SINTH)
1604              COSTH=COS(THET)
1605              FDOMIN=X*(COSTH+THET*SINTH)
1606              EXPAM=EXP(-A*PIHALF+FDOMIN)
1607              FAC=PI*EXPAM/A13
1608              FACD=2.D0*PI*EXPAM/Z/A23
1609            ENDIF
1610          ENDIF
1611          F4=A13**4
1612          IF (Z.LE.0.9D0) THEN
1613            ZDS=SQRT((1.D0-Z)*(1.D0+Z))
1614            PSI=(1.5D0*(LOG((1.D0+ZDS)/Z)-ZDS))**(2.D0/3.D0)
1615          ELSEIF (Z.GT.1.1D0) THEN
1616            ZDS=SQRT((Z-1.D0)*(1.D0+Z))
1617            PSI=-(1.5D0*(ZDS-ACOS(1.D0/Z)))**(2.D0/3.D0)
1618          ELSE
1619            Y=Z-1.D0
1620            ZMASF=COZMAS(12)
1621            DO 8 K=0,11
1622              J=11-K
1623              ZMASF=COZMAS(J)+Y*ZMASF
1624  8         CONTINUE
1625            ZMASF=ABS(Y)**(1.5D0)*ZMASF
1626            IF (Z.LT.1.D0) THEN
1627              PSI=(1.5D0*ZMASF)**(2.D0/3.D0)
1628            ELSE
1629              PSI=-(1.5D0*ZMASF)**(2.D0/3.D0)
1630            ENDIF
1631          ENDIF
1632          ETA=PSI/DF
1633          ARG=-A23*PSI
1634          PSI2=PSI*PSI
1635          PSI3=PSI2*PSI
1636          IF ((Z.GT.0.8D0).AND.(Z.LT.1.2D0)) THEN
1637            PHIS=0.D0
1638            CHI=0.D0
1639            SAS=0.D0
1640            SBS=0.D0
1641            SDS=0.D0
1642            SCS=1.D0
1643            BS=0.D0
1644            BSPO=0.D0
1645            DO 41 L=0,20
1646              IF (L.EQ.0) THEN
1647                ETAL=1.D0
1648              ELSE
1649                ETAL=ETAL*ETA
1650              ENDIF
1651              CHI=CHI+CHIN(L)*ETAL
1652  41        CONTINUE
1653            CHI=CHI/DF
1654            DO 10 K=0,20
1655              BSO=BS
1656              BSPO=BSP
1657              AS=0.D0
1658              BS=0.D0
1659              ASP=0.D0
1660              BSP=0.D0
1661              IF (K.EQ.0) THEN
1662                ETAK=1.D0
1663                A2K=1.D0
1664                SIG=1.D0
1665              ELSE
1666                ETAK=ETAK*ETA
1667                A2K=A2K*A2
1668                SIG=-1.D0*SIG
1669              ENDIF
1670              PHIS=PHIS+PHI(K)*ETAK
1671              F2=SIG/A2K
1672              IF (K.LE.5) THEN
1673                DO 20 J=0,20
1674                  IF (J.EQ.0) THEN
1675                    ETAJ=1.D0
1676                  ELSE
1677                    ETAJ=ETAJ*ETA
1678                  ENDIF
1679                  IF (J.LE.INDA(K)) THEN
1680                    AS=AS+AC(K,J)*ETAJ
1681                  ENDIF
1682                  IF (J.LE.INDB(K)) THEN
1683                    BS=BS+BC(K,J)*ETAJ
1684                  ENDIF
1685                  IF ((J+1).LE.INDA(K)) THEN
1686                    ASP=ASP+(J+1)*AC(K,J+1)*ETAJ
1687                  ENDIF
1688                  IF ((J+1).LE.INDB(K)) THEN
1689                    BSP=BSP+(J+1)*BC(K,J+1)*ETAJ
1690                  ENDIF
1691  20            CONTINUE
1692                ASP=1.D0/DF*ASP
1693                BS=BS*DF
1694                SAS=SAS+AS*F2
1695                SBS=SBS+BS*F2/DF
1696                SDS=SDS-(CHI*AS+ASP+PSI*BS)*F2
1697              ENDIF
1698              IF ((K.GT.0).AND.(K.LE.6)) THEN
1699                SCS=SCS+(AS+CHI*BSO+BSPO)*F2
1700              ENDIF
1701  10        CONTINUE
1702            PHIS=DF*PHIS
1703            SBS=DF*SBS
1704          ELSE
1705 CCCC EXACT VALUES CCCCCCCCCCCCCCCCCCCCC
1706            Z2=Z*Z
1707            Z21M=1.D0-Z2
1708            PHIEX=(4.D0*PSI/Z21M)**0.25D0
1709            PHIEX2=PHIEX*PHIEX
1710            A0EX=1.D0
1711            B0EX=-5.D0/48.D0/(PSI*PSI)+PHIEX2/48.D0/PSI*
1712      +    (5.D0/Z21M-3.D0)
1713            CHIEX=0.25D0/PSI*(1.D0-Z2*PHIEX2**3*0.25D0)
1714            CHI=CHIEX
1715            D0EX=-(7.D0/48.D0/PSI+PHIEX2/48.D0*
1716      +     (9.D0-7.D0/Z21M))
1717            IF (PSI.GT.0.D0) THEN
1718              PSI12=SQRT(PSI)
1719              DZZ=SQRT(Z21M)
1720            ELSE
1721              PSI12=SQRT(-PSI)
1722              DZZ=SQRT(ABS(Z21M))
1723            ENDIF
1724            B0PEX=5.D0/24.D0/PSI3+PHIEX2/48.D0*((2.D0*CHI*PSI-1.D0)/
1725      +       PSI2*(5.D0/Z21M-3.D0)-10.D0*Z2*PSI12/
1726      +       Z21M**2/DZZ/PSI)
1727            C0EX=1.D0
1728            SAS=A0EX
1729            SBS=B0EX
1730            SCS=C0EX
1731            SDS=D0EX
1732            BS=0.D0
1733            BSP=0.D0
1734            A2K=1.D0
1735            SIG=1.D0
1736            DO 30 K=1,6
1737              BSO=BS
1738              BSPO=BSP
1739              AS=0.D0
1740              BS=0.D0
1741              ASP=0.D0
1742              BSP=0.D0
1743              A2K=A2K*A2
1744              SIG=-1.D0*SIG
1745              F2=SIG/A2K
1746              IF (K.LE.5) THEN
1747                DO 40 J=0,20
1748                  IF (J.EQ.0) THEN
1749                    ETAJ=1.D0
1750                  ELSE
1751                    ETAJ=ETAJ*ETA
1752                  ENDIF
1753                  IF (J.LE.INDA(K)) THEN
1754                    AS=AS+AC(K,J)*ETAJ
1755                  ENDIF
1756                  IF (J.LE.INDB(K)) THEN
1757                    BS=BS+BC(K,J)*ETAJ
1758                  ENDIF
1759                  IF ((J+1).LE.INDA(K)) THEN
1760                    ASP=ASP+(J+1)*AC(K,J+1)*ETAJ
1761                  ENDIF
1762                  IF ((J+1).LE.INDB(K)) THEN
1763                    BSP=BSP+(J+1)*BC(K,J+1)*ETAJ
1764                  ENDIF
1765  40            CONTINUE
1766                BS=DF*BS
1767                ASP=ASP/DF
1768                SAS=SAS+AS*F2
1769                SBS=SBS+BS*F2
1770                SDS=SDS-(CHI*AS+ASP+PSI*BS)*F2
1771              ENDIF
1772              IF (K.GE.2) THEN
1773                SCS=SCS+(AS+CHI*BSO+BSPO)*F2
1774              ELSE
1775                SCS=SCS+(AS+CHI*B0EX+B0PEX)*F2
1776              ENDIF
1777  30        CONTINUE
1778            PHIS=PHIEX
1779          ENDIF
1780 CCCCC CALL THE AIRY ROUTINE CCCCCCCCC
1781          IFUN=1
1782          IFACA=1
1783          CALL AIZ(IFUN,IFACA,ARG,0.D0,AIR,AII,IERRO)
1784          IFUN=2
1785          IFACA=1
1786          CALL AIZ(IFUN,IFACA,ARG,0.D0,DAIR,DAII,IERRO)
1787          DKAI=FAC*PHIS*(AIR*SAS+DAIR*SBS/F4)
1788          DKAID=FACD/PHIS*(DAIR*SCS+AIR*SDS/A23)
1789        ELSE
1790          DKAI=0.D0
1791          DKAID=0.D0
1792        ENDIF
1793        RETURN
1794        END
1795
1796        SUBROUTINE DLIA(IFAC,X,A,DLI,DLID,IERRO)
1797 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1798 CCCC THE PURPOUSE OF THE ROUTINE IS THE CALCULATION OF THE
1799 CCCC MODIFIED BESSEL FUNCTIONS L_IA(X) AND L'_IA(X),
1800 CCCC WHERE I IS THE IMAGINARY UNIT AND X IS THE ARGUMENT OF THE
1801 CCCC FUNCTIONS. WE WILL REFER TO A AS THE ORDER OF THE FUNCTIONS.
1802 CCCC
1803 CCCC THE ROUTINE HAS THE OPTION OF COMPUTING SCALED FUNCTIONS.
1804 CCCC THIS SCALING CAN BE USED TO ENLARGE THE RANGE OF
1805 CCCC COMPUTATION.
1806 CCCC THE SCALED FUNCTIONS ARE DEFINED AS FOLLOWS:
1807 CCCC     (S STANDS FOR SCALED  AND
1808 CCCC      L=SQRT{X**2-A**2} + A*ARCSIN(A/X)):
1809 CCCC
1810 CCCC                EXP(-L)*L_IA(X),            IF X >=ABS(A)
1811 CCCC    SL_IA(X)  =
1812 CCCC                EXP(-ABS(A)*PI/2)*L_IA(X),  IF X < ABS(A)
1813 CCCC
1814 CCCC
1815 CCCC                EXP(-L)*L'_IA(X),           IF X >=ABS(A)
1816 CCCC    SL'_IA(X) =
1817 CCCC                EXP(-ABS(A)*PI/2)*L'_IA(X), IF X <ABS(A)
1818 CCCC
1819 CCCC THE RANGE OF THE PARAMETERS (X,A) FOR THE COMPUTATION IS:
1820 CCCC        0 < X <= 1500,  -1500 <= A <= 1500
1821 CCCC FOR FUNCTIONS WITHOUT SCALING, THE RANGE IS USUALLY LARGER
1822 CCCC THAN
1823 CCCC        0 < X <= 500,  -400 <= A <= 400,
1824 CCCC DEPENDING ON THE MACHINE OVERFLOW/UNDERFLOW PARAMETERS, WHICH
1825 CCCC ARE SET UP BY THE ROUTINE D1MACH.
1826 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1827 CCCC              METHODS OF COMPUTATION:
1828 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1829 CCCC      1) SERIES,
1830 CCCC          IF  ABS(A) > 0.044*ABS(X-3.1)**1.9+(X-3.1) OR
1831 CCCC              ABS(A) <= 25 AND X <= 60
1832 CCCC      2) ASYMPTOTIC EXPANSION, IF ABS(A) < 0.7*X-10
1833 CCCC      3) AIRY-TYPE ASYMPTOTIC EXPANSION,
1834 CCCC             IF ABS(A) < 3.7*X-10
1835 CCCC      4) QUADRATURES,
1836 CCCC           IN THE REST OF THE PLANE (X,A)
1837 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1838 CCCC        DESCRIPTION OF INPUT/OUTPUT VARIABLES:
1839 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1840 CCCC  INPUTS:
1841 CCCC    X      : ARGUMENT OF THE FUNCTIONS
1842 CCC     A      : ORDER OF THE FUNCTIONS
1843 CCCC    IFAC:   INTEGER FLAG FOR THE SCALING
1844 CCCC     * IFAC=1, THE CODE COMPUTES L_IA(X), L'_IA(X)
1845 CCCC     * OTHERWISE, THE CODE COMPUTES SCALED L_IA(X), L'_IA(X)
1846 CCCC  OUTPUTS:
1847 CCCC    DLI    :  LIA(X) FUNCTION
1848 CCCC    DLID   :  DERIVATIVE OF THE LIA(X) FUNCTION
1849 CCCC    IERRO: ERROR FLAG
1850 CCCC     * IF IERRO=0, COMPUTATION SUCCESSFUL.
1851 CCCC     * IF IERRO=1, COMPUTATION OUT OF RANGE.
1852 CCCC     * IF IERRO=2, ARGUMENT X AND/OR ORDER A  OUT OF RANGE.
1853 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1854 CCCC           ACCURACY:
1855 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1856 CCCC THE RELATIVE ACCURACY IS:
1857 CCCC  BETTER THAN 10**(-13)   FOR (X,A) IN (0,200]X[-200,200];
1858 CCCC  BETTER THAN 5.10**(-13) FOR (X,A) IN (0,500]X[-500,500];
1859 CCCC  CLOSE TO 10**(-12)      FOR (X,A) IN (0,1500]X[-1500,1500].
1860 CCCC NEAR THE ZEROS OF THE FUNCTIONS (THERE ARE INFINITELY
1861 CCCC MANY OF THEM IN THE OSCILLATORY REGION) RELATIVE PRECISION
1862 CCCC LOOSES MEANING AND ONLY ABSOLUTE PRECISION MAKES SENSE.
1863 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1864 C     AUTHORS:
1865 C        AMPARO GIL    (U. CANTABRIA, SANTANDER, SPAIN).
1866 C                      E-MAIL: AMPARO.GIL@UNICAN.ES
1867 C        JAVIER SEGURA (U. CANTABRIA, SANTANDER, SPAIN).
1868 C                      E-MAIL: SEGURAJJ@UNICAN.ES
1869 C        NICO M. TEMME (CWI, AMSTERDAM, THE NETHERLANDS).
1870 C                      E-MAIL: NICO.TEMME@CWI.NL
1871 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1872 C    REFERENCES:
1873 C     THIS IS THE COMPANION SOFTWARE OF THE ARTICLES
1874 C      1)'COMPUTING SOLUTIONS OF THE MODIFIED BESSEL DIFFERENTIAL
1875 C         EQUATION FOR IMAGINARY ORDERS AND POSITIVE ARGUMENTS',
1876 C         A. GIL, J. SEGURA, N.M. TEMME
1877 C         ACM TRANS. MATH. SOFT. (2004)
1878 C      2)'MODIFIED BESSEL FUNCTIONS OF IMAGINARY ORDER AND
1879 C         POSITIVE ARGUMENT',
1880 C         A. GIL, J. SEGURA, N.M. TEMME
1881 C         ACM TRANS. MATH. SOFT. (2004)
1882 C    ADDITIONAL REFERENCES:
1883 C     - 'COMPUTATION OF THE MODIFIED BESSEL FUNCTION OF THE
1884 C         THIRD KIND FOR IMAGINARY ORDERS'
1885 C         A. GIL, J. SEGURA, N.M. TEMME
1886 C         J. COMPUT. PHYS. 175 (2002) 398-411
1887 C     - 'COMPUTATION OF THE MODIFIED BESSEL FUNCTIONS OF THE
1888 C         THIRD KIND OF IMAGINARY ORDERS:
1889 C         UNIFORM AIRY-TYPE ASYMPTOTIC EXPANSION'
1890 C         A. GIL, J. SEGURA, N.M. TEMME, PROCEEDINGS OPSFA 2001
1891 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1892 CCC        LOCAL VARIABLES:
1893 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1894 C      D:   X-3.1
1895 C      DF1: 0.044*ABS(D)**1.9+(X-3.1)
1896 C      DF2: 0.7*X-10
1897 C      DF3: 3.7*X-10
1898 C      PNU: ORDER OF THE FUNCTION
1899 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1900        DOUBLE PRECISION A,D,DF1,DF2,DF3,DLI,DLID,PNU,X
1901        INTEGER IERRO,IFAC
1902        IERRO=0
1903        PNU=A
1904        IF ((X.GT.1500.D0).OR.(X.LE.0.D0)) THEN
1905          IERRO=2
1906          DLI=0.D0
1907          DLID=0.D0
1908        ENDIF
1909        IF (ABS(PNU).GT.1500.D0) THEN
1910          IERRO=2
1911          DLI=0.D0
1912          DLID=0.D0
1913        ENDIF
1914        IF (IERRO.EQ.0) THEN
1915 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1916 CCC THESE FUNCTIONS ARE EVEN FUNCTIONS IN THE  C
1917 CCC PARAMETER A (=PNU)                         C
1918 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1919          IF (PNU.LT.0.D0) PNU=-PNU
1920          D=X-3.1D0
1921          DF1=0.044D0*ABS(D)**1.9D0+D
1922          DF2=0.7D0*X-10.D0
1923          DF3=3.7D0*X-10.D0
1924          IF ((PNU.GT.DF1).OR.(((PNU.LE.25.D0).AND.(X.LE.60.D0))))
1925      +   THEN
1926            CALL SERLIA(IFAC,X,PNU,DLI,DLID,IERRO)
1927          ELSEIF (PNU.LT.DF2) THEN
1928            CALL EXPLIA(IFAC,X,PNU,DLI,DLID,IERRO)
1929          ELSEIF (PNU.LT.DF3) THEN
1930            CALL AIEXLI(IFAC,X,PNU,DLI,DLID,IERRO)
1931          ELSE
1932            CALL DLINT(IFAC,X,PNU,DLI,DLID,IERRO)
1933          ENDIF
1934        ENDIF
1935        RETURN
1936        END
1937
1938        SUBROUTINE DLINT(IFAC,XX,PNUA,DLINF,DLIND,IERRO)
1939 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1940 CCC  CALCULATION OF L, L' USING NON-OSCILLATING INTEGRAL       C
1941 CCC  REPRESENTATIONS                                           C
1942 CCC                                                            C
1943 CCC  INPUT:                                                    C
1944 CCC      XX:   ARGUMENT OF THE FUNCTIONS                       C
1945 CCC      PNUA: ORDER OF THE FUNCTIONS                          C
1946 CCC      IFAC:                                                 C
1947 CCC            =1,  NON SCALED FUNCTIONS                       C
1948 CCC            OTHERWISE,  SCALED FUNCTIONS                    C
1949 CCC  OUTPUT:                                                   C
1950 CCC      DLINF: L FUNCTION                                     C
1951 CCC      DLIND: DERIVATIVE OF THE L FUNCTION                   C
1952 CCC      IERRO: ERROR FLAG                                     C
1953 CCC      * IF IERRO=0, COMPUTATION SUCCESSFUL.                 C
1954 CCC      * IF IERRO=1, COMPUTATION OUT OF RANGE.               C
1955 CCC      * IF IERRO=2, ARGUMENT AND/OR ORDER,  OUT OF RANGE.   C
1956 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1957 CCC  METHOD OF COMPUTATION:
1958 CCC  * XX<PNUA:   THE NON-OSCILLATING INTEGRAL REPRESENTATIONS
1959 CCC               FOR THE OSCILLATORY REGION ARE USED
1960 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1961 CCC            LOCAL VARIABLES:
1962 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1963 C CHI    : X*SINH(MU)-PNU*MU
1964 C CONTR1 : CONTRIBUTION OF THE SEMI-INFINITE INTEGRAL
1965 C          IN THE OSCILLATORY CASE (INCLUDING ADDITIONAL
1966 C          FACTORS: CONTR1=FOS1*FACTORS).
1967 C COSCHI : COS(CHI)
1968 C COSH2M : COSH(2*MU)
1969 C COSHM  : COSH(MU)
1970 C DF1    : FACTOR (DEPENDING ON IFAC)
1971 C DFF    : (1-EXP(-2*PI*PNU))*EXP(-ETA)
1972 C DMAIN  : PI*PNU*0.5
1973 C DMU    : SOLUTION MU OF COSH(MU)=PNU/X
1974 C DMU2   : 2*MU
1975 C DMU3   : 3*MU
1976 C DMU5   : 5*MU
1977 C DMU7   : 7*MU
1978 C DMU9   : 9*MU
1979 C DMUFAC : MU*COSH(MU)-SINH(MU)
1980 C DMUTAN : MU-TANH(MU)
1981 C EXPO   : EXP(PI*PNU*0.5)
1982 C FOS1   : CONTRIBUTION OF THE SEMI-INFINITE INTEGRAL
1983 C          IN THE OSCILLATORY CASE (LIA(X)).
1984 C FOSD1  : CONTRIBUTION OF THE SEMI-INFINITE INTEGRAL
1985 C          IN THE OSCILLATORY CASE (LIA'(X)).
1986 C OVER   : OVERFLOW NUMBER
1987 C PI     : 3.1415...
1988 C PINU   : PI*PNU
1989 C PNUA   : SAME AS PNU (INPUT VARIABLE)
1990 C SINCHI : SIN(CHI)
1991 C SINH2M : SINH(2*MU)
1992 C SINHM  : SINH(MU)
1993 C UNDER  : UNDERFLOW NUMBER
1994 C X      : ARGUMENT OF THE FUNCTION
1995 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1996        DOUBLE PRECISION CHI,CONTR1,COSCHI,COSH2M,COSHM,
1997      + D1MACH,DF1,DFF,DLIND,DLINF,DMAIN,DMU,DMU2,DMU3,
1998      + DMU5,DMU7,DMU9,DMUFAC,DMUTAN,EXPO,FOS1,FOSD1,OVER,
1999      + PI,PINU,PNU,PNUA,SINCHI,SINH2M,SINHM,UNDER,X,XX
2000        INTEGER IERRO,IFAC
2001        COMMON/ARGU/X,PNU
2002        COMMON/PAROS1/DMU,COSHM,SINHM,DMUFAC,COSCHI,SINCHI
2003        COMMON/PAROS2/COSH2M,SINH2M,DMAIN
2004        COMMON/PAROS3/DMUTAN
2005 CCC MACHINE DEPENDENT CONSTANT (OVERFLOW NUMBER)
2006        OVER=D1MACH(2)*1.D-8
2007 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
2008        UNDER=D1MACH(1)*1.D+8
2009        X=XX
2010        PNU=PNUA
2011        PI=ACOS(-1.D0)
2012        IERRO=0
2013 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2014 CC  THE USE OF INTEGRALS FOR THE LIA FUNCTION IS   C
2015 CC  RESTRICTED TO THE CASE X < A, WHICH IS THE     C
2016 CC  OSCILLATORY CASE                               C
2017 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2018        PI=ACOS(-1.D0)
2019        IF (IFAC.EQ.1) THEN
2020          IF ((PI*PNU*0.5D0).GE.LOG(OVER)) IERRO=1
2021        ENDIF
2022        IF (IERRO.EQ.0) THEN
2023          COSHM=PNU/X
2024          DMU=LOG(COSHM+SQRT((COSHM-1.D0)*(COSHM+1.D0)))
2025          DMU2=2.D0*DMU
2026          COSH2M=COSH(DMU2)
2027          SINHM=SINH(DMU)
2028          SINH2M=SINH(DMU2)
2029          IF (DMU.GT.0.1D0) THEN
2030            CHI=X*SINHM-PNU*DMU
2031            DMUFAC=DMU*COSHM-SINHM
2032          ELSE
2033            DMU2=DMU*DMU
2034            DMU3=DMU2*DMU
2035            DMU5=DMU3*DMU2
2036            DMU7=DMU5*DMU2
2037            DMU9=DMU7*DMU2
2038            CHI=-2.D0*X*(1.D0/6.D0*DMU3+1.D0/60.D0*DMU5+
2039      +      3.D0/5040.D0*DMU7+4.D0/362880.D0*DMU9)
2040            DMUFAC=DMU3/3.D0+DMU5/30.D0+DMU7/840.D0+DMU9/45360.D0
2041          ENDIF
2042          DMUTAN=DMU-TANH(DMU)
2043          COSCHI=COS(CHI)
2044          SINCHI=SIN(CHI)
2045          PINU=PI*PNU
2046          DMAIN=PINU*0.5D0
2047 CCCCCCCCCCCCCCCCCC
2048 CCC  INTEGRALS  CC
2049 CCCCCCCCCCCCCCCCCC
2050          CALL TRAPRL(1,FOS1)
2051 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2052 CCC THEN, THE LIA FUNCTION IS GIVEN BY ... CC
2053 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2054          IF (IFAC.EQ.1) THEN
2055            EXPO=EXP(DMAIN)
2056            DFF=EXPO/PI
2057            IF ((-2.D0*PINU).LE.LOG(UNDER)) THEN
2058              DF1=DFF*0.5D0
2059            ELSE
2060              DF1=DFF*(1.D0-EXP(-2.D0*PINU))*0.5D0
2061            ENDIF
2062            CONTR1=DF1*FOS1
2063          ELSE
2064            DFF=1.D0/PI
2065            IF ((-2.D0*PINU).LE.LOG(UNDER)) THEN
2066              DF1=DFF*0.5D0
2067            ELSE
2068              DF1=DFF*(1.D0-EXP(-2.D0*PINU))*0.5D0
2069            ENDIF
2070            CONTR1=DF1*FOS1
2071          ENDIF
2072          DLINF=CONTR1
2073 CCCCCCCCCCCCCCCCCCCCCCCCCC
2074 CCC CALCULATION OF LIA' CC
2075 CCCCCCCCCCCCCCCCCCCCCCCCCC
2076 CCCCCCCCCCCCCCCCCC
2077 CCC  INTEGRALS  CC
2078 CCCCCCCCCCCCCCCCCC
2079          CALL TRAPRL(10,FOSD1)
2080 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2081 CCC THEN, THE LIA FUNCTION IS GIVEN BY ... CC
2082 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2083          IF (IFAC.EQ.1) THEN
2084            EXPO=EXP(DMAIN)
2085            DFF=EXPO/PI
2086            IF ((-2.D0*PINU).LE.LOG(UNDER)) THEN
2087              DF1=DFF*0.5D0
2088            ELSE
2089              DF1=DFF*(1.D0-EXP(-2.D0*PINU))*0.5D0
2090            ENDIF
2091            CONTR1=DF1*FOSD1
2092          ELSE
2093            DFF=1.D0/PI
2094            IF ((-2.D0*PINU).LE.LOG(UNDER)) THEN
2095              DF1=DFF*0.5D0
2096            ELSE
2097              DF1=DFF*(1.D0-EXP(-2.D0*PINU))*0.5D0
2098            ENDIF
2099            CONTR1=DF1*FOSD1
2100          ENDIF
2101          DLIND=CONTR1
2102        ELSE
2103          DLINF=0.D0
2104          DLIND=0.D0
2105        ENDIF
2106        RETURN
2107        END
2108
2109        DOUBLE PRECISION FUNCTION FSTAL2(T)
2110 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2111 CCC  COMPUTATION OF THE INTEGRAND FOR THE
2112 CCC  INTEGRAL FOR THE L FUNCTION IN EQ.(37) OF
2113 CCC  THE REFERENCE:
2114 CCC  'COMPUTING SOLUTIONS OF THE MODIFIED BESSEL
2115 CCC   DIFFERENTIAL EQUATION FOR IMAGINARY ORDERS
2116 CCC   AND POSITIVE ARGUMENTS',
2117 CCC   A. GIL, J. SEGURA, N.M. TEMME
2118 CCC   ACM TRANS. MATH. SOFT. (2004)
2119 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2120 CCC        INPUT VARIABLE:
2121 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2122 C    T    : ARGUMENT OF THE FUNCTION
2123 C           (VARIABLE OF INTEGRATION).
2124 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2125 CCC        LOCAL VARIABLES:
2126 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2127 C ARGU    : (COSH(MU)*U-DMUFAC)/SINH(U)
2128 C ARGU2   : ARGU**2
2129 C COSCHI  : COSH(CHI), CHI=X*SINH(MU)-PNU*MU
2130 C COSH2M  : COSH(2*MU)
2131 C COSHM   : COSH(MU)
2132 C D1      : COSH(MU)-DMUFAC
2133 C DERI    : (COSH(MU)/SINH(U)-D1*COSH(U)/SINH(U)**2)
2134 C           /SQRT(1-ARGU2)
2135 C DJACO   : COSH(T)*EXP(S)/SQRT(1+EXP(S)**2)
2136 C DMAIN   : PI*PNU*0.5
2137 C DMU     : SOLUTION MU OF COSH(MU)=PNU/X
2138 C DMU2    : 2*MU
2139 C DMUFAC  : MU*COSH(MU)-SINH(MU)
2140 C EXPON   : EXP(-(PHIB(U)-PI*PNU*0.5))
2141 C F1      : SINH(U-MU)/(U-MU)
2142 C G1      : Z/6+Z3/120+Z5/5040+Z7/362880, Z=2*Y
2143 C PHIB(U) : X*COSH(U)*COS(SIGMA(U))+PNU*SIGMA(U),
2144 C           WHERE X=ARGUMENT OF THE FUNCTION,
2145 C           SIGMA(U)=ARCSIN((COSH(MU)*U-DMUFAC)/SINH(U)
2146 C RESTO   : SIN(CHI)-COS(CHI)*DERI
2147 C S       : SINH(T)
2148 C SINCHI  : SIN(CHI)
2149 C SINH2M  : SINH(2*MU)
2150 C SINHM   : SINH(MU)
2151 C SINHU   : SINH(U)
2152 C SINHU2  : 2*SINH(U)
2153 C U       : MU+LOG(X+SQRT(X**2+1))
2154 C UNDER   : UNDERFLOW NUMBER
2155 C X       : EXP(S)
2156 C Y       : U-MU
2157 C Z       : 2*Y
2158 C Z2      : Z**2
2159 C Z3      : Z**3
2160 C Z5      : Z**5
2161 C Z7      : Z**7
2162 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2163        DOUBLE PRECISION ARGU,ARGU2,COSCHI,COSH2M,COSHM,
2164      + D1,D1MACH,DERI,DJACO,DMAIN,DMU,DMU2,DMUFAC,DMUTAN,
2165      + EXPON,F1,G1,PHIB,RESTO,S,SINCHI,SINH2M,SINHM,
2166      + SINHU,SINHU2,T,U,UNDER,X,Y,Z,Z2,Z3,Z5,Z7
2167        COMMON/PAROS1/DMU,COSHM,SINHM,DMUFAC,COSCHI,SINCHI
2168        COMMON/PAROS2/COSH2M,SINH2M,DMAIN
2169        COMMON/PAROS3/DMUTAN
2170 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
2171        UNDER=D1MACH(1)*1.D+8
2172        S=SINH(T)
2173        X=EXP(S)
2174        U=DMUTAN+LOG(X+SQRT(X*X+1.D0))
2175        Y=U-DMU
2176        IF (ABS(Y).LE.0.1D0) THEN
2177          IF (ABS(Y).GT.UNDER) THEN
2178            F1=SINH(Y)/Y
2179          ELSE
2180            F1=1.D0
2181          ENDIF
2182          Z=Y*2.D0
2183          Z2=Z*Z
2184          Z3=Z2*Z
2185          Z5=Z3*Z2
2186          Z7=Z5*Z2
2187          G1=Z/6.D0+Z3/120.D0+Z5/5040.D0+Z7/362880.D0
2188          DMU2=2.D0*DMU
2189          DERI=SINH(U)/(F1-COSHM*COSH(U))*
2190      +      SQRT(COSH(DMU2)*F1*F1+2.D0*SINH(DMU2)*G1-
2191      +      COSHM*COSHM)
2192          DERI=1.D0/DERI
2193        ELSE
2194          D1=COSHM*U-DMUFAC
2195          SINHU=SINH(U)
2196          ARGU=D1/SINHU
2197          ARGU2=ARGU*ARGU
2198          SINHU2=SINHU*SINHU
2199          DERI=1.D0/SQRT(1.D0-ARGU2)*
2200      +    (COSHM/SINHU-D1*COSH(U)/SINHU2)
2201          IF (U.LT.DMU) DERI=-DERI
2202        ENDIF
2203        DJACO=COSH(T)*X/SQRT(1.D0+X*X)
2204        IF ((-(PHIB(U)-DMAIN)).LE.LOG(UNDER)) THEN
2205          EXPON=0.D0
2206          FSTAL2=0.D0
2207        ELSE
2208          EXPON=EXP(-(PHIB(U)-DMAIN))
2209          RESTO=(SINCHI-COSCHI*DERI)
2210          FSTAL2=EXPON*RESTO*DJACO
2211        ENDIF
2212        RETURN
2213        END
2214
2215        DOUBLE PRECISION FUNCTION FDTAL2(T)
2216 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2217 CCC  COMPUTATION OF THE INTEGRAND FOR THE
2218 CCC  INTEGRAL FOR THE DERIVATIVE OF THE L FUNCTION
2219 CCC  IN EQ.(37) OF THE REFERENCE:
2220 CCC  'COMPUTING SOLUTIONS OF THE MODIFIED BESSEL
2221 CCC   DIFFERENTIAL EQUATION FOR IMAGINARY ORDERS
2222 CCC   AND POSITIVE ARGUMENTS',
2223 CCC   A. GIL, J. SEGURA, N.M. TEMME
2224 CCC   ACM TRANS. MATH. SOFT. (2004)
2225 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2226 CCC        INPUT VARIABLE:
2227 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2228 C    T    : ARGUMENT OF THE FUNCTION
2229 C           (VARIABLE OF INTEGRATION).
2230 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2231 CCC         LOCAL VARIABLES:
2232 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2233 C ARGU    : (COSH(MU)*U-DMUFAC)/SINH(U)
2234 C ARGU2   : ARGU**2
2235 C COSCHI  : COSH(CHI), CHI=X*SINH(MU)-PNU*MU
2236 C COSH2M  : COSH(2*MU)
2237 C COSHM   : COSH(MU)
2238 C COSHU   : COSH(U)
2239 C COSS    : COS(SIGMA(U))
2240 C D1      : COSH(MU)-DMUFAC
2241 C DELTA   : -SIN(CHI)*COS(SIGMA(U))*COSH(U)
2242 C           +COS(CHI)*SIN(SIGMA(U))*SINH(U)
2243 C DERI    : (COSH(MU)/SINH(U)-D1*COSH(U)/SINH(U)**2)
2244 C           /SQRT(1-ARGU2)
2245 C DJACO   : COSH(T)*EXP(S)/SQRT(1+EXP(S)**2)
2246 C DMAIN   : PI*PNU*0.5
2247 C DMU     : SOLUTION MU OF COSH(MU)=PNU/X
2248 C DMU2    : 2*MU
2249 C DMUFAC  : MU*COSH(MU)-SINH(MU)
2250 C EXPON   : EXP(-(PHIB(U)-PI*PNU*0.5))
2251 C F1      : SINH(U-MU)/(U-MU)
2252 C G1      : Z/6+Z3/120+Z5/5040+Z7/362880, Z=2*Y
2253 C GAMMA   : SIN(CHI)*SIN(SIGMA(U))*SINH(U)+
2254 C           COS(CHI)*COS(SIGMA(U))*COSH(U)
2255 C PHIB(U) : X*COSH(U)*COS(SIGMA(U))+PNU*SIGMA(U),
2256 C           WHERE X=ARGUMENT OF THE FUNCTION,
2257 C           SIGMA(U)=ARCSIN((COSH(MU)*U-DMUFAC)/SINH(U)
2258 C RESTO   : -SIN(CHI)*COS(SIGMA(U))*COSH(U)
2259 C           +COS(CHI)*SIN(SIGMA(U))*SINH(U)+
2260 C           (SIN(CHI)*SIN(SIGMA(U))*SINH(U)+
2261 C            COS(CHI)*COS(SIGMA(U))*COSH(U))*DERI
2262 C S       : SINH(T)
2263 C SIGMA(U): ASIN((COSH(MU)*U-DMUFAC)/SINH(U))
2264 C SIGMAU  : SIGMA(U)
2265 C SINCHI  : SIN(CHI)
2266 C SINH2M  : SINH(2*MU)
2267 C SINHM   : SINH(MU)
2268 C SINHU   : SINH(U)
2269 C SINHU2  : 2*SINH(U)
2270 C SINS    : SIN(SIGMA(U))
2271 C U       : MU+LOG(X+SQRT(X**2+1))
2272 C UNDER   : UNDERFLOW NUMBER
2273 C X       : EXP(S)
2274 C Y       : U-MU
2275 C Z       : 2*Y
2276 C Z2      : Z**2
2277 C Z3      : Z**3
2278 C Z5      : Z**5
2279 C Z7      : Z**7
2280 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2281        DOUBLE PRECISION ARGU,ARGU2,COSCHI,COSH2M,COSHM,
2282      + COSHU,COSS,D1,D1MACH,DELTA,DERI,DJACO,DMAIN,DMU,
2283      + DMU2,DMUFAC,DMUTAN,EXPON,F1,G1,GAMMA,PHIB,RESTO,
2284      + S,SIGMA,SIGMAU,SINCHI,SINH2M,SINHM,SINHU,SINHU2,
2285      + SINS,T,U,UNDER,X,Y,Z,Z2,Z3,Z5,Z7
2286        COMMON/PAROS1/DMU,COSHM,SINHM,DMUFAC,COSCHI,SINCHI
2287        COMMON/PAROS2/COSH2M,SINH2M,DMAIN
2288        COMMON/PAROS3/DMUTAN
2289 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
2290        UNDER=D1MACH(1)*1.D+8
2291        S=SINH(T)
2292        X=EXP(S)
2293        U=DMUTAN+LOG(X+SQRT(X*X+1.D0))
2294        Y=U-DMU
2295        COSHU=COSH(U)
2296        SINHU=SINH(U)
2297        IF (ABS(Y).LE.0.1) THEN
2298          IF (ABS(Y).GT.UNDER) THEN
2299            F1=SINH(Y)/Y
2300          ELSE
2301            F1=1.D0
2302          ENDIF
2303          Z=Y*2.D0
2304          Z2=Z*Z
2305          Z3=Z2*Z
2306          Z5=Z3*Z2
2307          Z7=Z5*Z2
2308          G1=Z/6.D0+Z3/120.D0+Z5/5040.D0+Z7/362880.D0
2309          DMU2=2.D0*DMU
2310          DERI=SINHU/(F1-COSHM*COSHU)*SQRT(COSH(DMU2)*F1*F1+
2311      +   2.D0*SINH(DMU2)*G1-COSHM*COSHM)
2312          DERI=1.D0/DERI
2313        ELSE
2314          D1=COSHM*U-DMUFAC
2315          ARGU=D1/SINHU
2316          ARGU2=ARGU*ARGU
2317          SINHU2=SINHU*SINHU
2318          DERI=1.D0/SQRT(1.D0-ARGU2)*(COSHM/SINHU-
2319      +   D1*COSHU/SINHU2)
2320          IF (U.LT.DMU) DERI=-DERI
2321        ENDIF
2322        DJACO=COSH(T)*X/SQRT(1.D0+X*X)
2323        IF ((-(PHIB(U)-DMAIN)).LE.LOG(UNDER)) THEN
2324          EXPON=0.D0
2325          FDTAL2=0.D0
2326        ELSE
2327          EXPON=EXP(-(PHIB(U)-DMAIN))
2328          SIGMAU=SIGMA(U)
2329          SINS=SIN(SIGMAU)
2330          COSS=COS(SIGMAU)
2331          GAMMA=SINCHI*SINS*SINHU+COSCHI*COSS*COSHU
2332          DELTA=-SINCHI*COSS*COSHU+COSCHI*SINS*SINHU
2333          RESTO=DELTA+GAMMA*DERI
2334          FDTAL2=EXPON*RESTO*DJACO
2335        ENDIF
2336        RETURN
2337        END
2338
2339       SUBROUTINE TRAPRL(IC,TI)
2340 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2341 CCC IMPLEMENTATION OF AN ADAPTIVE TRAPEZOIDAL RULE
2342 CCC FOR COMPUTING THE INTEGRAL REPRESENTATIONS OF
2343 CCC THE FUNCTIONS.
2344 CCC
2345 CCC   INPUT VARIABLE:
2346 CCC      IC: DEPENDING ON THE VALUES OF IC,
2347 CCC          DIFFERENT INTEGRALS ARE COMPUTED:
2348 CCC          *IC=1,  L FUNCTION, OSCILLATORY REGION
2349 CCC          *IC=10, DERIVATIVE OF THE L FUNCTION,
2350 CCC                  OSCILLATORY REGION.
2351 CCC   OUTPUT VARIABLE:
2352 CCC          TI,  RESULT OF THE INTEGRAL
2353 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2354 CCC            LOCAL VARIABLES
2355 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2356 C A        :  LOWER INTEGRATION LIMIT
2357 C B        :  UPPER INTEGRATION LIMIT
2358 C DELTA    :  CALCULATES THE RELATIVE PRECISION
2359 C EPS      :  RELATIVE PRECISION PARAMETER USED IN
2360 C             THE CALCULATION
2361 C H        :  INTEGRATION STEP
2362 C PNU      :  ORDER OF THE FUNCTION
2363 C SUM      :  ACCUMULATES THE ELEMENTARY
2364 C             CONTRIBUTIONS
2365 C TIN      :  EVALUATED INTEGRAL
2366 C X        :  ARGUMENT OF THE FUNCTION
2367 C XAC      :  INTEGRATION ABCISSA
2368 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2369        DOUBLE PRECISION A,B,D1MACH,DELTA,EPS,FINTEL,
2370      + H,PNU,SUM,TI,TIN,X,XAC,Z
2371        INTEGER I,IC,IFI,N
2372        COMMON/ARGU/X,PNU
2373        EPS=D1MACH(3)
2374        IF (EPS.LT.1.D-14) EPS=1.D-14
2375        N=0
2376 CCCCC   INTEGRATION LIMITS: A,B
2377        Z=X/PNU
2378        IF ((Z.GT.0.999D0).AND.(Z.LT.1.001D0)) THEN
2379          A=-4.5
2380        ELSE
2381          A=-2.5D0
2382        ENDIF
2383        B=-A
2384        H=B-A
2385        TI=0.5D0*H*(FINTEL(IC,A)+FINTEL(IC,B))
2386        DELTA=1.D0+EPS
2387 11     N=N+1
2388        H=0.5D0*H
2389        IF (N.EQ.1) THEN
2390          IFI=1
2391        ELSE
2392          IFI=2*IFI
2393        ENDIF
2394        SUM=0.D0
2395        DO 3 I=1,IFI
2396          XAC=A+DBLE(2*I-1)*H
2397          SUM=SUM+FINTEL(IC,XAC)
2398  3     CONTINUE
2399        TIN=0.5D0*TI+H*SUM
2400        IF ((TIN.NE.0.D0).AND.(N.GT.4)) THEN
2401          DELTA=ABS(1.D0-TI/TIN)
2402        ENDIF
2403        TI=TIN
2404        IF ((DELTA.GT.EPS).AND.(N.LT.9)) GOTO 11
2405        RETURN
2406        END
2407
2408        DOUBLE PRECISION FUNCTION FINTEL(IC,T)
2409 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2410 CCC  AUXILIARY FUNCTION FOR THE SUBROUTINE TRAPRL.
2411 CCC  THIS FUNCTION CALLS THE DIFFERENT FUNCTIONS
2412 CCC  CONTAINING THE INTEGRANDS WHICH ARE INTEGRATED
2413 CCC  BY TRAPRL.
2414 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2415 CCC   INPUT:                                      C
2416 CCC      IC: INTEGER PARAMETER. ITS ADMISSIBLE    C
2417 CCC          VALUES ARE THE SAME AS IN THE        C
2418 CCC          SUBROUTINE TRAPRL.                   C
2419 CCC      T:  INTEGRATION ABSCISSA                 C
2420 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2421 CCC       LOCAL VARIABLES:
2422 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2423 C FDTAL2  : OSCILLATORY PART: TAU CONTRIBUTION ROUTINE
2424 C          (LIA'(X) FUNCTION). SEMI-INFINITE INTEGRAL.
2425 C               USED FOR LARGE PNU
2426 C FSTAL2  : OSCILLATORY PART: TAU CONTRIBUTION ROUTINE
2427 C            (LIA(X) FUNCTION). SEMI-INFINITE INTEGRAL.
2428 C                  USED FOR LARGE PNU
2429 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2430        DOUBLE PRECISION FDTAL2,FSTAL2,T
2431        INTEGER IC
2432        IF (IC.EQ.1) THEN
2433          FINTEL=FSTAL2(T)
2434        ENDIF
2435        IF (IC.EQ.10) THEN
2436          FINTEL=FDTAL2(T)
2437        ENDIF
2438        RETURN
2439        END
2440
2441        SUBROUTINE SERLIA(IFAC,X,PNU,PSER,PSERD,IERRO)
2442 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2443 CCC   CALCULATION OF POWER SERIES FOR L, L'
2444 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2445 CCC   INPUT VARIABLES:                              C
2446 CCC      X:   ARGUMENT OF THE FUNCTIONS             C
2447 CCC      PNU: ORDER OF THE FUNCTIONS                C
2448 CCC      IFAC:                                      C
2449 CCC               =1,  NON SCALED FUNCTIONS         C
2450 CCC               OTHERWISE,  SCALED FUNCTIONS      C
2451 CCC   OUTPUT VARIABLES:                             C
2452 CCC      PSER: L FUNCTION                           C
2453 CCC      PSERD: DERIVATIVE OF THE L FUNCTION        C
2454 CCC      IERRO: ERROR FLAG                          C
2455 CCC      * IF IERRO=0, COMPUTATION SUCCESSFUL.      C
2456 CCC      * IF IERRO=1, COMPUTATION OUT OF RANGE.    C
2457 CCC      * IF IERRO=2, ARGUMENT AND/OR ORDER,       C
2458 CCC                    OUT OF RANGE.                C
2459 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2460 CCC         LOCAL VARIABLES:
2461 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2462 C A       : ORDER OF THE FUNCTIONS
2463 C A2      : A**2
2464 C A2H     : A**2/2
2465 C A2N     : A**(2*N)
2466 C ACCP    : ACCUMULATES THE P COEFFICIENTS
2467 C ACCQ    : ACCUMULATES THE Q COEFFICIENTS
2468 C ARGU    : SIG(0)-A*LOG(X/2)
2469 C C       : X**(2*K)/K
2470 C COCI    : 1/(K**2+A**2)
2471 C COSTH   : COS(THET)
2472 C DELTAK  : ACCUMULATES THE CONTRIBUTION FOR THE
2473 C           LIA(X) FUNCTION
2474 C DELTKP  : ACCUMULATES THE CONTRIBUTION FOR THE
2475 C           LIA'(X) FUNCTION
2476 C DF1     : FACTOR (DEPENDING ON IFAC)
2477 C ETA0    : PARAMETER FOR THE CALCULATION OF THE
2478 C           COULOMB PHASE SHIFT
2479 C ETA02   : ETA0**2
2480 C F(K)    : SIN(PHI(A,K)-A*LOG(X/2))
2481 C           /(A**2*(1+A**2)...(K**2+A**2))**1/2,
2482 C           WHERE PHI(A,K)=PHASE(GAMMA(1+K+IA))
2483 C FDOMIN  : X*(COS(THET)+THET*SIN(THET))
2484 C K       : CONTRIBUTION TO THE LIA(X) FUNCTION
2485 C KP      : CONTRIBUTION TO THE LIA'(X) FUNCTION
2486 C OVER    : OVERFLOW NUMBER
2487 C P0      : PARAMETERS FOR THE CALCULATION OF THE COULOMB
2488 C           PHASE SHIFT
2489 C P1      : PARAMETERS FOR THE CALCULATION OF THE COULOMB
2490 C           PHASE SHIFT
2491 C P2      : PARAMETERS FOR THE CALCULATION OF THE COULOMB
2492 C           PHASE SHIFT
2493 C PI      : 3.1415...
2494 C PIA     : PI*A
2495 C PIA2    : 2*PI*A
2496 C PRECI   : RELATIVE PRECISION USED IN THE CALCULATION
2497 C Q0      : PARAMETERS FOR THE CALCULATION OF THE COULOMB
2498 C           PHASE SHIFT
2499 C Q1      : PARAMETERS FOR THE CALCULATION OF THE COULOMB
2500 C           PHASE SHIFT
2501 C Q2      : PARAMETERS FOR THE CALCULATION OF THE COULOMB
2502 C           PHASE SHIFT
2503 C R(K)    : F(K)*A/TAN(PHI(A,K)-A*LOG(X/2))
2504 C SIG0    : COULOMB PHASE SHIFT
2505 C SINTH   : SIN(THET)
2506 C THET    : ASIN(A/X)
2507 C UNDER   : UNDERFLOW NUMBER
2508 C X2      : X*X
2509 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2510        DOUBLE PRECISION A,A2,A2H,A2N,ACCP,ACCQ,ARGU,C,COCI,
2511      + COSTH,D1MACH,DDS,DEE,DELTAK,DELTKP,DF1,DSMALL,ETA0,
2512      + ETA02,EULER,F(0:500),FDOMIN,K,KP,OVER,P0(0:9),P1(0:9),
2513      + P2(0:6),PI,PIA,PIA2,PNU,PRECI,PSER,PSERD,Q0(0:9),
2514      + Q1(0:9),Q2(0:6),R(0:500),SIG0,SINTH,THET,UNDER,X,X2
2515        INTEGER IERRO,IFAC,L,M,N
2516        SAVE P0,Q0,P1,Q1,P2,Q2
2517        DATA P0/1.08871504904797411683D5,3.64707573081160914640D5,
2518      +        4.88801471582878013158D5,3.36275736298197324009D5,
2519      +        1.26899226277838479804D5,
2520      +        2.60795543527084582682D4,2.73352480554497990544D3,
2521      +        1.26447543569902963184D2,
2522      +        1.85446022125533909390D0,1.90716219990037648146D-3/
2523        DATA Q0/6.14884786346071135090D5,2.29801588515708014282D6,
2524      +        3.50310844128424021934D6,
2525      +        2.81194990286041080264D6,1.28236441994358406742D6,
2526      +        3.35209348711803753154D5,
2527      +        4.84319580247948701171D4,3.54877039006873206531D3,
2528      +        1.11207201299804390166D2,1.D0/
2529        DATA P1/-1.044100987526487618670D10,-1.508574107180079913696D10,
2530      +        -5.582652833355901160542D9,4.052529174369477275446D8,
2531      +         5.461712273118594275192D8,
2532      +         9.510404403068169395714D7,6.281126609997342119416D6,
2533      +         1.651178048950518520416D5,
2534      +         1.498824421329341285521D3,2.974686506595477984776D0/
2535        DATA Q1/1.808868161493543887787D10,3.869142051704700267785D10,
2536      +        3.003264575147162634046D10,1.075554651494601843525D10,
2537      +        1.901298501823290694245D9,
2538      +        1.665999832151229472632D8,6.952188089169487375936D6,
2539      +        1.253235080625688652718D5,7.904420414560291396996D2,1.D0/
2540        DATA P2/7.08638611024520906826D-3,-6.54026368947801591128D-2,
2541      +        2.92684143106158043933D-1,4.66821392319665609167D0,
2542      +        -3.43943790382690949054D0,
2543      +        -7.72786486869252994370D0,-9.88841771200290647461D-01/
2544        DATA Q2/-7.08638611024520908189D-3,6.59931690706339630254D-2,
2545      +        -2.98754421632058618922D-1,-4.63752355513412248006D0,
2546      +        3.79700454098863541593D0,7.06184065426336718524D0,1.D0/
2547 CCC MACHINE DEPENDENT CONSTANT (OVERFLOW NUMBER)
2548        OVER=D1MACH(2)*1.D-8
2549 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
2550        UNDER=D1MACH(1)*1.D+8
2551        PI=ACOS(-1.D0)
2552        IERRO=0
2553        IF (IFAC.EQ.1) THEN
2554          IF (X.GE.PNU) THEN
2555            SINTH=PNU/X
2556            THET=ASIN(SINTH)
2557            COSTH=COS(THET)
2558            FDOMIN=X*(COSTH+THET*SINTH)
2559            IF (FDOMIN.GE.LOG(OVER)) IERRO=1
2560          ELSE
2561            IF ((PI*PNU*0.5D0).GE.LOG(OVER)) IERRO=1
2562          ENDIF
2563        ENDIF
2564        IF (IERRO.EQ.0) THEN
2565          ETA0=1.8055470716051069198764D0
2566          EULER=0.577215664901532860606512D0
2567          PRECI=D1MACH(3)*10
2568 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2569 CCC COEFFICIENTS FOR THE CALCULATION OF THE COULOMB PHASE SHIFT
2570 CCC FROM CODY & HILLSTROM, MATH. COMPUT. 24(111) 1970
2571 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2572          A=PNU
2573          ETA02=ETA0*ETA0
2574          A2=A*A
2575          N=0
2576          ACCP=0.D0
2577          ACCQ=0.D0
2578          IF (A.LE.2.D0) THEN
2579  33        A2N=A2**N
2580            ACCP=ACCP+P0(N)*A2N
2581            ACCQ=ACCQ+Q0(N)*A2N
2582            N=N+1
2583            IF (N.LE.9) GOTO 33
2584            SIG0=A*(A2-ETA02)*ACCP/ACCQ
2585          ELSE
2586            IF ((A.GT.2.D0).AND.(A.LE.4.D0)) THEN
2587  44          A2N=A2**N
2588              ACCP=ACCP+P1(N)*A2N
2589              ACCQ=ACCQ+Q1(N)*A2N
2590              N=N+1
2591              IF (N.LE.9) GOTO 44
2592              SIG0=A*ACCP/ACCQ
2593            ELSE
2594  55          A2N=A2**N
2595              ACCP=ACCP+P2(N)/A2N
2596              ACCQ=ACCQ+Q2(N)/A2N
2597              N=N+1
2598              IF (N.LE.6) GOTO 55
2599              SIG0=ATAN(A)*0.5D0+A*(LOG(1.D0+A2)*0.5D0+ACCP/ACCQ)
2600            ENDIF
2601          ENDIF
2602 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2603 CCC EVALUATION OF F(0), R(0), R(1)
2604 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2605          PIA=PI*A
2606          A2H=A2*0.5D0
2607          ARGU=SIG0-A*LOG(X*0.5D0)
2608          R(0)=COS(ARGU)
2609          R(1)=1.D0/(1.D0+A2)*(R(0)-A*SIN(ARGU))
2610          IF (A.LT.UNDER) THEN
2611            F(0)=-(EULER+LOG(X*0.5D0))
2612          ELSE
2613            F(0)=1.D0/A*SIN(ARGU)
2614          ENDIF
2615          C=1.D0
2616          K=R(0)
2617          KP=A2H*F(0)
2618 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2619 CCC RECURSION FOR F(K), R(K), C(K)
2620 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2621          X2=0.25D0*X*X
2622          M=1
2623          COCI=1.D0/(M*M+A2)
2624          DELTAK=K*10
2625          DELTKP=KP*10
2626  66      F(M)=(M*F(M-1)+R(M-1))*COCI
2627          C=X2*C/M
2628          DELTAK=R(M)*C
2629          K=K+DELTAK
2630          DELTKP=(M*R(M)+A2H*F(M))*C
2631          KP=KP+DELTKP
2632          M=M+1
2633          COCI=1.D0/(M*M+A2)
2634          R(M)=((2.D0*M-1.D0)*R(M-1)-R(M-2))*COCI
2635          IF (M.LT.500) THEN
2636            IF ((ABS(DELTAK/K).GT.PRECI).OR.
2637      +       (ABS(DELTKP/KP).GT.PRECI)) GOTO 66
2638          ENDIF
2639          PIA2=2.D0*PIA
2640          IF (IFAC.EQ.1) THEN
2641            IF (-PIA2.LE.LOG(UNDER)) THEN
2642              DEE=0.D0
2643            ELSE
2644              DEE=EXP(-PIA2)
2645            ENDIF
2646            IF (A.LT.1.D-1) THEN
2647              IF (A.LT.UNDER) THEN
2648                DF1=1.D0
2649              ELSE
2650                L=0
2651                DDS=1.D0
2652                DSMALL=1.D0
2653  47            L=L+1
2654                DDS=-DDS*PIA2/(L+1.D0)
2655                DSMALL=DSMALL+DDS
2656                IF (ABS(DDS/DSMALL).GT.PRECI) GOTO 47
2657                DF1=EXP(PIA*0.5D0)*SQRT(DSMALL)
2658              ENDIF
2659            ELSE
2660              DF1=EXP(PIA*0.5D0)*SQRT((1.D0-DEE)/PIA2)
2661            ENDIF
2662            PSER=K*DF1
2663            PSERD=DF1*KP*2.D0/X
2664          ELSE
2665            IF (X.LT.A) THEN
2666              IF (-PIA2.LE.LOG(UNDER)) THEN
2667                DEE=0.D0
2668              ELSE
2669                DEE=EXP(-PIA2)
2670              ENDIF
2671              IF (A.LT.1.D-1) THEN
2672                IF (A.LT.UNDER) THEN
2673                  DF1=1.D0
2674                ELSE
2675                  L=0
2676                  DDS=1.D0
2677                  DSMALL=1.D0
2678  48              L=L+1
2679                  DDS=-DDS*PIA2/(L+1.D0)
2680                  DSMALL=DSMALL+DDS
2681                  IF (ABS(DDS/DSMALL).GT.PRECI) GOTO 48
2682                  DF1=SQRT(DSMALL)
2683                ENDIF
2684              ELSE
2685                DF1=SQRT((1.D0-DEE)/PIA2)
2686              ENDIF
2687            ELSE
2688              SINTH=A/X
2689              THET=ASIN(SINTH)
2690              COSTH=COS(THET)
2691              FDOMIN=X*(COSTH+THET*SINTH)
2692              IF (-PIA2.LE.LOG(UNDER)) THEN
2693                DEE=0.D0
2694              ELSE
2695                DEE=EXP(-PIA2)
2696              ENDIF
2697              IF (A.LT.1.D-1) THEN
2698                IF (A.LT.UNDER) THEN
2699                  DF1=1.D0
2700                ELSE
2701                  L=0
2702                  DDS=1.D0
2703                  DSMALL=1.D0
2704  49              L=L+1
2705                  DDS=-DDS*PIA2/(L+1.D0)
2706                  DSMALL=DSMALL+DDS
2707                  IF (ABS(DDS/DSMALL).GT.PRECI) GOTO 49
2708                  DF1=EXP(PIA*0.5D0-FDOMIN)*SQRT(DSMALL)
2709                ENDIF
2710              ELSE
2711                DF1=EXP(PIA*0.5D0-FDOMIN)*SQRT((1.D0-DEE)/PIA2)
2712              ENDIF
2713            ENDIF
2714            PSER=K*DF1
2715            PSERD=DF1*KP*2.D0/X
2716          ENDIF
2717        ELSE
2718          PSER=0.D0
2719          PSERD=0.D0
2720        ENDIF
2721        RETURN
2722        END
2723
2724        SUBROUTINE EXPLIA(IFAC,X,PNU,PEXP,PEXPP,IERRO)
2725 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2726 CCC CALCULATION OF THE ASYMPTOTIC EXPANSIONS FOR       C
2727 CCC THE L, L' FUNCTIONS                                C
2728 CCC                                                    C
2729 CCC   INPUT VARIABLES:                                 C
2730 CCC      X:   ARGUMENT OF THE FUNCTIONS                C
2731 CCC      PNU:   ORDER OF THE FUNCTIONS                 C
2732 CCC      IFAC:                                         C
2733 CCC               =1,  NON SCALED FUNCTIONS            C
2734 CCC               OTHERWISE,  SCALED FUNCTIONS         C
2735 CCC   OUTPUT VARIABLES:                                C
2736 CCC      PEXP: L FUNCTION                              C
2737 CCC      PEXPP: DERIVATIVE OF THE L FUNCTION           C
2738 CCC      IERRO: ERROR FLAG                             C
2739 CCC      * IF IERRO=0, COMPUTATION SUCCESSFUL.         C
2740 CCC      * IF IERRO=1, COMPUTATION OUT OF RANGE .      C
2741 CCC      * IF IERRO=2, ARGUMENT AND/OR ORDER,          C
2742 CCC                     OUT OF RANGE.                  C
2743 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2744 CCC              LOCAL VARIABLES:
2745 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2746 C A        : ORDER OF THE FUNCTIONS
2747 C A2       : A**2
2748 C COSTH    : COS(THET)
2749 C DELTAK   : ACCUMULATES THE CONTRIBUTION FOR THE
2750 C            LIA(X) FUNCTION
2751 C DELTKP   : ACCUMULATES THE CONTRIBUTION FOR THE
2752 C            LIA'(X) FUNCTION
2753 C DF       : FACTOR (DEPENDING ON IFAC)
2754 C FDOMIN   : X*(COS(THET)+THET*SIN(THET))
2755 C IND      : (-1)**M
2756 C K        : CONTRIBUTION TO THE LIA(X) FUNCTION
2757 C KP       : CONTRIBUTION TO THE LIA'(X) FUNCTION
2758 C OVER     : OVERFLOW NUMBER
2759 C PI       : 3.1415...
2760 C PIA      : PI*A
2761 C PRECI    : RELATIVE PRECISION USED IN THE CALCULATION
2762 C RAX(K+1) : -((K+1/2)**2+A**2)/(K+1)*RAX(K)
2763 C SINTH    : SIN(THET)
2764 C THET     : ASIN(PNU/X)
2765 C XI       : 1/X
2766 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2767        DOUBLE PRECISION A,A2,COSTH,D1MACH,DELTAK,
2768      + DELTKP,DF,FDOMIN,K,KP,OVER,PEXP,PEXPP,PI,PIA,
2769      + PNU,PRECI,RAX,SINTH,THET,X,XI
2770        INTEGER IERRO,IFAC,IND,M
2771 CCC MACHINE DEPENDENT CONSTANT (OVERFLOW NUMBER)
2772        OVER=D1MACH(2)*1.D-8
2773        PI=ACOS(-1.D0)
2774        IERRO=0
2775         IF (IFAC.EQ.1) THEN
2776          IF (X.GE.PNU) THEN
2777            SINTH=PNU/X
2778            THET=ASIN(SINTH)
2779            COSTH=COS(THET)
2780            FDOMIN=X*(COSTH+THET*SINTH)
2781            IF (FDOMIN.GE.LOG(OVER)) IERRO=1
2782          ELSE
2783            IF ((PI*PNU*0.5D0).GE.LOG(OVER)) IERRO=1
2784          ENDIF
2785        ENDIF
2786        IF (IERRO.EQ.0) THEN
2787          PRECI=D1MACH(3)*10
2788          A=PNU
2789          XI=1.D0/X
2790          A2=A*A
2791          RAX=1.D0
2792          M=1
2793          DELTAK=RAX
2794          K=DELTAK
2795          DELTKP=RAX*(0.5D0*XI-1.D0)
2796          KP=DELTKP
2797          IND=-1
2798  44      RAX=-((M-0.5D0)*(1.D0-0.5D0/M)+A2/M)*0.5D0/X*RAX
2799          DELTAK=RAX*IND
2800          K=K+DELTAK
2801          IF (K.GE.OVER) M=2000
2802          M=M+1
2803          IND=-IND
2804          IF ((ABS(DELTAK/K).GT.PRECI).AND.(M.LT.1000)) GOTO 44
2805          RAX=1.D0
2806          M=1
2807          IND=-1
2808  45      RAX=-((M-0.5D0)*(1.D0-0.5D0/M)+A2/M)*0.5D0/X*RAX
2809          DELTAK=RAX*IND
2810          DELTKP=DELTAK*(XI*(0.5D0+M)-1.D0)
2811          KP=KP+DELTKP
2812          IF (KP.GE.OVER) M=2000
2813          M=M+1
2814          IND=-IND
2815          IF ((ABS(DELTKP/KP).GT.PRECI).AND.(M.LT.1000)) GOTO 45
2816          IF (IFAC.EQ.1) THEN
2817            DF=SQRT(0.5D0/(PI*X))*EXP(X)
2818            K=K*DF
2819            KP=-KP*DF
2820          ELSE
2821            IF (X.LT.A) THEN
2822              PIA=PI*A
2823              DF=SQRT(0.5D0/(PI*X))*EXP(X-0.5D0*PIA)
2824              K=K*DF
2825              KP=-KP*DF
2826            ELSE
2827              SINTH=A/X
2828              THET=ASIN(SINTH)
2829              COSTH=COS(THET)
2830              FDOMIN=X*(COSTH+THET*SINTH)
2831              DF=SQRT(0.5D0/(PI*X))*EXP(X-FDOMIN)
2832              K=K*DF
2833              KP=-KP*DF
2834            ENDIF
2835          ENDIF
2836          PEXP=K
2837          PEXPP=KP
2838        ELSE
2839          PEXP=0.D0
2840          PEXPP=0.D0
2841        ENDIF
2842        RETURN
2843        END
2844
2845        SUBROUTINE AIEXLI(IFAC,X,A,DLAI,DLAID,IERROL)
2846 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2847 CCC AIRY-TYPE ASYMPTOTIC EXPANSION FOR THE K AND K'   C
2848 CCC FUNCTIONS                                         C
2849 CCC                                                   C
2850 CCC   INPUT VARIABLES:                                C
2851 CCC      X:   ARGUMENT OF THE FUNCTIONS               C
2852 CCC      A:   ORDER OF THE FUNCTIONS                  C
2853 CCC      IFAC:                                        C
2854 CCC               =1,  NON SCALED FUNCTIONS           C
2855 CCC               OTHERWISE,  SCALED FUNCTIONS        C
2856 CCC   OUTPUT VARIABLES:                               C
2857 CCC      DLAI: L FUNCTION                             C
2858 CCC      DLAID: DERIVATIVE OF THE L FUNCTION          C
2859 CCC      IERROL: ERROR FLAG                           C
2860 CCC      * IF IERROL=0, COMPUTATION SUCCESSFUL.       C
2861 CCC      * IF IERROL=1, COMPUTATION OUT OF RANGE.     C
2862 CCC      * IF IERROL=2, ARGUMENT AND/OR ORDER,        C
2863 CCC                     OUT OF RANGE.                 C
2864 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2865 CCC          LOCAL VARIABLES:
2866 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2867 C A0EX    : EXACT VALUE OF THE COEFFICIENT A_0
2868 C A13     : A**(1/3)
2869 C A2      : A**2
2870 C A23     : A**(2/3)
2871 C A2K     : A**(2*K)
2872 C AC      : COEFFICIENTS A_S(PSI) USED IN THE EXPANSIONS
2873 C APIHAL  : A*PI/2
2874 C ARG     : -A**(2/3)*PSI
2875 C AS      : ACCUMULATES THE CONTRIBUTION OF THE AC COEFFICIENTS
2876 C           (FOR THE LIA(X) FUNCTION)
2877 C ASP     : ACCUMULATES THE CONTRIBUTION OF THE AC COEFFICIENTS
2878 C           (FOR THE LIA'(X) FUNCTION)
2879 C B0EX    : EXACT VALUE OF THE COEFFICIENT B_0
2880 C B0PEX   : EXACT VALUE OF THE COEFFICIENT B'_0
2881 C BC      : COEFFICIENTS B_S(PSI) USED IN THE EXPANSIONS
2882 C BII     : IMAGINARY PART OF THE AIRY FUNCTION BI(Z)
2883 C BIR     : REAL PART OF THE AIRY FUNCTION BI(Z)
2884 C BS      : ACCUMULATES THE CONTRIBUTION OF THE BC COEFFICIENTS
2885 C           (FOR THE LIA(X) FUNCTION)
2886 C BSO     : ACCUMULATES THE OLD VALUES OF BS
2887 C BSP     : ACCUMULATES THE CONTRIBUTION OF THE BC COEFFICIENTS
2888 C           (FOR THE LIA'(X) FUNCTION)
2889 C BSPO    : ACCUMULATES THE OLD VALUES OF BSP
2890 C C0EX    : EXACT VALUE OF THE COEFFICIENT C_0
2891 C CHI     : ACCUMULATES THE CONTRIBUTION OF THE CHIN COEFFICIENTS
2892 C CHIEX   : EXACT VALUE OF THE FUNCTION CHI(PSI)
2893 C CHIN    : COEFFICIENTS FOR THE EXPANSION OF THE FUNCTION
2894 C           CHI(PSI)
2895 C COSTH   : COS(THET)
2896 C D0EX    : EXACT VALUE OF THE COEFFICIENT D_0
2897 C DBII    : IMAGINARY PART OF THE AIRY FUNCTION BI'(Z)
2898 C DBIR    : REAL PART OF THE AIRY FUNCTION BI'(Z)
2899 C DF      : 2**(1/3)
2900 C DZZ     : (ABS(1-Z**2))**(1/2)
2901 C ETA     : PSI/2**(1/3)
2902 C ETAJ    : ETA**J
2903 C ETAK    : ETA**K
2904 C ETAL    : ETA**L
2905 C EXPAM   : EXP(A*PI/2-FDOMIN)
2906 C EXPAPI  : EXP(A*PI/2)
2907 C F2      : (-)**K/A**(2*K)
2908 C F4      : (A**(1/3))**4
2909 C FAC     : FACTOR FOR THE LIA(X) FUNCTION
2910 C FACD    : FACTOR FOR THE LIA'(X) FUNCTION
2911 C FDOMIN  : X*(COS(THET)+THET*SIN(THET))
2912 C IERRO   : ERROR FLAG FOR THE AIRY FUNCTIONS
2913 C IFACA   :
2914 C           *IF IFACA=1, COMPUTATION OF UNSCALED AIRY BI(Z),BI'(Z)
2915 C           *IF IFACA=2, COMPUTATION OF SCALED AIRY BI(Z),BI'(Z)
2916 C IFUN    :
2917 C           *IF IFUN=1, COMPUTATION OF AIRY BI(Z)
2918 C           *IF IFUN=2, COMPUTATION OF AIRY BI'(Z)
2919 C OVER    : OVERFLOW NUMBER
2920 C PHI     : COEFFICIENTS FOR THE EXPANSION OF THE FUNCTION
2921 C           PHI(PSI)
2922 C PHIEX   : EXACT VALUE OF THE FUNCTION PHI(PSI)
2923 C PHIEX2  : PHI(PSI)**2
2924 C PHIS    : VALUE OF THE FUNCTION PHI(PSI)
2925 C PI      : 3.1415...
2926 C PIHALF  : PI/2
2927 C PSI     :
2928 C          *IF 0<=Z<=1, 2/3*PSI**3/2=LOG((1+SQRT(1-Z**2))/Z)-SQRT(1-Z**2)
2929 C          *IF Z>=1, 2/3*(-PSI)**3/2=SQRT(Z**2-1)-ARCCOS(1/Z)
2930 C PSI12   : SQRT(PSI)
2931 C PSI2    : PSI*PSI
2932 C PSI3    : PSI**3
2933 C SAS     : ACCUMULATES THE CONTRIBUTION OF A_S(PSI)
2934 C SBS     : ACCUMULATES THE CONTRIBUTION OF B_S(PSI)
2935 C SCS     : ACCUMULATES THE CONTRIBUTION OF C_S(PSI)
2936 C SDS     : ACCUMULATES THE CONTRIBUTION OF D_S(PSI)
2937 C SIG     : (-)**K
2938 C SINTH   : SIN(THET)
2939 C THET    : ASIN(A/X)
2940 C Y       : Z-1
2941 C Z       : X/A
2942 C Z2      : Z**2
2943 C Z21M    : 1-Z**2
2944 C ZDS     : SQRT(1-Z**2)
2945 C ZMASF   : EXPANSION IN TERMS OF (Z-1) FOR THE
2946 C           CALCULATION OF PSI(Z)
2947 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2948        DOUBLE PRECISION A,A0EX,A13,A2,A23,A2K,AC(0:5,0:20),
2949      + APIHAL,ARG,AS,ASP,B0EX,B0PEX,BC(0:5,0:20),BII,BIR,
2950      + BS,BSO,BSP,BSPO,C0EX,CHI,CHIEX,CHIN(0:26),COSTH,
2951      + COZMAS(0:12),D0EX,D1MACH,DBII,DBIR,DF,DLAI,DLAID
2952        DOUBLE PRECISION DZZ,ETA,ETAJ,ETAK,ETAL,EXPAM,EXPAPI,
2953      + F2,F4,FAC,FACD,FDOMIN,OVER,PHI(0:35),PHIEX,PHIEX2,
2954      + PHIS,PI,PIHALF,PSI,PSI12,PSI2,PSI3,SAS,SBS,SCS,SDS,
2955      + SIG,SINTH,THET,X,Y,Z,Z2,Z21M,ZDS,ZMASF
2956        INTEGER IERRO,IERROL,IFAC,IFACA,IFUN,INDA(0:5),INDB(0:5),
2957      + J,K,L
2958        SAVE PHI,CHIN,AC,BC
2959 CCCC VALUES OF THE COEFFICIENTS
2960        DATA PHI/1.D0,0.2D0,.25714285714285714286D-1,
2961      + -.56507936507936507937D-2,-.39367965367965367965D-2,
2962      + -.5209362066504924D-3,.3708541264731741D-3,
2963      + .2123827840293627D-3,.2150629788753145D-4,
2964      + -.2636904062981799D-4,-.1405469826493129D-4,
2965      + -.1149328899029441D-5,.1972641193938624D-5,
2966      + .1014324305593329D-5,.7034331100192200D-7,
2967      + -.1525044777392676D-6,-.7677866256900572D-7,
2968      + -.4669842638693018D-8,.1206673645965329D-7,
2969      + .5990877668092344D-8,.3269102150077715D-9,
2970      + -.97138350244428335095D-9,-.47745934295232233834D-9,
2971      + -.23750318642839155779D-10,.79244598109106655567D-10,
2972      + .38653584230817865528D-10,.17734105846426807873D-11,
2973      + -.65332110030864751956D-11,-.31673512094772575686D-11,
2974      + -.13524195177121030660D-12,.54324103217903338951D-12,
2975      + .26204918647967626464D-12,.10488134973664584898D-13,
2976      + -.45490420121539001435D-13,-.21851238232690420878D-13,
2977      + -.82456080260379042800D-15/
2978        DATA CHIN/.2D0,.1142857142857142857142857D-1,
2979      + -.2438095238095238095238095D-1,-.1003471449185734900020614D-1,
2980      + .8811404468547325690182833D-3,.2318339655809043564145605D-2,
2981      + .7794494413564441575646057D-3,-.1373952504077228744949558D-3,
2982      + -.2162301322540308393022684D-3,-.6585004634375583559042795D-4,
2983      + .1502851367097217578058824D-4,.1991904617871647675455262D-4,
2984      + .5731972706719818525615427D-5,-.1496427612747891044606885D-5,
2985      + -.1821231830428939670133992D-5,-.5054254875930882534538217D-6,
2986      + .1433283859497625931203415D-6,.1657681319078678321113361D-6,
2987      + .4485592642702540575627044D-7,-.1346138242826094098161839D-7,
2988      + -.1504601862773535117708677D-7,-.3995660198654805921651406D-8,
2989      + .1250124931952495738882300D-8,.1363187174221864073749532D-8,
2990      + .3567608459777506132029204D-9,-.1152749290139859167119863D-9,
2991      + -.1233547289799408517691696D-9/
2992        DATA COZMAS/.9428090415820634D0,-.4242640687119285D0,
2993      + .2904188565587606D0,-.2234261009999160D0,
2994      + .1821754153944745D0,-.1539625154198624D0,
2995      + .1333737583085222D0,-.1176617834148007D0,
2996      + .1052687194772381D0,-.9524025714638822D-1,
2997      + .8695738197073783D-1,-.8000034897653656D-1,
2998      + .7407421797273086D-1/
2999        DATA AC(0,0)/1.D0/
3000        DATA AC(1,0)/-.44444444444444444445D-2/
3001        DATA AC(1,1)/-.18441558441558441558D-2/
3002        DATA AC(1,2)/.11213675213675213675D-2/
3003        DATA AC(1,3)/.13457752124418791086D-2/
3004        DATA AC(1,4)/.0003880626562979504D0/
3005        DATA AC(1,5)/-.0001830686723781799D0/
3006        DATA AC(1,6)/-.0001995460887806733D0/
3007        DATA AC(1,7)/-.00005256191234041587D0/
3008        DATA AC(1,8)/.00002460619652459158D0/
3009        DATA AC(1,9)/.00002519246780924541D0/
3010        DATA AC(1,10)/.6333157376533242D-5/
3011        DATA AC(1,11)/-.2957485733830202D-5/
3012        DATA AC(1,12)/-.2925255920564838D-5/
3013        DATA AC(1,13)/-.7159702610502009D-6/
3014        DATA AC(1,14)/.3331510720390949D-6/
3015        DATA AC(1,15)/.3227670475692310D-6/
3016        DATA AC(1,16)/.7767729381664199D-7/
3017        DATA AC(1,17)/-.3600954237921120D-7/
3018        DATA AC(1,18)/-.3441724449034226D-7/
3019        DATA AC(1,19)/-.8188194356398772D-8/
3020        DATA AC(1,20)/.3783148485152038D-8/
3021        DATA AC(2,0)/.69373554135458897365D-3/
3022        DATA AC(2,1)/.46448349036584330703D-3/
3023        DATA AC(2,2)/-.42838130171535112460D-3/
3024        DATA AC(2,3)/-.0007026702868771135D0/
3025        DATA AC(2,4)/-.0002632580046778811D0/
3026        DATA AC(2,5)/.0001663853666288703D0/
3027        DATA AC(2,6)/.0002212087687818584D0/
3028        DATA AC(2,7)/.00007020345615329662D0/
3029        DATA AC(2,8)/-.00004000421782540614D0/
3030        DATA AC(2,9)/-.00004786324966453962D0/
3031        DATA AC(2,10)/-.00001394600741473631D0/
3032        DATA AC(2,11)/.7536186591273727D-5/
3033        DATA AC(2,12)/.8478502161067667D-5/
3034        DATA AC(2,13)/.2345355228453912D-5/
3035        DATA AC(2,14)/-.1225943294710883D-5/
3036        DATA AC(2,15)/-.1325082343401027D-5/
3037        DATA AC(2,16)/-.3539954776569997D-6/
3038        DATA AC(2,17)/.1808291719376674D-6/
3039        DATA AC(2,18)/.1900383515233655D-6/
3040        DATA AC(3,0)/-.35421197145774384076D-3/
3041        DATA AC(3,1)/-.31232252789031883276D-3/
3042        DATA AC(3,2)/.3716442237502298D-3/
3043        DATA AC(3,3)/.0007539269155977733D0/
3044        DATA AC(3,4)/.0003408300059444739D0/
3045        DATA AC(3,5)/-.0002634968172069594D0/
3046        DATA AC(3,6)/-.0004089275726648432D0/
3047        DATA AC(3,7)/-.0001501108759563460D0/
3048        DATA AC(3,8)/.00009964015205538056D0/
3049        DATA AC(3,9)/.0001352492955751283D0/
3050        DATA AC(3,10)/.00004443117087272903D0/
3051        DATA AC(3,11)/-.00002713205071914117D0/
3052        DATA AC(3,12)/-.00003396796969771860D0/
3053        DATA AC(3,13)/-.00001040708865273043D0/
3054        DATA AC(3,14)/.6024639065414414D-5/
3055        DATA AC(3,15)/.7143919607846883D-5/
3056        DATA AC(4,0)/.378194199201773D-3/
3057        DATA AC(4,1)/.000404943905523634D0/
3058        DATA AC(4,2)/-.000579130526946453D0/
3059        DATA AC(4,3)/-.00138017901171011D0/
3060        DATA AC(4,4)/-.000722520056780091D0/
3061        DATA AC(4,5)/.000651265924036825D0/
3062        DATA AC(4,6)/.00114674563328389D0/
3063        DATA AC(4,7)/.000474423189340405D0/
3064        DATA AC(4,8)/-.000356495172735468D0/
3065        DATA AC(4,9)/-.000538157791035111D0/
3066        DATA AC(4,10)/-.000195687390661225D0/
3067        DATA AC(4,11)/.000132563525210293D0/
3068        DATA AC(4,12)/.000181949256267291D0/
3069        DATA AC(5,0)/-.69114139728829416760D-3/
3070        DATA AC(5,1)/-.00085995326611774383285D0/
3071        DATA AC(5,2)/.0014202335568143511489D0/
3072        DATA AC(5,3)/.0038535426995603052443D0/
3073        DATA AC(5,4)/.0022752811642901374595D0/
3074        DATA AC(5,5)/-.0023219572034556988366D0/
3075        DATA AC(5,6)/-.0045478643394434635622D0/
3076        DATA AC(5,7)/-.0020824431758272449919D0/
3077        DATA AC(5,8)/.0017370443573195808719D0/
3078        DATA BC(0,0)/.14285714285714285714D-1/
3079        DATA BC(0,1)/.88888888888888888889D-2/
3080        DATA BC(0,2)/.20482374768089053803D-2/
3081        DATA BC(0,3)/-.57826617826617826618D-3/
3082        DATA BC(0,4)/-.60412089799844901886D-3/
3083        DATA BC(0,5)/-.0001472685745626922D0/
3084        DATA BC(0,6)/.00005324102148009784D0/
3085        DATA BC(0,7)/.00005206561006583416D0/
3086        DATA BC(0,8)/.00001233115050894939D0/
3087        DATA BC(0,9)/-.4905932728531366D-5/
3088        DATA BC(0,10)/-.4632230987136350D-5/
3089        DATA BC(0,11)/-.1077174523455235D-5/
3090        DATA BC(0,12)/.4475963978932822D-6/
3091        DATA BC(0,13)/.4152586188464624D-6/
3092        DATA BC(0,14)/.9555819293589234D-7/
3093        DATA BC(0,15)/-.4060599208403059D-7/
3094        DATA BC(0,16)/-.3731367187988482D-7/
3095        DATA BC(0,17)/-.8532670645553778D-8/
3096        DATA BC(0,18)/.3673017245573624D-8/
3097        DATA BC(0,19)/.3355960460784536D-8/
3098        DATA BC(0,20)/.7643107095110475D-9/
3099        DATA BC(1,0)/-.11848595848595848596D-2/
3100        DATA BC(1,1)/-.13940630797773654917D-2/
3101        DATA BC(1,2)/-.48141005586383737645D-3/
3102        DATA BC(1,3)/.26841705366016142958D-3/
3103        DATA BC(1,4)/.0003419706982709903D0/
3104        DATA BC(1,5)/.0001034548234902078D0/
3105        DATA BC(1,6)/-.00005418191982095504D0/
3106        DATA BC(1,7)/-.00006202184830690167D0/
3107        DATA BC(1,8)/-.00001724885886056087D0/
3108        DATA BC(1,9)/.8744675992887053D-5/
3109        DATA BC(1,10)/.9420684216180929D-5/
3110        DATA BC(1,11)/.2494922112085850D-5/
3111        DATA BC(1,12)/-.1238458608836357D-5/
3112        DATA BC(1,13)/-.1285461713809769D-5/
3113        DATA BC(1,14)/-.3299710862537507D-6/
3114        DATA BC(1,15)/.1613441105788315D-6/
3115        DATA BC(1,16)/.1633623194402374D-6/
3116        DATA BC(1,17)/.4104252949605779D-7/
3117        DATA BC(1,18)/-.1984317042326989D-7/
3118        DATA BC(1,19)/-.1973948142769706D-7/
3119        DATA BC(1,20)/-.4882194808588752D-8/
3120        DATA BC(2,0)/.43829180944898810994D-3/
3121        DATA BC(2,1)/.71104865116708668943D-3/
3122        DATA BC(2,2)/.31858383945387580576D-3/
3123        DATA BC(2,3)/-.0002404809426804458D0/
3124        DATA BC(2,4)/-.0003722966038621536D0/
3125        DATA BC(2,5)/-.0001352752059595618D0/
3126        DATA BC(2,6)/.00008691694372704142D0/
3127        DATA BC(2,7)/.0001158750753591377D0/
3128        DATA BC(2,8)/.00003724965927846447D0/
3129        DATA BC(2,9)/-.00002198334949606935D0/
3130        DATA BC(2,10)/-.00002686449633870452D0/
3131        DATA BC(2,11)/-.8023061612032524D-5/
3132        DATA BC(2,12)/.4494756592180126D-5/
3133        DATA BC(2,13)/.5193504763856015D-5/
3134        DATA BC(2,14)/.1477156191529617D-5/
3135        DATA BC(2,15)/-.7988793826096784D-6/
3136        DATA BC(3,0)/-.37670439477105454219D-3/
3137        DATA BC(3,1)/-.75856271658798642365D-3/
3138        DATA BC(3,2)/-.0004103253968775048D0/
3139        DATA BC(3,3)/.0003791263310429010D0/
3140        DATA BC(3,4)/.0006850981673903450D0/
3141        DATA BC(3,5)/.0002878310571932216D0/
3142        DATA BC(3,6)/-.0002157010636115705D0/
3143        DATA BC(3,7)/-.0003260863991373500D0/
3144        DATA BC(3,8)/-.0001181317008748678D0/
3145        DATA BC(3,9)/.00007887526841582582D0/
3146        DATA BC(3,10)/.0001072081833420685D0/
3147        DATA BC(3,11)/.00003544595251288735D0/
3148        DATA BC(3,12)/-.00002201447920733824D0/
3149        DATA BC(3,13)/-.00002789336359620813D0/
3150        DATA BC(4,0)/.00058453330122076187255D0/
3151        DATA BC(4,1)/.0013854690422372401251D0/
3152        DATA BC(4,2)/.00086830374184946900245D0/
3153        DATA BC(4,3)/-.00093502904801345951693D0/
3154        DATA BC(4,4)/-.0019175486005525492059D0/
3155        DATA BC(4,5)/-.00090795047113308137941D0/
3156        DATA BC(4,6)/.00077050429806392235104D0/
3157        DATA BC(4,7)/.0012953100128255983488D0/
3158        DATA BC(4,8)/.00051933869471899550762D0/
3159        DATA BC(4,9)/-.00038482631948759834653D0/
3160        DATA BC(4,10)/-.00057335393099012476502D0/
3161        DATA BC(5,0)/-.0014301070053470410656D0/
3162        DATA BC(5,1)/-.0038637811942002539408D0/
3163        DATA BC(5,2)/-.0027322816261168328245D0/
3164        DATA BC(5,3)/.0033294980346743452748D0/
3165        DATA BC(5,4)/.0075972237660887795911D0/
3166        DATA BC(5,5)/.0039816655673062060620D0/
3167        DATA BC(5,6)/-.0037510180460986006595D0/
3168        DATA INDA/0,20,18,15,12,8/
3169        DATA INDB/20,20,15,13,10,6/
3170 CCCC CONSTANTS CCCCCCCCCCCCCC
3171        PI=ACOS(-1.D0)
3172 CCC MACHINE DEPENDENT CONSTANT (OVERFLOW NUMBER)
3173        OVER=D1MACH(2)*1.D-8
3174        IERROL=0
3175        IF (IFAC.EQ.1) THEN
3176          IF (X.GE.A) THEN
3177            SINTH=A/X
3178            THET=ASIN(SINTH)
3179            COSTH=COS(THET)
3180            FDOMIN=X*(COSTH+THET*SINTH)
3181            IF (FDOMIN.GE.LOG(OVER)) IERROL=1
3182          ELSE
3183            IF ((PI*A*0.5D0).GE.LOG(OVER)) IERROL=1
3184          ENDIF
3185        ENDIF
3186        IF (IERROL.EQ.0) THEN
3187          PIHALF=PI*0.5D0
3188          DF=2.D0**(1.D0/3.D0)
3189 CCCC VARIABLES CCCCCCCCCCCCCCC
3190          Z=X/A
3191          A2=A*A
3192          A13=A**(1.D0/3.D0)
3193          A23=A13*A13
3194          IF (IFAC.EQ.1) THEN
3195            APIHAL=A*PIHALF
3196            EXPAPI=EXP(APIHAL)
3197            FAC=EXPAPI*0.5D0/A13
3198            FACD=EXPAPI/Z/A23
3199          ELSE
3200            IF (X.LT.A) THEN
3201              FAC=0.5D0/A13
3202              FACD=1.D0/Z/A23
3203            ELSE
3204              SINTH=A/X
3205              THET=ASIN(SINTH)
3206              COSTH=COS(THET)
3207              FDOMIN=X*(COSTH+THET*SINTH)
3208              EXPAM=EXP(A*PIHALF-FDOMIN)
3209              FAC=EXPAM*0.5D0/A13
3210              FACD=EXPAM/Z/A23
3211            ENDIF
3212          ENDIF
3213          F4=A13**4
3214          IF (Z.LE.0.9D0) THEN
3215            ZDS=SQRT((1.D0-Z)*(1.D0+Z))
3216            PSI=(1.5D0*(LOG((1.D0+ZDS)/Z)-ZDS))**(2.D0/3.D0)
3217          ELSEIF (Z.GT.1.1D0) THEN
3218            ZDS=SQRT((Z-1.D0)*(1.D0+Z))
3219            PSI=-(1.5D0*(ZDS-ACOS(1.D0/Z)))**(2.D0/3.D0)
3220          ELSE
3221            Y=Z-1.D0
3222            ZMASF=COZMAS(12)
3223            DO 8 K=0,11
3224              J=11-K
3225              ZMASF=COZMAS(J)+Y*ZMASF
3226  8         CONTINUE
3227            ZMASF=ABS(Y)**(1.5D0)*ZMASF
3228            IF (Z.LT.1.D0) THEN
3229              PSI=(1.5D0*ZMASF)**(2.D0/3.D0)
3230            ELSE
3231              PSI=-(1.5D0*ZMASF)**(2.D0/3.D0)
3232            ENDIF
3233          ENDIF
3234          ETA=PSI/DF
3235          ARG=-A23*PSI
3236          PSI2=PSI*PSI
3237          PSI3=PSI2*PSI
3238          IF ((Z.GT.0.8D0).AND.(Z.LT.1.2D0)) THEN
3239            PHIS=0.D0
3240            CHI=0.D0
3241            SAS=0.D0
3242            SBS=0.D0
3243            SDS=0.D0
3244            SCS=1.D0
3245            BS=0.D0
3246            BSPO=0.D0
3247            DO 41 L=0,20
3248              IF (L.EQ.0) THEN
3249                ETAL=1.D0
3250              ELSE
3251                ETAL=ETAL*ETA
3252              ENDIF
3253              CHI=CHI+CHIN(L)*ETAL
3254  41        CONTINUE
3255            CHI=CHI/DF
3256            DO 10 K=0,20
3257              BSO=BS
3258              BSPO=BSP
3259              AS=0.D0
3260              BS=0.D0
3261              ASP=0.D0
3262              BSP=0.D0
3263              IF (K.EQ.0) THEN
3264                ETAK=1.D0
3265                A2K=1.D0
3266                SIG=1.D0
3267              ELSE
3268                ETAK=ETAK*ETA
3269                A2K=A2K*A2
3270                SIG=-1.D0*SIG
3271              ENDIF
3272              PHIS=PHIS+PHI(K)*ETAK
3273              F2=SIG/A2K
3274              IF (K.LE.5) THEN
3275                DO 20 J=0,20
3276                  IF (J.EQ.0) THEN
3277                    ETAJ=1.D0
3278                  ELSE
3279                    ETAJ=ETAJ*ETA
3280                  ENDIF
3281                  IF (J.LE.INDA(K)) THEN
3282                    AS=AS+AC(K,J)*ETAJ
3283                  ENDIF
3284                  IF (J.LE.INDB(K)) THEN
3285                    BS=BS+BC(K,J)*ETAJ
3286                  ENDIF
3287                  IF ((J+1).LE.INDA(K)) THEN
3288                    ASP=ASP+(J+1)*AC(K,J+1)*ETAJ
3289                  ENDIF
3290                  IF ((J+1).LE.INDB(K)) THEN
3291                    BSP=BSP+(J+1)*BC(K,J+1)*ETAJ
3292                  ENDIF
3293  20            CONTINUE
3294                ASP=1.D0/DF*ASP
3295                BS=BS*DF
3296                SAS=SAS+AS*F2
3297                SBS=SBS+BS*F2/DF
3298                SDS=SDS-(CHI*AS+ASP+PSI*BS)*F2
3299              ENDIF
3300              IF ((K.GT.0).AND.(K.LE.6)) THEN
3301                SCS=SCS+(AS+CHI*BSO+BSPO)*F2
3302              ENDIF
3303  10        CONTINUE
3304            PHIS=DF*PHIS
3305            SBS=DF*SBS
3306          ELSE
3307 CCCC EXACT VALUES CCCCCCCCCCCCCCCCCCCCC
3308            Z2=Z*Z
3309            Z21M=1.D0-Z2
3310            PHIEX=(4.D0*PSI/Z21M)**0.25D0
3311            PHIEX2=PHIEX*PHIEX
3312            A0EX=1.D0
3313            B0EX=-5.D0/48.D0/(PSI*PSI)+PHIEX2/48.D0/PSI*
3314      +      (5.D0/Z21M-3.D0)
3315            CHIEX=0.25D0/PSI*(1.D0-Z2*PHIEX2**3*0.25D0)
3316            CHI=CHIEX
3317            D0EX=-(7.D0/48.D0/PSI+PHIEX2/48.D0*
3318      +      (9.D0-7.D0/Z21M))
3319            IF (PSI.GT.0.D0) THEN
3320              PSI12=SQRT(PSI)
3321              DZZ=SQRT(Z21M)
3322            ELSE
3323              PSI12=SQRT(-PSI)
3324              DZZ=SQRT(ABS(Z21M))
3325            ENDIF
3326            B0PEX=5.D0/24.D0/PSI3+PHIEX2/48.D0*((2.D0*CHI*PSI-1.D0)/
3327      +       PSI2*(5.D0/Z21M-3.D0)-10.D0*Z2*PSI12/
3328      +       Z21M**2/DZZ/PSI)
3329            C0EX=1.D0
3330            SAS=A0EX
3331            SBS=B0EX
3332            SCS=C0EX
3333            SDS=D0EX
3334            BS=0.D0
3335            BSP=0.D0
3336            A2K=1.D0
3337            SIG=1.D0
3338            DO 30 K=1,6
3339              BSO=BS
3340              BSPO=BSP
3341              AS=0.D0
3342              BS=0.D0
3343              ASP=0.D0
3344              BSP=0.D0
3345              A2K=A2K*A2
3346              SIG=-1.D0*SIG
3347              F2=SIG/A2K
3348              IF (K.LE.5) THEN
3349                DO 40 J=0,20
3350                  IF (J.EQ.0) THEN
3351                    ETAJ=1.D0
3352                  ELSE
3353                    ETAJ=ETAJ*ETA
3354                  ENDIF
3355                  IF (J.LE.INDA(K)) THEN
3356                    AS=AS+AC(K,J)*ETAJ
3357                  ENDIF
3358                  IF (J.LE.INDB(K)) THEN
3359                    BS=BS+BC(K,J)*ETAJ
3360                  ENDIF
3361                  IF ((J+1).LE.INDA(K)) THEN
3362                    ASP=ASP+(J+1)*AC(K,J+1)*ETAJ
3363                  ENDIF
3364                  IF ((J+1).LE.INDB(K)) THEN
3365                    BSP=BSP+(J+1)*BC(K,J+1)*ETAJ
3366                  ENDIF
3367  40            CONTINUE
3368                BS=DF*BS
3369                ASP=ASP/DF
3370                SAS=SAS+AS*F2
3371                SBS=SBS+BS*F2
3372                SDS=SDS-(CHI*AS+ASP+PSI*BS)*F2
3373              ENDIF
3374              IF (K.GE.2) THEN
3375                SCS=SCS+(AS+CHI*BSO+BSPO)*F2
3376              ELSE
3377                SCS=SCS+(AS+CHI*B0EX+B0PEX)*F2
3378              ENDIF
3379  30        CONTINUE
3380            PHIS=PHIEX
3381          ENDIF
3382 CCCCC CALL THE AIRY ROUTINE CCCCCCCCC
3383          IFUN=1
3384          IFACA=1
3385          CALL BIZ(IFUN,IFACA,ARG,0.D0,BIR,BII,IERRO)
3386          IFUN=2
3387          IFACA=1
3388          CALL BIZ(IFUN,IFACA,ARG,0.D0,DBIR,DBII,IERRO)
3389          DLAI=FAC*PHIS*(BIR*SAS+DBIR*SBS/F4)
3390          DLAID=FACD/PHIS*(DBIR*SCS+BIR*SDS/A23)
3391        ELSE
3392          DLAI=0.D0
3393          DLAID=0.D0
3394        ENDIF
3395        RETURN
3396        END
3397
3398
3399        DOUBLE PRECISION FUNCTION V(U)
3400 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3401 CCC AUXILIARY FUNCTION:
3402 CCC THIS FUNCTION COMPUTES THE LAST EXPRESSION
3403 CCC IN EQ.(21) OF THE REFERENCE
3404 CCC  'COMPUTING SOLUTIONS OF THE MODIFIED BESSEL
3405 CCC   DIFFERENTIAL EQUATION FOR IMAGINARY ORDERS
3406 CCC   AND POSITIVE ARGUMENTS',
3407 CCC   A. GIL, J. SEGURA, N.M. TEMME
3408 CCC   ACM TRANS. MATH. SOFT. (2004)
3409 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3410 CCC        INPUT VARIABLE:
3411 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3412 C    U    : ARGUMENT OF THE FUNCTION
3413 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3414 CCC         LOCAL VARIABLES:
3415 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3416 C COSTH    : COS(THET)
3417 C FDOMIN   : X*(COS(THET)+THET*SIN(THET))
3418 C SINTH    : SIN(THET)
3419 C THET     : ASIN(PNU/X)
3420 C UNDER    : UNDERFLOW NUMBER
3421 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3422        DOUBLE PRECISION COSTH,D1MACH,FDOMIN,SINTH,
3423      + THET,U,UNDER
3424        COMMON/PARMON/THET,SINTH,COSTH,FDOMIN
3425 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
3426        UNDER=D1MACH(1)*1.D+6
3427        IF (ABS(U).LT.UNDER) THEN
3428          V=ASIN(SINTH)
3429        ELSE
3430          V=ASIN(U/SINH(U)*SINTH)
3431        ENDIF
3432        RETURN
3433        END
3434
3435        DOUBLE PRECISION FUNCTION PHIR(U)
3436 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3437 CCC AUXILIARY FUNCTION:
3438 CCC THIS FUNCTION COMPUTES THE EXPONENT IN THE
3439 CCC INTEGRAND IN EQ.(20) OF THE REFERENCE
3440 CCC  'COMPUTING SOLUTIONS OF THE MODIFIED BESSEL
3441 CCC   DIFFERENTIAL EQUATION FOR IMAGINARY ORDERS
3442 CCC   AND POSITIVE ARGUMENTS',
3443 CCC   A. GIL, J. SEGURA, N.M. TEMME
3444 CCC   ACM TRANS. MATH. SOFT. (2004)
3445 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3446 CCC        INPUT VARIABLE:
3447 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3448 C    U    : ARGUMENT OF THE FUNCTION
3449 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3450 CCC       LOCAL VARIABLES:
3451 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3452 C COSTH   : COS(THET)
3453 C COSVU   : COS(V(U))
3454 C FDOMIN  : X*(COS(THET)+THET*SIN(THET))
3455 C FU1     : 2*SINH(U/2)**2
3456 C FU2     : -2*SIN(FU3/2)*SIN(0.5*(V(U)+THET))
3457 C FU3     : ASIN(-SIN(THET)/(COS(THET)*U/SINH(U)
3458 C           +COS(V(U))*(SINH(U)+U)*FUAC/(SINH(U)*SINH(U)))
3459 C FUAC    : U**3/6+U**5/120+U**7/5040
3460 C OVER    : OVERFLOW NUMBER
3461 C PNU     : ORDER OF THE FUNCTION
3462 C SINHU   : SINH(U)
3463 C SINHUH  : SINH(U/2)
3464 C SINTH   : SIN(THET)
3465 C THET    : ASIN(PNU/X)
3466 C U2      : U**2
3467 C U3      : U**3
3468 C U5      : U**5
3469 C U7      : U**7
3470 C UH      : U/2
3471 C UNDER   : UNDERFLOW NUMBER
3472 C V(U)    : ASIN(U/SINH(U)*SIN(THET))
3473 C VU      : V(U)
3474 C X       : ARGUMENT OF THE FUNCTIONS
3475 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3476        DOUBLE PRECISION COSTH,COSVU,D1MACH,FDOMIN,FU1,FU2,
3477      + FU3,FUAC,OVER,PNU,SINHU,SINHUH,SINTH,THET,U,U2,
3478      + U3,U5,U7,UH,UNDER,V,VU,X
3479        COMMON/ARGU/X,PNU
3480        COMMON/PARMON/THET,SINTH,COSTH,FDOMIN
3481 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
3482        UNDER=D1MACH(1)*1.D+6
3483 CCC MACHINE DEPENDENT CONSTANT (OVERFLOW NUMBER)
3484        OVER=D1MACH(2)*1.D-6
3485        IF (U.LT.200.D0) THEN
3486          IF (ABS(U).LT.0.1D0) THEN
3487            IF (ABS(U).LT.UNDER) THEN
3488              PHIR=0.D0
3489            ELSE
3490              UH=U*0.5D0
3491              SINHUH=SINH(UH)
3492              FU1=2.D0*SINHUH*SINHUH
3493              U2=U*U
3494              U3=U2*U
3495              U5=U3*U2
3496              U7=U5*U2
3497              FUAC=U3/6.D0+U5/120.D0+U7/5040.D0
3498              SINHU=SINH(U)
3499              VU=V(U)
3500              COSVU=COS(VU)
3501              FU3=ASIN(-SINTH/(COSTH*U/SINHU+COSVU)*
3502      +          (SINHU+U)*FUAC/(SINHU*SINHU))
3503              FU2=-2.D0*SIN(FU3*0.5D0)*SIN(0.5D0*(VU+THET))
3504              PHIR=X*(FU1*COSVU+FU2+SINTH*FU3)
3505            ENDIF
3506          ELSE
3507            VU=V(U)
3508            COSVU=COS(VU)
3509            PHIR=X*COSH(U)*COSVU+PNU*VU-FDOMIN
3510          ENDIF
3511        ELSE
3512          PHIR=OVER
3513        ENDIF
3514        RETURN
3515        END
3516
3517        DOUBLE PRECISION FUNCTION SIGMA(U)
3518 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3519 CCC AUXILIARY FUNCTION:
3520 CCC THIS FUNCTION COMPUTES THE FUNCTION
3521 CCC SIGMA FROM EQ.(31) OF THE REFERENCE
3522 CCC  'COMPUTING SOLUTIONS OF THE MODIFIED BESSEL
3523 CCC   DIFFERENTIAL EQUATION FOR IMAGINARY ORDERS
3524 CCC   AND POSITIVE ARGUMENTS',
3525 CCC   A. GIL, J. SEGURA, N.M. TEMME
3526 CCC   ACM TRANS. MATH. SOFT. (2004)
3527 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3528 CCC        INPUT VARIABLE:
3529 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3530 C    U    : ARGUMENT OF THE FUNCTION
3531 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3532 CCC     LOCAL VARIABLES:
3533 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3534 C ARGU   : (COSH(MU)*U-DMUFAC)/SINH(U)
3535 C COSCHI : COS(CHI)
3536 C COSHM  : COSH(MU)
3537 C D1     : COSH(MU)*U-DMUFAC
3538 C DMU    : SOLUTION MU OF COSH(MU)=PNU/X
3539 C DMUFAC : MU*COSH(MU)-SINH(MU)
3540 C SINCHI : SIN(CHI)
3541 C SINHM  : SINH(MU)
3542 C SINHU  : SINH(U)
3543 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3544        DOUBLE PRECISION ARGU,COSCHI,COSHM,D1,DMU,DMUFAC,
3545      + PI,SINCHI,SINHM,SINHU,U
3546        COMMON/PAROS1/DMU,COSHM,SINHM,DMUFAC,COSCHI,SINCHI
3547        PI=ACOS(-1.D0)
3548        D1=COSHM*U-DMUFAC
3549        SINHU=SINH(U)
3550        ARGU=D1/SINH(U)
3551        IF (ABS(ARGU).GT.1.D0) THEN
3552          ARGU=1.D0
3553        ENDIF
3554        SIGMA=ASIN(ARGU)
3555        IF (U.LT.DMU) SIGMA=PI-SIGMA
3556        RETURN
3557        END
3558
3559        DOUBLE PRECISION FUNCTION PHIB(U)
3560 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3561 CCC THIS FUNCTION COMPUTES THE FIRST EXPRESSION
3562 CCC OF EQ.(30) OF THE REFERENCE
3563 CCC  'COMPUTING SOLUTIONS OF THE MODIFIED BESSEL
3564 CCC   DIFFERENTIAL EQUATION FOR IMAGINARY ORDERS
3565 CCC   AND POSITIVE ARGUMENTS',
3566 CCC   A. GIL, J. SEGURA, N.M. TEMME
3567 CCC   ACM TRANS. MATH. SOFT. (2004)
3568 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3569 CCC        INPUT VARIABLE:
3570 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3571 C    U    : ARGUMENT OF THE FUNCTION
3572 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3573 CCC         LOCAL VARIABLES:
3574 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3575 C OVER      : OVERFLOW NUMBER
3576 C PNU       : ORDER OF THE FUNCTIONS
3577 C SIGMA(U)  : ASIN((COSH(MU)*U-DMUFAC)/SINH(U))
3578 C SIGMAU    : SIGMA(U)
3579 C X         : ARGUMENT OF THE FUNCTIONS
3580 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3581        DOUBLE PRECISION U,SIGMA,SIGMAU,OVER,D1MACH,
3582      + X,PNU
3583        COMMON/ARGU/X,PNU
3584 CCC MACHINE DEPENDENT CONSTANT (OVERFLOW NUMBER)
3585        OVER=D1MACH(2)*1.D-8
3586        IF (U.LT.200.D0) THEN
3587          SIGMAU=SIGMA(U)
3588          PHIB=X*COSH(U)*COS(SIGMAU)+PNU*SIGMAU
3589        ELSE
3590          PHIB=OVER
3591        ENDIF
3592        RETURN
3593        END
3594