!'''
!@File    :   CC_python.py Ec_version solved by modified Numerove or Discrete Basis method 
!@Time    :   2025/06/22 
!@Author  :   Liao ZeHong, Kouichi Hagino
!@Version :   1.0
!@Contact :   liaozh26@mail2.sysu.edu.cn
! gfortran  ccfull.f90  -llapack -lblas -o solve
!'''


Module ccfull_initialization_mod
    Implicit None   
  !  """Principal global variables
  !  P                   - penetrability 
  !  Sigma               - fusion cross section, unit mb
  !  Spin                - mean angular momentum
  !    
  !  Apro, Zpro, Rpro    - atomic #, proton # and radius of the projectile
  !  Atar, Ztar, Rtar    - those of the target
  !  ReduceMass          - reduced mass, unit MeV/c**2
  !  Amu                 - nucleon mass, unit MeV/c**2
  !  E                   - bombarding energy in the center of mass frame, unit MeV
  !    
  !  V0,R0,A0            - depth, range, and dissuseness parameters of uncoupled 
  !                        nuclear potential, which is assumed to be a Woods-Saxon form
  !    
  !  IVIBROTT (IVIBROTP) - option for intrinsic degree of freedom 
  !                        = -1; no excitation (inert)
  !                        =  0; vibrational coupling 
  !                        =  1; rotational coupling
  !
  !  Ntar (Npro)         - the number of levels to be included    
  !
  !  BetaTn (BetaPn)     - i am not very clear about it (Liao)
  !  BetaT (BetaP)       - defeormation parameter
  !  OmegaT (OmegaP)     - excitation energy of the oscillator
  !  LambdaT (LambdaP)   - multipolarity
  !  NphononT (NphononP) - the number of phonon to be included
  !  Beta2T (Beta2P)     - static quadrupole deformation parameter
  !  Beta4T (Beta4P)     - static hexadecapole deformation parameter
  !  BetaT2n             - is the same as BetaTn but for the second mode
  !  BetaT2              - is the same as BETAT but for the second mode
  !  OmegaT2             - is the same as OmegaT but for the second mode
  !  LambdaT2            - is the same as LambdaT but for the second mode
  !  NphononT2           - is the same as NphononT but for the second mode
  !
  !  E2T (E2P)           - excitation energy of 2+ state in a rotational band
  !  NrotT (NrotP)       - the number of levels in the rotational band to be included 
  !                        (up to I^pi=2*NROT+ states are included)
  !
  !  L                   - angular momentum of the relative motion
  !  CPOT                - coupling matrix
  !  """
  Real(8), parameter :: Hbar = 197.329d0
  Real(8), parameter :: Pi   = 3.141592653d0
  Real(8), parameter :: Amu  = 938.0
  Real(8), parameter :: e2   = 1.44

  integer, parameter :: Nlevelmax = 30
  integer, Public ::  Apro, Zpro, Atar, Ztar, L_i
  Real(8), Public ::  R0P, R0T, Rpro, Rtar, ReduceMass, h2m
  Real(8), Public ::  OmegaT, BetaT, OmegaT2, BetaT2, E2T, Beta2T, Beta4T
  Real(8), Public ::  OmegaP, BetaP, E2P, Beta2P, Beta4P
  Real(8), Public ::  V0, R0, A0, Emin, Emax, dE, R_max, R_min, dr
  Real(8), Public ::  BetaTn, BetaPn, BetaT2n, BetaP2n
  Real(8), Public ::  R_barrier, V_barrier, curv, R_bottom, V_bottom
  Real(8), Public ::  R_barrier_l, V_barrier_l, curv_l, R_bottom_l, V_bottom_l
  Real(8), Public ::  t
  Real,    Public ::  Ftr, Qtrans
  Real,    Public ::  E
  Integer, Public ::  IVIBROTP, IVIBROTT, LambdaT, NphononT, LambdaT2, NphononT2
  Integer, Public ::  NrotT, Ntar, LambdaP, NphononP, NrotP, Npro
  Integer, Public ::  Ntrans, iq
  Integer, Public ::  Nlevel, R_iterat
  Real(8), Allocatable, Public :: Pot_para1(:), Pot_para2(:), Pot_para3(:)
  Real(8), Allocatable, Public :: Ev(:), eps(:)
  Real(8), Allocatable, Public :: omeahv(:), omeahv2(:), omeahvp(:), erott(:), erotp(:)
  Real(8), Allocatable, Public :: betnahv(:,:), betcahv(:,:), betnahv2(:,:), betcahv2(:,:)
  Real(8), Allocatable, Public :: betnahvp(:,:), betcahvp(:,:), bett(:,:), betp(:,:)
  Real(8), Allocatable, Public :: Evec(:,:), CPOT0(:,:), CPOTH(:,:), CPOT(:,:,:)
  
  Integer, Allocatable, Public :: imutual(:,:)


  Integer, Public ::  NumEmulator, Nx, Nx_cc, Nx1, Nx2
  Real(8),    Allocatable, Public ::  Echannel(:)
  complex*16, Allocatable, Public :: wf(:)
  complex*16, Allocatable, Public :: wf_base(:,:), wf_base2(:,:), wf_base3_dagger(:,:)
  complex*16, Allocatable, Public :: Hamilton(:,:), unit(:,:)
  complex*16, Allocatable, Public :: Ec_psi(:,:), Ec_psi0(:), Ec_psi0_2(:)
  complex*16, Allocatable, Public :: cc_ec(:,:), cc_ec_inv(:,:),dd_ec(:), dd2_ec(:), a_tra(:)
  
  character(len=2), dimension(0:109) :: Element

  data Element / &
       'O ','H ','He','Li','Be','B ','C ','N ','O ','F ','Ne', &
       'Na','Mg','Al','Si','P ','S ','Cl','Ar','K ','Ca','Sc','Ti','V ', &
       'Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As','Se','Br','Kr', &
       'Rb','Sr','Y ','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In', &
       'Sn','Sb','Te','I ','Xe','Cs','Ba','La','Ce','Pr','Nd','Pm','Sm', &
       'Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb','Lu','Hf','Ta','W ','Re', &
       'Os','Ir','Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn','Fr','Ra', &
       'Ac','Th','Pa','U ','Np','Pu','Am','Cm','Bk','Cf','Es','Fm','Md', &
       'No','Lr','XX','X1','X2','X3','X4','04' /

