program findmax
! J. Loreau
! Goal: compute cross sections based on statistical method.
! proba = 0 or 1 if channel open or closed.
! Takes into account the centrifugal barrier in the adiabats.
! Developed for adiabats of COb- + COa PES
! version called from script, generates cross sections for all j_i
!***** To check: np_f, file rot_levels.dat, nlevtot, mu (red mass), max number
!                adiabats

  implicit none

  integer, parameter :: np_f=142  !number of points for function, input
!      integer, parameter :: nps=10000  !number of points we work with in the splint subroutine.
  integer, parameter :: nlevtot=165  ! number of levels, some appear multiple times
  integer, parameter :: npar=2  ! number of parities from molscat
  integer, parameter :: Jmin=0  ! Jtot min
  integer, parameter :: Jmax=5  ! Jtot max
!      integer, parameter :: inist=1  ! initial state, # refers to levels.dat
  integer, parameter :: nE=500        ! number of energies
!      integer, parameter :: nf=10  !number of functions in input
!      double precision, parameter :: c = 137.035999     !in au,
!      c=1/alpha
!      double precision, parameter :: c = 299792458  !in m/s
!      double precision, parameter :: amu=1.66d-27  !in kg
!      double precision, parameter :: kb=1.38d-23 !in J/K
!      double precision, parameter :: kb2=8.617d-05 !in eV/K
  integer :: i,j,k,m,idx,nf,nlevuniq,Jtot,idx_p1,idx_p2
  integer :: l1,l2,l1min,l1max,l2min,l2max,iE,inist,openi(Jmax+1)
  character(2) :: x1,x2,x3
  double precision :: de,pi,bid,mu,wavek2
  double precision :: R(np_f),Deltaf(np_f)
!      double precision :: func(np_f,nf)
  double precision :: Rm(5),fm(5) ! 5=max number of maxima here.
  double precision :: Rmax,fmax,ener(nE),Emin,Emax
  double precision, allocatable :: func(:,:),proba(:)
  double precision :: probalev(nlevtot,npar),probalevtot(nlevtot)
  double precision :: probaJ(Jmax+1,nlevtot)  ! Jmax+1 as index J=0 wont work
  double precision :: E_lev(nlevtot),sigidx(nlevtot),erel
  integer :: jCO(nlevtot),jH2O(nlevtot),ja(nlevtot),jb(nlevtot)
  integer :: tau(nlevtot),tau2(nlevtot),taui,taumax,tauf
  integer :: jamax,jbmax,jai,jaf,jbi,jbf,deltakr
  double precision :: Elev(nlevtot),numlev(nlevtot) ! this is too large, but to be safe
  double precision :: openadiab(npar),openadiabtot(Jmax+1)
  double precision :: crossJ(Jmax+1,nlevtot),crosstot(nlevtot)
  double precision :: Etot, Einf(Jmax+1,npar,15000)   ! 500 = max number of adiabats in file
  double precision :: Ebarrier(Jmax+1,npar,15000)   ! 500 = max number of adiabats in file
  integer :: n_adiab(Jmax+1,npar)
  character*20 name_CS
  character(len=8) :: fmt ! format descriptor
  pi=dacos(-1d0)

  fmt = '(I2.2)' ! an integer of width 2 with zeros at the left

!****** read functions, extract energy of centrifugal barrier ****
  Ebarrier=0d0
  do Jtot=Jmin,Jmax       ! begin loop over Jtot
    write(400,*) "Jtot = ", Jtot
    do m=1,npar  ! loop over parities
      idx_p1=Jtot+m*1000
      open(idx_p1)
      read(idx_p1,*) nf
      n_adiab(Jtot+1,m)=nf    ! store the number of functions for given J, parity
      if (nf.eq.0) then  ! do nothing, only for J=0 in some cases
      else
        allocate(func(np_f,nf))
