!>
!> @file tmassmat.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> SPClibs 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.
!>
!> SPClibs 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 General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
PROGRAM main
!
!   FT of mass matrix
!
  USE bsplines
  USE matrix
!
  IMPLICIT NONE
  INTEGER :: nx, nidbas
  INTEGER :: ngauss, nrank, kl, ku
  INTEGER :: i, k, kmin, kmax
  DOUBLE PRECISION :: pi, xlenght, dx, arg0, arg
  DOUBLE PRECISION, ALLOCATABLE :: xgrid(:), arow(:)
  DOUBLE PRECISION, ALLOCATABLE :: fftmassm1(:), fftmassm2(:), fftmassm3(:)
  DOUBLE PRECISION, ALLOCATABLE :: fftmass_shifted(:)
  TYPE(spline1d) :: splx
  TYPE(periodic_mat) :: massm
!
  NAMELIST /newrun/ nx, nidbas, xlenght
!================================================================================
!              1.0 Prologue
!
  pi = 4.0d0*ATAN(1.0d0)
!
  nx = 8
  nidbas = 3
  xlenght = 2.0d0*pi
!
  READ(*,newrun)
  WRITE(*,newrun)
!
  ngauss = nidbas+1  ! Exact integration for polynomials of degree 2*nidbas
!
  ALLOCATE(xgrid(0:nx))
  dx = xlenght/REAL(nx)
  xgrid = (/ (i*dx,i=0,nx) /)
  WRITE(*,'(a/(10f8.3))') 'xgrid', xgrid
!
  CALL set_spline(nidbas, ngauss, xgrid, splx, .TRUE.)
!===========================================================================
!              2.0  Mass matrix
  nrank = nx
  kl = nidbas
  ku = kl
  CALL init(kl, ku, nrank, 1, massm)
  CALL dismat(splx, massm)
!
  ALLOCATE(arow(nrank))
!!$  WRITE(*,'(/a)') 'Mass matrix'
!!$  DO i=1,nrank
!!$     CALL getrow(massm, i, arow)
!!$     WRITE(*,'(10(1pe12.4))') arow
!!$  END DO
!===========================================================================
!              3.0  Fourier transform of Mass matrix
!
  ALLOCATE(fftmassm1(0:nx-1))
  ALLOCATE(fftmassm2(0:nx-1))
  ALLOCATE(fftmassm3(0:nx-1))
  IF(nidbas.LE.3) THEN
     CALL analytic(nidbas, fftmassm1)
     fftmassm1 = dx*fftmassm1
     WRITE(*,'(/a/(10(1pe12.4)))') 'Analytic', fftmassm1
  END IF
!
  CALL calc_fftmass_old(splx, fftmassm2)
  WRITE(*,'(/a/(10(1pe12.4)))') 'Old version', fftmassm2
  IF(nidbas.LE.3) THEN
     WRITE(*,'(a,1pe12.4)') 'error =', MAXVAL(ABS(fftmassm2-fftmassm1))
  END IF
!
!   Init DFT
  kmin = -nx/2
  kmax = nx/2-1
  CALL init_dft(splx, kmin, kmax)
!
  ALLOCATE(fftmass_shifted(kmin:kmax))
  CALL calc_fftmass(splx, fftmass_shifted)
  DO k=kmin, kmax
     fftmassm3(MODULO(k+nx,nx))=fftmass_shifted(k)
  END DO
  WRITE(*,'(/a/(10(1pe12.4)))') 'New version', fftmassm3
  WRITE(*,'(a,1pe12.4)') 'error =', MAXVAL(ABS(fftmassm3-fftmassm2))
!
!   Check dftcoefs
  WRITE(*,'(/a)') 'Check DFT of splines'
  PRINT*, 'dims of dftcoefs', SHAPE(splx%dft%coefs)
  DO i=0,nidbas
     WRITE(*,'(a,i3,2(1pe12.4))') 'Sum of coefs for spline', i, &
          &                        SUM(splx%dft%coefs(:,i))
  END DO
!===========================================================================
!              9.0  Clean up
!
  DEALLOCATE(xgrid)
  DEALLOCATE(fftmassm1)
  DEALLOCATE(fftmassm2)
  DEALLOCATE(fftmassm3)
  DEALLOCATE(fftmass_shifted)
  DEALLOCATE(arow)
  CALL destroy(massm)
  CALL destroy_sp(splx)
CONTAINS
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE dismat(spl, mat)
    TYPE(spline1d) :: spl
    TYPE(periodic_mat) :: mat
    INTEGER :: dim, nx, nidbas, ngauss
    DOUBLE PRECISION, ALLOCATABLE :: fun(:,:), xgauss(:), wgauss(:)
    INTEGER :: i, igauss, iw, jt, irow, jcol
    DOUBLE PRECISION :: contrib
!
    CALL get_dim(spl, dim, nx, nidbas)
    ALLOCATE(fun(0:nidbas,1)) ! Spline
    CALL get_gauss(spl, ngauss)
    ALLOCATE(xgauss(ngauss), wgauss(ngauss))
!
    DO i=1,nx
       CALL get_gauss(spl, ngauss, i, xgauss, wgauss)
       DO igauss=1,ngauss
          CALL basfun(xgauss(igauss), spl, fun, i)
          DO jt=0,nidbas
             DO iw=0,nidbas
                contrib = fun(jt,1) * fun(iw,1) * wgauss(igauss)
                irow=MODULO(i+iw-1,nx) + 1   ! Periodic BC
                jcol=MODULO(i+jt-1,nx) + 1
                CALL updtmat(mat, irow, jcol, contrib)
             END DO
          END DO
       END DO
    END DO
!
    DEALLOCATE(fun)
    DEALLOCATE(xgauss, wgauss)
  END SUBROUTINE dismat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE analytic(nidbas, mat)
!
!   Analytic form for nidbas .le. 3
!
    INTEGER, INTENT(in) :: nidbas
    DOUBLE PRECISION, INTENT(out) :: mat(0:)
    DOUBLE PRECISION :: arg0, arg, cosk
    INTEGER :: n, k
!
    n = SIZE(mat)
    arg0 = 2.0d0*pi/REAL(n,8)
    DO k=0,n-1
       arg = k*arg0
       cosk = COS(arg)
       SELECT CASE (nidbas)
       CASE (1)
          mat(k) = (2.0d0 + cosk)/3.0d0
       CASE (2)
          mat(k) = (16.0d0 + cosk*(13.0d0+cosk))/30.0d0
       CASE (3)
          mat(k) = (272.0d0 + cosk*(297.0d0+cosk*(60.0d0+cosk)))/630.0d0
       CASE default
          WRITE(*,'(a,i4,a)') 'ANALYTIC: nidbas =', nidbas, ' is not implemented!'
          STOP
       END SELECT
    END DO
  END SUBROUTINE analytic
END PROGRAM main