Update by Ludmila Malinina: Code rearranged for easier comparison with original
authorskowron <skowron@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 31 Oct 2002 13:38:30 +0000 (13:38 +0000)
committerskowron <skowron@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 31 Oct 2002 13:38:30 +0000 (13:38 +0000)
HBTAN/fsiini.F
HBTAN/fsiw.F
HBTAN/ltran12.F

index 810e717..f3c1421 100644 (file)
@@ -1,5 +1,5 @@
-          
-          SUBROUTINE fsiini
+
+          SUBROUTINE FSIINI
 C---Note:
 C-- ICH= 0 (1) if the Coulomb interaction is absent (present);
 C-- ISPIN= JJ= 1,2,..,MSPIN denote increasing values of the pair
@@ -14,7 +14,7 @@ C-- MSPIN= MSPINH(LL)= number of the values of the total pair spin S;
 C-- FD= FDH(LL,JJ), RD= RDH(LL,JJ)= scattering length and effective
 C-- radius for each value of the total pair spin S, JJ= 1,..,MSPIN;     ;
 C-- the corresponding square well parameters EB= EBH(LL,JJ), RB=
-C-- RBH(LL,JJ) (required if NS=1) may be calculated by sear.f;
+C-- RBH(LL,JJ) (required if NS=1) may be calculated by SEAR;
 C-- if the effective range approximation is not valid (as is the case,
 C-- e.g., for two-pion system) a code for calculation of the scattering
 C-- amplitude should be supplemented;
@@ -35,28 +35,13 @@ C-----------------------------------------------------------------------
       COMMON/FSI_FD/FD(10),RD(10)
       COMMON/FSI_C/C(10),AM,AMS,DM
       COMMON/FSI_CONS/PI,PI2,SPI,DR,W
-           COMPLEX*16 C
+      COMPLEX*16 C
       COMMON/FSI_AA/AA
       COMMON/FSI_AAPI/AAPI(20,2)/FSI_AAND/AAND(20,4)
-      COMMON/FSI_AAPIN/AAPIN(20,2)
       COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10),
      1              SBKRB(10),SDKK(10)
-      COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS 
-
-c      Include 'common_fsi_poc.inc'
-c      Include 'common_fsi_prf.inc'
-c      Include 'common_fsi_spin.inc'
-c      Include 'common_fsi_ach.inc'
-c      Include 'common_fsi_ns.inc'
-c      Include 'common_fsi_fd.inc'
-c-mlv      Include 'common_fsi_c.inc'
-c      Include 'common_fsi_cons.inc'
-c      Include 'common_fsi_aa.inc'
-c      Include 'common_fsi_aapi.inc'
-c      Include 'common_fsi_aapin.inc'
-c      Include 'common_fsi_sw.inc'
-c      Include 'common_fsi_aand.inc'
 
+      COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS 
 
       DIMENSION FDH(30,10),RDH(30,10),EBH(30,10),RBH(30,10)
       DIMENSION RHOH(30,10)
@@ -69,89 +54,63 @@ C============= declarations pour l'appel de READ_FILE()============
       INTEGER*4    IERR
 C
 C--- mass of the first and second particle
-      DATA AM1H/.93956563D0,.93827231D0,.93956563D0,3.72737978D0,
-     C          .13957D0,.13498D0,.13957D0, .93956563D0, .93827231D0,
+      DATA AM1H/.9395656D0,.9382723D0,.9395656D0,3.7294D0,.13957D0,
+     C          .13498D0,.13957D0, .9395656D0, .9382723D0,
      C          4*.13957D0,4*.493677D0,
-     C          2*1.87561339D0,2*2.80892165D0,2*.497672D0,
-     C          1.87561339D0,3*.93827231D0,.93956563D0, 2*0.D0/
-      DATA AM2H/.93956563D0,.93827231D0,.93827231D0,3.72737978D0,
-     C          .13957D0,.13498D0,.13957D0, 2*1.87561339D0,
-     C          2*.493677D0,2*.93827231D0,
-     C          2*.493677D0,2*.93827231D0,
-     C          1.87561339D0,3.72737978D0,2.80892165D0,3.72737978D0,
-     C          2*.497672D0,2*2.80892165D0,3.72737978D0,
-     C          2*1.115684D0,2*0.D0/
-c--------|---------|---------|---------|---------|---------|---------|----------
+     C          2*1.875613D0,2*2.808D0,2*.497672D0,
+     C          1.875613D0,2*.9382723D0, 4*0.D0/
+      DATA AM2H/.9395656D0,.9382723D0,.9382723D0,3.7294D0,.13957D0,
+     C          .13498D0,.13957D0, 2*1.875613D0,
+     C          2*.493677D0,2*.9382723D0,
+     C          2*.493677D0,2*.9382723D0,
+     C          1.875613D0,3.7294D0,2.808D0,3.7294D0,
+     C          2*.497672D0,2*2.808D0,3.7294D0, 4*0.D0/
 C---  charge of the first and second particle
       DATA C1H/0.D0,1.D0,0.D0,2.D0, 1.D0,0.D0,1.D0,0.D0,1.D0,
      C         3*1.D0,-1.D0,3*1.D0,-1.D0,
-     C         4*1.D0,2*0.D0,4*1.D0,0.D0, 2*0.D0/
+     C         4*1.D0,2*0.D0,3*1.D0, 4*0.D0/
       DATA C2H/0.D0,1.D0,1.D0,2.D0,-1.D0,0.D0,3*1.D0,
      C         -1.D0,3*1.D0,-1.D0,3*1.D0,
-     C         1.D0,2.D0,1.D0,2.D0,2*0.D0,2*1.D0,2.D0,2*0.D0,2*0.D0/
+     C         1.D0,2.D0,1.D0,2.D0,2*0.D0,2*1.D0,2.D0, 4*0.D0/
 C---MSPIN vs (LL)
-      DATA MSPINH/3*2,4*1,2*2,8*1,3,1,2,1,2*1,2*2,1,2*2, 2*0/
+      DATA MSPINH/3*2,4*1,2*2,8*1,3,1,2,1,2*1,2*2,1, 4*0/
 C---Spin factors RHO vs (LL,ISPIN)
       DATA RHOH/3*.25D0, 4*1.D0, 2*.3333D0, 8*1.D0, 
-     1          .1111D0,1.D0,.25D0,1.D0,2*1.D0,
-     1          .3333D0,.25D0,1.D0,2*.25D0, 2*0.D0,
-     2          3*.75D0, 4*0.D0, 2*.6667D0, 8*0.D0, 
-     2          .3333D0,.0D0,.75D0,.0D0,2*0.D0,
-     2          .6667D0,.75D0,0.D0,2*.75D0, 2*0.D0, 
-     3          17*.0D0,.5556D0,3*0.D0, 7*0.D0,2*0.D0,210*0.D0/
+     C          .1111D0,1.D0,.25D0,1.D0,2*1.D0,
+     1          .3333D0,.25D0,1.D0, 4*0.D0,
+     C          3*.75D0, 4*0.D0, 2*.6667D0, 8*0.D0, 
+     C          .3333D0,.0D0,.75D0,.0D0,2*0.D0,
+     2          .6667D0,.75D0,0.D0, 4*0.D0, 
+     C          17*.0D0,.5556D0,3*0.D0, 5*0.D0,4*0.D0,210*0.D0/
 C---Scattering length FD and effective radius RD in fm vs (LL,ISPIN)
       DATA FDH/17.0D0,7.8D0,23.7D0,2230.1218D0,.225D0,.081D0,-.063D0,
-     1        -.65D0,-2.73D0,
-     1        .137D0,-.071D0,-.148D0,.112D0,2*1.D-6,-.360D0,
-     1        2*1.D-6,1.344D0,6*1.D-6,-5.628D0,2.18D0,2.40D0, 2*0.D0,
-cc     2 -10.8D0,2*-5.4D0,4*0.D0,-6.35D0,-11.88D0,8*0.D0,9*0.D0,
-     2        3*-5.4D0,4*0.D0,-6.35D0,-11.88D0,8*0.D0,9*0.D0,
-     2        1.93D0,1.84D0,2*0.D0,
-     3        240*0.D0/
-c--------|---------|---------|---------|---------|---------|---------|----------     
+     C        -.65D0,-2.73D0,
+     C        .137D0,-.071D0,-.148D0,.112D0,2*1.D-6,-.360D0,
+     1        2*1.D-6,1.344D0,6*1.D-6,-5.628D0, 4*0.D0,
+     C -10.8D0,2*-5.4D0,4*0.D0,-6.35D0,-11.88D0,8*0.D0,9*0.D0,4*0.D0,
+     C 240*0.D0/
       DATA RDH/2.7D0,2.8D0,2.7D0,1.12139906D0,-44.36D0,64.0D0,784.9D0,