!       allocate(proba(nf))

        do i=np_f,1,-1  ! read in reverse to start from Rmax
          read(idx_p1,*) R(i),(func(i,j),j=1,nf)
        enddo
!       rewind(idx_p1)
        close(idx_p1)
!       do i=1,nlevuniq
!         print*, numlev(i),Elev(i)
!       enddo

        do j=1,nf      ! loop over functions
          Rmax=0d0
          fmax=0d0
          Deltaf=0d0
    !      idx=0
          do i=3,np_f  ! loop over distances, start from Rmax
            Deltaf(i)=func(i,j)-func(i-1,j)
    !        print*, i,Deltaf(i)
            if ((Deltaf(i) .lt. 0d0) .and. (Deltaf(i-1)) .gt. 0d0) then !find maxima
    !          idx=idx+1
    !          Rmax=(R(i-1)+R(i))/2
    !          fmax=(func(i-1,1)+func(i,1))/2
    !           print*, j,R(i-1),func(i-1,j)
               if (func(i-1,j) .gt. fmax) then
                 fmax=func(i-1,j)
                 Rmax=R(i-1)
               endif
            endif
          enddo     ! end loop over distances
          Ebarrier(Jtot+1,m,j)=fmax            ! store max value
    !      print*, Ebarrier(Jtot+1,m,j)
          Einf(Jtot+1,m,j)=func(1,j)
    !      print*, Rmax,fmax
          if ((Rmax.eq.0d0).and.(Deltaf(3).lt.0d0)) then
            Ebarrier(Jtot+1,m,j)=func(1,j)       ! attractive -> barrier=E_asy
            write(400,*) "a", j,  "attractive"
          elseif ((Rmax.eq.0d0).and.(Deltaf(3).gt.0d0)) then
            Ebarrier(Jtot+1,m,j)=1d10         ! arbitrary large value, proba=0
            write(400,*) j,  "repulsive"
          elseif (fmax.lt.func(1,j)) then   
            Ebarrier(Jtot+1,m,j)=func(1,j)    ! takes care of local min in the well
            write(400,*) "b", j,  "attractive"
          endif
    ! note: the order in the "if" is important! Do not invert 2 and 3.

          write(400,*) "Ebarrier = ", Ebarrier(Jtot+1,m,j)

        enddo     ! end loop over functions
        deallocate(func)

      endif  ! end if over nf=0 (to treat the J=0 case)

    enddo   ! end loop over parities m

  enddo     !end loop over Jtot

  ! print*, Ebarrier(1,1,:)

  !******
  !* we now know Ebarrier and Einf for each J, parity, # adiabat
  !* we can start the loop over the energies.

  ! read(5,*) jai, taui, jbi         ! read initial state. ja=H2O, jb=CO

  open(9,file='level.dat',status='old')
  read(9,*)         ! determine the number of unique levels, nlevuniq
  nlevuniq=0
  do i=1,nlevtot
    write (6,*) 'i = ', i
    read(9,*) sigidx(i),jH2O(i),tau2(i),bid, bid,jCO(i),E_lev(i)  
    write (6,*) 'sigidx = ', sigidx(i)          
    if (sigidx(i).ne.sigidx(i-1)) then
      nlevuniq=nlevuniq+1
      numlev(nlevuniq)=sigidx(i)
      Elev(nlevuniq)=E_lev(i)
      ja(nlevuniq)=jH2O(i)
      jb(nlevuniq)=jCO(i)
      tau(nlevuniq)=tau2(i)
