      PROGRAM NFDEC
      IMPLICIT REAL*8 (A-H,O-Z)   
      COMPLEX*16 FCL,FTT,FNR,FFR,SMTR,COUL,QMINU,QPLUS,FACTL
      COMMON/HELP/PLG(500),QLG(500)
      COMMON/HELP1/ETA,PI,SIGMA0,SK,THETAR,FCLR,FCLI 
      DIMENSION SMTR(500),COUL(500),DELTA(500),SABS(500)
      PI = 3.141592653589793238462643D0
      RADF=PI/180.0D0
      open(5,file='SMATR.PLT',form='formatted',status='old') 
      open(6,file='oc115nf.out',form='formatted')
	open(12,file='test.dat',form='formatted')
      READ(5,10)LMIN,LMAX,LSTP,AGMAX,AGSTP,Z12,SM12,SK
c     	write(*,*) LMIN,LMAX,LSTP,AGMAX,AGSTP,Z12,SM12,SK
  10  FORMAT(3I5,5F10.4)
C----Z12=Z1*Z2, SM12 is the reduced mass in MeV/c**2, SK is the wave 
C----number in 1/fm. They can be found in PTOLEMY output.-----------
      ETA=1.43986D0*Z12*SM12/(197.32858D0)**2/SK
      LMI=LMIN+1
      LMA=LMAX+1
      SIGMAR=0.0
      DO 2 L=LMI,LMA,LSTP
      READ(5,20) L1,L2,LX,SR,SI,TEMP,PHASE,PCB,TL     
  20  FORMAT (3I5,G16.8,2X,G16.8,2X,4G15.6)
ccmoj        READ(5,2000) L1,L2,LX,SR,SI,TEMP,PHASE,TL
	write(12,*) L1,L2,LX,SR,SI,TEMP,PHASE,TL
c  20  FORMAT (3I5,G16.8,2X,G16.8,2X,3G15.6)
  22  format(3I5,G16.8,1X,G13.8,3X,3G15.6)
2000  FORMAT (I3,2I5,G17.8,2X,G12.8,3X,G16.10,2X,F10.4,3X,G13.10)
      IF(L1.NE.L-1) THEN
      WRITE(6,30)
      	write(12,2000)L,L1,LMIN,LMAX,LSTP,AGMAX,AGSTP,Z12,SM12,SK 
  30  FORMAT(1X,' CHECK CONSISTENCY IN L-VALUES')
      STOP
      ELSE
      SMTR(L)=DCMPLX(SR,SI)
      SABSL=DSQRT(SR**2+SI**2)
      PHASEN=DATAN2D(SI,SR)
      ETAL=1.0D0-(SR**2+SI**2)
      SIGMAR=SIGMAR+DFLOAT(2*L-1)*ETAL
C----Calculate the point Coulomb phase---------------
      ACOU=2.0D0*SIGMAL(ETA,L-1)
      COUL(L)=DCMPLX(DCOS(ACOU),DSIN(ACOU))
      END IF
C----Check if we get the right input for S-matrix-----
c      write(6,20) L1,L2,LX,SR,SI,SABSL,0.5*phasen,0.5*acou,TL
      SABS(L)=SABSL
   2  DELTA(L)=PHASEN
      NTMAX=NINT(AGMAX/AGSTP)
      SIGMA0=SIGMAL(ETA,0)
      DO 3 NTHETA=1,NTMAX
      THETA=DFLOAT(NTHETA)*AGSTP
      THETAR=RADF*THETA
      X=DCOS(THETAR)
      Z=DSIN(0.5D0*THETAR) 
      IF(LSTP.EQ.1) THEN
C----Calculate Rutherford amplitude-----------------
      FACTR=-ETA/(2.0D0*SK*Z**2)
      ARGC=2.0D0*(SIGMA0-ETA*DLOG(Z))
      FCLR=FACTR*DCOS(ARGC)
      FCLI=FACTR*DSIN(ARGC)
      ELSE
C----Calculate Mott amplitude-----------------------
      Z1=DSIN(0.5D0*(PI-THETAR)) 
      FACTR=-ETA/(2.0D0*SK*Z**2)
      FACTR1=-ETA/(2.0D0*SK*Z1**2)
      ARGC=2.0D0*(SIGMA0-ETA*DLOG(Z))
      ARGC1=2.0D0*(SIGMA0-ETA*DLOG(Z1))
      FCLR=FACTR*DCOS(ARGC)+FACTR1*DCOS(ARGC1)
      FCLI=FACTR*DSIN(ARGC)+FACTR1*DSIN(ARGC1)
      END IF
      FCL=DCMPLX(FCLR,FCLI)
      CALL FRUTH(FCNR,FCNI,FCFR,FCFI)
      CALL PQL(X,LMI,LMA) 
      LL=LMI
      FTT=0.0
      FNR=0.0
      FFR=0.0
   4  IF(LL.GT.LMA) GO TO 5
      PL=PLG(LL)
      QL=QLG(LL)
      QMINU=DCMPLX(0.5*PL,1.0D0/PI*QL)  
      QPLUS=DCMPLX(0.5*PL,-1.0D0/PI*QL)  
      FACTL=DCMPLX(DFLOAT(2*LL-1),0.0D0)
      FACTL=FACTL*COUL(LL)*(SMTR(LL)-DCMPLX(1.0D0,0.0D0))
      FTT=FTT+FACTL*DCMPLX(PL,0.0D0)/DCMPLX(0.0D0,2.0D0*SK)
      FNR=FNR+FACTL*QMINU/DCMPLX(0.0D0,2.0D0*SK)
      FFR=FFR+FACTL*QPLUS/DCMPLX(0.0D0,2.0D0*SK)
      LL=LL+LSTP
      GO TO 4
   5  FACTL=DCMPLX(DFLOAT(LSTP),0.0D0)
      FTT=FACTL*FTT+FCL
      SIGC=10.0D0*CDABS(FCL)**2
      SIGT=10.0D0*CDABS(FTT)**2/SIGC
