| fe4da5cc |
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 |
git://git.uio.no / u / mrichter / AliRoot.git / blame