-     1  477.9D0, 2.27D0, 9*0.D0,-69.973D0, 6*0.D0,3.529D0,
-     1  3.19D0,3.15D0, 2*0.D0,
-     2  3*1.7D0,4*0.D0,2.0D0,2.63D0, 17*0.D0,3.35D0,3.37D0, 2*0.D0, 
-     3  240*0.D0/
+c--------|---------|---------|---------|---------|---------|---------|----------      
+     C  477.9D0, 2.27D0, 9*0.D0,-69.973D0, 6*0.D0,3.529D0, 4*0.D0,
+     C       3*1.7D0,4*0.D0,2.0D0,2.63D0, 17*0.D0, 4*0.D0, 240*0.D0/
 C---Corresponding square well parameters RB (width in fm) and
 C-- EB =SQRT(-AM*U) (in GeV/c); U is the well height
       DATA RBH/2.545739D0,   2.779789D0, 2.585795D0, 5.023544D0,
-     1 .124673D0, .3925180D0,.09D0, 2.D0, 4.058058D0, 17*0.D0, 
-     1  2.252623D0, 2.278575D0, 2*0.D0,
-     2         3*2.003144D0,
-     2         4*0.D0, 2.D0, 4.132163D0, 17*0.D0, 
-     2  2.272703D0, 2.256355D0, 2*0.D0, 
-     3         240*0.D0/
+     C .124673D0, .3925180D0,.09D0, 2.D0, 4.058058D0, 17*0.D0, 4*0.D0,
+     C         3*2.003144D0,
+     C         4*0.D0, 2.D0, 4.132163D0, 17*0.D0, 4*0.D0, 240*0.D0/
       DATA EBH/.1149517D0,    .1046257D0,   .1148757D0, .1186010D0,
-     1    .7947389D0,2.281208D0,8.7D0,.4D0,.1561219D0,17*0.D0,
-     1    .1013293D0, .1020966D0, 2*0.D0,
-     2         3*.1847221D0,
-     2         4*0.D0, .4D0, .1150687D0, 17*0.D0, 
-     2         .09736083D0, .09708310D0, 2*0.D0, 
-     3         240*0.D0/
+     C    .7947389D0,2.281208D0,8.7D0,.4D0,.1561219D0,17*0.D0,4*0.D0,
+     C         3*.1847221D0,
+     C         4*0.D0, .4D0, .1150687D0, 17*0.D0, 4*0.D0, 240*0.D0/
 C=======< constants >========================
       W=1/.1973D0    ! from fm to 1/GeV
       PI=4*DATAN(1.D0)
       PI2=2*PI
       SPI=DSQRT(PI)
       DR=180.D0/PI   ! from radian to degree
-
-c    WRITE(*,*)'from C++ to fortran W PI PI2 SPI DR',W,PI,PI2,SPI,DR
-
         AC1=1.D10
         AC2=1.D10
-C=======< condition de calculs >=============
-      NUNIT=11 ! for IBM or HP
-C      NUNIT=4 ! for SUN in Prague
-c-mlv      CALL readint4(NUNIT,'ITEST     ',ITEST)      
-c-mlv      CALL readint4(NUNIT,'LL        ',LL)        ! Two-particle system
-c-mlv      CALL readint4(NUNIT,'NS        ',NS) 
-c      CALL READ_FILE(NUNIT,'ITEST     ',CHAR,ITEST,REAL8,IERR)
-c      CALL READ_FILE(NUNIT,'LL        ',CHAR,LL,REAL8,IERR)
-c      CALL READ_FILE(NUNIT,'NS        ',CHAR,NS,REAL8,IERR)
-
-
  
 C---setting particle masses and charges
       AM1=AM1H(LL)
@@ -170,38 +129,22 @@ C-  to switch off the Coulomb force between the two particles
 C   put ICH=0 and substitute the strong amplitude parameters by
 C   the ones not affected by Coulomb interaction
  
-       IF(ITEST.EQ.0)THEN
+      IF(ITEST.EQ.0)THEN
+
        ICH=0
        IF(C1*C2.NE.0.D0) ICH=1
        IQS=0
        IF(C1+AM1.EQ.C2+AM2) IQS=1
        I3S=0  ! only this option is available
        ISI=1
-       I3C=0
-       IF(CN*ICH.NE.0.D0) I3C=1
+       I3C=1
+       
        ENDIF
-
        
-c      IF(ITEST.EQ.1)THEN
-c      NS=4 ! SPHER. WAVE
-c      ICH=0
-c      IQS=1
-c      ISI=0
-c      I3C=0
-      
-      
-c-mlv      CALL readint4(NUNIT,'ICH       ',ICH)
-c-mlv      CALL readint4(NUNIT,'IQS       ',IQS)
-c-mlv      CALL readint4(NUNIT,'ISI       ',ISI)
-c-mlv      CALL readint4(NUNIT,'I3C       ',I3C)
-c      CALL READ_FILE(NUNIT,'ICH     ',CHAR,ICH,REAL8,IERR)
-c      CALL READ_FILE(NUNIT,'IQS     ',CHAR,IQS,REAL8,IERR)
-c      CALL READ_FILE(NUNIT,'ISI     ',CHAR,ISI,REAL8,IERR)
-c      CALL READ_FILE(NUNIT,'I3C     ',CHAR,I3C,REAL8,IERR)
-c      ENDIF
-      
-      write(*,*)'====itest ll ich iqs isi i3c===',itest,ll,ich, iqs, isi, i3c
-      
+c23456       
+       write(*,*)'FSIINI ITEST ich iqs i3s isi i3c',ITEST, 
+     +  ICH,IQS,I3S,ISI,I3C
 C==================================================================
 C---fm to 1/GeV
       DO 3 J1=1,30
@@ -234,57 +177,48 @@ C-- also the corresp. square well parameters (EB, RB)
       RD(JJ)=RDH(LL,JJ)
       EB(JJ)=EBH(LL,JJ)
       RB(JJ)=RBH(LL,JJ)
+      IF(LL.NE.8.AND.LL.NE.9)GOTO25
 C---Resets FD and RD for a nucleon-deuteron system (LL=8,9)
-       IF(LL.EQ.8.OR.LL.EQ.9)THEN
-       JH=LL-7+2*JJ-2
-       FD(JJ)=AAND(1,JH)
-       RD(JJ)=AAND(2,JH)-2*AAND(3,JH)/AAND(1,JH)
-       ENDIF
+      JH=LL-7+2*JJ-2
+      FD(JJ)=AAND(1,JH)
+      RD(JJ)=AAND(2,JH)-2*AAND(3,JH)/AAND(1,JH)
 C---Resets FD and RD for a pion-pion system (LL=5,6,7)
-       IF(LL.EQ.5.OR.LL.EQ.6.OR.LL.EQ.7)THEN
-       IF(LL.EQ.7)FD(JJ)=AAPI(1,2)/AM
-       IF(LL.EQ.5)FD(JJ)=(.6667D0*AAPI(1,1)+.3333D0*AAPI(1,2))/AM
-       IF(LL.EQ.6)FD(JJ)=(.3333D0*AAPI(1,1)+.6667D0*AAPI(1,2))/AM
-       AKS=0.D0
-       DAKS=1.D-5
-       AKSH=AKS+DAKS
-       AKH=DSQRT(AKSH)
-       GPI1H=GPIPI(AKSH,1)
-       GPI2H=GPIPI(AKSH,2)
-       H=1/FD(JJ)
-       IF(LL.EQ.7)C(JJ)=1/DCMPLX(GPI2H,-AKH)
-       IF(LL.EQ.5)
-     + C(JJ)=.6667D0/DCMPLX(GPI1H,-AKH)+.3333D0/DCMPLX(GPI2H,-AKH)
-       IF(LL.EQ.6)
-     + C(JJ)=.3333D0/DCMPLX(GPI1H,-AKH)+.6667D0/DCMPLX(GPI2H,-AKH)
-       HH=DREAL(1/C(JJ))
-       RD(JJ)=2*(HH-H)/DAKS
-      ENDIF
-C---Resets FD and RD for a pion-nucleon system (LL=12,13) 
-      IF(LL.EQ.12.OR.LL.EQ.13)THEN
-       IF(LL.EQ.12)FD(JJ)=AAPIN(1,2)
-       IF(LL.EQ.13)FD(JJ)=(.6667D0*AAPIN(1,1)+.3333D0*AAPIN(1,2))
-       AKS=0.D0
-       DAKS=1.D-5
-       AKSH=AKS+DAKS
-       AKH=DSQRT(AKSH)
-       GPI1H=GPIN(AKSH,1)
-       GPI2H=GPIN(AKSH,2)
-       H=1/FD(JJ)
-       IF(LL.EQ.12)C(JJ)=1/DCMPLX(GPI2H,-AKH)
-       IF(LL.EQ.13)
-     + C(JJ)=.6667D0/DCMPLX(GPI1H,-AKH)+.3333D0/DCMPLX(GPI2H,-AKH)
-       HH=DREAL(1/C(JJ))
-       RD(JJ)=2*(HH-H)/DAKS
-      ENDIF
+  25  IF(LL.GT.7.OR.LL.LT.5)GOTO 24
+      IF(LL.EQ.7)FD(JJ)=AAPI(1,2)/AM
+      IF(LL.EQ.5)FD(JJ)=(.6667D0*AAPI(1,1)+.3333D0*AAPI(1,2))/AM
+      IF(LL.EQ.6)FD(JJ)=(.3333D0*AAPI(1,1)+.6667D0*AAPI(1,2))/AM
+      AKS=0.D0
+      DAKS=1.D-5
+      AKSH=AKS+DAKS
+      AKH=DSQRT(AKSH)
+      GPI1H=GPIPI(AKSH,1)
+      GPI2H=GPIPI(AKSH,2)
+      H=1/FD(JJ)
+      IF(LL.EQ.7)C(JJ)=1/DCMPLX(GPI2H,-AKH)
+      IF(LL.EQ.5)
+     +C(JJ)=.6667D0/DCMPLX(GPI1H,-AKH)+.3333D0/DCMPLX(GPI2H,-AKH)
+      IF(LL.EQ.6)
+     +C(JJ)=.3333D0/DCMPLX(GPI1H,-AKH)+.6667D0/DCMPLX(GPI2H,-AKH)
+      HH=DREAL(1/C(JJ))
+      RD(JJ)=2*(HH-H)/DAKS
+ 24   CONTINUE
 C---Calculation continues for any system (any LL)
  55   CONTINUE
-      RETURN
       END
+C=======================================================
+
+
+
+
+
+
+
+
+
 
 
       
-       FUNCTION GPIPI(X,J)                                                                    
+       FUNCTION GPIPIold(X,J)                                                                    
 C--- GPIPI = k*COTG(DELTA), X=k^2                                                            
 C--  J=1(2) corresponds to isospin=0(2)                                                      
        IMPLICIT REAL*8 (A-H,O-Z)                                                              
index e9bbbb0..15823c3 100644 (file)
@@ -122,32 +122,30 @@ C=======================================================================
      1                X2,Y2,Z2,T2,R2
       COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
       COMMON/FSI_FFPN/FF12,FF21
-      COMPLEX*16 FF12,FF21
       COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS
-       
+      COMPLEX*16 FF12,FF21
 C------------------------------------------------------------------
-c      write(*,*)'in fsiw C1 C2 CN',C1,C2,CN 
-
 C==> AC1,2 = "relativistic" Bohr radii for particle-nucleus systems
-
       C1N=C1*CN
       IF(C1N.NE.0.D0)AC1=137.036D0/(C1N*E1) !m1-->E1
       C2N=C2*CN
       IF(C2N.NE.0.D0)AC2=137.036D0/(C2N*E2) !m2-->E2
  
 C-----------------------------------------------------------
-      CALL FSIPN  !(WEIF) !weight due to particle-nucleus FSI
+         J=1
+      CALL FSIPN !weight due to particle-nucleus FSI
       JRAT=0
-      CALL FSIWF  !(WEI)  !weight due to particle-particle-nucleus FSI
+      CALL FSIWF  !weight due to particle-particle-nucleus FSI
       WEIN=WEI
-       IF(I3C*J.NE.0) THEN
+         IF(I3C*J.NE.0) THEN
       FF12=DCMPLX(1.D0,0.D0)
       FF21=DCMPLX(1.D0,0.D0)
       JRAT=1
-      CALL VZ    !(WEIN) ! weight due to particle-particle FSI
-       ENDIF
+      CALL VZ ! weight due to particle-particle FSI
+         ENDIF
       RETURN
       END
+
 C=======================================================================
  
       SUBROUTINE LTRAN(P0,P,PS)
@@ -190,20 +188,6 @@ C-----------------------------------------------------------------------
       RETURN
       END
  
-      FUNCTION GND(X,J)
-C--- GND = k*COTG(DELTA), X=k^2
-C--- J=1(2) corresp. to nd(pd), S=1/2,
-C--- J=3(4) corresp. to nd(pd), S=3/2
-      IMPLICIT REAL*8 (A-H,O-Z)
-      COMMON/FSI_AAND/AAND(20,4)
-      XX=X
-      GND=1/AAND(1,J)+.5D0*AAND(2,J)*X
-      DO 1 I=4,4
-      XX=XX*X
-   1  GND=GND+AAND(I,J)*XX
-      GND=GND/(1+AAND(3,J)*X)
-      RETURN
-      END
       
       FUNCTION GDD(X,J)
 C--- GDD = k*COTG(DELTA), X=k^2
@@ -235,7 +219,96 @@ C--- J=1,2,3 corresp. to S=0,1,2
       RETURN
       END
       
-      
+
+      SUBROUTINE FSIPN
+C  calculating particle-nucleus Coulomb Wave functions FFij
+      IMPLICIT REAL*8 (A-H,O-Z)
+      COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
+      COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1,  !part. momenta in NRF
+     1               P2X,P2Y,P2Z,E2,P2
+      COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emis. points in NRF
+     1                X2,Y2,Z2,T2,R2
+      COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
+      COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
+      COMMON/FSI_ICH1/ICH1
+      COMMON/FSI_ETA/ETA
+      COMMON/FSI_FFPN/FF12,FF21
+      COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS
+
+      COMPLEX*16 FF1,FF12,FF21
+      FF12=DCMPLX(1.D0,0.D0)
+      FF21=DCMPLX(1.D0,0.D0)
+      IF(I3C.EQ.0)RETURN
+      ICH1=IDINT(C1)
+      IF(ICH1.EQ.0)GOTO 11
+      XH=AC1*P1
+      ACH=ACP(XH)
+      ACHR=DSQRT(ACH)
+      ETA=0.D0
+      IF(XH.NE.0.D0)ETA=1/XH
+      RHOS=P1*R1
+      HS=X1*P1X+Y1*P1Y+Z1*P1Z
+      FF12=FF12*FF1(RHOS,HS)
+      IF(IQS.EQ.0)GOTO 11
+      RHOS=P1*R2
+      HS=X2*P1X+Y2*P1Y+Z2*P1Z
+      FF21=FF21*FF1(RHOS,HS)
+ 11   ICH1=IDINT(C2)
+      IF(ICH1.EQ.0)GOTO 10
+      XH=AC2*P2
+      ACH=ACP(XH)
+      ACHR=DSQRT(ACH)
+      ETA=0.D0
+      IF(XH.NE.0.D0)ETA=1/XH
+      RHOS=P2*R2
+      HS=X2*P2X+Y2*P2Y+Z2*P2Z
+      FF12=FF12*FF1(RHOS,HS)
+CW      WRITE(6,41)'AC2 ',AC2,'ACH ',ACH,'ETA ',ETA,'RHOS ',RHOS,'HS ',HS
+41    FORMAT(5(A5,E11.4))
+CW      WRITE(6,40)'FF12 ',DREAL(FF12),DIMAG(FF12)
+      IF(IQS.EQ.0)GOTO 10
+      RHOS=P2*R1
+      HS=X1*P2X+Y1*P2Y+Z1*P2Z
+      FF21=FF21*FF1(RHOS,HS)
+CW      WRITE(6,41)'AC1 ',AC1,'ACH ',ACH,'ETA ',ETA,'RHOS ',RHOS,'HS ',HS
+CW      WRITE(6,40)'FF21 ',DREAL(FF21),DIMAG(FF21)
+40    FORMAT(A7,2E12.4)
+ 10   CONTINUE
+C  WEIF = the weight due to the Coulomb particle-nucleus interaction
+      WEIF=DREAL(FF12)**2+DIMAG(FF12)**2
+      IF(IQS.EQ.1)WEIF=0.5D0*(WEIF+DREAL(FF21)**2+DIMAG(FF21)**2)
+      RETURN
+      END
+      FUNCTION GPIPI(X,J)
+C--- GPIPI = k*COTG(DELTA), X=k^2
+C--  J=1(2) corresponds to isospin=0(2)
+      IMPLICIT REAL*8 (A-H,O-Z)
+      COMMON/FSI_AAPI/AAPI(20,2)
+      COMMON/FSI_C/HELP(20),AM,AMS,DM
+      OM=DSQRT(X+AMS)
+      XX=X/AMS
+      GPIPI=OM/AAPI(1,J)
+      GPIPI=GPIPI*(1+(AAPI(3,J)-AAPI(1,J)**2)*XX+AAPI(4,J)*XX*XX)
+      GPIPI=GPIPI/(1+(AAPI(3,J)+AAPI(2,J)/AAPI(1,J))*XX)
+      RETURN
+      END
+      FUNCTION GND(X,J)
+C--- GND = k*COTG(DELTA), X=k^2
+C--- J=1(2) corresp. to nd(pd), S=1/2,
+C--- J=3(4) corresp. to nd(pd), S=3/2
+      IMPLICIT REAL*8 (A-H,O-Z)
+      COMMON/FSI_AAND/AAND(20,4)
+      XX=X
+      GND=1/AAND(1,J)+.5D0*AAND(2,J)*X
+      DO 1 I=4,4
+      XX=XX*X
+   1  GND=GND+AAND(I,J)*XX
+      GND=GND/(1+AAND(3,J)*X)
+      RETURN
+      END
+
       SUBROUTINE CKKB  ! calculates KK-b scattering amplitude,
                        ! saturated by S*(980) and delta(982) resonances
       IMPLICIT REAL*8 (A-H,O-Z)
@@ -254,9 +327,136 @@ C--- J=1,2,3 corresp. to S=0,1,2
       C(1)=C(1)+AAKK(8)/2/DCMPLX(AAKK(5)-S,
      ,-AK*AAKK(8)-AKPIETA*AAKK(9))
       RETURN
-      END 
-          
-      SUBROUTINE VZ  !(WEI)
+      END      
+
+
+C
+      SUBROUTINE FSIWF
+C==> Prepares necessary quantities and call VZ(WEI) to calculate
+C    the weight due to FSI
+      IMPLICIT REAL*8 (A-H,O-Z)
+      COMMON/FSI_CVK/V,CVK
+      COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1,  !part. momenta in NRF
+     1               P2X,P2Y,P2Z,E2,P2
+      COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
+     1               X,Y,Z,T,RP,RPS
+      COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emis. points in NRF
+     1                X2,Y2,Z2,T2,R2
+      COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
+      COMMON/FSI_SPIN/RHO(10)
+      COMMON/FSI_BP/B,P
+      COMMON/FSI_ETA/ETA
+      COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
+      COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10),
+     1              SBKRB(10),SDKK(10)
+      COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
+      COMMON/FSI_RR/F(10)
+      COMMON/FSI_FD/FD(10),RD(10)
+      COMMON/FSI_C/C(10),AM,AMS,DM
+      COMPLEX*16 C,F
+      COMMON/FSI_AA/AA
+      COMMON/FSI_SHH/SH,CHH
+      COMMON/FSI_AAPI/AAPI(20,2)/FSI_AAND/AAND(20,4)
+      COMMON/FSI_P12/P12X,P12Y,P12Z,E12,P12,AM12,EPM
+      COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS
+
+C==>calculating relative 4-coordinates of the particles in PRF
+C-  {T,X,Y,Z} from the relative coordinates {TS,XS,YS,ZS} in NRF
+      XS=X1-X2
+      YS=Y1-Y2
+      ZS=Z1-Z2
+      TS=T1-T2
+      RS12=XS*P12X+YS*P12Y+ZS*P12Z
+      H1=(RS12/EPM-TS)/AM12
+      
+c-mlv      X=XS+P12X*H1
+c-mlv      Y=YS+P12Y*H1
+c-mlv      Z=ZS+P12Z*H1
+c-mlv      T=(E12*TS-RS12)/AM12
+c-mlv      RPS=X*X+Y*Y+Z*Z
+c-mlv      RP=DSQRT(RPS)
+c      WRITE(6,38)'RP ',RP,'X ',X,Y,Z,T
+c38    FORMAT(A7,E11.4,A7,4E11.4)
+      CVK=(P12X*PPX+P12Y*PPY+P12Z*PPZ)/(P12*AK)
+      V=P12/E12
+      IF(ICH.EQ.0)GOTO 21
+      XH=AC*AK
+      ACH=ACP(XH)
+      ACHR=DSQRT(ACH)
+      ETA=0.D0
+      IF(XH.NE.0.D0)ETA=1/XH
+C---HCP, HPR needed (e.g. in GST) if ICH=1
+      HCP=HC(XH)
+      HPR=HCP+.1544313298D0
+  21  CONTINUE
+      MSP=MSPIN
+      DO 30 JJ=1,MSP
+      IF(NS.NE.1)GOTO22
+C---Calc. quantities for the square well potential;
+C-- for LL > 5 the square well potential is not possible or available
+cc      GAK=G(AK)
+      BK(JJ)=DSQRT(EB(JJ)**2+AKS)
+      XRA=2*RB(JJ)/AC
+      HRA=BK(JJ)*RB(JJ)
+      CALL SEQ(XRA,HRA)
+      SBKRB(JJ)=HRA*B
+      HRA=AK*RB(JJ)
+      CALL GST(XRA,HRA)
+      SDK(JJ)=SH
+      CDK(JJ)=CHH
+      SDKK(JJ)=RB(JJ)
+      IF(AK.NE.0.D0)SDKK(JJ)=SH/AK
+      IF(ICH.EQ.1)SDK(JJ)=ACH*SDK(JJ)
+  22  CONTINUE
+C-----------------------------------------------------------------------
+C---Calc. the strong s-wave scattering amplitude = C(JJ)
+C-- divided by Coulomb penetration factor squared (if ICH=1)
+      IF(NS.NE.1)GOTO 230
+      IF(LL.NE.4)GOTO 230 ! SW scat. amplitude used for alfa-alfa only
+      AKACH=AK
+      IF(ICH.EQ.1)AKACH=AK*ACH
+      C(JJ)=1/DCMPLX(GAK,-AKACH) ! amplitude for the SW-potential
+      GOTO 30
+ 230  IF(LL.EQ.5.OR.LL.EQ.6.OR.LL.EQ.7)GOTO20 ! pipi      
+      IF(LL.EQ.8.OR.LL.EQ.9)GOTO20   ! Nd
+      IF(LL.EQ.14.OR.LL.EQ.17.OR.LL.EQ.23)GOTO27 ! K+K-, K-p, K0K0-b
+      A1=RD(JJ)*FD(JJ)*AKS
+      A2=1+.5D0*A1
+      IF(ICH.EQ.1)A2=A2-2*HCP*FD(JJ)/AC
+      AKF=AK*FD(JJ)
+      IF(ICH.EQ.1)AKF=AKF*ACH
+      C(JJ)=FD(JJ)/DCMPLX(A2,-AKF)
+      GOTO30
+ 20   CONTINUE
+C---Calc. scatt. ampl. C(JJ) for pipi and Nd
+      JH=LL-7+2*JJ-2
+      IF(LL.GT.7)GPI2=GND(AKS,JH)
+      IF(LL.LE.7)GPI2=GPIPI(AKS,2)
+      C(JJ)=1.D0/DCMPLX(GPI2,-AK)
+      IF(LL.GE.7)GOTO27
+      GPI1=GPIPI(AKS,1)
+      IF(LL.EQ.5)C(JJ)=.6667D0/DCMPLX(GPI1,-AK)+.3333D0*C(JJ)
+      IF(LL.EQ.6)C(JJ)=.3333D0/DCMPLX(GPI1,-AK)+.6667D0*C(JJ)
+ 27   CONTINUE
+C---Calc. K+K-, K0K0-b or K-p s-wave scatt. ampl.
+      IF(LL.EQ.14.OR.LL.EQ.23)CALL CKKB
+      IF(LL.EQ.17)C(JJ)=DCMPLX(3.29D0,3.55D0)
+C---Calc. pi+pi-, pi+pi+, pd, K+K- or K-p s-wave scatt. amplitude
+C-- divided by Coulomb penetration factor squared (if ICH=1)
+      IF(ICH.EQ.0)GOTO 30
+      AAK=ACH*AK
+      HCP2=2*HCP/AC
+      C(JJ)=1/(1/C(JJ)-HCP2+DCMPLX(0.D0,AK-AAK))
+ 30   CONTINUE
+C***********************************************************************
+      CALL VZ
+      RETURN
+      END
+      
+      SUBROUTINE VZ
 C==> Calculates the weight WEI due to FSI
       IMPLICIT REAL*8 (A-H,O-Z)
       COMMON/FSI_JR/JRAT
@@ -273,50 +473,28 @@ C==> Calculates the weight WEI due to FSI
       COMMON/FSI_C/C(10),AM,AMS,DM
       COMMON/FSI_COULPH/EIDC
       COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS
+
       COMPLEX*16 F,C,G,PSI12,PSI21
       COMPLEX*16 F12,F21
       COMPLEX*16 EIDC
       COMPLEX*8 Z8,CGAMMA
       COMMON/FSI_FFPN/FF12,FF21
       COMPLEX*16 FF12,FF21
-c--mlv------------------------------------------------------
-c      write(*,*)'------- WE are in VZ-----------'
-c      write(*,*)'JRAT=', JRAT
-c      write(*,*)'AK=', AK,'RP=',RP
-c      write(*,*)'X Y Z', X,Y,Z
-c       write(*,*)'PPX PPY PPZ',PPX,PPY,PPZ        
-c      write(*,*)'ETA',ETA
-c      write(*,*)'MSPIN',MSPIN
-c      write(*,*)'ISI',ISI
-c       write(*,*)'RP RPS',RP,RPS
-c      write(*,*)'AA- BEFORE FIRT------',AA
-c      write(*,*)'IQS',IQS
-c      write(*,*)'ICH',ICH
-c      write(*,*)'MSPIN',MSPIN
-c      write(*,*)'C',(C(I),I=1,MSPIN)
-c      write(*,*)'ACHR',ACHR
-c      write(*,*)'F12',F12,'F21',F21,'FF12',FF12,'FF21',FF21
-
-c--mlv------------------------------------------------------
       WEI=0.D0
       IF(JRAT.EQ.1)GOTO 11
       RHOS=AK*RP
       HS=X*PPX+Y*PPY+Z*PPZ
-c      write(*,*)'rhos',rhos,'hs',hs
       IF(RHOS.LT.15.D0.AND.RHOS+DABS(HS).LT.20.D0)GOTO 2
 C---Calc. EIDC=exp(i*Coul.Ph.);
 C-- used in calc. of hypergeom. f-s in SEQA, FAS at k*R > 15, 20
-c      write(*,*)'1111'
       Z8=CMPLX(1.,SNGL(ETA))
-      Z8=CGAMMA(Z8) 
+      Z8=CGAMMA(Z8)
       EIDC=Z8/CABS(Z8)
-c      write(*,*)'1111 z8 eidc',z8,eidc
+      
+c      write(*,*)'Z8 EIDC',Z8,EIDC
       
  2    CALL FF(RHOS,HS)
  11   MSP=MSPIN
-      
-c      write(*,*)'before go to 4', wei
       IF(ISI.EQ.0)GOTO 4  ! the strong interaction ON (OFF) if ISI=1(0)
       IF(RP.LT.AA)GOTO 4
       IF(JRAT.NE.1) CALL FIRT
@@ -329,10 +507,7 @@ c      write(*,*)'before go to 4', wei
       PSI12=FF12*(F12+G)
       PSI21=FF21*(F21+G)
       G=PSI12+JSIGN*PSI21
- 1    WEI=WEI+RHO(JJ)*(DREAL(G)**2+DIMAG(G)**2) 
-c      WRITE(*,*)'WEI111=',WEI
-
+ 1    WEI=WEI+RHO(JJ)*(DREAL(G)**2+DIMAG(G)**2)
       GOTO 8
  5    DO 6 JJ=1,MSP
       G=F(JJ)*C(JJ)
@@ -346,8 +521,6 @@ CW      WRITE(6,38)'JJ ',JJ,'F21+G ',DREAL(F21+G),DIMAG(F21+G)
       PSI12=FF12*(F12+G)
  6    WEI=WEI+RHO(JJ)*(DREAL(PSI12)**2+DIMAG(PSI12)**2)
       RETURN
-      
-c      write(*,*)'before 4 wei',wei
  4    PSI12=FF12*F12
       IF(IQS.EQ.0)GOTO 50 ! the quantum statistics ON (OFF) if IQS=1(0)
       PSI21=FF21*F21
@@ -355,20 +528,18 @@ c      write(*,*)'before 4 wei',wei
       DO 3 JJ=1,MSP
       JSIGN=-JSIGN
       G=PSI12+JSIGN*PSI21
-c      write(*,*)'msp jsign psi12 psi21',msp,jsign,psi12,psi21
  3    WEI=WEI+RHO(JJ)*(DREAL(G)**2+DIMAG(G)**2)
-c      write(*,*)'wei rho(1)G dreal(G) dimag(G)',wei,rho(1),G,
-c     + dreal(G),dimag(G)
       GOTO 8
-      
  50   WEI=DREAL(PSI12)**2+DIMAG(PSI12)**2
-c      write(*,*)'WEI ---in VZ 50 before returne', WEI
       RETURN
  8    WEI=WEI/2
-c      write(*,*)'WEI ---in VZ 8 before returne and END', WEI
+c      [M %/WRITE(*,*)WEI
+
       RETURN
       END
-      
+
+
       SUBROUTINE FIRT
 C---CALC. THE F(JJ)
 C-- F(JJ)*C(JJ)= DEVIATION OF THE BETHE-SALPETER AMPL. FROM PLANE WAVE
@@ -425,10 +596,8 @@ C---F= ASYMPTOTIC FORMULA (T= 0 NOT REQUIRED); NS= 2
       IF(AK.NE.0.D0)GOTO4
       HM=H*RPS
       IF(HM.GE.3.D15)GOTO6
-      
       FS1=DFRSIN(HM)
       FC1=DFRCOS(HM)
-      
       FC2=FC1
       FS2=FS1
       GOTO5
@@ -437,12 +606,10 @@ C---F= ASYMPTOTIC FORMULA (T= 0 NOT REQUIRED); NS= 2
       HM=H*(RSS-H1)**2
       HP=H*(RSS+H1)**2
       IF(HP.GE.3.D15)GOTO6
-      
       FS1=DFRSIN(HM)
       FC1=DFRCOS(HM)
       FS2=DFRSIN(HP)
       FC2=DFRCOS(HP)
-      
       GOTO 5
   11  CONTINUE
       FS1=0.D0
@@ -456,10 +623,8 @@ C---F= ASYMPTOTIC FORMULA (T= 0 NOT REQUIRED); NS= 2
       IF(AK.NE.0.D0)GOTO 12
       HM=H*RPS
       IF(HM.GE.3.D15)GOTO6
-      
       FS1=FS1+DFRSIN(HM)/2
       FC1=FC1+DFRCOS(HM)/2
-      
       IF(I.EQ.1)GOTO 13
       FC2=FC1
       FS2=FS1
@@ -469,12 +634,10 @@ C---F= ASYMPTOTIC FORMULA (T= 0 NOT REQUIRED); NS= 2
       HM=H*(RSS-H1)**2
       HP=H*(RSS+H1)**2
       IF(HP.GE.3.D15)GOTO6
-      
       FS1=FS1+DFRSIN(HM)/2
       FC1=FC1+DFRCOS(HM)/2
       FS2=FS2+DFRSIN(HP)/2
       FC2=FC2+DFRCOS(HP)/2
-      
   13  CONTINUE
   5   C12=FC1+FS2
       S12=FS1+FC2
@@ -488,7 +651,6 @@ C---F= ASYMPTOTIC FORMULA (T= 0 NOT REQUIRED); NS= 2
  10   CONTINUE
       RETURN
       END
-      
       FUNCTION EXF(X)
       IMPLICIT REAL*8 (A-H,O-Z)
       IF(X.LT.-15.D0) GO TO 1
@@ -497,7 +659,6 @@ C---F= ASYMPTOTIC FORMULA (T= 0 NOT REQUIRED); NS= 2
   1   EXF=.0D0
       RETURN
       END
-      
       SUBROUTINE SEQ(X,H)
 C---CALC. FUNCTIONS B, P (EQS. (17) OF G-K-L-L);
 C-- NEEDED TO CALC. THE CONFLUENT HYPERGEOMETRIC FUNCTION GST.
@@ -526,7 +687,6 @@ C-- NEEDED TO CALC. THE CONFLUENT HYPERGEOMETRIC FUNCTION GST.
       PH(2)=PH(3)
       GOTO 2
       END
-      
       SUBROUTINE SEQA(X,H)
 C---CALC. FUNCTIONS CHH=REAL(GST), SH=IMAG(GST)/ACH, B=SH/H
 C-- IN THE ASYMPTOTIC REGION H=K*R >> 1.
@@ -545,7 +705,6 @@ C-- IN THE ASYMPTOTIC REGION H=K*R >> 1.
       B=SH/H
       RETURN
       END
-      
       SUBROUTINE FF(RHO,H)
 C---Calc. F12, F21;
 C-- F12= FF0* plane wave,  FF0=F*ACHR,
@@ -568,7 +727,6 @@ C-- ACHR=sqrt(ACH), where ACH is the Coulomb factor
       F21=FF0(RHO,-H)*F21
       RETURN
       END
-      
       FUNCTION FAS(RKS)
 C-- FAS=F*ACHR
 C---F is the confluent hypergeometric function at k*r >> 1
@@ -585,7 +743,6 @@ C-- ACHR=sqrt(ACH), where ACH is the Coulomb factor
       FAS=FAS-DCMPLX(DCOS(RKS),DSIN(RKS))*ETA/RKS/ZZ1
       RETURN
       END
-      
       FUNCTION FF0(RHO,H)
 C-- FF0=F*ACHR
 C-- F is the confluent hypergeometric function
@@ -618,7 +775,6 @@ CC      GOTO 5 ! rejects the approx. calcul. of hyperg. f-ion F
       FF0=FF0*ACHR
       RETURN
       END
-      
       FUNCTION HC(XA)
 C---HC = h-function of Landau-Lifshitz: h(x)=Re[psi(1-i/x)]+ln(x)
 C-- psi(x) is the digamma function (the logarithmic derivative of
@@ -654,7 +810,6 @@ CC      RETURN
    4  XP=X2*XP
       RETURN
       END
-      
       FUNCTION ACP(X)
 C--- ACP = COULOMB PENETRATION FACTOR
       IMPLICIT REAL*8 (A-H,O-Z)
@@ -665,7 +820,6 @@ C--- ACP = COULOMB PENETRATION FACTOR
    1  ACP=1.D-6
       RETURN
       END
-      
       SUBROUTINE GST(X,H)
 C---CALC. THE CONFL. HYPERGEOM. F-N = CHH+i*SH
 C-- AND THE COULOMB F-S B, P (CALLS SEQ OR SEQA).
@@ -691,7 +845,7 @@ C-- AND THE COULOMB F-S B, P (CALLS SEQ OR SEQA).
   4   CALL SEQA(X,H)
       RETURN
       END
-      
+
       FUNCTION FF1(RHO,H)
 C---FF1=FF0; used for particle-nucleus system
 C-- FF0=F12*ACHR
@@ -712,19 +866,15 @@ C-- ACHR=sqrt(ACH), where ACH is the Coulomb factor
       IF(RHO.LT.15.D0.AND.RHO+H.LT.20.D0)GOTO 2
 C---Calc. EIDC=exp(i*Coul.Ph.);
 C-- used in calc. of hypergeom. f-s in SEQA, FAS at k*R > 15, 20
-
       Z8=CMPLX(1.,SNGL(ETA))
       Z8=CGAMMA(Z8)
-      
       EIDC=Z8/CABS(Z8)
-      
  2    FF1=FF0(RHO,H)
       RETURN
       END
  
       FUNCTION G(AK)
 C---Used to calculate SW scattering amplitude for alpa-alpha system
-C-- and for sear.f (square well potential search)
 C---NOTE THAT SCATT. AMPL.= 1/CMPLX(G(AK),-AK*ACH)
       IMPLICIT REAL*8 (A-H,O-Z)
       COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10),
@@ -732,7 +882,6 @@ C---NOTE THAT SCATT. AMPL.= 1/CMPLX(G(AK),-AK*ACH)
       COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
       COMMON/FSI_ACH/HPR,AC,ACH,ACHR,JJ,MSPIN
       COMMON/FSI_BP/B,P/FSI_DERIV/BPR,PPR/FSI_SHH/SH,CHH
-      COMMON/FSI_DAK/DAK,HAC,IFUN
       HAC=.0D0
       IF(ICH.EQ.0)GOTO 1
       XH=AC*AK
@@ -770,7 +919,6 @@ C   GG= [B'(kp,d)*ReG(k,d)-B(kp,d)*ZZ]/d
 C    G= GG/[B(kp,d)*B'(k,d)-B'(kp,d)*B(k,d)]
       RETURN
       END
-      
       SUBROUTINE DERIW(X,H)
 C---CALLED BY F-N G(AK)
       IMPLICIT REAL*8 (A-H,O-Z)
@@ -793,15 +941,22 @@ C---CALLED BY F-N G(AK)
       PPR=PPR+X*(P-Q1)/HHH
       RETURN
       END
+C================================================================
+
 
 c-------------#include "gen/pilot.h"
       FUNCTION DFRSIN(X)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)      
 
       DIMENSION A(0:16),B(0:15),C1(0:25),C2(0:28)
+
       PARAMETER (Z1 = 1, R8 = Z1/8, R32 = Z1/32)
+
       DATA C0 /1.25331 41373 15500 3D0/
+
       DATA NA,NB,NC1,NC2 /16,15,25,28/
+
       DATA A( 0) / 0.76435 13866 41860 002D0/
       DATA A( 1) /-0.43135 54754 76601 793D0/
       DATA A( 2) / 0.43288 19997 97266 531D0/
@@ -966,14 +1121,17 @@ c-------------#include "gen/pilot.h"
       END
 
 c--#include "gen/pilot.h"
-      FUNCTION CGAMMA(Z) 
+      FUNCTION CGAMMA(Z)
       
       COMPLEX*8 CGAMMA
       COMPLEX*8 Z,U,V,F,H,S       
       CHARACTER NAME*(*)
       CHARACTER*80 ERRTXT      
       PARAMETER (NAME = 'CGAMMA')
+
       DIMENSION C(0:15)
+
       PARAMETER (Z1 = 1, HF = Z1/2)
 
 c--#if defined(CERNLIB_QF2C)
@@ -1048,202 +1206,3 @@ c----#endif
   101 FORMAT('ARGUMENT EQUALS NON-POSITIVE INTEGER = ',1P,E15.1)
       END
 
-
-      SUBROUTINE FSIPN !(WEIF)
-C  calculating particle-nucleus Coulomb Wave functions FFij
-      IMPLICIT REAL*8 (A-H,O-Z)
-      COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
-      COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1,  !part. momenta in NRF
-     1               P2X,P2Y,P2Z,E2,P2
-      COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emis. points in NRF
-     1                X2,Y2,Z2,T2,R2
-      COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
-      COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
-      COMMON/FSI_ICH1/ICH1
-      COMMON/FSI_ETA/ETA
-      COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS
-      COMMON/FSI_FFPN/FF12,FF21
-      COMPLEX*16 FF1,FF12,FF21
-      FF12=DCMPLX(1.D0,0.D0)
-      FF21=DCMPLX(1.D0,0.D0)
-      IF(I3C.EQ.0)RETURN
-      ICH1=IDINT(C1)
-      IF(ICH1.EQ.0)GOTO 11
-      XH=AC1*P1
-      ACH=ACP(XH)
-      ACHR=DSQRT(ACH)
-      ETA=0.D0
-      IF(XH.NE.0.D0)ETA=1/XH
-      RHOS=P1*R1
-      HS=X1*P1X+Y1*P1Y+Z1*P1Z
-      FF12=FF12*FF1(RHOS,HS)
-      IF(IQS.EQ.0)GOTO 11
-      RHOS=P1*R2
-      HS=X2*P1X+Y2*P1Y+Z2*P1Z
-      FF21=FF21*FF1(RHOS,HS)
- 11   ICH1=IDINT(C2)
-      IF(ICH1.EQ.0)GOTO 10
-      XH=AC2*P2
-      ACH=ACP(XH)
-      ACHR=DSQRT(ACH)
-      ETA=0.D0
-      IF(XH.NE.0.D0)ETA=1/XH
-      RHOS=P2*R2
-      HS=X2*P2X+Y2*P2Y+Z2*P2Z
-      FF12=FF12*FF1(RHOS,HS)
-CW      WRITE(6,41)'AC2 ',AC2,'ACH ',ACH,'ETA ',ETA,'RHOS ',RHOS,'HS ',HS
-41    FORMAT(5(A5,E11.4))
-CW      WRITE(6,40)'FF12 ',DREAL(FF12),DIMAG(FF12)
-      IF(IQS.EQ.0)GOTO 10
-      RHOS=P2*R1
-      HS=X1*P2X+Y1*P2Y+Z1*P2Z
-      FF21=FF21*FF1(RHOS,HS)
-CW      WRITE(6,41)'AC1 ',AC1,'ACH ',ACH,'ETA ',ETA,'RHOS ',RHOS,'HS ',HS
-CW      WRITE(6,40)'FF21 ',DREAL(FF21),DIMAG(FF21)
-40    FORMAT(A7,2E12.4)
- 10   CONTINUE
-C  WEIF = the weight due to the Coulomb particle-nucleus interaction
-      WEIF=DREAL(FF12)**2+DIMAG(FF12)**2
-      IF(IQS.EQ.1)WEIF=0.5D0*(WEIF+DREAL(FF21)**2+DIMAG(FF21)**2)
-      RETURN
-      END
-C=======================================================
-C
-      SUBROUTINE FSIWF !(WEI)
-C==> Prepares necessary quantities and call VZ(WEI) to calculate
-C    the weight due to FSI
-      IMPLICIT REAL*8 (A-H,O-Z)
-      COMMON/FSI_CVK/V,CVK
-      COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1,  !part. momenta in NRF
-     1               P2X,P2Y,P2Z,E2,P2
-      COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
-     1               X,Y,Z,T,RP,RPS
-      COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emis. points in NRF
-     1                X2,Y2,Z2,T2,R2
-      COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
-      COMMON/FSI_SPIN/RHO(10)
-      COMMON/FSI_BP/B,P
-      COMMON/FSI_ETA/ETA
-      COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
-      COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10),
-     1              SBKRB(10),SDKK(10)
-      COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
-      COMMON/FSI_RR/F(10)
-      COMMON/FSI_FD/FD(10),RD(10)
-      COMMON/FSI_C/C(10),AM,AMS,DM
-      COMPLEX*16 C,F
-      COMMON/FSI_AA/AA
-      COMMON/FSI_SHH/SH,CHH
-      COMMON/FSI_AAPI/AAPI(20,2)/FSI_AAND/AAND(20,4)
-      COMMON/FSI_AAPIN/AAPIN(20,2)
-      COMMON/FSI_P12/P12X,P12Y,P12Z,E12,P12,AM12,EPM
-      COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS
-      
-      
-cmlv
-        IF(IRANPOS.EQ.0)THEN
-C==>calculating relative 4-coordinates of the particles in PRF
-C-  {T,X,Y,Z} from the relative coordinates {TS,XS,YS,ZS} in NRF
-        XS=X1-X2
-        YS=Y1-Y2
-        ZS=Z1-Z2
-        TS=T1-T2      
-        RS12=XS*P12X+YS*P12Y+ZS*P12Z
-        H1=(RS12/EPM-TS)/AM12
-        X=XS+P12X*H1
-        Y=YS+P12Y*H1
-        Z=ZS+P12Z*H1
-        T=(E12*TS-RS12)/AM12
-        RPS=X*X+Y*Y+Z*Z
-        RP=DSQRT(RPS)            
-c        WRITE(6,38)'RP ',RP,'X ',X,Y,Z,T
-c        FORMAT(A7,E11.4,A7,4E11.4)
-        ENDIF 
-
-      CVK=(P12X*PPX+P12Y*PPY+P12Z*PPZ)/(P12*AK)
-      V=P12/E12
-      IF(ICH.EQ.0)GOTO 21
-      XH=AC*AK
-      ACH=ACP(XH)
-      ACHR=DSQRT(ACH)
-      ETA=0.D0
-      IF(XH.NE.0.D0)ETA=1/XH
-C---HCP, HPR needed (e.g. in GST) if ICH=1
-      HCP=HC(XH)
-      HPR=HCP+.1544313298D0
-  21  CONTINUE
-      MSP=MSPIN
-      DO 30 JJ=1,MSP
-      ISPIN=JJ
-      IF(NS.NE.1)GOTO22
-C---Calc. quantities for the square well potential;
-C-- for LL > 5 the square well potential is not possible or available
-      IF(LL.EQ.4)GOTO 22
-      BK(JJ)=DSQRT(EB(JJ)**2+AKS)
-      XRA=2*RB(JJ)/AC
-      HRA=BK(JJ)*RB(JJ)
-      CALL SEQ(XRA,HRA)
-      SBKRB(JJ)=HRA*B
-      HRA=AK*RB(JJ)
-      CALL GST(XRA,HRA)
-      SDK(JJ)=SH
-      CDK(JJ)=CHH
-      SDKK(JJ)=RB(JJ)
-      IF(AK.NE.0.D0)SDKK(JJ)=SH/AK
-      IF(ICH.EQ.1)SDK(JJ)=ACH*SDK(JJ)
-  22  CONTINUE
-C-----------------------------------------------------------------------
-C---Calc. the strong s-wave scattering amplitude = C(JJ)
-C-- divided by Coulomb penetration factor squared (if ICH=1)
-      IF(NS.NE.1)GOTO 230
-      IF(LL.NE.4)GOTO 230 ! SW scat. amplitude used for alfa-alfa only
-      GAK=G(AK)
-      AKACH=AK
-      IF(ICH.EQ.1)AKACH=AK*ACH
-      C(JJ)=1/DCMPLX(GAK,-AKACH) ! amplitude for the SW-potential
-      GOTO 30
- 230  IF(LL.EQ.5.OR.LL.EQ.6.OR.LL.EQ.7)GOTO20    ! pipi  
-      IF(LL.EQ.12.OR.LL.EQ.13)GOTO20             ! piN 
-      IF(LL.EQ.8.OR.LL.EQ.9.OR.LL.EQ.18)GOTO20   ! Nd, dd
-      IF(LL.EQ.14.OR.LL.EQ.17.OR.LL.EQ.23)GOTO27 ! K+K-, K-p, K0K0-b
-      A1=RD(JJ)*FD(JJ)*AKS
-      A2=1+.5D0*A1
-      IF(ICH.EQ.1)A2=A2-2*HCP*FD(JJ)/AC
-      AKF=AK*FD(JJ)
-      IF(ICH.EQ.1)AKF=AKF*ACH
-      C(JJ)=FD(JJ)/DCMPLX(A2,-AKF)
-      GOTO30
- 20   CONTINUE
-C---Calc. scatt. ampl. C(JJ) for pipi, piN and Nd, dd
-      JH=LL-7+2*JJ-2
-      IF(LL.EQ.8.OR.LL.EQ.9)GPI2=GND(AKS,JH)
-      IF(LL.EQ.18)GPI2=GDD(AKS,JJ)
-      IF(LL.EQ.5.OR.LL.EQ.6.OR.LL.EQ.7)GPI2=GPIPI(AKS,2)
-      IF(LL.EQ.12.OR.LL.EQ.13)GPI2=GPIN(AKS,2)
-      C(JJ)=1.D0/DCMPLX(GPI2,-AK) !pi+pi+, nd, pd, pi+p, dd
-      IF(LL.NE.5.AND.LL.NE.6.AND.LL.NE.13)GOTO27
-      IF(LL.EQ.5.OR.LL.EQ.6)GPI1=GPIPI(AKS,1)
-      IF(LL.EQ.13)GPI1=GPIN(AKS,1)
-      IF(LL.EQ.5.OR.LL.EQ.13)
-     c           C(JJ)=.6667D0/DCMPLX(GPI1,-AK)+.3333D0*C(JJ) !pi+pi-,pi-p
-      IF(LL.EQ.6)C(JJ)=.3333D0/DCMPLX(GPI1,-AK)+.6667D0*C(JJ) !pi0pi0
- 27   CONTINUE
-C---Calc. K+K-, K0K0-b or K-p s-wave scatt. ampl.
-      IF(LL.EQ.14.OR.LL.EQ.23)CALL CKKB
-      IF(LL.EQ.17)C(JJ)=DCMPLX(3.29D0,3.55D0)
-C---Calc. pi+pi-, pi+pi+, pd, pi+p, pi-p, K+K- or K-p s-wave scatt. ampl.
-C-- divided by Coulomb penetration factor squared (if ICH=1)
-      IF(ICH.EQ.0)GOTO 30
-      AAK=ACH*AK
-      HCP2=2*HCP/AC
-      C(JJ)=1/(1/C(JJ)-HCP2+DCMPLX(0.D0,AK-AAK))
-c      write(*,*)'c(jj)',c(jj)
- 30   CONTINUE
-C***********************************************************************
-c      write(*,*)'before call vz in fsiwf wei',wei
-      CALL VZ  !(WEI)
-      RETURN
-      END
index d071e2e..de694d7 100644 (file)
@@ -1,58 +1,84 @@
-      SUBROUTINE ltran12
-C==>calculating particle momentum in PRF {EE,PPX,PPY,PPZ} from
-C-  the momentum of the first particle {E1,P1X,P1Y,P1Z) in NRF
+
+       SUBROUTINE ltran12
+C
+
+c      SUBROUTINE TRANS(icrf,irot)
+C==>  TRANSformation to the Co-moving frame  (icrf>0) and
+C     Rotation to the system where (Pt || X),(irot=1).
+C-FSI ***************************************************
       IMPLICIT REAL*8 (A-H,O-Z)
