subroutine lagcheby1_interp_1d ( nd, xd, yd, ni, xi, yi ) !*****************************************************************************80 ! !! lagcheby1_interp_1d() evaluates the Lagrange Chebyshev 1 interpolant. ! ! Discussion: ! ! The weight vector WD computed below is only valid if the data points ! XD are, as expected, the Chebyshev Type 1 points for [-1,+1], or a linearly ! mapped version for [A,B]. The XD values may be computed by: ! ! xd = r8vec_cheby1space ( nd, a, b ) ! ! for instance. ! ! Thanks to John Ferrier for point out that DENOM and NUMER were not ! being initialized, 16 September 2013. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 September 2021 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Jean-Paul Berrut, Lloyd Trefethen, ! Barycentric Lagrange Interpolation, ! SIAM Review, ! Volume 46, Number 3, September 2004, pages 501-517. ! ! Input: ! ! integer ND, the number of data points. ! ND must be at least 1. ! ! real ( kind = rk ) XD(ND), the data points. ! ! real ( kind = rk ) YD(ND), the data values. ! ! integer NI, the number of interpolation points. ! ! real ( kind = rk ) XI(NI), the interpolation points. ! ! Output: ! ! real ( kind = rk ) YI(NI), the interpolated values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer nd integer ni real ( kind = rk ) denom(ni) integer exact(ni) integer i integer j real ( kind = rk ) numer(ni) real ( kind = rk ) r8_mop real ( kind = rk ), parameter :: r8_pi = 3.141592653589793D+00 real ( kind = rk ) t real ( kind = rk ) theta real ( kind = rk ) wd real ( kind = rk ) xd(nd) real ( kind = rk ) xi(ni) real ( kind = rk ) yd(nd) real ( kind = rk ) yi(ni) exact(1:ni) = 0 numer(1:ni) = 0.0D+00 denom(1:ni) = 0.0D+00 do j = 1, nd theta = real ( 2 * j - 1, kind = rk ) * r8_pi / real ( 2 * nd, kind = rk ) wd = r8_mop ( j + 1 ) * sin ( theta ) do i = 1, ni if ( xi(i) == xd(j) ) then exact(i) = j numer(i) = yd(j) denom(i) = 1.0D+00 end if if ( exact(i) == 0 ) then t = wd / ( xi(i) - xd(j) ) numer(i) = numer(i) + t * yd(j) denom(i) = denom(i) + t end if end do end do yi(1:ni) = numer(1:ni) / denom(1:ni) return end subroutine lagcheby2_interp_1d ( nd, xd, yd, ni, xi, yi ) !*****************************************************************************80 ! !! lagcheby2_interp_1d() evaluates the Lagrange Chebyshev 2 interpolant. ! ! Discussion: ! ! The weight vector WD computed below is only valid if the data points ! XD are, as expected, the Chebyshev Type 2 points for [-1,+1], or a ! linearly mapped version for [A,B]. The XD values may be computed by: ! ! xd = r8vec_cheby2space ( nd, a, b ) ! ! for instance. ! ! Thanks to John Ferrier for point out that DENOM and NUMER were not ! being initialized, 16 September 2013. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 September 2021 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Jean-Paul Berrut, Lloyd Trefethen, ! Barycentric Lagrange Interpolation, ! SIAM Review, ! Volume 46, Number 3, September 2004, pages 501-517. ! ! Input: ! ! integer ND, the number of data points. ! ND must be at least 1. ! ! real ( kind = rk ) XD(ND), the data points. ! ! real ( kind = rk ) YD(ND), the data values. ! ! integer NI, the number of interpolation points. ! ! real ( kind = rk ) XI(NI), the interpolation points. ! ! Output: ! ! real ( kind = rk ) YI(NI), the interpolated values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer nd integer ni real ( kind = rk ) denom(ni) integer exact(ni) integer i integer j real ( kind = rk ) numer(ni) real ( kind = rk ) r8_mop real ( kind = rk ) t real ( kind = rk ) wd real ( kind = rk ) xd(nd) real ( kind = rk ) xi(ni) real ( kind = rk ) yd(nd) real ( kind = rk ) yi(ni) exact(1:ni) = 0 numer(1:ni) = 0.0D+00 denom(1:ni) = 0.0D+00 do j = 1, nd wd = r8_mop ( j + 1 ) if ( j == 1 .or. j == nd ) then wd = 0.5D+00 * wd end if do i = 1, ni if ( xi(i) == xd(j) ) then exact(i) = j numer(i) = yd(j) denom(i) = 1.0D+00 end if if ( exact(i) == 0 ) then t = wd / ( xi(i) - xd(j) ) numer(i) = numer(i) + t * yd(j) denom(i) = denom(i) + t end if end do end do yi(1:ni) = numer(1:ni) / denom(1:ni) return end subroutine lageven_interp_1d ( nd, xd, yd, ni, xi, yi ) !*****************************************************************************80 ! !! lageven_value_1d() evaluates the Lagrange evenly-spaced interpolant. ! ! Discussion: ! ! The weight vector WD computed below is only valid if the data points ! XD are, as expected, evenly spaced in an interval [A,B] with ! spacing (B-A)/N. The XD values might be computed by: ! ! xd(i) = ( ( 2 * nd - 2 * i + 1 ) * a ! + ( 2 * i - 1 ) * b ) ! / ( 2 * nd ) ! ! for instance. ! ! Thanks to John Ferrier for point out that DENOM and NUMER were not ! being initialized, 16 September 2013. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 September 2021 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Jean-Paul Berrut, Lloyd Trefethen, ! Barycentric Lagrange Interpolation, ! SIAM Review, ! Volume 46, Number 3, September 2004, pages 501-517. ! ! Input: ! ! integer ND, the number of data points. ! ND must be at least 1. ! ! real ( kind = rk ) XD(ND), the data points. ! ! real ( kind = rk ) YD(ND), the data values. ! ! integer NI, the number of interpolation points. ! ! real ( kind = rk ) XI(NI), the interpolation points. ! ! Output: ! ! real ( kind = rk ) YI(NI), the interpolated values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer nd integer ni real ( kind = rk ) denom(ni) integer exact(ni) integer i integer j real ( kind = rk ) numer(ni) real ( kind = rk ) r8_choose real ( kind = rk ) r8_mop real ( kind = rk ) t real ( kind = rk ) wd real ( kind = rk ) xd(nd) real ( kind = rk ) xi(ni) real ( kind = rk ) yd(nd) real ( kind = rk ) yi(ni) exact(1:ni) = 0 numer(1:ni) = 0.0D+00 denom(1:ni) = 0.0D+00 do j = 1, nd wd = r8_mop ( j ) * r8_choose ( nd, j ) do i = 1, ni if ( xi(i) == xd(j) ) then exact(i) = j numer(i) = yd(j) denom(i) = 1.0D+00 end if if ( exact(i) == 0 ) then t = wd / ( xi(i) - xd(j) ) numer(i) = numer(i) + t * yd(j) denom(i) = denom(i) + t end if end do end do yi(1:ni) = numer(1:ni) / denom(1:ni) return end