!
!    Copyright 2020, 2021 Guy Munhoven
!
!    This file is part of SolveSAPHE v. 2

!    SolveSAPHE is free software: you can redistribute it and/or modify
!    it under the terms of the GNU Lesser General Public License as published by
!    the Free Software Foundation, either version 3 of the License, or
!    (at your option) any later version.
!
!    SolveSAPHE is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU Lesser General Public License for more details.
!
!    You should have received a copy of the GNU Lesser General Public License
!    along with SolveSAPHE.  If not, see <http://www.gnu.org/licenses/>.
!

program sdt


implicit none

integer, parameter :: ilevel_verbose = 0
                                    ! Set the following to true
logical, parameter :: lflag_printzeroes = .true.
logical, parameter :: lflag_printtexcode = .false.

integer :: i, j, k
integer :: n, nn

integer, allocatable, dimension(:)   :: iarr_m, iarr_f
integer, allocatable, dimension(:,:) :: iarr_prodpi
integer,              dimension(3)   :: iswap3
integer                              :: iswap

integer :: ijk, idx, idx_sorted, idx1, idx2
integer :: m, m_max, m_min, m_prev
integer, allocatable, dimension(:) :: n_m


integer, allocatable, dimension(:)   :: ijk_sorted

write(*,'(A)', advance="no") 'Please enter number of protons in the acid: '
read*, n

nn = (n+1) * ((n+1)*n) / 2 ! k : n+1
                           ! i, j: n + (n-1) + ... + 1
m_max = 3*n - 1  ! for (i,j,k) = (0,1,0)
m_min = 1        ! for (i,j,k) = (n-1, n, n)



allocate(iarr_m(nn))
allocate(iarr_f(nn))
allocate(iarr_prodpi(3,nn))

allocate(ijk_sorted(nn))

allocate(n_m(m_min:m_max))

ijk = 1
do k = 0, n
  do i = 0, n-1
    do j = i + 1, n
      iarr_m(ijk) = 3*n - i - j - k   ! exponent
      iarr_f(ijk) = (i - j)**2 * (i + j + 1 - 2*k)
      iarr_prodpi(1,ijk) = i
      iarr_prodpi(2,ijk) = j
      iarr_prodpi(3,ijk) = k
      ijk_sorted(ijk) = ijk
      ijk = ijk + 1
    enddo
  enddo
enddo


if (ilevel_verbose > 1) THEN
  do ijk = 1, nn
    write(*,'(3I3, " : ", i5, " * H^", I0 )') iarr_prodpi(:,ijk), iarr_f(ijk), iarr_m(ijk)
  enddo
endif


! order iarr_prodpi(:, ijk) for each ijk

do ijk = 1, nn

  if (iarr_prodpi(1,ijk) >  iarr_prodpi(2,ijk)) then
    iswap              = iarr_prodpi(1,ijk)
    iarr_prodpi(1,ijk) = iarr_prodpi(2,ijk)
    iarr_prodpi(2,ijk) = iswap
  endif

  if (iarr_prodpi(2,ijk) >  iarr_prodpi(3,ijk)) then
    iswap              = iarr_prodpi(2,ijk)
    iarr_prodpi(2,ijk) = iarr_prodpi(3,ijk)
    iarr_prodpi(3,ijk) = iswap
  endif

  if (iarr_prodpi(1,ijk) >  iarr_prodpi(2,ijk)) then
    iswap              = iarr_prodpi(1,ijk)
    iarr_prodpi(1,ijk) = iarr_prodpi(2,ijk)
    iarr_prodpi(2,ijk) = iswap
  endif

enddo



ijk = 1

n_m(:) = 0

do m = m_max, m_min, -1

  do idx_sorted = ijk, nn

    idx = ijk_sorted(idx_sorted)

    if (iarr_m(idx) == m) then

      if (idx == ijk) then

        n_m(m) = n_m(m) + 1
        ijk = ijk +1
        cycle

      else

        ijk_sorted(idx_sorted) = ijk_sorted(ijk)
        ijk_sorted(ijk)        = idx

        n_m(m) = n_m(m) + 1
        ijk = ijk +1

      endif

    endif      

  enddo

enddo