Contains

  Subroutine Read_Input()

    Open(10, File='ccfull.inp', Status = 'Unknown')
    Open(20, File='sigma.dat',  Status = 'Unknown')
    write(20,'(4A15)') 'E',  'sigma',  '<L>',  'P (L = 0)'
    Open(30, File='Output.dat', Status = 'Unknown')
    Open(40, File='Wavefunction_Db.dat', Status = 'Unknown')
    Read(10,*)Apro, Zpro, Atar, Ztar
    Read(10,*)R0P, IVIBROTP, R0T, IVIBROTT

    Rpro = R0P * Apro ** (1.d0 / 3.d0)
    Rtar = R0T * Atar ** (1.d0 / 3.d0)
    ReduceMass = Apro * Atar * 1.d0 / (Apro + Atar) * Amu

    h2m = Hbar**2 / (2 * ReduceMass)
    
    If (IVIBROTT == 0) Then

      Read(10,*)OmegaT, BetaT, LambdaT, NphononT
      Read(10,*)OmegaT2, BetaT2, LambdaT2, NphononT2
      Ntar = NphononT

    Else If (IVIBROTT == 1) Then

      Read(10,*)E2T, Beta2T, Beta4T, NrotT 
      Read(10,*)
      Ntar = NrotT
    Else
      Read(10,*)
      Read(10,*)
      Ntar = 0
    End If

    If (IVIBROTP == 0) Then

      Read(10,*)OmegaP, BetaP, LambdaP, NphononP
      Npro = NphononP

    Else If (IVIBROTP == 1) Then

      Read(10,*)E2P, Beta2P, Beta4P, NrotP
      Npro = NrotP
    Else
      Read(10,*)
      Npro = 0
    End If
    Nlevel=(Npro + 1)*(Ntar + 1)

    Read(10,*)Qtrans, Ftr
    Read(10,*)V0, R0, A0
    Read(10,*)Emin, Emax, dE
    Read(10,*)R_max, dr
  End Subroutine Read_Input

  Subroutine Output_information()
      Character(len=1) :: ans


      print '(A, I0 ,A, A,I0 ,A)', 'System: ', Int(Apro), Element(Zpro),'+', Int(Atar), Element(Ztar)

      print '(A, F8.3, A, F8.3, A, F8.3, A, F8.3, A)', &
            "Simulation range E from ", Emin, " MeV To ", Emax, " MeV, dE = ", dE, " MeV"

      print '(A, F8.4, A, F8.4, A, F8.3, A, F8.3, A)', &
            "Simulation range R from ", R_min, " fm To ", R_max, " fm, dR = ", dr, " fm"

      print '(A, F8.3, A, F8.3, A, F8.3, A, F8.3, A)', &
            "Potential parameters: V0= ", V0, "(MeV), a= ", A0, "(fm), r0= ", R0, "(fm), R= ",&
             R0*(Apro**(1.0/3.0) + Atar**(1.0/3.0)), "(fm)"

      print '(A, F8.4, A)', "Coulomb barrier position :", R_barrier, "fm"

      print '(A, F8.4, A)', "Coulomb barrier energy   :", V_barrier, "MeV"

      print '(A, F8.4, A)', "Coulomb barrier curv     :", curv, "MeV"

      print '(A, F8.4, A)', "Coulomb bottom position  :", R_bottom, "fm"

      print '(A, F10.4, A)', "Coulomb bottom energy    :", V_bottom, "MeV"

      write(30,'(A, I0 ,A, A,I0 ,A)') 'System: ', Int(Apro), Element(Zpro),'+', Int(Atar), Element(Ztar)

      write(30,'(A, F8.3, A, F8.3, A, F8.3, A, F8.3, A)') &
            "Simulation range E from ", Emin, " MeV To ", Emax, " MeV, dE = ", dE, " MeV"

      write(30,'(A, F8.4, A, F8.4, A, F8.3, A, F8.3, A)') &
            "Simulation range R from ", R_min, " fm To ", R_max, " fm, dR = ", dr, " fm"

      write(30,'(A, F8.3, A, F8.3, A, F8.3, A, F8.3, A)') &
            "Potential parameters: V0= ", V0, "(MeV), a= ", A0, "(fm), r0= ", R0, "(fm), R= ",&
             R0*(Apro**(1.0/3.0) + Atar**(1.0/3.0)), "(fm)"

      write(30,'(A, F8.4, A)') "Coulomb barrier position :", R_barrier, "fm"

      write(30,'(A, F8.4, A)') "Coulomb barrier energy   :", V_barrier, "MeV"

      write(30,'(A, F8.4, A)') "Coulomb barrier curv     :", curv, "MeV"

      write(30,'(A, F8.4, A)') "Coulomb bottom position  :", R_bottom, "fm"

      write(30,'(A, F10.4, A)') "Coulomb bottom energy    :", V_bottom, "MeV"

      if (Ntar /= 0) then
          if (IVIBROTT == 0) then
              write(*,'(A, F6.3, A, F6.3, A, I0, A, I0)') "Phonon Excitation in the targ.: beta=",&
                         BetaT, ", omega=", OmegaT, " (MeV), Lambda=", LambdaT, ", Nph=", NphononT
              BetaTn = BetaT
              Write(*,*)' Different beta_N from beta_C for this mode(n/y)?'
              Read(*,*)ans
              IF(ans .Eq. 'Y' .Or. ans .Eq. 'y') Then
                Write(*,*)'beta_N=?'
                Read(*,*)BetaTn
              End If
              write(*,'(A, F6.3, A, F6.3, A, F6.3, A)') "Phonon Excitation in the targ.: beta_N=",&
                                                 BetaTn, ", beta_C=", BetaT, ", r0=", R0T, " (fm),"
              write(*,'(A, F6.3, A, I0, A, I0)') "                              omega=", OmegaT, "&
                                                 (MeV), Lambda=", LambdaT, ", Nph=", NphononT
              write(30,'(A, F6.3, A, F6.3, A, F6.3, A)') "Phonon Excitation in the targ.: beta_N=",&
                                                   BetaTn, ", beta_C=", BetaT, ", r0=", R0T, " (fm),"
              write(30,'(A, F6.3, A, I0, A, I0)') "                              omega=", OmegaT, &
                                                      " (MeV), Lambda=", LambdaT, ", Nph=", NphononT
          else if (IVIBROTT == 1) then
              write(*,'(A, F6.3, A, F6.3, A, F6.3, A)') "Rotational Excitation in the targ.: beta2=", Beta2T, &
                                                        ", beta4=", Beta4T, ", r0=", R0T, " (fm),"
              write(*,'(A, F6.3, A, I0)') "                                   E2=", E2T, " (MeV), Nrot=", NrotT
              write(30,'(A, F6.3, A, F6.3, A, F6.3, A)') "Rotational Excitation in the targ.: beta2=", Beta2T,&
                                                         ", beta4=", Beta4T, ", r0=", R0T, " (fm),"
              write(30,'(A, F6.3, A, I0)') "                                   E2=", E2T, " (MeV), Nrot=", NrotT
          end if
      end if

      if (NphononT2 /= 0) then
          write(*,'(A, F6.3, A, F6.3, A, I0, A, I0)') "Phonon Excitation in the targ.: beta=", BetaT2, ", omega=",&
                                                         OmegaT2, " (MeV), Lambda=", LambdaT2, ", Nph=", NphononT2
          BetaT2n = BetaT2
          Write(*,*)' Different beta_N from beta_C for this mode(n/y)?'
          Read(*,*)ans
          IF(ans .Eq. 'Y' .Or. ans .Eq. 'y') Then
            Write(*,*)'beta_N=?'
            Read(*,*)BetaT2n
          End If
          write(*,'(A, F6.3, A, F6.3, A, F6.3, A)') "Phonon Excitation in the targ.: beta_N=", BetaT2n,&
                                                     ", beta_C=", BetaT2, ", r0=", R0T, " (fm),"
          write(*,'(A, F6.3, A, I0, A, I0)') "                              omega=", OmegaT2, " (MeV), Lambda=",&
                                             LambdaT2, ", Nph=", NphononT2
          write(30,'(A, F6.3, A, F6.3, A, F6.3, A)') "Phonon Excitation in the targ.: beta_N=", BetaT2n,&
                                                     ", beta_C=", BetaT2, ", r0=", R0T, " (fm),"
          write(30,'(A, F6.3, A, I0, A, I0)') "                              omega=", OmegaT2, " (MeV), Lambda=",&
                                               LambdaT2, ", Nph=", NphononT2
          Call Mutual
      End if

      print *, "------------------------------"
      print *, "Mode of excitation for projectile"

      if (Npro /= 0) then
          if (IVIBROTP == 0) then
              write(*,'(A, F6.3, A, F6.3, A, I0, A, I0)') "Phonon Excitation in the proj.: beta=", BetaP, ", omega=",&
                                                           OmegaP, " (MeV), Lambda=", LambdaP, ", Nph=", NphononP
              BetaPn = BetaP
              Write(*,*)' Different beta_N from beta_C for this mode(n/y)?'
              Read(*,*)ans
              IF(ans .Eq. 'Y' .Or. ans .Eq. 'y') Then
                Write(*,*)'beta_N=?'
                Read(*,*)BetaPn
              End If
              write(*,'(A, F6.3, A, F6.3, A, F6.3, A)') "Phonon Excitation in the proj.: beta_N=", BetaPn, ", beta_C=",&
                                                         BetaP, ", r0=", R0P, " (fm),"
              write(*,'(A, F6.3, A, I0, A, I0)') "                              omega=", OmegaP, " (MeV), Lambda=",&
                                                                                       LambdaP, ", Nph=", NphononP
              write(30,'(A, F6.3, A, F6.3, A, F6.3, A)') "Phonon Excitation in the proj.: beta_N=", BetaPn,&
                                                         ", beta_C=", BetaP, ", r0=", R0P, " (fm),"
              write(30,'(A, F6.3, A, I0, A, I0)') "                              omega=", OmegaP, " (MeV), Lambda=",&
                                                   LambdaP, ", Nph=", NphononP
          else if (IVIBROTP == 1) then
              write(*,'(A, F6.3, A, F6.3, A, F6.3, A)') "Rotational Excitation in the proj.: beta2=", Beta2P,&
                                                         ", beta4=", Beta4P, ", r0=", R0P, " (fm),"
              write(*,'(A, F6.3, A, I0)') "                                   E2=", E2P, " (MeV), Nrot=", NrotP
              write(30,'(A, F6.3, A, F6.3, A, F6.3, A)') "Rotational Excitation in the proj.: beta2=", Beta2P,&
                                                         ", beta4=", Beta4P, ", r0=", R0P, " (fm),"
              write(30,'(A, F6.3, A, I0)') "                                   E2=", E2P, " (MeV), Nrot=", NrotP
          end if
      end if
      Call Anharmonicity
      Call grotation

      print *, "------------------------------"
      print *, "Transfer coupled mode"

      if (Ntrans /= 0) then
          write(*,'(A, F6.3, A, F6.3, A)') "Transfer channel: Strength= ", Ftr, ", Q = ", Qtrans, " MeV"
          write(30,'(A, F6.3, A, F6.3, A)') "Transfer channel: Strength= ", Ftr, ", Q = ", Qtrans, " MeV"
      end if
      print *, "------------------------------"
  End Subroutine Output_information


 
  Subroutine Initialize_CCFull()

    Implicit None
    Allocate(imutual(0:Nlevelmax,0:Nlevelmax))
    Allocate(Ev(Nlevelmax))
    Allocate(Evec(Nlevelmax,Nlevelmax))
    Allocate(betnahv(0:Nlevelmax,0:Nlevelmax))
    Allocate(betcahv(0:Nlevelmax,0:Nlevelmax))
    Allocate(omeahv(0:Nlevelmax))
    Allocate(betnahv2(0:Nlevelmax,0:Nlevelmax))
    Allocate(betcahv2(0:Nlevelmax,0:Nlevelmax))
    Allocate(omeahv2(0:Nlevelmax+1))
    Allocate(betnahvp(0:Nlevelmax,0:Nlevelmax))
    Allocate(betcahvp(0:Nlevelmax,0:Nlevelmax))
    Allocate(omeahvp(0:Nlevelmax))
    Allocate(bett(0:Nlevelmax,0:Nlevelmax))
    Allocate(erott(0:Nlevelmax+1))
    Allocate(betp(0:Nlevelmax,0:Nlevelmax))
    Allocate(erotp(0:Nlevelmax))
    Allocate(eps(0:Nlevelmax))






    Call Read_Input()

    Call PotShape(0, R_barrier, V_barrier, curv, R_bottom, V_bottom)
    R_iterat = int((R_max - R_bottom) / dr)
    R_min = R_bottom
    R_max = R_bottom + dr * R_iterat
    t = Hbar**2 / (2 * ReduceMass * dr**2)
    
    Nx = R_iterat
    Nx_cc = Nlevel * Nx
    Nx1 = Nx_cc + Nlevel
    Nx2 = Nx_cc + 2 * Nlevel
    NumEmulator = 2
    Call Output_information()
    Allocate(CPOT(Nlevelmax,Nlevelmax,0:R_iterat+2))
    Allocate(CPOT0(Nlevelmax,Nlevelmax))
    Allocate(CPOTH(Nlevelmax,Nlevelmax))
    Allocate(wf(Nx2))
    Allocate(wf_base(Nx2, 1:NumEmulator+Nlevel))
    Allocate(wf_base2(Nx1, 1:NumEmulator+Nlevel))
    Allocate(wf_base3_dagger(1:NumEmulator+Nlevel, Nx1))
    Allocate(Hamilton(Nx1, Nx2))
    Allocate(Ec_psi(Nx2, Nx1))
    Allocate(Ec_psi0(Nx2))
    Allocate(Ec_psi0_2(Nx1))
    Allocate(Echannel(1:Nlevelmax))
    Allocate(a_tra(Nlevel))

    



  End Subroutine Initialize_CCFull


End Module ccfull_initialization_mod


Subroutine PotShape(l, rb, vb, curv_out, rmin, vmin)
  Use ccfull_initialization_mod
  
  Implicit None
  Integer :: n, i
  Integer, Intent(In) :: l
  Real(8), Intent(Out) :: rb, vb, curv_out, rmin, vmin
  Real(8) :: r, u0, u1, ra, rb_tmp, tolk, u, ddv00
  Real(8) :: dV_dr, V
  external dV_dr, V
  r = 50.5D0
  u0 = dV_dr(r, l)

  Do
    r = r - 1.0D0
    If (r < 0.0D0) Then
      rb = -5.0D0
      vb = -5.0D0
      curv_out = -1.0D0
      rmin = -1.0D0
      vmin = -1.0D0
      Return
    End If
    u1 = dV_dr(r, l)
    If (u0 * u1 <= 0.0D0) Exit
    u0 = u1
  End Do

  ra = r + 1.0D0
  rb_tmp = r
  tolk = 1.0D-6
  n = Int(Log10(Abs(rb_tmp - ra) / tolk) / Log10(2.0D0) + 0.5D0)

  Do i = 1, n
    r = (ra + rb_tmp) / 2.0D0
    u = dV_dr(r, l)
    If (u0 * u < 0.0D0) Then
      rb_tmp = r
    Else
      ra = r
    End If
  End Do

  rb = r
  ddv00 = (dV_dr(rb + 1.0D-5, l) - dV_dr(rb - 1.0D-5, l)) / (2.0D-5)

  If (ddv00 > 0.0D0) Then
    Print *, "Error: Second derivative positive at barrier top."
    Stop
  End If

  curv_out = Hbar * Sqrt(Abs(ddv00) / ReduceMass)
  vb = V(rb, l)

  ra = rb - 0.5D0
  u0 = dV_dr(ra, l)
  rb_tmp = 0.5D0
  u1 = dV_dr(rb_tmp, l)
  n = Int(Log10(Abs(rb - ra) / tolk) / Log10(2.0D0) + 0.5D0)

  Do i = 1, n
    r = (ra + rb_tmp) / 2.0D0
    u = dV_dr(r, l)
    If (u0 * u < 0.0D0) Then
      rb_tmp = r
    Else
      ra = r
    End If
  End Do

  rmin = r
  vmin = V(rmin, l)

End Subroutine PotShape

Real(8) Function Vn(r)
  Use ccfull_initialization_mod
  Real(8), Intent(In) :: r
  Real(8) :: R12
  R12 = R0 * (Apro**(1.0D0/3.0D0) + Atar**(1.0D0/3.0D0))
  Vn = -V0 / (1.0D0 + Exp((r - R12) / A0))
End Function Vn

Real(8) Function W(r)
  Use ccfull_initialization_mod
  Real(8), Intent(In) :: r
  Real(8) :: rw
  rw = 1.0D0 * (Apro**(1.0D0/3.0D0) + Atar**(1.0D0/3.0D0))
  W = 50.0D0 / (1.0D0 + Exp((r - rw) / 0.3D0))
End Function W

Real(8) Function dVn_dr(r)
  Use ccfull_initialization_mod
  Real(8), Intent(In) :: r
  Real(8) :: R12, exp_term
  R12 = R0 * (Apro**(1.0D0/3.0D0) + Atar**(1.0D0/3.0D0))
  exp_term = Exp((r - R12) / A0)
  dVn_dr = V0 / A0 * exp_term / (1.0D0 + exp_term)**2
End Function dVn_dr

Real(8) Function dVcent_dr(r, l)
  Use ccfull_initialization_mod
  Real(8), Intent(In) :: r
  Integer, Intent(In) :: l
  dVcent_dr = -2.0D0 * l * (l + 1.0D0) * Hbar**2 / (2.0D0 * ReduceMass * r**3)
End Function dVcent_dr

Real(8) Function dVc_dr(r)
  Use ccfull_initialization_mod
  Real(8), Intent(In) :: r
  dVc_dr = -Zpro * Ztar * Hbar / 137.0D0 / r**2
End Function dVc_dr

Real(8) Function dV_dr(r, l)
  Use ccfull_initialization_mod
  Real(8), Intent(In) :: r
  Integer, Intent(In) :: l
  Real(8)             :: dVc_dr, dVn_dr, dVcent_dr
  external dVc_dr, dVn_dr, dVcent_dr
  dV_dr = dVn_dr(r) + dVc_dr(r) + dVcent_dr(r, l)
End Function dV_dr

Real(8) Function Vc(r)
  Use ccfull_initialization_mod
  Real(8), Intent(In) :: r
  Vc = Zpro * Ztar * Hbar / 137.0D0 / r
End Function Vc

Real(8) Function Vr(r, l)
  Use ccfull_initialization_mod
  Real(8), Intent(In) :: r
  Integer, Intent(In) :: l
  Vr = h2m * l * (l + 1.0D0) / r**2
End Function Vr

