HBTAN/fsiini.F

   1
   2           SUBROUTINE fsiini
   3 C---Note:
   4 C-- ICH= 0 (1) if the Coulomb interaction is absent (present);
   5 C-- ISPIN= JJ= 1,2,..,MSPIN denote increasing values of the pair
   6 C-- total spin S.
   7 C-- To calculate the CF of two particles (with masses m1, m2 and
   8 C-- charges C1, C2) the following information is required:
   9 C-- AM= twice the reduced mass= 2*m1*m2/(m1+m2) in GeV/c^2,
  10 C-- DM= (m1-m2)/(m1+m2), required if NS=2;
  11 C-- AC= Bohr radius= 2*137.036*0.1973/(C1*C2*AMH) in fm;
  12 C-- AC > 1.D9 if C1*C2= 0, AC < 0 if C1*C2 < 0;
  13 C-- MSPIN= MSPINH(LL)= number of the values of the total pair spin S;
  14 C-- FD= FDH(LL,JJ), RD= RDH(LL,JJ)= scattering length and effective
  15 C-- radius for each value of the total pair spin S, JJ= 1,..,MSPIN;     ;
  16 C-- the corresponding square well parameters EB= EBH(LL,JJ), RB=
  17 C-- RBH(LL,JJ) (required if NS=1) may be calculated by sear.f;
  18 C-- if the effective range approximation is not valid (as is the case,
  19 C-- e.g., for two-pion system) a code for calculation of the scattering
  20 C-- amplitude should be supplemented;
  21 C-- RHO= RHOH(LL,JJ), SF= SFH(LL,JJ), SE= SEH(LL) are spin factors;
  22 C-- RHO= the probability that the spins j1 and j2 of the two particles
  23 C-- will combine in a total spin S;
  24 C-- RHO= (2*S+1)/[(2j1+1)*(2j2+1)] for unpolarized particles;
  25 C-- RHO= (1-P1*P2)/4 and (3+P1*P2)/4 correspond to S=0 and 1 in the
  26 C-- case of spin-1/2 particles with polarizations P1 and P2;
  27 C-----------------------------------------------------------------------
  28       IMPLICIT REAL*8 (A-H,O-Z)
  29       COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
  30       COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
  31      1               X,Y,Z,T,RP,RPS
  32       COMMON/FSI_SPIN/RHO(10)
  33       COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
  34       COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
  35       COMMON/FSI_FD/FD(10),RD(10)
  36       COMMON/FSI_C/C(10),AM,AMS,DM
  37       COMMON/FSI_CONS/PI,PI2,SPI,DR,W
  38            COMPLEX*16 C
  39       COMMON/FSI_AA/AA
  40       COMMON/FSI_AAPI/AAPI(20,2)/FSI_AAND/AAND(20,4)
  41       COMMON/FSI_AAPIN/AAPIN(20,2)
  42       COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10),
  43      1              SBKRB(10),SDKK(10)
  44       COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS
  45
  46 c      Include 'common_fsi_poc.inc'
  47 c      Include 'common_fsi_prf.inc'
  48 c      Include 'common_fsi_spin.inc'
  49 c      Include 'common_fsi_ach.inc'
  50 c      Include 'common_fsi_ns.inc'
  51 c      Include 'common_fsi_fd.inc'
  52 c-mlv      Include 'common_fsi_c.inc'
  53 c      Include 'common_fsi_cons.inc'
  54 c      Include 'common_fsi_aa.inc'
  55 c      Include 'common_fsi_aapi.inc'
  56 c      Include 'common_fsi_aapin.inc'
  57 c      Include 'common_fsi_sw.inc'
  58 c      Include 'common_fsi_aand.inc'
  59
  60
  61       DIMENSION FDH(30,10),RDH(30,10),EBH(30,10),RBH(30,10)
  62       DIMENSION RHOH(30,10)
  63       DIMENSION AM1H(30),AM2H(30),C1H(30),C2H(30),MSPINH(30)
  64 C============= declarations pour l'appel de READ_FILE()============
  65       CHARACTER*10 KEY
  66       CHARACTER*8  CH8
  67       INTEGER*4    INT4
  68       REAL*8       REAL8
  69       INTEGER*4    IERR
  70 C
  71 C--- mass of the first and second particle
  72       DATA AM1H/.93956563D0,.93827231D0,.93956563D0,3.72737978D0,
  73      C          .13957D0,.13498D0,.13957D0, .93956563D0, .93827231D0,
  74      C          4*.13957D0,4*.493677D0,
  75      C          2*1.87561339D0,2*2.80892165D0,2*.497672D0,
  76      C          1.87561339D0,3*.93827231D0,.93956563D0, 2*0.D0/
  77       DATA AM2H/.93956563D0,.93827231D0,.93827231D0,3.72737978D0,
  78      C          .13957D0,.13498D0,.13957D0, 2*1.87561339D0,
  79      C          2*.493677D0,2*.93827231D0,
  80      C          2*.493677D0,2*.93827231D0,
  81      C          1.87561339D0,3.72737978D0,2.80892165D0,3.72737978D0,
  82      C          2*.497672D0,2*2.80892165D0,3.72737978D0,
  83      C          2*1.115684D0,2*0.D0/
  84 c--------|---------|---------|---------|---------|---------|---------|----------
  85 C---  charge of the first and second particle
  86       DATA C1H/0.D0,1.D0,0.D0,2.D0, 1.D0,0.D0,1.D0,0.D0,1.D0,
  87      C         3*1.D0,-1.D0,3*1.D0,-1.D0,
  88      C         4*1.D0,2*0.D0,4*1.D0,0.D0, 2*0.D0/
  89       DATA C2H/0.D0,1.D0,1.D0,2.D0,-1.D0,0.D0,3*1.D0,
  90      C         -1.D0,3*1.D0,-1.D0,3*1.D0,
  91      C         1.D0,2.D0,1.D0,2.D0,2*0.D0,2*1.D0,2.D0,2*0.D0,2*0.D0/
  92 C---MSPIN vs (LL)
  93       DATA MSPINH/3*2,4*1,2*2,8*1,3,1,2,1,2*1,2*2,1,2*2, 2*0/
  94 C---Spin factors RHO vs (LL,ISPIN)
  95       DATA RHOH/3*.25D0, 4*1.D0, 2*.3333D0, 8*1.D0,
  96      1          .1111D0,1.D0,.25D0,1.D0,2*1.D0,
  97      1          .3333D0,.25D0,1.D0,2*.25D0, 2*0.D0,
  98      2          3*.75D0, 4*0.D0, 2*.6667D0, 8*0.D0,
  99      2          .3333D0,.0D0,.75D0,.0D0,2*0.D0,
 100      2          .6667D0,.75D0,0.D0,2*.75D0, 2*0.D0,
 101      3          17*.0D0,.5556D0,3*0.D0, 7*0.D0,2*0.D0,210*0.D0/
 102 C---Scattering length FD and effective radius RD in fm vs (LL,ISPIN)
 103       DATA FDH/17.0D0,7.8D0,23.7D0,2230.1218D0,.225D0,.081D0,-.063D0,
 104      1        -.65D0,-2.73D0,
 105      1        .137D0,-.071D0,-.148D0,.112D0,2*1.D-6,-.360D0,
 106      1        2*1.D-6,1.344D0,6*1.D-6,-5.628D0,2.18D0,2.40D0, 2*0.D0,
 107 cc     2 -10.8D0,2*-5.4D0,4*0.D0,-6.35D0,-11.88D0,8*0.D0,9*0.D0,
 108      2        3*-5.4D0,4*0.D0,-6.35D0,-11.88D0,8*0.D0,9*0.D0,
 109      2        1.93D0,1.84D0,2*0.D0,
 110      3        240*0.D0/
 111 c--------|---------|---------|---------|---------|---------|---------|----------
 112       DATA RDH/2.7D0,2.8D0,2.7D0,1.12139906D0,-44.36D0,64.0D0,784.9D0,
 113      1  477.9D0, 2.27D0, 9*0.D0,-69.973D0, 6*0.D0,3.529D0,
 114      1  3.19D0,3.15D0, 2*0.D0,
 115      2  3*1.7D0,4*0.D0,2.0D0,2.63D0, 17*0.D0,3.35D0,3.37D0, 2*0.D0,
 116      3  240*0.D0/
 117 C---Corresponding square well parameters RB (width in fm) and
 118 C-- EB =SQRT(-AM*U) (in GeV/c); U is the well height
 119       DATA RBH/2.545739D0,   2.779789D0, 2.585795D0, 5.023544D0,
 120      1 .124673D0, .3925180D0,.09D0, 2.D0, 4.058058D0, 17*0.D0,
 121      1  2.252623D0, 2.278575D0, 2*0.D0,
 122      2         3*2.003144D0,
 123      2         4*0.D0, 2.D0, 4.132163D0, 17*0.D0,
 124      2  2.272703D0, 2.256355D0, 2*0.D0,
 125      3         240*0.D0/
 126       DATA EBH/.1149517D0,    .1046257D0,   .1148757D0, .1186010D0,
 127      1    .7947389D0,2.281208D0,8.7D0,.4D0,.1561219D0,17*0.D0,
 128      1    .1013293D0, .1020966D0, 2*0.D0,
 129      2         3*.1847221D0,
 130      2         4*0.D0, .4D0, .1150687D0, 17*0.D0,
 131      2         .09736083D0, .09708310D0, 2*0.D0,
 132      3         240*0.D0/
 133 C=======< constants >========================
 134       W=1/.1973D0    ! from fm to 1/GeV
 135       PI=4*DATAN(1.D0)
 136       PI2=2*PI
 137       SPI=DSQRT(PI)
 138       DR=180.D0/PI   ! from radian to degree
 139
 140 c    WRITE(*,*)'from C++ to fortran W PI PI2 SPI DR',W,PI,PI2,SPI,DR
 141
 142         AC1=1.D10
 143         AC2=1.D10
 144 C=======< condition de calculs >=============
 145       NUNIT=11 ! for IBM or HP
 146 C      NUNIT=4 ! for SUN in Prague
 147 c-mlv      CALL readint4(NUNIT,'ITEST     ',ITEST)
 148 c-mlv      CALL readint4(NUNIT,'LL        ',LL)        ! Two-particle system
 149 c-mlv      CALL readint4(NUNIT,'NS        ',NS)
 150 c      CALL READ_FILE(NUNIT,'ITEST     ',CHAR,ITEST,REAL8,IERR)
 151 c      CALL READ_FILE(NUNIT,'LL        ',CHAR,LL,REAL8,IERR)
 152 c      CALL READ_FILE(NUNIT,'NS        ',CHAR,NS,REAL8,IERR)
 153
 154
 155
 156 C---setting particle masses and charges
 157       AM1=AM1H(LL)
 158       AM2=AM2H(LL)
 159       C1=C1H(LL)
 160       C2=C2H(LL)
 161
 162 C-- Switches:
 163 C   ISI=1(0)  the strong interaction between the two particles ON (OFF)
 164 C   IQS=1(0)  the quantum statistics ON (OFF);
 165 C             should be OFF for nonidentical particles
 166 C   I3C=1(0)  the Coulomb interaction with the nucleus ON (OFF)
 167 C   I3S=1(0)  the strong interaction with the nucleus ON (OFF)
 168 C   ICH=1(0)  if C1*C2 is different from 0 (is equal to 0)
 169 C-  to switch off the Coulomb force between the two particles
 170 C   put ICH=0 and substitute the strong amplitude parameters by
 171 C   the ones not affected by Coulomb interaction
 172
 173        IF(ITEST.EQ.0)THEN
 174        ICH=0
 175        IF(C1*C2.NE.0.D0) ICH=1
 176        IQS=0
 177        IF(C1+AM1.EQ.C2+AM2) IQS=1
 178        I3S=0  ! only this option is available
 179        ISI=1
 180        I3C=0
 181        IF(CN*ICH.NE.0.D0) I3C=1
 182        ENDIF
 183
 184
 185 c      IF(ITEST.EQ.1)THEN
 186 c      NS=4 ! SPHER. WAVE
 187 c      ICH=0
 188 c      IQS=1
 189 c      ISI=0
 190 c      I3C=0
 191
 192
 193 c-mlv      CALL readint4(NUNIT,'ICH       ',ICH)
 194 c-mlv      CALL readint4(NUNIT,'IQS       ',IQS)
 195 c-mlv      CALL readint4(NUNIT,'ISI       ',ISI)
 196 c-mlv      CALL readint4(NUNIT,'I3C       ',I3C)
 197 c      CALL READ_FILE(NUNIT,'ICH     ',CHAR,ICH,REAL8,IERR)
 198 c      CALL READ_FILE(NUNIT,'IQS     ',CHAR,IQS,REAL8,IERR)
 199 c      CALL READ_FILE(NUNIT,'ISI     ',CHAR,ISI,REAL8,IERR)
 200 c      CALL READ_FILE(NUNIT,'I3C     ',CHAR,I3C,REAL8,IERR)
 201 c      ENDIF
 202
 203       write(*,*)'====itest ll ich iqs isi i3c===',itest,ll,ich, iqs, isi, i3c
 204
 205 C==================================================================
 206 C---fm to 1/GeV
 207       DO 3 J1=1,30
 208       DO 3 J2=1,10
 209       FDH(J1,J2)=FDH(J1,J2)*W
 210       RDH(J1,J2)=RDH(J1,J2)*W
 211  3    RBH(J1,J2)=RBH(J1,J2)*W
 212 C---calcul. twice the reduced mass (AM), the relative mass difference
 213 C-- (DM) and the Bohr radius (AC)
 214       AM=2*AM1*AM2/(AM1+AM2)
 215       AMS=AM*AM
 216       DM=(AM1-AM2)/(AM1+AM2)
 217       AC=1.D10
 218       C12=C1*C2
 219       IF(C12.NE.0.D0)AC=2*137.036D0/(C12*AM)
 220 C---Setting spin factors
 221       MSPIN=MSPINH(LL)
 222       MSP=MSPIN
 223       DO 91 ISPIN=1,10
 224   91  RHO(ISPIN)=RHOH(LL,ISPIN)
 225 C---Integration limit AA in the spherical wave approximation
 226       AA=0.D0
 227 cc      IF(NS.EQ.2.OR.NS.EQ.4)AA=.5D0 !!in 1/GeV --> 0.1 fm
 228       IF(NS.EQ.2.OR.NS.EQ.4)AA=6.D0 !!in 1/GeV --> 1.2 fm
 229 C---Setting scatt. length (FD), eff. radius (RD) and, if possible,
 230 C-- also the corresp. square well parameters (EB, RB)
 231       DO 55 JJ=1,MSP
 232       ISPIN=JJ
 233       FD(JJ)=FDH(LL,JJ)
 234       RD(JJ)=RDH(LL,JJ)
 235       EB(JJ)=EBH(LL,JJ)
 236       RB(JJ)=RBH(LL,JJ)
 237 C---Resets FD and RD for a nucleon-deuteron system (LL=8,9)
 238        IF(LL.EQ.8.OR.LL.EQ.9)THEN
 239        JH=LL-7+2*JJ-2
 240        FD(JJ)=AAND(1,JH)
 241        RD(JJ)=AAND(2,JH)-2*AAND(3,JH)/AAND(1,JH)
 242        ENDIF
 243 C---Resets FD and RD for a pion-pion system (LL=5,6,7)
 244        IF(LL.EQ.5.OR.LL.EQ.6.OR.LL.EQ.7)THEN
 245        IF(LL.EQ.7)FD(JJ)=AAPI(1,2)/AM
 246        IF(LL.EQ.5)FD(JJ)=(.6667D0*AAPI(1,1)+.3333D0*AAPI(1,2))/AM
 247        IF(LL.EQ.6)FD(JJ)=(.3333D0*AAPI(1,1)+.6667D0*AAPI(1,2))/AM
 248        AKS=0.D0
 249        DAKS=1.D-5
 250        AKSH=AKS+DAKS
 251        AKH=DSQRT(AKSH)
 252        GPI1H=GPIPI(AKSH,1)
 253        GPI2H=GPIPI(AKSH,2)
 254        H=1/FD(JJ)
 255        IF(LL.EQ.7)C(JJ)=1/DCMPLX(GPI2H,-AKH)
 256        IF(LL.EQ.5)
 257      + C(JJ)=.6667D0/DCMPLX(GPI1H,-AKH)+.3333D0/DCMPLX(GPI2H,-AKH)
 258        IF(LL.EQ.6)
 259      + C(JJ)=.3333D0/DCMPLX(GPI1H,-AKH)+.6667D0/DCMPLX(GPI2H,-AKH)
 260        HH=DREAL(1/C(JJ))
 261        RD(JJ)=2*(HH-H)/DAKS
 262       ENDIF
 263 C---Resets FD and RD for a pion-nucleon system (LL=12,13)
 264       IF(LL.EQ.12.OR.LL.EQ.13)THEN
 265        IF(LL.EQ.12)FD(JJ)=AAPIN(1,2)
 266        IF(LL.EQ.13)FD(JJ)=(.6667D0*AAPIN(1,1)+.3333D0*AAPIN(1,2))
 267        AKS=0.D0
 268        DAKS=1.D-5
 269        AKSH=AKS+DAKS
 270        AKH=DSQRT(AKSH)
 271        GPI1H=GPIN(AKSH,1)
 272        GPI2H=GPIN(AKSH,2)
 273        H=1/FD(JJ)
 274        IF(LL.EQ.12)C(JJ)=1/DCMPLX(GPI2H,-AKH)
 275        IF(LL.EQ.13)
 276      + C(JJ)=.6667D0/DCMPLX(GPI1H,-AKH)+.3333D0/DCMPLX(GPI2H,-AKH)
 277        HH=DREAL(1/C(JJ))
 278        RD(JJ)=2*(HH-H)/DAKS
 279       ENDIF
 280 C---Calculation continues for any system (any LL)
 281  55   CONTINUE
 282       RETURN
 283       END
 284
 285
 286
 287        FUNCTION GPIPI(X,J)
 288 C--- GPIPI = k*COTG(DELTA), X=k^2
 289 C--  J=1(2) corresponds to isospin=0(2)
 290        IMPLICIT REAL*8 (A-H,O-Z)
 291 c--       Include 'common_fsi_aapi.inc'
 292 c--       Include 'common_fsi_c.inc'
 293        COMMON/FSI_AAPI/AAPI(20,2)
 294        COMMON/FSI_C/HELP(20),AM,AMS,DM
 295        OM=DSQRT(X+AMS)
 296        XX=X/AMS
 297        GPIPI=OM/AAPI(1,J)
 298        GPIPI=GPIPI*(1+(AAPI(3,J)-AAPI(1,J)**2)*XX+AAPI(4,J)*XX*XX)
 299        GPIPI=GPIPI/(1+(AAPI(3,J)+AAPI(2,J)/AAPI(1,J))*XX)
 300
 301        END
 302
 303        FUNCTION GPIN(X,J)
 304 C--- GPIN = k*COTG(DELTA), X=k^2
 305 C--  J=1(2) corresponds to piN isospin=1/2(3/2)
 306        IMPLICIT REAL*8 (A-H,O-Z)
 307 c--       Include 'common_fsi_aapin.inc'
 308        COMMON/FSI_AAPIN/AAPIN(20,2)
 309        GPIN=1/AAPIN(1,J)+.5D0*AAPIN(2,J)*X
 310
 311        END