!>
!> @file driv1.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
!
! Basis splines on a 2d grid.
!
  USE bsplines
  USE futils
  IMPLICIT NONE
  TYPE(spline1d) :: spx, spy
  INTEGER :: nx=10, ny=8, nidbas=2, ngauss=4, npts=1000
  DOUBLE PRECISION :: a, b, coefx(5), coefy(5)
  DOUBLE PRECISION, ALLOCATABLE :: xgrid(:), ygrid(:), fun(:,:)
  DOUBLE PRECISION, ALLOCATABLE :: xp(:), funxp(:,:), yp(:), funyp(:,:)
  DOUBLE PRECISION :: dx, dy
  INTEGER :: i, j, left
  CHARACTER(len=256) :: title
  INTEGER :: fid
  NAMELIST /newrun/ nx, ny, nidbas, ngauss, a, b, coefx, coefy
!===========================================================================
  nidbas = 3
  ngauss = 4
  nx = 10    ! Number of intevals in x
  ny = 8     ! Number of intevals in y
  a = 0.0d0  ! Left boundary of interval
  b = 1.0d0  ! Right boundary of interval
  coefx(1:5) = (/1.0d0, 0.d0, 0.d0, 0.d0, 1.d0/) ! Mesh point distribution function
  coefy(1:5) = (/1.0d0, 0.d0, 0.d0, 0.d0, 1.d0/) ! Mesh point distribution function
  READ(*,newrun)
  WRITE(*,newrun)
!===========================================================================
!                  1.0 Set up grids
!
  ALLOCATE( xgrid(0:nx) )
  xgrid(0) = a
  xgrid(nx) = b
  CALL meshdist(coefx, xgrid, nx)
!
!!$  dy = 2.d0*pi/REAL(ny)
  dy = 1.0d0
  ALLOCATE( ygrid(0:ny) )
  ygrid(0) = a
  ygrid(ny) = b
  CALL meshdist(coefy, ygrid, ny)
!===========================================================================
!                  2.0 Set up splines on (x,y)
!
  CALL set_spline(nidbas, ngauss, xgrid, spx)
  CALL set_spline(nidbas, ngauss, ygrid, spy, period=.TRUE.)
  WRITE(*,'(/a,i6,a,i6/(10(f8.3)))') 'KNOTS of x-splines', LBOUND(spx%knots), &
       &                       ':',UBOUND(spx%knots), spx%knots
  WRITE(*,'(2(a,i5, 2x))') 'NX =', nx, 'DIM =', spx%dim
  WRITE(*,'(/a,i6,a,i6/(10(f8.3)))') 'KNOTS of y-splines', LBOUND(spy%knots), &
       &                       ':',UBOUND(spy%knots), spy%knots
  WRITE(*,'(2(a,i5, 2x))') 'NY =', ny, 'DIM =', spy%dim
!===========================================================================
!                  3.0 Graph the splines on (x,y)
!
  ALLOCATE( fun(nidbas+1,1) )         ! Only 0-th derivative
!!$  ALLOCATE( fun(nidbas+1,0:1) )         ! Only 0-th derivative
  ALLOCATE( xp(npts), funxp(npts,0:spx%dim-1) )
  ALLOCATE( yp(npts-1), funyp(npts-1,0:spy%dim-1) )
!
!   Splines in X (non-peridic)
!
  WRITE(*,'(a)') 'Splines in x'
  dx = (xgrid(nx)-xgrid(0)) / REAL(NPTS-1)
  DO i=1,npts
     xp(i) = xgrid(0) + (i-1)*dx
     CALL locintv(spx, xp(i), left)
     CALL basfun(xp(i), spx, fun, left+1)
     funxp(i,left:left+nidbas) = fun(:,1)
  END DO
!
!   Splines in Y (periodic)
!
  WRITE(*,'(a)') 'Splines in y'
  dy = (ygrid(ny)-ygrid(0)) / REAL(NPTS-1)
  DO i=1,npts-1
     yp(i) = ygrid(0) + (i-1)*dy
     CALL locintv(spy, yp(i), left)
     CALL basfun(yp(i), spy, fun, left+1)
     funyp(i,left:left+nidbas) = fun(:,1)
  END DO
!
!  Create hdf5 file
!
  CALL creatf('driv1.h5', fid, real_prec='d')
!
  WRITE(title,'(a,i3,5x,a,i6)') 'Splines of degree =', nidbas, 'NX =', nx
  CALL putarr(fid, 'X', xp)
  CALL putarr(fid, 'KNOTSX', spx%knots)
  CALL putarr(fid, 'splinesx', funxp, TRIM(title))
  CALL putarr(fid, 'KNOTSY', spy%knots)
!
  WRITE(title,'(a,i3,5x,a,i6)') 'Periodic splines of degree =', nidbas, 'NY =', ny
  CALL putarr(fid, 'Y', yp)
  CALL putarr(fid, 'splinesy', funyp, TRIM(title))
  CALL closef(fid)
!===========================================================================
!                  9.0 Epilogue
!
  DEALLOCATE(fun)
  DEALLOCATE(xgrid, ygrid)
  DEALLOCATE(xp, funxp)
  CALL destroy_sp(spx)
  CALL destroy_sp(spy)
END PROGRAM main
!+++
SUBROUTINE meshdist(c, x, nx)
!
!   Construct an 1d  non-equidistant mesh given a
!   mesh distribution function.
!
  IMPLICIT NONE
  DOUBLE PRECISION, INTENT(in) :: c(5)
  INTEGER, INTENT(iN) :: nx
  DOUBLE PRECISION, INTENT(inout) :: x(0:nx)
  INTEGER :: nintg
  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xint, fint
  DOUBLE PRECISION :: a, b, dx, f0, f1, scal
  INTEGER :: i, k
!
  a=x(0)
  b=x(nx)
  nintg = 10*nx
  ALLOCATE(xint(0:nintg), fint(0:nintg))
!
!  Mesh distribution
!
  dx = (b-a)/REAL(nintg)
  xint(0) = a
  fint(0) = 0.0d0
  f1 = fdist(xint(0))
  DO i=1,nintg
     f0 = f1
     xint(i) = xint(i-1) + dx
     f1 = fdist(xint(i))
     fint(i) = fint(i-1) + 0.5*(f0+f1)
  END DO
!
!  Normalization
!
  scal = REAL(nx) / fint(nintg)
  fint(0:nintg) = fint(0:nintg) * scal
!!$  WRITE(*,'(a/(10f8.3))') 'FINT', fint
!
!  Obtain mesh point by (inverse) interpolation
!
  k = 1
  DO i=1,nintg-1
     IF( fint(i) .GE. REAL(k) ) THEN
        x(k) = xint(i) + (xint(i+1)-xint(i))/(fint(i+1)-fint(i)) * &
             &   (k-fint(i))
        k = k+1
     END IF
  END DO
!
  DEALLOCATE(xint, fint)
CONTAINS
  DOUBLE PRECISION FUNCTION fdist(x)
    DOUBLE PRECISION, INTENT(in) :: x
    fdist = c(1) + c(2)*x + c(3)*EXP(-((x-c(4))/c(5))**2)
  END FUNCTION fdist
END SUBROUTINE meshdist
!+++