!>
!> @file gyro.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
!
!   Gyro-average using splines.
!   F(x)=cos(x) => \bar F(x,rho) = J_0(rho) cos(x)
!
  USE bsplines
  USE futils
!
  IMPLICIT NONE
  INTEGER, PARAMETER :: nx=10, nidbas=3, dim=nx+nidbas, npt=100, &
       &                nrho=21, nnq=5
  DOUBLE PRECISION :: xgrid(0:nx), fgrid(0:nx), coefs(dim)
  DOUBLE PRECISION :: xpt(npt), fcalc(npt), fexact(npt)
  DOUBLE PRECISION :: averf(0:nx), averfexact(0:nx)
  DOUBLE PRECISION, POINTER :: splines(:,:)
  TYPE(spline1d) :: spl
  DOUBLE PRECISION :: pi, twopi, dx, lperiod, dth
  DOUBLE PRECISION :: drho, rho(nrho), erraver(nrho,nnq)
  INTEGER :: i, j, nq(nnq)
  DOUBLE PRECISION :: dbesj0
!
  CHARACTER(len=128) :: file='gyro.h5'
  INTEGER :: fid
!
  INTERFACE
     SUBROUTINE gyro(spl, xgrid, coefs, rho, nq, averf)
       USE bsplines
       TYPE(spline1d) :: spl
       DOUBLE PRECISION, INTENT(in) :: xgrid(0:), coefs(:), rho
       INTEGER, INTENT(in) :: nq
       DOUBLE PRECISION, INTENT(out) :: averf(0:)
     END SUBROUTINE gyro
  END INTERFACE
!
  pi = 4.0d0*ATAN(1.0d0)
  twopi = 2.0d0*pi
!
!   Create hdf5 file
  CALL creatf(file, fid, 'gyro Result File')
  CALL attach(fid, '/', 'NX', nx)
  CALL attach(fid, '/', 'NIDBAS', nidbas)
!
!   Grid and function values
  dx = twopi/nx
  xgrid(0) = 0.0d0
  DO i=1,nx
     xgrid(i) = xgrid(0) + i*dx
  END DO
  lperiod = xgrid(nx)-xgrid(0)
  fgrid = func(xgrid)
!
!   Spline interpolation
  CALL set_splcoef(nidbas, xgrid, spl, period=.TRUE.)
  CALL get_splcoef(spl, fgrid, coefs)
  WRITE(*,'(a)') 'Spline coefficients'
  WRITE(*,'(i5,f12.4)') (i-1, coefs(i), i=1,dim)
!
!   Error of interpolation
  CALL RANDOM_NUMBER(xpt)
  xpt = twopi*xpt
  CALL gridval(spl, xpt, fcalc, 0, coefs)
  fexact = func(xpt)
!!$  WRITE(*,'(a)') 'Interpolated and exact f'
!!$  WRITE(*,'(3(1pe12.3))') (xpt(i), fcalc(i), fexact(i), i=1,npt)
  WRITE(*,'(a,1pe12.3)') 'Interpolation error', norm2(fcalc-fexact)
  CALL putarr(fid, '/X', xpt)
  CALL putarr(fid, '/FEXACT', fexact, 'Exact values')
  CALL putarr(fid, '/FCALC', fcalc, 'Interpolated values')
!
!   Gyro-averaged of F at grid points xgrid(0:nx-1)
  drho = 5.0d0/nrho
  DO i=1,nrho
     rho(i) = i*drho
     DO j=1,nnq
        nq(j) = 2**(j+1)
        CALL gyro(spl, xgrid, coefs, rho(i), nq(j), averf)
        averfexact = dbesj0(rho(i))*COS(xgrid)
        erraver(i,j) = norm2(averfexact-averf)
     END DO
  END DO
!
  WRITE(*,'(a,f8.3,i5)') 'averaged f at rho, nq =', rho(nrho), nq(nnq)
  WRITE(*,'(3(1pe12.3))') (xgrid(i),averf(i),averfexact(i), i=0,nx)
  CALL putarr(fid, '/XGRID', xgrid)
  CALL putarr(fid, '/AVERF', averf, 'Averaged F')
  CALL putarr(fid, '/AVERFEXACT', averfexact, 'Averaged F exact')
!
  CALL putarr(fid, '/RHO', rho)
  CALL putarr(fid, '/NQ', nq)
  CALL putarr(fid, '/ERRAVER', erraver)
!
!   Clean up
  CALL destroy_sp(spl)
  CALL closef(fid)
CONTAINS
  FUNCTION func(x)
    DOUBLE PRECISION, INTENT(in) :: x(:)
    DOUBLE PRECISION :: func(SIZE(x))
    func = COS(x)
  END FUNCTION func
  FUNCTION norm2(x)
    DOUBLE PRECISION :: norm2
    DOUBLE PRECISION, INTENT(in) :: x(:)
    norm2 = SQRT(DOT_PRODUCT(x,x))
  END FUNCTION norm2
END PROGRAM main
!
SUBROUTINE gyro(spl, xgrid, coefs, rho, nq, averf)
!
!  Gyro-average using spline SPL and NQ point quadratire for
!  theta-integration.
!
  USE bsplines
  IMPLICIT NONE
  TYPE(spline1d) :: spl
  DOUBLE PRECISION, INTENT(in) :: xgrid(0:), coefs(:), rho
  INTEGER, INTENT(in) :: nq
  DOUBLE PRECISION, INTENT(out) :: averf(0:)
!
  DOUBLE PRECISION :: th(nq), wth(nq), xq(nq)
  DOUBLE PRECISION, ALLOCATABLE :: avermat(:,:)
  DOUBLE PRECISION :: pi, twopi, lperiod, dth
  DOUBLE PRECISION, POINTER :: splines(:,:)
  INTEGER :: i, j, iq, nx, dim
!
  pi = 4.0d0*ATAN(1.0d0)
  twopi = 2.0d0*pi
!
!   Quadrature in theta
  dth = twopi/nq
  th(1) = -pi + dth/2.0d0
  DO iq=2,nq
     th(iq) = th(iq-1)+dth
  END DO
  wth = dth
!
!   Gyro-averaging matrix
  CALL get_dim(spl, dim, nx)
  lperiod = xgrid(nx)-xgrid(0)
  ALLOCATE(avermat(0:nx,dim))
  DO i=0,nx
     xq = xgrid(i) + rho*COS(th)
     xq = xgrid(0) + MODULO(xq-xgrid(0), lperiod)
     CALL allsplines(spl, xq, splines)
     DO j=1,dim
        avermat(i,j) = DOT_PRODUCT(wth, splines(:,j))/twopi
     END DO
  END DO
!
!   Gyro-averaged of F at grid points xgrid(0:nx)
  averf = MATMUL(avermat, coefs)
!
  DEALLOCATE(avermat)
END SUBROUTINE gyro