]>
Commit | Line | Data |
---|---|---|
f5ab1a71 | 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 | |
f5ab1a71 | 125 | COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS |
ff4431bb | 126 | COMPLEX*16 FF12,FF21 |
f5ab1a71 | 127 | C------------------------------------------------------------------ |
f5ab1a71 | 128 | C==> AC1,2 = "relativistic" Bohr radii for particle-nucleus systems |
f5ab1a71 | 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----------------------------------------------------------- | |
ff4431bb | 135 | J=1 |
136 | CALL FSIPN !weight due to particle-nucleus FSI | |
f5ab1a71 | 137 | JRAT=0 |
ff4431bb | 138 | CALL FSIWF !weight due to particle-particle-nucleus FSI |
f5ab1a71 | 139 | WEIN=WEI |
ff4431bb | 140 | IF(I3C*J.NE.0) THEN |
f5ab1a71 | 141 | FF12=DCMPLX(1.D0,0.D0) |
142 | FF21=DCMPLX(1.D0,0.D0) | |
143 | JRAT=1 | |
ff4431bb | 144 | CALL VZ ! weight due to particle-particle FSI |
145 | ENDIF | |
f5ab1a71 | 146 | RETURN |
147 | END | |
ff4431bb | 148 | |
f5ab1a71 | 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 | ||
f5ab1a71 | 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 | ||
ff4431bb | 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 | ||
f5ab1a71 | 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 | |
ff4431bb | 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 | |
f5ab1a71 | 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 | |
ff4431bb | 476 | |
f5ab1a71 | 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 | |
f5ab1a71 | 483 | WEI=0.D0 |
484 | IF(JRAT.EQ.1)GOTO 11 | |
485 | RHOS=AK*RP | |
486 | HS=X*PPX+Y*PPY+Z*PPZ | |
f5ab1a71 | 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 | |
f5ab1a71 | 490 | Z8=CMPLX(1.,SNGL(ETA)) |
ff4431bb | 491 | Z8=CGAMMA(Z8) |
f5ab1a71 | 492 | EIDC=Z8/CABS(Z8) |
ff4431bb | 493 | |
494 | c write(*,*)'Z8 EIDC',Z8,EIDC | |
f5ab1a71 | 495 | |
496 | 2 CALL FF(RHOS,HS) | |
497 | 11 MSP=MSPIN | |
f5ab1a71 | 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 | |
ff4431bb | 510 | 1 WEI=WEI+RHO(JJ)*(DREAL(G)**2+DIMAG(G)**2) |
f5ab1a71 | 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 | |
f5ab1a71 | 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 | |
f5ab1a71 | 531 | 3 WEI=WEI+RHO(JJ)*(DREAL(G)**2+DIMAG(G)**2) |
f5ab1a71 | 532 | GOTO 8 |
f5ab1a71 | 533 | 50 WEI=DREAL(PSI12)**2+DIMAG(PSI12)**2 |
f5ab1a71 | 534 | RETURN |
535 | 8 WEI=WEI/2 | |
ff4431bb | 536 | |
537 | c [M %/WRITE(*,*)WEI | |
538 | ||
f5ab1a71 | 539 | RETURN |
540 | END | |
ff4431bb | 541 | |
542 | ||
f5ab1a71 | 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 | |
f5ab1a71 | 599 | FS1=DFRSIN(HM) |
600 | FC1=DFRCOS(HM) | |
f5ab1a71 | 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 | |
f5ab1a71 | 609 | FS1=DFRSIN(HM) |
610 | FC1=DFRCOS(HM) | |
611 | FS2=DFRSIN(HP) | |
612 | FC2=DFRCOS(HP) | |
f5ab1a71 | 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 | |
f5ab1a71 | 626 | FS1=FS1+DFRSIN(HM)/2 |
627 | FC1=FC1+DFRCOS(HM)/2 | |
f5ab1a71 | 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 | |
f5ab1a71 | 637 | FS1=FS1+DFRSIN(HM)/2 |
638 | FC1=FC1+DFRCOS(HM)/2 | |
639 | FS2=FS2+DFRSIN(HP)/2 | |
640 | FC2=FC2+DFRCOS(HP)/2 | |
f5ab1a71 | 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 | |
f5ab1a71 | 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 | |
f5ab1a71 | 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 | |
f5ab1a71 | 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 | |
f5ab1a71 | 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 | |
f5ab1a71 | 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 | |
f5ab1a71 | 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 | |
f5ab1a71 | 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 | |
f5ab1a71 | 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 | |
f5ab1a71 | 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 | |
ff4431bb | 848 | |
f5ab1a71 | 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 | |
f5ab1a71 | 869 | Z8=CMPLX(1.,SNGL(ETA)) |
870 | Z8=CGAMMA(Z8) | |
f5ab1a71 | 871 | EIDC=Z8/CABS(Z8) |
f5ab1a71 | 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 | |
f5ab1a71 | 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 | |
f5ab1a71 | 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 | |
f5ab1a71 | 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 | |
ff4431bb | 944 | C================================================================ |
945 | ||
f5ab1a71 | 946 | |
947 | c-------------#include "gen/pilot.h" | |
948 | FUNCTION DFRSIN(X) | |
ff4431bb | 949 | |
950 | IMPLICIT DOUBLE PRECISION (A-H,O-Z) | |
f5ab1a71 | 951 | |
952 | DIMENSION A(0:16),B(0:15),C1(0:25),C2(0:28) | |
ff4431bb | 953 | |
f5ab1a71 | 954 | PARAMETER (Z1 = 1, R8 = Z1/8, R32 = Z1/32) |
ff4431bb | 955 | |
f5ab1a71 | 956 | DATA C0 /1.25331 41373 15500 3D0/ |
ff4431bb | 957 | |
f5ab1a71 | 958 | DATA NA,NB,NC1,NC2 /16,15,25,28/ |
ff4431bb | 959 | |
f5ab1a71 | 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" | |
ff4431bb | 1124 | FUNCTION CGAMMA(Z) |
1125 | ||
f5ab1a71 | 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') | |
ff4431bb | 1132 | |
f5ab1a71 | 1133 | DIMENSION C(0:15) |
ff4431bb | 1134 | |
f5ab1a71 | 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 |