Real(8) Function V(r, l)
  Use ccfull_initialization_mod
  Real(8), Intent(In) :: r
  Integer, Intent(In) :: l
  Real(8)             :: Vn, Vc, Vr
  external Vn, Vc, Vr
  V = Vn(r) + Vc(r) + Vr(r, l)
End Function V

Real(8) Function VnCC(r, Xt)
  Use ccfull_initialization_mod
  Real(8), Intent(In) :: r, Xt
  Real(8) :: R12
  R12 = R0 * (Apro**(1.0D0/3.0D0) + Atar**(1.0D0/3.0D0))
  VnCC = -V0 / (1.0D0 + Exp((r - R12 - Xt) / A0))
End Function VnCC

Real(8) Function Fct(r)
  Use ccfull_initialization_mod
  IMPLICIT NONE
  Real(8), Intent(In) :: r
  Real(8) ::  result, Lambda

  Lambda = LambdaT
  If (r > Rtar) Then
     result = 3.D0 / (2.D0*Lambda + 1.D0) * Zpro * Ztar / 137.D0 * Hbar / r * (Rtar / r)**Lambda
  Else
     result = 3.D0 / (2.D0*Lambda + 1.D0) * Zpro * Ztar / 137.D0 * Hbar / Rtar * (r / Rtar)**Lambda
  End If

  result = result * BetaT / Sqrt(4.D0 * Pi)

  Fct = result
  Return
End Function Fct

Real(8) Function Fct2(r)
  Use ccfull_initialization_mod
  Implicit None
  Real(8), Intent(In) :: r
  Real(8) :: result
  Integer, Parameter :: Lambda = 2

  If (r > Rtar) Then
     result = 3.D0 / (2.D0*Lambda + 1.D0) * Zpro * Ztar / 137.D0 * Hbar / r * (Rtar / r) ** Lambda
  Else
     result = 3.D0 / (2.D0*Lambda + 1.D0) * Zpro * Ztar / 137.D0 * Hbar / Rtar * (r / Rtar) ** Lambda
  End If

  Fct2 = result
  Return
End Function Fct2

Real(8) Function Fct3(r)
  Use ccfull_initialization_mod
  Implicit None
  Real(8), Intent(In) :: r
  Real(8) :: result
  Integer, Parameter :: Lambda = 3

  If (r > Rtar) Then
     result = 3.D0 / (2.D0*Lambda + 1.D0) * Zpro * Ztar / 137.D0 * Hbar / r * (Rtar / r) ** Lambda
  Else
     result = 3.D0 / (2.D0*Lambda + 1.D0) * Zpro * Ztar / 137.D0 * Hbar / Rtar * (r / Rtar) ** Lambda
  End If

  result = result * (Beta4T + 7.D0 * Beta2T**2 / 7.D0 / Sqrt(Pi))

  Fct3 = result
  Return
End Function Fct3

Real(8) Function Fct4(r)
  Use ccfull_initialization_mod
  Implicit None
  Real(8), Intent(In) :: r
  Real(8) :: result
  Integer, Parameter :: Lambda = 4

  If (r > Rtar) Then
     result = 3.D0 / (2.D0*Lambda + 1.D0) * Zpro * Ztar / 137.D0 * Hbar / r * (Rtar / r) ** Lambda
  Else
     result = 3.D0 / (2.D0*Lambda + 1.D0) * Zpro * Ztar / 137.D0 * Hbar / Rtar * (r / Rtar) ** Lambda
  End If

  result = result * (Beta4T + 9.D0 * Beta2T**2 / 7.D0 / Sqrt(Pi))

  Fct4 = result
  Return
End Function Fct4

Real(8) Function Fcp(r)
  Use ccfull_initialization_mod
  Implicit None
  Real(8), Intent(In) :: r
  Real(8) :: result, Lambda

  Lambda = LambdaP

  If (r > Rpro) Then
     result = 3.D0 / (2.D0*Lambda + 1.D0) * Zpro * Ztar / 137.D0 * Hbar / r * (Rpro / r) ** Lambda
  Else
     result = 3.D0 / (2.D0*Lambda + 1.D0) * Zpro * Ztar / 137.D0 * Hbar / Rpro * (r / Rpro) ** Lambda
  End If

  result = result * BetaP / SQRT(4.D0 * Pi)

  Fcp = result
  Return
End Function Fcp

Real(8) Function Fcp2(r)
  Use ccfull_initialization_mod
  Implicit None
  Real(8), Intent(In) :: r
  Real(8) :: result
  Integer :: Lambda

  Lambda = 2

  If (r > Rpro) Then
     result = 3.D0 / (2.D0 * Lambda + 1.D0) * Zpro * Ztar / 137.D0 * Hbar / r * (Rpro / r)**Lambda
  Else
     result = 3.D0 / (2.D0 * Lambda + 1.D0) * Zpro * Ztar / 137.D0 * Hbar / Rpro * (r / Rpro)**Lambda
  End If

  Fcp2 = result
  Return
End Function Fcp2

Real(8) Function Fcp3(r)
  Use ccfull_initialization_mod
  Implicit None
  Real(8), Intent(In) :: r
  Real(8) :: result
  Integer :: Lambda

  Lambda = 3
  result = 0.D0

  If (r > Rpro) Then
     result = 3.D0 / (2.D0 * Lambda + 1.D0) * Zpro * Ztar / 137.D0 * Hbar / r * (Rpro / r)**Lambda
  Else
     result = 3.D0 / (2.D0 * Lambda + 1.D0) * Zpro * Ztar / 137.D0 * Hbar / Rpro * (r / Rpro)**Lambda
  End If

  result = result * (Beta4P + 7.D0 * Beta2P**2 / 7.D0 / Sqrt(Pi))

  Fcp3 = result
  Return
End Function Fcp3

Real(8) Function Fcp4(r)
  Use ccfull_initialization_mod
  Implicit None
  Real(8), Intent(In) :: r
  Real(8) :: result
  Integer :: Lambda

  Lambda = 4
  result = 0.D0

  If (r > Rpro) Then
     result = 3.D0 / (2.D0 * Lambda + 1.D0) * Zpro * Ztar / 137.D0 * Hbar / r * (Rpro / r)**Lambda
  Else
     result = 3.D0 / (2.D0 * Lambda + 1.D0) * Zpro * Ztar / 137.D0 * Hbar / Rpro * (r / Rpro)**Lambda
  End If

  result = result * (Beta4P + 9.D0 * Beta2P**2 / 7.D0 / Sqrt(pi))

  Fcp4 = result
  Return
End Function Fcp4

Real(8) Function Fctt(r)
  Use ccfull_initialization_mod
  Implicit None
  Real(8), Intent(In) :: r
  Real(8) :: result
  Integer :: Lambda

  Lambda = LambdaT2
  result = 0.D0

  If (r > Rtar) Then
     result = 3.D0 / (2.D0 * Lambda + 1.D0) * Zpro * Ztar / 137.D0 * Hbar / r * (Rtar / r)**Lambda
  Else
     result = 3.D0 / (2.D0 * Lambda + 1.D0) * Zpro * Ztar / 137.D0 * Hbar / Rtar * (r / Rtar)**Lambda
  End If

  result = result * BetaT2 / SQRT(4.D0 * Pi)

  Fctt = result
  Return
End Function Fctt

Real(8) Function Fct2v(r)
  Use ccfull_initialization_mod
  Implicit None
  Real(8), Intent(In) :: r
  Real(8) :: result
  Integer :: Lambda

  Lambda = 2
  result = 0.D0

  If (r > Rtar) Then
     result = 3.D0 / (2.D0 * Lambda + 1.D0) * Zpro * Ztar / 137.D0 * Hbar / r * (Rtar / r)**Lambda
  Else
     result = 3.D0 / (2.D0 * Lambda + 1.D0) * Zpro * Ztar / 137.D0 * Hbar / Rtar * (r / Rtar)**Lambda
  End If

  Fct2v = result
  Return
End Function Fct2v

Real(8) Function Fcp2v(r)
  Use ccfull_initialization_mod
  Implicit None
  Real(8), Intent(In) :: r
  Real(8) :: result
  Integer :: Lambda

  Lambda = 2
  result = 0.D0

  If (r > Rpro) Then
     result = 3.D0 / (2.D0 * Lambda + 1.D0) * Zpro * Ztar / 137.D0 * Hbar / r * (Rpro / r)**Lambda
  Else
     result = 3.D0 / (2.D0 * Lambda + 1.D0) * Zpro * Ztar / 137.D0 * Hbar / Rpro * (r / Rpro)**Lambda
  End If

  Fcp2v = result
  Return
End Function Fcp2v

Real(8) Function Ftrans(r)
  Use ccfull_initialization_mod
  Implicit None
  Real(8), Intent(In) :: r
  Real(8) :: result, dVn_dr
  external dVn_dr

  result = 0.D0
  result = Ftr * dVn_dr(r)

  Ftrans = result
  Return
End Function Ftrans

Real(8) Function CG(j1, m1, j2, m2, j3, m3)
    ! Clebsch-Gordan coefficient <j1 m1 j2 m2 | j3 m3>
    implicit real*8 (a-h, o-z)
    external fact

    if (m1 + m2 .ne. m3) then
        cg = 0.d0
        return
    end if

    if (j3 .lt. abs(j1 - j2)) then
        cg = 0.d0
        return
    end if

    if (j3 .gt. j1 + j2) then
        cg = 0.d0
        return
    end if

    ka = j1 + j2 - j3
    kb = j3 + j1 - j2
    kc = j2 + j3 - j1
    kd = j1 + j2 + j3 + 1

    del = sqrt(fact(ka) * fact(kb) * fact(kc) / fact(kd))

    cg = 0.d0
    do n = 0, max(j1 + j2 - j3, j1 - m1, j2 + m2)
        ka1 = j1 + j2 - j3 - n
        if (ka1 .lt. 0.d0) cycle

        ka2 = j3 - j2 + m1 + n
        if (ka2 .lt. 0.d0) cycle

        ka3 = j3 - j1 - m2 + n
        if (ka3 .lt. 0.d0) cycle

        ka4 = j1 - m1 - n
        if (ka4 .lt. 0.d0) cycle

        ka5 = j2 + m2 - n
        if (ka5 .lt. 0.d0) cycle

        cg = cg + (-1.d0)**n / &
                 (fact(n) * fact(ka1) * fact(ka2) * fact(ka3) * fact(ka4) * fact(ka5))
    end do

    cg = cg * sqrt(fact(j1 + m1) * fact(j1 - m1))
    cg = cg * sqrt(fact(j2 + m2) * fact(j2 - m2))
    cg = cg * sqrt(fact(j3 + m3) * fact(j3 - m3))

    cg = cg * sqrt(2.d0 * j3 + 1.d0) * del

    return
End Function CG

Subroutine mdiag(a,n,d,v)

	implicit real*8(a-h,o-z)
	parameter(nmax=100)
  parameter (nlevelmax=30)
  dimension a(nlevelmax,nlevelmax)
	dimension b(nmax),z(nmax),d(nlevelmax),v(nlevelmax,nlevelmax)

	do 12 ip=1,n
	do 11 iq=1,n
	v(ip,iq)=0.d0
 11	continue
	v(ip,ip)=1.d0
 12	continue

	do 13 ip=1,n
	b(ip)=a(ip,ip)
	d(ip)=b(ip)
	z(ip)=0.d0
 13	continue

	nrot=0

	do 24 i=1,50
	sm=0.d0

	do 15 ip=1,n-1
	do 14 iq=ip+1,n
	sm=sm+abs(a(ip,iq))
 14	continue
 15     continue

	if(sm.eq.0.d0) return

	if(i.lt.4) then
	tresh=0.2d0*sm/n**2
	else
	tresh=0.d0
	endif

	do 22 ip=1,n-1
	do 21 iq=ip+1,n
	g=100.d0*abs(a(ip,iq))
	
	if((i.gt.4).and.(abs(d(ip))+g.eq.abs(d(ip))).and.(abs(d(iq))+g.eq.abs(d(iq)))) then
	   a(ip,iq)=0.d0
	elseif(abs(a(ip,iq)).gt.tresh) then
	   h=d(iq)-d(ip)

	   if(abs(h)+g.eq.abs(h)) then
	     t=a(ip,iq)/h
	   else
	     theta=0.5d0*h/a(ip,iq)
	     t=1.d0/(abs(theta)+sqrt(1.d0+theta**2))
	     if(theta.lt.0.d0) t=-t
	   endif

	c=1.d0/sqrt(1.d0+t**2)
	s=t*c
	tau=s/(1.d0+c)
	h=t*a(ip,iq)
	z(ip)=z(ip)-h
	z(iq)=z(iq)+h
	d(ip)=d(ip)-h
	d(iq)=d(iq)+h
	a(ip,iq)=0.d0

	do 16 j=1,ip-1
	g=a(j,ip)
	h=a(j,iq)
	a(j,ip)=g-s*(h+g*tau)
	a(j,iq)=h+s*(g-h*tau)
 16	continue

	do 17 j=ip+1,iq-1
	g=a(ip,j)
	h=a(j,iq)
	a(ip,j)=g-s*(h+g*tau)
	a(j,iq)=h+s*(g-h*tau)
 17	continue

	do 18 j=iq+1,n
	g=a(ip,j)
	h=a(iq,j)
	a(ip,j)=g-s*(h+g*tau)
	a(iq,j)=h+s*(g-h*tau)
 18	continue

	do 19 j=1,n
	g=v(j,ip)
	h=v(j,iq)
	v(j,ip)=g-s*(h+g*tau)
	v(j,iq)=h+s*(g-h*tau)
 19	continue
	
	nrot=nrot+1
	endif

 21	continue
 22	continue

	do 23 ip=1,n
	b(ip)=b(ip)+z(ip)
	d(ip)=b(ip)
	z(ip)=0.d0

 23	continue
 24	continue

	return

