| 37e4484e |
1 | c 1 2 3 4 5 6 7 8 |
| 2 | c---------|---------|---------|---------|---------|---------|---------|---------| |
| 3 | *-- Author : R.Lednicky 20/01/95 |
| 4 | SUBROUTINE FSIW(J,WEIF,WEI,WEIN) |
| 5 | |
| 6 | C======================================================================= |
| 7 | C Calculates final state interaction (FSI) weights |
| 8 | C WEIF = weight due to particle - (effective) nucleus FSI (p-N) |
| 9 | C WEI = weight due to p-p-N FSI |
| 10 | C WEIN = weight due to p-p FSI; note that WEIN=WEI if I3C=0; |
| 11 | C note that if I3C=1 the calculation of |
| 12 | C WEIN can be skipped by putting J=0 |
| 13 | C....................................................................... |
| 14 | C Correlation Functions: |
| 15 | C CF(p-p-N) = sum(WEI)/sum(WEIF) |
| 16 | C CF(p-p) = sum(WEIN)/sum(1); here the nucleus is completely |
| 17 | C inactive |
| 18 | C CF(p-p-"N") = sum(WEIN*WEIF')/sum(WEIF'), where WEIN and WEIF' |
| 19 | C are not correlated (calculated at different emission |
| 20 | C points, e.g., for different events); |
| 21 | C thus here the nucleus affects one-particle |
| 22 | C spectra but not the correlation |
| 23 | C....................................................................... |
| 24 | C User must supply data file <fn> on unit NUNIT (e.g. =11) specifying |
| 25 | C LL : particle pair |
| 26 | C NS : approximation used to calculate Bethe-Salpeter amplitude |
| 27 | C ITEST: test switch |
| 28 | C If ITEST=1 then also following parameters are required |
| 29 | C ICH : 1(0) Coulomb interaction between the two particles ON (OFF) |
| 30 | C IQS : 1(0) quantum statistics for the two particles ON (OFF) |
| 31 | C ISI : 1(0) strong interaction between the two particles ON (OFF) |
| 32 | C I3C : 1(0) Coulomb interaction with residual nucleus ON (OFF) |
| 33 | C This data file can contain other information useful for the user. |
| 34 | C It is read by subroutines READINT4 and READREA8(4) (or READ_FILE). |
| 35 | C---------------------------------------------------------------------- |
| 36 | C- LL 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
| 37 | C- part. 1: n p n a pi+ pi0 pi+ n p pi+ pi+ pi+ pi- K+ K+ K+ K- |
| 38 | C- part. 2: n p p a pi- pi0 pi+ d d K- K+ p p K- K+ p p |
| 39 | C NS=1 y/n: + + + + + - - - - - - - - - - - - |
| 40 | C---------------------------------------------------------------------- |
| 41 | C- LL 18 19 20 21 22 23 24 25 26 27 28 29 30 |
| 42 | C- part. 1: d d t t K0 K0 d p p p n lam p |
| 43 | C- part. 2: d a t a K0 K0b t t a lam lam lam pb |
| 44 | C NS=1 y/n: - - - - - - - - - + + + - |
| 45 | C---------------------------------------------------------------------- |
| 46 | C NS=1 Square well potential, |
| 47 | C NS=3 not used |
| 48 | C NS=4 scattered wave approximated by the spherical wave, |
| 49 | C NS=2 same as NS=4 but the approx. of equal emission times in PRF |
| 50 | C not required (t=0 approx. used in all other cases). |
| 51 | C Note: if NS=2,4, the B-S amplitude diverges at zero distance r* in |
| 52 | C the two-particle c.m.s.; user can specify a cutoff AA in |
| 53 | C SUBROUTINE FSIINI, for example: |
| 54 | C IF(NS.EQ.2.OR.NS.EQ.4)AA=5.D0 !! in 1/GeV --> AA=1. fm |
| 55 | C--------------------------------------------------------------------- |
| 56 | C ITEST=1 any values of parameters ICH, IQS, ISI, I3C are allowed |
| 57 | C and should be given in data file <fn> |
| 58 | C ITEST=0 physical values of these parameters are put automatically |
| 59 | C in FSIINI (their values are not required in data file) |
| 60 | C===================================================================== |
| 61 | C At the beginning of calculation user should call FSIINI, |
| 62 | C which reads LL, NS, ITEST (and eventually ICH, IQS, ISI, I3C) |
| 63 | C and initializes various parameters. |
| 64 | C In particular the constants in |
| 65 | C COMMON/FSI_CONS/PI,PI2,SPI,DR,W |
| 66 | C may be useful for the user: |
| 67 | C W=1/.1973D0 ! from fm to 1/GeV |
| 68 | C PI=4*DATAN(1.D0) |
| 69 | C PI2=2*PI |
| 70 | C SPI=DSQRT(PI) |
| 71 | C DR=180.D0/PI ! from radian to degree |
| 72 | C _______________________________________________________ |
| 73 | C !! |Important note: all real quantities are assumed REAL*8 | !! |
| 74 | C ------------------------------------------------------- |
| 75 | C For each event user should fill in the following information |
| 76 | C in COMMONs (all COMMONs in FSI calculation start with FSI_): |
| 77 | C ................................................................... |
| 78 | C COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2 |
| 79 | C Only |
| 80 | C AMN = mass of the effective nucleus [GeV/c**2] |
| 81 | C CN = charge of the effective nucleus [elem. charge units] |
| 82 | C are required |
| 83 | C ................................................................... |
| 84 | C COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, !part. momenta in the rest frame |
| 85 | C 1 P2X,P2Y,P2Z,E2,P2 !of effective nucleus (NRF) |
| 86 | C Only the components |
| 87 | C PiX,PiY,PiZ [GeV/c] |
| 88 | C in NRF are required. |
| 89 | C To make the corresponding Lorentz transformation user can use the |
| 90 | C subroutines LTRAN and LTRANB |
| 91 | C ................................................................... |
| 92 | C COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emission |
| 93 | C 1 X2,Y2,Z2,T2,R2 ! points in NRF |
| 94 | C The componets |
| 95 | C Xi,Yi,Zi [fm] |
| 96 | C and emission times |
| 97 | C Ti [fm/c] |
| 98 | C should be given in NRF with the origin assumed at the center |
| 99 | C of the effective nucleus. If the effect of residual nucleus is |
| 100 | C not calculated within FSIW, the NRF can be any fixed frame. |
| 101 | C----------------------------------------------------------------------- |
| 102 | C Before calling FSIW the user must call |
| 103 | C CALL LTRAN12 |
| 104 | C Besides Lorentz transformation to pair rest frame: |
| 105 | C (p1-p2)/2 --> k* it also transforms 4-coordinates of |
| 106 | C emission points from fm to 1/GeV and calculates Ei,Pi and Ri. |
| 107 | C Note that |k*|=AK in COMMON/FSI_PRF/ |
| 108 | C----------------------------------------------------------------------- |
| 109 | C After making some additional filtering using k* (say k* < k*max) |
| 110 | C or direction of vector k*, |
| 111 | C user can finally call FSIW to calculate the FSI weights |
| 112 | C to be used to construct the correlation function |
| 113 | C======================================================================= |
| 114 | |
| 115 | IMPLICIT REAL*8 (A-H,O-Z) |
| 116 | COMMON/FSI_JR/JRAT |
| 117 | COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2 |
| 118 | COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, ! particle momenta in NRF |
| 119 | 1 P2X,P2Y,P2Z,E2,P2 |
| 120 | COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS, ! k*=(p1-p2)/2 and x1-x2 |
| 121 | 1 X,Y,Z,T,RP,RPS ! in pair rest frame (PRF) |
| 122 | COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, !4-coord. of emis. points in NRF |
| 123 | 1 X2,Y2,Z2,T2,R2 |
| 124 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| 125 | COMMON/FSI_FFPN/FF12,FF21 |
| 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 | CALL FSIPN(WEIF) !weight due to particle-nucleus FSI |
| 136 | JRAT=0 |
| 137 | call BoostToPrf() |
| 138 | CALL FSIWF(WEI) !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(WEIN) ! weight due to particle-particle FSI |
| 145 | ENDIF |
| 146 | RETURN |
| 147 | END |
| 148 | C======================================================================= |
| 149 | |
| 150 | SUBROUTINE LTRAN(P0,P,PS) |
| 151 | C==>calculating particle 4-momentum PS={PSX,PSY,PSZ,ES} |
| 152 | C in rest frame of a system 0 with 4-momentum P0={P0X,P0Y,P0Z,E0} |
| 153 | C from its 4-momentum P={PX,PY,PZ,E} |
| 154 | |
| 155 | IMPLICIT REAL*8 (A-H,O-Z) |
| 156 | DIMENSION P0(4),P(4),PS(4) |
| 157 | C----------------------------------------------------------------------- |
| 158 | P0S=P0(1)**2+P0(2)**2+P0(3)**2 |
| 159 | AM0=DSQRT(P0(4)**2-P0S) |
| 160 | EPM=P0(4)+AM0 |
| 161 | PP0=P(1)*P0(1)+P(2)*P0(2)+P(3)*P0(3) |
| 162 | H=(PP0/EPM-P(4))/AM0 |
| 163 | PS(1)=P(1)+P0(1)*H |
| 164 | PS(2)=P(2)+P0(2)*H |
| 165 | PS(3)=P(3)+P0(3)*H |
| 166 | PS(4)=(P0(4)*P(4)-PP0)/AM0 |
| 167 | RETURN |
| 168 | END |
| 169 | |
| 170 | SUBROUTINE LTRANB(P0,PS,P) |
| 171 | C==>calculating particle 4-momentum P={PX,PY,PZ,E} |
| 172 | C from its 4-momentum PS={PSX,PSY,PSZ,ES} |
| 173 | C in rest frame of a system 0 with 4-momentum P0={P0X,P0Y,P0Z,E0} |
| 174 | |
| 175 | IMPLICIT REAL*8 (A-H,O-Z) |
| 176 | DIMENSION P0(4),P(4),PS(4) |
| 177 | C----------------------------------------------------------------------- |
| 178 | P0S=P0(1)**2+P0(2)**2+P0(3)**2 |
| 179 | AM0=DSQRT(P0(4)**2-P0S) |
| 180 | EPM=P0(4)+AM0 |
| 181 | PSP0=PS(1)*P0(1)+PS(2)*P0(2)+PS(3)*P0(3) |
| 182 | HS=(PSP0/EPM+PS(4))/AM0 |
| 183 | P(1)=PS(1)+P0(1)*HS |
| 184 | P(2)=PS(2)+P0(2)*HS |
| 185 | P(3)=PS(3)+P0(3)*HS |
| 186 | P(4)=(P0(4)*PS(4)+PSP0)/AM0 |
| 187 | RETURN |
| 188 | END |
| 189 | |
| 190 | SUBROUTINE LTRAN12 |
| 191 | C==>calculating particle momentum in PRF {EE,PPX,PPY,PPZ} from |
| 192 | C- the momentum of the first particle {E1,P1X,P1Y,P1Z) in NRF |
| 193 | IMPLICIT REAL*8 (A-H,O-Z) |
| 194 | COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, !part. momenta in NRF |
| 195 | 1 P2X,P2Y,P2Z,E2,P2 |
| 196 | COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS, |
| 197 | 1 X,Y,Z,T,RP,RPS |
| 198 | COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2 |
| 199 | COMMON/FSI_P12/P12X,P12Y,P12Z,E12,P12,AM12,EPM |
| 200 | COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, !4-coord. of emis. points in NRF |
| 201 | 1 X2,Y2,Z2,T2,R2 |
| 202 | COMMON/FSI_CONS/PI,PI2,SPI,DR,W |
| 203 | |
| 204 | C fm --> 1/GeV |
| 205 | X1=X1*W |
| 206 | Y1=Y1*W |
| 207 | Z1=Z1*W |
| 208 | T1=T1*W |
| 209 | X2=X2*W |
| 210 | Y2=Y2*W |
| 211 | Z2=Z2*W |
| 212 | T2=T2*W |
| 213 | C calculating Ri, Pi and Ei |
| 214 | R1=DSQRT(X1*X1+Y1*Y1+Z1*Z1) |
| 215 | R2=DSQRT(X2*X2+Y2*Y2+Z2*Z2) |
| 216 | P1S=P1X*P1X+P1Y*P1Y+P1Z*P1Z |
| 217 | P2S=P2X*P2X+P2Y*P2Y+P2Z*P2Z |
| 218 | P1=DSQRT(P1S) |
| 219 | P2=DSQRT(P2S) |
| 220 | E1=DSQRT(AM1*AM1+P1S) |
| 221 | E2=DSQRT(AM2*AM2+P2S) |
| 222 | C----------------------------------------------------------------------- |
| 223 | E12=E1+E2 |
| 224 | P12X=P1X+P2X |
| 225 | P12Y=P1Y+P2Y |
| 226 | P12Z=P1Z+P2Z |
| 227 | P12S=P12X**2+P12Y**2+P12Z**2 |
| 228 | AM12=DSQRT(E12**2-P12S) |
| 229 | EPM=E12+AM12 |
| 230 | P12=DSQRT(P12S) |
| 231 | P112=P1X*P12X+P1Y*P12Y+P1Z*P12Z |
| 232 | H1=(P112/EPM-E1)/AM12 |
| 233 | PPX=P1X+P12X*H1 |
| 234 | PPY=P1Y+P12Y*H1 |
| 235 | PPZ=P1Z+P12Z*H1 |
| 236 | EE=(E12*E1-P112)/AM12 |
| 237 | AKS=EE**2-AM1**2 |
| 238 | AK=DSQRT(AKS) |
| 239 | |
| 240 | CW WRITE(6,38)'AK ',AK,'K ',PPX,PPY,PPZ,EE |
| 241 | 38 FORMAT(A7,E11.4,A7,4E11.4) |
| 242 | RETURN |
| 243 | END |
| 244 | |
| 245 | SUBROUTINE FSIPN(WEIF) |
| 246 | C calculating particle-nucleus Coulomb Wave functions FFij |
| 247 | IMPLICIT REAL*8 (A-H,O-Z) |
| 248 | COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2 |
| 249 | COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, !part. momenta in NRF |
| 250 | 1 P2X,P2Y,P2Z,E2,P2 |
| 251 | COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emis. points in NRF |
| 252 | 1 X2,Y2,Z2,T2,R2 |
| 253 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| 254 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| 255 | COMMON/FSI_ICH1/ICH1 |
| 256 | COMMON/FSI_ETA/ETA |
| 257 | COMMON/FSI_FFPN/FF12,FF21 |
| 258 | COMPLEX*16 FF1,FF12,FF21 |
| 259 | FF12=DCMPLX(1.D0,0.D0) |
| 260 | FF21=DCMPLX(1.D0,0.D0) |
| 261 | C ACH=1.D0 |
| 262 | C WEIF=1.D0 |
| 263 | IF(I3C.EQ.0)RETURN |
| 264 | ICH1=IDINT(C1) |
| 265 | IF(ICH1.EQ.0)GOTO 11 |
| 266 | XH=AC1*P1 |
| 267 | ACH=ACP(XH) |
| 268 | ACHR=DSQRT(ACH) |
| 269 | ETA=0.D0 |
| 270 | IF(XH.NE.0.D0)ETA=1/XH |
| 271 | RHOS=P1*R1 |
| 272 | HS=X1*P1X+Y1*P1Y+Z1*P1Z |
| 273 | FF12=FF12*FF1(RHOS,HS) |
| 274 | IF(IQS.EQ.0)GOTO 11 |
| 275 | RHOS=P1*R2 |
| 276 | HS=X2*P1X+Y2*P1Y+Z2*P1Z |
| 277 | FF21=FF21*FF1(RHOS,HS) |
| 278 | 11 ICH1=IDINT(C2) |
| 279 | IF(ICH1.EQ.0)GOTO 10 |
| 280 | XH=AC2*P2 |
| 281 | ACH=ACP(XH) |
| 282 | ACHR=DSQRT(ACH) |
| 283 | ETA=0.D0 |
| 284 | IF(XH.NE.0.D0)ETA=1/XH |
| 285 | RHOS=P2*R2 |
| 286 | HS=X2*P2X+Y2*P2Y+Z2*P2Z |
| 287 | FF12=FF12*FF1(RHOS,HS) |
| 288 | CW WRITE(6,41)'AC2 ',AC2,'ACH ',ACH,'ETA ',ETA,'RHOS ',RHOS,'HS ',HS |
| 289 | 41 FORMAT(5(A5,E11.4)) |
| 290 | CW WRITE(6,40)'FF12 ',DREAL(FF12),DIMAG(FF12) |
| 291 | IF(IQS.EQ.0)GOTO 10 |
| 292 | RHOS=P2*R1 |
| 293 | HS=X1*P2X+Y1*P2Y+Z1*P2Z |
| 294 | FF21=FF21*FF1(RHOS,HS) |
| 295 | CW WRITE(6,41)'AC1 ',AC1,'ACH ',ACH,'ETA ',ETA,'RHOS ',RHOS,'HS ',HS |
| 296 | CW WRITE(6,40)'FF21 ',DREAL(FF21),DIMAG(FF21) |
| 297 | 40 FORMAT(A7,2E12.4) |
| 298 | 10 CONTINUE |
| 299 | |
| 300 | C WEIF = the weight due to the Coulomb particle-nucleus interaction |
| 301 | WEIF=DREAL(FF12)**2+DIMAG(FF12)**2 |
| 302 | IF(IQS.EQ.1)WEIF=0.5D0*(WEIF+DREAL(FF21)**2+DIMAG(FF21)**2) |
| 303 | RETURN |
| 304 | END |
| 305 | |
| 306 | FUNCTION GPIPI(X,J) |
| 307 | C--- GPIPI = k*COTG(DELTA), X=k^2 |
| 308 | C-- J=1(2) corresponds to isospin=0(2) |
| 309 | IMPLICIT REAL*8 (A-H,O-Z) |
| 310 | COMMON/FSI_AAPI/AAPI(20,2) |
| 311 | COMMON/FSI_C/HELP(20),AM,AMS,DM |
| 312 | OM=DSQRT(X+AMS) |
| 313 | XX=X/AMS |
| 314 | GPIPI=OM/AAPI(1,J) |
| 315 | GPIPI=GPIPI*(1+(AAPI(3,J)-AAPI(1,J)**2)*XX+AAPI(4,J)*XX*XX) |
| 316 | GPIPI=GPIPI/(1+(AAPI(3,J)+AAPI(2,J)/AAPI(1,J))*XX) |
| 317 | RETURN |
| 318 | END |
| 319 | |
| 320 | FUNCTION GPIN(X,J) |
| 321 | C--- GPIN = k*COTG(DELTA), X=k^2 |
| 322 | C-- J=1(2) corresponds to piN isospin=1/2(3/2) |
| 323 | IMPLICIT REAL*8 (A-H,O-Z) |
| 324 | COMMON/FSI_AAPIN/AAPIN(20,2) |
| 325 | GPIN=1/AAPIN(1,J)+.5D0*AAPIN(2,J)*X |
| 326 | RETURN |
| 327 | END |
| 328 | |
| 329 | FUNCTION GND(X,J) |
| 330 | C--- GND = k*COTG(DELTA), X=k^2 |
| 331 | C--- J=1(2) corresp. to nd(pd), S=1/2, |
| 332 | C--- J=3(4) corresp. to nd(pd), S=3/2 |
| 333 | IMPLICIT REAL*8 (A-H,O-Z) |
| 334 | COMMON/FSI_AAND/AAND(20,4) |
| 335 | XX=X |
| 336 | GND=1/AAND(1,J)+.5D0*AAND(2,J)*X |
| 337 | DO 1 I=4,4 |
| 338 | XX=XX*X |
| 339 | 1 GND=GND+AAND(I,J)*XX |
| 340 | GND=GND/(1+AAND(3,J)*X) |
| 341 | RETURN |
| 342 | END |
| 343 | FUNCTION GDD(X,J) |
| 344 | C--- GDD = k*COTG(DELTA), X=k^2 |
| 345 | C--- J=1,2,3 corresp. to S=0,1,2 |
| 346 | IMPLICIT REAL*8 (A-H,O-Z) |
| 347 | COMMON/FSI_AADD/AADD(20,3) |
| 348 | COMMON/FSI_C/C(10),AM,AMS,DM |
| 349 | COMMON/FSI_CONS/PI,PI2,SPI,DR,W |
| 350 | COMPLEX*16 C |
| 351 | E=X/2/AM |
| 352 | ER=DSQRT(E) |
| 353 | IF(J.EQ.1)THEN |
| 354 | GDD=ER*(AADD(1,1)*DEXP(-E/AADD(2,1))-AADD(3,1)) |
| 355 | GDD=GDD/DR ! from degree to radian |
| 356 | TAND=DTAN(GDD) |
| 357 | IF(TAND.EQ.0.D0)TAND=1.D-10 |
| 358 | GDD=DSQRT(X)/TAND |
| 359 | END IF |
| 360 | IF(J.EQ.2)THEN |
| 361 | GDD=1.D10 |
| 362 | END IF |
| 363 | IF(J.EQ.3)THEN |
| 364 | GDD=ER*(AADD(1,3)+AADD(2,3)*E) |
| 365 | GDD=GDD/DR ! from degree to radian |
| 366 | TAND=DTAN(GDD) |
| 367 | IF(TAND.EQ.0.D0)TAND=1.D-10 |
| 368 | GDD=DSQRT(X)/TAND |
| 369 | END IF |
| 370 | RETURN |
| 371 | END |
| 372 | BLOCK DATA |
| 373 | IMPLICIT REAL*8 (A-H,O-Z) |
| 374 | COMMON/FSI_AAPI/AAPI(20,2)/FSI_AAND/AAND(20,4) |
| 375 | COMMON/FSI_AADD/AADD(20,3)/FSI_AAPIN/AAPIN(20,2) |
| 376 | COMMON/FSI_AAKK/AAKK(9)/FSI_AAPAP/AAPAPR(3,2),AAPAPI(3,2) |
| 377 | C---Parameters for GPIPI (I,J), J=1,2 -> isospin=0,2 |
| 378 | DATA AAPI/.2600D00, .2500D00, .3480D00,-.0637D00, 16*0.D0, |
| 379 | 1 -.0280D00,-.0820D00, .2795D00,-.0086D00, 16*0.D0/ |
| 380 | C---Parameters for GPIN (I,J), J=1,2 -> piN isospin=1/2,3/2 |
| 381 | DATA AAPIN/ .12265D1, .1563D2, 18*0.D0, |
| 382 | 1 -.750D0, .7688D2, 18*0.D0/ |
| 383 | C---Parameters for GND (I,J), J=1(2) corresp. to nd(pd), S=1/2, |
| 384 | C--- J=3(4) corresp. to nd(pd), S=3/2 |
| 385 | DATA AAND/-.3295D1,-.8819D3, .4537D4, .28645D5, 16*0.D0, |
| 386 | 1 -.13837D2, .11505D2, .0D0, .10416D2, 16*0.D0, |
| 387 | 2 -.32180D2, .1014D2, .0D0, .0D0, 16*0.D0, |
| 388 | 3 -.60213D2, .1333D2, .0D0,-.70309D2, 16*0.D0/ |
| 389 | DATA AADD/ .10617D4, .3194D-2, .56849D3, 17*0.D0, |
| 390 | 1 20*0.D0, |
| 391 | 2 -.1085D4, .1987D5, 18*0.D0/ |
| 392 | C--- AAKK= m_K^2, m_pi^2, m_eta^2, m_S*^2, m_delta^2, |
| 393 | C--- gam(S*-->KK-b), gam(S*-->pipi), gam(delta-->KK-b), |
| 394 | C--- gam(delta-->pi eta) |
| 395 | DATA AAKK/.247677D00,.01947977D00,.2997015D00,.9604D00, |
| 396 | 1 .96511D00, |
| 397 | cc 2 .792D00, .199D00, .333D00, .222D00/ ! Martin (77) |
| 398 | 2 .094D00, .110D00, .333D00, .222D00/ ! Morgan (93) |
| 399 | C---Parameters for PAP (I,J), j=1,2 -> isospin I=0,2 |
| 400 | C--- i=1-3 -> a_singlet, a_triplet, d [fm] |
| 401 | C--- Im a_IS (I=isospin, S=spin) are fixed by atomic data and |
| 402 | C n-bar survival time up to one free parameter, e.g. Im a_00 |
| 403 | C--- Batty (89), Kerbikov (93): |
| 404 | C--- Ima_10=1.96-Ima_00, Ima_01=0.367-Ima_00/3, Ima_11=0.453+Ima_00/3 |
| 405 | C--- In DATA we put Ima_00=0.3. |
| 406 | C--- Re a_IS are fixed by atomic data up to three free parameters |
| 407 | C--- Batty (89): |
| 408 | C--- Rea_aver(pp-bar)=Re[(a_00+a_10)+3(a_01+a_11)]/8=-0.9 |
| 409 | C--- In DATA we used Rea_IS from Paris potential Pignone (94) |
| 410 | C--- rescaled by 1.67 to satisfy the atomic constraint. |
| 411 | C--- Effective radius is is taken independent of IS from the phase |
| 412 | C--- shift fit by Pirner et al. (91). |
| 413 | DATA AAPAPR/-0.94D0, -1.98D0, .1D0, |
| 414 | 1 -1.40D0, 0.37D0, .1D0/ ! Re |
| 415 | DATA AAPAPI/ 0.3 D0, .267D0,-.01D0, |
| 416 | 1 1.66D0, .553D0,-.01D0/ ! Im |
| 417 | END |
| 418 | SUBROUTINE CKKB ! calculates KK-b scattering amplitude, |
| 419 | ! saturated by S*(980) and delta(982) resonances |
| 420 | IMPLICIT REAL*8 (A-H,O-Z) |
| 421 | COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS, |
| 422 | 1 X,Y,Z,T,RP,RPS |
| 423 | COMMON/FSI_AAKK/AAKK(9) |
| 424 | COMMON/FSI_C/C(10),AM,AMS,DM |
| 425 | COMPLEX*16 C |
| 426 | S4=AKS+AAKK(1) |
| 427 | S=4*S4 |
| 428 | AKPIPI=DSQRT(S4-AAKK(2)) |
| 429 | EETA2=(S+AAKK(3)-AAKK(2))**2/4/S |
| 430 | AKPIETA=DSQRT(EETA2-AAKK(3)) |
| 431 | C(1)=AAKK(6)/2/DCMPLX(AAKK(4)-S, |
| 432 | ,-AK*AAKK(6)-AKPIPI*AAKK(7)) |
| 433 | C(1)=C(1)+AAKK(8)/2/DCMPLX(AAKK(5)-S, |
| 434 | ,-AK*AAKK(8)-AKPIETA*AAKK(9)) |
| 435 | RETURN |
| 436 | END |
| 437 | SUBROUTINE CPAP ! calculates pp-bar scattering amplitude |
| 438 | ! accounting for nn-bar->pp-bar channel |
| 439 | IMPLICIT REAL*8 (A-H,O-Z) |
| 440 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| 441 | COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS, |
| 442 | 1 X,Y,Z,T,RP,RPS |
| 443 | COMMON/FSI_AAPAP/AAPAPR(3,2),AAPAPI(3,2) |
| 444 | COMMON/FSI_C/C(10),AM,AMS,DM |
| 445 | COMMON/FSI_2CHA/AK2,AK2S,AAK,HCP2,AMU2_AMU1 ! k* (kappa) for 2-nd channel |
| 446 | COMPLEX*16 C |
| 447 | DATA AM2/.93956563D0/ |
| 448 | AMU2_AMU1=AM2/AM ! AM2=2*mu(nn-bar), AM=2*mu(pp-bar) |
| 449 | AK2S0=AMU2_AMU1*AKS |
| 450 | AK2S =AK2S0-2*AM2*(AM2-AM) |
| 451 | IF(AK2S.GE.0.D0)THEN |
| 452 | AK2=DCMPLX(DSQRT(AK2S),0.D0) ! k2 |
| 453 | ELSE |
| 454 | AK2=DCMPLX(0.D0,DSQRT(-AK2S)) ! kappa2 |
| 455 | ENDIF |
| 456 | C(10)=C(6+(ISPIN-1)*2)+ |
| 457 | + DCMPLX(AAPAPR(3,ISPIN)*AKS/2-0.016D0-HCP2, |
| 458 | , AAPAPI(3,ISPIN)*AKS/2-AAK) ! (1/f)11 |
| 459 | C(5)=C(6+(ISPIN-1)*2)+ |
| 460 | + DCMPLX(AAPAPR(3,ISPIN)*AK2S0/2, |
| 461 | , AAPAPI(3,ISPIN)*AK2S0/2) |
| 462 | IF(AK2S.GE.0.D0)THEN |
| 463 | C(5)=C(5)-DCMPLX(0.D0,AK2) |
| 464 | ELSE |
| 465 | C(5)=C(5)+DCMPLX(AK2,0.D0) ! (1/f)22 |
| 466 | ENDIF |
| 467 | C(10)=C(10)*C(5)-C(7+(ISPIN-1)*2)*C(7+(ISPIN-1)*2) |
| 468 | C(ISPIN)=C(5)/C(10) ! f11 |
| 469 | C(ISPIN+2)=-C(7+(ISPIN-1)*2)/C(10) ! f12 |
| 470 | RETURN |
| 471 | END |
| 472 | |
| 473 | SUBROUTINE FSIIN(I_ITEST,I_ICH,I_IQS,I_ISI,I_I3C) |
| 474 | C SUBROUTINE FSIINI |
| 475 | C---Note: |
| 476 | C-- ICH= 0 (1) if the Coulomb interaction is absent (present); |
| 477 | C-- ISPIN= JJ= 1,2,..,MSPIN denote increasing values of the pair |
| 478 | C-- total spin S. |
| 479 | C-- To calculate the CF of two particles (with masses m1, m2 and |
| 480 | C-- charges C1, C2) the following information is required: |
| 481 | C-- AM= twice the reduced mass= 2*m1*m2/(m1+m2) in GeV/c^2, |
| 482 | C-- DM= (m1-m2)/(m1+m2), required if NS=2; |
| 483 | C-- AC= Bohr radius= 2*137.036*0.1973/(C1*C2*AMH) in fm; |
| 484 | C-- AC > 1.D9 if C1*C2= 0, AC < 0 if C1*C2 < 0; |
| 485 | C-- MSPIN= MSPINH(LL)= number of the values of the total pair spin S; |
| 486 | C-- FD= FDH(LL,JJ), RD= RDH(LL,JJ)= scattering length and effective |
| 487 | C-- radius for each value of the total pair spin S, JJ= 1,..,MSPIN; ; |
| 488 | C-- the corresponding square well parameters EB= EBH(LL,JJ), RB= |
| 489 | C-- RBH(LL,JJ) (required if NS=1) may be calculated by sear.f; |
| 490 | C-- if the effective range approximation is not valid (as is the case, |
| 491 | C-- e.g., for two-pion system) a code for calculation of the scattering |
| 492 | C-- amplitude should be supplemented; |
| 493 | C-- RHO= RHOH(LL,JJ), SF= SFH(LL,JJ), SE= SEH(LL) are spin factors; |
| 494 | C-- RHO= the probability that the spins j1 and j2 of the two particles |
| 495 | C-- will combine in a total spin S; |
| 496 | C-- RHO= (2*S+1)/[(2j1+1)*(2j2+1)] for unpolarized particles; |
| 497 | C-- RHO= (1-P1*P2)/4 and (3+P1*P2)/4 correspond to S=0 and 1 in the |
| 498 | C-- case of spin-1/2 particles with polarizations P1 and P2; |
| 499 | C----------------------------------------------------------------------- |
| 500 | IMPLICIT REAL*8 (A-H,O-Z) |
| 501 | COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2 |
| 502 | COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS, |
| 503 | 1 X,Y,Z,T,RP,RPS |
| 504 | COMMON/FSI_SPIN/RHO(10) |
| 505 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| 506 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| 507 | COMMON/FSI_FD/FD(10),RD(10) |
| 508 | COMMON/FSI_C/C(10),AM,AMS,DM |
| 509 | COMMON/FSI_CONS/PI,PI2,SPI,DR,W |
| 510 | COMPLEX*16 C |
| 511 | COMMON/FSI_AA/AA |
| 512 | COMMON/FSI_AAPI/AAPI(20,2)/FSI_AAND/AAND(20,4) |
| 513 | COMMON/FSI_AAPIN/AAPIN(20,2) |
| 514 | COMMON/FSI_AAPAP/AAPAPR(3,2),AAPAPI(3,2) |
| 515 | COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10), |
| 516 | 1 SBKRB(10),SDKK(10) |
| 517 | COMMON/FSI_FDH/FDH(30,10),RDH(30,10),EBH(30,10),RBH(30,10) |
| 518 | COMMON/FSI_RHOH/RHOH(30,10) |
| 519 | COMMON/FSI_AMCH/AM1H(30),AM2H(30),C1H(30),C2H(30),MSPINH(30) |
| 520 | C============= declarations pour l'appel de READ_FILE()============ |
| 521 | c CHARACTER*10 KEY |
| 522 | c CHARACTER*8 CH8 |
| 523 | c INTEGER*4 INT4 |
| 524 | c REAL*8 REAL8 |
| 525 | c INTEGER*4 IERR |
| 526 | C |
| 527 | C--- mass of the first and second particle |
| 528 | DATA AM1H/.93956563D0,.93827231D0,.93956563D0,3.72737978D0, |
| 529 | C .13957D0,.13498D0,.13957D0, .93956563D0, .93827231D0, |
| 530 | C 4*.13957D0,4*.493677D0, |
| 531 | C 2*1.87561339D0,2*2.80892165D0,2*.497672D0, |
| 532 | C 1.87561339D0,3*.93827231D0,.93956563D0, |
| 533 | C 1.115684D0,.93827231D0/ |
| 534 | DATA AM2H/.93956563D0,.93827231D0,.93827231D0,3.72737978D0, |
| 535 | C .13957D0,.13498D0,.13957D0, 2*1.87561339D0, |
| 536 | C 2*.493677D0,2*.93827231D0, |
| 537 | C 2*.493677D0,2*.93827231D0, |
| 538 | C 1.87561339D0,3.72737978D0,2.80892165D0,3.72737978D0, |
| 539 | C 2*.497672D0,2*2.80892165D0,3.72737978D0, |
| 540 | C 3*1.115684D0,.93827231D0/ |
| 541 | c--------|---------|---------|---------|---------|---------|---------|---------- |
| 542 | C--- charge of the first and second particle |
| 543 | DATA C1H/0.D0,1.D0,0.D0,2.D0, 1.D0,0.D0,1.D0,0.D0,1.D0, |
| 544 | C 3*1.D0,-1.D0,3*1.D0,-1.D0, |
| 545 | C 4*1.D0,2*0.D0,4*1.D0,2*0.D0, 1.D0/ |
| 546 | DATA C2H/0.D0,1.D0,1.D0,2.D0,-1.D0,0.D0,3*1.D0, |
| 547 | C -1.D0,3*1.D0,-1.D0,3*1.D0, |
| 548 | C 1.D0,2.D0,1.D0,2.D0,2*0.D0,2*1.D0,2.D0,3*0.D0,-1.D0/ |
| 549 | C---MSPIN vs (LL) |
| 550 | DATA MSPINH/3*2,4*1,2*2,8*1,3,1,2,1,2*1,2*2,1,3*2, 2/ |
| 551 | C---Spin factors <RHO vs (LL,ISPIN) |
| 552 | DATA RHOH/3*.25D0, 4*1.D0, 2*.3333D0, 8*1.D0, |
| 553 | 1 .1111D0,1.D0,.25D0,1.D0,2*1.D0, |
| 554 | 1 .3333D0,.25D0,1.D0,3*.25D0, .25D0, |
| 555 | 2 3*.75D0, 4*0.D0, 2*.6667D0, 8*0.D0, |
| 556 | 2 .3333D0,.0D0,.75D0,.0D0,2*0.D0, |
| 557 | 2 .6667D0,.75D0,0.D0,3*.75D0, .75D0, |
| 558 | 3 17*.0D0,.5556D0,3*0.D0, 8*0.D0,1*0.D0,210*0.D0/ |
| 559 | C---Scattering length FD and effective radius RD in fm vs (LL,ISPIN) |
| 560 | DATA FDH/ |
| 561 | 1 17.0D0,7.8D0,23.7D0,2230.1218D0,.225D0,.081D0,-.063D0, |
| 562 | 1 -.65D0,-2.73D0, |
| 563 | 1 .137D0,-.071D0,-.148D0,.112D0,2*1.D-6,-.360D0, |
| 564 | 1 2*1.D-6,1.344D0,6*1.D-6,-5.628D0,2.18D0,2.40D0, |
| 565 | 1 2.81D0, ! ND potential |
| 566 | C 1 0.50D0, ! NSC97e potential lam-lam |
| 567 | 1 1*0.001D0, |
| 568 | C c 2 -10.8D0,2*-5.4D0,4*0.D0,-6.35D0,-11.88D0,8*0.D0,9*0.D0, |
| 569 | 2 3*-5.4D0,4*0.D0,-6.35D0,-11.88D0,8*0.D0,9*0.D0, |
| 570 | 2 1.93D0,1.84D0, |
| 571 | 2 0.50D0, ! triplet f0 lam-lam=singlet f0 ND |
| 572 | ! not contributing in s-wave FSI approx. |
| 573 | 2 1*0.001D0, |
| 574 | 3 240*0.D0/ |
| 575 | c--------|---------|---------|---------|---------|---------|---------|---------- |
| 576 | DATA RDH/ |
| 577 | 1 2.7D0,2.8D0,2.7D0,1.12139906D0,-44.36D0,64.0D0,784.9D0, |
| 578 | 1 477.9D0, 2.27D0, 9*0.D0,-69.973D0, 6*0.D0,3.529D0, |
| 579 | 1 3.19D0,3.15D0, |
| 580 | 1 2.95D0, ! ND potential lam-lam |
| 581 | C 1 10.6D0, ! NSC97e potential lam-lam |
| 582 | 1 1*0.D0, |
| 583 | 2 3*1.7D0,4*0.D0,2.0D0,2.63D0, 17*0.D0,3.35D0,3.37D0, |
| 584 | 2 2.95D0, ! triplet d0 lam-lam=singlet d0 ND |
| 585 | ! not contributing in s-wave approx. |
| 586 | 2 1*0.D0, |
| 587 | 3 240*0.D0/ |
| 588 | C--- Corresponding square well parameters RB (width in fm) and |
| 589 | C-- EB =SQRT(-AM*U) (in GeV/c); U is the well height |
| 590 | DATA RBH/2.545739D0, 2.779789D0, 2.585795D0, 5.023544D0, |
| 591 | 1 .124673D0, .3925180D0,.09D0, 2.D0, 4.058058D0, 17*0.D0, |
| 592 | 1 2.252623D0, 2.278575D0, |
| 593 | 1 2.234089D0, ! ND potential lam-lam |
| 594 | C 1 3.065796D0, ! NSC97e potential lam-lam |
| 595 | 1 1*0.001D0, |
| 596 | 2 3*2.003144D0, |
| 597 | 2 4*0.D0, 2.D0, 4.132163D0, 17*0.D0, |
| 598 | 2 2.272703D0, 2.256355D0, |
| 599 | 2 2.234089D0, ! triplet potential lam-lam=singlet ND |
| 600 | ! not contributing in s-wave FSI approx. |
| 601 | 2 1*0.001D0, |
| 602 | 3 240*0.D0/ |
| 603 | DATA EBH/.1149517D0, .1046257D0, .1148757D0, .1186010D0, |
| 604 | 1 .7947389D0,2.281208D0,8.7D0,.4D0,.1561219D0,17*0.D0, |
| 605 | 1 .1013293D0, .1020966D0, |
| 606 | 1 .1080476D0, ! ND potential lam-lam |
| 607 | C 1 .04115994D0, ! NSC97e potential lam-lam |
| 608 | 1 1*0.001D0, |
| 609 | 2 3*.1847221D0, |
| 610 | 2 4*0.D0, .4D0, .1150687D0, 17*0.D0, |
| 611 | 2 .09736083D0, .09708310D0, |
| 612 | 2 .1080476D0, ! triplet potential lam-lam= singlet ND |
| 613 | ! not contributing in s-wave FSI approx. |
| 614 | 2 1*0.001D0, |
| 615 | 3 240*0.D0/ |
| 616 | |
| 617 | C++----- add to be able to call several time------- |
| 618 | integer ifirst |
| 619 | data ifirst/0/ |
| 620 | ifirst=ifirst+1 |
| 621 | |
| 622 | C=======< constants >======================== |
| 623 | W=1/.1973D0 ! from fm to 1/GeV |
| 624 | PI=4*DATAN(1.D0) |
| 625 | PI2=2*PI |
| 626 | SPI=DSQRT(PI) |
| 627 | DR=180.D0/PI ! from radian to degree |
| 628 | AC1=1.D10 |
| 629 | AC2=1.D10 |
| 630 | C=======< condition de calculs >============= |
| 631 | NUNIT=11 ! for IBM in Nantes |
| 632 | C NUNIT=4 ! for SUN in Prague |
| 633 | C++ CALL readint4(NUNIT,'ITEST ',ITEST) |
| 634 | C++ CALL readint4(NUNIT,'NS ',NS) |
| 635 | ITEST=I_ITEST |
| 636 | IF(ITEST.EQ.1)THEN |
| 637 | C++ CALL readint4(NUNIT,'ICH ',ICH) |
| 638 | C++ CALL readint4(NUNIT,'IQS ',IQS) |
| 639 | C++ CALL readint4(NUNIT,'ISI ',ISI) |
| 640 | C++ CALL readint4(NUNIT,'I3C ',I3C) |
| 641 | ICH=I_ICH |
| 642 | IQS=I_IQS |
| 643 | ISI=I_ISI |
| 644 | I3C=I_I3C |
| 645 | ENDIF |
| 646 | C============================================ |
| 647 | IF (IFIRST.LE.1) THEN |
| 648 | DO 3 J1=1,30 |
| 649 | DO 3 J2=1,10 |
| 650 | FDH(J1,J2)=FDH(J1,J2)*W |
| 651 | RDH(J1,J2)=RDH(J1,J2)*W |
| 652 | 3 RBH(J1,J2)=RBH(J1,J2)*W |
| 653 | C print *,"FD,RD,EB,RB: ",FDH(J1,J2),RDH(J1,J2),EBH(J1,J2),RBH(J1,J2) |
| 654 | |
| 655 | ENDIF |
| 656 | C=================================== |
| 657 | RETURN |
| 658 | END |
| 659 | C |
| 660 | SUBROUTINE LLINI(lll,I_NS,I_ITEST) |
| 661 | C===> Initialisation for a given LL value. |
| 662 | C =========================================== |
| 663 | IMPLICIT REAL*8 (A-H,O-Z) |
| 664 | COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2 |
| 665 | COMMON/FSI_SPIN/RHO(10) |
| 666 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| 667 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| 668 | COMMON/FSI_FD/FD(10),RD(10) |
| 669 | COMMON/FSI_C/C(10),AM,AMS,DM |
| 670 | COMMON/FSI_CONS/PI,PI2,SPI,DR,W |
| 671 | COMPLEX*16 C |
| 672 | COMMON/FSI_AA/AA |
| 673 | COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS, |
| 674 | 1 X,Y,Z,T,RP,RPS |
| 675 | COMMON/FSI_AAPI/AAPI(20,2)/FSI_AAND/AAND(20,4) |
| 676 | COMMON/FSI_AAPIN/AAPIN(20,2) |
| 677 | COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10), |
| 678 | 1 SBKRB(10),SDKK(10) |
| 679 | COMMON/FSI_FDH/FDH(30,10),RDH(30,10),EBH(30,10),RBH(30,10) |
| 680 | COMMON/FSI_RHOH/RHOH(30,10) |
| 681 | COMMON/FSI_AMCH/AM1H(30),AM2H(30),C1H(30),C2H(30),MSPINH(30) |
| 682 | COMMON/FSI_AAKK/AAKK(9)/FSI_AAPAP/AAPAPR(3,2),AAPAPI(3,2) |
| 683 | |
| 684 | C++----- add to be able to call several time------- |
| 685 | integer ifirst |
| 686 | data ifirst/0/ |
| 687 | ifirst=ifirst+1 |
| 688 | |
| 689 | C===> LL - Initialisation ======================================== |
| 690 | C---- Setting particle masses and charges |
| 691 | LL=lll |
| 692 | NS=I_NS |
| 693 | AM1=AM1H(LL) |
| 694 | AM2=AM2H(LL) |
| 695 | C1=C1H(LL) |
| 696 | C2=C2H(LL) |
| 697 | C print *,"llini: ",AM1,AM2,C1,C2 |
| 698 | C---> Switches: |
| 699 | C ISI=1(0) the strong interaction between the two particles ON (OFF) |
| 700 | C IQS=1(0) the quantum statistics ON (OFF); |
| 701 | C should be OFF for nonidentical particles |
| 702 | C I3C=1(0) the Coulomb interaction with the nucleus ON (OFF) |
| 703 | C I3S=1(0) the strong interaction with the nucleus ON (OFF) |
| 704 | C ICH=1(0) if C1*C2 is different from 0 (is equal to 0) |
| 705 | C- To switch off the Coulomb force between the two particles |
| 706 | C put ICH=0 and substitute the strong amplitude parameters by |
| 707 | C the ones not affected by Coulomb interaction |
| 708 | C---- -------------------------------------------------------------------- |
| 709 | IF (I_ITEST.NE.1) THEN |
| 710 | ICH=0 |
| 711 | IF(C1*C2.NE.0.D0) ICH=1 |
| 712 | IQS=0 |
| 713 | IF(C1+AM1.EQ.C2+AM2) IQS=1 |
| 714 | I3S=0 ! only this option is available |
| 715 | ISI=1 |
| 716 | I3C=1 |
| 717 | ENDIF |
| 718 | C---> Calcul. twice the reduced mass (AM), the relative |
| 719 | C mass difference (DM) and the Bohr radius (AC) |
| 720 | AM=2*AM1*AM2/(AM1+AM2) |
| 721 | AMS=AM*AM |
| 722 | DM=(AM1-AM2)/(AM1+AM2) |
| 723 | AC=1.D10 |
| 724 | C12=C1*C2 |
| 725 | IF(C12.NE.0.D0)AC=2*137.036D0/(C12*AM) |
| 726 | C---Setting spin factors |
| 727 | MSPIN=MSPINH(LL) |
| 728 | C print *,"MSPIN: ",MSPIN |
| 729 | MSP=MSPIN |
| 730 | DO 91 ISPIN=1,10 |
| 731 | 91 RHO(ISPIN)=RHOH(LL,ISPIN) |
| 732 | C 91 print *,"RHO: ",ISPIN,RHO(ISPIN) |
| 733 | C---> Integration limit AA in the spherical wave approximation |
| 734 | AA=0.D0 |
| 735 | cc IF(NS.EQ.2.OR.NS.EQ.4)AA=.5D0 !!in 1/GeV --> 0.1 fm |
| 736 | IF(NS.EQ.2.OR.NS.EQ.4)AA=6.D0 !!in 1/GeV --> 1.2 fm |
| 737 | C---> Setting scatt. length (FD), eff. radius (RD) and, if possible, |
| 738 | C-- also the corresp. square well parameters (EB, RB) |
| 739 | DO 55 JJ=1,MSP |
| 740 | ISPIN=JJ |
| 741 | FD(JJ)=FDH(LL,JJ) |
| 742 | RD(JJ)=RDH(LL,JJ) |
| 743 | EB(JJ)=EBH(LL,JJ) |
| 744 | RB(JJ)=RBH(LL,JJ) |
| 745 | C print *,"FD,RD,EB,RB: ",FD(JJ),RD(JJ),EB(JJ),RB(JJ) |
| 746 | C---Resets FD and RD for a nucleon-deuteron system (LL=8,9) |
| 747 | IF(LL.EQ.8.OR.LL.EQ.9)THEN |
| 748 | JH=LL-7+2*JJ-2 |
| 749 | FD(JJ)=AAND(1,JH) |
| 750 | RD(JJ)=AAND(2,JH)-2*AAND(3,JH)/AAND(1,JH) |
| 751 | ENDIF |
| 752 | C---Resets FD and RD for a pion-pion system (LL=5,6,7) |
| 753 | IF(LL.EQ.5.OR.LL.EQ.6.OR.LL.EQ.7)THEN |
| 754 | IF(LL.EQ.7)FD(JJ)=AAPI(1,2)/AM |
| 755 | IF(LL.EQ.5)FD(JJ)=(.6667D0*AAPI(1,1)+.3333D0*AAPI(1,2))/AM |
| 756 | IF(LL.EQ.6)FD(JJ)=(.3333D0*AAPI(1,1)+.6667D0*AAPI(1,2))/AM |
| 757 | AKS=0.D0 |
| 758 | DAKS=1.D-5 |
| 759 | AKSH=AKS+DAKS |
| 760 | AKH=DSQRT(AKSH) |
| 761 | GPI1H=GPIPI(AKSH,1) |
| 762 | GPI2H=GPIPI(AKSH,2) |
| 763 | H=1/FD(JJ) |
| 764 | IF(LL.EQ.7)C(JJ)=1/DCMPLX(GPI2H,-AKH) |
| 765 | IF(LL.EQ.5) |
| 766 | + C(JJ)=.6667D0/DCMPLX(GPI1H,-AKH)+.3333D0/DCMPLX(GPI2H,-AKH) |
| 767 | IF(LL.EQ.6) |
| 768 | + C(JJ)=.3333D0/DCMPLX(GPI1H,-AKH)+.6667D0/DCMPLX(GPI2H,-AKH) |
| 769 | HH=DREAL(1/C(JJ)) |
| 770 | RD(JJ)=2*(HH-H)/DAKS |
| 771 | ENDIF |
| 772 | C---Resets FD and RD for a pion-nucleon system (LL=12,13) |
| 773 | IF(LL.EQ.12.OR.LL.EQ.13)THEN |
| 774 | IF(LL.EQ.12)FD(JJ)=AAPIN(1,2) |
| 775 | IF(LL.EQ.13)FD(JJ)=(.6667D0*AAPIN(1,1)+.3333D0*AAPIN(1,2)) |
| 776 | AKS=0.D0 |
| 777 | DAKS=1.D-5 |
| 778 | AKSH=AKS+DAKS |
| 779 | AKH=DSQRT(AKSH) |
| 780 | GPI1H=GPIN(AKSH,1) |
| 781 | GPI2H=GPIN(AKSH,2) |
| 782 | H=1/FD(JJ) |
| 783 | IF(LL.EQ.12)C(JJ)=1/DCMPLX(GPI2H,-AKH) |
| 784 | IF(LL.EQ.13) |
| 785 | + C(JJ)=.6667D0/DCMPLX(GPI1H,-AKH)+.3333D0/DCMPLX(GPI2H,-AKH) |
| 786 | HH=DREAL(1/C(JJ)) |
| 787 | RD(JJ)=2*(HH-H)/DAKS |
| 788 | ENDIF |
| 789 | C---fm to 1/GeV for pp-bar system |
| 790 | IF(LL.EQ.30)THEN |
| 791 | IF(IFIRST.LE.1)THEN |
| 792 | DO 4 I3=1,3 |
| 793 | AAPAPR(I3,JJ)=AAPAPR(I3,JJ)*W |
| 794 | 4 AAPAPI(I3,JJ)=AAPAPI(I3,JJ)*W |
| 795 | C 4 print *,"AAPAPR,AAPAPI: ",AAPAPR(I3,JJ),AAPAPI(I3,JJ) |
| 796 | C--- Calculates complex elements M11=M22=C(6), M12=M21=C(7) for I=0 |
| 797 | C--- at k*=0 M11=M22=C(8), M12=M21=C(9) for I=1 |
| 798 | C(7+(JJ-1)*2)=2*DCMPLX(AAPAPR(1,JJ),AAPAPI(1,JJ))* |
| 799 | * DCMPLX(AAPAPR(2,JJ),AAPAPI(2,JJ)) ! 2a_0Sa_1S |
| 800 | C(6+(JJ-1)*2)=DCMPLX(AAPAPR(1,JJ)+AAPAPR(2,JJ), |
| 801 | , AAPAPI(1,JJ)+AAPAPI(2,JJ))/ |
| 802 | / C(7+(JJ-1)*2) ! M11=M22 |
| 803 | C(7+(JJ-1)*2)=-DCMPLX(AAPAPR(1,JJ)-AAPAPR(2,JJ), |
| 804 | , AAPAPI(1,JJ)-AAPAPI(2,JJ))/ |
| 805 | / C(7+(JJ-1)*2) ! M12=M21 |
| 806 | ENDIF |
| 807 | ENDIF |
| 808 | C---Calculation continues for any system (any LL) |
| 809 | 55 CONTINUE |
| 810 | RETURN |
| 811 | END |
| 812 | C======================================================= |
| 813 | C |
| 814 | |
| 815 | c++ This routine is used to init mass and charge of the nucleus. |
| 816 | |
| 817 | SUBROUTINE FSINUCL(R_AMN,R_CN) |
| 818 | |
| 819 | IMPLICIT REAL*8 (A-H,O-Z) |
| 820 | COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2 |
| 821 | |
| 822 | AMN=R_AMN |
| 823 | CN=R_CN |
| 824 | |
| 825 | RETURN |
| 826 | END |
| 827 | |
| 828 | C====================================================== |
| 829 | C |
| 830 | |
| 831 | SUBROUTINE FSIMOMENTUM(PP1,PP2) |
| 832 | |
| 833 | IMPLICIT REAL*8 (A-H,O-Z) |
| 834 | COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, ! particle momenta in NRF |
| 835 | 1 P2X,P2Y,P2Z,E2,P2 |
| 836 | |
| 837 | |
| 838 | REAL*8 PP1(3),PP2(3) |
| 839 | c Print *,"momentum",pp1,pp2 |
| 840 | P1X=PP1(1) |
| 841 | P1Y=PP1(2) |
| 842 | P1Z=PP1(3) |
| 843 | P2X=PP2(1) |
| 844 | P2Y=PP2(2) |
| 845 | P2Z=PP2(3) |
| 846 | RETURN |
| 847 | END |
| 848 | |
| 849 | |
| 850 | |
| 851 | C====================================================== |
| 852 | C |
| 853 | |
| 854 | SUBROUTINE FSIPOSITION(XT1,XT2) |
| 855 | |
| 856 | IMPLICIT REAL*8 (A-H,O-Z) |
| 857 | COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, !4-coord. of emis. points in NRF |
| 858 | 1 X2,Y2,Z2,T2,R2 |
| 859 | |
| 860 | REAL*8 XT1(4),XT2(4) |
| 861 | clc print *,'fsi',xt1,xt2 |
| 862 | X1=XT1(1) |
| 863 | Y1=XT1(2) |
| 864 | Z1=XT1(3) |
| 865 | T1=XT1(4) |
| 866 | X2=XT2(1) |
| 867 | Y2=XT2(2) |
| 868 | Z2=XT2(3) |
| 869 | T2=XT2(4) |
| 870 | RETURN |
| 871 | END |
| 872 | |
| 873 | |
| 874 | C====================================================== |
| 875 | C====================================================== |
| 876 | C |
| 877 | subroutine BoostToPrf() |
| 878 | IMPLICIT REAL*8 (A-H,O-Z) |
| 879 | COMMON/FSI_CVK/V,CVK |
| 880 | COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, !part. momenta in NRF |
| 881 | 1 P2X,P2Y,P2Z,E2,P2 |
| 882 | COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS, |
| 883 | 1 X,Y,Z,T,RP,RPS |
| 884 | COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emis. points in NRF |
| 885 | 1 X2,Y2,Z2,T2,R2 |
| 886 | COMMON/FSI_P12/P12X,P12Y,P12Z,E12,P12,AM12,EPM |
| 887 | |
| 888 | XS=X1-X2 |
| 889 | YS=Y1-Y2 |
| 890 | ZS=Z1-Z2 |
| 891 | TS=T1-T2 |
| 892 | RS12=XS*P12X+YS*P12Y+ZS*P12Z |
| 893 | H1=(RS12/EPM-TS)/AM12 |
| 894 | X=XS+P12X*H1 |
| 895 | Y=YS+P12Y*H1 |
| 896 | Z=ZS+P12Z*H1 |
| 897 | T=(E12*TS-RS12)/AM12 |
| 898 | RPS=X*X+Y*Y+Z*Z |
| 899 | RP=DSQRT(RPS) |
| 900 | CW WRITE(6,38)'RP ',RP,'X ',X,Y,Z,T |
| 901 | 38 FORMAT(A7,E11.4,A7,4E11.4) |
| 902 | |
| 903 | CVK=(P12X*PPX+P12Y*PPY+P12Z*PPZ)/(P12*AK) |
| 904 | V=P12/E12 |
| 905 | return |
| 906 | end |
| 907 | |
| 908 | SUBROUTINE FSIWF(WEI) |
| 909 | C==> Prepares necessary quantities and call VZ(WEI) to calculate |
| 910 | C the weight due to FSI |
| 911 | IMPLICIT REAL*8 (A-H,O-Z) |
| 912 | COMMON/FSI_CVK/V,CVK |
| 913 | COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, !part. momenta in NRF |
| 914 | 1 P2X,P2Y,P2Z,E2,P2 |
| 915 | COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS, |
| 916 | 1 X,Y,Z,T,RP,RPS |
| 917 | COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emis. points in NRF |
| 918 | 1 X2,Y2,Z2,T2,R2 |
| 919 | COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2 |
| 920 | COMMON/FSI_SPIN/RHO(10) |
| 921 | COMMON/FSI_BP/B,P |
| 922 | COMMON/FSI_ETA/ETA |
| 923 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| 924 | COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10), |
| 925 | 1 SBKRB(10),SDKK(10) |
| 926 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| 927 | COMMON/FSI_RR/F(10) |
| 928 | COMMON/FSI_FD/FD(10),RD(10) |
| 929 | COMMON/FSI_C/C(10),AM,AMS,DM |
| 930 | COMPLEX*16 C,F |
| 931 | COMMON/FSI_AA/AA |
| 932 | COMMON/FSI_SHH/SH,CHH |
| 933 | COMMON/FSI_AAPI/AAPI(20,2)/FSI_AAND/AAND(20,4) |
| 934 | COMMON/FSI_AAPIN/AAPIN(20,2) |
| 935 | COMMON/FSI_P12/P12X,P12Y,P12Z,E12,P12,AM12,EPM |
| 936 | COMMON/FSI_2CHA/AK2,AK2S,AAK,HCP2,AMU2_AMU1 ! k* (kappa) for 2-nd channel |
| 937 | C==>calculating relative 4-coordinates of the particles in PRF |
| 938 | C- {T,X,Y,Z} from the relative coordinates {TS,XS,YS,ZS} in NRF |
| 939 | XS=X1-X2 |
| 940 | YS=Y1-Y2 |
| 941 | ZS=Z1-Z2 |
| 942 | TS=T1-T2 |
| 943 | RS12=XS*P12X+YS*P12Y+ZS*P12Z |
| 944 | H1=(RS12/EPM-TS)/AM12 |
| 945 | X=XS+P12X*H1 |
| 946 | Y=YS+P12Y*H1 |
| 947 | Z=ZS+P12Z*H1 |
| 948 | T=(E12*TS-RS12)/AM12 |
| 949 | RPS=X*X+Y*Y+Z*Z |
| 950 | RP=DSQRT(RPS) |
| 951 | CW WRITE(6,38)'RP ',RP,'X ',X,Y,Z,T |
| 952 | 38 FORMAT(A7,E11.4,A7,4E11.4) |
| 953 | |
| 954 | CVK=(P12X*PPX+P12Y*PPY+P12Z*PPZ)/(P12*AK) |
| 955 | V=P12/E12 |
| 956 | |
| 957 | C ACH=1.D0 |
| 958 | IF(ICH.EQ.0)GOTO 21 |
| 959 | XH=AC*AK |
| 960 | ACH=ACP(XH) |
| 961 | ACHR=DSQRT(ACH) |
| 962 | ETA=0.D0 |
| 963 | IF(XH.NE.0.D0)ETA=1/XH |
| 964 | C---HCP, HPR needed (e.g. in GST) if ICH=1 |
| 965 | HCP=HC(XH) |
| 966 | HPR=HCP+.1544313298D0 |
| 967 | C AAK=ACH*AK ! |
| 968 | C HCP2=2*HCP/AC ! needed to calculate C(JJ) for charged particles |
| 969 | 21 CONTINUE |
| 970 | MSP=MSPIN |
| 971 | DO 30 JJ=1,MSP |
| 972 | ISPIN=JJ |
| 973 | IF(NS.NE.1)GOTO22 |
| 974 | C---Calc. quantities for the square well potential; |
| 975 | C-- for LL=6-26 the square well potential is not possible or available |
| 976 | IF(LL.EQ.4)GOTO 22 |
| 977 | BK(JJ)=DSQRT(EB(JJ)**2+AKS) |
| 978 | XRA=2*RB(JJ)/AC |
| 979 | HRA=BK(JJ)*RB(JJ) |
| 980 | CALL SEQ(XRA,HRA) |
| 981 | SBKRB(JJ)=HRA*B |
| 982 | HRA=AK*RB(JJ) |
| 983 | CALL GST(XRA,HRA) |
| 984 | SDK(JJ)=SH |
| 985 | CDK(JJ)=CHH |
| 986 | SDKK(JJ)=RB(JJ) |
| 987 | IF(AK.NE.0.D0)SDKK(JJ)=SH/AK |
| 988 | IF(ICH.EQ.1)SDK(JJ)=ACH*SDK(JJ) |
| 989 | 22 CONTINUE |
| 990 | C----------------------------------------------------------------------- |
| 991 | C---Calc. the strong s-wave scattering amplitude = C(JJ) |
| 992 | C-- divided by Coulomb penetration factor squared (if ICH=1) |
| 993 | IF(NS.NE.1)GOTO 230 |
| 994 | IF(LL.NE.4)GOTO 230 ! SW scat. amplitude used for alfa-alfa only |
| 995 | GAK=G(AK) |
| 996 | AKACH=AK |
| 997 | IF(ICH.EQ.1)AKACH=AK*ACH |
| 998 | C(JJ)=1/DCMPLX(GAK,-AKACH) ! amplitude for the SW-potential |
| 999 | GOTO 30 |
| 1000 | 230 IF(LL.EQ.5.OR.LL.EQ.6.OR.LL.EQ.7)GOTO20 ! pipi |
| 1001 | IF(LL.EQ.12.OR.LL.EQ.13)GOTO20 ! piN |
| 1002 | IF(LL.EQ.8.OR.LL.EQ.9.OR.LL.EQ.18)GOTO20 ! Nd, dd |
| 1003 | IF(LL.EQ.14.OR.LL.EQ.17.OR.LL.EQ.23)GOTO27 ! K+K-, K-p, K0K0-b |
| 1004 | IF(LL.EQ.30)GOTO 28 ! pp-bar |
| 1005 | A1=RD(JJ)*FD(JJ)*AKS |
| 1006 | A2=1+.5D0*A1 |
| 1007 | IF(ICH.EQ.1)A2=A2-2*HCP*FD(JJ)/AC |
| 1008 | AKF=AK*FD(JJ) |
| 1009 | IF(ICH.EQ.1)AKF=AKF*ACH |
| 1010 | C(JJ)=FD(JJ)/DCMPLX(A2,-AKF) |
| 1011 | GOTO30 |
| 1012 | 20 CONTINUE |
| 1013 | C---Calc. scatt. ampl. C(JJ) for pipi, piN and Nd, dd |
| 1014 | JH=LL-7+2*JJ-2 |
| 1015 | IF(LL.EQ.8.OR.LL.EQ.9)GPI2=GND(AKS,JH) |
| 1016 | IF(LL.EQ.18)GPI2=GDD(AKS,JJ) |
| 1017 | IF(LL.EQ.5.OR.LL.EQ.6.OR.LL.EQ.7)GPI2=GPIPI(AKS,2) |
| 1018 | IF(LL.EQ.12.OR.LL.EQ.13)GPI2=GPIN(AKS,2) |
| 1019 | C(JJ)=1.D0/DCMPLX(GPI2,-AK) !pi+pi+, nd, pd, pi+p, dd |
| 1020 | IF(LL.NE.5.AND.LL.NE.6.AND.LL.NE.13)GOTO27 |
| 1021 | IF(LL.EQ.5.OR.LL.EQ.6)GPI1=GPIPI(AKS,1) |
| 1022 | IF(LL.EQ.13)GPI1=GPIN(AKS,1) |
| 1023 | IF(LL.EQ.5.OR.LL.EQ.13) |
| 1024 | c C(JJ)=.6667D0/DCMPLX(GPI1,-AK)+.3333D0*C(JJ) !pi+pi-,pi-p |
| 1025 | IF(LL.EQ.6)C(JJ)=.3333D0/DCMPLX(GPI1,-AK)+.6667D0*C(JJ) !pi0pi0 |
| 1026 | 27 CONTINUE |
| 1027 | C---Calc. K+K-, K0K0-b or K-p s-wave scatt. ampl. |
| 1028 | IF(LL.EQ.14.OR.LL.EQ.23)CALL CKKB |
| 1029 | c IF(LL.EQ.17)C(JJ)=DCMPLX(-3.29D0,3.55D0) ! Martin'76 (K-p) |
| 1030 | c IF(LL.EQ.17)C(JJ)=DCMPLX(3.29D0,3.55D0) ! Martin'76 (K-p) WRONG SIGN!!! |
| 1031 | IF(LL.EQ.17)C(JJ)=DCMPLX(-2.585D0,4.156D0) ! Borasoy'04 (K-p) |
| 1032 | c IF(LL.EQ.17)C(JJ)=DCMPLX(-3.371D0,3.244D0) ! Martin'81 (K-p) |
| 1033 | C---Calc. pi+pi-, pi+pi+, pd, pi+p, pi-p, K+K- or K-p s-wave scatt. ampl. |
| 1034 | C-- divided by Coulomb penetration factor squared (if ICH=1) |
| 1035 | IF(ICH.EQ.0)GOTO 30 |
| 1036 | AAK=ACH*AK ! |
| 1037 | HCP2=2*HCP/AC ! needed to calculate C(JJ) for charged particles |
| 1038 | C(JJ)=1/(1/C(JJ)-HCP2+DCMPLX(0.D0,AK-AAK)) |
| 1039 | GOTO 30 |
| 1040 | 28 CONTINUE |
| 1041 | C---Calc. pp-bar s-wave scatt. ampl. |
| 1042 | CALL CPAP |
| 1043 | 30 CONTINUE |
| 1044 | C*********************************************************************** |
| 1045 | CALL VZ(WEI) |
| 1046 | RETURN |
| 1047 | END |
| 1048 | SUBROUTINE VZ(WEI) |
| 1049 | C==> Calculates the weight WEI due to FSI |
| 1050 | IMPLICIT REAL*8 (A-H,O-Z) |
| 1051 | COMMON/FSI_JR/JRAT |
| 1052 | COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS, |
| 1053 | 1 X,Y,Z,T,RP,RPS |
| 1054 | COMMON/FSI_SPIN/RHO(10) |
| 1055 | COMMON/FSI_ETA/ETA |
| 1056 | COMMON/FSI_AA/AA |
| 1057 | COMMON/FSI_FFF/F12,F21 |
| 1058 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| 1059 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| 1060 | COMMON/FSI_FD/FD(10),RD(10) |
| 1061 | COMMON/FSI_RR/F(10) |
| 1062 | COMMON/FSI_C/C(10),AM,AMS,DM |
| 1063 | COMMON/FSI_COULPH/EIDC |
| 1064 | COMPLEX*16 F,C,G,PSI12,PSI21 |
| 1065 | COMPLEX*16 F12,F21 |
| 1066 | COMPLEX*16 EIDC |
| 1067 | COMPLEX*8 Z8,CGAMMA |
| 1068 | COMMON/FSI_FFPN/FF12,FF21 |
| 1069 | COMPLEX*16 FF12,FF21 |
| 1070 | COMMON/FSI_2CHA/AK2,AK2S,AAK,HCP2,AMU2_AMU1 ! k* (kappa) for 2-nd channel |
| 1071 | WEI=0.D0 |
| 1072 | IF(JRAT.EQ.1)GOTO 11 |
| 1073 | RHOS=AK*RP |
| 1074 | HS=X*PPX+Y*PPY+Z*PPZ |
| 1075 | IF(RHOS.LT.15.D0.AND.RHOS+DABS(HS).LT.20.D0)GOTO 2 |
| 1076 | C---Calc. EIDC=exp(i*Coul.Ph.); |
| 1077 | C-- used in calc. of hypergeom. f-s in SEQA, FAS at k*R > 15, 20 |
| 1078 | Z8=CMPLX(1.,SNGL(ETA)) |
| 1079 | Z8=CGAMMA(Z8) |
| 1080 | EIDC=Z8/CABS(Z8) |
| 1081 | 2 CALL FF(RHOS,HS) |
| 1082 | 11 MSP=MSPIN |
| 1083 | IF(ISI.EQ.0)GOTO 4 ! the strong interaction ON (OFF) if ISI=1(0) |
| 1084 | IF(RP.LT.AA)GOTO 4 |
| 1085 | IF(JRAT.NE.1) CALL FIRT |
| 1086 | IF(IQS.EQ.0)GOTO 5 ! the quantum statistics ON (OFF) if IQS=1(0) |
| 1087 | JSIGN=-1 |
| 1088 | DO 1 JJ=1,MSP |
| 1089 | JSIGN=-JSIGN |
| 1090 | G=F(JJ)*C(JJ) |
| 1091 | IF(ICH.EQ.1)G=G*ACHR |
| 1092 | PSI12=FF12*(F12+G) |
| 1093 | PSI21=FF21*(F21+G) |
| 1094 | G=PSI12+JSIGN*PSI21 |
| 1095 | 1 WEI=WEI+RHO(JJ)*(DREAL(G)**2+DIMAG(G)**2) |
| 1096 | GOTO 8 |
| 1097 | 5 DO 6 JJ=1,MSP |
| 1098 | G=F(JJ)*C(JJ) |
| 1099 | IF(ICH.EQ.1)G=G*ACHR |
| 1100 | CW WRITE(6,38)'JJ ',JJ,'F ',DREAL(F(JJ)),DIMAG(F(JJ)) |
| 1101 | CW WRITE(6,38)'JJ ',JJ,'C ',DREAL(C(JJ)),DIMAG(C(JJ)) |
| 1102 | CW WRITE(6,38)'JJ ',JJ,'G ',DREAL(G),DIMAG(G) |
| 1103 | CW WRITE(6,38)'JJ ',JJ,'F12+G ',DREAL(F12+G),DIMAG(F12+G) |
| 1104 | CW WRITE(6,38)'JJ ',JJ,'F21+G ',DREAL(F21+G),DIMAG(F21+G) |
| 1105 | 38 FORMAT(A7,I3,A7,2E11.4) |
| 1106 | PSI12=FF12*(F12+G) |
| 1107 | 6 WEI=WEI+RHO(JJ)*(DREAL(PSI12)**2+DIMAG(PSI12)**2) |
| 1108 | c--- Account for nn-bar->pp-bar channel --------------------------- |
| 1109 | IF(LL.EQ.30)THEN |
| 1110 | DO 61 JJ=1,MSP |
| 1111 | HH=RHO(JJ)*(DREAL(C(JJ+2))**2+DIMAG(C(JJ+2))**2)* |
| 1112 | * AMU2_AMU1*ACH/RPS |
| 1113 | IF(AK2S.LT.0)HH=HH*DEXP(-2*RP*AK2) |
| 1114 | 61 WEI=WEI+HH |
| 1115 | ENDIF |
| 1116 | c------------------------------------------------------------------ |
| 1117 | RETURN |
| 1118 | 4 PSI12=FF12*F12 |
| 1119 | IF(IQS.EQ.0)GOTO 50 ! the quantum statistics ON (OFF) if IQS=1(0) |
| 1120 | PSI21=FF21*F21 |
| 1121 | JSIGN=-1 |
| 1122 | DO 3 JJ=1,MSP |
| 1123 | JSIGN=-JSIGN |
| 1124 | G=PSI12+JSIGN*PSI21 |
| 1125 | 3 WEI=WEI+RHO(JJ)*(DREAL(G)**2+DIMAG(G)**2) |
| 1126 | GOTO 8 |
| 1127 | 50 WEI=DREAL(PSI12)**2+DIMAG(PSI12)**2 |
| 1128 | RETURN |
| 1129 | 8 WEI=WEI/2 |
| 1130 | RETURN |
| 1131 | END |
| 1132 | SUBROUTINE FIRT |
| 1133 | C---CALC. THE F(JJ) |
| 1134 | C-- F(JJ)*C(JJ)= DEVIATION OF THE BETHE-SALPETER AMPL. FROM PLANE WAVE |
| 1135 | IMPLICIT REAL*8 (A-H,O-Z) |
| 1136 | COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS, |
| 1137 | 1 X,Y,Z,T,RP,RPS |
| 1138 | COMMON/FSI_SHH/SH,CHH |
| 1139 | COMMON/FSI_BP/B,P |
| 1140 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| 1141 | COMMON/FSI_C/C(10),AM,AMS,DM |
| 1142 | COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10), |
| 1143 | 1 SBKRB(10),SDKK(10) |
| 1144 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| 1145 | COMMON/FSI_RR/F(10) |
| 1146 | EQUIVALENCE(RSS,RP),(TSS,T) |
| 1147 | COMPLEX*16 F,C,CH1 |
| 1148 | MSP=MSPIN |
| 1149 | DO 10 JJ=1,MSP |
| 1150 | IF(JJ.GT.1)GOTO 3 |
| 1151 | XRA=2*RSS/AC |
| 1152 | IF(AK.NE.0.D0)GOTO2 |
| 1153 | SHK=1.D0 |
| 1154 | SH=.0D0 |
| 1155 | SHH=SH |
| 1156 | CHH=1/RSS |
| 1157 | GOTO3 |
| 1158 | 2 H=AK*RSS |
| 1159 | CALL GST(XRA,H) |
| 1160 | SH=SH/RSS |
| 1161 | CHH=CHH/RSS |
| 1162 | SHH=SH |
| 1163 | IF(ICH.EQ.1) SHH=ACH*SH |
| 1164 | 3 IF(NS.EQ.2)GOTO1 |
| 1165 | C---F= ASYMPTOTIC FORMULA (T= 0 APPROX.); NS= 4 |
| 1166 | 6 F(JJ)=DCMPLX(CHH,SHH) |
| 1167 | IF(NS.NE.1)GOTO 10 |
| 1168 | C---F INSIDE THE SQUARE-WELL (T= 0 APPROX.); NS= 1 |
| 1169 | IF(RSS.GE.RB(JJ)) GOTO 10 |
| 1170 | IF(AK.NE.0.D0.AND.JJ.EQ.1)SHK=SH/AK |
| 1171 | H=BK(JJ)*RSS |
| 1172 | CALL GST(XRA,H) |
| 1173 | SKR=B*BK(JJ) |
| 1174 | F(JJ)=DCMPLX(CDK(JJ),SDK(JJ))*SKR |
| 1175 | CH1=(SDKK(JJ)*SKR-SHK*SBKRB(JJ))/C(JJ) |
| 1176 | F(JJ)=(F(JJ)+CH1)/SBKRB(JJ) |
| 1177 | GOTO 10 |
| 1178 | 1 CONTINUE |
| 1179 | C---F= ASYMPTOTIC FORMULA (T= 0 NOT REQUIRED); NS= 2 |
| 1180 | IF(JJ.GT.1)GOTO 8 |
| 1181 | IF(TSS.EQ.0.D0)GOTO6 |
| 1182 | TSSA=DABS(TSS) |
| 1183 | IF(DM.NE.0.D0)GOTO 11 |
| 1184 | H=AM*.5D0/TSSA |
| 1185 | IF(AK.NE.0.D0)GOTO4 |
| 1186 | HM=H*RPS |
| 1187 | IF(HM.GE.3.D15)GOTO6 |
| 1188 | FS1=DFRSIN(HM) |
| 1189 | FC1=DFRCOS(HM) |
| 1190 | FC2=FC1 |
| 1191 | FS2=FS1 |
| 1192 | GOTO5 |
| 1193 | 4 CONTINUE |
| 1194 | H1=AK*TSSA/AM |
| 1195 | HM=H*(RSS-H1)**2 |
| 1196 | HP=H*(RSS+H1)**2 |
| 1197 | IF(HP.GE.3.D15)GOTO6 |
| 1198 | FS1=DFRSIN(HM) |
| 1199 | FC1=DFRCOS(HM) |
| 1200 | FS2=DFRSIN(HP) |
| 1201 | FC2=DFRCOS(HP) |
| 1202 | GOTO 5 |
| 1203 | 11 CONTINUE |
| 1204 | FS1=0.D0 |
| 1205 | FS2=0.D0 |
| 1206 | FC1=0.D0 |
| 1207 | FC2=0.D0 |
| 1208 | DO 13 I=1,2 |
| 1209 | IF(I.EQ.1)TSSH=TSSA*(1+DM) |
| 1210 | IF(I.EQ.2)TSSH=TSSA*(1-DM) |
| 1211 | H=AM*.5D0/TSSH |
| 1212 | IF(AK.NE.0.D0)GOTO 12 |
| 1213 | HM=H*RPS |
| 1214 | IF(HM.GE.3.D15)GOTO6 |
| 1215 | FS1=FS1+DFRSIN(HM)/2 |
| 1216 | FC1=FC1+DFRCOS(HM)/2 |
| 1217 | IF(I.EQ.1)GOTO 13 |
| 1218 | FC2=FC1 |
| 1219 | FS2=FS1 |
| 1220 | GOTO 13 |
| 1221 | 12 CONTINUE |
| 1222 | H1=AK*TSSH/AM |
| 1223 | HM=H*(RSS-H1)**2 |
| 1224 | HP=H*(RSS+H1)**2 |
| 1225 | IF(HP.GE.3.D15)GOTO6 |
| 1226 | FS1=FS1+DFRSIN(HM)/2 |
| 1227 | FC1=FC1+DFRCOS(HM)/2 |
| 1228 | FS2=FS2+DFRSIN(HP)/2 |
| 1229 | FC2=FC2+DFRCOS(HP)/2 |
| 1230 | 13 CONTINUE |
| 1231 | 5 C12=FC1+FS2 |
| 1232 | S12=FS1+FC2 |
| 1233 | A12=FS1-FC2 |
| 1234 | A21=FS2-FC1 |
| 1235 | A2=.5D0*(CHH*(A12+A21)+SH*(A12-A21))+SHH |
| 1236 | A1=.5D0*(CHH*(C12+S12)+SH*(C12-S12)) |
| 1237 | F(JJ)=.3989422D0*DCMPLX(A1,A2) |
| 1238 | GOTO 10 |
| 1239 | 8 F(JJ)=F(1) |
| 1240 | 10 CONTINUE |
| 1241 | RETURN |
| 1242 | END |
| 1243 | FUNCTION EXF(X) |
| 1244 | IMPLICIT REAL*8 (A-H,O-Z) |
| 1245 | IF(X.LT.-15.D0) GO TO 1 |
| 1246 | EXF=DEXP(X) |
| 1247 | RETURN |
| 1248 | 1 EXF=.0D0 |
| 1249 | RETURN |
| 1250 | END |
| 1251 | SUBROUTINE SEQ(X,H) |
| 1252 | C---CALC. FUNCTIONS B, P (EQS. (17) OF G-K-L-L); |
| 1253 | C-- NEEDED TO CALC. THE CONFLUENT HYPERGEOMETRIC FUNCTION GST. |
| 1254 | IMPLICIT REAL*8 (A-H,O-Z) |
| 1255 | COMMON/FSI_BP/B,P |
| 1256 | DIMENSION BH(3),PH(3) |
| 1257 | DATA ERR/1.D-7/ |
| 1258 | BH(1)=1.D0 |
| 1259 | PH(1)=1.D0 |
| 1260 | PH(2)=.0D0 |
| 1261 | BH(2)=.5D0*X |
| 1262 | B=1+BH(2) |
| 1263 | P=1.D0 |
| 1264 | HS=H*H |
| 1265 | J=0 |
| 1266 | 2 J=J+1 |
| 1267 | BH(3)=(X*BH(2)-HS*BH(1))/((J+1)*(J+2)) |
| 1268 | PH(3)=(X*PH(2)-HS*PH(1)-(2*J+1)*X*BH(2))/(J*(J+1)) |
| 1269 | B=B+BH(3) |
| 1270 | P=P+PH(3) |
| 1271 | Z=DABS(BH(2))+DABS(BH(3))+DABS(PH(2))+DABS(PH(3)) |
| 1272 | IF(Z.LT.ERR)RETURN |
| 1273 | BH(1)=BH(2) |
| 1274 | BH(2)=BH(3) |
| 1275 | PH(1)=PH(2) |
| 1276 | PH(2)=PH(3) |
| 1277 | GOTO 2 |
| 1278 | END |
| 1279 | SUBROUTINE SEQA(X,H) |
| 1280 | C---CALC. FUNCTIONS CHH=REAL(GST), SH=IMAG(GST)/ACH, B=SH/H |
| 1281 | C-- IN THE ASYMPTOTIC REGION H=K*R >> 1. |
| 1282 | IMPLICIT REAL*8 (A-H,O-Z) |
| 1283 | COMMON/FSI_BP/B,P |
| 1284 | COMMON/FSI_SHH/SH,CHH |
| 1285 | COMMON/FSI_ETA/ETA |
| 1286 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| 1287 | COMMON/FSI_COULPH/EIDC |
| 1288 | COMPLEX*16 EIDC,GST |
| 1289 | ARG=H-ETA*DLOG(2*H) |
| 1290 | GST=DCMPLX(DCOS(ARG),DSIN(ARG)) |
| 1291 | GST=ACHR*EIDC*GST |
| 1292 | CHH=DREAL(GST) |
| 1293 | SH=DIMAG(GST)/ACH |
| 1294 | B=SH/H |
| 1295 | RETURN |
| 1296 | END |
| 1297 | SUBROUTINE FF(RHO,H) |
| 1298 | C---Calc. F12, F21; |
| 1299 | C-- F12= FF0* plane wave, FF0=F*ACHR, |
| 1300 | C---F is the confluent hypergeometric function, |
| 1301 | C-- ACHR=sqrt(ACH), where ACH is the Coulomb factor |
| 1302 | IMPLICIT REAL*8 (A-H,O-Z) |
| 1303 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| 1304 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| 1305 | COMMON/FSI_ETA/ETA |
| 1306 | COMMON/FSI_FFF/F12,F21 |
| 1307 | COMPLEX*16 FF0,F12,F21 |
| 1308 | C=DCOS(H) |
| 1309 | S=DSIN(H) |
| 1310 | F12=DCMPLX(C,-S) |
| 1311 | F21=DCMPLX(C,S) |
| 1312 | IF(ICH.EQ.0)RETURN |
| 1313 | RHOP=RHO+H |
| 1314 | RHOM=RHO-H |
| 1315 | F12=FF0(RHO,H)*F12 |
| 1316 | F21=FF0(RHO,-H)*F21 |
| 1317 | RETURN |
| 1318 | END |
| 1319 | FUNCTION FAS(RKS) |
| 1320 | C-- FAS=F*ACHR |
| 1321 | C---F is the confluent hypergeometric function at k*r >> 1 |
| 1322 | C-- ACHR=sqrt(ACH), where ACH is the Coulomb factor |
| 1323 | IMPLICIT REAL*8 (A-H,O-Z) |
| 1324 | COMPLEX*16 FAS,EIDC,ZZ1 |
| 1325 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| 1326 | COMMON/FSI_ETA/ETA |
| 1327 | COMMON/FSI_COULPH/EIDC |
| 1328 | D1=DLOG(RKS)*ETA |
| 1329 | D2=ETA*ETA/RKS |
| 1330 | ZZ1=DCMPLX(DCOS(D1),DSIN(D1))/EIDC |
| 1331 | FAS=DCMPLX(1.D0,-D2)*ZZ1 |
| 1332 | FAS=FAS-DCMPLX(DCOS(RKS),DSIN(RKS))*ETA/RKS/ZZ1 |
| 1333 | RETURN |
| 1334 | END |
| 1335 | FUNCTION FF0(RHO,H) |
| 1336 | C-- FF0=F*ACHR |
| 1337 | C-- F is the confluent hypergeometric function |
| 1338 | C-- (Eq. (15) of G-K-L-L), F= 1 at r* << AC |
| 1339 | C-- ACHR=sqrt(ACH), where ACH is the Coulomb factor |
| 1340 | IMPLICIT REAL*8 (A-H,O-Z) |
| 1341 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| 1342 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| 1343 | COMMON/FSI_ETA/ETA |
| 1344 | COMPLEX*16 ALF,ALF1,Z,S,A,FF0,FAS |
| 1345 | DATA ERR/1.D-5/ |
| 1346 | S=DCMPLX(1.D0,0.D0) |
| 1347 | FF0=S |
| 1348 | RHOP=RHO+H |
| 1349 | CC GOTO 5 ! rejects the approx. calcul. of hyperg. f-ion F |
| 1350 | IF(RHOP.LT.20.D0)GOTO5 |
| 1351 | FF0=FAS(RHOP) ! approx. calc. |
| 1352 | RETURN |
| 1353 | 5 ALF=DCMPLX(.0D0,-ETA) |
| 1354 | ALF1=ALF-1 |
| 1355 | Z=DCMPLX(.0D0,RHOP) |
| 1356 | J=0 |
| 1357 | 3 J=J+1 |
| 1358 | A=(ALF1+J)/(J*J) |
| 1359 | S=S*A*Z |
| 1360 | FF0=FF0+S |
| 1361 | ZR=DABS(DREAL(S)) |
| 1362 | ZI=DABS(DIMAG(S)) |
| 1363 | IF((ZR+ZI).GT.ERR)GOTO3 |
| 1364 | FF0=FF0*ACHR |
| 1365 | RETURN |
| 1366 | END |
| 1367 | FUNCTION HC(XA) |
| 1368 | C---HC = h-function of Landau-Lifshitz: h(x)=Re[psi(1-i/x)]+ln(x) |
| 1369 | C-- psi(x) is the digamma function (the logarithmic derivative of |
| 1370 | C-- the gamma function) |
| 1371 | IMPLICIT REAL*8 (A-H,O-Z) |
| 1372 | DIMENSION BN(15) |
| 1373 | DATA BN/.8333333333D-1,.8333333333D-2,.396825396825D-2, |
| 1374 | 1 .4166666667D-2,.7575757576D-2,.2109279609D-1, |
| 1375 | 2 .8333333333D-1,.4432598039D0 ,.305395433D1, |
| 1376 | 3 .2645621212D2, .2814601449D3, .3607510546D4, |
| 1377 | 4 .5482758333D5, .9749368235D6, .200526958D8/ |
| 1378 | X=DABS(XA) |
| 1379 | IF(X.LT..33D0) GOTO 1 |
| 1380 | CC IF(X.GE.3.5D0) GO TO 2 |
| 1381 | S=.0D0 |
| 1382 | N=0 |
| 1383 | 3 N=N+1 |
| 1384 | DS=1.D0/N/((N*X)**2+1) |
| 1385 | S=S+DS |
| 1386 | IF(DS.GT.0.1D-12) GOTO 3 |
| 1387 | C---Provides 7 digit accuracy |
| 1388 | HC=S-.5772156649D0+DLOG(X) |
| 1389 | RETURN |
| 1390 | CC 2 HC=1.2D0/X**2+DLOG(X)-.5772156649 D0 |
| 1391 | CC RETURN |
| 1392 | 1 X2=X*X |
| 1393 | XP=X2 |
| 1394 | HC=0.D0 |
| 1395 | IMA=9 |
| 1396 | IF(X.LT.0.1D0)IMA=3 |
| 1397 | DO 4 I=1,IMA |
| 1398 | HC=HC+XP*BN(I) |
| 1399 | 4 XP=X2*XP |
| 1400 | RETURN |
| 1401 | END |
| 1402 | FUNCTION ACP(X) |
| 1403 | C--- ACP = COULOMB PENETRATION FACTOR |
| 1404 | IMPLICIT REAL*8 (A-H,O-Z) |
| 1405 | IF(X.LT.0.05D0.AND.X.GE.0.D0) GO TO 1 |
| 1406 | Y=6.2831853D0/X |
| 1407 | ACP=Y/(EXF(Y)-1) |
| 1408 | RETURN |
| 1409 | 1 ACP=1.D-6 |
| 1410 | RETURN |
| 1411 | END |
| 1412 | SUBROUTINE GST(X,H) |
| 1413 | C---CALC. THE CONFL. HYPERGEOM. F-N = CHH+i*SH |
| 1414 | C-- AND THE COULOMB F-S B, P (CALLS SEQ OR SEQA). |
| 1415 | IMPLICIT REAL*8 (A-H,O-Z) |
| 1416 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| 1417 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| 1418 | COMMON/FSI_SHH/SH,CHH |
| 1419 | COMMON/FSI_BP/B,P |
| 1420 | 1 IF(ICH.EQ.1)GOTO 2 |
| 1421 | 3 SH=DSIN(H) |
| 1422 | CHH=DCOS(H) |
| 1423 | B=1.D0 |
| 1424 | IF(H.NE.0.D0)B=SH/H |
| 1425 | P=CHH |
| 1426 | RETURN |
| 1427 | 2 CONTINUE |
| 1428 | IF(H.GT.15.D0)GOTO4 ! comment out if you want to reject |
| 1429 | ! the approximate calculation of hyperg. f-ion G |
| 1430 | CALL SEQ(X,H) ! exact calculation |
| 1431 | SH=H*B |
| 1432 | CHH=P+B*X*(DLOG(DABS(X))+HPR) |
| 1433 | RETURN |
| 1434 | 4 CALL SEQA(X,H) |
| 1435 | RETURN |
| 1436 | END |
| 1437 | FUNCTION FF1(RHO,H) |
| 1438 | C---FF1=FF0; used for particle-nucleus system |
| 1439 | C-- FF0=F12*ACHR |
| 1440 | C-- F12 is the confluent hypergeometric function |
| 1441 | C-- (Eq. (15) of G-K-L-L), F12= 1 at r* << AC |
| 1442 | C-- ACHR=sqrt(ACH), where ACH is the Coulomb factor |
| 1443 | IMPLICIT REAL*8 (A-H,O-Z) |
| 1444 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN |
| 1445 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| 1446 | COMMON/FSI_ETA/ETA |
| 1447 | COMMON/FSI_COULPH/EIDC |
| 1448 | COMMON/FSI_ICH1/ICH1 |
| 1449 | COMPLEX*16 FF0,FF1 |
| 1450 | COMPLEX*16 EIDC |
| 1451 | COMPLEX*8 Z8,CGAMMA |
| 1452 | FF1=DCMPLX(1.D0,0.D0) |
| 1453 | IF(ICH1.EQ.0)GOTO 2 |
| 1454 | IF(RHO.LT.15.D0.AND.RHO+H.LT.20.D0)GOTO 2 |
| 1455 | C---Calc. EIDC=exp(i*Coul.Ph.); |
| 1456 | C-- used in calc. of hypergeom. f-s in SEQA, FAS at k*R > 15, 20 |
| 1457 | Z8=CMPLX(1.,SNGL(ETA)) |
| 1458 | Z8=CGAMMA(Z8) |
| 1459 | EIDC=Z8/CABS(Z8) |
| 1460 | 2 FF1=FF0(RHO,H) |
| 1461 | RETURN |
| 1462 | END |
| 1463 | |
| 1464 | FUNCTION G(AK) |
| 1465 | C---Used to calculate SW scattering amplitude for alpa-alpha system |
| 1466 | C-- and for sear.f (square well potential search) |
| 1467 | C---NOTE THAT SCATT. AMPL.= 1/CMPLX(G(AK),-AK*ACH) |
| 1468 | IMPLICIT REAL*8 (A-H,O-Z) |
| 1469 | COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10), |
| 1470 | 1 SBKRB(10),SDKK(10) |
| 1471 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| 1472 | COMMON/FSI_ACH/HPR,AC,ACH,ACHR,JJ,MSPIN |
| 1473 | COMMON/FSI_BP/B,P/FSI_DERIV/BPR,PPR/FSI_SHH/SH,CHH |
| 1474 | COMMON/FSI_DAK/DAK,HAC,IFUN |
| 1475 | HAC=.0D0 |
| 1476 | ACH=1.D0 |
| 1477 | IF(ICH.EQ.0)GOTO 1 |
| 1478 | XH=AC*AK |
| 1479 | HCP=HC(XH) |
| 1480 | HPR=HCP+.1544313298D0 |
| 1481 | ACH=ACP(XH) |
| 1482 | HAC=2*HCP/AC |
| 1483 | 1 AKS=AK**2 |
| 1484 | BK(JJ)=DSQRT(AKS+EB(JJ)**2) ! kappa=kp |
| 1485 | X=2*RB(JJ)/AC |
| 1486 | H=BK(JJ)*RB(JJ) ! kp*d |
| 1487 | CALL GST(X,H) |
| 1488 | BRHO=B ! B(kp,d) |
| 1489 | SBKRB(JJ)=SH ! kp*d*B(kp,d) |
| 1490 | CALL DERIW(X,H) |
| 1491 | BRHOP=BPR ! B'(kp,d)= dB(kp,r)/dln(r) at r=d |
| 1492 | H=AK*RB(JJ) |
| 1493 | CALL GST(X,H) |
| 1494 | CDK(JJ)=CHH ! ReG(k,d) |
| 1495 | BRHOS=B ! B(k,d) |
| 1496 | PRHOS=P ! P(k,d) |
| 1497 | SDK(JJ)=SH |
| 1498 | IF(ICH.EQ.0)GOTO 2 |
| 1499 | SDK(JJ)=ACH*SH ! ImG(k,d) |
| 1500 | IF(AK.EQ.0.D0.AND.AC.LT.0.D0)SDK(JJ)=3.14159*X*B |
| 1501 | 2 SDKK(JJ)=RB(JJ) |
| 1502 | IF(AK.NE.0.D0)SDKK(JJ)=SH/AK ! d*B(k,d) |
| 1503 | CALL DERIW(X,H) ! PPR=P'(k,d)= dP(k,r)/dln(r) at r=d |
| 1504 | ZZ=PPR-PRHOS |
| 1505 | IF(ICH.EQ.1)ZZ=ZZ+X*(BRHOS+BPR*(DLOG(DABS(X))+HPR)) |
| 1506 | C ZZ= P'(k,d)-P(k,d)+x*{B(k,d)+B'(k,d)*[ln!x!+2*C-1+h(k*ac)]} |
| 1507 | GG=(BRHOP*CDK(JJ)-BRHO*ZZ)/RB(JJ) |
| 1508 | C GG= [B'(kp,d)*ReG(k,d)-B(kp,d)*ZZ]/d |
| 1509 | G=GG/(BRHO*BPR-BRHOP*BRHOS) |
| 1510 | C G= GG/[B(kp,d)*B'(k,d)-B'(kp,d)*B(k,d)] |
| 1511 | RETURN |
| 1512 | END |
| 1513 | SUBROUTINE DERIW(X,H) |
| 1514 | C---CALLED BY F-N G(AK) |
| 1515 | IMPLICIT REAL*8 (A-H,O-Z) |
| 1516 | COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S |
| 1517 | COMMON/FSI_BP/B,P/FSI_DERIV/BPR,PPR |
| 1518 | HH=.1D-3 |
| 1519 | CALL GST(X,H-HH) |
| 1520 | Q1=P |
| 1521 | B1=B |
| 1522 | CALL GST(X,H+HH) |
| 1523 | HHH=HH+HH |
| 1524 | BPR=H*(B-B1)/HHH |
| 1525 | PPR=H*(P-Q1)/HHH |
| 1526 | IF(ICH.EQ.0)RETURN |
| 1527 | CALL GST(X-HH,H) |
| 1528 | Q1=P |
| 1529 | B1=B |
| 1530 | CALL GST(X+HH,H) |
| 1531 | BPR=BPR+X*(B-B1)/HHH |
| 1532 | PPR=PPR+X*(P-Q1)/HHH |
| 1533 | RETURN |
| 1534 | END |
| 1535 | C================================================================ |
| 1536 | |
| 1537 | SUBROUTINE READ_FILE(KEY,CH8,INT4,REAL8,IERR,NUNIT) |
| 1538 | C ========== |
| 1539 | C |
| 1540 | C Routine to read one parameter of the program in the file |
| 1541 | C DATA NUNIT defined in FSI3B EXEC |
| 1542 | C NUNIT=11 for IBM in Nantes, 4 for SUN in Prague |
| 1543 | C |
| 1544 | C INPUT : KEY (CHARACTER*10) : |
| 1545 | C OUTPUT : case of KEY : CH8 : (CHARACTER*8) |
| 1546 | C INT4 : (INTEGER*4) |
| 1547 | C REAL8 : (REAL*8) |
| 1548 | C (only one of them) |
| 1549 | C IERR (INTEGER) : 0 : no error |
| 1550 | C 1 : key not found |
| 1551 | |
| 1552 | CHARACTER*10 KEY,TEST |
| 1553 | CHARACTER*4 TYPE |
| 1554 | CHARACTER*8 CH8 |
| 1555 | INTEGER*4 INT4 |
| 1556 | REAL*8 REAL8 |
| 1557 | INTEGER*4 IERR |
| 1558 | |
| 1559 | IERR=0 |
| 1560 | REWIND(NUNIT) |
| 1561 | 1 READ(NUNIT,FMT='(A10,2X,A4)')TEST,TYPE |
| 1562 | IF (TEST.EQ.KEY) THEN |
| 1563 | BACKSPACE(NUNIT) |
| 1564 | IF (TYPE.EQ.'CHAR') READ(NUNIT,FMT='(18X,A8,54X)')CH8 |
| 1565 | IF (TYPE.EQ.'INT4') READ(NUNIT,FMT='(18X,I8,54X)')INT4 |
| 1566 | IF (TYPE.EQ.'REA8') READ(NUNIT,FMT='(18X,F10.5,52X)')REAL8 |
| 1567 | ELSE |
| 1568 | IF (TEST.NE.'* E.O.F. *') THEN |
| 1569 | GOTO 1 |
| 1570 | ELSE |
| 1571 | IERR=1 |
| 1572 | ENDIF |
| 1573 | ENDIF |
| 1574 | c IF(IERR.EQ.1)STOP |
| 1575 | RETURN |
| 1576 | END |
| 1577 | |
git://git.uio.no / u / mrichter / AliRoot.git / blame