-      COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1,  !part. momenta in NRF
+      COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1,  ! momenta in NRF
      1               P2X,P2Y,P2Z,E2,P2
-      COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
+      COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emis. 
+     1                X2,Y2,Z2,T2,R2  ! points in NRF
+      COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS, ! momenta in PRF
      1               X,Y,Z,T,RP,RPS
       COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
       COMMON/FSI_P12/P12X,P12Y,P12Z,E12,P12,AM12,EPM
-      COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, !4-coord. of emis. points in NRF
-     1                X2,Y2,Z2,T2,R2
-      COMMON/FSI_CONS/PI,PI2,SPI,DR,W
-      COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS
+      COMMON/FSI_CVK/V,CVK
+C-FSI ***************************************************
+      COMMON /PAIR/P12T,V12Z,GAMZ,V12T,CPHI,SPHI
       
-C   fm --> 1/GeV
-c      write(*,*)'in LTRAN12 W', W
-c      write(*,*)'in LTRAN12 p1x p1y p1z', P1X,P1Y,P1Z
-c      write(*,*)'in LTRAN12 p2x p2y p2z', P2X,P2Y,P2Z
-c       write(*,*)'in LTRAN12 x y z t', X,Y,Z,T
-c      write(*,*)'in LTRAN12 x2 y2 z2', X2,Y2,Z2,T2
-
-CMLV
-       IF(IRANPOS.EQ.0)THEN      
-      X1=X1*W
-      Y1=Y1*W
-      Z1=Z1*W
-      T1=T1*W
-      X2=X2*W
-      Y2=Y2*W
-      Z2=Z2*W
-      T2=T2*W
-C   calculating Ri, Pi and Ei
-      R1=DSQRT(X1*X1+Y1*Y1+Z1*Z1)
-      R2=DSQRT(X2*X2+Y2*Y2+Z2*Z2)
-      ENDIF
-CMLV
-
+      icrf=1
+      irot=1
+      
+C---> Particle energies ---------
       P1S=P1X*P1X+P1Y*P1Y+P1Z*P1Z
       P2S=P2X*P2X+P2Y*P2Y+P2Z*P2Z
