!>
!> @file fit1d_cmpl.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 a 1d complex function
!
  USE bsplines
  IMPLICIT NONE
  INTEGER,PARAMETER :: NX=10, NIDBAS=3, NPTS=40
  DOUBLE PRECISION :: pi, dx, xgrid(0:NX), xpt(NPTS), err
  DOUBLE COMPLEX, ALLOCATABLE :: coefs(:)
  DOUBLE COMPLEX :: fgrid(0:NX), fexact(NPTS), fcalc(NPTS)
  INTEGER :: dim, i
  TYPE(spline1d) :: spl
!================================================================================
!
!   Define grid and function values on grid
!
  pi = 4.0d0*ATAN(1.0d0)
  xgrid(0) = 0.0d0
  dx = 2.0d0*pi/NX
  DO i=1,NX
     xgrid(i) = xgrid(0) + i*dx
  END DO
!
  fgrid = func(xgrid)
!
  WRITE(*,'(2a10)') 'x', 'f'
  DO i=0,NX
     WRITE(*,'(3f10.4)') xgrid(i), fgrid(i)
  END DO
!
!   Set up spline
!
  CALL set_splcoef(NIDBAS, xgrid, spl, period=.TRUE.)
  CALL get_dim(spl, dim)
  ALLOCATE(coefs(dim))
!
!   Get Spline coefficients
!
  CALL get_splcoef(spl, fgrid, coefs)
  WRITE(*,'(a)') 'Interpolation coefs'
  DO i=1,dim
     WRITE(*,'(2(1pe12.3))')  coefs(i)
  END DO
!
!   Check interpolation
!
  CALL RANDOM_NUMBER(xpt)
  xpt = (2.0d0*pi) * xpt
  fexact = func(xpt)
!
  CALL gridval(spl, xpt, fcalc, 0, coefs)
!
  WRITE(*,'(a10,2a20)') 'x', 'fexact', 'fcacl'
  DO i=1,NPTS
     WRITE(*,'(5f10.4)') xpt(i), fexact(i), fcalc(i)
  END DO
  err = norm2(fcalc-fexact)
  WRITE(*,'(a,1pe12.3)') 'error', err
!
!   Clean up
!
  DEALLOCATE(coefs)
  CALL destroy_sp(spl)  
!================================================================================
CONTAINS
  FUNCTION func(x)
    DOUBLE PRECISION, INTENT(in) :: x(:)
    DOUBLE COMPLEX :: func(size(x))
    func = EXP( CMPLX(0.0d0, x))
  END FUNCTION func
!
  FUNCTION norm2(x)
!
!  Compute the 2-norm of vector x
!
    DOUBLE PRECISION :: norm2
    DOUBLE COMPLEX, INTENT(in) :: x(:)
!
    norm2 = SQRT(DOT_PRODUCT(x,x))
  END FUNCTION norm2
END PROGRAM main