| | 1 | c 1 2 3 4 5 6 7 8 |
| | 2 | c---------|---------|---------|---------|---------|---------|---------|---------| |
| | 3 | *-- Author : R.Lednicky 20/01/95 |
| | 4 | SUBROUTINE fsiw |
| | 5 | C======================================================================= |
| | 6 | C Calculates final state interaction (FSI) weights |
| | 7 | C WEIF = weight due to particle - (effective) nucleus FSI (p-N) |
| | 8 | C WEI = weight due to p-p-N FSI |
| | 9 | C WEIN = weight due to p-p FSI; note that WEIN=WEI if I3C=0; |
| | 10 | C note that if I3C=1 the calculation of |
| | 11 | C WEIN can be skipped by putting J=0 |
| | 12 | C....................................................................... |
| | 13 | C Correlation Functions: |
| | 14 | C CF(p-p-N) = sum(WEI)/sum(WEIF) |
| | 15 | C CF(p-p) = sum(WEIN)/sum(1); here the nucleus is completely |
| | 16 | C inactive |
| | 17 | C CF(p-p-"N") = sum(WEIN*WEIF')/sum(WEIF'), where WEIN and WEIF' |
| | 18 | C are not correlated (calculated at different emission |
| | 19 | C points, e.g., for different events); |
| | 20 | C thus here the nucleus affects one-particle |
| | 21 | C spectra but not the correlation |
| | 22 | C....................................................................... |
| | 23 | C User must supply data file <fn> on unit NUNIT (e.g. =11) specifying |
| | 24 | C LL : particle pair |
| | 25 | C NS : approximation used to calculate Bethe-Salpeter amplitude |
| | 26 | C ITEST: test switch |
| | 27 | C If ITEST=1 then also following parameters are required |
| | 28 | C ICH : 1(0) Coulomb interaction between the two particles ON (OFF) |
| | 29 | C IQS : 1(0) quantum statistics for the two particles ON (OFF) |
| | 30 | C ISI : 1(0) strong interaction between the two particles ON (OFF) |
| | 31 | C I3C : 1(0) Coulomb interaction with residual nucleus ON (OFF) |
| | 32 | C This data file can contain other information useful for the user. |
| | 33 | C It is read by subroutines READINT4 and READREA8(4) (or READ_FILE). |
| | 34 | C---------------------------------------------------------------------- |
| | 35 | C- LL 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
| | 36 | C- part. 1: n p n alfa pi+ pi0 pi+ n p pi+ pi+ pi+ pi- K+ K+ K+ K- |
| | 37 | C- part. 2: n p p alfa pi- pi0 pi+ d d K- K+ p p K- K+ p p |
| | 38 | C NS=1 y/n: + + + + + - - - - - - - - - - - - |
| | 39 | C---------------------------------------------------------------------- |
| | 40 | C- LL 18 19 20 21 22 23 24 25 26 27 28 |
| | 41 | C- part. 1: d d t t K0 K0 d p p p n |
| | 42 | C- part. 2: d alfa t alfa K0 K0b t t alfa lambda lambda |
| | 43 | C NS=1 y/n: - - - - - - - - - + + |
| | 44 | C---------------------------------------------------------------------- |
| | 45 | C NS=1 Square well potential, |
| | 46 | C NS=3 not used |
| | 47 | C NS=4 scattered wave approximated by the spherical wave, |
| | 48 | C NS=2 same as NS=4 but the approx. of equal emission times in PRF |
| | 49 | C not required (t=0 approx. used in all other cases). |
| | 50 | C Note: if NS=2,4, the B-S amplitude diverges at zero distance r* in |
| | 51 | C the two-particle c.m.s.; user can specify a cutoff AA in |
| | 52 | C SUBROUTINE FSIINI, for example: |
| | 53 | C IF(NS.EQ.2.OR.NS.EQ.4)AA=5.D0 !! in 1/GeV --> AA=1. fm |
| | 54 | C--------------------------------------------------------------------- |
| | 55 | C ITEST=1 any values of parameters ICH, IQS, ISI, I3C are allowed |
| | 56 | C and should be given in data file <fn> |
| | 57 | C ITEST=0 physical values of these parameters are put automatically |
| | 58 | C in FSIINI (their values are not required in data file) |
| | 59 | C===================================================================== |
| | 60 | C At the beginning of calculation user should call FSIINI, |
| | 61 | C which reads LL, NS, ITEST (and eventually ICH, IQS, ISI, I3C) |
| | 62 | C and initializes various parameters. |
| | 63 | C In particular the constants in |
| | 64 | C COMMON/FSI_CONS/PI,PI2,SPI,DR,W |
| | 65 | C may be useful for the user: |
| | 66 | C W=1/.1973D0 ! from fm to 1/GeV |
| | 67 | C PI=4*DATAN(1.D0) |
| | 68 | C PI2=2*PI |
| | 69 | C SPI=DSQRT(PI) |
| | 70 | C DR=180.D0/PI ! from radian to degree |
| | 71 | C _______________________________________________________ |
| | 72 | C !! |Important note: all real quantities are assumed REAL*8 | !! |
| | 73 | C ------------------------------------------------------- |
| | 74 | C For each event user should fill in the following information |
| | 75 | C in COMMONs (all COMMONs in FSI calculation start with FSI_): |
| | 76 | C ................................................................... |
| | 77 | C COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2 |
| | 78 | C Only |
| | 79 | C AMN = mass of the effective nucleus [GeV/c**2] |
| | 80 | C CN = charge of the effective nucleus [elem. charge units] |
| | 81 | C are required |
| | 82 | C ................................................................... |
| | 83 | C COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, !part. momenta in the rest frame |
| | 84 | C 1 P2X,P2Y,P2Z,E2,P2 !of effective nucleus (NRF) |
| | 85 | C Only the components |
| | 86 | C PiX,PiY,PiZ [GeV/c] |
| | 87 | C in NRF are required. |
| | 88 | C To make the corresponding Lorentz transformation user can use the |
| | 89 | C subroutines LTRAN and LTRANB |
| | 90 | C ................................................................... |
| | 91 | C COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emission |
| | 92 | C 1 X2,Y2,Z2,T2,R2 ! points in NRF |
| | 93 | C The componets |
| | 94 | C Xi,Yi,Zi [fm] |
| | 95 | C and emission times |
| | 96 | C Ti [fm/c] |
| | 97 | C should be given in NRF with the origin assumed at the center |
| | 98 | C of the effective nucleus. If the effect of residual nucleus is |
| | 99 | C not calculated within FSIW, the NRF can be any fixed frame. |
| | 100 | C----------------------------------------------------------------------- |
| | 101 | C Before calling FSIW the user must call |
| | 102 | C CALL LTRAN12 |
| | 103 | C Besides Lorentz transformation to pair rest frame: |
| | 104 | C (p1-p2)/2 --> k* it also transforms 4-coordinates of |
| | 105 | C emission points from fm to 1/GeV and calculates Ei,Pi and Ri. |
| | 106 | C Note that |k*|=AK in COMMON/FSI_PRF/ |
| | 107 | C----------------------------------------------------------------------- |
| | 108 | C After making some additional filtering using k* (say k* < k*max) |
| | 109 | C or direction of vector k*, |
| | 110 | C user can finally call FSIW to calculate the FSI weights |
| | 111 | C to be used to construct the correlation function |
| | 112 | C======================================================================= |
| | 113 | |
| | 114 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 115 | COMMON/FSI_JR/JRAT |
| | 116 | COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2 |
| | 117 | COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, ! particle momenta in NRF |
| | 118 | 1 P2X,P2Y,P2Z,E2,P2 |
| | 119 | COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS, ! k*=(p1-p2)/2 and x1-x2 |
| | 120 | 1 X,Y,Z,T,RP,RPS ! in pair rest frame (PRF) |
| | 121 | COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, !4-coord. of emis. points in NRF |
| | 122 | 1 X2,Y2,Z2,T2,R2 |
| | 123 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| | 124 | COMMON/FSI_FFPN/FF12,FF21 |
| | 125 | COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS |
| | 126 | COMPLEX*16 FF12,FF21 |
| | 127 | C------------------------------------------------------------------ |
| | 128 | C==> AC1,2 = "relativistic" Bohr radii for particle-nucleus systems |
| | 129 | C1N=C1*CN |
| | 130 | IF(C1N.NE.0.D0)AC1=137.036D0/(C1N*E1) !m1-->E1 |
| | 131 | C2N=C2*CN |
| | 132 | IF(C2N.NE.0.D0)AC2=137.036D0/(C2N*E2) !m2-->E2 |
| | 133 | |
| | 134 | C----------------------------------------------------------- |
| | 135 | J=1 |
| | 136 | CALL FSIPN !weight due to particle-nucleus FSI |
| | 137 | JRAT=0 |
| | 138 | CALL FSIWF !weight due to particle-particle-nucleus FSI |
| | 139 | WEIN=WEI |
| | 140 | IF(I3C*J.NE.0) THEN |
| | 141 | FF12=DCMPLX(1.D0,0.D0) |
| | 142 | FF21=DCMPLX(1.D0,0.D0) |
| | 143 | JRAT=1 |
| | 144 | CALL VZ ! weight due to particle-particle FSI |
| | 145 | ENDIF |
| | 146 | RETURN |
| | 147 | END |
| | 148 | |
| | 149 | C======================================================================= |
| | 150 | |
| | 151 | SUBROUTINE LTRAN(P0,P,PS) |
| | 152 | C==>calculating particle 4-momentum PS={PSX,PSY,PSZ,ES} |
| | 153 | C in rest frame of a system 0 with 4-momentum P0={P0X,P0Y,P0Z,E0} |
| | 154 | C from its 4-momentum P={PX,PY,PZ,E} |
| | 155 | |
| | 156 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 157 | DIMENSION P0(4),P(4),PS(4) |
| | 158 | C----------------------------------------------------------------------- |
| | 159 | P0S=P0(1)**2+P0(2)**2+P0(3)**2 |
| | 160 | AM0=DSQRT(P0(4)**2-P0S) |
| | 161 | EPM=P0(4)+AM0 |
| | 162 | PP0=P(1)*P0(1)+P(2)*P0(2)+P(3)*P0(3) |
| | 163 | H=(PP0/EPM-P(4))/AM0 |
| | 164 | PS(1)=P(1)+P0(1)*H |
| | 165 | PS(2)=P(2)+P0(2)*H |
| | 166 | PS(3)=P(3)+P0(3)*H |
| | 167 | PS(4)=(P0(4)*P(4)-PP0)/AM0 |
| | 168 | RETURN |
| | 169 | END |
| | 170 | |
| | 171 | SUBROUTINE LTRANB(P0,PS,P) |
| | 172 | C==>calculating particle 4-momentum P={PX,PY,PZ,E} |
| | 173 | C from its 4-momentum PS={PSX,PSY,PSZ,ES} |
| | 174 | C in rest frame of a system 0 with 4-momentum P0={P0X,P0Y,P0Z,E0} |
| | 175 | |
| | 176 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 177 | DIMENSION P0(4),P(4),PS(4) |
| | 178 | C----------------------------------------------------------------------- |
| | 179 | P0S=P0(1)**2+P0(2)**2+P0(3)**2 |
| | 180 | AM0=DSQRT(P0(4)**2-P0S) |
| | 181 | EPM=P0(4)+AM0 |
| | 182 | PSP0=PS(1)*P0(1)+PS(2)*P0(2)+PS(3)*P0(3) |
| | 183 | HS=(PSP0/EPM+PS(4))/AM0 |
| | 184 | P(1)=PS(1)+P0(1)*HS |
| | 185 | P(2)=PS(2)+P0(2)*HS |
| | 186 | P(3)=PS(3)+P0(3)*HS |
| | 187 | P(4)=(P0(4)*PS(4)+PSP0)/AM0 |
| | 188 | RETURN |
| | 189 | END |
| | 190 | |
| | 191 | |
| | 192 | FUNCTION GDD(X,J) |
| | 193 | C--- GDD = k*COTG(DELTA), X=k^2 |
| | 194 | C--- J=1,2,3 corresp. to S=0,1,2 |
| | 195 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 196 | COMMON/FSI_AADD/AADD(20,3) |
| | 197 | COMMON/FSI_C/C(10),AM,AMS,DM |
| | 198 | COMMON/FSI_CONS/PI,PI2,SPI,DR,W |
| | 199 | COMPLEX*16 C |
| | 200 | E=X/2/AM |
| | 201 | ER=DSQRT(E) |
| | 202 | IF(J.EQ.1)THEN |
| | 203 | GDD=ER*(AADD(1,1)*DEXP(-E/AADD(2,1))-AADD(3,1)) |
| | 204 | GDD=GDD/DR ! from degree to radian |
| | 205 | TAND=DTAN(GDD) |
| | 206 | IF(TAND.EQ.0.D0)TAND=1.D-10 |
| | 207 | GDD=DSQRT(X)/TAND |
| | 208 | END IF |
| | 209 | IF(J.EQ.2)THEN |
| | 210 | GDD=1.D10 |
| | 211 | END IF |
| | 212 | IF(J.EQ.3)THEN |
| | 213 | GDD=ER*(AADD(1,3)+AADD(2,3)*E) |
| | 214 | GDD=GDD/DR ! from degree to radian |
| | 215 | TAND=DTAN(GDD) |
| | 216 | IF(TAND.EQ.0.D0)TAND=1.D-10 |
| | 217 | GDD=DSQRT(X)/TAND |
| | 218 | END IF |
| | 219 | RETURN |
| | 220 | END |
| | 221 | |
| | 222 | |
| | 223 | SUBROUTINE FSIPN |
| | 224 | C calculating particle-nucleus Coulomb Wave functions FFij |
| | 225 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 226 | COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2 |
| | 227 | COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, !part. momenta in NRF |
| | 228 | 1 P2X,P2Y,P2Z,E2,P2 |
| | 229 | COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emis. points in NRF |
| | 230 | 1 X2,Y2,Z2,T2,R2 |
| | 231 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| | 232 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| | 233 | COMMON/FSI_ICH1/ICH1 |
| | 234 | COMMON/FSI_ETA/ETA |
| | 235 | COMMON/FSI_FFPN/FF12,FF21 |
| | 236 | COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS |
| | 237 | |
| | 238 | COMPLEX*16 FF1,FF12,FF21 |
| | 239 | FF12=DCMPLX(1.D0,0.D0) |
| | 240 | FF21=DCMPLX(1.D0,0.D0) |
| | 241 | IF(I3C.EQ.0)RETURN |
| | 242 | ICH1=IDINT(C1) |
| | 243 | IF(ICH1.EQ.0)GOTO 11 |
| | 244 | XH=AC1*P1 |
| | 245 | ACH=ACP(XH) |
| | 246 | ACHR=DSQRT(ACH) |
| | 247 | ETA=0.D0 |
| | 248 | IF(XH.NE.0.D0)ETA=1/XH |
| | 249 | RHOS=P1*R1 |
| | 250 | HS=X1*P1X+Y1*P1Y+Z1*P1Z |
| | 251 | FF12=FF12*FF1(RHOS,HS) |
| | 252 | IF(IQS.EQ.0)GOTO 11 |
| | 253 | RHOS=P1*R2 |
| | 254 | HS=X2*P1X+Y2*P1Y+Z2*P1Z |
| | 255 | FF21=FF21*FF1(RHOS,HS) |
| | 256 | 11 ICH1=IDINT(C2) |
| | 257 | IF(ICH1.EQ.0)GOTO 10 |
| | 258 | XH=AC2*P2 |
| | 259 | ACH=ACP(XH) |
| | 260 | ACHR=DSQRT(ACH) |
| | 261 | ETA=0.D0 |
| | 262 | IF(XH.NE.0.D0)ETA=1/XH |
| | 263 | RHOS=P2*R2 |
| | 264 | HS=X2*P2X+Y2*P2Y+Z2*P2Z |
| | 265 | FF12=FF12*FF1(RHOS,HS) |
| | 266 | CW WRITE(6,41)'AC2 ',AC2,'ACH ',ACH,'ETA ',ETA,'RHOS ',RHOS,'HS ',HS |
| | 267 | 41 FORMAT(5(A5,E11.4)) |
| | 268 | CW WRITE(6,40)'FF12 ',DREAL(FF12),DIMAG(FF12) |
| | 269 | IF(IQS.EQ.0)GOTO 10 |
| | 270 | RHOS=P2*R1 |
| | 271 | HS=X1*P2X+Y1*P2Y+Z1*P2Z |
| | 272 | FF21=FF21*FF1(RHOS,HS) |
| | 273 | CW WRITE(6,41)'AC1 ',AC1,'ACH ',ACH,'ETA ',ETA,'RHOS ',RHOS,'HS ',HS |
| | 274 | CW WRITE(6,40)'FF21 ',DREAL(FF21),DIMAG(FF21) |
| | 275 | 40 FORMAT(A7,2E12.4) |
| | 276 | 10 CONTINUE |
| | 277 | |
| | 278 | C WEIF = the weight due to the Coulomb particle-nucleus interaction |
| | 279 | WEIF=DREAL(FF12)**2+DIMAG(FF12)**2 |
| | 280 | IF(IQS.EQ.1)WEIF=0.5D0*(WEIF+DREAL(FF21)**2+DIMAG(FF21)**2) |
| | 281 | RETURN |
| | 282 | END |
| | 283 | |
| | 284 | FUNCTION GPIPI(X,J) |
| | 285 | C--- GPIPI = k*COTG(DELTA), X=k^2 |
| | 286 | C-- J=1(2) corresponds to isospin=0(2) |
| | 287 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 288 | COMMON/FSI_AAPI/AAPI(20,2) |
| | 289 | COMMON/FSI_C/HELP(20),AM,AMS,DM |
| | 290 | OM=DSQRT(X+AMS) |
| | 291 | XX=X/AMS |
| | 292 | GPIPI=OM/AAPI(1,J) |
| | 293 | GPIPI=GPIPI*(1+(AAPI(3,J)-AAPI(1,J)**2)*XX+AAPI(4,J)*XX*XX) |
| | 294 | GPIPI=GPIPI/(1+(AAPI(3,J)+AAPI(2,J)/AAPI(1,J))*XX) |
| | 295 | RETURN |
| | 296 | END |
| | 297 | FUNCTION GND(X,J) |
| | 298 | C--- GND = k*COTG(DELTA), X=k^2 |
| | 299 | C--- J=1(2) corresp. to nd(pd), S=1/2, |
| | 300 | C--- J=3(4) corresp. to nd(pd), S=3/2 |
| | 301 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 302 | COMMON/FSI_AAND/AAND(20,4) |
| | 303 | XX=X |
| | 304 | GND=1/AAND(1,J)+.5D0*AAND(2,J)*X |
| | 305 | DO 1 I=4,4 |
| | 306 | XX=XX*X |
| | 307 | 1 GND=GND+AAND(I,J)*XX |
| | 308 | GND=GND/(1+AAND(3,J)*X) |
| | 309 | RETURN |
| | 310 | END |
| | 311 | |
| | 312 | SUBROUTINE CKKB ! calculates KK-b scattering amplitude, |
| | 313 | ! saturated by S*(980) and delta(982) resonances |
| | 314 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 315 | COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS, |
| | 316 | 1 X,Y,Z,T,RP,RPS |
| | 317 | COMMON/FSI_AAKK/AAKK(9) |
| | 318 | COMMON/FSI_C/C(10),AM,AMS,DM |
| | 319 | COMPLEX*16 C |
| | 320 | S4=AKS+AAKK(1) |
| | 321 | S=4*S4 |
| | 322 | AKPIPI=DSQRT(S4-AAKK(2)) |
| | 323 | EETA2=(S+AAKK(3)-AAKK(2))**2/4/S |
| | 324 | AKPIETA=DSQRT(EETA2-AAKK(3)) |
| | 325 | C(1)=AAKK(6)/2/DCMPLX(AAKK(4)-S, |
| | 326 | ,-AK*AAKK(6)-AKPIPI*AAKK(7)) |
| | 327 | C(1)=C(1)+AAKK(8)/2/DCMPLX(AAKK(5)-S, |
| | 328 | ,-AK*AAKK(8)-AKPIETA*AAKK(9)) |
| | 329 | RETURN |
| | 330 | END |
| | 331 | |
| | 332 | |
| | 333 | C |
| | 334 | SUBROUTINE FSIWF |
| | 335 | C==> Prepares necessary quantities and call VZ(WEI) to calculate |
| | 336 | C the weight due to FSI |
| | 337 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 338 | COMMON/FSI_CVK/V,CVK |
| | 339 | COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, !part. momenta in NRF |
| | 340 | 1 P2X,P2Y,P2Z,E2,P2 |
| | 341 | COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS, |
| | 342 | 1 X,Y,Z,T,RP,RPS |
| | 343 | COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emis. points in NRF |
| | 344 | 1 X2,Y2,Z2,T2,R2 |
| | 345 | COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2 |
| | 346 | COMMON/FSI_SPIN/RHO(10) |
| | 347 | COMMON/FSI_BP/B,P |
| | 348 | COMMON/FSI_ETA/ETA |
| | 349 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| | 350 | COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10), |
| | 351 | 1 SBKRB(10),SDKK(10) |
| | 352 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| | 353 | COMMON/FSI_RR/F(10) |
| | 354 | COMMON/FSI_FD/FD(10),RD(10) |
| | 355 | COMMON/FSI_C/C(10),AM,AMS,DM |
| | 356 | COMPLEX*16 C,F |
| | 357 | COMMON/FSI_AA/AA |
| | 358 | COMMON/FSI_SHH/SH,CHH |
| | 359 | COMMON/FSI_AAPI/AAPI(20,2)/FSI_AAND/AAND(20,4) |
| | 360 | COMMON/FSI_P12/P12X,P12Y,P12Z,E12,P12,AM12,EPM |
| | 361 | COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS |
| | 362 | |
| | 363 | |
| | 364 | C==>calculating relative 4-coordinates of the particles in PRF |
| | 365 | C- {T,X,Y,Z} from the relative coordinates {TS,XS,YS,ZS} in NRF |
| | 366 | XS=X1-X2 |
| | 367 | YS=Y1-Y2 |
| | 368 | ZS=Z1-Z2 |
| | 369 | TS=T1-T2 |
| | 370 | RS12=XS*P12X+YS*P12Y+ZS*P12Z |
| | 371 | H1=(RS12/EPM-TS)/AM12 |
| | 372 | |
| | 373 | c-mlv X=XS+P12X*H1 |
| | 374 | c-mlv Y=YS+P12Y*H1 |
| | 375 | c-mlv Z=ZS+P12Z*H1 |
| | 376 | c-mlv T=(E12*TS-RS12)/AM12 |
| | 377 | c-mlv RPS=X*X+Y*Y+Z*Z |
| | 378 | c-mlv RP=DSQRT(RPS) |
| | 379 | c WRITE(6,38)'RP ',RP,'X ',X,Y,Z,T |
| | 380 | c38 FORMAT(A7,E11.4,A7,4E11.4) |
| | 381 | |
| | 382 | CVK=(P12X*PPX+P12Y*PPY+P12Z*PPZ)/(P12*AK) |
| | 383 | V=P12/E12 |
| | 384 | |
| | 385 | IF(ICH.EQ.0)GOTO 21 |
| | 386 | XH=AC*AK |
| | 387 | ACH=ACP(XH) |
| | 388 | ACHR=DSQRT(ACH) |
| | 389 | ETA=0.D0 |
| | 390 | IF(XH.NE.0.D0)ETA=1/XH |
| | 391 | C---HCP, HPR needed (e.g. in GST) if ICH=1 |
| | 392 | HCP=HC(XH) |
| | 393 | HPR=HCP+.1544313298D0 |
| | 394 | 21 CONTINUE |
| | 395 | MSP=MSPIN |
| | 396 | DO 30 JJ=1,MSP |
| | 397 | IF(NS.NE.1)GOTO22 |
| | 398 | C---Calc. quantities for the square well potential; |
| | 399 | C-- for LL > 5 the square well potential is not possible or available |
| | 400 | cc GAK=G(AK) |
| | 401 | BK(JJ)=DSQRT(EB(JJ)**2+AKS) |
| | 402 | XRA=2*RB(JJ)/AC |
| | 403 | HRA=BK(JJ)*RB(JJ) |
| | 404 | CALL SEQ(XRA,HRA) |
| | 405 | SBKRB(JJ)=HRA*B |
| | 406 | HRA=AK*RB(JJ) |
| | 407 | CALL GST(XRA,HRA) |
| | 408 | SDK(JJ)=SH |
| | 409 | CDK(JJ)=CHH |
| | 410 | SDKK(JJ)=RB(JJ) |
| | 411 | IF(AK.NE.0.D0)SDKK(JJ)=SH/AK |
| | 412 | IF(ICH.EQ.1)SDK(JJ)=ACH*SDK(JJ) |
| | 413 | 22 CONTINUE |
| | 414 | C----------------------------------------------------------------------- |
| | 415 | C---Calc. the strong s-wave scattering amplitude = C(JJ) |
| | 416 | C-- divided by Coulomb penetration factor squared (if ICH=1) |
| | 417 | IF(NS.NE.1)GOTO 230 |
| | 418 | IF(LL.NE.4)GOTO 230 ! SW scat. amplitude used for alfa-alfa only |
| | 419 | AKACH=AK |
| | 420 | IF(ICH.EQ.1)AKACH=AK*ACH |
| | 421 | C(JJ)=1/DCMPLX(GAK,-AKACH) ! amplitude for the SW-potential |
| | 422 | GOTO 30 |
| | 423 | 230 IF(LL.EQ.5.OR.LL.EQ.6.OR.LL.EQ.7)GOTO20 ! pipi |
| | 424 | IF(LL.EQ.8.OR.LL.EQ.9)GOTO20 ! Nd |
| | 425 | IF(LL.EQ.14.OR.LL.EQ.17.OR.LL.EQ.23)GOTO27 ! K+K-, K-p, K0K0-b |
| | 426 | A1=RD(JJ)*FD(JJ)*AKS |
| | 427 | A2=1+.5D0*A1 |
| | 428 | IF(ICH.EQ.1)A2=A2-2*HCP*FD(JJ)/AC |
| | 429 | AKF=AK*FD(JJ) |
| | 430 | IF(ICH.EQ.1)AKF=AKF*ACH |
| | 431 | C(JJ)=FD(JJ)/DCMPLX(A2,-AKF) |
| | 432 | GOTO30 |
| | 433 | 20 CONTINUE |
| | 434 | C---Calc. scatt. ampl. C(JJ) for pipi and Nd |
| | 435 | JH=LL-7+2*JJ-2 |
| | 436 | IF(LL.GT.7)GPI2=GND(AKS,JH) |
| | 437 | IF(LL.LE.7)GPI2=GPIPI(AKS,2) |
| | 438 | C(JJ)=1.D0/DCMPLX(GPI2,-AK) |
| | 439 | IF(LL.GE.7)GOTO27 |
| | 440 | GPI1=GPIPI(AKS,1) |
| | 441 | IF(LL.EQ.5)C(JJ)=.6667D0/DCMPLX(GPI1,-AK)+.3333D0*C(JJ) |
| | 442 | IF(LL.EQ.6)C(JJ)=.3333D0/DCMPLX(GPI1,-AK)+.6667D0*C(JJ) |
| | 443 | 27 CONTINUE |
| | 444 | C---Calc. K+K-, K0K0-b or K-p s-wave scatt. ampl. |
| | 445 | IF(LL.EQ.14.OR.LL.EQ.23)CALL CKKB |
| | 446 | IF(LL.EQ.17)C(JJ)=DCMPLX(3.29D0,3.55D0) |
| | 447 | C---Calc. pi+pi-, pi+pi+, pd, K+K- or K-p s-wave scatt. amplitude |
| | 448 | C-- divided by Coulomb penetration factor squared (if ICH=1) |
| | 449 | IF(ICH.EQ.0)GOTO 30 |
| | 450 | AAK=ACH*AK |
| | 451 | HCP2=2*HCP/AC |
| | 452 | C(JJ)=1/(1/C(JJ)-HCP2+DCMPLX(0.D0,AK-AAK)) |
| | 453 | 30 CONTINUE |
| | 454 | C*********************************************************************** |
| | 455 | CALL VZ |
| | 456 | RETURN |
| | 457 | END |
| | 458 | |
| | 459 | SUBROUTINE VZ |
| | 460 | C==> Calculates the weight WEI due to FSI |
| | 461 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 462 | COMMON/FSI_JR/JRAT |
| | 463 | COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS, |
| | 464 | 1 X,Y,Z,T,RP,RPS |
| | 465 | COMMON/FSI_SPIN/RHO(10) |
| | 466 | COMMON/FSI_ETA/ETA |
| | 467 | COMMON/FSI_AA/AA |
| | 468 | COMMON/FSI_FFF/F12,F21 |
| | 469 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| | 470 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| | 471 | COMMON/FSI_FD/FD(10),RD(10) |
| | 472 | COMMON/FSI_RR/F(10) |
| | 473 | COMMON/FSI_C/C(10),AM,AMS,DM |
| | 474 | COMMON/FSI_COULPH/EIDC |
| | 475 | COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS |
| | 476 | |
| | 477 | COMPLEX*16 F,C,G,PSI12,PSI21 |
| | 478 | COMPLEX*16 F12,F21 |
| | 479 | COMPLEX*16 EIDC |
| | 480 | COMPLEX*8 Z8,CGAMMA |
| | 481 | COMMON/FSI_FFPN/FF12,FF21 |
| | 482 | COMPLEX*16 FF12,FF21 |
| | 483 | WEI=0.D0 |
| | 484 | IF(JRAT.EQ.1)GOTO 11 |
| | 485 | RHOS=AK*RP |
| | 486 | HS=X*PPX+Y*PPY+Z*PPZ |
| | 487 | IF(RHOS.LT.15.D0.AND.RHOS+DABS(HS).LT.20.D0)GOTO 2 |
| | 488 | C---Calc. EIDC=exp(i*Coul.Ph.); |
| | 489 | C-- used in calc. of hypergeom. f-s in SEQA, FAS at k*R > 15, 20 |
| | 490 | Z8=CMPLX(1.,SNGL(ETA)) |
| | 491 | Z8=CGAMMA(Z8) |
| | 492 | EIDC=Z8/CABS(Z8) |
| | 493 | |
| | 494 | c write(*,*)'Z8 EIDC',Z8,EIDC |
| | 495 | |
| | 496 | 2 CALL FF(RHOS,HS) |
| | 497 | 11 MSP=MSPIN |
| | 498 | IF(ISI.EQ.0)GOTO 4 ! the strong interaction ON (OFF) if ISI=1(0) |
| | 499 | IF(RP.LT.AA)GOTO 4 |
| | 500 | IF(JRAT.NE.1) CALL FIRT |
| | 501 | IF(IQS.EQ.0)GOTO 5 ! the quantum statistics ON (OFF) if IQS=1(0) |
| | 502 | JSIGN=-1 |
| | 503 | DO 1 JJ=1,MSP |
| | 504 | JSIGN=-JSIGN |
| | 505 | G=F(JJ)*C(JJ) |
| | 506 | IF(ICH.EQ.1)G=G*ACHR |
| | 507 | PSI12=FF12*(F12+G) |
| | 508 | PSI21=FF21*(F21+G) |
| | 509 | G=PSI12+JSIGN*PSI21 |
| | 510 | 1 WEI=WEI+RHO(JJ)*(DREAL(G)**2+DIMAG(G)**2) |
| | 511 | GOTO 8 |
| | 512 | 5 DO 6 JJ=1,MSP |
| | 513 | G=F(JJ)*C(JJ) |
| | 514 | IF(ICH.EQ.1)G=G*ACHR |
| | 515 | CW WRITE(6,38)'JJ ',JJ,'F ',DREAL(F(JJ)),DIMAG(F(JJ)) |
| | 516 | CW WRITE(6,38)'JJ ',JJ,'C ',DREAL(C(JJ)),DIMAG(C(JJ)) |
| | 517 | CW WRITE(6,38)'JJ ',JJ,'G ',DREAL(G),DIMAG(G) |
| | 518 | CW WRITE(6,38)'JJ ',JJ,'F12+G ',DREAL(F12+G),DIMAG(F12+G) |
| | 519 | CW WRITE(6,38)'JJ ',JJ,'F21+G ',DREAL(F21+G),DIMAG(F21+G) |
| | 520 | 38 FORMAT(A7,I3,A7,2E11.4) |
| | 521 | PSI12=FF12*(F12+G) |
| | 522 | 6 WEI=WEI+RHO(JJ)*(DREAL(PSI12)**2+DIMAG(PSI12)**2) |
| | 523 | RETURN |
| | 524 | 4 PSI12=FF12*F12 |
| | 525 | IF(IQS.EQ.0)GOTO 50 ! the quantum statistics ON (OFF) if IQS=1(0) |
| | 526 | PSI21=FF21*F21 |
| | 527 | JSIGN=-1 |
| | 528 | DO 3 JJ=1,MSP |
| | 529 | JSIGN=-JSIGN |
| | 530 | G=PSI12+JSIGN*PSI21 |
| | 531 | 3 WEI=WEI+RHO(JJ)*(DREAL(G)**2+DIMAG(G)**2) |
| | 532 | GOTO 8 |
| | 533 | 50 WEI=DREAL(PSI12)**2+DIMAG(PSI12)**2 |
| | 534 | RETURN |
| | 535 | 8 WEI=WEI/2 |
| | 536 | |
| | 537 | c [M %/WRITE(*,*)WEI |
| | 538 | |
| | 539 | RETURN |
| | 540 | END |
| | 541 | |
| | 542 | |
| | 543 | SUBROUTINE FIRT |
| | 544 | C---CALC. THE F(JJ) |
| | 545 | C-- F(JJ)*C(JJ)= DEVIATION OF THE BETHE-SALPETER AMPL. FROM PLANE WAVE |
| | 546 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 547 | COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS, |
| | 548 | 1 X,Y,Z,T,RP,RPS |
| | 549 | COMMON/FSI_SHH/SH,CHH |
| | 550 | COMMON/FSI_BP/B,P |
| | 551 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| | 552 | COMMON/FSI_C/C(10),AM,AMS,DM |
| | 553 | COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10), |
| | 554 | 1 SBKRB(10),SDKK(10) |
| | 555 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| | 556 | COMMON/FSI_RR/F(10) |
| | 557 | EQUIVALENCE(RSS,RP),(TSS,T) |
| | 558 | COMPLEX*16 F,C,CH1 |
| | 559 | MSP=MSPIN |
| | 560 | DO 10 JJ=1,MSP |
| | 561 | IF(JJ.GT.1)GOTO 3 |
| | 562 | XRA=2*RSS/AC |
| | 563 | IF(AK.NE.0.D0)GOTO2 |
| | 564 | SHK=1.D0 |
| | 565 | SH=.0D0 |
| | 566 | SHH=SH |
| | 567 | CHH=1/RSS |
| | 568 | GOTO3 |
| | 569 | 2 H=AK*RSS |
| | 570 | CALL GST(XRA,H) |
| | 571 | SH=SH/RSS |
| | 572 | CHH=CHH/RSS |
| | 573 | SHH=SH |
| | 574 | IF(ICH.EQ.1) SHH=ACH*SH |
| | 575 | 3 IF(NS.EQ.2)GOTO1 |
| | 576 | C---F= ASYMPTOTIC FORMULA (T= 0 APPROX.); NS= 4 |
| | 577 | 6 F(JJ)=DCMPLX(CHH,SHH) |
| | 578 | IF(NS.NE.1)GOTO 10 |
| | 579 | C---F INSIDE THE SQUARE-WELL (T= 0 APPROX.); NS= 1 |
| | 580 | IF(RSS.GE.RB(JJ)) GOTO 10 |
| | 581 | IF(AK.NE.0.D0.AND.JJ.EQ.1)SHK=SH/AK |
| | 582 | H=BK(JJ)*RSS |
| | 583 | CALL GST(XRA,H) |
| | 584 | SKR=B*BK(JJ) |
| | 585 | F(JJ)=DCMPLX(CDK(JJ),SDK(JJ))*SKR |
| | 586 | CH1=(SDKK(JJ)*SKR-SHK*SBKRB(JJ))/C(JJ) |
| | 587 | F(JJ)=(F(JJ)+CH1)/SBKRB(JJ) |
| | 588 | GOTO 10 |
| | 589 | 1 CONTINUE |
| | 590 | C---F= ASYMPTOTIC FORMULA (T= 0 NOT REQUIRED); NS= 2 |
| | 591 | IF(JJ.GT.1)GOTO 8 |
| | 592 | IF(TSS.EQ.0.D0)GOTO6 |
| | 593 | TSSA=DABS(TSS) |
| | 594 | IF(DM.NE.0.D0)GOTO 11 |
| | 595 | H=AM*.5D0/TSSA |
| | 596 | IF(AK.NE.0.D0)GOTO4 |
| | 597 | HM=H*RPS |
| | 598 | IF(HM.GE.3.D15)GOTO6 |
| | 599 | FS1=DFRSIN(HM) |
| | 600 | FC1=DFRCOS(HM) |
| | 601 | FC2=FC1 |
| | 602 | FS2=FS1 |
| | 603 | GOTO5 |
| | 604 | 4 CONTINUE |
| | 605 | H1=AK*TSSA/AM |
| | 606 | HM=H*(RSS-H1)**2 |
| | 607 | HP=H*(RSS+H1)**2 |
| | 608 | IF(HP.GE.3.D15)GOTO6 |
| | 609 | FS1=DFRSIN(HM) |
| | 610 | FC1=DFRCOS(HM) |
| | 611 | FS2=DFRSIN(HP) |
| | 612 | FC2=DFRCOS(HP) |
| | 613 | GOTO 5 |
| | 614 | 11 CONTINUE |
| | 615 | FS1=0.D0 |
| | 616 | FS2=0.D0 |
| | 617 | FC1=0.D0 |
| | 618 | FC2=0.D0 |
| | 619 | DO 13 I=1,2 |
| | 620 | IF(I.EQ.1)TSSH=TSSA*(1+DM) |
| | 621 | IF(I.EQ.2)TSSH=TSSA*(1-DM) |
| | 622 | H=AM*.5D0/TSSH |
| | 623 | IF(AK.NE.0.D0)GOTO 12 |
| | 624 | HM=H*RPS |
| | 625 | IF(HM.GE.3.D15)GOTO6 |
| | 626 | FS1=FS1+DFRSIN(HM)/2 |
| | 627 | FC1=FC1+DFRCOS(HM)/2 |
| | 628 | IF(I.EQ.1)GOTO 13 |
| | 629 | FC2=FC1 |
| | 630 | FS2=FS1 |
| | 631 | GOTO 13 |
| | 632 | 12 CONTINUE |
| | 633 | H1=AK*TSSH/AM |
| | 634 | HM=H*(RSS-H1)**2 |
| | 635 | HP=H*(RSS+H1)**2 |
| | 636 | IF(HP.GE.3.D15)GOTO6 |
| | 637 | FS1=FS1+DFRSIN(HM)/2 |
| | 638 | FC1=FC1+DFRCOS(HM)/2 |
| | 639 | FS2=FS2+DFRSIN(HP)/2 |
| | 640 | FC2=FC2+DFRCOS(HP)/2 |
| | 641 | 13 CONTINUE |
| | 642 | 5 C12=FC1+FS2 |
| | 643 | S12=FS1+FC2 |
| | 644 | A12=FS1-FC2 |
| | 645 | A21=FS2-FC1 |
| | 646 | A2=.5D0*(CHH*(A12+A21)+SH*(A12-A21))+SHH |
| | 647 | A1=.5D0*(CHH*(C12+S12)+SH*(C12-S12)) |
| | 648 | F(JJ)=.3989422D0*DCMPLX(A1,A2) |
| | 649 | GOTO 10 |
| | 650 | 8 F(JJ)=F(1) |
| | 651 | 10 CONTINUE |
| | 652 | RETURN |
| | 653 | END |
| | 654 | FUNCTION EXF(X) |
| | 655 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 656 | IF(X.LT.-15.D0) GO TO 1 |
| | 657 | EXF=DEXP(X) |
| | 658 | RETURN |
| | 659 | 1 EXF=.0D0 |
| | 660 | RETURN |
| | 661 | END |
| | 662 | SUBROUTINE SEQ(X,H) |
| | 663 | C---CALC. FUNCTIONS B, P (EQS. (17) OF G-K-L-L); |
| | 664 | C-- NEEDED TO CALC. THE CONFLUENT HYPERGEOMETRIC FUNCTION GST. |
| | 665 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 666 | COMMON/FSI_BP/B,P |
| | 667 | DIMENSION BH(3),PH(3) |
| | 668 | DATA ERR/1.D-7/ |
| | 669 | BH(1)=1.D0 |
| | 670 | PH(1)=1.D0 |
| | 671 | PH(2)=.0D0 |
| | 672 | BH(2)=.5D0*X |
| | 673 | B=1+BH(2) |
| | 674 | P=1.D0 |
| | 675 | HS=H*H |
| | 676 | J=0 |
| | 677 | 2 J=J+1 |
| | 678 | BH(3)=(X*BH(2)-HS*BH(1))/((J+1)*(J+2)) |
| | 679 | PH(3)=(X*PH(2)-HS*PH(1)-(2*J+1)*X*BH(2))/(J*(J+1)) |
| | 680 | B=B+BH(3) |
| | 681 | P=P+PH(3) |
| | 682 | Z=DABS(BH(2))+DABS(BH(3))+DABS(PH(2))+DABS(PH(3)) |
| | 683 | IF(Z.LT.ERR)RETURN |
| | 684 | BH(1)=BH(2) |
| | 685 | BH(2)=BH(3) |
| | 686 | PH(1)=PH(2) |
| | 687 | PH(2)=PH(3) |
| | 688 | GOTO 2 |
| | 689 | END |
| | 690 | SUBROUTINE SEQA(X,H) |
| | 691 | C---CALC. FUNCTIONS CHH=REAL(GST), SH=IMAG(GST)/ACH, B=SH/H |
| | 692 | C-- IN THE ASYMPTOTIC REGION H=K*R >> 1. |
| | 693 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 694 | COMMON/FSI_BP/B,P |
| | 695 | COMMON/FSI_SHH/SH,CHH |
| | 696 | COMMON/FSI_ETA/ETA |
| | 697 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| | 698 | COMMON/FSI_COULPH/EIDC |
| | 699 | COMPLEX*16 EIDC,GST |
| | 700 | ARG=H-ETA*DLOG(2*H) |
| | 701 | GST=DCMPLX(DCOS(ARG),DSIN(ARG)) |
| | 702 | GST=ACHR*EIDC*GST |
| | 703 | CHH=DREAL(GST) |
| | 704 | SH=DIMAG(GST)/ACH |
| | 705 | B=SH/H |
| | 706 | RETURN |
| | 707 | END |
| | 708 | SUBROUTINE FF(RHO,H) |
| | 709 | C---Calc. F12, F21; |
| | 710 | C-- F12= FF0* plane wave, FF0=F*ACHR, |
| | 711 | C---F is the confluent hypergeometric function, |
| | 712 | C-- ACHR=sqrt(ACH), where ACH is the Coulomb factor |
| | 713 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 714 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| | 715 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| | 716 | COMMON/FSI_ETA/ETA |
| | 717 | COMMON/FSI_FFF/F12,F21 |
| | 718 | COMPLEX*16 FF0,F12,F21 |
| | 719 | C=DCOS(H) |
| | 720 | S=DSIN(H) |
| | 721 | F12=DCMPLX(C,-S) |
| | 722 | F21=DCMPLX(C,S) |
| | 723 | IF(ICH.EQ.0)RETURN |
| | 724 | RHOP=RHO+H |
| | 725 | RHOM=RHO-H |
| | 726 | F12=FF0(RHO,H)*F12 |
| | 727 | F21=FF0(RHO,-H)*F21 |
| | 728 | RETURN |
| | 729 | END |
| | 730 | FUNCTION FAS(RKS) |
| | 731 | C-- FAS=F*ACHR |
| | 732 | C---F is the confluent hypergeometric function at k*r >> 1 |
| | 733 | C-- ACHR=sqrt(ACH), where ACH is the Coulomb factor |
| | 734 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 735 | COMPLEX*16 FAS,EIDC,ZZ1 |
| | 736 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| | 737 | COMMON/FSI_ETA/ETA |
| | 738 | COMMON/FSI_COULPH/EIDC |
| | 739 | D1=DLOG(RKS)*ETA |
| | 740 | D2=ETA*ETA/RKS |
| | 741 | ZZ1=DCMPLX(DCOS(D1),DSIN(D1))/EIDC |
| | 742 | FAS=DCMPLX(1.D0,-D2)*ZZ1 |
| | 743 | FAS=FAS-DCMPLX(DCOS(RKS),DSIN(RKS))*ETA/RKS/ZZ1 |
| | 744 | RETURN |
| | 745 | END |
| | 746 | FUNCTION FF0(RHO,H) |
| | 747 | C-- FF0=F*ACHR |
| | 748 | C-- F is the confluent hypergeometric function |
| | 749 | C-- (Eq. (15) of G-K-L-L), F= 1 at r* << AC |
| | 750 | C-- ACHR=sqrt(ACH), where ACH is the Coulomb factor |
| | 751 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 752 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| | 753 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| | 754 | COMMON/FSI_ETA/ETA |
| | 755 | COMPLEX*16 ALF,ALF1,Z,S,A,FF0,FAS |
| | 756 | DATA ERR/1.D-5/ |
| | 757 | S=DCMPLX(1.D0,0.D0) |
| | 758 | FF0=S |
| | 759 | RHOP=RHO+H |
| | 760 | CC GOTO 5 ! rejects the approx. calcul. of hyperg. f-ion F |
| | 761 | IF(RHOP.LT.20.D0)GOTO5 |
| | 762 | FF0=FAS(RHOP) ! approx. calc. |
| | 763 | RETURN |
| | 764 | 5 ALF=DCMPLX(.0D0,-ETA) |
| | 765 | ALF1=ALF-1 |
| | 766 | Z=DCMPLX(.0D0,RHOP) |
| | 767 | J=0 |
| | 768 | 3 J=J+1 |
| | 769 | A=(ALF1+J)/(J*J) |
| | 770 | S=S*A*Z |
| | 771 | FF0=FF0+S |
| | 772 | ZR=DABS(DREAL(S)) |
| | 773 | ZI=DABS(DIMAG(S)) |
| | 774 | IF((ZR+ZI).GT.ERR)GOTO3 |
| | 775 | FF0=FF0*ACHR |
| | 776 | RETURN |
| | 777 | END |
| | 778 | FUNCTION HC(XA) |
| | 779 | C---HC = h-function of Landau-Lifshitz: h(x)=Re[psi(1-i/x)]+ln(x) |
| | 780 | C-- psi(x) is the digamma function (the logarithmic derivative of |
| | 781 | C-- the gamma function) |
| | 782 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 783 | DIMENSION BN(15) |
| | 784 | DATA BN/.8333333333D-1,.8333333333D-2,.396825396825D-2, |
| | 785 | 1 .4166666667D-2,.7575757576D-2,.2109279609D-1, |
| | 786 | 2 .8333333333D-1,.4432598039D0 ,.305395433D1, |
| | 787 | 3 .2645621212D2, .2814601449D3, .3607510546D4, |
| | 788 | 4 .5482758333D5, .9749368235D6, .200526958D8/ |
| | 789 | X=DABS(XA) |
| | 790 | IF(X.LT..33D0) GOTO 1 |
| | 791 | CC IF(X.GE.3.5D0) GO TO 2 |
| | 792 | S=.0D0 |
| | 793 | N=0 |
| | 794 | 3 N=N+1 |
| | 795 | DS=1.D0/N/((N*X)**2+1) |
| | 796 | S=S+DS |
| | 797 | IF(DS.GT.0.1D-12) GOTO 3 |
| | 798 | C---Provides 7 digit accuracy |
| | 799 | HC=S-.5772156649D0+DLOG(X) |
| | 800 | RETURN |
| | 801 | CC 2 HC=1.2D0/X**2+DLOG(X)-.5772156649 D0 |
| | 802 | CC RETURN |
| | 803 | 1 X2=X*X |
| | 804 | XP=X2 |
| | 805 | HC=0.D0 |
| | 806 | IMA=9 |
| | 807 | IF(X.LT.0.1D0)IMA=3 |
| | 808 | DO 4 I=1,IMA |
| | 809 | HC=HC+XP*BN(I) |
| | 810 | 4 XP=X2*XP |
| | 811 | RETURN |
| | 812 | END |
| | 813 | FUNCTION ACP(X) |
| | 814 | C--- ACP = COULOMB PENETRATION FACTOR |
| | 815 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 816 | IF(X.LT.0.05D0.AND.X.GE.0.D0) GO TO 1 |
| | 817 | Y=6.2831853D0/X |
| | 818 | ACP=Y/(EXF(Y)-1) |
| | 819 | RETURN |
| | 820 | 1 ACP=1.D-6 |
| | 821 | RETURN |
| | 822 | END |
| | 823 | SUBROUTINE GST(X,H) |
| | 824 | C---CALC. THE CONFL. HYPERGEOM. F-N = CHH+i*SH |
| | 825 | C-- AND THE COULOMB F-S B, P (CALLS SEQ OR SEQA). |
| | 826 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 827 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| | 828 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| | 829 | COMMON/FSI_SHH/SH,CHH |
| | 830 | COMMON/FSI_BP/B,P |
| | 831 | 1 IF(ICH.EQ.1)GOTO 2 |
| | 832 | 3 SH=DSIN(H) |
| | 833 | CHH=DCOS(H) |
| | 834 | B=1.D0 |
| | 835 | IF(H.NE.0.D0)B=SH/H |
| | 836 | P=CHH |
| | 837 | RETURN |
| | 838 | 2 CONTINUE |
| | 839 | IF(H.GT.15.D0)GOTO4 ! comment out if you want to reject |
| | 840 | ! the approximate calculation of hyperg. f-ion G |
| | 841 | CALL SEQ(X,H) ! exact calculation |
| | 842 | SH=H*B |
| | 843 | CHH=P+B*X*(DLOG(DABS(X))+HPR) |
| | 844 | RETURN |
| | 845 | 4 CALL SEQA(X,H) |
| | 846 | RETURN |
| | 847 | END |
| | 848 | |
| | 849 | FUNCTION FF1(RHO,H) |
| | 850 | C---FF1=FF0; used for particle-nucleus system |
| | 851 | C-- FF0=F12*ACHR |
| | 852 | C-- F12 is the confluent hypergeometric function |
| | 853 | C-- (Eq. (15) of G-K-L-L), F12= 1 at r* << AC |
| | 854 | C-- ACHR=sqrt(ACH), where ACH is the Coulomb factor |
| | 855 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 856 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| | 857 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| | 858 | COMMON/FSI_ETA/ETA |
| | 859 | COMMON/FSI_COULPH/EIDC |
| | 860 | COMMON/FSI_ICH1/ICH1 |
| | 861 | COMPLEX*16 FF0,FF1 |
| | 862 | COMPLEX*16 EIDC |
| | 863 | COMPLEX*8 Z8,CGAMMA |
| | 864 | FF1=DCMPLX(1.D0,0.D0) |
| | 865 | IF(ICH1.EQ.0)GOTO 2 |
| | 866 | IF(RHO.LT.15.D0.AND.RHO+H.LT.20.D0)GOTO 2 |
| | 867 | C---Calc. EIDC=exp(i*Coul.Ph.); |
| | 868 | C-- used in calc. of hypergeom. f-s in SEQA, FAS at k*R > 15, 20 |
| | 869 | Z8=CMPLX(1.,SNGL(ETA)) |
| | 870 | Z8=CGAMMA(Z8) |
| | 871 | EIDC=Z8/CABS(Z8) |
| | 872 | 2 FF1=FF0(RHO,H) |
| | 873 | RETURN |
| | 874 | END |
| | 875 | |
| | 876 | FUNCTION G(AK) |
| | 877 | C---Used to calculate SW scattering amplitude for alpa-alpha system |
| | 878 | C---NOTE THAT SCATT. AMPL.= 1/CMPLX(G(AK),-AK*ACH) |
| | 879 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 880 | COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10), |
| | 881 | 1 SBKRB(10),SDKK(10) |
| | 882 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| | 883 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,JJ,MSPIN |
| | 884 | COMMON/FSI_BP/B,P/FSI_DERIV/BPR,PPR/FSI_SHH/SH,CHH |
| | 885 | HAC=.0D0 |
| | 886 | IF(ICH.EQ.0)GOTO 1 |
| | 887 | XH=AC*AK |
| | 888 | HCP=HC(XH) |
| | 889 | HPR=HCP+.1544313298D0 |
| | 890 | ACH=ACP(XH) |
| | 891 | HAC=2*HCP/AC |
| | 892 | 1 AKS=AK**2 |
| | 893 | BK(JJ)=DSQRT(AKS+EB(JJ)**2) ! kappa=kp |
| | 894 | X=2*RB(JJ)/AC |
| | 895 | H=BK(JJ)*RB(JJ) ! kp*d |
| | 896 | CALL GST(X,H) |
| | 897 | BRHO=B ! B(kp,d) |
| | 898 | SBKRB(JJ)=SH ! kp*d*B(kp,d) |
| | 899 | CALL DERIW(X,H) |
| | 900 | BRHOP=BPR ! B'(kp,d)= dB(kp,r)/dln(r) at r=d |
| | 901 | H=AK*RB(JJ) |
| | 902 | CALL GST(X,H) |
| | 903 | CDK(JJ)=CHH ! ReG(k,d) |
| | 904 | BRHOS=B ! B(k,d) |
| | 905 | PRHOS=P ! P(k,d) |
| | 906 | SDK(JJ)=SH |
| | 907 | IF(ICH.EQ.0)GOTO 2 |
| | 908 | SDK(JJ)=ACH*SH ! ImG(k,d) |
| | 909 | IF(AK.EQ.0.D0.AND.AC.LT.0.D0)SDK(JJ)=3.14159*X*B |
| | 910 | 2 SDKK(JJ)=RB(JJ) |
| | 911 | IF(AK.NE.0.D0)SDKK(JJ)=SH/AK ! d*B(k,d) |
| | 912 | CALL DERIW(X,H) ! PPR=P'(k,d)= dP(k,r)/dln(r) at r=d |
| | 913 | ZZ=PPR-PRHOS |
| | 914 | IF(ICH.EQ.1)ZZ=ZZ+X*(BRHOS+BPR*(DLOG(DABS(X))+HPR)) |
| | 915 | C ZZ= P'(k,d)-P(k,d)+x*{B(k,d)+B'(k,d)*[ln!x!+2*C-1+h(k*ac)]} |
| | 916 | GG=(BRHOP*CDK(JJ)-BRHO*ZZ)/RB(JJ) |
| | 917 | C GG= [B'(kp,d)*ReG(k,d)-B(kp,d)*ZZ]/d |
| | 918 | G=GG/(BRHO*BPR-BRHOP*BRHOS) |
| | 919 | C G= GG/[B(kp,d)*B'(k,d)-B'(kp,d)*B(k,d)] |
| | 920 | RETURN |
| | 921 | END |
| | 922 | SUBROUTINE DERIW(X,H) |
| | 923 | C---CALLED BY F-N G(AK) |
| | 924 | IMPLICIT REAL*8 (A-H,O-Z) |
| | 925 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| | 926 | COMMON/FSI_BP/B,P/FSI_DERIV/BPR,PPR |
| | 927 | HH=.1D-3 |
| | 928 | CALL GST(X,H-HH) |
| | 929 | Q1=P |
| | 930 | B1=B |
| | 931 | CALL GST(X,H+HH) |
| | 932 | HHH=HH+HH |
| | 933 | BPR=H*(B-B1)/HHH |
| | 934 | PPR=H*(P-Q1)/HHH |
| | 935 | IF(ICH.EQ.0)RETURN |
| | 936 | CALL GST(X-HH,H) |
| | 937 | Q1=P |
| | 938 | B1=B |
| | 939 | CALL GST(X+HH,H) |
| | 940 | BPR=BPR+X*(B-B1)/HHH |
| | 941 | PPR=PPR+X*(P-Q1)/HHH |
| | 942 | RETURN |
| | 943 | END |
| | 944 | C================================================================ |
| | 945 | |
| | 946 | |
| | 947 | c-------------#include "gen/pilot.h" |
| | 948 | FUNCTION DFRSIN(X) |
| | 949 | |
| | 950 | IMPLICIT DOUBLE PRECISION (A-H,O-Z) |
| | 951 | |
| | 952 | DIMENSION A(0:16),B(0:15),C1(0:25),C2(0:28) |
| | 953 | |
| | 954 | PARAMETER (Z1 = 1, R8 = Z1/8, R32 = Z1/32) |
| | 955 | |
| | 956 | DATA C0 /1.25331 41373 15500 3D0/ |
| | 957 | |
| | 958 | DATA NA,NB,NC1,NC2 /16,15,25,28/ |
| | 959 | |
| | 960 | DATA A( 0) / 0.76435 13866 41860 002D0/ |
| | 961 | DATA A( 1) /-0.43135 54754 76601 793D0/ |
| | 962 | DATA A( 2) / 0.43288 19997 97266 531D0/ |
| | 963 | DATA A( 3) /-0.26973 31033 83871 110D0/ |
| | 964 | DATA A( 4) / 0.08416 04532 08769 354D0/ |
| | 965 | DATA A( 5) /-0.01546 52448 44613 820D0/ |
| | 966 | DATA A( 6) / 0.00187 85542 34398 220D0/ |
| | 967 | DATA A( 7) /-0.00016 26497 76188 875D0/ |
| | 968 | DATA A( 8) / 0.00001 05739 76563 833D0/ |
| | 969 | DATA A( 9) /-0.00000 05360 93398 892D0/ |
| | 970 | DATA A(10) / 0.00000 00218 16584 549D0/ |
| | 971 | DATA A(11) /-0.00000 00007 29016 212D0/ |
| | 972 | DATA A(12) / 0.00000 00000 20373 325D0/ |
| | 973 | DATA A(13) /-0.00000 00000 00483 440D0/ |
| | 974 | DATA A(14) / 0.00000 00000 00009 865D0/ |
| | 975 | DATA A(15) /-0.00000 00000 00000 175D0/ |
| | 976 | DATA A(16) / 0.00000 00000 00000 003D0/ |
| | 977 | |
| | 978 | DATA B( 0) / 0.63041 40431 45705 392D0/ |
| | 979 | DATA B( 1) /-0.42344 51140 57053 335D0/ |
| | 980 | DATA B( 2) / 0.37617 17264 33436 566D0/ |
| | 981 | DATA B( 3) /-0.16249 48915 45095 674D0/ |
| | 982 | DATA B( 4) / 0.03822 25577 86330 087D0/ |
| | 983 | DATA B( 5) /-0.00564 56347 71321 909D0/ |
| | 984 | DATA B( 6) / 0.00057 45495 19768 974D0/ |
| | 985 | DATA B( 7) /-0.00004 28707 15321 020D0/ |
| | 986 | DATA B( 8) / 0.00000 24512 07499 233D0/ |
| | 987 | DATA B( 9) /-0.00000 01109 88418 409D0/ |
| | 988 | DATA B(10) / 0.00000 00040 82497 317D0/ |
| | 989 | DATA B(11) /-0.00000 00001 24498 302D0/ |
| | 990 | DATA B(12) / 0.00000 00000 03200 484D0/ |
| | 991 | DATA B(13) /-0.00000 00000 00070 324D0/ |
| | 992 | DATA B(14) / 0.00000 00000 00001 336D0/ |
| | 993 | DATA B(15) /-0.00000 00000 00000 022D0/ |
| | 994 | |
| | 995 | DATA C1( 0) / 0.99056 04793 73497 549D0/ |
| | 996 | DATA C1( 1) /-0.01218 35098 31478 997D0/ |
| | 997 | DATA C1( 2) /-0.00248 27428 23113 060D0/ |
| | 998 | DATA C1( 3) / 0.00026 60949 52647 247D0/ |
| | 999 | DATA C1( 4) /-0.00000 10790 68987 406D0/ |
| | 1000 | DATA C1( 5) /-0.00000 48836 81753 933D0/ |
| | 1001 | DATA C1( 6) / 0.00000 09990 55266 368D0/ |
| | 1002 | DATA C1( 7) /-0.00000 00750 92717 372D0/ |
| | 1003 | DATA C1( 8) /-0.00000 00190 79487 573D0/ |
| | 1004 | DATA C1( 9) / 0.00000 00090 90797 293D0/ |
| | 1005 | DATA C1(10) /-0.00000 00019 66236 033D0/ |
| | 1006 | DATA C1(11) / 0.00000 00001 64772 911D0/ |
| | 1007 | DATA C1(12) / 0.00000 00000 63079 714D0/ |
| | 1008 | DATA C1(13) /-0.00000 00000 36432 219D0/ |
| | 1009 | DATA C1(14) / 0.00000 00000 10536 930D0/ |
| | 1010 | DATA C1(15) /-0.00000 00000 01716 438D0/ |
| | 1011 | DATA C1(16) /-0.00000 00000 00107 124D0/ |
| | 1012 | DATA C1(17) / 0.00000 00000 00204 099D0/ |
| | 1013 | DATA C1(18) /-0.00000 00000 00090 064D0/ |
| | 1014 | DATA C1(19) / 0.00000 00000 00025 506D0/ |
| | 1015 | DATA C1(20) /-0.00000 00000 00004 036D0/ |
| | 1016 | DATA C1(21) /-0.00000 00000 00000 570D0/ |
| | 1017 | DATA C1(22) / 0.00000 00000 00000 762D0/ |
| | 1018 | DATA C1(23) /-0.00000 00000 00000 363D0/ |
| | 1019 | DATA C1(24) / 0.00000 00000 00000 118D0/ |
| | 1020 | DATA C1(25) /-0.00000 00000 00000 025D0/ |
| | 1021 | |
| | 1022 | DATA C2( 0) / 0.04655 77987 37516 4561D0/ |
| | 1023 | DATA C2( 1) / 0.04499 21302 01239 4140D0/ |
| | 1024 | DATA C2( 2) /-0.00175 42871 39651 4532D0/ |
| | 1025 | DATA C2( 3) /-0.00014 65340 02581 0678D0/ |
| | 1026 | DATA C2( 4) / 0.00003 91330 40863 0159D0/ |
| | 1027 | DATA C2( 5) /-0.00000 34932 28659 7731D0/ |
| | 1028 | DATA C2( 6) /-0.00000 03153 53003 2345D0/ |
| | 1029 | DATA C2( 7) / 0.00000 01876 58200 8529D0/ |
| | 1030 | DATA C2( 8) /-0.00000 00377 55280 4930D0/ |
| | 1031 | DATA C2( 9) / 0.00000 00026 65516 5010D0/ |
| | 1032 | DATA C2(10) / 0.00000 00010 88144 8122D0/ |
| | 1033 | DATA C2(11) /-0.00000 00005 35500 7671D0/ |
| | 1034 | DATA C2(12) / 0.00000 00001 31576 5447D0/ |
| | 1035 | DATA C2(13) /-0.00000 00000 15286 0881D0/ |
| | 1036 | DATA C2(14) /-0.00000 00000 03394 7646D0/ |
| | 1037 | DATA C2(15) / 0.00000 00000 02702 0267D0/ |
| | 1038 | DATA C2(16) /-0.00000 00000 00946 3142D0/ |
| | 1039 | DATA C2(17) / 0.00000 00000 00207 1565D0/ |
| | 1040 | DATA C2(18) /-0.00000 00000 00012 6931D0/ |
| | 1041 | DATA C2(19) /-0.00000 00000 00013 9756D0/ |
| | 1042 | DATA C2(20) / 0.00000 00000 00008 5929D0/ |
| | 1043 | DATA C2(21) /-0.00000 00000 00003 1070D0/ |
| | 1044 | DATA C2(22) / 0.00000 00000 00000 7515D0/ |
| | 1045 | DATA C2(23) /-0.00000 00000 00000 0648D0/ |
| | 1046 | DATA C2(24) /-0.00000 00000 00000 0522D0/ |
| | 1047 | DATA C2(25) / 0.00000 00000 00000 0386D0/ |
| | 1048 | DATA C2(26) /-0.00000 00000 00000 0165D0/ |
| | 1049 | DATA C2(27) / 0.00000 00000 00000 0050D0/ |
| | 1050 | DATA C2(28) /-0.00000 00000 00000 0009D0/ |
| | 1051 | |
| | 1052 | V=ABS(X) |
| | 1053 | IF(V .LT. 8) THEN |
| | 1054 | Y=R8*V |
| | 1055 | H=2*Y**2-1 |
| | 1056 | ALFA=H+H |
| | 1057 | B1=0 |
| | 1058 | B2=0 |
| | 1059 | DO 4 I = NB,0,-1 |
| | 1060 | B0=B(I)+ALFA*B1-B2 |
| | 1061 | B2=B1 |
| | 1062 | 4 B1=B0 |
| | 1063 | H=SQRT(V)*Y*(B0-B2) |
| | 1064 | ELSE |
| | 1065 | R=1/V |
| | 1066 | H=10*R-1 |
| | 1067 | ALFA=H+H |
| | 1068 | B1=0 |
| | 1069 | B2=0 |
| | 1070 | DO 5 I = NC1,0,-1 |
| | 1071 | B0=C1(I)+ALFA*B1-B2 |
| | 1072 | B2=B1 |
| | 1073 | 5 B1=B0 |
| | 1074 | S=B0-H*B2 |
| | 1075 | B1=0 |
| | 1076 | B2=0 |
| | 1077 | DO 6 I = NC2,0,-1 |
| | 1078 | B0=C2(I)+ALFA*B1-B2 |
| | 1079 | B2=B1 |
| | 1080 | 6 B1=B0 |
| | 1081 | H=C0-SQRT(R)*(S*COS(V)+(B0-H*B2)*SIN(V)) |
| | 1082 | END IF |
| | 1083 | IF(X .LT. 0) H=-H |
| | 1084 | GO TO 9 |
| | 1085 | |
| | 1086 | ENTRY DFRCOS(X) |
| | 1087 | |
| | 1088 | V=ABS(X) |
| | 1089 | IF(V .LT. 8) THEN |
| | 1090 | H=R32*V**2-1 |
| | 1091 | ALFA=H+H |
| | 1092 | B1=0 |
| | 1093 | B2=0 |
| | 1094 | DO 1 I = NA,0,-1 |
| | 1095 | B0=A(I)+ALFA*B1-B2 |
| | 1096 | B2=B1 |
| | 1097 | 1 B1=B0 |
| | 1098 | H=SQRT(V)*(B0-H*B2) |
| | 1099 | ELSE |
| | 1100 | R=1/V |
| | 1101 | H=10*R-1 |
| | 1102 | ALFA=H+H |
| | 1103 | B1=0 |
| | 1104 | B2=0 |
| | 1105 | DO 2 I = NC1,0,-1 |
| | 1106 | B0=C1(I)+ALFA*B1-B2 |
| | 1107 | B2=B1 |
| | 1108 | 2 B1=B0 |
| | 1109 | S=B0-H*B2 |
| | 1110 | B1=0 |
| | 1111 | B2=0 |
| | 1112 | DO 3 I = NC2,0,-1 |
| | 1113 | B0=C2(I)+ALFA*B1-B2 |
| | 1114 | B2=B1 |
| | 1115 | 3 B1=B0 |
| | 1116 | H=C0-SQRT(R)*((B0-H*B2)*COS(V)-S*SIN(V)) |
| | 1117 | END IF |
| | 1118 | IF(X .LT. 0) H=-H |
| | 1119 | 9 DFRSIN=H |
| | 1120 | RETURN |
| | 1121 | END |
| | 1122 | |
| | 1123 | c--#include "gen/pilot.h" |
| | 1124 | FUNCTION CGAMMA(Z) |
| | 1125 | |
| | 1126 | |
| | 1127 | COMPLEX*8 CGAMMA |
| | 1128 | COMPLEX*8 Z,U,V,F,H,S |
| | 1129 | CHARACTER NAME*(*) |
| | 1130 | CHARACTER*80 ERRTXT |
| | 1131 | PARAMETER (NAME = 'CGAMMA') |
| | 1132 | |
| | 1133 | DIMENSION C(0:15) |
| | 1134 | |
| | 1135 | PARAMETER (Z1 = 1, HF = Z1/2) |
| | 1136 | |
| | 1137 | c--#if defined(CERNLIB_QF2C) |
| | 1138 | c--#include "gen/gcmpfun.inc" |
| | 1139 | c--#endif |
| | 1140 | |
| | 1141 | DATA PI /3.14159 26535 89793 24D0/ |
| | 1142 | DATA C1 /2.50662 82746 31000 50D0/ |
| | 1143 | |
| | 1144 | DATA C( 0) / 41.62443 69164 39068D0/ |
| | 1145 | DATA C( 1) /-51.22424 10223 74774D0/ |
| | 1146 | DATA C( 2) / 11.33875 58134 88977D0/ |
| | 1147 | DATA C( 3) / -0.74773 26877 72388D0/ |
| | 1148 | DATA C( 4) / 0.00878 28774 93061D0/ |
| | 1149 | DATA C( 5) / -0.00000 18990 30264D0/ |
| | 1150 | DATA C( 6) / 0.00000 00019 46335D0/ |
| | 1151 | DATA C( 7) / -0.00000 00001 99345D0/ |
| | 1152 | DATA C( 8) / 0.00000 00000 08433D0/ |
| | 1153 | DATA C( 9) / 0.00000 00000 01486D0/ |
| | 1154 | DATA C(10) / -0.00000 00000 00806D0/ |
| | 1155 | DATA C(11) / 0.00000 00000 00293D0/ |
| | 1156 | DATA C(12) / -0.00000 00000 00102D0/ |
| | 1157 | DATA C(13) / 0.00000 00000 00037D0/ |
| | 1158 | DATA C(14) / -0.00000 00000 00014D0/ |
| | 1159 | DATA C(15) / 0.00000 00000 00006D0/ |
| | 1160 | |
| | 1161 | c----#if !defined(CERNLIB_QF2C) |
| | 1162 | c----#include "gen/gcmpfun.inc" |
| | 1163 | c----#endif |
| | 1164 | |
| | 1165 | COMPLEX*8 GREAL,GIMAG,XARG,YARG |
| | 1166 | |
| | 1167 | COMPLEX*8 ZARG,GCONJG,GCMPLX |
| | 1168 | GREAL( ZARG)=REAL( ZARG) |
| | 1169 | GIMAG( ZARG)=AIMAG(ZARG) |
| | 1170 | GCONJG(ZARG)=CONJG(ZARG) |
| | 1171 | c-- GCMPLX(XARG,YARG)= CMPLX(XARG,YARG) |
| | 1172 | |
| | 1173 | U=Z |
| | 1174 | X=U |
| | 1175 | IF(GIMAG(U) .EQ. 0 .AND. -ABS(X) .EQ. INT(X)) THEN |
| | 1176 | F=0 |
| | 1177 | H=0 |
| | 1178 | WRITE(ERRTXT,101) X |
| | 1179 | c- CALL MTLPRT(NAME,'C305.1',ERRTXT) |
| | 1180 | ELSE |
| | 1181 | IF(X .GE. 1) THEN |
| | 1182 | F=1 |
| | 1183 | V=U |
| | 1184 | ELSEIF(X .GE. 0) THEN |
| | 1185 | F=1/U |
| | 1186 | V=1+U |
| | 1187 | ELSE |
| | 1188 | F=1 |
| | 1189 | V=1-U |
| | 1190 | END IF |
| | 1191 | H=1 |
| | 1192 | S=C(0) |
| | 1193 | DO 1 K = 1,15 |
| | 1194 | H=((V-K)/(V+(K-1)))*H |
| | 1195 | 1 S=S+C(K)*H |
| | 1196 | H=V+(4+HF) |
| | 1197 | H=C1*EXP((V-HF)*LOG(H)-H)*S |
| | 1198 | IF(X .LT. 0) H=PI/(SIN(PI*U)*H) |
| | 1199 | ENDIF |
| | 1200 | |
| | 1201 | c----#if !defined(CERNLIB_DOUBLE) |
| | 1202 | CGAMMA=F*H |
| | 1203 | c----#endif |
| | 1204 | |
| | 1205 | RETURN |
| | 1206 | 101 FORMAT('ARGUMENT EQUALS NON-POSITIVE INTEGER = ',1P,E15.1) |
| | 1207 | END |
| | 1208 | |
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental