C************************************************************
C*****                                                  *****
C*****                   MAIN PROGRAM                   *****
C*****                                                  *****
C*****     SCATTERING WAVE FUNCTION CALCULATED WITH     *****
C*****              EIKONAL APPROXIMATION               *****
C*****                                                  *****
C************************************************************
C
C=======================================================================
      PROGRAM WFEIKONAL02
      IMPLICIT REAL*8(A-H,O-Z)
C=======================================================================
C
      PARAMETER (HC0=197.3269601D0,AC0=931.494013D0,EC0=137.03599976D0)
      PARAMETER(NBDIM=1001)
C
      COMPLEX*16 SEK(NBDIM)
C
      COMMON /PHYSCONST/HC,AC,EC
      COMMON /MACHINE/EPSLIM,EXPMAX
CCC---------------------------------------------------------------------
      HC=HC0
      AC=AC0
      EC=EC0
C
C-------------------------------------
C---  LIMIT VALUE FOR COMPUTATION  ---
C-------------------------------------
C
      EXPMAX=300.0D0
      EPSLIM=1.0D-14
C
C
C-------------------
C---  FILE OPEN  ---
C-------------------
C
      CALL FOPEN
C
C
C-------------------------
C---  PARAMETER INPUT  ---
C-------------------------
C
      CALL INPUT(AP,AT,ZP,ZT,ELAB,RMAX,DR,NRMAX,V0,RV,AV
     1          ,W0,RW,AW,WD,RD,AD,BMIN,BMAX,DB,NBMAX,KIBCS
     2          ,THMIN,THMAX,DTH,NTHMAX,ICM,IREL)
C
C
C-----------------------------------
C---  SCATTERING WAVE FUNCTION  ---
C-----------------------------------
C
      CALL WF(AP,AT,ZP,ZT,ELAB,RMAX,DR,NRMAX,V0,RV,AV
     1       ,W0,RW,AW,WD,RD,AD,BMIN,BMAX,DB,NBMAX,NBDIM,SEK,FK
     2       ,ICM,IREL)
C
C
C-----------------------
C---  CROSS SECTION  ---
C-----------------------
C
      CALL CS(FK,SEK,BMIN,BMAX,DB,NBMAX,NBDIM,KIBCS
     1       ,THMIN,THMAX,DTH,NTHMAX)
C
      STOP 0
      END
C
C
C
Cver0.2
C=======================================================================
      SUBROUTINE CS(FK,SEK,BMIN,BMAX,DB,NBMAX,NBDIM,KIBCS
     1             ,THMIN,THMAX,DTH,NTHMAX)
      IMPLICIT REAL*8 (A-H,O-Z)
C=======================================================================
C
      COMPLEX*16 CSUM,SEK(NBDIM)
CCC---------------------------------------------------------------------
C
      PI=ACOS(-1.0D0)
      PIAD=PI/180.0D0
C
      WRITE(KIBCS,*) ' THETA    CS[mb/sr]'
C
      DO ITH=1,NTHMAX
       THT=(ITH-1)*DTH+THMIN
       RAD=THT*PIAD
C
       Q=2.0D0*FK*SIN(RAD/2.0D0)
C
       CSUM=0.0D0
       SUME=0.0D0
       SUMR=0.0D0
       SUMT=0.0D0
       DO IB=1,NBMAX
        B=(IB-1)*DB+BMIN
        ARG=B*FK*RAD
        CSUM=CSUM+BESSJ0(ARG)*(SEK(IB)-1.0D0)*B
        SUME=SUME+ABS(1.0D0-SEK(IB))**2*B
        SUMR=SUMR+(1.0D0-ABS(SEK(IB))**2)*B
        SUMT=SUMT+2.0D0*(1.0D0-DBLE(SEK(IB)))*B
       END DO
       CSUM=CSUM*DB
C
       SIG=FK**2*ABS(CSUM)**2*10.0D0
C
       WRITE(KIBCS,600) THT,SIG
  600  FORMAT(F8.4,E13.5)
C
      END DO
C
      SIGE=SUME*2.0D0*PI*DB*10.0D0
      SIGR=SUMR*2.0D0*PI*DB*10.0D0
      SIGT=SUMT*2.0D0*PI*DB*10.0D0
C
C-----OUTPUT
C
      WRITE(6,601) SIGE,SIGR,SIGT
  601 FORMAT(/3X,'-- INTEGRATED CROSS SECTION --'
     1       /8X,'NUCLEAR ELASTIC      :',F10.3,' [mb]'
     2       /8X,'REACTION             :',F10.3,' [mb]'
     3       /8X,'TOTAL                :',F10.3,' [mb]')
C
      RETURN
C
      END
Cver0.2
C
C
C
C=======================================================================
      SUBROUTINE WF(AP,AT,ZP,ZT,ELAB,RMAX,DR,NRMAX,V0,RV,AV
     1             ,W0,RW,AW,WD,RD,AD,BMIN,BMAX,DB,NBMAX,NBDIM,SEK,FK
     2             ,ICM,IREL)
      IMPLICIT REAL*8 (A-H,O-Z)
C=======================================================================
C
      COMPLEX*16 CSUM,WI
      COMPLEX*16 SEK(NBDIM)
C
      COMMON /PHYSCONST/HC,AC,EC
      COMMON /MACHINE/EPSLIM,EXPMAX
CCC---------------------------------------------------------------------
C
C-----DIMENSION CHECK FOR B
C
      IF(NBMAX.GT.NBDIM) GO TO 910
C
C-----INITIALIZATION
C
      SEK(:)=0.0D0
C
      WI=(0.0D0,1.0D0)
      RMAX=(NRMAX-1)*DR
      FM=AP*AC
      FK=SQRT(2.0D0*FM*ELAB)/HC
      HVEL=FK*HC**2/FM
C
C-----KINEMATICAL CORRECTION
C
      IF(IREL.EQ.1) THEN
        FMP=AP*AC
        FMT=AT*AC
        EPL=ELAB+FMP
        ETL=FMT
        FKPL=SQRT(ELAB*(ELAB+2.0D0*FMP))/HC
        IF(ICM.EQ.0) THEN
          FK=FKPL
          HVEL=FKPL*HC**2/EPL            ! invariant
         ELSE
          FK=FKPL*FMT/SQRT(FMP**2+FMT**2+2.0D0*FMT*EPL)
          HVEL=FKPL*HC**2/EPL            ! invariant
        END IF
       ELSE IF(ICM.EQ.1) THEN
        ECM=ELAB*AT/(AP+AT)
        FM=AP*AT/(AP+AT)*AC
        FK=SQRT(2.0D0*FM*ECM)/HC
        HVEL=FK*HC**2/FM
      END IF
C
C-----OUTPUT
C
      WRITE(6,600) FK,HVEL
  600 FORMAT(/3X,'-- KINEMATICS --'
     1       /8X,'FK (WAVE NUMBER)     :',F10.4,' [1/fm]'
     2       /8X,'HBAR*C               :',F10.4,' [MeV fm]')
C
      WRITE(6,601)
  601 FORMAT(/3X,'-- S-MATRIX --'
     1       /7X,'b      realS        imagS')
C
C
      DO IB=1,NBMAX
       B=(IB-1)*DB+BMIN
       IF(B.GT.RMAX) THEN
         NBMAX=IB-1
         RETURN
       END IF
C
       DZ=DR
       ZMAX=SQRT(RMAX**2-B**2)
       NZMAX0=ZMAX/DZ+1.01
       NZMAX=2*(NZMAX0-1)+1
       ZMAX=(NZMAX0-1)*DZ
C
       CSUM=0.0D0
       DO IZ=1,NZMAX
        Z=(IZ-1)*DZ-ZMAX
        R=SQRT(Z**2+B**2)
        ARGV=(R-RV)/AV
        IF(ARGV.GT.EXPMAX) THEN
          VOPT=0.0D0
         ELSE
          VOPT=-V0/(1.0D0+EXP(ARGV))
        END IF
        ARGW=(R-RW)/AW
        IF(ARGW.GT.EXPMAX) THEN
          WOPTV=0.0D0
         ELSE
          WOPTV=-W0/(1.0D0+EXP(ARGW))
        END IF
        ARGD=(R-RD)/AD
        IF(2.0D0*ARGD.GT.EXPMAX) THEN
          WOPTS=0.0D0
         ELSE
          WOPTS=-4.0D0*WD*EXP(ARGD)/(1.0D0+EXP(ARGD))**2
        END IF
        WOPT=WOPTV+WOPTS
        CSUM=CSUM+VOPT+WI*WOPT
       END DO
C
       SEK(IB)=EXP(CSUM*DZ/(HVEL*WI))
C
       WRITE(6,602) B,SEK(IB)
  602  FORMAT(F10.4,2E13.5)
C
      END DO
C
      RETURN
C
  910 WRITE(6,6910) NBMAX,NBDIM
 6910 FORMAT(/3X,'-- DIMENSION ERROR IN WF --'
     1       /8X,'NBMAX > NBDIM:  NBMAX,NBDIM=',2I8)
      STOP
C
      END
C
C
C
C=======================================================================
      SUBROUTINE INPUT(AP,AT,ZP,ZT,ELAB,RMAX,DR,NRMAX,V0,RV,AV
     1                ,W0,RW,AW,WD,RD,AD,BMIN,BMAX,DB,NBMAX,KIBCS
     2                ,THMIN,THMAX,DTH,NTHMAX,ICM,IREL)
      IMPLICIT REAL*8 (A-H,O-Z)
C=======================================================================
C
      CHARACTER*2 NAMEP,NAMET
      CHARACTER*50 COMMENT
C
      COMMON /MACHINE/EPSLIM,EXPMAX
CCC---------------------------------------------------------------------
C
      KIBAN=5
C
      AV=0.6D0
      AW=0.6D0
      AD=0.6D0
C
      READ(KIBAN,500) COMMENT
  500 FORMAT(A50)
      READ(KIBAN,501) AP,ZP,AT,ZT
  501 FORMAT(5F10.0)
      READ(KIBAN,501) ELAB,RMAX,DR
      READ(KIBAN,501) BMIN,BMAX,DB
      READ(KIBAN,501) FICM,FIREL
      READ(KIBAN,501) V0,RV,AV
      READ(KIBAN,501) W0,RW,AW
      READ(KIBAN,501) WD,RD,AD
      READ(KIBAN,500) DUMMY
      READ(KIBAN,501) FKIBCS,THMIN,THMAX,DTH
C
      NRMAX=RMAX/DR+1.01
      NBMAX=(BMAX-BMIN)/DB+1.01
      BMAX=BMIN+(NBMAX-1)*DB
C
      ICM=FICM
      IREL=FIREL
      KIBCS=FKIBCS
      NTHMAX=(THMAX-THMIN)/DTH+1.01
C
      IF(RV.LT.-EPSLIM) THEN
        RV=-RV
       ELSE
        RV=RV*AT**(1.0D0/3.0D0)
      END IF
      IF(RW.LT.-EPSLIM) THEN
        RW=-RW
       ELSE
        RW=RW*AT**(1.0D0/3.0D0)
      END IF
      IF(RD.LT.-EPSLIM) THEN
        RD=-RD
       ELSE
        RD=RD*AT**(1.0D0/3.0D0)
      END IF
C

C----------------------------------------
C---  NAMES OF PROJECTILE AND TARGET  ---
C----------------------------------------
C
      CALL ELEMNT(AP,ZP,MASSP,NAMEP)
      CALL ELEMNT(AT,ZT,MASST,NAMET)
C
C
C-------------------------------
C---  OUTPUT FOR THE INPUTS  ---
C-------------------------------
C
      WRITE(6,600)
  600 FORMAT(5X,'***************************************************'
     1         ,'**************'
     2      /5X,'*****                                              '
     3         ,'         *****'
     4      /5X,'*****                    PROGRAM WFEK              '
     5         ,'         *****'
     6      /5X,'*****                                              '
     7         ,'         *****'
     8      /5X,'***************************************************'
     9         ,'**************')
C
      WRITE(6,601) MASSP,NAMEP,MASST,NAMET,NINT(ELAB)
  601 FORMAT(//3X,'==================================================='
     1       /8X,'REACTION:',3X,I3,A2,' + ',I3,A2,5X,'at ',I5,'MeV'
     2       /3X,'===================================================')
C
      WRITE(6,602) COMMENT
  602 FORMAT(/5X,'USER''S COMMENT:'/10X,A50)
C
      WRITE(6,603) AP,ZP,AT,ZT, ELAB,RMAX,DR, BMIN,BMAX,DB ,ICM,IREL
  603 FORMAT(/3X,'-- GENERAL INPUTS --'
     1       /8X,'AP ZP AT ZT          :',2F10.5,2F10.2
     2       /8X,'ELAB RMAX DR         :',3F10.5
     3       /8X,'BMIN BMAX DB         :',3F10.5
     4       /8X,'ICM  IREL            :',2I5)
C
      IF(ABS(ZP*ZT).GT.EPSLIM) WRITE(6,604)
  604 FORMAT(/8X,'NOTE: BOTH PARTICLES ARE CHARGED BUT NO COULOMB'
     1          ,' INTERACTION IS INCLUDED')
C
      WRITE(6,605)
  605 FORMAT(/3X,'-- OPTICAL POTENTIAL PARAMETERS --')
C
      IF(ABS(V0).GT.EPSLIM) THEN
        WRITE(6,606) V0,RV,AV
       ELSE
        WRITE(6,607)
      END IF
  606 FORMAT(8X,'V0   RV   AV  :',3F10.5)
  607 FORMAT(8X,'V0   RV   AV  : NO REAL POT.')
      IF(ABS(W0).GT.EPSLIM) THEN
        WRITE(6,608) W0,RW,AW
       ELSE
        WRITE(6,609)
      END IF
  608 FORMAT(8X,'W0   RW   AW  :',3F10.5)
  609 FORMAT(8X,'W0   RW   AW  : NO VOLUME IMAGINARY POT.')
      IF(ABS(W0).GT.EPSLIM) THEN
        WRITE(6,610) WD,RD,AD
       ELSE
        WRITE(6,611)
      END IF
  610 FORMAT(8X,'WD   RD   AD  :',3F10.5)
  611 FORMAT(8X,'WD   RD   AD  : NO SURFACE IMAGINARY POT.')
C
      WRITE(6,612) KIBCS,THMIN,THMAX,DTH
  612 FORMAT(/3X,'-- CONTROL OF OUTPUT --'
     1       /8X,'KIBCS  THMIN THMAX DTH        :',I5,3F10.5)
C
      RETURN
C
      END
C
C
C
C=======================================================================
      FUNCTION BESSJ0(Z)
      IMPLICIT REAL*8 (A-H,O-Z)
C=======================================================================
C
CCC---------------------------------------------------------------------
C
      IF(ABS(Z).LE.8.0D0) THEN
        P00= 0.5756849057376692D11
        P01=-0.1336259035395191D11
        P02= 0.6516196407500873D09
        P03=-0.1121442418227296D08
        P04= 0.7739233017240460D05
        P05=-0.1849052456430540D03
C
        Q00= 0.5756849041067652D11
        Q01= 0.1029532985766340D10
        Q02= 0.9494680718778989D07
        Q03= 0.5927264853423237D05
        Q04= 0.2678532712126684D03
        Q05= 1.0D0
C
        Y=Z**2
        BESSJ0=(P00+Y*(P01+Y*(P02+Y*(P03+Y*(P04+Y*P05)))))
     1        /(Q00+Y*(Q01+Y*(Q02+Y*(Q03+Y*(Q04+Y*Q05)))))
C
       ELSE
        Y=(8.0D0/Z)**2
C
        P00= 0.9999999999947955D0
        P01=-0.109863241367192D-2
        P02= 0.2737589862975D-4
        P03=-0.216060245228D-5
        P04= 0.30994056326D-6
        P05=-0.4048673907D-7
C
        P0=P00+Y*(P01+Y*(P02+Y*(P03+Y*(P04+Y*P05))))
C
        P00=-0.1562499999676469D-1
        P01= 0.14305089740097D-3
        P02=-0.692763042016D-5
        P03= 0.80888732887D-6
        P04=-0.14748826511D-6
        P05= 0.2176282453D-7
C
        Q0=8.0D0/ABS(Z)*(P00+Y*(P01+Y*(P02+Y*(P03+Y*(P04+Y*P05)))))
C
        PI=ACOS(-1.0D0)
        ARG=ABS(Z)-PI/4.0D0
        BESSJ0=SQRT(2.0D0/PI/ABS(Z))*(P0*COS(ARG)-Q0*SIN(ARG))
      END IF
C
      RETURN
      END
C
C
C
C-----------------------------------------------------------------------
C---  Subroutines below are made by Y. Iseri (Chiba-Keizai College)  ---
C-----------------------------------------------------------------------
CCC
CCC        ***** FOPEN *****
CCC
      SUBROUTINE FOPEN
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER FNAME*50,STA*8,COMM*60,OFF*1
C ==========================================
C --- UNIT-5 ---
C     WRITE(*,100)
C 100 FORMAT(//1x,'>>> Input File-Name of Unit 5 ==> ',$)
C     READ(*,'(A)') FNAME

C     OPEN(5,FILE=FNAME,STATUS='OLD')

C --- OTHER UNITS ---
      READ(5,'(A)') COMM
      WRITE(*,200) COMM
  200 FORMAT(/1x,'---- (Comment in File) ----'
     1       /1x,' << ',A,' >>'
     2      //1x,'---- (  Open  Files  ) ----'
     3      //1x,'  Unit  Status      File-Name')

    1 READ(5,300) OFF,IUNIT,STA,FNAME
  300 FORMAT(A1,I3,1X,A8,2X,A50)
       IF(OFF.NE.' ') GO TO 1
       IF(IUNIT.GE.100 .OR. IUNIT.LT.0) GO TO 900
       WRITE(*,310) IUNIT,STA,FNAME
  310  FORMAT(1x,3X,I2,3X,A,3X,A)
       OPEN(IUNIT,FILE=FNAME,STATUS=STA)
      GO TO 1

C --- ONE MORE LINE FOR DISCRIMINATION ---
  900 READ(5,'(A)') COMM
      RETURN
      END
CCCCC
CCCCC      ***** ELEMNT *****
CCCCC
      SUBROUTINE ELEMNT(A,Z,MASS,NAME)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (MMAX=109)
      CHARACTER*2 NAME
      CHARACTER*2 NAT(MMAX),NAP(3),NAA(2),NEU,NON
      DATA  NEU,NAT / 'n ',
     1      'H ','He','Li','Be','B ','C ','N ','O ','F ','Ne','Na',
     2      'Mg','Al','Si','P ','S ','Cl','Ar','K ','Ca','Sc','Ti',
     3      'V ','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As',
     4      'Se','Br','Kr','Rb','Sr','Y ','Zr','Nb','Mo','Tc','Ru',
     5      'Rh','Pd','Ag','Cd','In','Sn','Sb','Te','I ','Xe','Cs',
     6      'Ba','La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy',
     7      'Ho','Er','Tm','Yb','Lu','Hf','Ta','W ','Re','Os','Ir',
     8      'Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn','Fr','Ra',
     9      'Ac','Th','Pa','U ','Np','Pu','Am','Cm','Bk','Cf','Es',
     A      'Fm','Md','No','Lr','Rf','Db','Sg','Bh','Hs','Mt'/
      DATA NAP /'p ','d ','t '/,  NAA/'-p','-n'/, NON/'??'/
CCC ============================================================
      MASS=NINT(A)
      IZ=NINT(Z)
C
      IF(IZ.EQ.-1) THEN
         NAME= NAA(1)
      ELSE IF(IZ.EQ. 0) THEN
         NAME= NEU
      ELSE IF(IZ.EQ. 1 .AND. MASS.LE.3) THEN
         NAME= NAP(MASS)
      ELSE IF(IZ.GE.1 .AND. IZ.LE.MMAX) THEN
         NAME= NAT(IZ)
      ELSE
         NAME= NON
      END IF
C
      RETURN
      END
