MINICERN/mathlib/gen/c/reli364.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:11  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if defined(CERNLIB_DOUBLE)
  11       FUNCTION DELI3(X,AKP,P)
  12 C
  13 #include "gen/imp64.inc"
  14 C
  15       CHARACTER*(*) NAME
  16       PARAMETER(NAME='RELI3/DELI3')
  17 #endif
  18 #if !defined(CERNLIB_DOUBLE)
  19       FUNCTION RELI3(X,AKP,P)
  20 C
  21       CHARACTER*(*) NAME
  22       PARAMETER(NAME='RELI3')
  23 #endif
  24 C
  25 C     Translation of Algol procedure el3(x,kc,p) in
  26 C      R. BULIRSCH Numerical Calculation of Elliptic Integrals and
  27 C      Elliptic Functions III.,  Numer. Math. 13 (1969) 305-315
  28 C
  29       LOGICAL LBO,LBK
  30
  31       PARAMETER (ID = 16, IB = 4)
  32       PARAMETER (PI  = 3.14159 26535 89793 24D0)
  33       PARAMETER (AL2 = 0.69314 71805 59945 31D0)
  34       PARAMETER (ALB = IB*AL2, RLB = 1/ALB)
  35       PARAMETER (Z1 = 1, Z10 = 10, HF = Z1/2, C1 = Z1/10)
  36       PARAMETER (ND = ID-2, CA = Z10**(-ID/2), CB = Z10**(-(ID+2)))
  37       PARAMETER (ZD = HF/(ND+1))
  38
  39       CHARACTER*80 ERRTXT
  40
  41       DIMENSION RA(2:ND),RB(2:ND),RR(2:ND)
  42
  43       DATA (RB(K),RA(K),K=2,ND)
  44      2/2.50000000000000000D-01, 7.50000000000000000D-01,
  45      3 1.66666666666666667D-01, 8.33333333333333333D-01,
  46      4 1.25000000000000000D-01, 8.75000000000000000D-01,
  47      5 1.00000000000000000D-01, 9.00000000000000000D-01,
  48      6 8.33333333333333333D-02, 9.16666666666666667D-01,
  49      7 7.14285714285714286D-02, 9.28571428571428571D-01,
  50      8 6.25000000000000000D-02, 9.37500000000000000D-01,
  51      9 5.55555555555555556D-02, 9.44444444444444445D-01,
  52      A 5.00000000000000000D-02, 9.50000000000000000D-01,
  53      B 4.54545454545454545D-02, 9.54545454545454545D-01,
  54      C 4.16666666666666667D-02, 9.58333333333333333D-01,
  55      D 3.84615384615384615D-02, 9.61538461538461538D-01,
  56      E 3.57142857142857143D-02, 9.64285714285714286D-01/
  57
  58       IF(X .EQ. 0) THEN
  59        H=0
  60        GO TO 9
  61       ENDIF
  62       HH=X**2
  63       F=P*HH
  64       S=AKP
  65       IF(S .EQ. 0) S=CA/(1+ABS(X))
  66       T=S**2
  67       PM=HF*T
  68       E=HH*T
  69       Z=ABS(F)
  70       R=ABS(P)
  71       H1=HH+1
  72       IF(E .LT. C1 .AND. Z .LT. C1 .AND. T .LT. 1 .AND. R .LT. 1) THEN
  73        S=P+PM
  74        DO 1 K = 2,ND
  75        RR(K)=S
  76        PM=PM*T*RA(K)
  77     1  S=S*P+PM
  78        U=S*ZD
  79        S=U
  80        LBO=.FALSE.
  81        DO 2 K = ND,2,-1
  82        U=U+(RR(K)-U)*RB(K)
  83        LBO=.NOT.LBO
  84        V=U
  85        IF(LBO) V=-U
  86     2  S=S*HH+V
  87        IF(LBO) S=-S
  88        U=HF*(U+1)
  89 #if defined(CERNLIB_DOUBLE)
  90        H=(U-S*H1)*SQRT(H1)*X+U*DASINH(X)
  91 #endif
  92 #if !defined(CERNLIB_DOUBLE)
  93        H=(U-S*H1)*SQRT(H1)*X+U*ASINH(X)
  94 #endif
  95       ELSE
  96        W=1+F
  97        IF(W .EQ. 0) THEN
  98         H=0
  99         WRITE(ERRTXT,103) X,P
 100         CALL MTLPRT(NAME,'C346.3',ERRTXT)
 101         GO TO 9
 102        ENDIF
 103        PP=P
 104        IF(PP .EQ. 0) PP=CB/HH
 105        S=ABS(S)
 106        Y=ABS(X)
 107        G=PP-1
 108        IF(G .EQ. 0) G=CB
 109        F=PP-T
 110        IF(F .EQ. 0) F=CB*T
 111        AM=1-T
 112        AP=1+E
 113        R=PP*H1
 114        FA=G/(F*PP)
 115        LBO=FA .GT. 0
 116        FA=ABS(FA)
 117        PZ=ABS(G*F)
 118        DE=SQRT(PZ)
 119        Q=SQRT(ABS(PP))
 120        PM=PP-MIN(PM,HF)
 121        IF(PM .GE. 0) THEN
 122         U=SQRT(R*AP)
 123         V=Y*DE
 124         IF(G .LT. 0) V=-V
 125         D=1/Q
 126         C=1
 127        ELSE
 128         U=SQRT(H1*AP*PZ)
 129         YE=Y*Q
 130         V=AM*YE
 131         Q=-DE/G
 132         D=-AM/DE
 133         C=0
 134         PZ=AP-R
 135        ENDIF
 136        IF(LBO) THEN
 137         R=V/U
 138         Z=1
 139         K=1
 140         IF(PM .LT. 0) THEN
 141          H1=Y*SQRT(H1/(AP*FA))
 142          H1=1/H1-H1
 143          Z=H1-2*R
 144          R=2+R*H1
 145          IF(R .EQ. 0) R=CB
 146          IF(Z .EQ. 0) Z=H1*CB
 147          R=R/Z
 148          Z=R
 149          W=PZ
 150         ENDIF
 151         U=U/W
 152         V=V/W
 153        ELSE
 154         T=U+ABS(V)
 155         LBK=.TRUE.
 156         IF(PP .LT. 0) THEN
 157          DE=V/PZ
 158          YE=2*U*YE
 159          U=T/PZ
 160          V=-(F+G*E)/T
 161          T=PZ*ABS(W)
 162          Z=(HH*R*F-G*AP+YE)/T
 163          YE=YE/T
 164         ELSE
 165          DE=V/W
 166          YE=0
 167          U=(E+PP)/T
 168          V=T/W
 169          Z=1
 170         ENDIF
 171         IF(S .GT. 1) THEN
 172          H1=U
 173          U=V
 174          V=H1
 175         ENDIF
 176        ENDIF
 177        Y=1/Y
 178        E=S
 179        T=1
 180        N=1
 181        L=0
 182        M=0
 183     3  Y=Y-E/Y
 184        IF(Y .EQ. 0) Y =SQRT(E)*CB
 185        F=C
 186        C=D/Q+C
 187        G=E/Q
 188        D=2*(F*G+D)
 189        Q=G+Q
 190        G=T
 191        T=S+T
 192        N=2*N
 193        M=2*M
 194        IF(LBO) THEN
 195         IF(Z .LT. 0) M=K+M
 196         K=0
 197         IF(R .GT. 0) K=1
 198         IF(R .LT. 0) K=-1
 199         H1=E/(U**2+V**2)
 200         U=U*(1+H1)
 201         V=V*(1-H1)
 202        ELSE
 203         R=U/V
 204         H1=Z*R
 205         Z=H1*Z
 206         HH=E/V
 207         IF(LBK) THEN
 208          DE=DE/U
 209          YE=YE*(H1+1/H1)+DE*(1+R)
 210          DE=DE*(U-HH)
 211          LBK=ABS(YE) .LT. 1
 212         ELSE
 213          A=LOG(Z)
 214          K=RLB*A
 215          IF(Z .GE. 1) K=K+1
 216          Z=EXP(A-K*ALB)
 217          M=M+K
 218         ENDIF
 219        ENDIF
 220        IF(ABS(G-S) .GT. CA*G) THEN
 221         IF(LBO) THEN
 222          G=HF*(1/R-R)
 223          HH=U+V*G
 224          H1=G*U-V
 225          IF(HH .EQ. 0) HH=U*CB
 226          IF(H1 .EQ. 0) H1=V*CB
 227          Z=R*H1
 228          R=HH/H1
 229         ELSE
 230          U=U+E/U
 231          V=V+HH
 232         ENDIF
 233         S=2*SQRT(E)
 234         E=S*T
 235         L=2*L
 236         IF(Y .LT. 0) L=L+1
 237         GO TO 3
 238        ENDIF
 239        IF(Y .LT. 0) L=L+1
 240        E=(ATAN(T/Y)+PI*L)*(C*T+D)/(T*(T+Q))
 241        IF(LBO) THEN
 242         H1=V/(T+U)
 243         Z=1-R*H1
 244         H1=R+H1
 245         IF(Z .EQ. 0) Z=CB
 246         IF(Z .LT. 0) THEN
 247          IF(H1 .GT. 0) M=M+1
 248          IF(H1 .LT. 0) M=M-1
 249         ENDIF
 250         S=ATAN(H1/Z)+M*PI
 251        ELSE
 252         IF(LBK) THEN
 253 #if defined(CERNLIB_DOUBLE)
 254          S=HF*DASINH(YE)
 255 #endif
 256 #if !defined(CERNLIB_DOUBLE)
 257          S=HF*ASINH(YE)
 258 #endif
 259         ELSE
 260          S=HF*(LOG(Z)+M*ALB)
 261         ENDIF
 262        ENDIF
 263        H=(E+SQRT(FA)*S)/N
 264        IF(X .LT. 0) H=-H
 265       ENDIF
 266 #if defined(CERNLIB_DOUBLE)
 267     9 DELI3=H
 268 #endif
 269 #if !defined(CERNLIB_DOUBLE)
 270     9 RELI3=H
 271 #endif
 272       RETURN
 273   103 FORMAT('FUNCTION SINGULAR FOR ',1P,
 274      1 'X = ',D15.8,3X,'P = ',D15.8,4X,'(P*X**2 = -1)')
 275       END