if (ilevel_verbose > 1) THEN
  write(*,'()')
  write(*,'()')
  
  m_prev = m_max
  do ijk = 1, nn
    idx = ijk_sorted(ijk)
    if (iarr_m(idx) /= m_prev) then
      write(*,'()')
      m_prev = iarr_m(idx)
    endif
    if ((iarr_f(idx) /= 0) .or. lflag_printzeroes) &
    write(*,'(3I3, " : ", i5, " * H^", I0 )') iarr_prodpi(:,idx), iarr_f(idx), iarr_m(idx)
  enddo
endif

! Within each exponent:
!   sort iarr_prodpi(3, :) ascending
!     then sort iarr_prodpi(2, :) ascending within each of the  iarr_prodpi(3, :)
!       then sort iarr_prodpi(1, :) ascending within each of the  iarr_prodpi(2, :)

ijk = 0

do m = m_max, m_min, -1

  do i = 1, n_m(m) - 1
    do j = ijk + 1, ijk + n_m(m) - i
      idx1 = ijk_sorted(j)
      idx2 = ijk_sorted(j+1)
      if     (iarr_prodpi(3,idx2)  < iarr_prodpi(3,idx1)) then
        ijk_sorted(j+1) = idx1
        ijk_sorted(j)   = idx2
      elseif (iarr_prodpi(3,idx2) == iarr_prodpi(3,idx1)) then
        if     (iarr_prodpi(2,idx2)  < iarr_prodpi(2,idx1)) then
          ijk_sorted(j+1) = idx1
          ijk_sorted(j)   = idx2
        elseif (iarr_prodpi(2,idx2) == iarr_prodpi(2,idx1)) then
          if      (iarr_prodpi(2,idx2) < iarr_prodpi(2,idx1)) then
            ijk_sorted(j+1) = idx1
            ijk_sorted(j)   = idx2
          endif
        endif
      endif
    enddo 
  enddo
  ijk = ijk + n_m(m)

enddo


if (ilevel_verbose > 0) THEN
  write(*,'()')
  write(*,'()')
  
  m_prev = m_max
  do ijk = 1, nn
    idx = ijk_sorted(ijk)
    if (iarr_m(idx) /= m_prev) then
      write(*,'()')
      m_prev = iarr_m(idx)
    endif
    if ((iarr_f(idx) /= 0) .or. lflag_printzeroes) &
    write(*,'(3I3, " : ", i5, " * H^", I0 )') iarr_prodpi(:,idx), iarr_f(idx), iarr_m(idx)
  enddo
endif


! Within each exponent: simplify the triplets

ijk = 0

do m = m_max, m_min, -1

  do i = ijk + 1, ijk + n_m(m) - 1
    idx1 = ijk_sorted(i)
    idx2 = ijk_sorted(i+1)
    if (ALL(iarr_prodpi(:,idx1) == iarr_prodpi(:,idx2))) then
      iarr_f(idx2) = iarr_f(idx2) + iarr_f(idx1)
      iarr_f(idx1) = 0   ! Set iarr_f(idx1) to zero
      iarr_m(idx1) = -1  ! and flag the term idx1 obsolete
    endif
  enddo
  ijk = ijk + n_m(m)

enddo


write(*,'()')
write(*,'()')
if (lflag_printtexcode) then
  write(*,'("$n = ", i0, "$")') n
else
  write(*,'("n = ", i0)') n
endif
write(*,'()')

m_prev = m_max
do ijk = 1, nn
  idx = ijk_sorted(ijk)
  if (iarr_m(idx) == -1) cycle      ! Skip the obsolete terms
  if (iarr_m(idx) /= m_prev) then
    write(*,'()')
    m_prev = iarr_m(idx)
  endif
  
  if ((iarr_f(idx) /= 0) .or. lflag_printzeroes) then
    if (lflag_printtexcode) then
      write(*,'("$", i0, 1x, 3("\Pi_{",i0,"} "), "H^{", I0, "}$")') iarr_f(idx), iarr_prodpi(:,idx), iarr_m(idx)
    else
      write(*,'(3I3, " : ", i5, " * H^", I0 )') iarr_prodpi(:,idx), iarr_f(idx), iarr_m(idx)
    endif
  endif
enddo

!~ write(*,'()')
!~ do m = m_max, m_min, -1
!~   write(*,'(i3, 2x, i3)') m, n_m(m)
!~ enddo

end
