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