End Subroutine  mdiag

SUBROUTINE matinv(nmax, c, d)
  IMPLICIT NONE
  
  ! Argument declarations
  INTEGER, INTENT(IN) :: nmax
  COMPLEX(8), DIMENSION(30,30), INTENT(INOUT) :: c, d
  
  ! Parameter and local variable declarations
  COMPLEX(8) :: u(30), v(30), t, a, b
  REAL(8) :: deter
  INTEGER :: m, n, j, k, l
  
  ! Initialize
  deter = 1.0_8
  d = (0.0_8, 0.0_8)
  DO m = 1, nmax
      d(m, m) = (1.0_8, 0.0_8)
  END DO
  
  ! Main inversion loop
  DO n = 1, nmax
      t = c(n, n)
      
      ! Check for zero pivot element
      IF (ABS(t) < 1.0E-10_8) THEN
          j = n
          DO
              j = j + 1
              IF (j > nmax) THEN
                  PRINT *, 'matrix not invertible'
                  RETURN
              END IF
              
              t = c(n, j)
              deter = -deter
              IF (ABS(t) > 1.0E-10_8) EXIT
          END DO
          
          ! Swap rows
          DO k = 1, nmax
              u(k) = c(n, k)
              v(k) = d(n, k)
              c(n, k) = c(j, k)
              d(n, k) = d(j, k)
              c(j, k) = u(k)
              d(j, k) = v(k)
          END DO
      END IF
      
      ! Eliminate column
      DO k = 1, nmax
          IF (k == n) CYCLE
          
          a = c(k, n) / c(n, n)
          DO l = 1, nmax
              c(k, l) = c(k, l) - a * c(n, l)
              d(k, l) = d(k, l) - a * d(n, l)
          END DO
      END DO
      
      ! Normalize row
      b = c(n, n)
      deter = b * deter
      DO m = 1, nmax
          c(n, m) = c(n, m) / b
          d(n, m) = d(n, m) / b
      END DO
  END DO
  
END SUBROUTINE matinv

subroutine matinv_lapack(n, A, Ainv)
  implicit none
  integer, intent(in) :: n
  complex*16, intent(in) :: A(n,n)
  complex*16, intent(out) :: Ainv(n,n)

  ! 局部变量
  integer :: info, lwork
  integer, allocatable :: ipiv(:)
  complex*16, allocatable :: work(:)
  complex*16 :: tmp(n,n)

  ! 拷贝 A 到临时变量
  tmp = A
  Ainv = (0.0d0, 0.0d0)

  allocate(ipiv(n))

  ! LU 分解：tmp 会被覆盖为 LU
  call zgetrf(n, n, tmp, n, ipiv, info)
  if (info /= 0) then
    print *, "Error in zgetrf: info = ", info
    stop
  end if

  ! 工作空间大小（可调优）
  lwork = n * n
  allocate(work(lwork))

  ! 计算逆矩阵
  call zgetri(n, tmp, n, ipiv, work, lwork, info)
  if (info /= 0) then
    print *, "Error in zgetri: info = ", info
    stop
  end if

  Ainv = tmp

  deallocate(ipiv)
  deallocate(work)
end subroutine matinv_lapack

Real(8) Function fact(n)
    Implicit none
    Integer n, i

    If (n .lt. 0) Then
        fact = 0.d0
        Return
    End If

    If (n .eq. 0) Then
        fact = 1.d0
        Return
    End If

    fact = 1.d0
    Do i = 1, n
        fact = fact * i * 1.d0
    End Do

    Return
End Function fact

Subroutine grotation()
    Use ccfull_initialization_mod
    Implicit None
    Integer :: it, jt, ip, jp, i
    Character(len=1) :: ans
    Real(8) :: erott0, erotp0

    Do it = 0, Ntar
      Do jt = 0, Ntar
        bett(it,jt) = 0.0D0
        If (it == jt - 1 .Or. it == jt + 1 .Or. it == jt) bett(it,jt) = Beta2T
      End Do
      erott(it) = 2.0D0 * it * (2.0D0 * it + 1.0D0) / 6.0D0 * E2T
    End Do

    If (IVIBROTT == 1 .And. Ntar > 1) Then
      !Print *, ''
      !Print *, 'Generalised E2 couplings for the target rotor (n/y)?'
      !Read(*,'(A)') ans
      ans = 'n'
      If (ans == 'y' .Or. ans == 'Y') Then
        Print *, 'Beta2_targ=', Beta2T
        Do it = 1, Ntar - 1
          jt = it + 1
          Print *, 'Modify beta2 for transition from', 2*it, '+ to', 2*jt, '+ ? (n/y)'
          Read(*,'(A)') ans
          If (ans == 'y' .Or. ans == 'Y') Then
            Print *, 'BETA2=?'
            Read(*,*) bett(it,jt)
            bett(jt,it) = bett(it,jt)
            Print *, 'Beta2 set to: ', bett(it,jt)
          End If
        End Do

        Print *, 'Reorientation terms:'
        Do it = 1, Ntar
          jt = it
          Print *, 'Modify beta2 for transition from', 2*it, '+ to', 2*jt, '+ ? (n/y)'
          Read(*,'(A)') ans
          If (ans == 'y' .Or. ans == 'Y') Then
            Print *, 'BETA2=?'
            Read(*,*) bett(it,jt)
            bett(jt,it) = bett(it,jt)
          End If
        End Do

        Print *, 'Excitation energies:'
        Do i = 2, Ntar
          Print *, 'Energy of the', 2*i, '+ state for pure rotor =', erott(i)
          Print *, 'Modify this energy? (n/y)'
          Read(*,'(A)') ans
          If (ans == 'y' .Or. ans == 'Y') Then
            erott0 = erott(i)
            Print *, 'Energy=?'
            Read(*,*) erott(i)
            Print *, 'Modified from', erott0, ' to ', erott(i)
          End If
        End Do
      End If
    End If

    Do ip = 0, Npro
      Do jp = 0, Npro
        betp(ip,jp) = 0.0D0
        If (ip == jp - 1 .Or. ip == jp + 1 .Or. ip == jp) betp(ip,jp) = Beta2P
      End Do
      erotp(ip) = 2.0D0 * ip * (2.0D0 * ip + 1.0D0) / 6.0D0 * E2P
    End Do

    If (IVIBROTP == 1 .And. Npro > 1) Then
      !Print *, 'Generalised E2 couplings for the projectile rotor (n/y)?'
      !Read(*,'(A)') ans
      ans = 'n'
      If (ans == 'y' .Or. ans == 'Y') Then
        Print *, 'Beta2_proj=', Beta2P
        Do ip = 1, Npro - 1
          jp = ip + 1
          Print *, 'Modify beta2 for transition from', 2*ip, '+ to', 2*jp, '+ ? (n/y)'
          Read(*,'(A)') ans
          If (ans == 'y' .Or. ans == 'Y') Then
            Print *, 'BETA2=?'
            Read(*,*) betp(ip,jp)
            betp(jp,ip) = betp(ip,jp)
          End If
        End Do

        Print *, 'Reorientation terms:'
        Do ip = 1, Npro
          jp = ip
          Print *, 'Modify beta2 for transition from', 2*ip, '+ to', 2*jp, '+ ? (n/y)'
          Read(*,'(A)') ans
          If (ans == 'y' .Or. ans == 'Y') Then
            Print *, 'BETA2=?'
            Read(*,*) betp(ip,jp)
            betp(jp,ip) = betp(ip,jp)
          End If
        End Do

        Print *, 'Excitation energies:'
        Do i = 2, Npro
          Print *, 'Energy of the', 2*i, '+ state for pure rotor =', erotp(i)
          Print *, 'Modify this energy? (n/y)'
          Read(*,'(A)') ans
          If (ans == 'y' .Or. ans == 'Y') Then
            erotp0 = erotp(i)
            Print *, 'Energy=?'
            Read(*,*) erotp(i)
            Print *, 'Modified from', erotp0, ' to ', erotp(i)
          End If
        End Do
      End If
    End If

End Subroutine grotation

Subroutine Mutual()
    Use ccfull_initialization_mod
    Implicit None
    Integer :: i, j
    Character(len=1) :: ans, sign1, sign2
    Integer :: imut

    Do i = 0, Nlevelmax
      Do j = 0, Nlevelmax
        imutual(i,j) = 0
      End Do
    End Do

    Print *, ''
    Print *, 'Mutual excitations in the *target* nucleus'
    Print *, 'Include the mutual excitations (y/n)?'
    Read(*,'(A)') ans

    If (ans == 'n' .Or. ans == 'N') Then
      imut = 0
      Nlevel = (Ntar + NphononT2 + 1) * (Npro + 1)
      Do i = 0, Ntar
        Do j = 0, NphononT2
          If (i == 0 .Or. j == 0) imutual(i,j) = 1
        End Do
      End Do
      Return
    End If

    Print *, 'All the possible mutual excitation channels (n/y)?'
    Read(*,'(A)') ans

    If (ans == 'y' .Or. ans == 'Y') Then
      imut = 1
      Nlevel = (Ntar + 1) * (NphononT2 + 1) * (Npro + 1)
      Do i = 0, Ntar
        Do j = 0, NphononT2
          imutual(i,j) = 1
        End Do
      End Do
      Return
    End If

    imut = 2
    If (Mod(LambdaT, 2) == 0) Then
      sign1 = '+'
    Else
      sign1 = '-'
    End If
    If (Mod(LambdaT2, 2) == 0) Then
      sign2 = '+'
    Else
      sign2 = '-'
    End If

    Do i = 0, Ntar
      Do j = 0, NphononT2
        If (i == 0 .Or. j == 0) Then
          imutual(i,j) = 1
        Else
          Print *, 'Include (', LambdaT, sign1, '^', i, ',', LambdaT2, sign2, '^', j, ') state ? (y/n)'
          Read(*,'(A)') ans
          If (ans == 'n' .Or. ans == 'N') Then
            imutual(i,j) = 0
          Else
            imutual(i,j) = 1
          End If
        End If
      End Do
    End Do

    Print *, 'Excited states in the target to be included:'
    Nlevel = 0

    Do i = 0, Ntar
      Do j = 0, NphononT2
        If (imutual(i,j) == 0) Cycle
        Print *, '(', LambdaT, sign1, '^', i, ',', LambdaT2, sign2, '^', j, ') state'
        Nlevel = Nlevel + 1
      End Do
    End Do

    Nlevel = Nlevel * (Npro + 1)

End Subroutine Mutual

