Homec4science
Diffusion SPClibs (feature/fft-split) Browse

spclibs/examples/fit1d_cmpl.f90
03540012fedefeature/fft-split

Last Change
Raw File

fit1d_cmpl.f90
View Options

!>
!> @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/>.
!>
!> @authors
!> (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
c4science · Help