]> 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
Extra Control-Ms removed
[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)
148aadde 799 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
f5ab1a71 800
801 DIMENSION A(0:16),B(0:15),C1(0:25),C2(0:28)
f5ab1a71 802 PARAMETER (Z1 = 1, R8 = Z1/8, R32 = Z1/32)
f5ab1a71 803 DATA C0 /1.25331 41373 15500 3D0/
f5ab1a71 804 DATA NA,NB,NC1,NC2 /16,15,25,28/
f5ab1a71 805 DATA A( 0) / 0.76435 13866 41860 002D0/
806 DATA A( 1) /-0.43135 54754 76601 793D0/
807 DATA A( 2) / 0.43288 19997 97266 531D0/
808 DATA A( 3) /-0.26973 31033 83871 110D0/
809 DATA A( 4) / 0.08416 04532 08769 354D0/
810 DATA A( 5) /-0.01546 52448 44613 820D0/
811 DATA A( 6) / 0.00187 85542 34398 220D0/
812 DATA A( 7) /-0.00016 26497 76188 875D0/
813 DATA A( 8) / 0.00001 05739 76563 833D0/
814 DATA A( 9) /-0.00000 05360 93398 892D0/
815 DATA A(10) / 0.00000 00218 16584 549D0/
816 DATA A(11) /-0.00000 00007 29016 212D0/
817 DATA A(12) / 0.00000 00000 20373 325D0/
818 DATA A(13) /-0.00000 00000 00483 440D0/
819 DATA A(14) / 0.00000 00000 00009 865D0/
820 DATA A(15) /-0.00000 00000 00000 175D0/
821 DATA A(16) / 0.00000 00000 00000 003D0/
822
823 DATA B( 0) / 0.63041 40431 45705 392D0/
824 DATA B( 1) /-0.42344 51140 57053 335D0/
825 DATA B( 2) / 0.37617 17264 33436 566D0/
826 DATA B( 3) /-0.16249 48915 45095 674D0/
827 DATA B( 4) / 0.03822 25577 86330 087D0/
828 DATA B( 5) /-0.00564 56347 71321 909D0/
829 DATA B( 6) / 0.00057 45495 19768 974D0/
830 DATA B( 7) /-0.00004 28707 15321 020D0/
831 DATA B( 8) / 0.00000 24512 07499 233D0/
832 DATA B( 9) /-0.00000 01109 88418 409D0/
833 DATA B(10) / 0.00000 00040 82497 317D0/
834 DATA B(11) /-0.00000 00001 24498 302D0/
835 DATA B(12) / 0.00000 00000 03200 484D0/
836 DATA B(13) /-0.00000 00000 00070 324D0/
837 DATA B(14) / 0.00000 00000 00001 336D0/
838 DATA B(15) /-0.00000 00000 00000 022D0/
839
840 DATA C1( 0) / 0.99056 04793 73497 549D0/
841 DATA C1( 1) /-0.01218 35098 31478 997D0/
842 DATA C1( 2) /-0.00248 27428 23113 060D0/
843 DATA C1( 3) / 0.00026 60949 52647 247D0/
844 DATA C1( 4) /-0.00000 10790 68987 406D0/
845 DATA C1( 5) /-0.00000 48836 81753 933D0/
846 DATA C1( 6) / 0.00000 09990 55266 368D0/
847 DATA C1( 7) /-0.00000 00750 92717 372D0/
848 DATA C1( 8) /-0.00000 00190 79487 573D0/
849 DATA C1( 9) / 0.00000 00090 90797 293D0/
850 DATA C1(10) /-0.00000 00019 66236 033D0/
851 DATA C1(11) / 0.00000 00001 64772 911D0/
852 DATA C1(12) / 0.00000 00000 63079 714D0/
853 DATA C1(13) /-0.00000 00000 36432 219D0/
854 DATA C1(14) / 0.00000 00000 10536 930D0/
855 DATA C1(15) /-0.00000 00000 01716 438D0/
856 DATA C1(16) /-0.00000 00000 00107 124D0/
857 DATA C1(17) / 0.00000 00000 00204 099D0/
858 DATA C1(18) /-0.00000 00000 00090 064D0/
859 DATA C1(19) / 0.00000 00000 00025 506D0/
860 DATA C1(20) /-0.00000 00000 00004 036D0/
861 DATA C1(21) /-0.00000 00000 00000 570D0/
862 DATA C1(22) / 0.00000 00000 00000 762D0/
863 DATA C1(23) /-0.00000 00000 00000 363D0/
864 DATA C1(24) / 0.00000 00000 00000 118D0/
865 DATA C1(25) /-0.00000 00000 00000 025D0/
866
867 DATA C2( 0) / 0.04655 77987 37516 4561D0/
868 DATA C2( 1) / 0.04499 21302 01239 4140D0/
869 DATA C2( 2) /-0.00175 42871 39651 4532D0/
870 DATA C2( 3) /-0.00014 65340 02581 0678D0/
871 DATA C2( 4) / 0.00003 91330 40863 0159D0/
872 DATA C2( 5) /-0.00000 34932 28659 7731D0/
873 DATA C2( 6) /-0.00000 03153 53003 2345D0/
874 DATA C2( 7) / 0.00000 01876 58200 8529D0/
875 DATA C2( 8) /-0.00000 00377 55280 4930D0/
876 DATA C2( 9) / 0.00000 00026 65516 5010D0/
877 DATA C2(10) / 0.00000 00010 88144 8122D0/
878 DATA C2(11) /-0.00000 00005 35500 7671D0/
879 DATA C2(12) / 0.00000 00001 31576 5447D0/
880 DATA C2(13) /-0.00000 00000 15286 0881D0/
881 DATA C2(14) /-0.00000 00000 03394 7646D0/
882 DATA C2(15) / 0.00000 00000 02702 0267D0/
883 DATA C2(16) /-0.00000 00000 00946 3142D0/
884 DATA C2(17) / 0.00000 00000 00207 1565D0/
885 DATA C2(18) /-0.00000 00000 00012 6931D0/
886 DATA C2(19) /-0.00000 00000 00013 9756D0/
887 DATA C2(20) / 0.00000 00000 00008 5929D0/
888 DATA C2(21) /-0.00000 00000 00003 1070D0/
889 DATA C2(22) / 0.00000 00000 00000 7515D0/
890 DATA C2(23) /-0.00000 00000 00000 0648D0/
891 DATA C2(24) /-0.00000 00000 00000 0522D0/
892 DATA C2(25) / 0.00000 00000 00000 0386D0/
893 DATA C2(26) /-0.00000 00000 00000 0165D0/
894 DATA C2(27) / 0.00000 00000 00000 0050D0/
895 DATA C2(28) /-0.00000 00000 00000 0009D0/
896
897 V=ABS(X)
898 IF(V .LT. 8) THEN
899 Y=R8*V
900 H=2*Y**2-1
901 ALFA=H+H
902 B1=0
903 B2=0
904 DO 4 I = NB,0,-1
905 B0=B(I)+ALFA*B1-B2
906 B2=B1
907 4 B1=B0
908 H=SQRT(V)*Y*(B0-B2)
909 ELSE
910 R=1/V
911 H=10*R-1
912 ALFA=H+H
913 B1=0
914 B2=0
915 DO 5 I = NC1,0,-1
916 B0=C1(I)+ALFA*B1-B2
917 B2=B1
918 5 B1=B0
919 S=B0-H*B2
920 B1=0
921 B2=0
922 DO 6 I = NC2,0,-1
923 B0=C2(I)+ALFA*B1-B2
924 B2=B1
925 6 B1=B0
926 H=C0-SQRT(R)*(S*COS(V)+(B0-H*B2)*SIN(V))
927 END IF
928 IF(X .LT. 0) H=-H
929 GO TO 9
930
931 ENTRY DFRCOS(X)
932
933 V=ABS(X)
934 IF(V .LT. 8) THEN
935 H=R32*V**2-1
936 ALFA=H+H
937 B1=0
938 B2=0
939 DO 1 I = NA,0,-1
940 B0=A(I)+ALFA*B1-B2
941 B2=B1
942 1 B1=B0
943 H=SQRT(V)*(B0-H*B2)
944 ELSE
945 R=1/V
946 H=10*R-1
947 ALFA=H+H
948 B1=0
949 B2=0
950 DO 2 I = NC1,0,-1
951 B0=C1(I)+ALFA*B1-B2
952 B2=B1
953 2 B1=B0
954 S=B0-H*B2
955 B1=0
956 B2=0
957 DO 3 I = NC2,0,-1
958 B0=C2(I)+ALFA*B1-B2
959 B2=B1
960 3 B1=B0
961 H=C0-SQRT(R)*((B0-H*B2)*COS(V)-S*SIN(V))
962 END IF
963 IF(X .LT. 0) H=-H
964 9 DFRSIN=H
965 RETURN
966 END
967
968c--#include "gen/pilot.h"
148aadde 969 FUNCTION CGAMMA(Z)
f5ab1a71 970
971 COMPLEX*8 CGAMMA
972 COMPLEX*8 Z,U,V,F,H,S
973 CHARACTER NAME*(*)
974 CHARACTER*80 ERRTXT
975 PARAMETER (NAME = 'CGAMMA')
f5ab1a71 976 DIMENSION C(0:15)
f5ab1a71 977 PARAMETER (Z1 = 1, HF = Z1/2)
978
979c--#if defined(CERNLIB_QF2C)
980c--#include "gen/gcmpfun.inc"
981c--#endif
982
983 DATA PI /3.14159 26535 89793 24D0/
984 DATA C1 /2.50662 82746 31000 50D0/
985
986 DATA C( 0) / 41.62443 69164 39068D0/
987 DATA C( 1) /-51.22424 10223 74774D0/
988 DATA C( 2) / 11.33875 58134 88977D0/
989 DATA C( 3) / -0.74773 26877 72388D0/
990 DATA C( 4) / 0.00878 28774 93061D0/
991 DATA C( 5) / -0.00000 18990 30264D0/
992 DATA C( 6) / 0.00000 00019 46335D0/
993 DATA C( 7) / -0.00000 00001 99345D0/
994 DATA C( 8) / 0.00000 00000 08433D0/
995 DATA C( 9) / 0.00000 00000 01486D0/
996 DATA C(10) / -0.00000 00000 00806D0/
997 DATA C(11) / 0.00000 00000 00293D0/
998 DATA C(12) / -0.00000 00000 00102D0/
999 DATA C(13) / 0.00000 00000 00037D0/
1000 DATA C(14) / -0.00000 00000 00014D0/
1001 DATA C(15) / 0.00000 00000 00006D0/
1002
1003c----#if !defined(CERNLIB_QF2C)
1004c----#include "gen/gcmpfun.inc"
1005c----#endif
1006
1007 COMPLEX*8 GREAL,GIMAG,XARG,YARG
1008
1009 COMPLEX*8 ZARG,GCONJG,GCMPLX
1010 GREAL( ZARG)=REAL( ZARG)
1011 GIMAG( ZARG)=AIMAG(ZARG)
1012 GCONJG(ZARG)=CONJG(ZARG)
1013c-- GCMPLX(XARG,YARG)= CMPLX(XARG,YARG)
1014
1015 U=Z
1016 X=U
1017 IF(GIMAG(U) .EQ. 0 .AND. -ABS(X) .EQ. INT(X)) THEN
1018 F=0
1019 H=0
1020 WRITE(ERRTXT,101) X
1021c- CALL MTLPRT(NAME,'C305.1',ERRTXT)
1022 ELSE
1023 IF(X .GE. 1) THEN
1024 F=1
1025 V=U
1026 ELSEIF(X .GE. 0) THEN
1027 F=1/U
1028 V=1+U
1029 ELSE
1030 F=1
1031 V=1-U
1032 END IF
1033 H=1
1034 S=C(0)
1035 DO 1 K = 1,15
1036 H=((V-K)/(V+(K-1)))*H
1037 1 S=S+C(K)*H
1038 H=V+(4+HF)
1039 H=C1*EXP((V-HF)*LOG(H)-H)*S
1040 IF(X .LT. 0) H=PI/(SIN(PI*U)*H)
1041 ENDIF
1042
1043c----#if !defined(CERNLIB_DOUBLE)
1044 CGAMMA=F*H
1045c----#endif
1046
1047 RETURN
1048 101 FORMAT('ARGUMENT EQUALS NON-POSITIVE INTEGER = ',1P,E15.1)
1049 END
1050
1051
1052 SUBROUTINE FSIPN !(WEIF)
1053C calculating particle-nucleus Coulomb Wave functions FFij
1054 IMPLICIT REAL*8 (A-H,O-Z)
1055 COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
1056 COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, !part. momenta in NRF
1057 1 P2X,P2Y,P2Z,E2,P2
1058 COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emis. points in NRF
1059 1 X2,Y2,Z2,T2,R2
1060 COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
1061 COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
1062 COMMON/FSI_ICH1/ICH1
1063 COMMON/FSI_ETA/ETA
1064 COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS
1065 COMMON/FSI_FFPN/FF12,FF21
1066 COMPLEX*16 FF1,FF12,FF21
1067 FF12=DCMPLX(1.D0,0.D0)
1068 FF21=DCMPLX(1.D0,0.D0)
1069 IF(I3C.EQ.0)RETURN
1070 ICH1=IDINT(C1)
1071 IF(ICH1.EQ.0)GOTO 11
1072 XH=AC1*P1
1073 ACH=ACP(XH)
1074 ACHR=DSQRT(ACH)
1075 ETA=0.D0
1076 IF(XH.NE.0.D0)ETA=1/XH
1077 RHOS=P1*R1
1078 HS=X1*P1X+Y1*P1Y+Z1*P1Z
1079 FF12=FF12*FF1(RHOS,HS)
1080 IF(IQS.EQ.0)GOTO 11
1081 RHOS=P1*R2
1082 HS=X2*P1X+Y2*P1Y+Z2*P1Z
1083 FF21=FF21*FF1(RHOS,HS)
1084 11 ICH1=IDINT(C2)
1085 IF(ICH1.EQ.0)GOTO 10
1086 XH=AC2*P2
1087 ACH=ACP(XH)
1088 ACHR=DSQRT(ACH)
1089 ETA=0.D0
1090 IF(XH.NE.0.D0)ETA=1/XH
1091 RHOS=P2*R2
1092 HS=X2*P2X+Y2*P2Y+Z2*P2Z
1093 FF12=FF12*FF1(RHOS,HS)
1094CW WRITE(6,41)'AC2 ',AC2,'ACH ',ACH,'ETA ',ETA,'RHOS ',RHOS,'HS ',HS
109541 FORMAT(5(A5,E11.4))
1096CW WRITE(6,40)'FF12 ',DREAL(FF12),DIMAG(FF12)
1097 IF(IQS.EQ.0)GOTO 10
1098 RHOS=P2*R1
1099 HS=X1*P2X+Y1*P2Y+Z1*P2Z
1100 FF21=FF21*FF1(RHOS,HS)
1101CW WRITE(6,41)'AC1 ',AC1,'ACH ',ACH,'ETA ',ETA,'RHOS ',RHOS,'HS ',HS
1102CW WRITE(6,40)'FF21 ',DREAL(FF21),DIMAG(FF21)
110340 FORMAT(A7,2E12.4)
1104 10 CONTINUE
1105
1106C WEIF = the weight due to the Coulomb particle-nucleus interaction
1107 WEIF=DREAL(FF12)**2+DIMAG(FF12)**2
1108 IF(IQS.EQ.1)WEIF=0.5D0*(WEIF+DREAL(FF21)**2+DIMAG(FF21)**2)
1109 RETURN
1110 END
1111
1112C=======================================================
1113C
1114 SUBROUTINE FSIWF !(WEI)
1115C==> Prepares necessary quantities and call VZ(WEI) to calculate
1116C the weight due to FSI
1117 IMPLICIT REAL*8 (A-H,O-Z)
1118 COMMON/FSI_CVK/V,CVK
1119 COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, !part. momenta in NRF
1120 1 P2X,P2Y,P2Z,E2,P2
1121 COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
1122 1 X,Y,Z,T,RP,RPS
1123 COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emis. points in NRF
1124 1 X2,Y2,Z2,T2,R2
1125 COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
1126 COMMON/FSI_SPIN/RHO(10)
1127 COMMON/FSI_BP/B,P
1128 COMMON/FSI_ETA/ETA
1129 COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
1130 COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10),
1131 1 SBKRB(10),SDKK(10)
1132 COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
1133 COMMON/FSI_RR/F(10)
1134 COMMON/FSI_FD/FD(10),RD(10)
1135 COMMON/FSI_C/C(10),AM,AMS,DM
1136 COMPLEX*16 C,F
1137 COMMON/FSI_AA/AA
1138 COMMON/FSI_SHH/SH,CHH
1139 COMMON/FSI_AAPI/AAPI(20,2)/FSI_AAND/AAND(20,4)
1140 COMMON/FSI_AAPIN/AAPIN(20,2)
1141 COMMON/FSI_P12/P12X,P12Y,P12Z,E12,P12,AM12,EPM
1142 COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS
1143
1144
1145cmlv
1146 IF(IRANPOS.EQ.0)THEN
1147C==>calculating relative 4-coordinates of the particles in PRF
1148C- {T,X,Y,Z} from the relative coordinates {TS,XS,YS,ZS} in NRF
1149 XS=X1-X2
1150 YS=Y1-Y2
1151 ZS=Z1-Z2
1152 TS=T1-T2
1153 RS12=XS*P12X+YS*P12Y+ZS*P12Z
1154 H1=(RS12/EPM-TS)/AM12
1155 X=XS+P12X*H1
1156 Y=YS+P12Y*H1
1157 Z=ZS+P12Z*H1
1158 T=(E12*TS-RS12)/AM12
1159 RPS=X*X+Y*Y+Z*Z
1160 RP=DSQRT(RPS)
1161c WRITE(6,38)'RP ',RP,'X ',X,Y,Z,T
1162c FORMAT(A7,E11.4,A7,4E11.4)
1163 ENDIF
1164
1165 CVK=(P12X*PPX+P12Y*PPY+P12Z*PPZ)/(P12*AK)
1166 V=P12/E12
1167
1168 IF(ICH.EQ.0)GOTO 21
1169 XH=AC*AK
1170 ACH=ACP(XH)
1171 ACHR=DSQRT(ACH)
1172 ETA=0.D0
1173 IF(XH.NE.0.D0)ETA=1/XH
1174C---HCP, HPR needed (e.g. in GST) if ICH=1
1175 HCP=HC(XH)
1176 HPR=HCP+.1544313298D0
1177 21 CONTINUE
1178 MSP=MSPIN
1179 DO 30 JJ=1,MSP
1180 ISPIN=JJ
1181 IF(NS.NE.1)GOTO22
1182C---Calc. quantities for the square well potential;
1183C-- for LL > 5 the square well potential is not possible or available
1184 IF(LL.EQ.4)GOTO 22
1185 BK(JJ)=DSQRT(EB(JJ)**2+AKS)
1186 XRA=2*RB(JJ)/AC
1187 HRA=BK(JJ)*RB(JJ)
1188 CALL SEQ(XRA,HRA)
1189 SBKRB(JJ)=HRA*B
1190 HRA=AK*RB(JJ)
1191 CALL GST(XRA,HRA)
1192 SDK(JJ)=SH
1193 CDK(JJ)=CHH
1194 SDKK(JJ)=RB(JJ)
1195 IF(AK.NE.0.D0)SDKK(JJ)=SH/AK
1196 IF(ICH.EQ.1)SDK(JJ)=ACH*SDK(JJ)
1197 22 CONTINUE
1198C-----------------------------------------------------------------------
1199C---Calc. the strong s-wave scattering amplitude = C(JJ)
1200C-- divided by Coulomb penetration factor squared (if ICH=1)
1201 IF(NS.NE.1)GOTO 230
1202 IF(LL.NE.4)GOTO 230 ! SW scat. amplitude used for alfa-alfa only
1203 GAK=G(AK)
1204 AKACH=AK
1205 IF(ICH.EQ.1)AKACH=AK*ACH
1206 C(JJ)=1/DCMPLX(GAK,-AKACH) ! amplitude for the SW-potential
1207 GOTO 30
1208 230 IF(LL.EQ.5.OR.LL.EQ.6.OR.LL.EQ.7)GOTO20 ! pipi
1209 IF(LL.EQ.12.OR.LL.EQ.13)GOTO20 ! piN
1210 IF(LL.EQ.8.OR.LL.EQ.9.OR.LL.EQ.18)GOTO20 ! Nd, dd
1211 IF(LL.EQ.14.OR.LL.EQ.17.OR.LL.EQ.23)GOTO27 ! K+K-, K-p, K0K0-b
1212 A1=RD(JJ)*FD(JJ)*AKS
1213 A2=1+.5D0*A1
1214 IF(ICH.EQ.1)A2=A2-2*HCP*FD(JJ)/AC
1215 AKF=AK*FD(JJ)
1216 IF(ICH.EQ.1)AKF=AKF*ACH
1217 C(JJ)=FD(JJ)/DCMPLX(A2,-AKF)
1218 GOTO30
1219 20 CONTINUE
1220C---Calc. scatt. ampl. C(JJ) for pipi, piN and Nd, dd
1221 JH=LL-7+2*JJ-2
1222 IF(LL.EQ.8.OR.LL.EQ.9)GPI2=GND(AKS,JH)
1223 IF(LL.EQ.18)GPI2=GDD(AKS,JJ)
1224 IF(LL.EQ.5.OR.LL.EQ.6.OR.LL.EQ.7)GPI2=GPIPI(AKS,2)
1225 IF(LL.EQ.12.OR.LL.EQ.13)GPI2=GPIN(AKS,2)
1226 C(JJ)=1.D0/DCMPLX(GPI2,-AK) !pi+pi+, nd, pd, pi+p, dd
1227 IF(LL.NE.5.AND.LL.NE.6.AND.LL.NE.13)GOTO27
1228 IF(LL.EQ.5.OR.LL.EQ.6)GPI1=GPIPI(AKS,1)
1229 IF(LL.EQ.13)GPI1=GPIN(AKS,1)
1230 IF(LL.EQ.5.OR.LL.EQ.13)
1231 c C(JJ)=.6667D0/DCMPLX(GPI1,-AK)+.3333D0*C(JJ) !pi+pi-,pi-p
1232 IF(LL.EQ.6)C(JJ)=.3333D0/DCMPLX(GPI1,-AK)+.6667D0*C(JJ) !pi0pi0
1233 27 CONTINUE
1234C---Calc. K+K-, K0K0-b or K-p s-wave scatt. ampl.
1235 IF(LL.EQ.14.OR.LL.EQ.23)CALL CKKB
1236 IF(LL.EQ.17)C(JJ)=DCMPLX(3.29D0,3.55D0)
1237C---Calc. pi+pi-, pi+pi+, pd, pi+p, pi-p, K+K- or K-p s-wave scatt. ampl.
1238C-- divided by Coulomb penetration factor squared (if ICH=1)
1239 IF(ICH.EQ.0)GOTO 30
1240 AAK=ACH*AK
1241 HCP2=2*HCP/AC
1242 C(JJ)=1/(1/C(JJ)-HCP2+DCMPLX(0.D0,AK-AAK))
1243c write(*,*)'c(jj)',c(jj)
1244 30 CONTINUE
1245C***********************************************************************
1246c write(*,*)'before call vz in fsiwf wei',wei
1247 CALL VZ !(WEI)
1248 RETURN
1249 END