Subroutine Anharmonicity()
  Use ccfull_initialization_mod
  Implicit None

  Character(1) :: ans, sign
  Integer :: it, jt, i

  ! Target: AHV Couplings (1st mode)
  Do it = 0, Ntar
    Do jt = 0, Ntar
      betnahv(it, jt) = 0.0
      betcahv(it, jt) = 0.0

      If (it == jt - 1) Then
        betnahv(it, jt) = BetaTn * Sqrt(Real(jt))
        betcahv(it, jt) = BetaT  * Sqrt(Real(jt))
      Else If (jt == it - 1) Then
        betnahv(it, jt) = BetaTn * Sqrt(Real(it))
        betcahv(it, jt) = BetaT  * Sqrt(Real(it))
      End If
    End Do
    omeahv(it) = it * OmegaT
  End Do

  If (IVIBROTT == 0 .And. Ntar > 1) Then
    Write(*, '(A)') 'AHV couplings for the first mode in the target phonon (n/y)?'
    Read(*, '(A1)') ans
    If (ans == 'Y' .Or. ans == 'y') Then
      Write(*, '(A)') '**** AHV Couplings in the target (the 1st mode)'
      If ((-1.)**LambdaT ==  1) sign = '+'
      If ((-1.)**LambdaT == -1) sign = '-'

      Do it = 0, Ntar
        Do jt = 0, Ntar
          If (it > jt)Cycle
          If ((it == 0 .And. jt <= 1)) Cycle
          Write(*, *) ' '
          Write(*, '(A,I2,2A,I2,A,I2,2A,I2,A)') 'Transition from the', LambdaT, sign, '^', it, ' to the',&
                                                 LambdaT, sign, '^', jt, 'state:'
          Write(*, '(A,2G15.5)') '   beta_N and beta_C in the HO limit=', betnahv(it, jt), betcahv(it, jt)
          Write(*, '(2A)') '    Modify these beta_N a/o beta_C (n/y)?'
          Read(*, '(A1)') ans
          If (ans == 'Y' .Or. ans == 'y') Then
            Write(*, *) '   beta_N and beta_C =?'
            Read(*, *) betnahv(it, jt), betcahv(it, jt)
            Write(*, '(G15.5, G15.5)') betnahv(it, jt), betcahv(it, jt)
          End If
        End Do
      End Do

      Write(*, *) 'Excitation energy for the first mode:'
      Do i = 2, Ntar
        Write(*, '(A,I2,A,G15.5)') '    Energy of the', i, '-phonon state in the HO=', i * OmegaT
        Write(*, *) '   Modify this energy(n/y)?'
        Read(*, '(A1)') ans
        If (ans == 'Y' .Or. ans == 'y') Then
          Write(*, *) '   Energy=?'
          Read(*, *) omeahv(i)
          Write(*, '(A,I2,A,G15.5,A,G12.5,A)') 'Energy of the', i, '-phonon state=', omeahv(i), '(HO: ', i * OmegaT, ')'
        End If
      End Do
    End If
  End If


  ! Target: AHV Couplings (2st mode)
  Do it = 0, NphononT2
    Do jt = 0, NphononT2
      betnahv2(it, jt) = 0.0
      betcahv2(it, jt) = 0.0

      If (it == jt - 1) Then
        betnahv2(it, jt) = BetaT2n * Sqrt(Real(jt))
        betcahv2(it, jt) = BetaT2  * Sqrt(Real(jt))
      Else If (jt == it - 1) Then
        betnahv(it, jt) = BetaT2n * Sqrt(Real(it))
        betcahv(it, jt) = BetaT2  * Sqrt(Real(it))
      End If
    End Do
    omeahv2(it) = it * OmegaT2
  End Do

  If (NphononT2 > 1) Then
    Write(*, '(A)') 'AHV couplings for the second mode in the target phonon (n/y)?'
    Read(*, '(A1)') ans
    If (ans == 'Y' .Or. ans == 'y') Then
      Write(*, '(A)') '**** AHV Couplings in the target (the 1st mode)'
      If ((-1.)**LambdaT2 ==  1) sign = '+'
      If ((-1.)**LambdaT2 == -1) sign = '-'

      Do it = 0, NphononT2
        Do jt = 0, NphononT2
          If (it > jt)Cycle
          If ((it == 0 .And. jt <= 1)) Cycle
          Write(*, *) ' '
          Write(*, '(A,I2,2A,I2,A,I2,2A,I2,A)') 'Transition from the', LambdaT2, sign, '^', it, ' to the',&
                                               LambdaT2, sign, '^', jt, 'state:'
          Write(*, '(A,2G15.5)') '   beta_N and beta_C in the HO limit=', betnahv2(it, jt), betcahv2(it, jt)
          Write(*, '(2A)') '    Modify these beta_N a/o beta_C (n/y)?'
          Read(*, '(A1)') ans
          If (ans == 'Y' .Or. ans == 'y') Then
            Write(*, *) '   beta_N and beta_C =?'
            Read(*, *) betnahv2(it, jt), betcahv2(it, jt)
            Write(*, '(G15.5, G15.5)') betnahv2(it, jt), betcahv2(it, jt)
          End If
        End Do
      End Do

      Write(*, *) 'Excitation energy for the second mode:'
      Do i = 2, Ntar
        Write(*, '(A,I2,A,G15.5)') '    Energy of the', i, '-phonon state in the HO=', i * OmegaT2
        Write(*, *) '   Modify this energy(n/y)?'
        Read(*, '(A1)') ans
        If (ans == 'Y' .Or. ans == 'y') Then
          Write(*, *) '   Energy=?'
          Read(*, *) omeahv2(i)
          Write(*, '(A,I2,A,G15.5,A,G12.5,A)') 'Energy of the', i, '-phonon state=', omeahv2(i), '(HO: ', i * OmegaT2, ')'
        End If
      End Do
    End If
  End If


  Do it = 0, Npro
    Do jt = 0, Npro
      betnahvp(it, jt) = 0.0
      betcahvp(it, jt) = 0.0

      If (it == jt - 1) Then
        betnahvp(it, jt) = BetaPn * Sqrt(Real(jt))
        betcahvp(it, jt) = BetaP  * Sqrt(Real(jt))
      Else If (jt == it - 1) Then
        betnahvp(it, jt) = BetaPn * Sqrt(Real(it))
        betcahvp(it, jt) = BetaP  * Sqrt(Real(it))
      End If
    End Do
    omeahvp(it) = it * OmegaP
  End Do

  If (IVIBROTP == 0 .And. Npro > 1) Then
      Write(*, '(A)') 'AHV couplings for the first mode in the projectile phonon (n/y)?'
      Read(*, '(A1)') ans
      If (ans == 'Y' .Or. ans == 'y') Then
        Write(*, '(A)') '**** AHV Couplings in the projectile'
        If ((-1.)**LambdaP ==  1) sign = '+'
        If ((-1.)**LambdaP == -1) sign = '-'

        Do it = 0, Npro
          Do jt = 0, Npro
            If (it > jt)Cycle
            If ((it == 0 .And. jt <= 1)) Cycle
            Write(*, *) ' '
            Write(*, '(A,I2,2A,I2,A,I2,2A,I2,A)') 'Transition from the', LambdaP, sign, '^', it, ' to the',&
                                                 LambdaP, sign, '^', jt, 'state:'
            Write(*, '(A,2G15.5)') '   beta_N and beta_C in the HO limit=', betnahvp(it, jt), betcahvp(it, jt)
            Write(*, '(2A)') '    Modify these beta_N a/o beta_C (n/y)?'
            Read(*, '(A1)') ans
            If (ans == 'Y' .Or. ans == 'y') Then
              Write(*, *) '   beta_N and beta_C =?'
              Read(*, *) betnahvp(it, jt), betcahvp(it, jt)
              Write(*, '(G15.5, G15.5)') betnahvp(it, jt), betcahvp(it, jt)
            End If
          End Do
        End Do

        Write(*, *) 'Excitation energy for the first mode:'
        Do i = 2, Npro
          Write(*, '(A,I2,A,G15.5)') '    Energy of the', i, '-phonon state in the HO=', i * OmegaP
          Write(*, *) '   Modify this energy(n/y)?'
          Read(*, '(A1)') ans
          If (ans == 'Y' .Or. ans == 'y') Then
            Write(*, *) '   Energy=?'
            Read(*, *) omeahvp(i)
            Write(*, '(A,I2,A,G15.5,A,G12.5,A)') 'Energy of the', i, '-phonon state=', omeahvp(i), '(HO: ', i * OmegaP, ')'
          End If
        End Do
      End If
    End If

End Subroutine Anharmonicity

Subroutine coupled_matrix0()
  Use ccfull_initialization_mod
  implicit none
  real(8) :: A(Nlevelmax, Nlevelmax)
  Real(8) :: C, CG
  integer :: i, j, ip, it, it2, jp, jt, jt2
  external CG
  
  !initialization of A matrix
  Do i = 1, Nlevelmax
    Do j = 1, Nlevelmax
     A(i,j) = 0.0d0
    End Do
  End Do

  If((Nlevel - Ntrans) == 1)Return
  i = 0
  Do ip = 0, Npro
    Do it = 0, Ntar
      Do it2 = 0, NphononT2

        if (NphononT2 .Ne. 0)Then
          if(imutual(it, it2) == 0)Cycle
        End if

        i = i + 1 
        j = 0
        Do jp = 0, Npro
          Do jt = 0, Ntar
            Do jt2 = 0, NphononT2

              if (NphononT2 .Ne. 0)Then
                if(imutual(jt, jt2) == 0)Cycle
              End if

              j = j + 1 

              If (i > j)Then
                A(i,j) = A(j,i)
                Cycle
              End If

              C = 0.0d0

              if (ip == jp .and. it2 == jt2) then
                  if (IVIBROTT == 0) then
                      C = Rtar * betnahv(it, jt) / sqrt(4.0d0 * pi)

                  else
                      C = Rtar * bett(it, jt) * sqrt((2 * 2 * it + 1) * 5.0d0 * (2 * 2 * jt + 1) / (4.0d0 * pi)) &
                          * CG(2 * it, 0, 2, 0, 2 * jt, 0)**2 / (2.0d0 * 2 * jt + 1)

                      C = C + Rtar * Beta4T * sqrt((2 * 2 * it + 1) * 9.0d0 * (2 * 2 * jt + 1) / (4.0d0 * pi)) &
                          * CG(2 * it, 0, 4, 0, 2 * jt, 0)**2 / (2.0d0 * 2 * jt + 1)
                  end if
              end if

              if (ip == jp .and. it == jt) then
                  C = C + Rtar * betnahv2(it2, jt2) / sqrt(4.0d0 * pi)
              end if

              if (it == jt .and. it2 == jt2) then
                  if (IVIBROTP == 0) then
                      C = C + Rpro * betnahvp(ip, jp) / sqrt(4.0d0 * pi)
                  else
                      C = C + Rpro * betp(ip, jp) * sqrt((2 * 2 * ip + 1) * 5.0d0 * (2 * 2 * jp + 1) / (4.0d0 * pi)) &
                          * CG(2 * ip, 0, 2, 0, 2 * jp, 0)**2 / (2.0d0 * 2 * jp + 1)

                      C = C + Rpro * Beta4P * sqrt((2 * 2 * ip + 1) * 9.0d0 * (2 * 2 * jp + 1) / (4.0d0 * pi)) &
                          * CG(2 * ip, 0, 4, 0, 2 * jp, 0)**2 / (2.0d0 * 2 * jp + 1)
                  end if
              end if

              A(i,j) = C


            End Do
          End Do
        End Do

      End Do
    End Do
  End Do

  Call mdiag(A, Nlevel - Ntrans, Ev, Evec)

  Return

End Subroutine coupled_matrix0

