MINICERN/mathlib/gen/b/rsrtnt64.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:01:51  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if defined(CERNLIB_DOUBLE)
  11       SUBROUTINE DSRTNT(K,N,A,B,C,U1,V1,RES,LRL)
  12 #endif
  13 #if !defined(CERNLIB_DOUBLE)
  14       SUBROUTINE RSRTNT(K,N,A,B,C,U1,V1,RES,LRL)
  15 #endif
  16 #include "gen/imp64.inc"
  17       LOGICAL LRL,LLL
  18       CHARACTER NAME*(*)
  19       CHARACTER*80 ERRTXT
  20 #if defined(CERNLIB_DOUBLE)
  21       PARAMETER (NAME = 'DSRTNT')
  22 #endif
  23 #if !defined(CERNLIB_DOUBLE)
  24       PARAMETER (NAME = 'RSRTNT')
  25 #endif
  26       PARAMETER (R1 = 1, HF = R1/2)
  27
  28       DIMENSION BK(0:4,0:4),SGN(0:4)
  29
  30       DATA BK(0,0) /1/, (BK(1,J),J=0,1) /1,1/, (BK(2,J),J=0,2) /1,2,1/
  31       DATA (BK(3,J),J=0,3) /1,3,3,1/, (BK(4,J),J=0,4) /1,4,6,4,1/
  32       DATA (SGN(J),J=0,4) /1,-1,1,-1,1/
  33
  34       P(X)=A+B*X+C*X**2
  35       RT(X)=SQRT(P(X))
  36
  37       H=0
  38       LLL=ABS(A)+ABS(B)+ABS(C) .GT. 0
  39       IF(.NOT.LLL) THEN
  40        CALL MTLPRT(NAME,'B300.1','A = B = C = 0')
  41        GO TO 9
  42       ENDIF
  43       LLL=ABS(K) .LE. 3 .AND. (N .EQ. 1 .OR. N .EQ. 3)
  44       IF(.NOT.LLL) THEN
  45        WRITE(ERRTXT,102) K,N
  46        CALL MTLPRT(NAME,'B300.2',ERRTXT)
  47        GO TO 9
  48       ENDIF
  49
  50       U=MIN(U1,V1)
  51       V=MAX(U1,V1)
  52       LLL=U .EQ. V
  53       IF(LLL) GO TO 9
  54
  55       LLL=K .GE. 0 .OR. K .EQ. -1 .AND. U*V .NE. 0
  56      1             .OR. K .LE. -2 .AND. U*V .GT. 0
  57       IF(.NOT.LLL) GO TO 9
  58
  59       DELTA=4*A*C-B**2
  60       RTD=SQRT(ABS(DELTA))
  61       IF(C .EQ. 0) THEN
  62        IF(B .LT. 0) THEN
  63         X0=-A/B
  64         LLL=U .LE. X0 .AND. V .LE. X0
  65        ELSEIF(B .EQ. 0) THEN
  66         LLL=A .GT. 0
  67        ELSE
  68         X0=-A/B
  69         LLL=U .GE. X0 .AND. V .GE. X0
  70        ENDIF
  71       ELSE
  72        IF(DELTA .GT. 0) THEN
  73         LLL=C .GT. 0
  74        ELSEIF(DELTA .EQ. 0) THEN
  75         IF(C .LT. 0) THEN
  76          LLL=.FALSE.
  77         ELSE
  78          X0=-B/(2*C)
  79          LLL=U .LT. X0 .AND. V .LT. X0 .OR. U .GT. X0 .AND. V .GT. X0
  80         ENDIF
  81        ELSE
  82         A1=(-B+RTD)/(2*C)
  83         A2=(-B-RTD)/(2*C)
  84         W1=MIN(A1,A2)
  85         W2=MAX(A1,A2)
  86         IF(C .LT. 0) THEN
  87          LLL=W1 .LE. U .AND. U .LE. W2 .AND. W1 .LE. V .AND. V .LE. W2
  88         ELSE
  89          LLL=U .LE. W1 .AND. V .LE. W1 .OR. U .GE. W2 .AND. V .GE. W2
  90         ENDIF
  91        ENDIF
  92       ENDIF
  93       IF(.NOT.LLL) GO TO 9
  94
  95       IF(K .GE. 0) THEN
  96        IF(C .EQ. 0) THEN
  97         IF(B .EQ. 0) THEN
  98          K1=K+1
  99          H=(V**K1-U**K1)/(K1*SQRT(A)**N)
 100         ELSE
 101          N1=2-N
 102          XV=A+B*V
 103          XU=A+B*U
 104          IF(A .EQ. 0) THEN
 105           K1=2*K+N1
 106           H=2*(SQRT(XV)**K1-SQRT(XU)**K1)/(K1*B**(K+1))
 107          ELSE
 108           HV=SQRT(XV)**N1
 109           HU=SQRT(XU)**N1
 110           H1=-XV/A
 111           H2=-XU/A
 112           S=N1*(HV-HU)
 113           DO 1 J = 1,K
 114     1     S=S+BK(K,J)*(H1**J*HV-H2**J*HU)/(2*J+N1)
 115           H=2*(-A/B)**K*S/B
 116          ENDIF
 117         ENDIF
 118        ELSE
 119         IF(N .EQ. 1) THEN
 120          ASSIGN 11 TO JMP1
 121          GO TO 10
 122    11    IF(K .EQ. 0) THEN
 123           H=H
 124          ELSEIF(K .EQ. 1) THEN
 125           H=(RT(V)-RT(U)-HF*B*H)/C
 126          ELSEIF(K .EQ. 2) THEN
 127           H1=4*C
 128           H2=6*B
 129           H=((H1*V-H2)*RT(V)-(H1*U-H2)*RT(U)+(2*B**2-DELTA)*H)/(8*C**2)
 130          ELSEIF(K .EQ. 3) THEN
 131           H1=8*C**2
 132           H2=10*B*C
 133           G1=15*B**2
 134           G2=A*C
 135           H3=G1-16*G2
 136           H=(((H1*V-H2)*V+H3)*RT(V)-((H1*U-H2)*U+H3)*RT(U)-
 137      1       (HF*G1-18*G2)*B*H)/(24*C**3)
 138          ENDIF
 139         ELSE
 140          IF(DELTA .EQ. 0) THEN
 141           IF(B .EQ. 0) THEN
 142            IF(K .EQ. 2) THEN
 143             H=LOG(ABS(V/U))/SQRT(C)**3
 144            ELSE
 145             K1=K-2
 146             H=(V**K1-U**K1)/(K1*SQRT(C)**3)
 147            ENDIF
 148            H=SIGN(R1,U)*H
 149           ELSE
 150            X0=B/(2*C)
 151            XV=V+X0
 152            XU=U+X0
 153            HV=1/XV**2
 154            HU=1/XU**2
 155            H1=-XV/X0
 156            H2=-XU/X0
 157            S=HF*(HU-HV)
 158            DO 2 J = 1,K
 159            IF(J .NE. 2) THEN
 160             S=S+BK(K,J)*(H1**J*HV-H2**J*HU)/(J-2)
 161            ELSE
 162             S=S+BK(K,2)*LOG(ABS(XV/XU))/X0**2
 163            ENDIF
 164     2      CONTINUE
 165            H=(-X0)**K*S/SQRT(C)**3
 166            IF(U .LT. -X0) H=-H
 167           ENDIF
 168          ELSE
 169           IF(K .EQ. 0) THEN
 170            H1=2*C
 171            H=2*((H1*V+B)/RT(V)-(H1*U+B)/RT(U))/DELTA
 172           ELSEIF(K .EQ. 1) THEN
 173            H1=2*A
 174            H=2*((H1+B*U)/RT(U)-(H1+B*V)/RT(V))/DELTA
 175            LB1=11
 176           ELSEIF(K .EQ. 2) THEN
 177            ASSIGN 12 TO JMP1
 178            GO TO 10
 179    12      H1=DELTA-B**2
 180            H2=2*A*B
 181            H=(((H1*U-H2)/RT(U)-(H1*V-H2)/RT(V))/DELTA+H)/C
 182           ELSEIF(K .EQ. 3) THEN
 183            ASSIGN 13 TO JMP1
 184            GO TO 10
 185    13      H1=C*DELTA
 186            G1=A*C
 187            G2=3*B**2
 188            H2=B*(10*G1-G2)
 189            H3=A*(8*G1-G2)
 190            H=(2*(((H1*V+H2)*V+H3)/RT(V)-((H1*U+H2)*U+H3)/RT(U))/
 191      1       DELTA-3*B*H)/(2*C**2)
 192           ENDIF
 193          ENDIF
 194         ENDIF
 195        ENDIF
 196       ELSE
 197        IF(A .EQ. 0) THEN
 198         IF(B .EQ. 0) THEN
 199          K1=K-N+1
 200          H=SIGN(R1,U)*(V**K1-U**K1)/(K1*SQRT(C)**N)
 201         ELSE
 202          IF(C .EQ. 0) THEN
 203           K1=2*K-N+2
 204           H=2*(SQRT(B*V)**K1-SQRT(B*U)**K1)/(K1*B**(K+1))
 205          ELSE
 206           XV=SQRT(C+B/V)
 207           XU=SQRT(C+B/U)
 208           N1=N-2
 209           K1=-K+N1
 210           S=0
 211           DO 4 J = 0,K1
 212           KJ=2*J-N1
 213     4     S=S+SGN(J)*BK(K1,J)*(XU**KJ-XV**KJ)/(KJ*C**J)
 214           H=2*(-C/B)**K1*S/B
 215           IF(U .LT. 0 .AND. V .LT. 0) H=-H
 216          ENDIF
 217         ENDIF
 218        ELSE
 219         IF(N .EQ. 1) THEN
 220          ASSIGN 21 TO JMP2
 221          GO TO 20
 222    21    IF(K .EQ. -1) THEN
 223           H=H
 224          ELSEIF(K .EQ. -2) THEN
 225           H=(RT(U)/U-RT(V)/V-HF*B*H)/A
 226          ELSEIF(K .EQ. -3) THEN
 227           H1=6*B
 228           H2=4*A
 229           H=((H1*V-H2)*RT(V)/V**2-(H1*U-H2)*RT(U)/U**2+
 230      1       (3*B**2-H2*C)*H)/(8*A**2)
 231          ENDIF
 232         ELSE
 233          IF(DELTA .EQ. 0) THEN
 234           IF(C .EQ. 0) THEN
 235            IF(K .EQ. -1) THEN
 236             H=LOG(ABS(V/U))/SQRT(A)**3
 237            ELSE
 238             K1=K+1
 239             H=(V**K1-U**K1)/(K1*SQRT(A)**3)
 240            ENDIF
 241           ELSE
 242            X0=B/(2*C)
 243            XV=1+X0/V
 244            XU=1+X0/U
 245            K1=-K+1
 246            K2=-K-1
 247            S=0
 248            DO 3 J = 0,K1
 249            KJ=K2-J
 250            IF(KJ .NE. 0) THEN
 251             S=S+SGN(J)*BK(K1,J)*(XV**KJ-XU**KJ)/KJ
 252            ELSE
 253             S=S+SGN(K2)*BK(K1,K2)*LOG(ABS(XV/XU))
 254            ENDIF
 255     3      CONTINUE
 256            H=-S/(SQRT(C)**3*X0**(K1+1))
 257            IF(U. LT. -X0) H=-H
 258           ENDIF
 259          ELSE
 260           ASSIGN 22 TO JMP2
 261           GO TO 20
 262    22     IF(K .EQ. -1) THEN
 263            H1=B*C
 264            H2=B**2-2*A*C
 265            H=(2*((H1*U+H2)/RT(U)-(H1*V+H2)/RT(V))/DELTA+H)/A
 266           ELSEIF(K .EQ. -2) THEN
 267            G1=3*B**2
 268            G2=A*C
 269            H1=(G1-8*G2)*C
 270            H2=(G1-10*G2)*B
 271            H3=A*DELTA
 272            H=(((H1*V+H2-H3/V)/RT(V)-(H1*U+H2-H3/U)/RT(U))/DELTA
 273      1       -3*HF*B*H)/A**2
 274           ELSEIF(K .EQ. -3) THEN
 275            G1=A*DELTA
 276            G2=A*C
 277            G3=B**2
 278            G4=15*G3
 279            H1=2*A*G1
 280            H2=5*B*G1
 281            H3=(G4-62*G2)*G3+24*G2**2
 282            H4=B*C*(G4-52*G2)
 283            H=((((H2-H1/V)/V-H3-H4*V)/RT(V)-((H2-H1/U)/U-H3-H4*U)/RT(U))
 284      1        /DELTA+HF*(G4-12*G2)*H)/(4*A**3)
 285           ENDIF
 286          ENDIF
 287         ENDIF
 288        ENDIF
 289       ENDIF
 290       GO TO 9
 291
 292    10 C2=2*C
 293       IF(DELTA .GT. 0 .OR. DELTA .LT. 0 .AND. C .GT. 0) THEN
 294        H=LOG(ABS((2*SQRT(C*P(V))+C2*V+B)/
 295      1           (2*SQRT(C*P(U))+C2*U+B)))/SQRT(C)
 296       ELSEIF(DELTA .EQ. 0) THEN
 297        H=ABS(LOG(ABS((C2*V+B)/(C2*U+B))))/SQRT(C)
 298       ELSE
 299        H=(ASIN((C2*U+B)/RTD)-ASIN((C2*V+B)/RTD))/SQRT(-C)
 300       ENDIF
 301       GO TO JMP1, (11,12,13)
 302
 303    20 IF(C .EQ. 0) THEN
 304        IF(B .EQ. 0) THEN
 305         H=LOG(ABS(V/U))/SQRT(A)
 306        ELSE
 307         IF(A .LT. 0) THEN
 308          H=2*(ATAN(SQRT(-(A+B*V)/A))-ATAN(SQRT(-(A+B*U)/A)))/SQRT(-A)
 309         ELSE
 310          WA=SQRT(A)
 311          WU=SQRT(A+B*U)
 312          WV=SQRT(A+B*V)
 313          H=LOG(ABS((WV-WA)*(WU+WA)/((WV+WA)*(WU-WA))))/WA
 314         ENDIF
 315        ENDIF
 316       ELSE
 317        A2=2*A
 318        IF(DELTA .GT. 0 .OR. DELTA .LT. 0 .AND. A .GT. 0) THEN
 319         H=LOG(ABS((-2*SQRT(A*P(V))+B*V+A2)*U/
 320      1           ((-2*SQRT(A*P(U))+B*U+A2)*V)))/SQRT(A)
 321        ELSEIF(DELTA .EQ. 0) THEN
 322         H=LOG(ABS((B*U+A2)*V/((B*V+A2)*U)))/SQRT(A)
 323         IF(U*V .GT. 0) H=SIGN(H,U)
 324        ELSE
 325         H=(ASIN((B*V+A2)/(V*RTD))-ASIN((B*U+A2)/(U*RTD)))/SQRT(-A)
 326         IF(U .LT. 0 .AND. V .LT. 0) H=-H
 327        ENDIF
 328       ENDIF
 329       GO TO JMP2, (21,22)
 330
 331     9 RES=SIGN(R1,V1-U1)*H
 332       LRL=LLL
 333       RETURN
 334   102 FORMAT('ILLEGAL VALUE(S) K = ',I5,', N = ',I5)
 335       END