PWG2/FEMTOSCOPY/AliFemto/AliFemtoFsiWeightLednicky.F

   1 c         1         2         3         4         5         6         7         8
   2 c---------|---------|---------|---------|---------|---------|---------|---------|
   3 *-- Author :    R.Lednicky      20/01/95
   4       SUBROUTINE FSIW(J,WEIF,WEI,WEIN)
   5
   6 C=======================================================================
   7 C    Calculates final state interaction (FSI) weights
   8 C    WEIF = weight due to particle - (effective) nucleus FSI (p-N)
   9 C    WEI  = weight due to p-p-N FSI
  10 C    WEIN = weight due to p-p FSI; note that WEIN=WEI if I3C=0;
  11 C                                  note that if I3C=1 the calculation of
  12 C                                  WEIN can be skipped by putting J=0
  13 C.......................................................................
  14 C    Correlation Functions:
  15 C    CF(p-p-N)   = sum(WEI)/sum(WEIF)
  16 C    CF(p-p)     = sum(WEIN)/sum(1); here the nucleus is completely
  17 C                                    inactive
  18 C    CF(p-p-"N") = sum(WEIN*WEIF')/sum(WEIF'), where WEIN and WEIF'
  19 C                  are not correlated (calculated at different emission
  20 C                  points, e.g., for different events);
  21 C                  thus here the nucleus affects one-particle
  22 C                  spectra but not the correlation
  23 C.......................................................................
  24 C    User must supply data file <fn> on unit NUNIT (e.g. =11) specifying
  25 C    LL   : particle pair
  26 C    NS   : approximation used to calculate Bethe-Salpeter amplitude
  27 C    ITEST: test switch
  28 C           If ITEST=1 then also following parameters are required
  29 C    ICH  : 1(0) Coulomb interaction between the two particles ON (OFF)
  30 C    IQS  : 1(0) quantum statistics for the two particles ON (OFF)
  31 C    ISI  : 1(0) strong interaction between the two particles ON (OFF)
  32 C    I3C  : 1(0) Coulomb interaction with residual nucleus ON (OFF)
  33 C    This data file can contain other information useful for the user.
  34 C    It is read by subroutines READINT4 and READREA8(4) (or READ_FILE).
  35 C----------------------------------------------------------------------
  36 C-   LL       1  2  3  4  5   6   7   8  9 10  11  12  13  14 15 16 17
  37 C-   part. 1: n  p  n  a  pi+ pi0 pi+ n  p pi+ pi+ pi+ pi- K+ K+ K+ K-
  38 C-   part. 2: n  p  p  a  pi- pi0 pi+ d  d  K-  K+  p   p  K- K+ p  p
  39 C   NS=1 y/n: +  +  +  +  +   -   -   -  -  -   -   -   -  -  -  -  -
  40 C----------------------------------------------------------------------
  41 C-   LL       18 19 20 21 22 23  24 25 26 27 28 29 30
  42 C-   part. 1: d  d  t  t  K0 K0  d  p  p  p  n  lam p
  43 C-   part. 2: d  a  t  a  K0 K0b t  t  a  lam lam lam pb
  44 C   NS=1 y/n: -  -  -  -  -  -   -  -  -  +  +  +  -
  45 C----------------------------------------------------------------------
  46 C   NS=1  Square well potential,
  47 C   NS=3  not used
  48 C   NS=4  scattered wave approximated by the spherical wave,
  49 C   NS=2  same as NS=4 but the approx. of equal emission times in PRF
  50 C         not required (t=0 approx. used in all other cases).
  51 C   Note: if NS=2,4, the B-S amplitude diverges at zero distance r* in
  52 C         the two-particle c.m.s.; user can specify a cutoff AA in
  53 C         SUBROUTINE FSIINI, for example:
  54 C         IF(NS.EQ.2.OR.NS.EQ.4)AA=5.D0 !! in 1/GeV --> AA=1. fm
  55 C---------------------------------------------------------------------
  56 C    ITEST=1 any values of parameters ICH, IQS, ISI, I3C are allowed
  57 C            and should be given in data file <fn>
  58 C    ITEST=0 physical values of these parameters are put automatically
  59 C            in FSIINI (their values are not required in data file)
  60 C=====================================================================
  61 C    At the beginning of calculation user should call FSIINI,
  62 C    which reads LL, NS, ITEST (and eventually ICH, IQS, ISI, I3C)
  63 C    and initializes various parameters.
  64 C    In particular the constants in
  65 C      COMMON/FSI_CONS/PI,PI2,SPI,DR,W
  66 C    may be useful for the user:
  67 C     W=1/.1973D0    ! from fm to 1/GeV
  68 C     PI=4*DATAN(1.D0)
  69 C     PI2=2*PI
  70 C     SPI=DSQRT(PI)
  71 C     DR=180.D0/PI   ! from radian to degree
  72 C      _______________________________________________________
  73 C  !! |Important note: all real quantities are assumed REAL*8 | !!
  74 C      -------------------------------------------------------
  75 C    For each event user should fill in the following information
  76 C    in COMMONs (all COMMONs in FSI calculation start with FSI_):
  77 C    ...................................................................
  78 C     COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
  79 C    Only
  80 C         AMN  = mass of the effective nucleus   [GeV/c**2]
  81 C         CN   = charge of the effective nucleus [elem. charge units]
  82 C    are required
  83 C    ...................................................................
  84 C     COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, !part. momenta in the rest frame
  85 C    1               P2X,P2Y,P2Z,E2,P2  !of effective nucleus (NRF)
  86 C    Only the components
  87 C                        PiX,PiY,PiZ  [GeV/c]
  88 C    in NRF are required.
  89 C    To make the corresponding Lorentz transformation user can use the
  90 C    subroutines LTRAN and LTRANB
  91 C    ...................................................................
  92 C    COMMON/FSI_COOR/X1,Y1,Z1,T1,R1,     ! 4-coord. of emission
  93 C    1               X2,Y2,Z2,T2,R2      ! points in NRF
  94 C    The componets
  95 C                       Xi,Yi,Zi  [fm]
  96 C    and emission times
  97 C                          Ti   [fm/c]
  98 C    should be given in NRF with the origin assumed at the center
  99 C    of the effective nucleus. If the effect of residual nucleus is
 100 C    not calculated within FSIW, the NRF can be any fixed frame.
 101 C-----------------------------------------------------------------------
 102 C    Before calling FSIW the user must call
 103 C     CALL LTRAN12
 104 C    Besides Lorentz transformation to pair rest frame:
 105 C    (p1-p2)/2 --> k* it also transforms 4-coordinates of
 106 C    emission points from fm to 1/GeV and calculates Ei,Pi and Ri.
 107 C    Note that |k*|=AK in COMMON/FSI_PRF/
 108 C-----------------------------------------------------------------------
 109 C    After making some additional filtering using k* (say k* < k*max)
 110 C    or direction of vector k*,
 111 C    user can finally call FSIW to calculate the FSI weights
 112 C    to be used to construct the correlation function
 113 C=======================================================================
 114
 115       IMPLICIT REAL*8 (A-H,O-Z)
 116       COMMON/FSI_JR/JRAT
 117       COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
 118       COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1,  ! particle momenta in NRF
 119      1               P2X,P2Y,P2Z,E2,P2
 120       COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS, ! k*=(p1-p2)/2 and x1-x2
 121      1               X,Y,Z,T,RP,RPS      ! in pair rest frame (PRF)
 122       COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, !4-coord. of emis. points in NRF
 123      1                X2,Y2,Z2,T2,R2
 124       COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
 125       COMMON/FSI_FFPN/FF12,FF21
 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       CALL FSIPN(WEIF) !weight due to particle-nucleus FSI
 136       JRAT=0
 137       call BoostToPrf()
 138       CALL FSIWF(WEI)  !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(WEIN) ! weight due to particle-particle FSI
 145          ENDIF
 146       RETURN
 147       END
 148 C=======================================================================
 149
 150       SUBROUTINE LTRAN(P0,P,PS)
 151 C==>calculating particle 4-momentum PS={PSX,PSY,PSZ,ES}
 152 C   in rest frame of a system 0 with 4-momentum P0={P0X,P0Y,P0Z,E0}
 153 C   from its 4-momentum P={PX,PY,PZ,E}
 154
 155       IMPLICIT REAL*8 (A-H,O-Z)
 156       DIMENSION P0(4),P(4),PS(4)
 157 C-----------------------------------------------------------------------
 158       P0S=P0(1)**2+P0(2)**2+P0(3)**2
 159       AM0=DSQRT(P0(4)**2-P0S)
 160       EPM=P0(4)+AM0
 161       PP0=P(1)*P0(1)+P(2)*P0(2)+P(3)*P0(3)
 162       H=(PP0/EPM-P(4))/AM0
 163       PS(1)=P(1)+P0(1)*H
 164       PS(2)=P(2)+P0(2)*H
 165       PS(3)=P(3)+P0(3)*H
 166       PS(4)=(P0(4)*P(4)-PP0)/AM0
 167       RETURN
 168       END
 169
 170       SUBROUTINE LTRANB(P0,PS,P)
 171 C==>calculating particle 4-momentum P={PX,PY,PZ,E}
 172 C   from its 4-momentum PS={PSX,PSY,PSZ,ES}
 173 C   in rest frame of a system 0 with 4-momentum P0={P0X,P0Y,P0Z,E0}
 174
 175       IMPLICIT REAL*8 (A-H,O-Z)
 176       DIMENSION P0(4),P(4),PS(4)
 177 C-----------------------------------------------------------------------
 178       P0S=P0(1)**2+P0(2)**2+P0(3)**2
 179       AM0=DSQRT(P0(4)**2-P0S)
 180       EPM=P0(4)+AM0
 181       PSP0=PS(1)*P0(1)+PS(2)*P0(2)+PS(3)*P0(3)
 182       HS=(PSP0/EPM+PS(4))/AM0
 183       P(1)=PS(1)+P0(1)*HS
 184       P(2)=PS(2)+P0(2)*HS
 185       P(3)=PS(3)+P0(3)*HS
 186       P(4)=(P0(4)*PS(4)+PSP0)/AM0
 187       RETURN
 188       END
 189
 190       SUBROUTINE LTRAN12
 191 C==>calculating particle momentum in PRF {EE,PPX,PPY,PPZ} from
 192 C-  the momentum of the first particle {E1,P1X,P1Y,P1Z) in NRF
 193       IMPLICIT REAL*8 (A-H,O-Z)
 194       COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1,  !part. momenta in NRF
 195      1               P2X,P2Y,P2Z,E2,P2
 196       COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
 197      1               X,Y,Z,T,RP,RPS
 198       COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
 199       COMMON/FSI_P12/P12X,P12Y,P12Z,E12,P12,AM12,EPM
 200       COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, !4-coord. of emis. points in NRF
 201      1                X2,Y2,Z2,T2,R2
 202       COMMON/FSI_CONS/PI,PI2,SPI,DR,W
 203
 204 C   fm --> 1/GeV
 205       X1=X1*W
 206       Y1=Y1*W
 207       Z1=Z1*W
 208       T1=T1*W
 209       X2=X2*W
 210       Y2=Y2*W
 211       Z2=Z2*W
 212       T2=T2*W
 213 C   calculating Ri, Pi and Ei
 214       R1=DSQRT(X1*X1+Y1*Y1+Z1*Z1)
 215       R2=DSQRT(X2*X2+Y2*Y2+Z2*Z2)
 216       P1S=P1X*P1X+P1Y*P1Y+P1Z*P1Z
 217       P2S=P2X*P2X+P2Y*P2Y+P2Z*P2Z
 218       P1=DSQRT(P1S)
 219       P2=DSQRT(P2S)
 220       E1=DSQRT(AM1*AM1+P1S)
 221       E2=DSQRT(AM2*AM2+P2S)
 222 C-----------------------------------------------------------------------
 223       E12=E1+E2
 224       P12X=P1X+P2X
 225       P12Y=P1Y+P2Y
 226       P12Z=P1Z+P2Z
 227       P12S=P12X**2+P12Y**2+P12Z**2
 228       AM12=DSQRT(E12**2-P12S)
 229       EPM=E12+AM12
 230       P12=DSQRT(P12S)
 231       P112=P1X*P12X+P1Y*P12Y+P1Z*P12Z
 232       H1=(P112/EPM-E1)/AM12
 233       PPX=P1X+P12X*H1
 234       PPY=P1Y+P12Y*H1
 235       PPZ=P1Z+P12Z*H1
 236       EE=(E12*E1-P112)/AM12
 237       AKS=EE**2-AM1**2
 238       AK=DSQRT(AKS)
 239
 240 CW      WRITE(6,38)'AK ',AK,'K ',PPX,PPY,PPZ,EE
 241 38    FORMAT(A7,E11.4,A7,4E11.4)
 242       RETURN
 243       END
 244
 245       SUBROUTINE FSIPN(WEIF)
 246 C  calculating particle-nucleus Coulomb Wave functions FFij
 247       IMPLICIT REAL*8 (A-H,O-Z)
 248       COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
 249       COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1,  !part. momenta in NRF
 250      1               P2X,P2Y,P2Z,E2,P2
 251       COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emis. points in NRF
 252      1                X2,Y2,Z2,T2,R2
 253       COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
 254       COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
 255       COMMON/FSI_ICH1/ICH1
 256       COMMON/FSI_ETA/ETA
 257       COMMON/FSI_FFPN/FF12,FF21
 258       COMPLEX*16 FF1,FF12,FF21
 259       FF12=DCMPLX(1.D0,0.D0)
 260       FF21=DCMPLX(1.D0,0.D0)
 261 C      ACH=1.D0
 262 C      WEIF=1.D0
 263       IF(I3C.EQ.0)RETURN
 264       ICH1=IDINT(C1)
 265       IF(ICH1.EQ.0)GOTO 11
 266       XH=AC1*P1
 267       ACH=ACP(XH)
 268       ACHR=DSQRT(ACH)
 269       ETA=0.D0
 270       IF(XH.NE.0.D0)ETA=1/XH
 271       RHOS=P1*R1
 272       HS=X1*P1X+Y1*P1Y+Z1*P1Z
 273       FF12=FF12*FF1(RHOS,HS)
 274       IF(IQS.EQ.0)GOTO 11
 275       RHOS=P1*R2
 276       HS=X2*P1X+Y2*P1Y+Z2*P1Z
 277       FF21=FF21*FF1(RHOS,HS)
 278  11   ICH1=IDINT(C2)
 279       IF(ICH1.EQ.0)GOTO 10
 280       XH=AC2*P2
 281       ACH=ACP(XH)
 282       ACHR=DSQRT(ACH)
 283       ETA=0.D0
 284       IF(XH.NE.0.D0)ETA=1/XH
 285       RHOS=P2*R2
 286       HS=X2*P2X+Y2*P2Y+Z2*P2Z
 287       FF12=FF12*FF1(RHOS,HS)
 288 CW      WRITE(6,41)'AC2 ',AC2,'ACH ',ACH,'ETA ',ETA,'RHOS ',RHOS,'HS ',HS
 289 41    FORMAT(5(A5,E11.4))
 290 CW      WRITE(6,40)'FF12 ',DREAL(FF12),DIMAG(FF12)
 291       IF(IQS.EQ.0)GOTO 10
 292       RHOS=P2*R1
 293       HS=X1*P2X+Y1*P2Y+Z1*P2Z
 294       FF21=FF21*FF1(RHOS,HS)
 295 CW      WRITE(6,41)'AC1 ',AC1,'ACH ',ACH,'ETA ',ETA,'RHOS ',RHOS,'HS ',HS
 296 CW      WRITE(6,40)'FF21 ',DREAL(FF21),DIMAG(FF21)
 297 40    FORMAT(A7,2E12.4)
 298  10   CONTINUE
 299
 300 C  WEIF = the weight due to the Coulomb particle-nucleus interaction
 301       WEIF=DREAL(FF12)**2+DIMAG(FF12)**2
 302       IF(IQS.EQ.1)WEIF=0.5D0*(WEIF+DREAL(FF21)**2+DIMAG(FF21)**2)
 303       RETURN
 304       END
 305
 306       FUNCTION GPIPI(X,J)
 307 C--- GPIPI = k*COTG(DELTA), X=k^2
 308 C--  J=1(2) corresponds to isospin=0(2)
 309       IMPLICIT REAL*8 (A-H,O-Z)
 310       COMMON/FSI_AAPI/AAPI(20,2)
 311       COMMON/FSI_C/HELP(20),AM,AMS,DM
 312       OM=DSQRT(X+AMS)
 313       XX=X/AMS
 314       GPIPI=OM/AAPI(1,J)
 315       GPIPI=GPIPI*(1+(AAPI(3,J)-AAPI(1,J)**2)*XX+AAPI(4,J)*XX*XX)
 316       GPIPI=GPIPI/(1+(AAPI(3,J)+AAPI(2,J)/AAPI(1,J))*XX)
 317       RETURN
 318       END
 319
 320       FUNCTION GPIN(X,J)
 321 C--- GPIN = k*COTG(DELTA), X=k^2
 322 C--  J=1(2) corresponds to piN isospin=1/2(3/2)
 323       IMPLICIT REAL*8 (A-H,O-Z)
 324       COMMON/FSI_AAPIN/AAPIN(20,2)
 325       GPIN=1/AAPIN(1,J)+.5D0*AAPIN(2,J)*X
 326       RETURN
 327       END
 328
 329       FUNCTION GND(X,J)
 330 C--- GND = k*COTG(DELTA), X=k^2
 331 C--- J=1(2) corresp. to nd(pd), S=1/2,
 332 C--- J=3(4) corresp. to nd(pd), S=3/2
 333       IMPLICIT REAL*8 (A-H,O-Z)
 334       COMMON/FSI_AAND/AAND(20,4)
 335       XX=X
 336       GND=1/AAND(1,J)+.5D0*AAND(2,J)*X
 337       DO 1 I=4,4
 338       XX=XX*X
 339    1  GND=GND+AAND(I,J)*XX
 340       GND=GND/(1+AAND(3,J)*X)
 341       RETURN
 342       END
 343       FUNCTION GDD(X,J)
 344 C--- GDD = k*COTG(DELTA), X=k^2
 345 C--- J=1,2,3 corresp. to S=0,1,2
 346       IMPLICIT REAL*8 (A-H,O-Z)
 347       COMMON/FSI_AADD/AADD(20,3)
 348       COMMON/FSI_C/C(10),AM,AMS,DM
 349       COMMON/FSI_CONS/PI,PI2,SPI,DR,W
 350       COMPLEX*16 C
 351       E=X/2/AM
 352       ER=DSQRT(E)
 353       IF(J.EQ.1)THEN
 354        GDD=ER*(AADD(1,1)*DEXP(-E/AADD(2,1))-AADD(3,1))
 355        GDD=GDD/DR   ! from degree to radian
 356        TAND=DTAN(GDD)
 357        IF(TAND.EQ.0.D0)TAND=1.D-10
 358        GDD=DSQRT(X)/TAND
 359       END IF
 360       IF(J.EQ.2)THEN
 361        GDD=1.D10
 362       END IF
 363       IF(J.EQ.3)THEN
 364        GDD=ER*(AADD(1,3)+AADD(2,3)*E)
 365        GDD=GDD/DR    ! from degree to radian
 366        TAND=DTAN(GDD)
 367        IF(TAND.EQ.0.D0)TAND=1.D-10
 368        GDD=DSQRT(X)/TAND
 369       END IF
 370       RETURN
 371       END
 372       BLOCK DATA
 373       IMPLICIT REAL*8 (A-H,O-Z)
 374       COMMON/FSI_AAPI/AAPI(20,2)/FSI_AAND/AAND(20,4)
 375       COMMON/FSI_AADD/AADD(20,3)/FSI_AAPIN/AAPIN(20,2)
 376       COMMON/FSI_AAKK/AAKK(9)/FSI_AAPAP/AAPAPR(3,2),AAPAPI(3,2)
 377 C---Parameters for GPIPI (I,J), J=1,2 -> isospin=0,2
 378       DATA AAPI/.2600D00, .2500D00, .3480D00,-.0637D00, 16*0.D0,
 379      1         -.0280D00,-.0820D00, .2795D00,-.0086D00, 16*0.D0/
 380 C---Parameters for GPIN (I,J), J=1,2 -> piN isospin=1/2,3/2
 381       DATA AAPIN/ .12265D1, .1563D2,     18*0.D0,
 382      1            -.750D0,  .7688D2,     18*0.D0/
 383 C---Parameters for GND (I,J), J=1(2) corresp. to nd(pd), S=1/2,
 384 C---                          J=3(4) corresp. to nd(pd), S=3/2
 385       DATA AAND/-.3295D1,-.8819D3, .4537D4, .28645D5, 16*0.D0,
 386      1          -.13837D2, .11505D2, .0D0, .10416D2,  16*0.D0,
 387      2          -.32180D2, .1014D2,  .0D0, .0D0,      16*0.D0,
 388      3          -.60213D2, .1333D2,  .0D0,-.70309D2,  16*0.D0/
 389       DATA AADD/ .10617D4, .3194D-2, .56849D3, 17*0.D0,
 390      1           20*0.D0,
 391      2          -.1085D4, .1987D5, 18*0.D0/
 392 C--- AAKK= m_K^2, m_pi^2, m_eta^2, m_S*^2, m_delta^2,
 393 C---       gam(S*-->KK-b), gam(S*-->pipi), gam(delta-->KK-b),
 394 C---       gam(delta-->pi eta)
 395       DATA AAKK/.247677D00,.01947977D00,.2997015D00,.9604D00,
 396      1          .96511D00,
 397 cc     2          .792D00, .199D00, .333D00, .222D00/ ! Martin (77)
 398      2          .094D00, .110D00, .333D00, .222D00/ ! Morgan (93)
 399 C---Parameters for PAP (I,J), j=1,2 -> isospin I=0,2
 400 C---                          i=1-3 -> a_singlet, a_triplet, d [fm]
 401 C---    Im a_IS (I=isospin, S=spin) are fixed by atomic data and
 402 C       n-bar survival time up to one free parameter, e.g. Im a_00
 403 C---    Batty (89), Kerbikov (93):
 404 C--- Ima_10=1.96-Ima_00, Ima_01=0.367-Ima_00/3, Ima_11=0.453+Ima_00/3
 405 C---       In DATA we put Ima_00=0.3.
 406 C---    Re a_IS are fixed by atomic data up to three free parameters
 407 C---    Batty (89):
 408 C---         Rea_aver(pp-bar)=Re[(a_00+a_10)+3(a_01+a_11)]/8=-0.9
 409 C---       In DATA we used Rea_IS from Paris potential Pignone (94)
 410 C---       rescaled by 1.67 to satisfy the atomic constraint.
 411 C---    Effective radius is is taken independent of IS from the phase
 412 C---    shift fit by Pirner et al. (91).
 413       DATA AAPAPR/-0.94D0, -1.98D0,  .1D0,
 414      1            -1.40D0,  0.37D0,  .1D0/ ! Re
 415       DATA AAPAPI/ 0.3 D0,  .267D0,-.01D0,
 416      1             1.66D0,  .553D0,-.01D0/ ! Im
 417       END
 418       SUBROUTINE CKKB  ! calculates KK-b scattering amplitude,
 419                        ! saturated by S*(980) and delta(982) resonances
 420       IMPLICIT REAL*8 (A-H,O-Z)
 421       COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
 422      1               X,Y,Z,T,RP,RPS
 423       COMMON/FSI_AAKK/AAKK(9)
 424       COMMON/FSI_C/C(10),AM,AMS,DM
 425       COMPLEX*16 C
 426       S4=AKS+AAKK(1)
 427       S=4*S4
 428       AKPIPI=DSQRT(S4-AAKK(2))
 429       EETA2=(S+AAKK(3)-AAKK(2))**2/4/S
 430       AKPIETA=DSQRT(EETA2-AAKK(3))
 431       C(1)=AAKK(6)/2/DCMPLX(AAKK(4)-S,
 432      ,-AK*AAKK(6)-AKPIPI*AAKK(7))
 433       C(1)=C(1)+AAKK(8)/2/DCMPLX(AAKK(5)-S,
 434      ,-AK*AAKK(8)-AKPIETA*AAKK(9))
 435       RETURN
 436       END
 437       SUBROUTINE CPAP           ! calculates pp-bar scattering amplitude
 438                                 ! accounting for nn-bar->pp-bar channel
 439       IMPLICIT REAL*8 (A-H,O-Z)
 440       COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
 441       COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
 442      1     X,Y,Z,T,RP,RPS
 443       COMMON/FSI_AAPAP/AAPAPR(3,2),AAPAPI(3,2)
 444       COMMON/FSI_C/C(10),AM,AMS,DM
 445       COMMON/FSI_2CHA/AK2,AK2S,AAK,HCP2,AMU2_AMU1 ! k* (kappa) for 2-nd channel
 446       COMPLEX*16 C
 447       DATA AM2/.93956563D0/
 448       AMU2_AMU1=AM2/AM          ! AM2=2*mu(nn-bar), AM=2*mu(pp-bar)
 449       AK2S0=AMU2_AMU1*AKS
 450       AK2S =AK2S0-2*AM2*(AM2-AM)
 451       IF(AK2S.GE.0.D0)THEN
 452          AK2=DCMPLX(DSQRT(AK2S),0.D0) !  k2
 453       ELSE
 454          AK2=DCMPLX(0.D0,DSQRT(-AK2S)) !  kappa2
 455       ENDIF
 456       C(10)=C(6+(ISPIN-1)*2)+
 457      +     DCMPLX(AAPAPR(3,ISPIN)*AKS/2-0.016D0-HCP2,
 458      ,     AAPAPI(3,ISPIN)*AKS/2-AAK) ! (1/f)11
 459       C(5)=C(6+(ISPIN-1)*2)+
 460      +     DCMPLX(AAPAPR(3,ISPIN)*AK2S0/2,
 461      ,     AAPAPI(3,ISPIN)*AK2S0/2)
 462       IF(AK2S.GE.0.D0)THEN
 463          C(5)=C(5)-DCMPLX(0.D0,AK2)
 464       ELSE
 465          C(5)=C(5)+DCMPLX(AK2,0.D0) ! (1/f)22
 466       ENDIF
 467       C(10)=C(10)*C(5)-C(7+(ISPIN-1)*2)*C(7+(ISPIN-1)*2)
 468       C(ISPIN)=C(5)/C(10)       ! f11
 469       C(ISPIN+2)=-C(7+(ISPIN-1)*2)/C(10) ! f12
 470       RETURN
 471       END
 472
 473       SUBROUTINE FSIIN(I_ITEST,I_ICH,I_IQS,I_ISI,I_I3C)
 474 C          SUBROUTINE FSIINI
 475 C---Note:
 476 C-- ICH= 0 (1) if the Coulomb interaction is absent (present);
 477 C-- ISPIN= JJ= 1,2,..,MSPIN denote increasing values of the pair
 478 C-- total spin S.
 479 C-- To calculate the CF of two particles (with masses m1, m2 and
 480 C-- charges C1, C2) the following information is required:
 481 C-- AM= twice the reduced mass= 2*m1*m2/(m1+m2) in GeV/c^2,
 482 C-- DM= (m1-m2)/(m1+m2), required if NS=2;
 483 C-- AC= Bohr radius= 2*137.036*0.1973/(C1*C2*AMH) in fm;
 484 C-- AC > 1.D9 if C1*C2= 0, AC < 0 if C1*C2 < 0;
 485 C-- MSPIN= MSPINH(LL)= number of the values of the total pair spin S;
 486 C-- FD= FDH(LL,JJ), RD= RDH(LL,JJ)= scattering length and effective
 487 C-- radius for each value of the total pair spin S, JJ= 1,..,MSPIN;     ;
 488 C-- the corresponding square well parameters EB= EBH(LL,JJ), RB=
 489 C-- RBH(LL,JJ) (required if NS=1) may be calculated by sear.f;
 490 C-- if the effective range approximation is not valid (as is the case,
 491 C-- e.g., for two-pion system) a code for calculation of the scattering
 492 C-- amplitude should be supplemented;
 493 C-- RHO= RHOH(LL,JJ), SF= SFH(LL,JJ), SE= SEH(LL) are spin factors;
 494 C-- RHO= the probability that the spins j1 and j2 of the two particles
 495 C-- will combine in a total spin S;
 496 C-- RHO= (2*S+1)/[(2j1+1)*(2j2+1)] for unpolarized particles;
 497 C-- RHO= (1-P1*P2)/4 and (3+P1*P2)/4 correspond to S=0 and 1 in the
 498 C-- case of spin-1/2 particles with polarizations P1 and P2;
 499 C-----------------------------------------------------------------------
 500       IMPLICIT REAL*8 (A-H,O-Z)
 501       COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
 502       COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
 503      1               X,Y,Z,T,RP,RPS
 504       COMMON/FSI_SPIN/RHO(10)
 505       COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
 506       COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
 507       COMMON/FSI_FD/FD(10),RD(10)
 508       COMMON/FSI_C/C(10),AM,AMS,DM
 509       COMMON/FSI_CONS/PI,PI2,SPI,DR,W
 510       COMPLEX*16 C
 511       COMMON/FSI_AA/AA
 512       COMMON/FSI_AAPI/AAPI(20,2)/FSI_AAND/AAND(20,4)
 513       COMMON/FSI_AAPIN/AAPIN(20,2)
 514       COMMON/FSI_AAPAP/AAPAPR(3,2),AAPAPI(3,2)
 515       COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10),
 516      1              SBKRB(10),SDKK(10)
 517       COMMON/FSI_FDH/FDH(30,10),RDH(30,10),EBH(30,10),RBH(30,10)
 518       COMMON/FSI_RHOH/RHOH(30,10)
 519       COMMON/FSI_AMCH/AM1H(30),AM2H(30),C1H(30),C2H(30),MSPINH(30)
 520 C============= declarations pour l'appel de READ_FILE()============
 521 c      CHARACTER*10 KEY
 522 c      CHARACTER*8  CH8
 523 c      INTEGER*4    INT4
 524 c      REAL*8       REAL8
 525 c      INTEGER*4    IERR
 526 C
 527 C--- mass of the first and second particle
 528       DATA AM1H/.93956563D0,.93827231D0,.93956563D0,3.72737978D0,
 529      C          .13957D0,.13498D0,.13957D0, .93956563D0, .93827231D0,
 530      C          4*.13957D0,4*.493677D0,
 531      C          2*1.87561339D0,2*2.80892165D0,2*.497672D0,
 532      C          1.87561339D0,3*.93827231D0,.93956563D0,
 533      C          1.115684D0,.93827231D0/
 534       DATA AM2H/.93956563D0,.93827231D0,.93827231D0,3.72737978D0,
 535      C          .13957D0,.13498D0,.13957D0, 2*1.87561339D0,
 536      C          2*.493677D0,2*.93827231D0,
 537      C          2*.493677D0,2*.93827231D0,
 538      C          1.87561339D0,3.72737978D0,2.80892165D0,3.72737978D0,
 539      C          2*.497672D0,2*2.80892165D0,3.72737978D0,
 540      C          3*1.115684D0,.93827231D0/
 541 c--------|---------|---------|---------|---------|---------|---------|----------
 542 C---  charge of the first and second particle
 543       DATA C1H/0.D0,1.D0,0.D0,2.D0, 1.D0,0.D0,1.D0,0.D0,1.D0,
 544      C         3*1.D0,-1.D0,3*1.D0,-1.D0,
 545      C         4*1.D0,2*0.D0,4*1.D0,2*0.D0, 1.D0/
 546       DATA C2H/0.D0,1.D0,1.D0,2.D0,-1.D0,0.D0,3*1.D0,
 547      C         -1.D0,3*1.D0,-1.D0,3*1.D0,
 548      C         1.D0,2.D0,1.D0,2.D0,2*0.D0,2*1.D0,2.D0,3*0.D0,-1.D0/
 549 C---MSPIN vs (LL)
 550       DATA MSPINH/3*2,4*1,2*2,8*1,3,1,2,1,2*1,2*2,1,3*2, 2/
 551 C---Spin factors <RHO vs (LL,ISPIN)
 552       DATA RHOH/3*.25D0, 4*1.D0, 2*.3333D0, 8*1.D0,
 553      1          .1111D0,1.D0,.25D0,1.D0,2*1.D0,
 554      1          .3333D0,.25D0,1.D0,3*.25D0, .25D0,
 555      2          3*.75D0, 4*0.D0, 2*.6667D0, 8*0.D0,
 556      2          .3333D0,.0D0,.75D0,.0D0,2*0.D0,
 557      2          .6667D0,.75D0,0.D0,3*.75D0, .75D0,
 558      3          17*.0D0,.5556D0,3*0.D0, 8*0.D0,1*0.D0,210*0.D0/
 559 C---Scattering length FD and effective radius RD in fm vs (LL,ISPIN)
 560       DATA FDH/
 561      1     17.0D0,7.8D0,23.7D0,2230.1218D0,.225D0,.081D0,-.063D0,
 562      1     -.65D0,-2.73D0,
 563      1     .137D0,-.071D0,-.148D0,.112D0,2*1.D-6,-.360D0,
 564      1     2*1.D-6,1.344D0,6*1.D-6,-5.628D0,2.18D0,2.40D0,
 565      1     2.81D0,              ! ND potential
 566 C     1     0.50D0,              ! NSC97e potential lam-lam
 567      1     1*0.001D0,
 568 C     c     2 -10.8D0,2*-5.4D0,4*0.D0,-6.35D0,-11.88D0,8*0.D0,9*0.D0,
 569      2     3*-5.4D0,4*0.D0,-6.35D0,-11.88D0,8*0.D0,9*0.D0,
 570      2     1.93D0,1.84D0,
 571      2     0.50D0,              ! triplet f0 lam-lam=singlet f0 ND
 572                                 ! not contributing in s-wave FSI approx.
 573      2     1*0.001D0,
 574      3     240*0.D0/
 575 c--------|---------|---------|---------|---------|---------|---------|----------
 576       DATA RDH/
 577      1     2.7D0,2.8D0,2.7D0,1.12139906D0,-44.36D0,64.0D0,784.9D0,
 578      1     477.9D0, 2.27D0, 9*0.D0,-69.973D0, 6*0.D0,3.529D0,
 579      1     3.19D0,3.15D0,
 580      1     2.95D0,              ! ND potential lam-lam
 581 C     1  10.6D0, ! NSC97e potential lam-lam
 582      1     1*0.D0,
 583      2     3*1.7D0,4*0.D0,2.0D0,2.63D0, 17*0.D0,3.35D0,3.37D0,
 584      2     2.95D0,              ! triplet d0 lam-lam=singlet d0 ND
 585                                 ! not contributing in s-wave approx.
 586      2     1*0.D0,
 587      3     240*0.D0/
 588 C---  Corresponding square well parameters RB (width in fm) and
 589 C--   EB =SQRT(-AM*U) (in GeV/c); U is the well height
 590       DATA RBH/2.545739D0,   2.779789D0, 2.585795D0, 5.023544D0,
 591      1     .124673D0, .3925180D0,.09D0, 2.D0, 4.058058D0, 17*0.D0,
 592      1     2.252623D0, 2.278575D0,
 593      1     2.234089D0,          ! ND potential lam-lam
 594 C     1  3.065796D0, ! NSC97e potential lam-lam
 595      1     1*0.001D0,
 596      2     3*2.003144D0,
 597      2     4*0.D0, 2.D0, 4.132163D0, 17*0.D0,
 598      2     2.272703D0, 2.256355D0,
 599      2     2.234089D0,          ! triplet potential lam-lam=singlet ND
 600                                 ! not contributing in s-wave FSI approx.
 601      2     1*0.001D0,
 602      3     240*0.D0/
 603       DATA EBH/.1149517D0,    .1046257D0,   .1148757D0, .1186010D0,
 604      1     .7947389D0,2.281208D0,8.7D0,.4D0,.1561219D0,17*0.D0,
 605      1     .1013293D0, .1020966D0,
 606      1     .1080476D0,          ! ND potential lam-lam
 607 C     1    .04115994D0, ! NSC97e potential lam-lam
 608      1     1*0.001D0,
 609      2     3*.1847221D0,
 610      2     4*0.D0, .4D0, .1150687D0, 17*0.D0,
 611      2     .09736083D0, .09708310D0,
 612      2     .1080476D0,          ! triplet potential lam-lam= singlet ND
 613                                 ! not contributing in s-wave FSI approx.
 614      2     1*0.001D0,
 615      3     240*0.D0/
 616
 617 C++----- add to be able to call several time-------
 618       integer ifirst
 619       data ifirst/0/
 620       ifirst=ifirst+1
 621
 622 C=======< constants >========================
 623       W=1/.1973D0    ! from fm to 1/GeV
 624       PI=4*DATAN(1.D0)
 625       PI2=2*PI
 626       SPI=DSQRT(PI)
 627       DR=180.D0/PI   ! from radian to degree
 628       AC1=1.D10
 629       AC2=1.D10
 630 C=======< condition de calculs >=============
 631       NUNIT=11 ! for IBM in Nantes
 632 C      NUNIT=4 ! for SUN in Prague
 633 C++      CALL readint4(NUNIT,'ITEST     ',ITEST)
 634 C++      CALL readint4(NUNIT,'NS        ',NS)
 635       ITEST=I_ITEST
 636       IF(ITEST.EQ.1)THEN
 637 C++      CALL readint4(NUNIT,'ICH       ',ICH)
 638 C++      CALL readint4(NUNIT,'IQS       ',IQS)
 639 C++      CALL readint4(NUNIT,'ISI       ',ISI)
 640 C++      CALL readint4(NUNIT,'I3C       ',I3C)
 641        ICH=I_ICH
 642        IQS=I_IQS
 643        ISI=I_ISI
 644        I3C=I_I3C
 645       ENDIF
 646 C============================================
 647       IF (IFIRST.LE.1) THEN
 648          DO 3 J1=1,30
 649             DO 3 J2=1,10
 650                FDH(J1,J2)=FDH(J1,J2)*W
 651                RDH(J1,J2)=RDH(J1,J2)*W
 652  3             RBH(J1,J2)=RBH(J1,J2)*W
 653 C     print *,"FD,RD,EB,RB: ",FDH(J1,J2),RDH(J1,J2),EBH(J1,J2),RBH(J1,J2)
 654
 655       ENDIF
 656 C===================================
 657       RETURN
 658       END
 659 C
 660       SUBROUTINE LLINI(lll,I_NS,I_ITEST)
 661 C===> Initialisation for a given LL value.
 662 C     ===========================================
 663       IMPLICIT REAL*8 (A-H,O-Z)
 664       COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
 665       COMMON/FSI_SPIN/RHO(10)
 666       COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
 667       COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
 668       COMMON/FSI_FD/FD(10),RD(10)
 669       COMMON/FSI_C/C(10),AM,AMS,DM
 670       COMMON/FSI_CONS/PI,PI2,SPI,DR,W
 671       COMPLEX*16 C
 672       COMMON/FSI_AA/AA
 673       COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
 674      1               X,Y,Z,T,RP,RPS
 675       COMMON/FSI_AAPI/AAPI(20,2)/FSI_AAND/AAND(20,4)
 676       COMMON/FSI_AAPIN/AAPIN(20,2)
 677       COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10),
 678      1              SBKRB(10),SDKK(10)
 679       COMMON/FSI_FDH/FDH(30,10),RDH(30,10),EBH(30,10),RBH(30,10)
 680       COMMON/FSI_RHOH/RHOH(30,10)
 681       COMMON/FSI_AMCH/AM1H(30),AM2H(30),C1H(30),C2H(30),MSPINH(30)
 682       COMMON/FSI_AAKK/AAKK(9)/FSI_AAPAP/AAPAPR(3,2),AAPAPI(3,2)
 683
 684 C++----- add to be able to call several time-------
 685       integer ifirst
 686       data ifirst/0/
 687       ifirst=ifirst+1
 688
 689 C===> LL - Initialisation ========================================
 690 C---- Setting particle masses and charges
 691       LL=lll
 692       NS=I_NS
 693       AM1=AM1H(LL)
 694       AM2=AM2H(LL)
 695       C1=C1H(LL)
 696       C2=C2H(LL)
 697 C      print *,"llini: ",AM1,AM2,C1,C2
 698 C---> Switches:
 699 C     ISI=1(0)  the strong interaction between the two particles ON (OFF)
 700 C     IQS=1(0)  the quantum statistics ON (OFF);
 701 C               should be OFF for nonidentical particles
 702 C     I3C=1(0)  the Coulomb interaction with the nucleus ON (OFF)
 703 C     I3S=1(0)  the strong interaction with the nucleus ON (OFF)
 704 C     ICH=1(0)  if C1*C2 is different from 0 (is equal to 0)
 705 C-    To switch off the Coulomb force between the two particles
 706 C     put ICH=0 and substitute the strong amplitude parameters by
 707 C     the ones not affected by Coulomb interaction
 708 C---- --------------------------------------------------------------------
 709       IF (I_ITEST.NE.1) THEN
 710          ICH=0
 711          IF(C1*C2.NE.0.D0) ICH=1
 712          IQS=0
 713          IF(C1+AM1.EQ.C2+AM2) IQS=1
 714          I3S=0                  ! only this option is available
 715          ISI=1
 716          I3C=1
 717       ENDIF
 718 C---> Calcul. twice the reduced mass (AM), the relative
 719 C     mass difference (DM) and the Bohr radius (AC)
 720       AM=2*AM1*AM2/(AM1+AM2)
 721       AMS=AM*AM
 722       DM=(AM1-AM2)/(AM1+AM2)
 723       AC=1.D10
 724       C12=C1*C2
 725       IF(C12.NE.0.D0)AC=2*137.036D0/(C12*AM)
 726 C---Setting spin factors
 727       MSPIN=MSPINH(LL)
 728 C      print *,"MSPIN: ",MSPIN
 729       MSP=MSPIN
 730       DO 91 ISPIN=1,10
 731  91      RHO(ISPIN)=RHOH(LL,ISPIN)
 732 C 91      print *,"RHO: ",ISPIN,RHO(ISPIN)
 733 C---> Integration limit AA in the spherical wave approximation
 734       AA=0.D0
 735 cc      IF(NS.EQ.2.OR.NS.EQ.4)AA=.5D0 !!in 1/GeV --> 0.1 fm
 736       IF(NS.EQ.2.OR.NS.EQ.4)AA=6.D0 !!in 1/GeV --> 1.2 fm
 737 C---> Setting scatt. length (FD), eff. radius (RD) and, if possible,
 738 C--   also the corresp. square well parameters (EB, RB)
 739       DO 55 JJ=1,MSP
 740          ISPIN=JJ
 741          FD(JJ)=FDH(LL,JJ)
 742          RD(JJ)=RDH(LL,JJ)
 743          EB(JJ)=EBH(LL,JJ)
 744          RB(JJ)=RBH(LL,JJ)
 745 C         print *,"FD,RD,EB,RB: ",FD(JJ),RD(JJ),EB(JJ),RB(JJ)
 746 C---Resets FD and RD for a nucleon-deuteron system (LL=8,9)
 747       IF(LL.EQ.8.OR.LL.EQ.9)THEN
 748        JH=LL-7+2*JJ-2
 749        FD(JJ)=AAND(1,JH)
 750        RD(JJ)=AAND(2,JH)-2*AAND(3,JH)/AAND(1,JH)
 751       ENDIF
 752 C---Resets FD and RD for a pion-pion system (LL=5,6,7)
 753       IF(LL.EQ.5.OR.LL.EQ.6.OR.LL.EQ.7)THEN
 754        IF(LL.EQ.7)FD(JJ)=AAPI(1,2)/AM
 755        IF(LL.EQ.5)FD(JJ)=(.6667D0*AAPI(1,1)+.3333D0*AAPI(1,2))/AM
 756        IF(LL.EQ.6)FD(JJ)=(.3333D0*AAPI(1,1)+.6667D0*AAPI(1,2))/AM
 757        AKS=0.D0
 758        DAKS=1.D-5
 759        AKSH=AKS+DAKS
 760        AKH=DSQRT(AKSH)
 761        GPI1H=GPIPI(AKSH,1)
 762        GPI2H=GPIPI(AKSH,2)
 763        H=1/FD(JJ)
 764        IF(LL.EQ.7)C(JJ)=1/DCMPLX(GPI2H,-AKH)
 765        IF(LL.EQ.5)
 766      + C(JJ)=.6667D0/DCMPLX(GPI1H,-AKH)+.3333D0/DCMPLX(GPI2H,-AKH)
 767        IF(LL.EQ.6)
 768      + C(JJ)=.3333D0/DCMPLX(GPI1H,-AKH)+.6667D0/DCMPLX(GPI2H,-AKH)
 769        HH=DREAL(1/C(JJ))
 770        RD(JJ)=2*(HH-H)/DAKS
 771       ENDIF
 772 C---Resets FD and RD for a pion-nucleon system (LL=12,13)
 773       IF(LL.EQ.12.OR.LL.EQ.13)THEN
 774        IF(LL.EQ.12)FD(JJ)=AAPIN(1,2)
 775        IF(LL.EQ.13)FD(JJ)=(.6667D0*AAPIN(1,1)+.3333D0*AAPIN(1,2))
 776        AKS=0.D0
 777        DAKS=1.D-5
 778        AKSH=AKS+DAKS
 779        AKH=DSQRT(AKSH)
 780        GPI1H=GPIN(AKSH,1)
 781        GPI2H=GPIN(AKSH,2)
 782        H=1/FD(JJ)
 783        IF(LL.EQ.12)C(JJ)=1/DCMPLX(GPI2H,-AKH)
 784        IF(LL.EQ.13)
 785      + C(JJ)=.6667D0/DCMPLX(GPI1H,-AKH)+.3333D0/DCMPLX(GPI2H,-AKH)
 786        HH=DREAL(1/C(JJ))
 787        RD(JJ)=2*(HH-H)/DAKS
 788       ENDIF
 789 C---fm to 1/GeV for pp-bar system
 790       IF(LL.EQ.30)THEN
 791          IF(IFIRST.LE.1)THEN
 792             DO 4 I3=1,3
 793                AAPAPR(I3,JJ)=AAPAPR(I3,JJ)*W
 794  4             AAPAPI(I3,JJ)=AAPAPI(I3,JJ)*W
 795 C     4       print *,"AAPAPR,AAPAPI: ",AAPAPR(I3,JJ),AAPAPI(I3,JJ)
 796 C---  Calculates complex elements M11=M22=C(6), M12=M21=C(7) for I=0
 797 C---  at k*=0                  M11=M22=C(8), M12=M21=C(9) for I=1
 798           C(7+(JJ-1)*2)=2*DCMPLX(AAPAPR(1,JJ),AAPAPI(1,JJ))*
 799      *         DCMPLX(AAPAPR(2,JJ),AAPAPI(2,JJ)) ! 2a_0Sa_1S
 800           C(6+(JJ-1)*2)=DCMPLX(AAPAPR(1,JJ)+AAPAPR(2,JJ),
 801      ,         AAPAPI(1,JJ)+AAPAPI(2,JJ))/
 802      /         C(7+(JJ-1)*2)    ! M11=M22
 803           C(7+(JJ-1)*2)=-DCMPLX(AAPAPR(1,JJ)-AAPAPR(2,JJ),
 804      ,         AAPAPI(1,JJ)-AAPAPI(2,JJ))/
 805      /         C(7+(JJ-1)*2)    ! M12=M21
 806         ENDIF
 807       ENDIF
 808 C---Calculation continues for any system (any LL)
 809  55   CONTINUE
 810       RETURN
 811       END
 812 C=======================================================
 813 C
 814
 815 c++  This routine is used to init mass and charge of the nucleus.
 816
 817       SUBROUTINE FSINUCL(R_AMN,R_CN)
 818
 819       IMPLICIT REAL*8 (A-H,O-Z)
 820       COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
 821
 822       AMN=R_AMN
 823       CN=R_CN
 824
 825       RETURN
 826       END
 827
 828 C======================================================
 829 C
 830
 831       SUBROUTINE FSIMOMENTUM(PP1,PP2)
 832
 833       IMPLICIT REAL*8 (A-H,O-Z)
 834       COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1,  ! particle momenta in NRF
 835      1               P2X,P2Y,P2Z,E2,P2
 836
 837
 838       REAL*8 PP1(3),PP2(3)
 839 c      Print *,"momentum",pp1,pp2
 840       P1X=PP1(1)
 841       P1Y=PP1(2)
 842       P1Z=PP1(3)
 843       P2X=PP2(1)
 844       P2Y=PP2(2)
 845       P2Z=PP2(3)
 846       RETURN
 847       END
 848
 849
 850
 851 C======================================================
 852 C
 853
 854       SUBROUTINE FSIPOSITION(XT1,XT2)
 855
 856       IMPLICIT REAL*8 (A-H,O-Z)
 857       COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, !4-coord. of emis. points in NRF
 858      1                X2,Y2,Z2,T2,R2
 859
 860       REAL*8 XT1(4),XT2(4)
 861 clc      print *,'fsi',xt1,xt2
 862       X1=XT1(1)
 863       Y1=XT1(2)
 864       Z1=XT1(3)
 865       T1=XT1(4)
 866       X2=XT2(1)
 867       Y2=XT2(2)
 868       Z2=XT2(3)
 869       T2=XT2(4)
 870       RETURN
 871       END
 872
 873
 874 C======================================================
 875 C======================================================
 876 C
 877       subroutine BoostToPrf()
 878       IMPLICIT REAL*8 (A-H,O-Z)
 879       COMMON/FSI_CVK/V,CVK
 880       COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1,  !part. momenta in NRF
 881      1               P2X,P2Y,P2Z,E2,P2
 882       COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
 883      1               X,Y,Z,T,RP,RPS
 884       COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emis. points in NRF
 885      1                X2,Y2,Z2,T2,R2
 886       COMMON/FSI_P12/P12X,P12Y,P12Z,E12,P12,AM12,EPM
 887
 888       XS=X1-X2
 889       YS=Y1-Y2
 890       ZS=Z1-Z2
 891       TS=T1-T2
 892       RS12=XS*P12X+YS*P12Y+ZS*P12Z
 893       H1=(RS12/EPM-TS)/AM12
 894       X=XS+P12X*H1
 895       Y=YS+P12Y*H1
 896       Z=ZS+P12Z*H1
 897       T=(E12*TS-RS12)/AM12
 898       RPS=X*X+Y*Y+Z*Z
 899       RP=DSQRT(RPS)
 900 CW      WRITE(6,38)'RP ',RP,'X ',X,Y,Z,T
 901 38    FORMAT(A7,E11.4,A7,4E11.4)
 902
 903       CVK=(P12X*PPX+P12Y*PPY+P12Z*PPZ)/(P12*AK)
 904       V=P12/E12
 905       return
 906       end
 907
 908       SUBROUTINE FSIWF(WEI)
 909 C==> Prepares necessary quantities and call VZ(WEI) to calculate
 910 C    the weight due to FSI
 911       IMPLICIT REAL*8 (A-H,O-Z)
 912       COMMON/FSI_CVK/V,CVK
 913       COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1,  !part. momenta in NRF
 914      1               P2X,P2Y,P2Z,E2,P2
 915       COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
 916      1               X,Y,Z,T,RP,RPS
 917       COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emis. points in NRF
 918      1                X2,Y2,Z2,T2,R2
 919       COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
 920       COMMON/FSI_SPIN/RHO(10)
 921       COMMON/FSI_BP/B,P
 922       COMMON/FSI_ETA/ETA
 923       COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
 924       COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10),
 925      1              SBKRB(10),SDKK(10)
 926       COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
 927       COMMON/FSI_RR/F(10)
 928       COMMON/FSI_FD/FD(10),RD(10)
 929       COMMON/FSI_C/C(10),AM,AMS,DM
 930       COMPLEX*16 C,F
 931       COMMON/FSI_AA/AA
 932       COMMON/FSI_SHH/SH,CHH
 933       COMMON/FSI_AAPI/AAPI(20,2)/FSI_AAND/AAND(20,4)
 934       COMMON/FSI_AAPIN/AAPIN(20,2)
 935       COMMON/FSI_P12/P12X,P12Y,P12Z,E12,P12,AM12,EPM
 936       COMMON/FSI_2CHA/AK2,AK2S,AAK,HCP2,AMU2_AMU1 ! k* (kappa) for 2-nd channel
 937 C==>calculating relative 4-coordinates of the particles in PRF
 938 C-  {T,X,Y,Z} from the relative coordinates {TS,XS,YS,ZS} in NRF
 939       XS=X1-X2
 940       YS=Y1-Y2
 941       ZS=Z1-Z2
 942       TS=T1-T2
 943       RS12=XS*P12X+YS*P12Y+ZS*P12Z
 944       H1=(RS12/EPM-TS)/AM12
 945       X=XS+P12X*H1
 946       Y=YS+P12Y*H1
 947       Z=ZS+P12Z*H1
 948       T=(E12*TS-RS12)/AM12
 949       RPS=X*X+Y*Y+Z*Z
 950       RP=DSQRT(RPS)
 951 CW      WRITE(6,38)'RP ',RP,'X ',X,Y,Z,T
 952 38    FORMAT(A7,E11.4,A7,4E11.4)
 953
 954       CVK=(P12X*PPX+P12Y*PPY+P12Z*PPZ)/(P12*AK)
 955       V=P12/E12
 956
 957 C      ACH=1.D0
 958       IF(ICH.EQ.0)GOTO 21
 959       XH=AC*AK
 960       ACH=ACP(XH)
 961       ACHR=DSQRT(ACH)
 962       ETA=0.D0
 963       IF(XH.NE.0.D0)ETA=1/XH
 964 C---HCP, HPR needed (e.g. in GST) if ICH=1
 965       HCP=HC(XH)
 966       HPR=HCP+.1544313298D0
 967 C      AAK=ACH*AK    !
 968 C      HCP2=2*HCP/AC ! needed to calculate C(JJ) for charged particles
 969   21  CONTINUE
 970       MSP=MSPIN
 971       DO 30 JJ=1,MSP
 972       ISPIN=JJ
 973       IF(NS.NE.1)GOTO22
 974 C---Calc. quantities for the square well potential;
 975 C-- for LL=6-26 the square well potential is not possible or available
 976       IF(LL.EQ.4)GOTO 22
 977       BK(JJ)=DSQRT(EB(JJ)**2+AKS)
 978       XRA=2*RB(JJ)/AC
 979       HRA=BK(JJ)*RB(JJ)
 980       CALL SEQ(XRA,HRA)
 981       SBKRB(JJ)=HRA*B
 982       HRA=AK*RB(JJ)
 983       CALL GST(XRA,HRA)
 984       SDK(JJ)=SH
 985       CDK(JJ)=CHH
 986       SDKK(JJ)=RB(JJ)
 987       IF(AK.NE.0.D0)SDKK(JJ)=SH/AK
 988       IF(ICH.EQ.1)SDK(JJ)=ACH*SDK(JJ)
 989   22  CONTINUE
 990 C-----------------------------------------------------------------------
 991 C---Calc. the strong s-wave scattering amplitude = C(JJ)
 992 C-- divided by Coulomb penetration factor squared (if ICH=1)
 993       IF(NS.NE.1)GOTO 230
 994       IF(LL.NE.4)GOTO 230 ! SW scat. amplitude used for alfa-alfa only
 995       GAK=G(AK)
 996       AKACH=AK
 997       IF(ICH.EQ.1)AKACH=AK*ACH
 998       C(JJ)=1/DCMPLX(GAK,-AKACH) ! amplitude for the SW-potential
 999       GOTO 30