C----For a symmetric system (12C+12C etc.) the near-far decomposition
C----should be performed for 'unsymmetric' cross section only, i.e. 
C----PTOLEMY output with LSTP=1. See also Fuller, P.R. C12 (1975) 1561.
      IF(LSTP.EQ.1) THEN
      FNR=FACTL*FNR+DCMPLX(FCNR,FCNI)
      FFR=FACTL*FFR+DCMPLX(FCFR,FCFI)
      SIGN=10.0D0*CDABS(FNR)**2/SIGC
      SIGF=10.0D0*CDABS(FFR)**2/SIGC
      SIGNC=10.0D0*((FCNR)**2+(FCNI)**2)/SIGC
      SIGFC=10.0D0*((FCFR)**2+(FCFI)**2)/SIGC
      ELSE
      SIGN=0.0
      SIGF=0.0
      SIGNC=0.0     
      SIGFC=0.0
      END IF
      WRITE(6,40)THETA,SIGT,SIGN,SIGF,SIGC,SIGNC,SIGFC
  40  FORMAT(F10.2,6G15.6)
   3  CONTINUE
      SIGMAR=SIGMAR*10.0D0*DFLOAT(LSTP)*PI/SK**2
      WRITE(6,50)
  50  FORMAT(/,3X,' THETA',3X,' SIGT./RUTH.',4X,'SIGN./RUTH.',
     *4X,'SIGF./RUTH.',4X,'RUTH. [mb]',4X,'SIGNC./RUTH.',4X,
     *'SIGFC./RUTH.'/)
      DO 6 L=LMI+1,LMA,LSTP
      DEFLN=DELTA(L)-DELTA(L-1)
      DLH=DFLOAT(L)-0.5D0
      DEFLC=2.0D0*DATAND(ETA/DLH)
      if(defln.gt.0.0)defln=defln-360.0
      DEFLS=DEFLN+DEFLC
      DRK=DLH/SK
   6  WRITE(6,60) L-1,DRK,SABS(L),DEFLS,DEFLN,DEFLC
  60  FORMAT(3X,I3,5G16.8)
      WRITE(6,70)
  70  FORMAT(/,3X,' L',5X,' D [fm]',8X,' |S(L)|',
     *8X,' DEFS.[deg]',6X,' DEFN.[deg]',4X,' DEFC.[deg]'/)
      WRITE(6,80) ETA,SIGMAR
  80  FORMAT(/,3X,' ETA =',F10.6,3X,' REACTION SIGMA [mb]=',
     *F10.3,/,/,3X,' NFDEC results obtained by using PTOLEMY 
     *output.'/,3X,' With best regards from Dao T. Khoa !!!' )
      END
C---------------------------------------------------------------------------
      SUBROUTINE FRUTH(FNR,FNI,FFR,FFI)
      IMPLICIT REAL*8 (A-H,O-Z)   
      COMPLEX*16 BB,CZ,A,B,ARG,PSI,S,TERM,AMPN,AMPF,AMPT
      COMMON/HELP1/ETA,PI,SIGMA0,SK,AN,FCLR,FCLI 
C----Near-Far decomposition of the Rutherford amplitude using Eqs.(14) and
C----(15) from R.C. Fuller, Phys. Rev. C12 (1975) 1563---------------------
C      CALL LIB$FLT_UNDER(0)
      G=0.577215664901532860606512D0
      CONVG=1.0D-9
      Z=DSIN(0.5D0*AN)**2 
      CZ=DCMPLX(Z,0.0D0)
      IF(ETA.EQ.0.0D0) THEN
      AMPN=-EXP(2.*DCMPLX(0.0,1.0)*SIGMA0)/(4.*SK*PI*Z)
      AMPF=-AMPN
      ELSE
C
      BB = DCMPLX(1.D0,ETA)	
      IF(Z.GT.0.5D0) GOTO 100
      S=1.0/BB
      N=1
  20  TERM=CZ**N/(BB+FLOAT(N))
      S=S+TERM
      N=N+1
      IF(CDABS(TERM)/CDABS(S).GT.CONVG) GO TO 20
      GOTO 300
 100  Z =1.0D0-Z
      CZ=DCMPLX(Z,0.0D0)
      A=1.0
C----START BY EVALUATING PSI( 30 + i ETA)---------------------------
      ARG=DCMPLX(30.D0,ETA)
C----USING ASYMPTOTIC FORM A&S 6.3.18-------------------------------
      PSI=CDLOG(ARG)-0.5/ARG-1.0/12.0/(ARG**2)+
     & 1.0/120.0/(ARG**4)-1.0/252.0/(ARG**6)
C----NOW BRING IT BACK TO PSI( 1 + i ETA)---------------------------
      DO 110 I=1,29		
      DFI=DFLOAT(I)
 110  PSI=PSI-1.0/DCMPLX(DFI,ETA)
      B=-CDLOG(CZ)-PSI-DCMPLX(G,0.0D0)
      S=A*B			
      N=1
 120  DFN=DFLOAT(N)
      A=A*DCMPLX(DFN,ETA)/FLOAT(N)
      B=B+DCMPLX(0.D0,ETA)/FLOAT(N)/DCMPLX(DFN,ETA)
      TERM=A*B*(CZ**N)
      S=S+TERM
      N=N+1
      IF(CDABS(TERM)/CDABS(S).GT.CONVG) GO TO 120
      Z=1.0-Z
      CZ=DCMPLX(Z,0.0D0)
 300  S=DCMPLX(0.D0,0.5/PI)*(CZ**DCMPLX(1.D0,ETA))*S
      XKT=DEXP(-2.0D0*PI*ETA)
      XKTN=1.0D0/(1.0D0-XKT)
      XKTF=-XKT/(1.0D0-XKT)
C----THE NEAR AMPLITUDE--------------------------------------------
      AMPN=DCMPLX(XKTN,0.0D0)-S
C----THE FAR AMPLITUDE--------------------------------------------
      AMPF=DCMPLX(XKTF,0.0D0)+S
      AMPT=DCMPLX(FCLR,FCLI)
      AMPN=AMPN*AMPT
      AMPF=AMPF*AMPT
      END IF
C
      FNR=DREAL(AMPN)
      FFR=DREAL(AMPF)
      FNI=DIMAG(AMPN)
      FFI=DIMAG(AMPF)
      RETURN
      END
C--------------------------------------------------------------
      FUNCTION SIGMAL(ETA,L)
      IMPLICIT REAL*8 (A-H,O-Z)   
C---USING FORM. 6.1.27 AND 6.3.2 FROM ABRAMOWITZ & STEGUN
      GAMMA=0.577215664901532860606512D0
      CONVG=1.0D-9
      IF(L.EQ.0) THEN 
      SIG0=-ETA*GAMMA
      ELSE
      SK=0.0D0      
      DO 1 K=1,L
   1  SK=SK+1.0D0/DFLOAT(K)
      SIG0=ETA*(-GAMMA+SK)
      END IF
      N=0
      SN=0.0D0
   2  FACTN=1.0D0+DFLOAT(L+N)
      TERM=ETA/FACTN-DATAN(ETA/FACTN)
      SN=SN+TERM
      N=N+1
      IF(DABS(TERM)/DABS(SN).GT.CONVG) GO TO 2
      SIGMAL=SIG0+SN
      RETURN
      END
C-------------------------------------------------------------
      SUBROUTINE PQL(X,LMI,LMA) 
      IMPLICIT REAL*8 (A-H,O-Z)   
      COMMON/HELP/PLG(500),QLG(500)
C---X=COS(THETA), 0 < THETA < PI-----------------------------
C---LMI MUST BE EITHER 1 OR AN INTEGER > 3--------------------
      PL0=1.0D0
      PL1=X 
      QL0=0.5D0*DLOG((1.0D0+X)/(1.0D0-X))
      QL1=X*QL0-1.0D0 
      IF(LMI.EQ.1) THEN 
      PLG(1)=PL0 
      QLG(1)=QL0 
      PLG(2)=PL1 
      QLG(2)=QL1 
      L1=3
      GO TO 20
      ELSE 
      DO 10 L=3,LMI-1
      PPL=DFLOAT(2*L-3)*X*PL1-DFLOAT(L-2)*PL0
      QQL=DFLOAT(2*L-3)*X*QL1-DFLOAT(L-2)*QL0
      PP=PPL/DFLOAT(L-1)
      QQ=QQL/DFLOAT(L-1)
      PL0=PL1
      QL0=QL1
      PL1=PP
   10 QL1=QQ
      L1=LMI
      END IF
   20 DO 30 L=L1,LMA
      PPL=DFLOAT(2*L-3)*X*PL1-DFLOAT(L-2)*PL0
      QQL=DFLOAT(2*L-3)*X*QL1-DFLOAT(L-2)*QL0
      PP=PPL/DFLOAT(L-1)
      QQ=QQL/DFLOAT(L-1)
      PL0=PL1
      QL0=QL1
      PL1=PP
      QL1=QQ
      PLG(L)=PP
   30 QLG(L)=QQ
      RETURN
      END
