!>
!> @file moments.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>
!>
!
!   Compute moments of f(x), using its Spline representation
!
MODULE globals
  USE bsplines
  USE matrix
  IMPLICIT NONE
  DOUBLE PRECISION, PARAMETER   :: pi = 3.14159265359d0
  DOUBLE PRECISION, ALLOCATABLE :: xgrid(:), rhs(:), sol(:)
  DOUBLE PRECISION, ALLOCATABLE :: finteg(:,:), moms(:)
  TYPE(spline1d), SAVE :: splx
  TYPE(gbmat), SAVE    :: mat
CONTAINS
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE dismat(spl, mat)
!
!   Assembly FE matrix mat using spline spl
!
    TYPE(spline1d), INTENT(in) :: spl
    TYPE(gbmat), INTENT(inout) :: mat
    INTEGER :: nrank, nx, nidbas, ngauss
    INTEGER :: i, igauss, iterm, iw, jt, irow, jcol
    DOUBLE PRECISION, ALLOCATABLE :: fun(:,:)
    DOUBLE PRECISION, ALLOCATABLE :: xgauss(:), wgauss(:)
    DOUBLE PRECISION :: contrib
!
    INTEGER :: kterms                ! Number of terms in weak form
    INTEGER, ALLOCATABLE :: idert(:), iderw(:)  ! Derivative order
    DOUBLE PRECISION, ALLOCATABLE :: coefs(:)
!===========================================================================
!              1.0 Prologue
!
!   Properties of spline space
!
    CALL get_dim(spl, nrank, nx, nidbas)
    ALLOCATE(fun(0:nidbas,0:1)) ! Spline and its 1st derivative
!
!   Weak form
!
    kterms = mat%nterms
    ALLOCATE(idert(kterms), iderw(kterms), coefs(kterms))
!
!   Gauss quadature
!
    CALL get_gauss(spl, ngauss)
    ALLOCATE(xgauss(ngauss), wgauss(ngauss))
!===========================================================================
!              2.0 Assembly loop
!
    DO i=1,nx
       CALL get_gauss(spl, ngauss, i, xgauss, wgauss)
       DO igauss=1,ngauss
          CALL basfun(xgauss(igauss), spl, fun, i)
          CALL coefeq(xgauss(igauss), idert, iderw, coefs)
          DO iterm=1,kterms
             DO jt=0,nidbas
                DO iw=0,nidbas
                   contrib = fun(jt,idert(iterm)) * coefs(iterm) * &
                        &    fun(iw,iderw(iterm)) * wgauss(igauss)
                   irow=i+iw; jcol=i+jt
                   CALL updtmat(mat, irow, jcol, contrib)
                END DO
             END DO
          END DO
       END DO
    END DO
!===========================================================================
!              9.0  Epilogue
!
    DEALLOCATE(fun)
    DEALLOCATE(xgauss, wgauss)
    DEALLOCATE(iderw, idert, coefs)
!
  CONTAINS
    SUBROUTINE coefeq(x, idt, idw, c)
      DOUBLE PRECISION, INTENT(in) :: x
      INTEGER, INTENT(out) :: idt(:), idw(:)
      DOUBLE PRECISION, INTENT(out) :: c(SIZE(idt))
!
! Mass matrix
!
      c(1) = 1.d0
      idt(1) = 0
      idw(1) = 0
    END SUBROUTINE coefeq
  END SUBROUTINE dismat
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE disrhs(spl, rhs)
!
!   Assenbly the RHS using spline spl
!
    TYPE(spline1d), INTENT(in) :: spl
    DOUBLE PRECISION, INTENT(out) :: rhs(:)
    INTEGER :: nrank, nx, nidbas, ngauss
    INTEGER :: i, igauss
    DOUBLE PRECISION, ALLOCATABLE :: fun(:,:)
    DOUBLE PRECISION, ALLOCATABLE :: xgauss(:), wgauss(:)
    DOUBLE PRECISION:: contrib
!===========================================================================
!              1.0 Prologue
!
!   Properties of spline space
!
    CALL get_dim(spl, nrank, nx, nidbas)
!
    ALLOCATE(fun(nidbas+1,1)) ! needs only basis functions (no derivatives)
!
!   Gauss quadature
!
    CALL get_gauss(spl, ngauss)
    ALLOCATE(xgauss(ngauss), wgauss(ngauss))
!===========================================================================
!              2.0 Assembly loop
!
    rhs(1:nrank) = 0.0d0
!
    DO i=1,nx
       CALL get_gauss(spl, ngauss, i, xgauss, wgauss)
       DO igauss=1,ngauss
          CALL basfun(xgauss(igauss), spl, fun, i)
          contrib = wgauss(igauss)*rhseq(xgauss(igauss))
          rhs(i:i+nidbas) = rhs(i:i+nidbas) + contrib*fun(:,1)
       END DO
    END DO
!===========================================================================
!              9.0  Epilogue
!
    DEALLOCATE(fun)
    DEALLOCATE(xgauss, wgauss)
!
  CONTAINS
    DOUBLE PRECISION FUNCTION rhseq(v)
      DOUBLE PRECISION, INTENT(in) :: v
      rhseq = SQRT(1.0d0/(2.0d0*pi)) * EXP(-0.5d0*v**2)
    END FUNCTION rhseq
  END SUBROUTINE disrhs
END MODULE globals
!
!================================================================================
PROGRAM main
  USE bsplines
  USE globals
  IMPLICIT NONE
  INTEGER          :: nidbas, nx, nmoms, ngauss, rank, kl, ku
  INTEGER          :: i
  DOUBLE PRECISION :: a, b, dx
!
!   Input
!
  WRITE(*,'(a)') 'Enter, nidbas, a, b, nx, nmoms'
  READ(*,*) nidbas, a, b, nx, nmoms
!
!   Equidistant mesh
!
  ALLOCATE(xgrid(0:nx))
  dx = (b-a)/REAL(nx)
  DO i=0, nx
     xgrid(i) = a + i*dx
  END DO
  WRITE(*,'(a/(8(1pe12.4)))') 'XGRID', xgrid
!
!   Set up spline
!
  ngauss = nidbas+1
  CALL set_spline(nidbas, ngauss, xgrid, splx)
!
!    Mass matrix
!
  CALL get_dim(splx, rank)       ! Rank of the FE Mass matrix
  kl = nidbas
  ku = kl
  CALL init(kl, ku, rank, 1, mat)
  WRITE(*,'(a,3i6)') 'kl, ku, rank', kl, ku, rank
  CALL dismat(splx, mat)
!
!   RHS
!
  ALLOCATE(rhs(rank), sol(rank))
  CALL disrhs(splx, rhs)
!
!   Solve for Spline coefs
!
  CALL factor(mat)
  CALL bsolve(mat, rhs, sol)
  WRITE(*,'(a/(8(1pe12.4)))') 'SOL', sol
  WRITE(*,'(a,1pe20.12)') ' Integral of sol using FINTG', fintg(splx, sol)
!
!   Moments
!
  ALLOCATE(finteg(rank,0:nmoms), moms(0:nmoms))
  CALL calc_integ(splx, finteg)
  DO i=0,nmoms
     moms(i) = DOT_PRODUCT(sol(:), finteg(:,i))
  END DO
  WRITE(*,'(a,i3)') 'Moments of orders from 0 to', nmoms
  WRITE(*,'(8(1pe20.12))') moms
!
  DEALLOCATE(finteg, moms)
  DEALLOCATE(rhs, sol)
  DEALLOCATE(xgrid)
  CALL destroy(mat)
  CALL destroy_sp(splx)
END PROGRAM main