Updates in GRP Preprocessor (Ernesto)
[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
f5ab1a71 125 COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS
ff4431bb 126 COMPLEX*16 FF12,FF21
f5ab1a71 127C------------------------------------------------------------------
f5ab1a71 128C==> 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
134C-----------------------------------------------------------
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 149C=======================================================================
150
151 SUBROUTINE LTRAN(P0,P,PS)
152C==>calculating particle 4-momentum PS={PSX,PSY,PSZ,ES}
153C in rest frame of a system 0 with 4-momentum P0={P0X,P0Y,P0Z,E0}
154C 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)
158C-----------------------------------------------------------------------
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)
172C==>calculating particle 4-momentum P={PX,PY,PZ,E}
173C from its 4-momentum PS={PSX,PSY,PSZ,ES}
174C 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)
178C-----------------------------------------------------------------------
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)
193C--- GDD = k*COTG(DELTA), X=k^2
194C--- 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
224C 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)
266CW WRITE(6,41)'AC2 ',AC2,'ACH ',ACH,'ETA ',ETA,'RHOS ',RHOS,'HS ',HS
26741 FORMAT(5(A5,E11.4))
268CW 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)
273CW WRITE(6,41)'AC1 ',AC1,'ACH ',ACH,'ETA ',ETA,'RHOS ',RHOS,'HS ',HS
274CW WRITE(6,40)'FF21 ',DREAL(FF21),DIMAG(FF21)
27540 FORMAT(A7,2E12.4)
276 10 CONTINUE
277
278C 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)
285C--- GPIPI = k*COTG(DELTA), X=k^2
286C-- 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)
298C--- GND = k*COTG(DELTA), X=k^2
299C--- J=1(2) corresp. to nd(pd), S=1/2,
300C--- 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
333C
334 SUBROUTINE FSIWF
335C==> Prepares necessary quantities and call VZ(WEI) to calculate
336C 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
364C==>calculating relative 4-coordinates of the particles in PRF
365C- {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
373c-mlv X=XS+P12X*H1
374c-mlv Y=YS+P12Y*H1
375c-mlv Z=ZS+P12Z*H1
376c-mlv T=(E12*TS-RS12)/AM12
377c-mlv RPS=X*X+Y*Y+Z*Z
378c-mlv RP=DSQRT(RPS)
379c WRITE(6,38)'RP ',RP,'X ',X,Y,Z,T
380c38 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
391C---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
398C---Calc. quantities for the square well potential;
399C-- for LL > 5 the square well potential is not possible or available
400cc 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
414C-----------------------------------------------------------------------
415C---Calc. the strong s-wave scattering amplitude = C(JJ)
416C-- 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
434C---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
444C---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)
447C---Calc. pi+pi-, pi+pi+, pd, K+K- or K-p s-wave scatt. amplitude
448C-- 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
454C***********************************************************************
455 CALL VZ
456 RETURN
457 END
458
459 SUBROUTINE VZ
f5ab1a71 460C==> 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
488C---Calc. EIDC=exp(i*Coul.Ph.);
489C-- 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
494c 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
515CW WRITE(6,38)'JJ ',JJ,'F ',DREAL(F(JJ)),DIMAG(F(JJ))
516CW WRITE(6,38)'JJ ',JJ,'C ',DREAL(C(JJ)),DIMAG(C(JJ))
517CW WRITE(6,38)'JJ ',JJ,'G ',DREAL(G),DIMAG(G)
518CW WRITE(6,38)'JJ ',JJ,'F12+G ',DREAL(F12+G),DIMAG(F12+G)
519CW WRITE(6,38)'JJ ',JJ,'F21+G ',DREAL(F21+G),DIMAG(F21+G)
52038 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
537c [M %/WRITE(*,*)WEI
538
f5ab1a71 539 RETURN
540 END
ff4431bb 541
542
f5ab1a71 543 SUBROUTINE FIRT
544C---CALC. THE F(JJ)
545C-- 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
576C---F= ASYMPTOTIC FORMULA (T= 0 APPROX.); NS= 4
577 6 F(JJ)=DCMPLX(CHH,SHH)
578 IF(NS.NE.1)GOTO 10
579C---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
590C---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)
663C---CALC. FUNCTIONS B, P (EQS. (17) OF G-K-L-L);
664C-- 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)
691C---CALC. FUNCTIONS CHH=REAL(GST), SH=IMAG(GST)/ACH, B=SH/H
692C-- 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)
709C---Calc. F12, F21;
710C-- F12= FF0* plane wave, FF0=F*ACHR,
711C---F is the confluent hypergeometric function,
712C-- 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)
731C-- FAS=F*ACHR
732C---F is the confluent hypergeometric function at k*r >> 1
733C-- 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)
747C-- FF0=F*ACHR
748C-- F is the confluent hypergeometric function
749C-- (Eq. (15) of G-K-L-L), F= 1 at r* << AC
750C-- 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
760CC 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)
779C---HC = h-function of Landau-Lifshitz: h(x)=Re[psi(1-i/x)]+ln(x)
780C-- psi(x) is the digamma function (the logarithmic derivative of
781C-- 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
791CC 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
798C---Provides 7 digit accuracy
799 HC=S-.5772156649D0+DLOG(X)
800 RETURN
801CC 2 HC=1.2D0/X**2+DLOG(X)-.5772156649 D0
802CC 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)
814C--- 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)
824C---CALC. THE CONFL. HYPERGEOM. F-N = CHH+i*SH
825C-- 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)
850C---FF1=FF0; used for particle-nucleus system
851C-- FF0=F12*ACHR
852C-- F12 is the confluent hypergeometric function
853C-- (Eq. (15) of G-K-L-L), F12= 1 at r* << AC
854C-- 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
867C---Calc. EIDC=exp(i*Coul.Ph.);
868C-- 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)
877C---Used to calculate SW scattering amplitude for alpa-alpha system
f5ab1a71 878C---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))
915C 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)
917C GG= [B'(kp,d)*ReG(k,d)-B(kp,d)*ZZ]/d
918 G=GG/(BRHO*BPR-BRHOP*BRHOS)
919C G= GG/[B(kp,d)*B'(k,d)-B'(kp,d)*B(k,d)]
920 RETURN
921 END
f5ab1a71 922 SUBROUTINE DERIW(X,H)
923C---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 944C================================================================
945
f5ab1a71 946
947c-------------#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
1123c--#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
1137c--#if defined(CERNLIB_QF2C)
1138c--#include "gen/gcmpfun.inc"
1139c--#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
1161c----#if !defined(CERNLIB_QF2C)
1162c----#include "gen/gcmpfun.inc"
1163c----#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)
1171c-- 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
1179c- 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
1201c----#if !defined(CERNLIB_DOUBLE)
1202 CGAMMA=F*H
1203c----#endif
1204
1205 RETURN
1206 101 FORMAT('ARGUMENT EQUALS NON-POSITIVE INTEGER = ',1P,E15.1)
1207 END
1208