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