1000  230  IF(LL.EQ.5.OR.LL.EQ.6.OR.LL.EQ.7)GOTO20    ! pipi
1001       IF(LL.EQ.12.OR.LL.EQ.13)GOTO20             ! piN
1002       IF(LL.EQ.8.OR.LL.EQ.9.OR.LL.EQ.18)GOTO20   ! Nd, dd
1003       IF(LL.EQ.14.OR.LL.EQ.17.OR.LL.EQ.23)GOTO27 ! K+K-, K-p, K0K0-b
1004         IF(LL.EQ.30)GOTO 28                        ! pp-bar
1005       A1=RD(JJ)*FD(JJ)*AKS
1006       A2=1+.5D0*A1
1007       IF(ICH.EQ.1)A2=A2-2*HCP*FD(JJ)/AC
1008       AKF=AK*FD(JJ)
1009       IF(ICH.EQ.1)AKF=AKF*ACH
1010       C(JJ)=FD(JJ)/DCMPLX(A2,-AKF)
1011       GOTO30
1012  20   CONTINUE
1013 C---Calc. scatt. ampl. C(JJ) for pipi, piN and Nd, dd
1014       JH=LL-7+2*JJ-2
1015       IF(LL.EQ.8.OR.LL.EQ.9)GPI2=GND(AKS,JH)
1016       IF(LL.EQ.18)GPI2=GDD(AKS,JJ)
1017       IF(LL.EQ.5.OR.LL.EQ.6.OR.LL.EQ.7)GPI2=GPIPI(AKS,2)
1018       IF(LL.EQ.12.OR.LL.EQ.13)GPI2=GPIN(AKS,2)
1019       C(JJ)=1.D0/DCMPLX(GPI2,-AK) !pi+pi+, nd, pd, pi+p, dd
1020       IF(LL.NE.5.AND.LL.NE.6.AND.LL.NE.13)GOTO27
1021       IF(LL.EQ.5.OR.LL.EQ.6)GPI1=GPIPI(AKS,1)
1022       IF(LL.EQ.13)GPI1=GPIN(AKS,1)
1023       IF(LL.EQ.5.OR.LL.EQ.13)
1024      c           C(JJ)=.6667D0/DCMPLX(GPI1,-AK)+.3333D0*C(JJ) !pi+pi-,pi-p
1025       IF(LL.EQ.6)C(JJ)=.3333D0/DCMPLX(GPI1,-AK)+.6667D0*C(JJ) !pi0pi0
1026  27   CONTINUE
1027 C---Calc. K+K-, K0K0-b or K-p s-wave scatt. ampl.
1028       IF(LL.EQ.14.OR.LL.EQ.23)CALL CKKB
1029 c      IF(LL.EQ.17)C(JJ)=DCMPLX(-3.29D0,3.55D0)     ! Martin'76 (K-p)
1030 c      IF(LL.EQ.17)C(JJ)=DCMPLX(3.29D0,3.55D0)     ! Martin'76 (K-p) WRONG SIGN!!!
1031       IF(LL.EQ.17)C(JJ)=DCMPLX(-2.585D0,4.156D0)   ! Borasoy'04 (K-p)
1032 c      IF(LL.EQ.17)C(JJ)=DCMPLX(-3.371D0,3.244D0)   ! Martin'81 (K-p)
1033 C---Calc. pi+pi-, pi+pi+, pd, pi+p, pi-p, K+K- or K-p s-wave scatt. ampl.
1034 C-- divided by Coulomb penetration factor squared (if ICH=1)
1035       IF(ICH.EQ.0)GOTO 30
1036       AAK=ACH*AK    !
1037       HCP2=2*HCP/AC ! needed to calculate C(JJ) for charged particles
1038       C(JJ)=1/(1/C(JJ)-HCP2+DCMPLX(0.D0,AK-AAK))
1039       GOTO 30
1040  28   CONTINUE
1041 C---Calc. pp-bar s-wave scatt. ampl.
1042       CALL CPAP
1043  30   CONTINUE
1044 C***********************************************************************
1045       CALL VZ(WEI)
1046       RETURN
1047       END
1048       SUBROUTINE VZ(WEI)
1049 C==> Calculates the weight WEI due to FSI
1050       IMPLICIT REAL*8 (A-H,O-Z)
1051       COMMON/FSI_JR/JRAT
1052       COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
1053      1               X,Y,Z,T,RP,RPS
1054       COMMON/FSI_SPIN/RHO(10)
1055       COMMON/FSI_ETA/ETA
1056       COMMON/FSI_AA/AA
1057       COMMON/FSI_FFF/F12,F21
1058       COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
1059       COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
1060       COMMON/FSI_FD/FD(10),RD(10)
1061       COMMON/FSI_RR/F(10)
1062       COMMON/FSI_C/C(10),AM,AMS,DM
1063       COMMON/FSI_COULPH/EIDC
1064       COMPLEX*16 F,C,G,PSI12,PSI21
1065       COMPLEX*16 F12,F21
1066       COMPLEX*16 EIDC
1067       COMPLEX*8 Z8,CGAMMA
1068       COMMON/FSI_FFPN/FF12,FF21
1069       COMPLEX*16 FF12,FF21
1070       COMMON/FSI_2CHA/AK2,AK2S,AAK,HCP2,AMU2_AMU1 ! k* (kappa) for 2-nd channel
1071       WEI=0.D0
1072       IF(JRAT.EQ.1)GOTO 11
1073       RHOS=AK*RP
1074       HS=X*PPX+Y*PPY+Z*PPZ
1075       IF(RHOS.LT.15.D0.AND.RHOS+DABS(HS).LT.20.D0)GOTO 2
1076 C---Calc. EIDC=exp(i*Coul.Ph.);
1077 C-- used in calc. of hypergeom. f-s in SEQA, FAS at k*R > 15, 20
1078       Z8=CMPLX(1.,SNGL(ETA))
1079       Z8=CGAMMA(Z8)
1080       EIDC=Z8/CABS(Z8)
1081  2    CALL FF(RHOS,HS)
1082  11   MSP=MSPIN
1083       IF(ISI.EQ.0)GOTO 4  ! the strong interaction ON (OFF) if ISI=1(0)
1084       IF(RP.LT.AA)GOTO 4
1085       IF(JRAT.NE.1) CALL FIRT
1086       IF(IQS.EQ.0)GOTO 5  ! the quantum statistics ON (OFF) if IQS=1(0)
1087       JSIGN=-1
1088       DO 1 JJ=1,MSP
1089       JSIGN=-JSIGN
1090       G=F(JJ)*C(JJ)
1091       IF(ICH.EQ.1)G=G*ACHR
1092       PSI12=FF12*(F12+G)
1093       PSI21=FF21*(F21+G)
1094       G=PSI12+JSIGN*PSI21
1095  1    WEI=WEI+RHO(JJ)*(DREAL(G)**2+DIMAG(G)**2)
1096       GOTO 8
1097  5    DO 6 JJ=1,MSP
1098       G=F(JJ)*C(JJ)
1099       IF(ICH.EQ.1)G=G*ACHR
1100 CW      WRITE(6,38)'JJ ',JJ,'F ',DREAL(F(JJ)),DIMAG(F(JJ))
1101 CW      WRITE(6,38)'JJ ',JJ,'C ',DREAL(C(JJ)),DIMAG(C(JJ))
1102 CW      WRITE(6,38)'JJ ',JJ,'G ',DREAL(G),DIMAG(G)
1103 CW      WRITE(6,38)'JJ ',JJ,'F12+G ',DREAL(F12+G),DIMAG(F12+G)
1104 CW      WRITE(6,38)'JJ ',JJ,'F21+G ',DREAL(F21+G),DIMAG(F21+G)
1105 38    FORMAT(A7,I3,A7,2E11.4)
1106       PSI12=FF12*(F12+G)
1107  6    WEI=WEI+RHO(JJ)*(DREAL(PSI12)**2+DIMAG(PSI12)**2)
1108 c--- Account for nn-bar->pp-bar channel ---------------------------
1109       IF(LL.EQ.30)THEN
1110       DO 61 JJ=1,MSP
1111         HH=RHO(JJ)*(DREAL(C(JJ+2))**2+DIMAG(C(JJ+2))**2)*
1112      *  AMU2_AMU1*ACH/RPS
1113         IF(AK2S.LT.0)HH=HH*DEXP(-2*RP*AK2)
1114  61   WEI=WEI+HH
1115         ENDIF
1116 c------------------------------------------------------------------
1117       RETURN
1118  4    PSI12=FF12*F12
1119       IF(IQS.EQ.0)GOTO 50 ! the quantum statistics ON (OFF) if IQS=1(0)
1120       PSI21=FF21*F21
1121       JSIGN=-1
1122       DO 3 JJ=1,MSP
1123       JSIGN=-JSIGN
1124       G=PSI12+JSIGN*PSI21
1125  3    WEI=WEI+RHO(JJ)*(DREAL(G)**2+DIMAG(G)**2)
1126       GOTO 8
1127  50   WEI=DREAL(PSI12)**2+DIMAG(PSI12)**2
1128       RETURN
1129  8    WEI=WEI/2
1130       RETURN
1131       END
1132       SUBROUTINE FIRT
1133 C---CALC. THE F(JJ)
1134 C-- F(JJ)*C(JJ)= DEVIATION OF THE BETHE-SALPETER AMPL. FROM PLANE WAVE
1135       IMPLICIT REAL*8 (A-H,O-Z)
1136       COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
1137      1               X,Y,Z,T,RP,RPS
1138       COMMON/FSI_SHH/SH,CHH
1139       COMMON/FSI_BP/B,P
1140       COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
1141       COMMON/FSI_C/C(10),AM,AMS,DM
1142       COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10),
1143      1              SBKRB(10),SDKK(10)
1144       COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
1145       COMMON/FSI_RR/F(10)
1146       EQUIVALENCE(RSS,RP),(TSS,T)
1147       COMPLEX*16 F,C,CH1
1148       MSP=MSPIN
1149       DO 10 JJ=1,MSP
1150       IF(JJ.GT.1)GOTO 3
1151       XRA=2*RSS/AC
1152       IF(AK.NE.0.D0)GOTO2
1153       SHK=1.D0
1154       SH=.0D0
1155       SHH=SH
1156       CHH=1/RSS
1157       GOTO3
1158   2   H=AK*RSS
1159       CALL GST(XRA,H)
1160       SH=SH/RSS
1161       CHH=CHH/RSS
1162       SHH=SH
1163       IF(ICH.EQ.1) SHH=ACH*SH
1164   3   IF(NS.EQ.2)GOTO1
1165 C---F= ASYMPTOTIC FORMULA (T= 0 APPROX.); NS= 4
1166   6   F(JJ)=DCMPLX(CHH,SHH)
1167       IF(NS.NE.1)GOTO 10
1168 C---F INSIDE THE SQUARE-WELL (T= 0 APPROX.); NS= 1
1169       IF(RSS.GE.RB(JJ)) GOTO 10
1170       IF(AK.NE.0.D0.AND.JJ.EQ.1)SHK=SH/AK
1171       H=BK(JJ)*RSS
1172       CALL GST(XRA,H)
1173       SKR=B*BK(JJ)
1174       F(JJ)=DCMPLX(CDK(JJ),SDK(JJ))*SKR
1175       CH1=(SDKK(JJ)*SKR-SHK*SBKRB(JJ))/C(JJ)
1176       F(JJ)=(F(JJ)+CH1)/SBKRB(JJ)
1177       GOTO 10
1178   1   CONTINUE
1179 C---F= ASYMPTOTIC FORMULA (T= 0 NOT REQUIRED); NS= 2
1180       IF(JJ.GT.1)GOTO 8
1181       IF(TSS.EQ.0.D0)GOTO6
1182       TSSA=DABS(TSS)
1183       IF(DM.NE.0.D0)GOTO 11
1184       H=AM*.5D0/TSSA
1185       IF(AK.NE.0.D0)GOTO4
1186       HM=H*RPS
1187       IF(HM.GE.3.D15)GOTO6
1188       FS1=DFRSIN(HM)
1189       FC1=DFRCOS(HM)
1190       FC2=FC1
1191       FS2=FS1
1192       GOTO5
1193   4   CONTINUE
1194       H1=AK*TSSA/AM
1195       HM=H*(RSS-H1)**2
1196       HP=H*(RSS+H1)**2
1197       IF(HP.GE.3.D15)GOTO6
1198       FS1=DFRSIN(HM)
1199       FC1=DFRCOS(HM)
1200       FS2=DFRSIN(HP)
1201       FC2=DFRCOS(HP)
1202       GOTO 5
1203   11  CONTINUE
1204       FS1=0.D0
1205       FS2=0.D0
1206       FC1=0.D0
1207       FC2=0.D0
1208       DO 13 I=1,2
1209       IF(I.EQ.1)TSSH=TSSA*(1+DM)
1210       IF(I.EQ.2)TSSH=TSSA*(1-DM)
1211       H=AM*.5D0/TSSH
1212       IF(AK.NE.0.D0)GOTO 12
1213       HM=H*RPS
1214       IF(HM.GE.3.D15)GOTO6
1215       FS1=FS1+DFRSIN(HM)/2
1216       FC1=FC1+DFRCOS(HM)/2
1217       IF(I.EQ.1)GOTO 13
1218       FC2=FC1
1219       FS2=FS1
1220       GOTO 13
1221   12  CONTINUE
1222       H1=AK*TSSH/AM
1223       HM=H*(RSS-H1)**2
1224       HP=H*(RSS+H1)**2
1225       IF(HP.GE.3.D15)GOTO6
1226       FS1=FS1+DFRSIN(HM)/2
1227       FC1=FC1+DFRCOS(HM)/2
1228       FS2=FS2+DFRSIN(HP)/2
1229       FC2=FC2+DFRCOS(HP)/2
1230   13  CONTINUE
1231   5   C12=FC1+FS2
1232       S12=FS1+FC2
1233       A12=FS1-FC2
1234       A21=FS2-FC1
1235       A2=.5D0*(CHH*(A12+A21)+SH*(A12-A21))+SHH
1236       A1=.5D0*(CHH*(C12+S12)+SH*(C12-S12))
1237       F(JJ)=.3989422D0*DCMPLX(A1,A2)
1238       GOTO 10
1239   8   F(JJ)=F(1)
1240  10   CONTINUE
1241       RETURN
1242       END
1243       FUNCTION EXF(X)
1244       IMPLICIT REAL*8 (A-H,O-Z)
1245       IF(X.LT.-15.D0) GO TO 1
1246       EXF=DEXP(X)
1247       RETURN
1248   1   EXF=.0D0
1249       RETURN
1250       END
1251       SUBROUTINE SEQ(X,H)
1252 C---CALC. FUNCTIONS B, P (EQS. (17) OF G-K-L-L);
1253 C-- NEEDED TO CALC. THE CONFLUENT HYPERGEOMETRIC FUNCTION GST.
1254       IMPLICIT REAL*8 (A-H,O-Z)
1255       COMMON/FSI_BP/B,P
1256       DIMENSION BH(3),PH(3)
1257       DATA ERR/1.D-7/
1258       BH(1)=1.D0
1259       PH(1)=1.D0
1260       PH(2)=.0D0
1261       BH(2)=.5D0*X
1262       B=1+BH(2)
1263       P=1.D0
1264       HS=H*H
1265       J=0
1266   2   J=J+1
1267       BH(3)=(X*BH(2)-HS*BH(1))/((J+1)*(J+2))
1268       PH(3)=(X*PH(2)-HS*PH(1)-(2*J+1)*X*BH(2))/(J*(J+1))
1269       B=B+BH(3)
1270       P=P+PH(3)
1271       Z=DABS(BH(2))+DABS(BH(3))+DABS(PH(2))+DABS(PH(3))
1272       IF(Z.LT.ERR)RETURN
1273       BH(1)=BH(2)
1274       BH(2)=BH(3)
1275       PH(1)=PH(2)
1276       PH(2)=PH(3)
1277       GOTO 2
1278       END
1279       SUBROUTINE SEQA(X,H)
1280 C---CALC. FUNCTIONS CHH=REAL(GST), SH=IMAG(GST)/ACH, B=SH/H
1281 C-- IN THE ASYMPTOTIC REGION H=K*R >> 1.
1282       IMPLICIT REAL*8 (A-H,O-Z)
1283       COMMON/FSI_BP/B,P
1284       COMMON/FSI_SHH/SH,CHH
1285       COMMON/FSI_ETA/ETA
1286       COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
1287       COMMON/FSI_COULPH/EIDC
1288       COMPLEX*16 EIDC,GST
1289       ARG=H-ETA*DLOG(2*H)
1290       GST=DCMPLX(DCOS(ARG),DSIN(ARG))
1291       GST=ACHR*EIDC*GST
1292       CHH=DREAL(GST)
1293       SH=DIMAG(GST)/ACH
1294       B=SH/H
1295       RETURN
1296       END
1297       SUBROUTINE FF(RHO,H)
1298 C---Calc. F12, F21;
1299 C-- F12= FF0* plane wave,  FF0=F*ACHR,
1300 C---F is the confluent hypergeometric function,
1301 C-- ACHR=sqrt(ACH), where ACH is the Coulomb factor
1302       IMPLICIT REAL*8 (A-H,O-Z)
1303       COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
1304       COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
1305       COMMON/FSI_ETA/ETA
1306       COMMON/FSI_FFF/F12,F21
1307       COMPLEX*16 FF0,F12,F21
1308       C=DCOS(H)
1309       S=DSIN(H)
1310       F12=DCMPLX(C,-S)
1311       F21=DCMPLX(C,S)
1312       IF(ICH.EQ.0)RETURN
1313       RHOP=RHO+H
1314       RHOM=RHO-H
1315       F12=FF0(RHO,H)*F12
1316       F21=FF0(RHO,-H)*F21
1317       RETURN
1318       END
1319       FUNCTION FAS(RKS)
1320 C-- FAS=F*ACHR
1321 C---F is the confluent hypergeometric function at k*r >> 1
1322 C-- ACHR=sqrt(ACH), where ACH is the Coulomb factor
1323       IMPLICIT REAL*8 (A-H,O-Z)
1324       COMPLEX*16 FAS,EIDC,ZZ1
1325       COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
1326       COMMON/FSI_ETA/ETA
1327       COMMON/FSI_COULPH/EIDC
1328       D1=DLOG(RKS)*ETA
1329       D2=ETA*ETA/RKS
1330       ZZ1=DCMPLX(DCOS(D1),DSIN(D1))/EIDC
1331       FAS=DCMPLX(1.D0,-D2)*ZZ1
1332       FAS=FAS-DCMPLX(DCOS(RKS),DSIN(RKS))*ETA/RKS/ZZ1
1333       RETURN
1334       END
1335       FUNCTION FF0(RHO,H)
1336 C-- FF0=F*ACHR
1337 C-- F is the confluent hypergeometric function
1338 C-- (Eq. (15) of G-K-L-L), F= 1 at r* << AC
1339 C-- ACHR=sqrt(ACH), where ACH is the Coulomb factor
1340       IMPLICIT REAL*8 (A-H,O-Z)
1341       COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
1342       COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
1343       COMMON/FSI_ETA/ETA
1344       COMPLEX*16 ALF,ALF1,Z,S,A,FF0,FAS
1345       DATA ERR/1.D-5/
1346       S=DCMPLX(1.D0,0.D0)
1347       FF0=S
1348       RHOP=RHO+H
1349 CC      GOTO 5 ! rejects the approx. calcul. of hyperg. f-ion F
1350       IF(RHOP.LT.20.D0)GOTO5
1351       FF0=FAS(RHOP) ! approx. calc.
1352       RETURN
1353   5   ALF=DCMPLX(.0D0,-ETA)
1354       ALF1=ALF-1
1355       Z=DCMPLX(.0D0,RHOP)
1356       J=0
1357   3   J=J+1
1358       A=(ALF1+J)/(J*J)
1359       S=S*A*Z
1360       FF0=FF0+S
1361       ZR=DABS(DREAL(S))
1362       ZI=DABS(DIMAG(S))
1363       IF((ZR+ZI).GT.ERR)GOTO3
1364       FF0=FF0*ACHR
1365       RETURN
1366       END
1367       FUNCTION HC(XA)
1368 C---HC = h-function of Landau-Lifshitz: h(x)=Re[psi(1-i/x)]+ln(x)
1369 C-- psi(x) is the digamma function (the logarithmic derivative of
1370 C-- the gamma function)
1371       IMPLICIT REAL*8 (A-H,O-Z)
1372       DIMENSION BN(15)
1373       DATA BN/.8333333333D-1,.8333333333D-2,.396825396825D-2,
1374      1        .4166666667D-2,.7575757576D-2,.2109279609D-1,
1375      2        .8333333333D-1,.4432598039D0 ,.305395433D1,
1376      3        .2645621212D2, .2814601449D3, .3607510546D4,
1377      4        .5482758333D5, .9749368235D6, .200526958D8/
1378       X=DABS(XA)
1379       IF(X.LT..33D0) GOTO 1
1380 CC      IF(X.GE.3.5D0) GO TO 2
1381       S=.0D0
1382       N=0
1383    3  N=N+1
1384       DS=1.D0/N/((N*X)**2+1)
1385       S=S+DS
1386       IF(DS.GT.0.1D-12) GOTO 3
1387 C---Provides 7 digit accuracy
1388       HC=S-.5772156649D0+DLOG(X)
1389       RETURN
1390 CC   2  HC=1.2D0/X**2+DLOG(X)-.5772156649 D0
1391 CC      RETURN
1392    1  X2=X*X
1393       XP=X2
1394       HC=0.D0
1395       IMA=9
1396       IF(X.LT.0.1D0)IMA=3
1397       DO 4 I=1,IMA
1398       HC=HC+XP*BN(I)
1399    4  XP=X2*XP
1400       RETURN
1401       END
1402       FUNCTION ACP(X)
1403 C--- ACP = COULOMB PENETRATION FACTOR
1404       IMPLICIT REAL*8 (A-H,O-Z)
1405       IF(X.LT.0.05D0.AND.X.GE.0.D0) GO TO 1
1406       Y=6.2831853D0/X
1407       ACP=Y/(EXF(Y)-1)
1408       RETURN
1409    1  ACP=1.D-6
1410       RETURN
1411       END
1412       SUBROUTINE GST(X,H)
1413 C---CALC. THE CONFL. HYPERGEOM. F-N = CHH+i*SH
1414 C-- AND THE COULOMB F-S B, P (CALLS SEQ OR SEQA).
1415       IMPLICIT REAL*8 (A-H,O-Z)
1416       COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
1417       COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
1418       COMMON/FSI_SHH/SH,CHH
1419       COMMON/FSI_BP/B,P
1420   1   IF(ICH.EQ.1)GOTO 2
1421   3   SH=DSIN(H)
1422       CHH=DCOS(H)
1423       B=1.D0
1424       IF(H.NE.0.D0)B=SH/H
1425       P=CHH
1426       RETURN
1427   2   CONTINUE
1428       IF(H.GT.15.D0)GOTO4 ! comment out if you want to reject
1429                    ! the approximate calculation of hyperg. f-ion G
1430       CALL SEQ(X,H) ! exact calculation
1431       SH=H*B
1432       CHH=P+B*X*(DLOG(DABS(X))+HPR)
1433       RETURN
1434   4   CALL SEQA(X,H)
1435       RETURN
1436       END
1437       FUNCTION FF1(RHO,H)
1438 C---FF1=FF0; used for particle-nucleus system
1439 C-- FF0=F12*ACHR
1440 C-- F12 is the confluent hypergeometric function
1441 C-- (Eq. (15) of G-K-L-L), F12= 1 at r* << AC
1442 C-- ACHR=sqrt(ACH), where ACH is the Coulomb factor
1443       IMPLICIT REAL*8 (A-H,O-Z)
1444       COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
1445       COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
1446       COMMON/FSI_ETA/ETA
1447       COMMON/FSI_COULPH/EIDC
1448       COMMON/FSI_ICH1/ICH1
1449       COMPLEX*16 FF0,FF1
1450       COMPLEX*16 EIDC
1451       COMPLEX*8 Z8,CGAMMA
1452       FF1=DCMPLX(1.D0,0.D0)
1453       IF(ICH1.EQ.0)GOTO 2
1454       IF(RHO.LT.15.D0.AND.RHO+H.LT.20.D0)GOTO 2
1455 C---Calc. EIDC=exp(i*Coul.Ph.);
1456 C-- used in calc. of hypergeom. f-s in SEQA, FAS at k*R > 15, 20
1457       Z8=CMPLX(1.,SNGL(ETA))
1458       Z8=CGAMMA(Z8)
1459       EIDC=Z8/CABS(Z8)
1460  2    FF1=FF0(RHO,H)
1461       RETURN
1462       END
1463
1464       FUNCTION G(AK)
1465 C---Used to calculate SW scattering amplitude for alpa-alpha system
1466 C-- and for sear.f (square well potential search)
1467 C---NOTE THAT SCATT. AMPL.= 1/CMPLX(G(AK),-AK*ACH)
1468       IMPLICIT REAL*8 (A-H,O-Z)
1469       COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10),
1470      1              SBKRB(10),SDKK(10)
1471       COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
1472       COMMON/FSI_ACH/HPR,AC,ACH,ACHR,JJ,MSPIN
1473       COMMON/FSI_BP/B,P/FSI_DERIV/BPR,PPR/FSI_SHH/SH,CHH
1474       COMMON/FSI_DAK/DAK,HAC,IFUN
1475       HAC=.0D0
1476       ACH=1.D0
1477       IF(ICH.EQ.0)GOTO 1
1478       XH=AC*AK
1479       HCP=HC(XH)
1480       HPR=HCP+.1544313298D0
1481       ACH=ACP(XH)
1482       HAC=2*HCP/AC
1483    1  AKS=AK**2
1484       BK(JJ)=DSQRT(AKS+EB(JJ)**2) ! kappa=kp
1485       X=2*RB(JJ)/AC
1486       H=BK(JJ)*RB(JJ)             ! kp*d
1487       CALL GST(X,H)
1488       BRHO=B                      ! B(kp,d)
1489       SBKRB(JJ)=SH                ! kp*d*B(kp,d)
1490       CALL DERIW(X,H)
1491       BRHOP=BPR                   ! B'(kp,d)= dB(kp,r)/dln(r) at r=d
1492       H=AK*RB(JJ)
1493       CALL GST(X,H)
1494       CDK(JJ)=CHH                 ! ReG(k,d)
1495       BRHOS=B                     !  B(k,d)
1496       PRHOS=P                     !  P(k,d)
1497       SDK(JJ)=SH
1498       IF(ICH.EQ.0)GOTO 2
1499       SDK(JJ)=ACH*SH              ! ImG(k,d)
1500       IF(AK.EQ.0.D0.AND.AC.LT.0.D0)SDK(JJ)=3.14159*X*B
1501    2  SDKK(JJ)=RB(JJ)
1502       IF(AK.NE.0.D0)SDKK(JJ)=SH/AK ! d*B(k,d)
1503       CALL DERIW(X,H)              ! PPR=P'(k,d)= dP(k,r)/dln(r) at r=d
1504       ZZ=PPR-PRHOS
1505       IF(ICH.EQ.1)ZZ=ZZ+X*(BRHOS+BPR*(DLOG(DABS(X))+HPR))
1506 C   ZZ= P'(k,d)-P(k,d)+x*{B(k,d)+B'(k,d)*[ln!x!+2*C-1+h(k*ac)]}
1507       GG=(BRHOP*CDK(JJ)-BRHO*ZZ)/RB(JJ)
1508 C   GG= [B'(kp,d)*ReG(k,d)-B(kp,d)*ZZ]/d
1509       G=GG/(BRHO*BPR-BRHOP*BRHOS)
1510 C    G= GG/[B(kp,d)*B'(k,d)-B'(kp,d)*B(k,d)]
1511       RETURN
1512       END
1513       SUBROUTINE DERIW(X,H)
1514 C---CALLED BY F-N G(AK)
1515       IMPLICIT REAL*8 (A-H,O-Z)
1516       COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
1517       COMMON/FSI_BP/B,P/FSI_DERIV/BPR,PPR
1518       HH=.1D-3
1519       CALL GST(X,H-HH)
1520       Q1=P
1521       B1=B
1522       CALL GST(X,H+HH)
1523       HHH=HH+HH
1524       BPR=H*(B-B1)/HHH
1525       PPR=H*(P-Q1)/HHH
1526       IF(ICH.EQ.0)RETURN
1527       CALL GST(X-HH,H)
1528       Q1=P
1529       B1=B
1530       CALL GST(X+HH,H)
1531       BPR=BPR+X*(B-B1)/HHH
1532       PPR=PPR+X*(P-Q1)/HHH
1533       RETURN
1534       END
1535 C================================================================
1536
1537       SUBROUTINE READ_FILE(KEY,CH8,INT4,REAL8,IERR,NUNIT)
1538 C     ==========
1539 C
1540 C     Routine to read one parameter of the program in the file
1541 C     DATA NUNIT defined in FSI3B EXEC
1542 C     NUNIT=11 for IBM in Nantes, 4 for SUN in Prague
1543 C
1544 C     INPUT  : KEY (CHARACTER*10) :
1545 C     OUTPUT : case of KEY : CH8   : (CHARACTER*8)
1546 C                            INT4  : (INTEGER*4)
1547 C                            REAL8 : (REAL*8)
1548 C                     (only one of them)
1549 C              IERR (INTEGER) : 0 : no error
1550 C                               1 : key not found
1551
1552       CHARACTER*10 KEY,TEST
1553       CHARACTER*4  TYPE
1554       CHARACTER*8  CH8
1555       INTEGER*4    INT4
1556       REAL*8       REAL8
1557       INTEGER*4    IERR
1558
1559       IERR=0
1560       REWIND(NUNIT)
1561 1     READ(NUNIT,FMT='(A10,2X,A4)')TEST,TYPE
1562       IF (TEST.EQ.KEY) THEN
1563         BACKSPACE(NUNIT)
1564         IF (TYPE.EQ.'CHAR') READ(NUNIT,FMT='(18X,A8,54X)')CH8
1565         IF (TYPE.EQ.'INT4') READ(NUNIT,FMT='(18X,I8,54X)')INT4
1566         IF (TYPE.EQ.'REA8') READ(NUNIT,FMT='(18X,F10.5,52X)')REAL8
1567       ELSE
1568         IF (TEST.NE.'* E.O.F. *') THEN
1569           GOTO 1
1570         ELSE
1571           IERR=1
1572         ENDIF
1573       ENDIF
1574 c      IF(IERR.EQ.1)STOP
1575       RETURN
1576       END
1577