-      P1=DSQRT(P1S)
-      P2=DSQRT(P2S)
       E1=DSQRT(AM1*AM1+P1S)
       E2=DSQRT(AM2*AM2+P2S)
-C-----------------------------------------------------------------------
+C---> Pair parameters -----------
+      E12=E1+E2       ! Energy
+      P12X=P1X+P2X    ! Px
+      P12Y=P1Y+P2Y    ! Py
+      P12Z=P1Z+P2Z    ! Pz
+      P12S=P12X**2+P12Y**2+P12Z**2
+      P12 =DSQRT(P12S)! Momentum
+      V12 =P12/E12    ! Velocity
+      CTH =P12Z/P12   ! cos(theta)
+      STH =DSQRT(1.D0-CTH**2) !sin
+      V12Z=V12*CTH    ! Longit. V
+      GAMZ=1.D0/DSQRT(1.D0-V12Z**2)
+C--      V12T=V12*STH    ! Transv. V in CMS (not needed)
+      P12TS=P12X*P12X+P12Y*P12Y
+      P12T=DSQRT(P12TS) !Pt
+C===> Azimuthal rotation (Pt||X) ============
+      IF(V12T.NE.0.D0) THEN
+        CPHI=P12X/P12T  ! cos(phi)
+        SPHI=P12Y/P12T  ! sin(phi)
+        IF((irot.eq.1)) THEN 
+          CALL ROT8(P1X,P1Y,SPHI,CPHI,P1X,P1Y)       
+          CALL ROT8(P2X,P2Y,SPHI,CPHI,P2X,P2Y)
+          CALL ROT8(X1,Y1,SPHI,CPHI,X1,Y1)       
+          CALL ROT8(X2,Y2,SPHI,CPHI,X2,Y2) 
+        END IF             
+      ELSE ! Rotation impossible 
+       CPHI=2.D0 ! to avoid
+       SPHI=2.D0 ! using it !  
+      END IF             
+C===> Co-moving ref. frame       ============
+      IF(icrf.gt.0) THEN
+        CALL LTR8(P1Z,E1,V12Z,GAMZ,P1Z,E1a)
+        CALL LTR8(P2Z,E2,V12Z,GAMZ,P2Z,E2a)
+        P1S=P1X*P1X+P1Y*P1Y+P1Z*P1Z
+        P2S=P2X*P2X+P2Y*P2Y+P2Z*P2Z
+        E1=DSQRT(AM1*AM1+P1S)
+        E2=DSQRT(AM2*AM2+P2S)
+        CALL LTR8(Z1,T1,V12Z,GAMZ,Z1,T1)
+        CALL LTR8(Z2,T2,V12Z,GAMZ,Z2,T2)
+      END IF        
+C===> Pair reference frame       ============
+      P1=DSQRT(P1S)
+      P2=DSQRT(P2S)
       E12=E1+E2
       P12X=P1X+P2X
       P12Y=P1Y+P2Y
       P12Z=P1Z+P2Z
       P12S=P12X**2+P12Y**2+P12Z**2
-      AM12=DSQRT(E12**2-P12S)
+      P12 =DSQRT(P12S)
+      AM12S=E12*E12-P12S
+      AM12=DSQRT(AM12S)
       EPM=E12+AM12
-      P12=DSQRT(P12S)
       P112=P1X*P12X+P1Y*P12Y+P1Z*P12Z
       H1=(P112/EPM-E1)/AM12
       PPX=P1X+P12X*H1
@@ -61,8 +87,71 @@ C-----------------------------------------------------------------------
       EE=(E12*E1-P112)/AM12
       AKS=EE**2-AM1**2
       AK=DSQRT(AKS)
-CW      WRITE(6,38)'AK ',AK,'K ',PPX,PPY,PPZ,EE
-38    FORMAT(A7,E11.4,A7,4E11.4)
+      CVK=(P12X*PPX+P12Y*PPY+P12Z*PPZ)/(P12*AK)
+      V=P12/E12 
+      V12T=P12T/SQRT(AM12S+P12TS) ! transverse velocity in LCMS   
+C---> Coordinates -----------------------------
+      XS=X1-X2
+      YS=Y1-Y2
+      ZS=Z1-Z2
+      TS=T1-T2  
+      RS12=XS*P12X+YS*P12Y+ZS*P12Z
+      H1=(RS12/EPM-TS)/AM12
+
+cmlv      X=XS+P12X*H1
+cmlv      Y=YS+P12Y*H1
+cmlv      Z=ZS+P12Z*H1
+cmlv      T=(E12*TS-RS12)/AM12
+cmlv      RPS=X*X+Y*Y+Z*Z
+cmlv      RP=DSQRT(RPS)
+
+
+      RETURN
+      END
+C====       
+      SUBROUTINE LTR8(Z,T,BETA,GAMMA,ZT,TT)
+C===> Lorentz Transf. of Z(Pz) and T(E) to moving ref. frame.(REAL*8) 
+CInp: Z,T-Zcoord,Time before tr., BETA,GAMMA- velocity, Lor.fact.
+COut: ZT,TT- "    "   after transformation.
+C==== ===============================================================
+      IMPLICIT REAL*8 (A-H,O-Z)
+      ZH=GAMMA*(Z-BETA*T)     
+      TT=GAMMA*(T-BETA*Z)
+      ZT=ZH
       RETURN
       END
+C====       
+      SUBROUTINE LTR4(Z,T,BETA,GAMMA,ZT,TT)
+C===> Lorentz Transf. of Z(Pz) and T(E) to moving ref. frame.(real*4) 
+CInp: Z,T-Zcoord,Time before tr., BETA,GAMMA- velocity, Lor.fact.
+COut: ZT,TT- "    "   after transformation.
+C==== ===============================================================
+      ZH=GAMMA*(Z-BETA*T)     
+      TT=GAMMA*(T-BETA*Z)
+      ZT=ZH
+      RETURN
+      END
+C====           
+      SUBROUTINE ROT8(X,Y,SF,CF,XR,YR)
+C===> Rotation with the angle f. (REAL*8)
+CInp: X,Y-coord. before rotation; SF=sin(f), CF=cos(f),
+COut: XR,YR - coordinates after rotation.
+C==== =================================================
+      IMPLICIT REAL*8 (A-H,O-Z)
+      XH=X*CF+Y*SF !Y              
+      YR=Y*CF-X*SF !   _-X'   
+      XR=XH        ! _- f      
+      RETURN       !------>
+      END          !     X
+
+
+
+
+
+
+
+
+
+
+
  
\ No newline at end of file