!>
!> @file fit1dp.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
!
!  Fit grid value function to a spline of any order
!  Periodic case.
!
  USE bsplines
  USE futils
!
  IMPLICIT NONE
  INTEGER :: nx, nidbas
  DOUBLE PRECISION :: a, b, coefx(5)
  DOUBLE PRECISION, ALLOCATABLE ::  xgrid(:), fgrid(:,:), coefs(:,:)
  INTEGER :: i, dim
  TYPE(spline1d) :: spl
  DOUBLE PRECISION :: dx, x0, x1
  INTEGER :: npts
  DOUBLE PRECISION, ALLOCATABLE :: xpt(:), fcalc(:), fexact(:), err(:)
  DOUBLE PRECISION, POINTER :: splines(:,:) => null()
!
  CHARACTER(len=128) :: file='fit1d.h5'
  INTEGER :: fid
!
  NAMELIST /newrun/ nx, nidbas, a, b, coefx
!===========================================================================
!              1.0 Prologue
!
!   Read in data specific to run
!
  nx = 10    ! Number of intevals in x
  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
!
  READ(*,newrun)
  WRITE(*,newrun)
!
!   Define grid and function values
!
  ALLOCATE(xgrid(0:nx), fgrid(0:nx,1))
  xgrid(0) = a
  xgrid(nx) = b
  CALL meshdist(coefx, xgrid, nx)
  fgrid(:,1) = func(xgrid)
  WRITE(*,'(a/(10f8.3))') 'xgrid', xgrid
  WRITE(*,'(a/(10f8.3))') 'fgrid', fgrid
!
!   Create hdf5 file
!
  CALL creatf(file, fid, 'FIT1D Result File')
  CALL attach(fid, '/', 'NX', nx)
  CALL attach(fid, '/', 'NIDBAS', nidbas)
!===========================================================================
!              2.0 Spline interpolation
!
!   Set up the spline interpolation
  CALL set_splcoef(nidbas, xgrid, spl, period=.TRUE.)
  CALL get_dim(spl, dim)
  WRITE(*,'(2(a,i5, 2x))') 'NX =', nx, 'DIM =', dim
  WRITE(*,'(/a,i6,a,i6/(10(f8.3)))') 'KNOTS of splines', LBOUND(spl%knots), &
       &                       ':',UBOUND(spl%knots), spl%knots
!
  ALLOCATE(coefs(dim,1))
!
!   From given grid values fgrid, compute the spline coefs
  CALL get_splcoef(spl, fgrid, coefs)
  WRITE(*,'(a/(10f8.3))') 'coefs', coefs
!
!   Plot all splines
!
  npts = 100
  ALLOCATE(xpt(npts))
!!$  x0 = a
!!$  x1 = b
  x0 = spl%knots(0)
  x1 = spl%knots(nx)
  dx = (x1 -x0)/REAL(npts)  ! Last point b not inluded
  DO i=1,npts
     xpt(i) = x0 + (i-1)*dx
  END DO
  CALL allsplines(spl, xpt, splines)
  CALL putarr(fid, '/X', xpt)
  CALL putarr(fid,'/SPLINES', splines)
!
!   Check interpolation
!
  ALLOCATE(fcalc(npts), fexact(npts), err(npts))
  fexact = func(xpt)
!
  CALL gridval(spl, xpt, fcalc, 0, coefs(:,1))
  err = fexact - fcalc
  CALL putarr(fid, '/FEXACT', fexact, 'Exact values')
  CALL putarr(fid, '/FCALC', fcalc, 'Interpolated values')
  CALL putarr(fid, '/ERROR', err, 'Interpolation errors')
  WRITE(*,'(/a,1pe12.3)') 'Norm of error', norm2(err)/norm2(fexact)
!
!  Derivatives values
  CALL gridval(spl, xpt, fcalc, 1)
  fexact = func1(xpt)
  err = fexact - fcalc
  CALL putarr(fid, '/FEXACT1', fexact, 'Exact values of first derivative')
  CALL putarr(fid, '/FCALC1', fcalc, 'Interpolated first derivative')
  CALL putarr(fid, '/ERROR1', err, 'Interpolation errors on first derivative')
  WRITE(*,'(/a,1pe12.3)') 'Norm of error', norm2(err)/norm2(fexact)
!===========================================================================
!              9.0  Epilogue
!
  DEALLOCATE(xpt, splines, fexact, fcalc, err)
  DEALLOCATE(xgrid, fgrid)
  CALL destroy_sp(spl)
  CALL closef(fid)
!===========================================================================
!
CONTAINS
  FUNCTION func(x)
    DOUBLE PRECISION, INTENT(in) :: x(:)
    DOUBLE PRECISION :: func(SIZE(x))
    DOUBLE PRECISION :: pi
    pi = 4.0*ATAN(1.0d0)
    func = SIN(2.d0*pi*x) + 2.0d0*COS(8.d0*pi*x)
  END FUNCTION func
  FUNCTION func1(x)
    DOUBLE PRECISION, INTENT(in) :: x(:)
    DOUBLE PRECISION :: func1(SIZE(x))
    DOUBLE PRECISION :: pi
    pi = 4.0*ATAN(1.0d0)
    func1 = 2.d0*pi*COS(2.d0*pi*x) - 16.0d0*pi*SIN(8.d0*pi*x)
  END FUNCTION func1
!
  FUNCTION norm2(x)
!
!  Compute the 2-norm of array x
!
    DOUBLE PRECISION :: norm2
    DOUBLE PRECISION, INTENT(in) :: x(:)
    DOUBLE PRECISION :: sum2
    INTEGER :: i
!
    sum2 = 0.0d0
    DO i=1,SIZE(x,1)
       sum2 = sum2 + x(i)**2
    END DO
    norm2 = SQRT(sum2)
  END FUNCTION norm2
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
!+++