!>
!> @file fit2dbc.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 2d grid value function to a 2d spline of any order
!  BC using derivatives at both ends, in the non-periodic direction.
!
  USE bsplines
  USE futils
!
  IMPLICIT NONE
  CHARACTER(len=128) :: file='fit2d.h5'
  INTEGER :: fid
  INTEGER :: nx, ny, nidbas(2), mbes, dims(2)
  INTEGER, PARAMETER :: nptx=100, npty=100
  DOUBLE PRECISION :: pi, coefx(5), coefy(5)
  DOUBLE PRECISION, ALLOCATABLE :: xgrid(:),ygrid(:),fgrid(:,:),bcoefs(:,:)
  DOUBLE PRECISION :: dx, dy, xpt(nptx), ypt(npty)
  DOUBLE PRECISION, DIMENSION(nptx,npty) :: fcalc, fexact, errs
  TYPE(spline2d) :: splxy
  DOUBLE PRECISION :: mem
  INTEGER :: i, j, ii
  INTEGER :: ibc1(2,10), ibc2(2,10)
  DOUBLE PRECISION, ALLOCATABLE:: fbc(:,:,:)
!
  NAMELIST /newrun/ nx, ny, nidbas, mbes, coefx, coefy

!===========================================================================
!              1.0 Prologue
!
!   Read in data specific to run
!
  nx = 8              ! Number of intevals in x
  ny = 8              ! Number of intevals in y
  nidbas = (/3,3/)    ! Degree of splines
  mbes = 2
  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)
!
!   Define grid on x axis
!
  pi = 4.0d0*ATAN(1.0d0)
  ALLOCATE(xgrid(0:nx), ygrid(0:ny), fgrid(0:nx,0:ny))
  xgrid(0) = 0.0d0;   xgrid(nx) = 1.0d0
  CALL meshdist(coefx, xgrid, nx)
  ygrid(0) = 0.0d0;   ygrid(ny) = 1.0d0
  CALL meshdist(coefy, ygrid, ny)
  fgrid = func(xgrid,ygrid)
!
  WRITE(*,'(a/(12f8.3))') 'XGRID', xgrid(0:nx)
  WRITE(*,'(a/(12f8.3))') 'YGRID', ygrid(0:ny)
  IF( nx.LE.10 .AND. ny.LE.10 ) THEN
     WRITE(*,'(a)') 'FGRID'
     DO j=0,ny
        WRITE(*,'(12f8.3)') fgrid(:,j)
     END DO
     WRITE(*,'(a,f8.3)') 'Memory used so far (MB) =', mem()
  END IF
!
!   Create hdf5 file
!
  CALL creatf(file, fid, 'PDE2D Result File')
  CALL attach(fid, '/', 'NX', nx)
  CALL attach(fid, '/', 'NY', ny)
  CALL attach(fid, '/', 'NIDBAS1', nidbas(1))
  CALL attach(fid, '/', 'NIDBAS2', nidbas(2))
  CALL attach(fid, '/', 'MBESS', mbes)
!===========================================================================
!              2.0 Spline interpolation
!
!   Setup the spline interpolation
  ii=1           ! Start with first derivative
  DO i = 1, nidbas(1)/2
     ibc1(1,i) = ii+i-1
     ibc1(2,i) = ii+i-1
  END DO
  ibc2 = ibc1
  CALL set_splcoef(nidbas, xgrid, ygrid, splxy, (/.FALSE., .FALSE./),&
       &           ibc1=ibc1, ibc2=ibc2)
  WRITE(*,'(a,f8.3)') 'Memory used so far (MB) =', mem()
!
!   Compute spline interpolation coefficients
  CALL get_dim(splxy, dims)
  ALLOCATE(bcoefs(dims(1),dims(2)))
  WRITE(*,'(a,2i4)') 'Dims of spline', dims
!
  ALLOCATE(fbc(2, nidbas(1)/2, 0:ny))
  fbc=0.0d0
!
  WRITE(*,'(a/(10f8.3))') 'fbc(1)', fbc(1,1,:)
  WRITE(*,'(a/(10f8.3))') 'fbc(2)', fbc(2,1,:)
!
  CALL get_splcoef(splxy, fgrid, bcoefs, fbc1=fbc, fbc2=fbc)
!
  DEALLOCATE(fbc)
!===========================================================================
!              2.0 Check interpolation
!
  dx=(xgrid(nx)-xgrid(0))/REAL(nptx-1)
  dy=(ygrid(ny)-ygrid(0))/REAL(npty-1)
  DO i=1,nptx
     xpt(i) = xgrid(0) + (i-1)*dx
  END DO
  DO i=1,npty
     ypt(i) = ygrid(0) + (i-1)*dy
  END DO
  fexact = func(xpt,ypt)
  CALL gridval(splxy, xpt, ypt, fcalc, (/0,0/), bcoefs)
  errs = fcalc-fexact
  WRITE(*,'(a,2(1pe12.3))') 'Min max or errors', MINVAL(errs), MAXVAL(errs)
!
  CALL putarr(fid, '/xpt', xpt, 'r')
  CALL putarr(fid, '/ypt', ypt, '\theta')
  CALL putarr(fid, '/bcoefs', bcoefs, 'bcoefs')
  CALL putarr(fid, '/fcalc', fcalc, 'Interpolated')
  CALL putarr(fid, '/fexact', fexact, 'Exact')
  CALL putarr(fid, '/errs', errs, 'Errors')
!===========================================================================
!              9.0  Epilogue
!
  DEALLOCATE(bcoefs)
  DEALLOCATE(xgrid, ygrid, fgrid)
  CALL destroy_sp(splxy)
  CALL closef(fid)
!
!===========================================================================
CONTAINS
  FUNCTION func(x,y)
!
!   A function with zeo derivatives at both ends
!
    DOUBLE PRECISION, INTENT(in) :: x(:), y(:)
    DOUBLE PRECISION :: func(SIZE(x), SIZE(y))
    DOUBLE PRECISION :: zy
    INTEGER :: j
    DO j=1,SIZE(y)
       zy = y(j)*y(j)*(y(j)-1.5d0)
       func(:,j) = x(:)*x(:)*(x(:)-1.5d0) + zy
    END DO
  END FUNCTION func
!
  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