MINICERN/mathlib/gen/c/fconc64.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:04  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if !defined(CERNLIB_DOUBLE)
  11       FUNCTION RFCONC(X,TAU,M)
  12 #include "gen/defc64.inc"
  13      +           CGAMMA,CLOGAM
  14 #endif
  15 #if defined(CERNLIB_DOUBLE)
  16       FUNCTION DFCONC(X,TAU,M)
  17 #include "gen/imp64.inc"
  18 #include "gen/defc64.inc"
  19      +           WGAMMA,WLOGAM
  20 #endif
  21 C
  22 C     Computes the (real) Conical Function of the first kind
  23 C                 P M (-1/2 + I*TAU) (X)
  24 C     for M = 0 and M = 1. P m nu (x) is the Legendre function of
  25 C     the first kind.
  26 C     Based on K.S. Koelbig, A program for computing the conical
  27 C     functions of the first kind P ... for m = 0 and m = 1,
  28 C     Computer Phys. Comm. 23 (1981) 51-61.
  29 C
  30 #include "gen/defc64.inc"
  31      +     CGM,CLG,CRG,I,A,B,C,TI,R,RR,U(0:3),V(0:3),W(19)
  32       LOGICAL LM0,LM1,LTA
  33       CHARACTER NAME*(*)
  34       CHARACTER*80 ERRTXT
  35 #if !defined(CERNLIB_DOUBLE)
  36       PARAMETER (NAME = 'RFCONC')
  37 #endif
  38 #if defined(CERNLIB_DOUBLE)
  39       PARAMETER (NAME = 'RFCONC/DFCONC')
  40 #endif
  41       DIMENSION T(7),H(9),S(5),P(11),D(-1:6)
  42
  43
  44       PARAMETER (PI  = 3.14159 26535 89793 24D+0)
  45       PARAMETER (RPI = 1.77245 38509 05516 03D+0)
  46       PARAMETER (I = (0,1))
  47       PARAMETER (Z1 = 1, HF = Z1/2, TH = 1+HF, C1 = Z1/10, C2 = Z1/5)
  48       PARAMETER (RPH = 2/PI, RPW = 2/RPI, TW = 20, NMAX = 200)
  49
  50       DATA EPS /1D-14/
  51
  52 #if !defined(CERNLIB_DOUBLE)
  53       ENTRY FCONC(X,TAU,M)
  54 #endif
  55
  56 #if defined(CERNLIB_DOUBLE)
  57       FEK(ARG)=DELIKC(ARG)
  58       FEE(ARG)=DELIEC(ARG)
  59       FJ0(ARG)=DBESJ0(ARG)
  60       FJ1(ARG)=DBESJ1(ARG)
  61       FI0(ARG)=DBESI0(ARG)
  62       FI1(ARG)=DBESI1(ARG)
  63       CGM(CRG)=WGAMMA(CRG)
  64       CLG(CRG)=WLOGAM(CRG)
  65 #endif
  66 #if !defined(CERNLIB_DOUBLE)
  67       FEK(ARG)=RELIKC(ARG)
  68       FEE(ARG)=RELIEC(ARG)
  69       FJ0(ARG)= BESJ0(ARG)
  70       FJ1(ARG)= BESJ1(ARG)
  71       FI0(ARG)= BESI0(ARG)
  72       FI1(ARG)= BESI1(ARG)
  73       CGM(CRG)=CGAMMA(CRG)
  74       CLG(CRG)=CLOGAM(CRG)
  75 #endif
  76       LM0=M .EQ. 0
  77       LM1=M .EQ. 1
  78       IF(X .LE. -1 .OR. TAU .LT. 0 .OR. .NOT.(LM0 .OR. LM1)) THEN
  79        FC=0
  80        WRITE(ERRTXT,101) X,TAU,M
  81        CALL MTLPRT(NAME,'C331.1',ERRTXT)
  82        GO TO 99
  83       END IF
  84       FC=1-M
  85       IF(X .EQ. 1) GO TO 99
  86
  87       IF(TAU .EQ. 0) THEN
  88        IF(X .LT. 1) THEN
  89         Y=SQRT(HF*(1-X))
  90         IF(LM0) THEN
  91          FC=RPH*FEK(Y)
  92         ELSE
  93          FC=RPH*(FEE(Y)-HF*(1+X)*FEK(Y))/SQRT(1-X**2)
  94         END IF
  95        ELSE
  96         Y=SQRT((X-1)/(X+1))
  97         IF(LM0) THEN
  98          FC=RPH*FEK(Y)/SQRT(HF*(X+1))
  99         ELSE
 100          FC=RPH*SQRT(HF/(X-1))*(FEE(Y)-FEK(Y))
 101         END IF
 102        END IF
 103        GO TO 99
 104       END IF
 105
 106       TI=I*TAU
 107       FM=M
 108       IF(-1 .LT. X .AND. X .LE. 0) GO TO 11
 109       IF(0  .LT. X .AND. X .LE. C1 .AND. TAU .LE. 17) GO TO 11
 110       IF(C1 .LT. X .AND. X .LE. C2 .AND. TAU .LE.  5) GO TO 11
 111       IF(C1 .LT. X .AND. X .LE. C2 .AND. TAU .LE. 17) GO TO 12
 112       IF(C2 .LT. X .AND. X .LE. TH .AND. TAU .LE. 20) GO TO 12
 113       IF(TH .LT. X .AND. TAU .LE. MAX(TW,X)) GO TO 13
 114
 115       IF(X .LT. 1) THEN
 116        Y=SQRT(1-X**2)
 117        T(1)=ACOS(X)
 118        H(1)=TAU*T(1)
 119        B0=FI0(H(1))
 120        B1=FI1(H(1))
 121        Z=-1
 122       ELSE
 123        Y=SQRT(X**2-1)
 124        T(1)=LOG(X+Y)
 125        H(1)=TAU*T(1)
 126        B0=FJ0(H(1))
 127        B1=FJ1(H(1))
 128        Z=1
 129       END IF
 130       H(1)=T(1)*X/Y
 131       P(1)=1/TAU
 132       S(1)=1/T(1)
 133       DO 5 J = 2,7
 134       T(J)=T(J-1)*T(1)
 135     5 H(J)=H(J-1)*H(1)
 136       DO 6 J = 2,11
 137     6 P(J)=P(J-1)*P(1)
 138       DO 7 J = 2,5
 139     7 S(J)=S(J-1)*S(1)
 140 C
 141       IF(LM0) THEN
 142        D(-1)=0
 143        D(0)=1
 144        D(1)=(H(1)-1)/(8*T(1))
 145        D(2)=(9*H(2)+6*H(1)-15-Z*8*T(2))/(128*T(2))
 146        D(3)=5*(15*H(3)+27*H(2)+21*H(1)-63-Z*T(2)*(16*H(1)+24))/
 147      1  (1024*T(3))
 148        D(4)=7*(525*H(4)+1500*H(3)+2430*H(2)+1980*H(1)-6435
 149      1 +192*T(4)-Z*T(2)*(720*H(2)+1600*H(1)+2160))/(32768*T(4))
 150        D(5)=21*(2835*H(5)+11025*H(4)+24750*H(3)+38610*H(2)
 151      1  +32175*H(1)-109395+T(4)*(1984*H(1)+4032)
 152      2  -Z*T(2)*(4800*H(3)+15120*H(2)+26400*H(1)+34320))/
 153      3  (262144*T(5))
 154        D(6)=11*(218295*H(6)+1071630*H(5)+3009825*H(4)+6142500*
 155      1  H(3)+9398025*H(2)+7936110*H(1)-27776385+T(4)*(254016*H(2)
 156      2  +749952*H(1)+1100736)-Z*T(2)*(441000*H(4)+1814400*H(3)
 157      3  +4127760*H(2)+6552000*H(1)+8353800+31232*T(4)))/
 158      4  (4194304*T(6))
 159       ELSE
 160        D(-1)=-1
 161        D(0)=3*(1-H(1))/(8*T(1))
 162        D(1)=(-15*H(2)+6*H(1)+9+Z*8*T(2))/(128*T(2))
 163        D(2)=3*(-35*H(3)-15*H(2)+15*H(1)+35+Z*T(2)*(32*H(1)+8))/
 164      1  (1024*T(3))
 165        D(3)=(-4725*H(4)-6300*H(3)-3150*H(2)+3780*H(1)+10395
 166      1   -1216*T(4)+Z*T(2)*(6000*H(2)+5760*H(1)+1680))/(32768*T(4))
 167        D(4)=7*(-10395*H(5)-23625*H(4)-28350*H(3)-14850*H(2)
 168      1  +19305*H(1)+57915-T(4)*(6336*H(1)+6080)+Z*T(2)*
 169      2  (16800*H(3)+30000*H(2)+25920*H(1)+7920))/(262144*T(5))
 170        D(5)=(-2837835*H(6)-9168390*H(5)-16372125*H(4)-18918900*H(3)
 171      1  -10135125*H(2)+13783770*H(1)+43648605-T(4)*(3044160*H(2)
 172      2  +5588352*H(1)+4213440)+Z*T(2)*(5556600*H(4)+14817600*H(3)
 173      3  +20790000*H(2)+17297280*H(1)+5405400+323072*T(4)))/
 174      4  (4194304*T(6))
 175        D(6)=0
 176       END IF
 177       S0=D(0)+(-4*D(3)*S(1)+D(4))*P(4)+
 178      1 (-192*D(5)*S(3)+144*D(6)*S(2))*P(8)
 179      2 +Z*(-D(2)*P(2)+(-24*D(4)*S(2)+12*D(5)*S(1)-D(6))*P(6)
 180      3 +(-1920*D(6)*S(4))*P(10))
 181       S1=D(1)*P(1)+(8*(D(3)*S(2)-D(4)*S(1))+D(5))*P(5)
 182      1 +(384*D(5)*S(4)-768*D(6)*S(3))*P(9)
 183      2 +Z*(D(-1)*TAU+(2*D(2)*S(1)-D(3))*P(3)+(48*D(4)*S(3)
 184      3 -72*D(5)*S(2)+18*D(6)*S(1))*P(7)+(3840*D(6)*S(5))*P(11))
 185       FC=SQRT(T(1)/Y)*(B0*S0+B1*S1)
 186       GO TO 99
 187
 188    11 LTA=TAU .LE. 10
 189       X1=X**2
 190       A=HF*((HF-FM)-TI)
 191       B=HF*((HF-FM)+TI)
 192       C=HF
 193       ASSIGN 1 TO JP
 194       GO TO 20
 195     1 R1=R
 196       R1=R1/ABS(CGM(A+HF))**2
 197       A=HF*((TH-FM)-TI)
 198       B=HF*((TH-FM)+TI)
 199       C=TH
 200       ASSIGN 2 TO JP
 201       GO TO 20
 202     2 R2=R
 203       FC=RPI*(R1-2*X*R2/ABS(CGM(A-HF))**2)
 204       IF(LM1) FC=2*FC/SQRT(1-X1)
 205       GO TO 99
 206
 207    12 LTA=X .GT. 1 .OR. X .LE. 1 .AND. TAU .LE. 5
 208       X1=HF*(1-X)
 209       A=(HF+FM)-TI
 210       B=(HF+FM)+TI
 211       C=FM+1
 212       ASSIGN 3 TO JP
 213       GO TO 20
 214     3 FC=R
 215       IF(LM1) FC=SIGN(HF,1-X)*(TAU**2+HF**2)*SQRT(ABS(X**2-1))*FC
 216       GO TO 99
 217
 218    13 LTA=.TRUE.
 219       X1=1/X**2
 220       A=HF*((HF-FM)-TI)
 221       B=HF*((TH-FM)-TI)
 222       C=1-TI
 223       ASSIGN 4 TO JP
 224       GO TO 20
 225     4 R1=EXP((TI-HF)*LOG(X+X)+CLG(1+TI)-CLG((TH-FM)+TI))*
 226      1        R*((HF-FM)+TI)/TI
 227       FC=RPW*R1
 228       IF(LM1) FC=FC/SQRT(1-X1)
 229       GO TO 99
 230
 231    20 IF(LTA) THEN
 232        Y=-X1
 233        Y2=Y**2
 234        Y3=Y*Y2
 235        W(1)=A+1
 236        W(2)=A+2
 237        W(3)=B+1
 238        W(4)=B+2
 239        W(5)=C+1
 240        W(6)=C*W(5)
 241        W(7)=A+B
 242        W(8)=A*B
 243        W(9)=(W(8)/C)*Y
 244        W(10)=W(1)*W(3)
 245        W(11)=W(2)*W(4)
 246        W(12)=1+(W(11)/(W(5)+W(5)))*Y
 247        W(13)=W(7)-6
 248        W(14)=W(7)+6
 249        W(15)=2-W(8)
 250        W(16)=W(15)-W(7)-W(7)
 251
 252        V(0)=1
 253        V(1)=1+(W(10)/(C+C))*Y
 254        V(2)=W(12)+(W(10)*W(11)/(12*W(6)))*Y2
 255        U(0)=1
 256        U(1)=V(1)-W(9)
 257        U(2)=V(2)-W(9)*W(12)+(W(8)*W(10)/(W(6)+W(6)))*Y2
 258
 259        R=1
 260        DO 21 N = 3,NMAX
 261        FN=N
 262        RR=R
 263        H(1)=FN-1
 264        H(2)=FN-2
 265        H(3)=FN-3
 266        H(4)=FN+FN
 267        H(5)=H(4)-3
 268        H(6)=H(5)+H(5)
 269        H(7)=4*(H(4)-1)*H(5)
 270        H(8)=8*H(5)**2*(H(4)-5)
 271        H(9)=3*FN**2
 272        W(1)=A+H(1)
 273        W(2)=A+H(2)
 274        W(3)=B+H(1)
 275        W(4)=B+H(2)
 276        W(5)=C+H(1)
 277        W(6)=C+H(2)
 278        W(7)=C+H(3)
 279        W(8)=H(2)-A
 280        W(9)=H(2)-B
 281        W(10)=H(1)-C
 282        W(11)=W(1)*W(3)
 283        W(12)=W(5)*W(6)
 284
 285        W(17)=1+((H(9)+W(13)*FN+W(16))/(H(6)*W(5)))*Y
 286        W(18)=-((W(11)*W(10)/H(6)+(H(9)-W(14)*FN+W(15))*W(11)*Y/H(7))/
 287      1          W(12))*Y
 288        W(19)=(W(2)*W(11)*W(4)*W(8)*W(9)/(H(8)*W(7)*W(12)))*Y3
 289        V(3)=W(17)*V(2)+W(18)*V(1)+W(19)*V(0)
 290        U(3)=W(17)*U(2)+W(18)*U(1)+W(19)*U(0)
 291        R=U(3)/V(3)
 292        IF(ABS(R-RR) .LT. EPS) GO TO JP, (1,2,3,4)
 293        DO 22 J = 1,3
 294        V(J-1)=V(J)
 295    22  U(J-1)=U(J)
 296    21  CONTINUE
 297       ELSE
 298        W(1)=X1*A*B/C
 299        R=1+W(1)
 300        DO 23 N = 1,NMAX
 301        FN=N
 302        RR=R
 303        W(1)=W(1)*X1*(A+FN)*(B+FN)/((C+FN)*(FN+1))
 304        R=R+W(1)
 305        IF(ABS(R-RR) .LT. EPS) GO TO JP, (1,2,3,4)
 306    23  CONTINUE
 307       END IF
 308       FC=0
 309       WRITE(ERRTXT,102) X
 310       CALL MTLPRT(NAME,'C331.2',ERRTXT)
 311 #if defined(CERNLIB_DOUBLE)
 312    99 DFCONC=FC
 313 #endif
 314 #if !defined(CERNLIB_DOUBLE)
 315    99 RFCONC=FC
 316 #endif
 317       RETURN
 318
 319   101 FORMAT('ILLEGAL ARGUMENT(S)  X = ',D15.8,' TAU = ',D15.8,
 320      1       ' M = ',I3)
 321   102 FORMAT('CONVERGENCE PROBLEM FOR HYPERGEOMETRIC FUNCTION, X = ',
 322      1        D15.8)
 323       END