!>
!> @file fit2dbc_y.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/>.
!>
!> @authors
!> (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.
!
! Testing BC on derivative along second 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(npty,nptx) :: fcalc, fexact, errs
  TYPE(spline2d) :: splxy
  DOUBLE PRECISION :: mem
  INTEGER :: i, j, ii
  INTEGER :: ibc(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:ny,0:nx))
  xgrid(0) = 0.0d0;   xgrid(nx) = 1.0d0
  CALL meshdist(coefx, xgrid, nx)
  ygrid(0) = 0.0d0;   ygrid(ny) = 2.d0*pi
  CALL meshdist(coefy, ygrid, ny)
  fgrid = TRANSPOSE(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(y, x)'
     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
     ibc(1,i) = ii+i-1
     ibc(2,i) = ii+i-1
  END DO
  CALL set_splcoef((/nidbas(2), nidbas(1)/), ygrid, xgrid, splxy, (/.TRUE., .FALSE./),&
       &           ibc2=ibc)
  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
!
!!$!   Exact first derivatives at boundaries
!!$  fbc(1,1:1,:) = func1(xgrid(0:0), ygrid(0:ny))
!!$  fbc(2,1:1,:) = func1(xgrid(nx:nx), ygrid(0:ny))
!
!!$!   Derivatives at boundaries approximated with FD
!!$  DO j=0,ny
!!$     fbc(1,1,j+1) = fgrid(1,j)-fgrid(0,j)
!!$     fbc(2,1,j+1) = fgrid(nx,j)-fgrid(nx-1,j)
!!$  END DO
!!$  fbc(1,1,:) = fbc(1,1,:)/(xgrid(1)-xgrid(0))
!!$  fbc(2,1,:) = fbc(2,1,:)/(xgrid(nx)-xgrid(nx-1))
!
  WRITE(*,'(a/(10f8.3))') 'fbc(1)', fbc(1,1,:)
  WRITE(*,'(a/(10f8.3))') 'fbc(2)', fbc(2,1,:)
!
  CALL get_splcoef(splxy, fgrid, bcoefs, 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 = TRANSPOSE(func(xpt,ypt))
  CALL gridval(splxy, ypt, xpt, 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, '/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)
    DOUBLE PRECISION, INTENT(in) :: x(:), y(:)
    DOUBLE PRECISION :: func(SIZE(x), SIZE(y))
    DOUBLE PRECISION :: zy
    INTEGER :: j
    DO j=1,SIZE(y)
       zy = -mbes * SIN(mbes*y(j))
       func(:,j) =(1-x(:)**2) * x(:)**mbes * zy
    END DO
  END FUNCTION func
  FUNCTION func1(x,y)
    DOUBLE PRECISION, INTENT(in) :: x(:), y(:)
    DOUBLE PRECISION :: func1(SIZE(x), SIZE(y))
    DOUBLE PRECISION :: zy
    INTEGER :: j
    DO j=1,SIZE(y)
       zy = -mbes * SIN(mbes*y(j))
       func1(:,j) = (mbes - (mbes+2.0d0)*x(:)**2) * x(:)**(mbes-1) * zy
    END DO
  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