Subroutine coupled_matrix(r, cpot_matrix)
  Use ccfull_initialization_mod
  Implicit None
  Real(8), Intent(In) :: r
  real(8), intent(out) :: cpot_matrix(Nlevelmax, Nlevelmax)
  Real(8) :: C, A, C0
  Real(8) :: Fct, Fct2, Fct3, Fct4, Fct2v, Fctt
  Real(8) :: Ftrans, Fcp, Fcp2, Fcp3, Fcp4, Fcp2v
  Real(8) :: CG, Vn, VnCC
  integer :: i, j, k, ip, it, it2, jp, jt, jt2
  external CG, Fct, Fct2, Fct3, Fct4, Fct2v, Fctt
  external Ftrans, Fcp, Fcp2, Fcp3, Fcp4, Fcp2v, Vn, VnCC
  
  Do i = 1, nlevelmax
    eps(i) = 0.
  End Do

  Do i = 1, Nlevelmax
    Do j = 1, Nlevelmax
     cpot_matrix(i,j) = 0.0d0
    End Do
  End Do

  If((Nlevel - Ntrans) == 1)Return

  i = 0
  Do ip = 0, Npro
    Do it = 0, Ntar
      Do it2 = 0, NphononT2

        if (NphononT2 .Ne. 0)Then
          if(imutual(it, it2) == 0.)Cycle
        End if

        i = i + 1 
        j = 0
        Do jp = 0, Npro
          Do jt = 0, Ntar
            Do jt2 = 0, NphononT2

            if (NphononT2 .Ne. 0)Then
              if(imutual(jt, jt2) == 0.)Cycle
            End if

            j = j + 1 

            If (i > j)Then

              cpot_matrix(i,j) = cpot_matrix(j,i)
              Cycle
            End If

            C = 0.0d0
            ! nuclear coupling
            Do k = 1, Nlevel-Ntrans
              C = C + VnCC(r, Ev(k))*Evec(i, k)*Evec(j, k)

            End Do

            ! ---- Target contribution ----
            if (ip == jp .and. it2 == jt2) then
                if (IVIBROTT == 0) then
                    if (it /= jt) then
                        A = betcahv(it, jt) * Fct(r)
                        if (BetaT /= 0.d0) A = A / BetaT
                        
                    else
                        A = betcahv(it, jt) * Fct2v(r) / sqrt(4.d0 * pi)

                    end if
                    C = C + A
                else
                    C = C + sqrt((2 * 2 * it + 1) * 5.d0 * (2 * 2 * jt + 1) / (4.d0 * pi)) &
                          * CG(2 * it, 0, 2, 0, 2 * jt, 0)**2 / (2.d0 * 2 * jt + 1) * Fct2(r) &
                          * (bett(it, jt) + 2.d0 * sqrt(5.d0 / pi) * bett(it, jt)**2 / 7.d0)

                    C = C + sqrt((2 * 2 * it + 1) * 9.d0 * (2 * 2 * jt + 1) / (4.d0 * pi)) &
                          * CG(2 * it, 0, 4, 0, 2 * jt, 0)**2 / (2.d0 * 2 * jt + 1) * Fct4(r)
                end if
            end if

            ! ---- Projectile contribution ----
            if (it == jt .and. it2 == jt2) then
                if (IVIBROTP == 0) then
                    if (ip /= jp) then
                        A = betcahvp(ip, jp) * Fcp(r)
                        if (BetaP /= 0.d0) A = A / BetaP
                    else
                        A = betcahvp(ip, jp) * Fcp2v(r) / sqrt(4.d0 * pi)
                    end if
                    C = C + A
                else
                    C = C + sqrt((2 * 2 * ip + 1) * 5.d0 * (2 * 2 * jp + 1) / (4.d0 * pi)) &
                          * CG(2 * ip, 0, 2, 0, 2 * jp, 0)**2 / (2.d0 * 2 * jp + 1) * Fcp2(r) &
                          * (betp(ip, jp) + 2.d0 * sqrt(5.d0 / pi) * betp(ip, jp)**2 / 7.d0)

                    C = C + sqrt((2 * 2 * ip + 1) * 9.d0 * (2 * 2 * jp + 1) / (4.d0 * pi)) &
                          * CG(2 * ip, 0, 4, 0, 2 * jp, 0)**2 / (2.d0 * 2 * jp + 1) * Fcp4(r)
                end if
            end if

            ! ---- BetaT2 part (vibrational excitation of target) ----
            if (ip == jp .and. it == jt) then
                if (it2 /= jt2) then
                    A = betcahv2(it2, jt2) * Fctt(r)
                    if (BetaT2 /= 0.d0) A = A / BetaT2
                else
                    A = betcahv2(it2, jt2) * Fct2v(r) / sqrt(4.d0 * pi)
                end if
                C = C + A
            end if

            ! ---- Excitation energy shift ----
            If (it == jt .and. ip == jp .and. it2 == jt2) then
                if (IVIBROTT == 0) then
                    C = C + omeahv(it)
                    eps(i) = eps(i) + omeahv(it)
                else
                    C = C + erott(it)
                    eps(i) = eps(i) + erott(it)
                end if

                if (IVIBROTP == 0) then
                    C = C + omeahvp(ip)
                    eps(i) = eps(i) + omeahvp(ip)
                else
                    C = C + erotp(ip)
                    eps(i) = eps(i) + erotp(ip)
                end if

                C = C + omeahv2(it2)
                eps(i) = eps(i) + omeahv2(it2)
            End if

            cpot_matrix(i, j) = C


            End Do
          End Do
        End Do

      End Do
    End Do
  End Do

  C0 = cpot_matrix(1,1)
  Do i = 1, Nlevel-Ntrans
    cpot_matrix(i,i) = cpot_matrix(i,i) - Vn(r)
    !cpot_matrix(i,i) = cpot_matrix(i,i) - C0
  End Do

  ! transfer coupling
  If (Ntrans == 1) Then
      cpot_matrix(1, Nlevel) = Ftrans(r)
      cpot_matrix(Nlevel, 1) = Ftrans(r)
      cpot_matrix(Nlevel, Nlevel) = cpot_matrix(1, 1) - Qtrans + (Zpro + iq) * (Ztar - iq) &
                                        / r * Hbar / 137.0 - Zpro * Ztar / r * Hbar / 137.0
  End If

  Return

End Subroutine coupled_matrix

Subroutine rkutta00(psi0, phi0, psi1)
  Use ccfull_initialization_mod
  Implicit None
  complex*16, Intent(In)  ::  psi0(Nlevelmax), phi0(Nlevelmax)
  complex*16, intent(out) ::  psi1(Nlevelmax)
  complex*16 :: ai
  complex*16, dimension(Nlevelmax) :: ak1, ak2, ak3, ak4
  complex*16, dimension(Nlevelmax) :: bk1, bk2, bk3, bk4

  Real(8) ::  fac, r, rh, rpp, V
  Integer ::  j1, i0, ic, is
  external V

  ai = (0.d0, 1.d0)
  fac = dr * (2.d0 * ReduceMass / Hbar**2)

  r  = R_min
  rh = R_min + dr/2.d0
  rpp= R_min + dr
  Do j1 = 1, Nlevel
      ak1(j1) = 0.0D0
      ak2(j1) = 0.0D0
      ak3(j1) = 0.0D0
      ak4(j1) = 0.0D0
      bk1(j1) = 0.0D0
      bk2(j1) = 0.0D0
      bk3(j1) = 0.0D0
      bk4(j1) = 0.0D0
  End Do

  Do i0 = 1, Nlevel
      Do ic = 1, Nlevel
        ak1(i0) = ak1(i0) + fac * CPOT(i0, ic, 0) * psi0(ic)
      End Do
      ak1(i0) = ak1(i0) - fac * (E - V(r, L_i)) * psi0(i0)
      bk1(i0) = dr * phi0(i0)
  End Do

  Do i0 = 1, Nlevel
      Do ic = 1, Nlevel
        ak2(i0) = ak2(i0) + fac * CPOTH(i0, ic) * (psi0(ic) + 1.0D0 / 2.0D0 * bk1(ic))
      End Do
      ak2(i0) = ak2(i0) - fac * (E - V(rh, L_i)) * (psi0(i0) + 1.0D0 / 2.0D0 * bk1(i0))
      bk2(i0) = dr * (phi0(i0) + 1.0D0 / 2.0D0 * ak1(i0))
  End Do

  Do i0 = 1, Nlevel
      Do ic = 1, Nlevel
        ak3(i0) = ak3(i0) + fac * CPOTH(i0, ic) * (psi0(ic) + 1.0D0 / 2.0D0 * bk2(ic))
      End Do
      ak3(i0) = ak3(i0) - fac * (E - V(rh, L_i)) * (psi0(i0) + 1.0D0 / 2.0D0 * bk2(i0))
      bk3(i0) = dr * (phi0(i0) + 1.0D0 / 2.0D0 * ak2(i0))
  End Do

  Do i0 = 1, Nlevel
      Do ic = 1, Nlevel
        ak4(i0) = ak4(i0) + fac * CPOT(i0, ic, 1) * (psi0(ic) + bk3(ic))
      End Do
      ak4(i0) = ak4(i0) - fac * (E - V(rpp, L_i)) * (psi0(i0) + bk3(i0))
      bk4(i0) = dr * (phi0(i0) + ak3(i0))
  End Do

  Do is = 1, nlevel
      psi1(is) = psi0(is) + (1.0D0 / 6.0D0) * (bk1(is) + 2.0D0 * bk2(is) + 2.0D0 * bk3(is) + bk4(is))
  End Do

  Return

End Subroutine rkutta00

Subroutine stabilize(xi1, xi, aa, ir1)
  Use ccfull_initialization_mod
  Implicit None
  Integer, Intent(In)  ::ir1
  complex*16, Intent(Inout)  ::  aa(Nlevelmax, Nlevelmax)
  complex*16, Intent(Inout)  ::  xi(Nlevelmax, Nlevelmax)
  complex*16, Intent(Inout)  ::  xi1(Nlevelmax, Nlevelmax)
  complex*16 :: ai
  complex*16, dimension(Nlevelmax, Nlevelmax)  :: psi
  complex*16, dimension(Nlevelmax, Nlevelmax)  :: cc, cin
  complex*16, dimension(Nlevelmax, Nlevelmax)  :: aa0
  complex*16, dimension(Nlevelmax, Nlevelmax)  :: xid, xid1
  Real(8) ::  r, V, fac
  Integer ::  i, j, k, i0, ic, ich
  external  V

  Do i = 1, nlevel
      Do j = 1, nlevel
        aa0(i, j)  = aa(i, j)
        xid(i, j)  = xi(i, j)
        xid1(i, j) = xi1(i, j)
      End Do
  End Do

  ai = (0.d0, 1.d0)
  fac = dr**2 * (2.d0 * ReduceMass / Hbar**2)
  r = R_min + ir1*dr

  Do i0 = 1, Nlevel
    Do ic = 1, Nlevel
        cc(i0, ic) = -fac / 12.0D0 * CPOT(i0, ic, ir1)
        If (i0 == ic) Then
          cc(i0, ic) = cc(i0, ic) - fac / 12.0D0 * (V(r, L_i) - E) + 1.0D0
        End If
    End Do
  End Do

  Call matinv(Nlevel,cc,cin)

  Do ich = 1, nlevel
      Do i0 = 1, nlevel
        psi(i0, ich) = 0.0D0
        Do ic = 1, nlevel
            psi(i0, ich) = psi(i0, ich) + cin(i0, ic) * xi(ic, ich)
        End Do
      End Do
  End Do

  Do i = 1, nlevel
      Do j = 1, nlevel
        aa(i, j) = 0.0D0
        Do k = 1, nlevel
            aa(i, j) = aa(i, j) + psi(i, k) * aa0(k, j)
        End Do
      End Do
  End Do

  call matinv(Nlevel,psi,cin)

  Do i = 1, Nlevel
      Do j = 1, Nlevel
        xi(i, j)  = 0.0D0
        xi1(i, j) = 0.0D0
        Do k = 1, Nlevel
            xi(i, j)  = xi(i, j) + xid(i, k)  * cin(k, j)
            xi1(i, j) = xi1(i, j) + xid1(i, k) * cin(k, j)
        End Do
      End Do
  End Do

  Return
End Subroutine stabilize