!          print*,nlevuniq,tau(nlevuniq)
    endif
  enddo
  close(9)
  jamax=maxval(ja)
  jbmax=maxval(jb)

  do jai=0,5
    taumax=2*jai+1
    ! print*, "taumax=",taumax
    do taui=1,taumax
      do jbi = 0,10
        ! print*, jai,taui,jbi

        inist=0
        write(6,*) 'nlevuniq = ', nlevuniq
        do i=1,nlevuniq
          write(6,*) 'i = ', i
          if ((ja(i).eq.jai) .and. (jb(i).eq.jbi) .and. (tau(i).eq.taui)) then
            inist=i    ! find number of level corresponding to initial state
          endif
        enddo
        print*, "inist=", inist, "with energy", Elev(inist)
        if (inist .ne. 0) then  ! if inist =0, dont do anything (wrong parity)
          ! print*, "initial state does not exist"
          ! print*, "inist=", inist

          write (x1,fmt) ja(inist) ! converting integer to string using a 'internal file'
          write (x2,fmt) tau(inist) !
          write (x3,fmt) jb(inist) !
          ! print*,x1,x2,x3
          name_CS='cross_'//x1//'_'//x2//'__'//x3//'.dat'
          ! print*, name_CS
          open(26,file=name_CS,status='replace')
          open(400,file='output_adiab',status='replace')

          Emin=1d0
          Emax=8d2

          do iE=1,nE    ! loop over energies
            ener(iE)=Emin+(iE-1)*(Emax-Emin)/(nE-1)  !define scale for energy
            ! ener(1)=500d0
            write(400,*) 'energy: ',ener(iE)
            ! ener(iE)=10+(iE-1)*20d0  !define scale for energy
            ! ener(iE)=ener(iE)-Elev(inist)     ! energy in cm-1

            probaJ=0d0
            crosstot=0d0

            do Jtot=Jmin,Jmax       ! begin loop over Jtot
              ! write(6,*) "E = ", ener(iE), "Jtot = ", Jtot

              openi(Jtot+1)=0
              probalev=0d0
              probalevtot=0d0
              openadiab=0d0

              do m=1,npar  ! loop over parities

                nf=n_adiab(Jtot+1,m)

                if (nf.eq.0) then  ! do nothing, only for J=0 for some cases
                else

                  allocate(proba(nf))

                  do j=1,nf         ! loop over adiabats

                    do i=1,nlevuniq ! determine the level associated with the adiabat
                      ! if (Einf(Jtot+1,m,j).eq.Elev(i)) then
                      if (abs(Einf(Jtot+1,m,j)-Elev(i)).lt.1d-4) then
                        k=i
                      endif
                    enddo
                    Etot=ener(iE)+Elev(inist)

                    if (Etot.ge.Ebarrier(Jtot+1,m,j)) then
                      proba(j)=1d0
                      ! if (Einf(Jtot+1,m,j).eq.Elev(inist)) then
                      if (abs(Einf(Jtot+1,m,j)-Elev(inist)).lt.1d-3) then
                        openi(Jtot+1)=openi(Jtot+1)+1   ! count the number of open initial channels 
                      endif
                    else
                      proba(j)=0d0
                    endif

                    probalev(k,m)=probalev(k,m)+proba(j)  ! sum the probas for each unique level and parity
                    openadiab(m)=openadiab(m)+proba(j)   ! counts the number of open channels per parity
                  enddo     ! end loop over adiabats
                  write(400,*) "Jtot = ", Jtot, "open channels for Jtot=", Jtot, "and parity: ", m, "= ", openadiab(m)

                  ! print*, "J = ", Jtot
                  ! do i=1,nlevuniq
                  !   print*, i,probalev(i,m),probalev(i,m)/sum(proba)
                  ! enddo

                  deallocate(proba)
                  ! Ptot=sum(proba)
                  ! print*, Ptot
                  ! do i=1,nlevuniq ! determine the level associated with the adiabat
                  !   if (func(1,j).eq.Elev(i)) then
                  !     k=i
                  !   endif
                  ! enddo
                  !
                  ! do j=1,nf
                  !   if k(j)=i then
                  !     probalev=probalev+proba(j)
                  !   endif

                endif  ! end if over nf=0 (to treat the J=0 case)

              enddo   ! end loop over parities m
              openadiabtot(Jtot+1)=sum(openadiab)       ! sum over m
              write(400,*) "channel n°,       open channels,      probability:"
              do i=1,nlevuniq
                do m=1,npar
                  probalevtot(i)=probalevtot(i)+probalev(i,m)
                enddo
                ! print*, i,probalev(i,m),probalev(i,m)/sum(proba)

                if (openadiabtot(Jtot+1).eq.0d0) then
                  probaJ(Jtot+1,i)=0d0    ! to avoid infinity below
                  write(400,*) i,probalevtot(i),probaJ(Jtot+1,i)
                else
                  probaJ(Jtot+1,i)=probalevtot(i)/sum(openadiab)
                  write(400,*) i,probalevtot(i),probaJ(Jtot+1,i)
                endif
              enddo

            enddo     !end loop over Jtot


            ! call cross(probaJ,)

            ! *** Compute the cross section based on S-matrix (= probas)

            mu=10.956d0 ! reduced mass
            wavek2=2d0*mu*1822.8885d0*ener(iE)/219474.63d0  ! in a.u.
            crossJ=0d0
            ! jbi=jb(inist)  ! already defined
            ! jai=0  ! COa
            do Jtot=Jmin,Jmax
              ! print*, Jtot, openi(Jtot+1)
              ! print*, openadiabtot(Jtot+1)
              if (openadiabtot(Jtot+1).eq.0d0) then
                crossJ(Jtot+1,:)=0d0
              else
                do i=1,nlevuniq
                  jaf=ja(i)
                  jbf=jb(i)
                  tauf=tau(i)

                  ! this sum is not necessary, already taken into account by adiabats #
                  ! l1min=abs(Jtot-ji)       ! change  ji and jf here!
                  ! l1max=Jtot+ji
                  ! l2min=abs(Jtot-jf)
                  ! l2max=Jtot+jf
                  ! print*, 'l1...', l1min, l1max, l2min, l2max
                  ! do l1=l1min,l1max
                  !   do l2=l2min,l2max
                        crossJ(Jtot+1,i)=crossJ(Jtot+1,i) &
                          ! (deltakr(ji,jf)*deltakr(l1,l2)-probaJ(Jtot+1,i))**2  &
                          ! (deltakr(ji,jf)-probaJ(Jtot+1,i))**2  !no square!  &
                          + abs(deltakr(jbi,jbf)*deltakr(jai,jaf)*deltakr(taui,tauf) &
                          - probaJ(Jtot+1,i))
                  !   enddo
                  ! enddo

                  crossJ(Jtot+1,i)=openi(Jtot+1)*crossJ(Jtot+1,i)*(2d0*Jtot+1d0)
                  crosstot(i)=crosstot(i)+crossJ(Jtot+1,i)
                  if (crossJ(Jtot+1,i).ne.0d0) then
                    write(300+i,*) ener(iE), Jtot, crossJ(Jtot+1,i)
                  endif
                enddo
              endif
            enddo

            write(26,*)"    energy       jbi   jai   taui    jbf   jaf  tauf      CS"
            do i=1,nlevuniq
              crosstot(i)=crosstot(i)*pi/((2d0*jai+1d0)*(2d0*jbi+1d0)*wavek2)
              crosstot(i)=crosstot(i)*((5.291772d-1)**2)
              if (crosstot(i).ne.0d0) then
                write(26,111) ener(iE),jbi,jai,taui,jb(i),ja(i),tau(i),crosstot(i)
              endif
            enddo

          enddo   ! end loop over energies
          close(400)
          close(26)

          ! print*, probaJ
111       format(1(ES14.6),6I6,ES14.4)
112       format(17(E16.8,1h,))

        endif ! end if over inist=0
      enddo  ! end loops over jbi
    enddo  ! end loops over taui
  enddo  ! end loop over jai
end program 


function deltakr(a,b)
  implicit none
  integer :: a,b,deltakr
  if (a.eq.b) then
    deltakr=1d0
  else
    deltakr=0d0
  endif
end function

!***********************************************************************