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