Subroutine Prameterize_H(H2, dbpsi,dbpsi0)
  Use ccfull_initialization_mod
  implicit none
  !H2, dbpsi, dbpsi0
  complex*16 H1(Nx_cc,Nx_cc)
  complex*16 H2(Nx1, Nx2)
  complex*16 dbpsi(Nx2, Nx1)
  complex*16 dbpsi0(Nx2)
  complex*16 ai

  Integer ::  i, j, ir, ie, il, idx, jdx
  Real(8) ::  r,rh, V, eta, rho, k, ec
  Real(8), dimension(0:200) :: fcw, gcw, fpcw, gpcw
  Real(8), dimension(0:200) :: sigmad
  Integer, dimension(0:200) :: iexp
  external V
  H1 = (0.d0, 0.d0)
  H2 = (0.d0, 0.d0)
  ai = (0.d0,1.d0)

  !BetaT = para
  ! contruct the cpot
  Call coupled_matrix0

  rh = R_min + dr / 2.0
  Call coupled_matrix(rh, CPOTH)

  Do ir = 0, R_iterat + 1
    r = R_min + dr * ir

    Call coupled_matrix(r, CPOT0)
    Do  i = 1, Nlevel
      Do  j = 1, Nlevel
        CPOT(i, j, ir) = CPOT0(i, j)

        ! contruct the H, potential term
        If(ir .lt. R_iterat)Then

          idx = Nlevel * ir + i
          jdx = Nlevel * ir + j
          H1(idx, jdx) = H1(idx, jdx) + CPOT0(i, j)
          
          
          If(i == j)then
            H1(idx, jdx) = H1(idx, jdx) + 2*t - E
          End If

        End If 

      End Do
    End Do

    If(ir .lt. R_iterat)Then
      Do i = 1, Nlevel
          idx = i + ir* Nlevel
          H1(idx, idx) = H1(idx, idx) + V(r, L_i) 

      End Do
    End If
  End Do

  ! contruct the H, knetic term 
  Do i = 1, Nlevel*(R_iterat-1)
        H1(i, i+nlevel) = -t
        H1(i+nlevel, i) = -t
  End Do

  Do i = 1, Nlevel
    Echannel(i) = E - V(R_min, L_i) - CPOT(i,i,0)
    
    If(Echannel(i) > 0)Then
     
      !k = acos(1. - Echannel(i)/2. / t)/dr
      k = sqrt(2. * ReduceMass / Hbar**2 * Echannel(i))

      H1(i, i) = H1(i, i) + (-t * exp(ai * k * dr))

    Else
      !k = acos(1. - abs(Echannel(i))/2. / t)/dr
      k = sqrt(2. * ReduceMass / Hbar**2 * abs(Echannel(i)))
      H1(i, i) = H1(i, i) + (-t * exp(k * dr))
    End If 

  End do

  Do i = 1, Nlevel*R_iterat
      Do j = 1, Nlevel*R_iterat
        H2(i,j) = H1(i,j)
      End Do
  End Do

  Do i = 1, Nlevel
    H2(Nlevel*R_iterat+i, Nlevel*R_iterat+i) = 2.0d0*t - E + eps(i) + V(R_max,L_i) 
    H2(Nlevel*R_iterat + i - nlevel, Nlevel*R_iterat + i) = -t
    H2(Nlevel*R_iterat + i, Nlevel*R_iterat + i  - nlevel) = -t
    H2(Nlevel*R_iterat + i, Nlevel*R_iterat + i  + nlevel) = -t
  End Do

  Do i = 1, Nlevel*R_iterat
      dbpsi(i, i) = (1.0d0, 0.0d0)
  End Do


  Do i = 1, Nlevel
    ec = E - CPOT(i,i,R_iterat)
    rho = acos(1. - ec/2. / t)/dr * (R_max)
    rho = sqrt(2.d0 * ReduceMass * ec) / Hbar * (R_max)

    eta = (Zpro * Ztar / 137.d0) * Sqrt(ReduceMass / (2.d0 * ec))
    Call myCOULFG(L_i*1.d0, eta, rho, fcw(L_i), fpcw(L_i), gcw(L_i), gpcw(L_i), iexp(L_i))
    dbpsi(Nlevel*R_iterat + i, Nlevel*R_iterat + i) = gcw(L_i) + ai*fcw(L_i)

    ec = E - CPOT(i,i,R_iterat+1)
    !rho = acos(1. - ec/2. / t)/dr  * (R_max + dr)
    rho = sqrt(2.d0 * ReduceMass * ec) / Hbar * (R_max + dr)
    eta = (Zpro * Ztar / 137.d0) * Sqrt(ReduceMass / (2.d0 * ec))
    Call myCOULFG(L_i*1.d0, eta, rho, fcw(L_i), fpcw(L_i), gcw(L_i), gpcw(L_i), iexp(L_i))
    dbpsi(Nlevel*R_iterat + i + Nlevel, Nlevel*R_iterat + i) = gcw(L_i) + ai*fcw(L_i)
  End Do

  do i = 1, Nx2
      dbpsi0(i) = dconjg(dbpsi(i, Nx1 - Nlevel + 1))
  end do

End Subroutine Prameterize_H

Subroutine DiscreteBasis(P)
  Use ccfull_initialization_mod

  Implicit None
  real(8), intent(out) :: P
  Real(8) ::  k, kk
  Integer ::  i, j
  complex*16 Hpsi(Nx1, Nx1)
  complex*16 Hpsi0(Nx1, 1)
  complex*16 cc(Nx1)
  Call MATMUL_CPLX(Hamilton, Ec_psi,  Hpsi, Nx1,  Nx2, Nx1)
  Call MATMUL_CPLX(Hamilton, Ec_psi0, Hpsi0, Nx1, Nx2, 1)
  call CSOLVE_EQUATION(Nx1,Hpsi,Hpsi0, cc)

  P = 0.

  Do i = 1, Nlevel
    if(Echannel(i) .gt. 0)then
      !k  = acos(1. - Echannel(i)/2. / t)/dr
      !kk = acos(1. - E/2. / t)/dr

      k=sqrt((2.d0*ReduceMass/Hbar**2*Echannel(i)))
      kk=sqrt(2.d0*ReduceMass/hbar**2*E)

      P = P+(abs(cc(i)))**2*k/kk
    end if 
  end do

  Do i = 1, Nx2
    wf(i) = Ec_psi0(i)
  End Do

  Do i = 1, Nx1
    Do j = 1, Nx2
      wf(j) = wf(j) + cc(i) * Ec_psi(j, i)
    End Do
  End Do
  !wf = wf / wf(R_iterat)
End Subroutine

SUBROUTINE MATMUL_CPLX(H2, PSI, HPSI, M, N, P)
    IMPLICIT NONE
    INTEGER M, N, P
    COMPLEX*16 H2(M, N), PSI(N, P), HPSI(M, P)

    COMPLEX*16 ALPHA, BETA
    CHARACTER*1 TRANSA, TRANSB

    ALPHA = (1.0D0, 0.0D0)
    BETA  = (0.0D0, 0.0D0)
    TRANSA = 'N'
    TRANSB = 'N'
    CALL ZGEMM(TRANSA, TRANSB, M, P, N, ALPHA,H2,M,PSI, N,BETA, HPSI, M)
    RETURN
END

SUBROUTINE CSOLVE_EQUATION(N, HPSI, HPSI0, CC)
      IMPLICIT NONE
      INTEGER N, INFO
      COMPLEX*16 HPSI(N,N), HPSI0(N), CC(N)

      INTEGER IPIV(N)
      COMPLEX*16 A(N,N), B(N)
      INTEGER I, J

      DO 20 J = 1, N
          DO 10 I = 1, N
              A(I,J) = HPSI(I,J)
   10     CONTINUE
   20 CONTINUE

      DO 30 I = 1, N
          B(I) = -HPSI0(I)
   30 CONTINUE

      CALL ZGESV(N, 1, A, N, IPIV, B, N, INFO)

      DO 40 I = 1, N
          CC(I) = B(I)
   40 CONTINUE

      RETURN
END

Subroutine myCOULFG(angL,eta,rho,myf,mydf,myg,mydg,iv)
  double precision  angL,eta, rho, myf, mydf,myg,mydg
  integer iv,lmax 
  double precision XLMIN, XLMAX
  double precision XX, ETA1
  integer  MODE1, KFN, IFAIL
  double precision , allocatable, dimension(:):: FC,GC,FCP,GCP
  intent(in) angL
  intent(in) eta
  intent(in) rho
  intent(out) myf
  intent(out) mydf
  intent(out) myg
  intent(out) mydg
  intent(out) iv     
  XLMIN=angL
  XLMAX=angL
  lmax=IDINT(XLMAX)+1
  allocate (FC(lmax),GC(lmax),FCP(lmax),GCP(lmax))
  MODE1=1
  KFN=0 
  XX=rho
  ETA1=eta
  call COULFG(XX,ETA1,XLMIN,XLMAX,FC,GC,FCP,GCP,MODE1,KFN,IFAIL)
  myf=FC(lmax)
  mydf=FCP(lmax)
  myg=GC(lmax)
  mydg=GCP(lmax) 
  iv=IFAIL
  return 
End Subroutine

SUBROUTINE COULFG(XX, ETA1, XLMIN, XLMAX, FC, GC, FCP, GCP, MODE1, KFN, IFAIL)
    IMPLICIT NONE
    ! Argument type declarations
    DOUBLE PRECISION, INTENT(IN) :: XX, ETA1, XLMIN, XLMAX
    DOUBLE PRECISION, DIMENSION(*), INTENT(OUT) :: FC, GC, FCP, GCP
    INTEGER, INTENT(IN) :: MODE1, KFN
    INTEGER, INTENT(OUT) :: IFAIL

    ! Local variable declarations
    DOUBLE PRECISION :: ACCUR, ACC, ACC4, ACCH, ETA, GJWKB, PACCQ, X, XLM, E2MM1, DELL, XLL, XI, FCL, PK, PX, F, D, C
    DOUBLE PRECISION :: PK1, EK, RK2, TK, DF, FPL, XL, RL, EL, SL, FCL1, TA, WI, P, Q, AR, AI, BR, BI, DR, DI, DP, DQ
    DOUBLE PRECISION :: A, B, GAM, W, ALPHA, BETA, FCM, GCL, GPL, FJWKB, GCL1
    INTEGER :: MODE, NFP, NPQ, IEXP, M1, LXTRA, L1, LP, MAXL, L, LP1
    LOGICAL :: ETANE0, XLTURN

    ! Common block
    COMMON /STEE/ PACCQ, NFP, NPQ, IEXP, M1

    ! Constants
    DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, TEN2 = 1.0D2, ABORT = 2.0D4
    DOUBLE PRECISION, PARAMETER :: HALF = 0.5D0, TM30 = 1.0D-30, BIG = 1.0D+100
    DOUBLE PRECISION, PARAMETER :: RT2DPI = 0.7978845608286535587989211986876373D0

    ! Initialize accuracy and other variables
    ACCUR = 1.0D-16
    MODE = 1
    IF (MODE1 == 2 .OR. MODE1 == 3) MODE = MODE1
    IFAIL = 0
    IEXP = 1
    NPQ = 0
    ETA = ETA1
    GJWKB = ZERO
    PACCQ = ONE
    IF (KFN /= 0) ETA = ZERO
    ETANE0 = (ETA /= ZERO)
    ACC = ACCUR * 10.0D0
    ACC4 = ACC * TEN2 * TEN2
    ACCH = SQRT(ACC)

    ! Check for small XX
    IF (XX <= ACCH) THEN
        IFAIL = -1
        RETURN
    END IF

    X = XX
    XLM = XLMIN
    IF (KFN == 2) XLM = XLM - HALF
    IF (XLM <= -ONE .OR. XLMAX < XLMIN) THEN
        IFAIL = -2
        RETURN
    END IF

    E2MM1 = ETA*ETA + XLM*XLM + XLM
    XLTURN = (X*(X - TWO*ETA) < XLM*XLM + XLM)
    DELL = XLMAX - XLMIN + ACC
    LXTRA = INT(DELL)
    XLL = XLM + DBLE(LXTRA)

    M1 = MAX(INT(XLMIN + ACC), 0) + 1
    L1 = M1 + LXTRA

    XI = ONE / X
    FCL = ONE
    PK = XLL + ONE
    PX = PK + ABORT
    F = ETA/PK + PK*XI
    IF (ABS(F) < TM30) F = TM30
    D = ZERO
    C = F

    DO
        PK1 = PK + ONE
        EK = ETA / PK
        RK2 = ONE + EK*EK
        TK = (PK + PK1)*(XI + EK/PK1)
        D = TK - RK2 * D
        C = TK - RK2 / C
        IF (ABS(C) < TM30) C = TM30
        IF (ABS(D) < TM30) D = TM30
        D = ONE / D
        DF = D * C
        F = F * DF
        IF (D < ZERO) FCL = -FCL
        PK = PK1
        IF (PK > PX) THEN
            IFAIL = -3
            RETURN
        END IF
        IF (ABS(DF-ONE) < ACC) EXIT
    END DO

    NFP = PK - XLL - 1
    IF (LXTRA /= 0) THEN
        FCL = FCL / BIG
        FPL = FCL * F
        IF (MODE == 1) FCP(L1) = FPL
        FC(L1) = FCL
        XL = XLL
        RL = ONE
        EL = ZERO
        DO LP = 1, LXTRA
            IF (ETANE0) EL = ETA / XL
            IF (ETANE0) RL = SQRT(ONE + EL*EL)
            SL = EL + XL*XI
            L = L1 - LP
            FCL1 = (FCL * SL + FPL) / RL
            FPL = FCL1 * SL - FCL * RL
            FCL = FCL1
            FC(L) = FCL
            IF (MODE == 1) FCP(L) = FPL
            IF (MODE /= 3 .AND. ETANE0) GC(L+1) = RL
            IF (ABS(FCL) > BIG) THEN
                DO LP1 = L, M1+LXTRA
                    IF (MODE == 1) FCP(LP1) = FCP(LP1) * 1.0D-20
                    FC(LP1) = FC(LP1) * 1.0D-20
                END DO
                FCL = FC(L)
                FPL = FPL * 1.0D-20
            END IF
            XL = XL - ONE
        END DO
        IF (FCL == ZERO) FCL = ACC
        F = FPL / FCL
    END IF

    IF (XLTURN) CALL JWKB(X, ETA, MAX(XLM, ZERO), FJWKB, GJWKB, IEXP)
    IF (IEXP > 1 .OR. GJWKB > ONE/(ACCH*TEN2)) THEN
        W = FJWKB
        GAM = GJWKB * W
        P = F
        Q = ONE
    ELSE
        XLTURN = .FALSE.
        TA = TWO * ABORT
        PK = ZERO
        WI = ETA + ETA
        P = ZERO
        Q = ONE - ETA*XI
        AR = -E2MM1
        AI = ETA
        BR = TWO*(X - ETA)
        BI = TWO
        DR = BR/(BR*BR + BI*BI)
        DI = -BI/(BR*BR + BI*BI)
        DP = -XI*(AR*DI + AI*DR)
        DQ = XI*(AR*DR - AI*DI)
        DO
            P = P + DP
            Q = Q + DQ
            PK = PK + TWO
            AR = AR + PK
            AI = AI + WI
            BI = BI + TWO
            D = AR*DR - AI*DI + BR
            DI = AI*DR + AR*DI + BI
            C = ONE/(D*D + DI*DI)
            DR = C * D
            DI = -C * DI
            A = BR*DR - BI*DI - ONE
            B = BI*DR + BR*DI
            C = DP*A - DQ*B
            DQ = DP*B + DQ*A
            DP = C
            IF (PK > TA) THEN
                IFAIL = -4
                RETURN
            END IF
            IF (ABS(DP) + ABS(DQ) < (ABS(P) + ABS(Q)) * ACC) EXIT
        END DO
        NPQ = INT(PK / TWO)
        PACCQ = HALF * ACC / MIN(ABS(Q), ONE)
        IF (ABS(P) > ABS(Q)) PACCQ = PACCQ * ABS(P)
        GAM = (F - P) / Q
        IF (Q <= ACC4 * ABS(P)) THEN
            IFAIL = -5
            RETURN
        END IF
        W = ONE / SQRT((F - P)*GAM + Q)
    END IF

    ALPHA = ZERO
    IF (KFN == 1) ALPHA = XI
    IF (KFN == 2) ALPHA = XI * HALF
    BETA = ONE
    IF (KFN == 1) BETA = XI
    IF (KFN == 2) BETA = SQRT(XI) * RT2DPI
    FCM = SIGN(W, FCL) * BETA
    FC(M1) = FCM
    IF (MODE /= 3) THEN
        IF (.NOT. XLTURN) THEN
            GCL = FCM * GAM
        ELSE
            GCL = GJWKB * BETA
        END IF
        IF (KFN /= 0) GCL = -GCL
        GC(M1) = GCL
        GPL = GCL*(P - Q/GAM) - ALPHA*GCL
        IF (MODE /= 2) THEN
            GCP(M1) = GPL
            FCP(M1) = FCM*(F - ALPHA)
        END IF
    END IF

    IF (LXTRA == 0) RETURN

    W = BETA * W / ABS(FCL)
    MAXL = L1 - 1
    DO L = M1, MAXL
        IF (MODE /= 3) THEN
            XL = XL + ONE
            IF (ETANE0) EL = ETA / XL
            IF (ETANE0) RL = GC(L+1)
            SL = EL + XL*XI
            GCL1 = ((SL - ALPHA)*GCL - GPL) / RL
            IF (ABS(GCL1) > BIG) THEN
                IFAIL = L1 - L
                RETURN
            END IF
            GPL = RL*GCL - (SL + ALPHA)*GCL1
            GCL = GCL1
            GC(L+1) = GCL1
            IF (MODE /= 2) THEN
                GCP(L+1) = GPL
                FCP(L+1) = W*(FCP(L+1) - ALPHA*FC(L+1))
            END IF
        END IF
        FC(L+1) = W * FC(L+1)
    END DO

END SUBROUTINE COULFG

subroutine JWKB(XX, ETA1, XL, FJWKB, GJWKB, IEXP)
    implicit none
    ! ===== 输入输出参数 =====
    real(8), intent(in)  :: XX, ETA1, XL
    real(8), intent(out) :: FJWKB, GJWKB
    integer, intent(out) :: IEXP

    ! ===== 内部变量 =====
    real(8) :: ZERO, HALF, ONE, SIX, TEN
    real(8) :: DZ, RL35
    real(8) :: aloge, X, ETA, GH2, XLL1, HLL, HL, SL, RL2, GH
    real(8) :: PHI, PHI10

    ! ===== 常量赋值 =====
    ZERO = 0.0d0
    HALF = 0.5d0
    ONE  = 1.0d0
    SIX  = 6.0d0
    TEN  = 10.0d0
    DZ   = 0.0d0
    RL35 = 35.0d0

    aloge = log(TEN)

    ! ===== 核心计算 =====
    X   = XX
    ETA = ETA1

    GH2  = X * (2.0d0*ETA - X)
    XLL1 = max(XL*XL + XL, DZ)

    if (GH2 + XLL1 <= ZERO) then
        FJWKB = 0.0d0
        GJWKB = 0.0d0
        IEXP  = 0
        return
    end if

    HLL = XLL1 + SIX / RL35
    HL  = sqrt(HLL)
    SL  = ETA/HL + HL/X
    RL2 = ONE + ETA*ETA/HLL
    GH  = sqrt(GH2 + HLL) / X

    PHI = X*GH - HALF * ( HL*log((GH+SL)**2 / RL2) - log(GH) )
    if (ETA /= ZERO) PHI = PHI - ETA*atan2(X*GH, X - ETA)

    PHI10 = -PHI * aloge
    IEXP  = int(PHI10)

    if (IEXP > 70) then
        GJWKB = TEN**(PHI10 - dble(IEXP))
    else
        GJWKB = exp(-PHI)
        IEXP  = 0
    end if

    FJWKB = HALF / (GH * GJWKB)

end subroutine JWKB

Subroutine WHIT(HETA, R, XK, E, LL, F, FD, IE)
    implicit none
    ! ====== 输入输出参数 ======
    integer, intent(in) :: LL
    integer, intent(inout) :: IE
    real(8), intent(in) :: HETA, R, XK, E
    real(8), intent(out) :: F(LL+1), FD(LL+1)

    ! ====== 内部变量 ======
    integer :: L, LP1, LM, LMP1, IS, H
    integer :: I, M, M1, M2, JS, IFEQL
    real(8) :: fpmax, fpminl
    real(8) :: EE, AK, ETA, RHO, RHOA, PJE, A, B, C, D
    real(8) :: T(12), S(7)

    fpmax = 1.0d290
    L = LL + 1

    EE = -1.0d0
    AK = XK
    ETA = HETA
    LP1 = L + 1
    RHO = AK * R

    S(:) = 0.0d0

    ! ====== 设置 LM ======
    if (L <= 50) then
        LM = 60
    else
        LM = L + 10
    end if
    LMP1 = LM + 1
    IS = 7

    PJE = 30.0d0 * RHO + 1.0d0
    H = max(int(PJE), 4)
    H = RHO / H

    RHOA = 10.0d0 * (ETA + 1.0d0)
    if (RHOA <= RHO) then
        IFEQL = 1
        RHOA = RHO
    else
        IFEQL = 0
    end if

    PJE = RHOA / H + 0.5d0
    RHOA = H * int(PJE)

    if (IFEQL == 1) then
        if (RHOA <= RHO + 1.5d0 * H) then
            RHOA = RHO + 2.0d0 * H
        end if
    end if

    if (EE < 0.0d0) then
        ! ----------- 级数展开部分 ----------
        C = 1.0d0 / RHOA
        A = 1.0d0
        B = 1.0d0 - C * ETA
        F(1) = A
        FD(1) = B
        do M = 1, 26
            D = 0.5d0 * (ETA + dble(M - 1)) * (ETA + dble(M)) * C / dble(M)
            A = -A * D
            B = -B * D - A * C
            F(1) = F(1) + A
            FD(1) = FD(1) + B
        end do

        A = -ETA * log(2.0d0 * RHOA) - RHOA
        fpminl = -log(fpmax)
        if (IE == 0 .and. A < fpminl) IE = int(fpminl - A)
        A = exp(A + dble(IE))

        F(1) = A * F(1)
        FD(1) = A * FD(1) * (-1.0d0 - 2.0d0 * ETA / RHOA)

        if (IFEQL == 1) then
            S(IS) = F(1)
            if (IS == 7) then
                IS = 6
                RHOA = RHOA + H
                ! 再做一次展开
                C = 1.0d0 / RHOA
                A = 1.0d0
                B = 1.0d0 - C * ETA
                F(1) = A
                FD(1) = B
                do M = 1, 26
                    D = 0.5d0 * (ETA + dble(M - 1)) * (ETA + dble(M)) * C / dble(M)
                    A = -A * D
                    B = -B * D - A * C
                    F(1) = F(1) + A
                    FD(1) = FD(1) + B
                end do
                A = -ETA * log(2.0d0 * RHOA) - RHOA
                fpminl = -log(fpmax)
                if (IE == 0 .and. A < fpminl) IE = int(fpminl - A)
                A = exp(A + dble(IE))
                F(1) = A * F(1)
                FD(1) = A * FD(1) * (-1.0d0 - 2.0d0 * ETA / RHOA)
            end if
        else
            F(1) = T(1)
            FD(1) = T(2)
        end if

    else
        stop "WHIT1: EE >= 0 not implemented"
    end if

    ! ----------- 递推关系 ----------
    C = 1.0d0 / RHO
    do M = 1, L-1
        A = ETA / dble(M)
        B = A + C * dble(M)
        F(M+1) = (B * F(M) - FD(M)) / (A + 1.0d0)
        FD(M+1) = (A - 1.0d0) * F(M) - B * F(M+1)
    end do

    do M = 1, L
        FD(M) = AK * FD(M)
    end do

    return
End Subroutine



Program CCFull_Main
  Use ccfull_initialization_mod
  Implicit None
  Integer ::  i, j, ir, ie, il
  Integer ::  Iestep
  Real(8) ::  sigma, spin, P_0, P, s0, P1,P2
  Real(8) ::  V, r, rh, k, kk
  Real(8) ::  test_set
  real(8), dimension(2) :: training_set

  external V
  Call Initialize_CCFull()
  Write(*,*)'Nlevle = ', Nlevel

  Iestep = int((Emax - Emin)/dE)
  Do ie = 0, Iestep

    E = Emin + dE * ie
    sigma = 0.d0
    spin  = 0.d0
    P_0 = 0.d0


    Do il = 0, 200
      L_i = il
      s0 = sigma

      Call PotShape(L_i, R_barrier_l, V_barrier_l, curv_l, R_bottom_l, V_bottom_l)
      If( V(R_bottom_l, L_i) > E .or. R_bottom_l < 0 )Then
          P = 0.
          exit  
      End If

      if (E == 55 .and. L_i == 0)Then
        Call Prameterize_H(Hamilton, Ec_psi, Ec_psi0)
        Call DiscreteBasis(P)
        Do i = 0, R_iterat+1
          Write(40,*)R_min+i*dr, Real(wf(nlevel*i+1)),  AIMAG(wf(nlevel*i+1)), Real(wf(nlevel*i+2)),  AIMAG(wf(nlevel*i+2))
        End Do 
      End If

      Call Prameterize_H(Hamilton, Ec_psi, Ec_psi0)
      Call DiscreteBasis(P)
      sigma = sigma + (2.0D0 * L_i + 1.0D0) * P * pi * Hbar**2 / (2.0D0 * ReduceMass * E) * 10.0D0
      spin  = spin  + (2.0D0 * L_i + 1.0D0) * P * pi * Hbar**2 / (2.0D0 * ReduceMass * E) * 10.0D0 * L_i
      If (L_i == 0) Then
          P_0 = P
      End If
      If(sigma - s0 < s0*1.d-4)Exit 
    End Do


    if (sigma < 0.01D0) then
        write(*,'(A,F8.3,A,E15.5,A,F8.5,A,E15.5)') ' E = ', E, ' MeV,   Sigma = ', sigma,&
                                       ' mb, <L> = ', spin / sigma, ' hbar, <P_0> = ', P_0
        write(20,'(F15.5,E15.5,F15.5,E15.5)') E,  sigma,  spin / sigma,  P_0
    else
        write(*,'(A,F8.3,A,F15.5,A,F8.5,A,E15.5)') ' E = ', E, ' MeV,   Sigma = ', sigma,&
                                       ' mb, <L> = ', spin / sigma, ' hbar, <P_0> = ', P_0
        write(20,'(4F15.5)') E,  sigma,  spin / sigma,  P_0
    end if

  End Do

  Print *, "Program completed successfully."

End Program CCFull_Main
