module m_spline !! Implements various routines for cubic spline interpolation. implicit none contains !*****************************************************************************80 subroutine basis_function_b_val ( tdata, tval, yval ) !*****************************************************************************80 ! !! BASIS_FUNCTION_B_VAL evaluates the B spline basis function. ! ! Discussion: ! ! The B spline basis function is a piecewise cubic which ! has the properties that: ! ! * it equals 2/3 at TDATA(3), 1/6 at TDATA(2) and TDATA(4); ! * it is 0 for TVAL <= TDATA(1) or TDATA(5) <= TVAL; ! * it is strictly increasing from TDATA(1) to TDATA(3), ! and strictly decreasing from TDATA(3) to TDATA(5); ! * the function and its first two derivatives are continuous ! at each node TDATA(I). ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 06 April 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Alan Davies, Philip Samuels, ! An Introduction to Computational Geometry for Curves and Surfaces, ! Clarendon Press, 1996, ! ISBN: 0-19-851478-6, ! LC: QA448.D38. ! ! Input: ! ! real ( kind = 8 ) TDATA(5), the nodes associated with the ! basis function. The entries of TDATA are assumed to be distinct ! and increasing. ! ! real ( kind = 8 ) TVAL, a point at which the B spline basis ! function is to be evaluated. ! ! Output: ! ! real ( kind = 8 ) YVAL, the value of the function at TVAL. ! implicit none integer ( kind = 4 ), parameter :: ndata = 5 integer ( kind = 4 ) left integer ( kind = 4 ) right real ( kind = 8 ) tdata(ndata) real ( kind = 8 ) tval real ( kind = 8 ) u real ( kind = 8 ) yval if ( tval <= tdata(1) .or. tdata(ndata) <= tval ) then yval = 0.0D+00 return end if ! ! Find the interval [ TDATA(LEFT), TDATA(RIGHT) ] containing TVAL. ! call r8vec_bracket ( ndata, tdata, tval, left, right ) ! ! U is the normalized coordinate of TVAL in this interval. ! u = ( tval - tdata(left) ) / ( tdata(right) - tdata(left) ) ! ! Now evaluate the function. ! if ( tval < tdata(2) ) then yval = u**3 / 6.0D+00 else if ( tval < tdata(3) ) then yval = ( ( ( - 3.0D+00 & * u + 3.0D+00 ) & * u + 3.0D+00 ) & * u + 1.0D+00 ) / 6.0D+00 else if ( tval < tdata(4) ) then yval = ( ( ( + 3.0D+00 & * u - 6.0D+00 ) & * u + 0.0D+00 ) & * u + 4.0D+00 ) / 6.0D+00 else if ( tval < tdata(5) ) then yval = ( 1.0D+00 - u )**3 / 6.0D+00 end if return end subroutine basis_function_beta_val ( beta1, beta2, tdata, tval, yval ) !*****************************************************************************80 ! !! BASIS_FUNCTION_BETA_VAL evaluates the beta spline basis function. ! ! Discussion: ! ! With BETA1 = 1 and BETA2 = 0, the beta spline basis function ! equals the B spline basis function. ! ! With BETA1 large, and BETA2 = 0, the beta spline basis function ! skews to the right, that is, its maximum increases, and occurs ! to the right of the center. ! ! With BETA1 = 1 and BETA2 large, the beta spline becomes more like ! a linear basis function; that is, its value in the outer two intervals ! goes to zero, and its behavior in the inner two intervals approaches ! a piecewise linear function that is 0 at the second node, 1 at the ! third, and 0 at the fourth. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 09 April 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Alan Davies, Philip Samuels, ! An Introduction to Computational Geometry for Curves and Surfaces, ! Clarendon Press, 1996, ! ISBN: 0-19-851478-6, ! LC: QA448.D38. ! ! Parameters: ! ! Input, real ( kind = 8 ) BETA1, the skew or bias parameter. ! BETA1 = 1 for no skew or bias. ! ! Input, real ( kind = 8 ) BETA2, the tension parameter. ! BETA2 = 0 for no tension. ! ! Input, real ( kind = 8 ) TDATA(5), the nodes associated with the ! basis function. The entries of TDATA are assumed to be distinct ! and increasing. ! ! Input, real ( kind = 8 ) TVAL, a point at which the B spline ! basis function is to be evaluated. ! ! Output, real ( kind = 8 ) YVAL, the value of the function at TVAL. ! implicit none integer ( kind = 4 ), parameter :: ndata = 5 real ( kind = 8 ) a real ( kind = 8 ) b real ( kind = 8 ) beta1 real ( kind = 8 ) beta2 real ( kind = 8 ) c real ( kind = 8 ) d integer ( kind = 4 ) left integer ( kind = 4 ) right real ( kind = 8 ) tdata(ndata) real ( kind = 8 ) tval real ( kind = 8 ) u real ( kind = 8 ) yval if ( tval <= tdata(1) .or. tdata(ndata) <= tval ) then yval = 0.0D+00 return end if ! ! Find the interval [ TDATA(LEFT), TDATA(RIGHT) ] containing TVAL. ! call r8vec_bracket ( ndata, tdata, tval, left, right ) ! ! U is the normalized coordinate of TVAL in this interval. ! u = ( tval - tdata(left) ) / ( tdata(right) - tdata(left) ) ! ! Now evaluate the function. ! if ( tval < tdata(2) ) then yval = 2.0D+00 * u**3 else if ( tval < tdata(3) ) then a = beta2 + 4.0D+00 * beta1 + 4.0D+00 * beta1 * beta1 & + 6.0D+00 * ( 1.0D+00 - beta1 * beta1 ) & - 3.0D+00 * ( 2.0D+00 + beta2 + 2.0D+00 * beta1 ) & + 2.0D+00 * ( 1.0D+00 + beta2 + beta1 + beta1 * beta1 ) b = - 6.0D+00 * ( 1.0D+00 - beta1 * beta1 ) & + 6.0D+00 * ( 2.0D+00 + beta2 + 2.0D+00 * beta1 ) & - 6.0D+00 * ( 1.0D+00 + beta2 + beta1 + beta1 * beta1 ) c = - 3.0D+00 * ( 2.0D+00 + beta2 + 2.0D+00 * beta1 ) & + 6.0D+00 * ( 1.0D+00 + beta2 + beta1 + beta1 * beta1 ) d = - 2.0D+00 * ( 1.0D+00 + beta2 + beta1 + beta1 * beta1 ) yval = ( ( d * u + c ) * u + b ) * u + a else if ( tval < tdata(4) ) then a = beta2 + 4.0D+00 * beta1 + 4.0D+00 * beta1 * beta1 b = - 6.0D+00 * beta1 * ( 1.0D+00 - beta1 * beta1 ) c = - 3.0D+00 * ( beta2 + 2.0D+00 * beta1**2 + 2.0D+00 * beta1**3 ) d = 2.0D+00 * ( beta2 + beta1 + beta1**2 + beta1**3 ) yval = ( ( d * u + c ) * u + b ) * u + a else if ( tval < tdata(5) ) then yval = 2.0D+00 * beta1**3 * ( 1.0D+00 - u )**3 end if yval = yval / ( 2.0D+00 + beta2 + 4.0D+00 * beta1 + 4.0D+00 * beta1**2 & + 2.0D+00 * beta1**3 ) return end subroutine basis_matrix_b_uni ( mbasis ) !*****************************************************************************80 ! !! BASIS_MATRIX_B_UNI sets up the uniform B spline basis matrix. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 17 June 2006 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! James Foley, Andries vanDam, Steven Feiner, John Hughes, ! Computer Graphics, Principles and Practice, ! Second Edition, ! Addison Wesley, 1995, ! ISBN: 0201848406, ! LC: T385.C5735. ! ! Parameters: ! ! Output, real ( kind = 8 ) MBASIS(4,4), the basis matrix. ! implicit none real ( kind = 8 ) mbasis(4,4) ! ! In the following statement, the matrix appears as though it ! has been transposed. ! mbasis(1:4,1:4) = real ( reshape ( & (/ -1, 3, -3, 1, & 3, -6, 0, 4, & -3, 3, 3, 1, & 1, 0, 0, 0 /), & (/ 4, 4 /) ), kind = 8 ) / 6.0D+00 return end subroutine basis_matrix_beta_uni ( beta1, beta2, mbasis ) !*****************************************************************************80 ! !! BASIS_MATRIX_BETA_UNI sets up the uniform beta spline basis matrix. ! ! Discussion: ! ! If BETA1 = 1 and BETA2 = 0, then the beta spline reduces to ! the B spline. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 12 February 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! James Foley, Andries vanDam, Steven Feiner, John Hughes, ! Computer Graphics, Principles and Practice, ! Second Edition, ! Addison Wesley, 1995, ! ISBN: 0201848406, ! LC: T385.C5735. ! ! Parameters: ! ! Input, real ( kind = 8 ) BETA1, the skew or bias parameter. ! BETA1 = 1 for no skew or bias. ! ! Input, real ( kind = 8 ) BETA2, the tension parameter. ! BETA2 = 0 for no tension. ! ! Output, real ( kind = 8 ) MBASIS(4,4), the basis matrix. ! implicit none real ( kind = 8 ) beta1 real ( kind = 8 ) beta2 real ( kind = 8 ) delta real ( kind = 8 ) mbasis(4,4) mbasis(1,1) = - 2.0D+00 * beta1 * beta1 * beta1 mbasis(1,2) = 2.0D+00 * beta2 & + 2.0 * beta1 * ( beta1 * beta1 + beta1 + 1.0D+00 ) mbasis(1,3) = - 2.0D+00 * ( beta2 + beta1 * beta1 + beta1 + 1.0D+00 ) mbasis(1,4) = 2.0D+00 mbasis(2,1) = 6.0D+00 * beta1 * beta1 * beta1 mbasis(2,2) = - 3.0D+00 * beta2 & - 6.0D+00 * beta1 * beta1 * ( beta1 + 1.0D+00 ) mbasis(2,3) = 3.0D+00 * beta2 + 6.0D+00 * beta1 * beta1 mbasis(2,4) = 0.0D+00 mbasis(3,1) = - 6.0D+00 * beta1 * beta1 * beta1 mbasis(3,2) = 6.0D+00 * beta1 * ( beta1 - 1.0D+00 ) * ( beta1 + 1.0D+00 ) mbasis(3,3) = 6.0D+00 * beta1 mbasis(3,4) = 0.0D+00 mbasis(4,1) = 2.0D+00 * beta1 * beta1 * beta1 mbasis(4,2) = 4.0D+00 * beta1 * ( beta1 + 1.0D+00 ) + beta2 mbasis(4,3) = 2.0D+00 mbasis(4,4) = 0.0D+00 delta = ( ( 2.0D+00 & * beta1 + 4.0D+00 ) & * beta1 + 4.0D+00 ) & * beta1 + 2.0D+00 + beta2 mbasis(1:4,1:4) = mbasis(1:4,1:4) / delta return end subroutine basis_matrix_bezier ( mbasis ) !*****************************************************************************80 ! !! BASIS_MATRIX_BEZIER sets up the cubic Bezier spline basis matrix. ! ! Discussion: ! ! This basis matrix assumes that the data points are stored as ! ( P1, P2, P3, P4 ). P1 is the function value at T = 0, while ! P2 is used to approximate the derivative at T = 0 by ! dP/dt = 3 * ( P2 - P1 ). Similarly, P4 is the function value ! at T = 1, and P3 is used to approximate the derivative at T = 1 ! by dP/dT = 3 * ( P4 - P3 ). ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 05 April 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! James Foley, Andries vanDam, Steven Feiner, John Hughes, ! Computer Graphics, Principles and Practice, ! Second Edition, ! Addison Wesley, 1995, ! ISBN: 0201848406, ! LC: T385.C5735. ! ! Parameters: ! ! Output, real ( kind = 8 ) MBASIS(4,4), the basis matrix. ! implicit none real ( kind = 8 ) mbasis(4,4) mbasis(1,1) = -1.0D+00 mbasis(1,2) = 3.0D+00 mbasis(1,3) = -3.0D+00 mbasis(1,4) = 1.0D+00 mbasis(2,1) = 3.0D+00 mbasis(2,2) = -6.0D+00 mbasis(2,3) = 3.0D+00 mbasis(2,4) = 0.0D+00 mbasis(3,1) = -3.0D+00 mbasis(3,2) = 3.0D+00 mbasis(3,3) = 0.0D+00 mbasis(3,4) = 0.0D+00 mbasis(4,1) = 1.0D+00 mbasis(4,2) = 0.0D+00 mbasis(4,3) = 0.0D+00 mbasis(4,4) = 0.0D+00 return end subroutine basis_matrix_hermite ( mbasis ) !*****************************************************************************80 ! !! BASIS_MATRIX_HERMITE sets up the Hermite spline basis matrix. ! ! Discussion: ! ! This basis matrix assumes that the data points are stored as ! ( P1, P2, P1', P2' ), with P1 and P1' being the data value and ! the derivative dP/dT at T = 0, while P2 and P2' apply at T = 1. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 06 April 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! James Foley, Andries vanDam, Steven Feiner, John Hughes, ! Computer Graphics, Principles and Practice, ! Second Edition, ! Addison Wesley, 1995, ! ISBN: 0201848406, ! LC: T385.C5735. ! ! Parameters: ! ! Output, real ( kind = 8 ) MBASIS(4,4), the basis matrix. ! implicit none real ( kind = 8 ) mbasis(4,4) mbasis(1,1) = 2.0D+00 mbasis(1,2) = -2.0D+00 mbasis(1,3) = 1.0D+00 mbasis(1,4) = 1.0D+00 mbasis(2,1) = -3.0D+00 mbasis(2,2) = 3.0D+00 mbasis(2,3) = -2.0D+00 mbasis(2,4) = -1.0D+00 mbasis(3,1) = 0.0D+00 mbasis(3,2) = 0.0D+00 mbasis(3,3) = 1.0D+00 mbasis(3,4) = 0.0D+00 mbasis(4,1) = 1.0D+00 mbasis(4,2) = 0.0D+00 mbasis(4,3) = 0.0D+00 mbasis(4,4) = 0.0D+00 return end subroutine basis_matrix_overhauser_nonuni ( alpha, beta, mbasis ) !*****************************************************************************80 ! !! BASIS_MATRIX_OVERHAUSER_NONUNI: nonuniform Overhauser spline basis matrix. ! ! Discussion: ! ! This basis matrix assumes that the data points P1, P2, P3 and ! P4 are not uniformly spaced in T, and that P2 corresponds to T = 0, ! and P3 to T = 1. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 05 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) ALPHA, BETA. ! ALPHA = || P2 - P1 || / ( || P3 - P2 || + || P2 - P1 || ) ! BETA = || P3 - P2 || / ( || P4 - P3 || + || P3 - P2 || ). ! ! Output, real ( kind = 8 ) MBASIS(4,4), the basis matrix. ! implicit none real ( kind = 8 ) alpha real ( kind = 8 ) beta real ( kind = 8 ) mbasis(4,4) mbasis(1,1) = - ( 1.0D+00 - alpha ) * ( 1.0D+00 - alpha ) / alpha mbasis(1,2) = beta + ( 1.0D+00 - alpha ) / alpha mbasis(1,3) = alpha - 1.0D+00 / ( 1.0D+00 - beta ) mbasis(1,4) = beta * beta / ( 1.0D+00 - beta ) mbasis(2,1) = 2.0D+00 * ( 1.0D+00 - alpha ) * ( 1.0D+00 - alpha ) / alpha mbasis(2,2) = ( - 2.0D+00 * ( 1.0D+00 - alpha ) - alpha * beta ) / alpha mbasis(2,3) = ( 2.0D+00 * ( 1.0D+00 - alpha ) & - beta * ( 1.0D+00 - 2.0D+00 * alpha ) ) / ( 1.0D+00 - beta ) mbasis(2,4) = - beta * beta / ( 1.0D+00 - beta ) mbasis(3,1) = - ( 1.0D+00 - alpha ) * ( 1.0D+00 - alpha ) / alpha mbasis(3,2) = ( 1.0D+00 - 2.0D+00 * alpha ) / alpha mbasis(3,3) = alpha mbasis(3,4) = 0.0D+00 mbasis(4,1) = 0.0D+00 mbasis(4,2) = 1.0D+00 mbasis(4,3) = 0.0D+00 mbasis(4,4) = 0.0D+00 return end subroutine basis_matrix_overhauser_nul ( alpha, mbasis ) !*****************************************************************************80 ! !! BASIS_MATRIX_OVERHAUSER_NUL sets the nonuniform left Overhauser basis matrix. ! ! Discussion: ! ! This basis matrix assumes that the data points P1, P2, and ! P3 are not uniformly spaced in T, and that P1 corresponds to T = 0, ! and P2 to T = 1. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 27 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) ALPHA. ! ALPHA = || P2 - P1 || / ( || P3 - P2 || + || P2 - P1 || ) ! ! Output, real ( kind = 8 ) MBASIS(3,3), the basis matrix. ! implicit none real ( kind = 8 ) alpha real ( kind = 8 ) mbasis(3,3) mbasis(1,1) = 1.0D+00 / alpha mbasis(1,2) = - 1.0D+00 / ( alpha * ( 1.0D+00 - alpha ) ) mbasis(1,3) = 1.0D+00 / ( 1.0D+00 - alpha ) mbasis(2,1) = - ( 1.0D+00 + alpha ) / alpha mbasis(2,2) = 1.0D+00 / ( alpha * ( 1.0D+00 - alpha ) ) mbasis(2,3) = - alpha / ( 1.0D+00 - alpha ) mbasis(3,1) = 1.0D+00 mbasis(3,2) = 0.0D+00 mbasis(3,3) = 0.0D+00 return end subroutine basis_matrix_overhauser_nur ( beta, mbasis ) !*****************************************************************************80 ! !! BASIS_MATRIX_OVERHAUSER_NUR: the nonuniform right Overhauser basis matrix. ! ! Discussion: ! ! This basis matrix assumes that the data points PN-2, PN-1, and ! PN are not uniformly spaced in T, and that PN-1 corresponds to T = 0, ! and PN to T = 1. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 27 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) BETA. ! BETA = || P(N) - P(N-1) || ! / ( || P(N) - P(N-1) || + || P(N-1) - P(N-2) || ) ! ! Output, real ( kind = 8 ) MBASIS(3,3), the basis matrix. ! implicit none real ( kind = 8 ) beta real ( kind = 8 ) mbasis(3,3) mbasis(1,1) = 1.0D+00 / beta mbasis(1,2) = - 1.0D+00 / ( beta * ( 1.0D+00 - beta ) ) mbasis(1,3) = 1.0D+00 / ( 1.0D+00 - beta ) mbasis(2,1) = - ( 1.0D+00 + beta ) / beta mbasis(2,2) = 1.0D+00 / ( beta * ( 1.0D+00 - beta ) ) mbasis(2,3) = - beta / ( 1.0D+00 - beta ) mbasis(3,1) = 1.0D+00 mbasis(3,2) = 0.0D+00 mbasis(3,3) = 0.0D+00 return end subroutine basis_matrix_overhauser_uni ( mbasis ) !*****************************************************************************80 ! !! BASIS_MATRIX_OVERHAUSER_UNI sets the uniform Overhauser spline basis matrix. ! ! Discussion: ! ! This basis matrix assumes that the data points P1, P2, P3 and ! P4 are uniformly spaced in T, and that P2 corresponds to T = 0, ! and P3 to T = 1. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 05 April 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! James Foley, Andries vanDam, Steven Feiner, John Hughes, ! Computer Graphics, Principles and Practice, ! Second Edition, ! Addison Wesley, 1995, ! ISBN: 0201848406, ! LC: T385.C5735. ! ! Parameters: ! ! Output, real ( kind = 8 ) MBASIS(4,4), the basis matrix. ! implicit none real ( kind = 8 ) mbasis(4,4) mbasis(1,1) = - 1.0D+00 / 2.0D+00 mbasis(1,2) = 3.0D+00 / 2.0D+00 mbasis(1,3) = - 3.0D+00 / 2.0D+00 mbasis(1,4) = 1.0D+00 / 2.0D+00 mbasis(2,1) = 2.0D+00 / 2.0D+00 mbasis(2,2) = - 5.0D+00 / 2.0D+00 mbasis(2,3) = 4.0D+00 / 2.0D+00 mbasis(2,4) = - 1.0D+00 / 2.0D+00 mbasis(3,1) = - 1.0D+00 / 2.0D+00 mbasis(3,2) = 0.0D+00 mbasis(3,3) = 1.0D+00 / 2.0D+00 mbasis(3,4) = 0.0D+00 mbasis(4,1) = 0.0D+00 mbasis(4,2) = 2.0D+00 / 2.0D+00 mbasis(4,3) = 0.0D+00 mbasis(4,4) = 0.0D+00 return end subroutine basis_matrix_overhauser_uni_l ( mbasis ) !*****************************************************************************80 ! !! BASIS_MATRIX_OVERHAUSER_UNI_L sets the left uniform Overhauser basis matrix. ! ! Discussion: ! ! This basis matrix assumes that the data points P1, P2, and P3 ! are not uniformly spaced in T, and that P1 corresponds to T = 0, ! and P2 to T = 1. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 05 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) MBASIS(3,3), the basis matrix. ! implicit none real ( kind = 8 ) mbasis(3,3) mbasis(1,1) = 2.0D+00 mbasis(1,2) = - 4.0D+00 mbasis(1,3) = 2.0D+00 mbasis(2,1) = - 3.0D+00 mbasis(2,2) = 4.0D+00 mbasis(2,3) = - 1.0D+00 mbasis(3,1) = 1.0D+00 mbasis(3,2) = 0.0D+00 mbasis(3,3) = 0.0D+00 return end subroutine basis_matrix_overhauser_uni_r ( mbasis ) !*****************************************************************************80 ! !! BASIS_MATRIX_OVERHAUSER_UNI_R sets the right uniform Overhauser basis matrix. ! ! Discussion: ! ! This basis matrix assumes that the data points P(N-2), P(N-1), ! and P(N) are uniformly spaced in T, and that P(N-1) corresponds to ! T = 0, and P(N) to T = 1. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 05 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) MBASIS(3,3), the basis matrix. ! implicit none real ( kind = 8 ) mbasis(3,3) mbasis(1,1) = 2.0D+00 mbasis(1,2) = - 4.0D+00 mbasis(1,3) = 2.0D+00 mbasis(2,1) = - 3.0D+00 mbasis(2,2) = 4.0D+00 mbasis(2,3) = - 1.0D+00 mbasis(3,1) = 1.0D+00 mbasis(3,2) = 0.0D+00 mbasis(3,3) = 0.0D+00 return end subroutine basis_matrix_tmp ( left, n, mbasis, ndata, tdata, ydata, tval, yval ) !*****************************************************************************80 ! !! BASIS_MATRIX_TMP computes Q = T * MBASIS * P ! ! Discussion: ! ! YDATA is a vector of data values, most frequently the values of some ! function sampled at uniformly spaced points. MBASIS is the basis ! matrix for a particular kind of spline. T is a vector of the ! powers of the normalized difference between TVAL and the left ! endpoint of the interval. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 10 February 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) LEFT, indicats that TVAL is in the interval ! [ TDATA(LEFT), TDATA(LEFT+1) ], or that this is the "nearest" ! interval to TVAL. ! For TVAL < TDATA(1), use LEFT = 1. ! For TDATA(NDATA) < TVAL, use LEFT = NDATA - 1. ! ! Input, integer ( kind = 4 ) N, the order of the basis matrix. ! ! Input, real ( kind = 8 ) MBASIS(N,N), the basis matrix. ! ! Input, integer ( kind = 4 ) NDATA, the dimension of the vectors TDATA ! and YDATA. ! ! Input, real ( kind = 8 ) TDATA(NDATA), the abscissa values. This routine ! assumes that the TDATA values are uniformly spaced, with an ! increment of 1.0. ! ! Input, real ( kind = 8 ) YDATA(NDATA), the data values to be ! interpolated or approximated. ! ! Input, real ( kind = 8 ) TVAL, the value of T at which the spline is to be ! evaluated. ! ! Output, real ( kind = 8 ) YVAL, the value of the spline at TVAL. ! implicit none integer ( kind = 4 ), parameter :: maxn = 4 integer ( kind = 4 ) n integer ( kind = 4 ) ndata real ( kind = 8 ) arg integer ( kind = 4 ) first integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) left real ( kind = 8 ) mbasis(n,n) real ( kind = 8 ) tdata(ndata) real ( kind = 8 ) tval real ( kind = 8 ) tvec(maxn) real ( kind = 8 ) ydata(ndata) real ( kind = 8 ) yval if ( left == 1 ) then arg = 0.5D+00 * ( tval - tdata(left) ) first = left else if ( left < ndata - 1 ) then arg = tval - tdata(left) first = left - 1 else if ( left == ndata - 1 ) then arg = 0.5D+00 * ( 1.0D+00 + tval - tdata(left) ) first = left - 1 end if ! ! TVEC(I) = ARG^(N-I). ! tvec(n) = 1.0D+00 do i = n-1, 1, -1 tvec(i) = arg * tvec(i+1) end do yval = 0.0D+00 do j = 1, n yval = yval + dot_product ( tvec(1:n), mbasis(1:n,j) ) & * ydata(first - 1 + j) end do return end subroutine bc_val ( n, t, xcon, ycon, xval, yval ) !*****************************************************************************80 ! !! BC_VAL evaluates a parameterized N-th degree Bezier curve in 2D. ! ! Discussion: ! ! BC_VAL(T) is the value of a vector function of the form ! ! BC_VAL(T) = ( X(T), Y(T) ) ! ! where ! ! X(T) = sum ( 0 <= I <= N ) XCON(I) * BERN(I,N)(T) ! Y(T) = sum ( 0 <= I <= N ) YCON(I) * BERN(I,N)(T) ! ! BERN(I,N)(T) is the I-th Bernstein polynomial of order N ! defined on the interval [0,1], ! ! XCON(0:N) and YCON(0:N) are the coordinates of N+1 "control points". ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 12 February 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the degree of the Bezier curve. ! N must be at least 0. ! ! Input, real ( kind = 8 ) T, the point at which the Bezier curve should ! be evaluated. The best results are obtained within the interval ! [0,1] but T may be anywhere. ! ! Input, real ( kind = 8 ) XCON(0:N), YCON(0:N), the X and Y coordinates ! of the control points. The Bezier curve will pass through ! the points ( XCON(0), YCON(0) ) and ( XCON(N), YCON(N) ), but ! generally NOT through the other control points. ! ! Output, real ( kind = 8 ) XVAL, YVAL, the X and Y coordinates of the point ! on the Bezier curve corresponding to the given T value. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) bval(0:n) real ( kind = 8 ) t real ( kind = 8 ) xcon(0:n) real ( kind = 8 ) xval real ( kind = 8 ) ycon(0:n) real ( kind = 8 ) yval call bp01 ( n, t, bval ) xval = dot_product ( xcon(0:n), bval(0:n) ) yval = dot_product ( ycon(0:n), bval(0:n) ) return end function bez_val ( n, x, a, b, y ) !*****************************************************************************80 ! !! BEZ_VAL evaluates an N-th degree Bezier function at a point. ! ! Discussion: ! ! The Bezier function has the form: ! ! BEZ(X) = sum ( 0 <= I <= N ) Y(I) * BERN(N,I)( (X-A)/(B-A) ) ! ! BERN(N,I)(X) is the I-th Bernstein polynomial of order N ! defined on the interval [0,1], ! ! Y(0:N) is a set of coefficients, ! ! and if, for I = 0 to N, we define the N+1 points ! ! X(I) = ( (N-I)*A + I*B) / N, ! ! equally spaced in [A,B], the pairs ( X(I), Y(I) ) can be regarded as ! "control points". ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 12 February 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the degree of the Bezier function. ! N must be at least 0. ! ! Input, real ( kind = 8 ) X, the point at which the Bezier function should ! be evaluated. The best results are obtained within the interval ! [A,B] but X may be anywhere. ! ! Input, real ( kind = 8 ) A, B, the interval over which the Bezier function ! has been defined. This is the interval in which the control ! points have been set up. Note BEZ(A) = Y(0) and BEZ(B) = Y(N), ! although BEZ will not, in general pass through the other ! control points. A and B must not be equal. ! ! Input, real ( kind = 8 ) Y(0:N), a set of data defining the Y coordinates ! of the control points. ! ! Output, real ( kind = 8 ) BEZ_VAL, the value of the Bezier function at X. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a real ( kind = 8 ) b real ( kind = 8 ) bez_val real ( kind = 8 ) bval(0:n) real ( kind = 8 ) x real ( kind = 8 ) x01 real ( kind = 8 ) y(0:n) if ( b - a == 0.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BEZ_VAL - Fatal error!' write ( *, '(a,g14.6)' ) ' Null interval, A = B = ', a stop 1 end if ! ! X01 lies in [0,1], in the same relative position as X in [A,B]. ! x01 = ( x - a ) / ( b - a ) call bp01 ( n, x01, bval ) bez_val = dot_product ( y(0:n), bval(0:n) ) return end subroutine bp01 ( n, x, bern ) !*****************************************************************************80 ! !! BP01 evaluates the Bernstein basis polynomials for [0,1] at a point. ! ! Discussion: ! ! For any N greater than or equal to 0, there is a set of N+1 Bernstein ! basis polynomials, each of degree N, which form a basis for ! all polynomials of degree N on [0,1]. ! ! BERN(N,I,X) = [N!/(I!*(N-I)!)] * (1-X)^(N-I) * X^I ! ! N is the degree; ! ! 0 <= I <= N indicates which of the N+1 basis polynomials ! of degree N to choose; ! ! X is a point in [0,1] at which to evaluate the basis polynomial. ! ! First values: ! ! B(0,0,X) = 1 ! ! B(1,0,X) = 1-X ! B(1,1,X) = X ! ! B(2,0,X) = (1-X)^2 ! B(2,1,X) = 2 * (1-X) * X ! B(2,2,X) = X^2 ! ! B(3,0,X) = (1-X)^3 ! B(3,1,X) = 3 * (1-X)^2 * X ! B(3,2,X) = 3 * (1-X) * X^2 ! B(3,3,X) = X^3 ! ! B(4,0,X) = (1-X)^4 ! B(4,1,X) = 4 * (1-X)^3 * X ! B(4,2,X) = 6 * (1-X)^2 * X^2 ! B(4,3,X) = 4 * (1-X) * X^3 ! B(4,4,X) = X^4 ! ! Special values: ! ! B(N,I,1/2) = C(N,K) / 2^N ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 12 February 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the degree of the Bernstein basis ! polynomials. N must be at least 0. ! ! Input, real ( kind = 8 ) X, the evaluation point. ! ! Output, real ( kind = 8 ) BERN(0:N), the values of the N+1 Bernstein basis ! polynomials at X. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) bern(0:n) integer ( kind = 4 ) i integer ( kind = 4 ) j real ( kind = 8 ) x if ( n == 0 ) then bern(0) = 1.0D+00 else if ( 0 < n ) then bern(0) = 1.0D+00 - x bern(1) = x do i = 2, n bern(i) = x * bern(i-1) do j = i-1, 1, -1 bern(j) = x * bern(j-1) + ( 1.0D+00 - x ) * bern(j) end do bern(0) = ( 1.0D+00 - x ) * bern(0) end do end if return end subroutine bpab ( n, a, b, x, bern ) !*****************************************************************************80 ! !! BPAB evaluates the Bernstein basis polynomials for [A,B] at a point. ! ! Discussion: ! ! BERN(N,I,X) = [N!/(I!*(N-I)!)] * (B-X)^(N-I) * (X-A)^I / (B-A)^N ! ! B(0,0,X) = 1 ! ! B(1,0,X) = ( B-X ) / (B-A) ! B(1,1,X) = ( X-A ) / (B-A) ! ! B(2,0,X) = ( (B-X)^2 ) / (B-A)^2 ! B(2,1,X) = ( 2 * (B-X * (X-A) ) / (B-A)^2 ! B(2,2,X) = ( (X-A)^2 ) / (B-A)^2 ! ! B(3,0,X) = ( (B-X)^3 ) / (B-A)^3 ! B(3,1,X) = ( 3 * (B-X)^2 * (X-A) ) / (B-A)^3 ! B(3,2,X) = ( 3 * (B-X) * (X-A)^2 ) / (B-A)^3 ! B(3,3,X) = ( (X-A)^3 ) / (B-A)^3 ! ! B(4,0,X) = ( (B-X)^4 ) / (B-A)^4 ! B(4,1,X) = ( 4 * (B-X)^3 * (X-A) ) / (B-A)^4 ! B(4,2,X) = ( 6 * (B-X)^2 * (X-A)^2 ) / (B-A)^4 ! B(4,3,X) = ( 4 * (B-X) * (X-A)^3 ) / (B-A)^4 ! B(4,4,X) = ( (X-A)^4 ) / (B-A)^4 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 12 February 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the degree of the Bernstein basis ! polynomials. There is a set of N+1 Bernstein basis polynomials, each ! of degree N, which form a basis for polynomials of degree N on [A,B]. ! N must be at least 0. ! ! Input, real ( kind = 8 ) A, B, the endpoints of the interval on which the ! polynomials are to be based. A and B should not be equal. ! Input, real ( kind = 8 ) X, the point at which the polynomials are to be ! evaluated. X need not lie in the interval [A,B]. ! ! Output, real ( kind = 8 ) BERN(0:N), the values of the N+1 Bernstein basis ! polynomials at X. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a real ( kind = 8 ) b real ( kind = 8 ) bern(0:n) integer ( kind = 4 ) i integer ( kind = 4 ) j real ( kind = 8 ) x if ( b == a ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BPAB - Fatal error!' write ( *, '(a,g14.6)' ) ' A = B = ', a stop 1 end if if ( n == 0 ) then bern(0) = 1.0D+00 else if ( 0 < n ) then bern(0) = ( b - x ) / ( b - a ) bern(1) = ( x - a ) / ( b - a ) do i = 2, n bern(i) = ( x - a ) * bern(i-1) / ( b - a ) do j = i-1, 1, -1 bern(j) = ( ( b - x ) * bern(j) + ( x - a ) * bern(j-1) ) / ( b - a ) end do bern(0) = ( b - x ) * bern(0) / ( b - a ) end do end if return end subroutine bpab_approx ( n, a, b, ydata, xval, yval ) !*****************************************************************************80 ! !! BPAB_APPROX evaluates the Bernstein polynomial approximant to F(X) on [A,B]. ! ! Formula: ! ! BERN(F)(X) = sum ( 0 <= I <= N ) F(X(I)) * B_BASE(I,X) ! ! where ! ! X(I) = ( ( N - I ) * A + I * B ) / N ! B_BASE(I,X) is the value of the I-th Bernstein basis polynomial at X. ! ! Discussion: ! ! The Bernstein polynomial BERN(F) for F(X) is an approximant, not an ! interpolant; in other words, its value is not guaranteed to equal ! that of F at any particular point. However, for a fixed interval ! [A,B], if we let N increase, the Bernstein polynomial converges ! uniformly to F everywhere in [A,B], provided only that F is continuous. ! Even if F is not continuous, but is bounded, the polynomial converges ! pointwise to F(X) at all points of continuity. On the other hand, ! the convergence is quite slow compared to other interpolation ! and approximation schemes. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 12 February 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the degree of the Bernstein polynomial ! to be used. N must be at least 0. ! ! Input, real ( kind = 8 ) A, B, the endpoints of the interval on which the ! approximant is based. A and B should not be equal. ! ! Input, real ( kind = 8 ) YDATA(0:N), the data values at N+1 equally ! spaced points in [A,B]. If N = 0, then the evaluation point should ! be 0.5 * ( A + B). Otherwise, evaluation point I should be ! ( (N-I)*A + I*B ) / N ). ! ! Input, real ( kind = 8 ) XVAL, the point at which the Bernstein polynomial ! approximant is to be evaluated. XVAL does not have to lie in the ! interval [A,B]. ! ! Output, real ( kind = 8 ) YVAL, the value of the Bernstein polynomial ! approximant for F, based in [A,B], evaluated at XVAL. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a real ( kind = 8 ) b real ( kind = 8 ) bvec(0:n) real ( kind = 8 ) xval real ( kind = 8 ) ydata(0:n) real ( kind = 8 ) yval ! ! Evaluate the Bernstein basis polynomials at XVAL. ! call bpab ( n, a, b, xval, bvec ) ! ! Now compute the sum of YDATA(I) * BVEC(I). ! yval = dot_product ( ydata(0:n), bvec(0:n) ) return end subroutine chfev ( x1, x2, f1, f2, d1, d2, ne, xe, fe, next, ierr ) !*****************************************************************************80 ! !! CHFEV evaluates a cubic polynomial given in Hermite form. ! ! Discussion: ! ! This routine evaluates a cubic polynomial given in Hermite form at an ! array of points. While designed for use by SPLINE_PCHIP_VAL, it may ! be useful directly as an evaluator for a piecewise cubic ! Hermite function in applications, such as graphing, where ! the interval is known in advance. ! ! The cubic polynomial is determined by function values ! F1, F2 and derivatives D1, D2 on the interval [X1,X2]. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 25 June 2008 ! ! Author: ! ! Original FORTRAN77 version by Fred Fritsch. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Fred Fritsch, Ralph Carlson, ! Monotone Piecewise Cubic Interpolation, ! SIAM Journal on Numerical Analysis, ! Volume 17, Number 2, April 1980, pages 238-246. ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Parameters: ! ! Input, real ( kind = 8 ) X1, X2, the endpoints of the interval of ! definition of the cubic. X1 and X2 must be distinct. ! ! Input, real ( kind = 8 ) F1, F2, the values of the function at X1 and ! X2, respectively. ! ! Input, real ( kind = 8 ) D1, D2, the derivative values at X1 and ! X2, respectively. ! ! Input, integer ( kind = 4 ) NE, the number of evaluation points. ! ! Input, real ( kind = 8 ) XE(NE), the points at which the function is to ! be evaluated. If any of the XE are outside the interval ! [X1,X2], a warning error is returned in NEXT. ! ! Output, real ( kind = 8 ) FE(NE), the value of the cubic function ! at the points XE. ! ! Output, integer ( kind = 4 ) NEXT(2), indicates the number of ! extrapolation points: ! NEXT(1) = number of evaluation points to the left of interval. ! NEXT(2) = number of evaluation points to the right of interval. ! ! Output, integer ( kind = 4 ) IERR, error flag. ! 0, no errors. ! -1, NE < 1. ! -2, X1 == X2. ! implicit none integer ( kind = 4 ) ne real ( kind = 8 ) c2 real ( kind = 8 ) c3 real ( kind = 8 ) d1 real ( kind = 8 ) d2 real ( kind = 8 ) del1 real ( kind = 8 ) del2 real ( kind = 8 ) delta real ( kind = 8 ) f1 real ( kind = 8 ) f2 real ( kind = 8 ) fe(ne) real ( kind = 8 ) h integer ( kind = 4 ) i integer ( kind = 4 ) ierr integer ( kind = 4 ) next(2) real ( kind = 8 ) x real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) xe(ne) real ( kind = 8 ) xma real ( kind = 8 ) xmi if ( ne < 1 ) then ierr = -1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CHFEV - Fatal error!' write ( *, '(a)' ) ' Number of evaluation points is less than 1.' write ( *, '(a,i8)' ) ' NE = ', ne stop 1 end if h = x2 - x1 if ( h == 0.0D+00 ) then ierr = -2 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CHFEV - Fatal error!' write ( *, '(a)' ) ' The interval [X1,X2] is of zero length.' stop 1 end if ! ! Initialize. ! ierr = 0 next(1) = 0 next(2) = 0 xmi = min ( 0.0D+00, h ) xma = max ( 0.0D+00, h ) ! ! Compute cubic coefficients expanded about X1. ! delta = ( f2 - f1 ) / h del1 = ( d1 - delta ) / h del2 = ( d2 - delta ) / h c2 = -( del1 + del1 + del2 ) c3 = ( del1 + del2 ) / h ! ! Evaluation loop. ! do i = 1, ne x = xe(i) - x1 fe(i) = f1 + x * ( d1 + x * ( c2 + x * c3 ) ) ! ! Count the extrapolation points. ! if ( x < xmi ) then next(1) = next(1) + 1 end if if ( xma < x ) then next(2) = next(2) + 1 end if end do return end subroutine data_to_dif ( ntab, xtab, ytab, diftab ) !*****************************************************************************80 ! !! DATA_TO_DIF sets up a divided difference table from raw data. ! ! Discussion: ! ! Space can be saved by using a single array for both the DIFTAB and ! YTAB dummy parameters. In that case, the divided difference table will ! overwrite the Y data without interfering with the computation. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 04 September 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Carl deBoor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663. ! ! Parameters: ! ! Input, integer ( kind = 4 ) NTAB, the number of pairs of points ! (XTAB(I),YTAB(I)) which are to be used as data. The ! number of entries to be used in DIFTAB, XTAB and YTAB. ! ! Input, real ( kind = 8 ) XTAB(NTAB), the X values at which data was taken. ! These values must be distinct. ! ! Input, real ( kind = 8 ) YTAB(NTAB), the corresponding Y values. ! ! Output, real ( kind = 8 ) DIFTAB(NTAB), the divided difference coefficients ! corresponding to the input (XTAB,YTAB). ! implicit none integer ( kind = 4 ) ntab real ( kind = 8 ) diftab(ntab) integer ( kind = 4 ) i integer ( kind = 4 ) j ! logical r8vec_distinct real ( kind = 8 ) xtab(ntab) real ( kind = 8 ) ytab(ntab) if ( .not. r8vec_distinct ( ntab, xtab ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DATA_TO_DIF - Fatal error!' write ( *, '(a)' ) ' Two entries of XTAB are equal!' return end if ! ! Copy the data values into DIFTAB. ! diftab(1:ntab) = ytab(1:ntab) ! ! Compute the divided differences. ! do i = 2, ntab do j = ntab, i, -1 diftab(j) = ( diftab(j) - diftab(j-1) ) / ( xtab(j) - xtab(j+1-i) ) end do end do return end subroutine dif_val ( ntab, xtab, diftab, xval, yval ) !*****************************************************************************80 ! !! DIF_VAL evaluates a divided difference polynomial at a point. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 05 September 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Carl deBoor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663. ! ! Parameters: ! ! Input, integer ( kind = 4 ) NTAB, the number of divided difference ! coefficients in DIFTAB, and the number of points XTAB. ! ! Input, real ( kind = 8 ) XTAB(NTAB), the X values upon which the ! divided difference polynomial is based. ! ! Input, real ( kind = 8 ) DIFTAB(NTAB), the divided difference ! polynomial coefficients. ! ! Input, real ( kind = 8 ) XVAL, the value where the polynomial ! is to be evaluated. ! ! Output, real ( kind = 8 ) YVAL, the value of the polynomial at XVAL. ! implicit none integer ( kind = 4 ) ntab real ( kind = 8 ) diftab(ntab) integer ( kind = 4 ) i real ( kind = 8 ) xtab(ntab) real ( kind = 8 ) xval real ( kind = 8 ) yval yval = diftab(ntab) do i = 1, ntab-1 yval = diftab(ntab-i) + ( xval - xtab(ntab-i) ) * yval end do return end subroutine least_set_old ( ntab, xtab, ytab, ndeg, ptab, b, c, d, eps, ierror ) !*****************************************************************************80 ! !! LEAST_SET_OLD constructs the least squares polynomial approximation to data. ! ! Discussion: ! ! The least squares polynomial is not returned directly as a simple ! polynomial. Instead, it is represented in terms of a set of ! orthogonal polynomials appopriate for the given data. This makes ! the computation more accurate, but means that the user can not ! easily evaluate the computed polynomial. Instead, the routine ! LEAST_EVAL should be used to evaluate the least squares polynomial ! at any point. (However, the value of the least squares polynomial ! at each of the data points is returned as part of this computation.) ! ! ! A discrete unweighted inner product is used, so that ! ! ( F(X), G(X) ) = sum ( 1 <= I <= NTAB ) F(XTAB(I)) * G(XTAB(I)). ! ! The least squares polynomial is determined using a set of ! orthogonal polynomials PHI. These polynomials can be defined ! recursively by: ! ! PHI(0)(X) = 1 ! PHI(1)(X) = X - B(1) ! PHI(I)(X) = ( X - B(I) ) * PHI(I-1)(X) - D(I) * PHI(I-2)(X) ! ! The array B(1:NDEG) contains the values ! ! B(I) = ( X*PHI(I-1), PHI(I-1) ) / ( PHI(I-1), PHI(I-1) ) ! ! The array D(2:NDEG) contains the values ! ! D(I) = ( PHI(I-1), PHI(I-1) ) / ( PHI(I-2), PHI(I-2) ) ! ! Using this basis, the least squares polynomial can be represented as ! ! P(X)(I) = sum ( 0 <= I <= NDEG ) C(I) * PHI(I)(X) ! ! The array C(0:NDEG) contains the values ! ! C(I) = ( YTAB(I), PHI(I) ) / ( PHI(I), PHI(I) ) ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 15 May 2004 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Gisela Engeln-Muellges, Frank Uhlig, ! Numerical Algorithms with C, ! Springer, 1996, ! ISBN: 3-540-60530-4. ! ! Parameters: ! ! Input, integer ( kind = 4 ) NTAB, the number of data points. ! ! Input, real ( kind = 8 ) XTAB(NTAB), the X data. The values in XTAB ! should be distinct, and in increasing order. ! ! Input, real ( kind = 8 ) YTAB(NTAB), the Y data values corresponding ! to the X data in XTAB. ! ! Input, integer ( kind = 4 ) NDEG, the degree of the polynomial which the ! program is to use. NDEG must be at least 0, and less than or ! equal to NTAB-1. ! ! Output, real ( kind = 8 ) PTAB(NTAB), the value of the least ! squares polynomial at the points XTAB(1:NTAB). ! ! Output, real ( kind = 8 ) B(1:NDEG), C(0:NDEG), D(2:NDEG), arrays ! needed to evaluate the polynomial. ! ! Output, real ( kind = 8 ) EPS, the root-mean-square discrepancy of the ! polynomial fit. ! ! Output, integer ( kind = 4 ) IERROR, error flag. ! zero, no error occurred; ! nonzero, an error occurred, and the polynomial could not be computed. ! implicit none integer ( kind = 4 ) ndeg integer ( kind = 4 ) ntab real ( kind = 8 ) b(1:ndeg) real ( kind = 8 ) c(0:ndeg) real ( kind = 8 ) d(2:ndeg) real ( kind = 8 ) eps integer ( kind = 4 ) i integer ( kind = 4 ) i0l1 integer ( kind = 4 ) i1l1 integer ( kind = 4 ) ierror integer ( kind = 4 ) it integer ( kind = 4 ) k integer ( kind = 4 ) mdeg real ( kind = 8 ) ptab(ntab) real ( kind = 8 ) rn0 real ( kind = 8 ) rn1 real ( kind = 8 ) s real ( kind = 8 ) sum2 real ( kind = 8 ) xtab(ntab) real ( kind = 8 ) y_sum real ( kind = 8 ) ytab(ntab) real ( kind = 8 ) ztab(2*ntab) ierror = 0 ! ! Check NDEG. ! if ( ndeg < 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LEAST_SET_OLD - Fatal error!' write ( *, '(a)' ) ' NDEG < 0.' stop 1 end if if ( ntab <= ndeg ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LEAST_SET_OLD - Fatal error!' write ( *, '(a)' ) ' NTAB <= NDEG.' stop 1 end if ! ! Check that the abscissas are strictly increasing. ! do i = 1, ntab-1 if ( xtab(i+1) <= xtab(i) ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LEAST_SET_OLD - Fatal error!' write ( *, '(a)' ) ' XTAB must be strictly increasing, but' write ( *, '(a,i8,a,g14.6)' ) ' XTAB(', i, ') = ', xtab(i) write ( *, '(a,i8,a,g14.6)' ) ' XTAB(', i+1, ') = ', xtab(i+1) stop 1 end if end do i0l1 = 0 i1l1 = ntab ! ! The polynomial is of degree at least 0. ! y_sum = sum ( ytab(1:ntab) ) rn0 = ntab c(0) = y_sum / real ( ntab, kind = 8 ) ptab(1:ntab) = y_sum / real ( ntab, kind = 8 ) if ( ndeg == 0 ) then eps = sum ( ( ptab(1:ntab) - ytab(1:ntab) )**2 ) eps = sqrt ( eps / real ( ntab, kind = 8 ) ) return end if ! ! The polynomial is of degree at least 1. ! b(1) = sum ( xtab(1:ntab) ) / real ( ntab, kind = 8 ) s = 0.0D+00 sum2 = 0.0D+00 do i = 1, ntab ztab(i1l1+i) = xtab(i) - b(1) s = s + ztab(i1l1+i)**2 sum2 = sum2 + ztab(i1l1+i) * ( ytab(i) - ptab(i) ) end do rn1 = s c(1) = sum2 / s do i = 1, ntab ptab(i) = ptab(i) + c(1) * ztab(i1l1+i) end do if ( ndeg == 1 ) then eps = sum ( ( ptab(1:ntab) - ytab(1:ntab) )**2 ) eps = sqrt ( eps / real ( ntab, kind = 8 ) ) return end if ztab(1:ntab) = 1.0D+00 mdeg = 2 k = 2 do d(k) = rn1 / rn0 sum2 = 0.0D+00 do i = 1, ntab sum2 = sum2 + xtab(i) * ztab(i1l1+i) * ztab(i1l1+i) end do b(k) = sum2 / rn1 s = 0.0D+00 sum2 = 0.0D+00 do i = 1, ntab ztab(i0l1+i) = ( xtab(i) - b(k) ) * ztab(i1l1+i) & - d(k) * ztab(i0l1+i) s = s + ztab(i0l1+i) * ztab(i0l1+i) sum2 = sum2 + ztab(i0l1+i) * ( ytab(i) - ptab(i) ) end do rn0 = rn1 rn1 = s c(k) = sum2 / rn1 it = i0l1 i0l1 = i1l1 i1l1 = it do i = 1, ntab ptab(i) = ptab(i) + c(k) * ztab(i1l1+i) end do if ( ndeg <= mdeg ) then exit end if mdeg = mdeg + 1 k = k + 1 end do ! ! Compute the RMS error. ! eps = sum ( ( ptab(1:ntab) - ytab(1:ntab) )**2 ) eps = sqrt ( eps / real ( ntab, kind = 8 ) ) return end subroutine least_val_old ( x, ndeg, b, c, d, value ) !*****************************************************************************80 ! !! LEAST_VAL_OLD evaluates a least squares polynomial defined by LEAST_SET_OLD. ! ! Discussion: ! ! This is an "old" version of the routine. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 15 May 2004 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Gisela Engeln-Muellges, Frank Uhlig, ! Numerical Algorithms with C, ! Springer, 1996, ! ISBN: 3-540-60530-4. ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the point at which the polynomial is ! to be evaluated. ! ! Input, integer ( kind = 4 ) NDEG, the degree of the least squares ! polynomial. ! ! Input, real ( kind = 8 ) B(1:NDEG), C(0:NDEG), D(2:NDEG), arrays ! defined by LEAST_SET_OLD, and needed to evaluate the polynomial. ! ! Output, real ( kind = 8 ) VALUE, the value of the polynomial at X. ! implicit none integer ( kind = 4 ) ndeg real ( kind = 8 ) b(1:ndeg) real ( kind = 8 ) c(0:ndeg) real ( kind = 8 ) d(2:ndeg) integer ( kind = 4 ) k real ( kind = 8 ) sk real ( kind = 8 ) skp1 real ( kind = 8 ) skp2 real ( kind = 8 ) value real ( kind = 8 ) x if ( ndeg <= 0 ) then value = c(0) else if ( ndeg == 1 ) then value = c(0) + c(1) * ( x - b(1) ) else skp2 = c(ndeg) skp1 = c(ndeg-1) + ( x - b(ndeg) ) * skp2 do k = ndeg-2, 0, -1 sk = c(k) + ( x - b(k+1) ) * skp1 - d(k+2) * skp2 skp2 = skp1 skp1 = sk end do value = sk end if return end subroutine least_set ( point_num, x, f, w, nterms, b, c, d ) !*****************************************************************************80 ! !! LEAST_SET defines a least squares polynomial for given data. ! ! Discussion: ! ! This routine is based on ORTPOL by Conte and deBoor. ! ! The polynomial may be evaluated at any point X by calling LEAST_VAL. ! ! Thanks to Andrew Telford for pointing out a mistake in the form of ! the check that there are enough unique data points, 25 June 2008. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 25 June 2008 ! ! Author: ! ! Original FORTRAN77 version by Samuel Conte, Carl deBoor. ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Samuel Conte, Carl deBoor, ! Elementary Numerical Analysis, ! Second Edition, ! McGraw Hill, 1972, ! ISBN: 07-012446-4, ! LC: QA297.C65. ! ! Parameters: ! ! Input, integer ( kind = 4 ) POINT_NUM, the number of data values. ! ! Input, real ( kind = 8 ) X(POINT_NUM), the abscissas of the data points. ! At least NTERMS of the values in X must be distinct. ! ! Input, real ( kind = 8 ) F(POINT_NUM), the data values at the points X(*). ! ! Input, real ( kind = 8 ) W(POINT_NUM), the weights associated with ! the data points. Each entry of W should be positive. ! ! Input, integer ( kind = 4 ) NTERMS, the number of terms to use in the ! approximating polynomial. NTERMS must be at least 1. ! The degree of the polynomial is NTERMS-1. ! ! Output, real ( kind = 8 ) B(NTERMS), C(NTERMS), D(NTERMS), are quantities ! defining the least squares polynomial for the input data, ! which will be needed to evaluate the polynomial. ! implicit none integer ( kind = 4 ) point_num integer ( kind = 4 ) nterms real ( kind = 8 ) b(nterms) real ( kind = 8 ) c(nterms) real ( kind = 8 ) d(nterms) real ( kind = 8 ) f(point_num) integer ( kind = 4 ) i integer ( kind = 4 ) j real ( kind = 8 ) p real ( kind = 8 ) pj(point_num) real ( kind = 8 ) pjm1(point_num) real ( kind = 8 ) s(nterms) real ( kind = 8 ), parameter :: tol = 0.0D+00 integer ( kind = 4 ) unique_num real ( kind = 8 ) w(point_num) real ( kind = 8 ) x(point_num) ! ! Make sure at least NTERMS X values are unique. ! call r8vec_unique_count ( point_num, x, tol, unique_num ) if ( unique_num < nterms ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LEAST_SET - Fatal error!' write ( *, '(a)' ) ' The number of distinct X values must be' write ( *, '(a,i8)') ' at least NTERMS = ', nterms write ( *, '(a,i8)' ) ' but the input data has only ', unique_num write ( *, '(a)' ) ' distinct entries.' return end if ! ! Make sure all W entries are positive. ! do i = 1, point_num if ( w(i) <= 0.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LEAST_SET - Fatal error!' write ( *, '(a)' ) ' All weights W must be positive,' write ( *, '(a,i8)' ) ' but weight ', i write ( *, '(a,g14.6)' ) ' is ', w(i) return end if end do ! ! Start inner product summations at zero. ! b(1:nterms) = 0.0D+00 c(1:nterms) = 0.0D+00 d(1:nterms) = 0.0D+00 s(1:nterms) = 0.0D+00 ! ! Set the values of P(-1,X) and P(0,X) at all data points. ! pjm1(1:point_num) = 0.0D+00 pj(1:point_num) = 1.0D+00 ! ! Now compute the value of P(J,X(I)) as ! ! P(J,X(I)) = ( X(I) - B(J) ) * P(J-1,X(I)) - C(J) * P(J-2,X(I)) ! ! where ! ! S(J) = < P(J,X), P(J,X) > ! B(J) = < x*P(J,X), P(J,X) > / < P(J,X), P(J,X) > ! C(J) = S(J) / S(J-1) ! ! and the least squares coefficients are ! ! D(J) = < F(X), P(J,X) > / < P(J,X), P(J,X) > ! do j = 1, nterms d(j) = d(j) + sum ( w(1:point_num) * f(1:point_num) * pj(1:point_num) ) b(j) = b(j) + sum ( w(1:point_num) * x(1:point_num) * pj(1:point_num)**2 ) s(j) = s(j) + sum ( w(1:point_num) * pj(1:point_num)**2 ) d(j) = d(j) / s(j) if ( j == nterms ) then c(j) = 0.0D+00 return end if b(j) = b(j) / s(j) if ( j == 1 ) then c(j) = 0.0D+00 else c(j) = s(j) / s(j-1) end if do i = 1, point_num p = pj(i) pj(i) = ( x(i) - b(j) ) * pj(i) - c(j) * pjm1(i) pjm1(i) = p end do end do return end subroutine least_val ( nterms, b, c, d, x, px ) !*****************************************************************************80 ! !! LEAST_VAL evaluates a least squares polynomial defined by LEAST_SET. ! ! Discussion: ! ! The least squares polynomial is assumed to be defined as a sum ! ! P(X) = sum ( 1 <= I <= NTERMS ) D(I) * P(I-1,X) ! ! where the orthogonal basis polynomials P(I,X) satisfy the following ! three term recurrence: ! ! P(-1,X) = 0 ! P(0,X) = 1 ! P(I,X) = ( X - B(I-1) ) * P(I-1,X) - C(I-1) * P(I-2,X) ! ! Therefore, the least squares polynomial can be evaluated as follows: ! ! If NTERMS is 1, then the value of P(X) is D(1) * P(0,X) = D(1). ! ! Otherwise, P(X) is defined as the sum of NTERMS > 1 terms. We can ! reduce the number of terms by 1, because the polynomial P(NTERMS,X) ! can be rewritten as a sum of polynomials; Therefore, P(NTERMS,X) ! can be eliminated from the sum, and its coefficient merged in with ! those of other polynomials. Repeat this process for P(NTERMS-1,X) ! and so on until a single term remains. ! P(NTERMS,X) of P(NTERMS-1,X) ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 24 May 2005 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Samuel Conte, Carl deBoor, ! Elementary Numerical Analysis, ! Second Edition, ! McGraw Hill, 1972, ! ISBN: 07-012446-4, ! LC: QA297.C65. ! ! Parameters: ! ! Input, integer ( kind = 4 ) NTERMS, the number of terms in the least ! squares polynomial. NTERMS must be at least 1. The input value of NTERMS ! may be reduced from the value given to LEAST_SET. This will ! evaluate the least squares polynomial of the lower degree specified. ! ! Input, real ( kind = 8 ) B(NTERMS), C(NTERMS), D(NTERMS), the information ! computed by LEAST_SET. ! ! Input, real ( kind = 8 ) X, the point at which the least squares polynomial ! is to be evaluated. ! ! Output, real ( kind = 8 ) PX, the value of the least squares ! polynomial at X. ! implicit none integer ( kind = 4 ) nterms real ( kind = 8 ) b(nterms) real ( kind = 8 ) c(nterms) real ( kind = 8 ) d(nterms) integer ( kind = 4 ) i real ( kind = 8 ) prev real ( kind = 8 ) prev2 real ( kind = 8 ) px real ( kind = 8 ) x px = d(nterms) prev = 0.0D+00 do i = nterms - 1, 1, -1 prev2 = prev prev = px if ( i == nterms-1 ) then px = d(i) + ( x - b(i) ) * prev else px = d(i) + ( x - b(i) ) * prev - c(i+1) * prev2 end if end do return end subroutine least_val2 ( nterms, b, c, d, x, px, pxp ) !*****************************************************************************80 ! !! LEAST_VAL2 evaluates a least squares polynomial defined by LEAST_SET. ! ! Discussion: ! ! This routine also computes the derivative of the polynomial. ! ! The least squares polynomial is assumed to be defined as a sum ! ! P(X) = sum ( 1 <= I <= NTERMS ) D(I) * P(I-1,X) ! ! where the orthogonal basis polynomials P(I,X) satisfy the following ! three term recurrence: ! ! P(-1,X) = 0 ! P(0,X) = 1 ! P(I,X) = ( X - B(I-1) ) * P(I-1,X) - C(I-1) * P(I-2,X) ! ! Therefore, the least squares polynomial can be evaluated as follows: ! ! If NTERMS is 1, then the value of P(X) is D(1) * P(0,X) = D(1). ! ! Otherwise, P(X) is defined as the sum of NTERMS > 1 terms. We can ! reduce the number of terms by 1, because the polynomial P(NTERMS,X) ! can be rewritten as a sum of polynomials; Therefore, P(NTERMS,X) ! can be eliminated from the sum, and its coefficient merged in with ! those of other polynomials. Repeat this process for P(NTERMS-1,X) ! and so on until a single term remains. ! P(NTERMS,X) of P(NTERMS-1,X) ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 24 May 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NTERMS, the number of terms in the least ! squares polynomial. NTERMS must be at least 1. The value of NTERMS ! may be reduced from the value given to LEAST_SET. ! This will cause LEAST_VAL to evaluate the least squares polynomial ! of the lower degree specified. ! ! Input, real ( kind = 8 ) B(NTERMS), C(NTERMS), D(NTERMS), the information ! computed by LEAST_SET. ! ! Input, real ( kind = 8 ) X, the point at which the least squares polynomial ! is to be evaluated. ! ! Output, real ( kind = 8 ) PX, PXP, the value and derivative of the least ! squares polynomial at X. ! implicit none integer ( kind = 4 ) nterms real ( kind = 8 ) b(nterms) real ( kind = 8 ) c(nterms) real ( kind = 8 ) d(nterms) integer ( kind = 4 ) i real ( kind = 8 ) px real ( kind = 8 ) pxm1 real ( kind = 8 ) pxm2 real ( kind = 8 ) pxp real ( kind = 8 ) pxpm1 real ( kind = 8 ) pxpm2 real ( kind = 8 ) x px = d(nterms) pxp = 0.0D+00 pxm1 = 0.0D+00 pxpm1 = 0.0D+00 do i = nterms - 1, 1, -1 pxm2 = pxm1 pxpm2 = pxpm1 pxm1 = px pxpm1 = pxp if ( i == nterms - 1 ) then px = d(i) + ( x - b(i) ) * pxm1 pxp = pxm1 + ( x - b(i) ) * pxpm1 else px = d(i) + ( x - b(i) ) * pxm1 - c(i+1) * pxm2 pxp = pxm1 + ( x - b(i) ) * pxpm1 - c(i+1) * pxpm2 end if end do return end subroutine parabola_val2 ( dim_num, ndata, tdata, ydata, left, tval, yval ) !*****************************************************************************80 ! !! PARABOLA_VAL2 evaluates a parabolic interpolant through tabular data. ! ! Discussion: ! ! This routine is a utility routine used by OVERHAUSER_SPLINE_VAL. ! It constructs the parabolic interpolant through the data in ! 3 consecutive entries of a table and evaluates this interpolant ! at a given abscissa value. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 26 January 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) DIM_NUM, the dimension of a single data point. ! DIM_NUM must be at least 1. ! ! Input, integer ( kind = 4 ) NDATA, the number of data points. ! NDATA must be at least 3. ! ! Input, real ( kind = 8 ) TDATA(NDATA), the abscissas of the data ! points. The values in TDATA must be in strictly ascending order. ! ! Input, real ( kind = 8 ) YDATA(DIM_NUM,NDATA), the data points ! corresponding to the abscissas. ! ! Input, integer ( kind = 4 ) LEFT, the location of the first of the three ! consecutive data points through which the parabolic interpolant ! must pass. 1 <= LEFT <= NDATA - 2. ! ! Input, real ( kind = 8 ) TVAL, the value of T at which the parabolic ! interpolant is to be evaluated. Normally, TDATA(1) <= TVAL <= T(NDATA), ! and the data will be interpolated. For TVAL outside this range, ! extrapolation will be used. ! ! Output, real ( kind = 8 ) YVAL(DIM_NUM), the value of the parabolic ! interpolant at TVAL. ! implicit none integer ( kind = 4 ) ndata integer ( kind = 4 ) dim_num real ( kind = 8 ) dif1 real ( kind = 8 ) dif2 integer ( kind = 4 ) i integer ( kind = 4 ) left real ( kind = 8 ) t1 real ( kind = 8 ) t2 real ( kind = 8 ) t3 real ( kind = 8 ) tval real ( kind = 8 ) tdata(ndata) real ( kind = 8 ) ydata(dim_num,ndata) real ( kind = 8 ) y1 real ( kind = 8 ) y2 real ( kind = 8 ) y3 real ( kind = 8 ) yval(dim_num) ! ! Check. ! if ( left < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PARABOLA_VAL2 - Fatal error!' write ( *, '(a)' ) ' LEFT < 1.' write ( *, '(a,i8)' ) ' LEFT = ', left stop 1 end if if ( ndata-2 < left ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PARABOLA_VAL2 - Fatal error!' write ( *, '(a)' ) ' NDATA-2 < LEFT.' write ( *, '(a,i8)' ) ' NDATA = ', ndata write ( *, '(a,i8)' ) ' LEFT = ', left stop 1 end if if ( dim_num < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PARABOLA_VAL2 - Fatal error!' write ( *, '(a)' ) ' DIM_NUM < 1.' stop 1 end if ! ! Copy out the three abscissas. ! t1 = tdata(left) t2 = tdata(left+1) t3 = tdata(left+2) if ( t2 <= t1 .or. t3 <= t2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PARABOLA_VAL2 - Fatal error!' write ( *, '(a)' ) ' T2 <= T1 or T3 <= T2.' stop 1 end if ! ! Construct and evaluate a parabolic interpolant for the data ! in each dimension. ! do i = 1, dim_num y1 = ydata(i,left) y2 = ydata(i,left+1) y3 = ydata(i,left+2) dif1 = ( y2 - y1 ) / ( t2 - t1 ) dif2 = ( ( y3 - y1 ) / ( t3 - t1 ) & - ( y2 - y1 ) / ( t2 - t1 ) ) / ( t3 - t2 ) yval(i) = y1 + ( tval - t1 ) * ( dif1 + ( tval - t2 ) * dif2 ) end do return end function pchst ( arg1, arg2 ) !*****************************************************************************80 ! !! PCHST: PCHIP sign-testing routine. ! ! Discussion: ! ! This routine essentially computes the sign of ARG1 * ARG2. ! ! The object is to do this without multiplying ARG1 * ARG2, to avoid ! possible over/underflow problems. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 25 June 2008 ! ! Author: ! ! Original FORTRAN77 version by Fred Fritsch. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Fred Fritsch, Ralph Carlson, ! Monotone Piecewise Cubic Interpolation, ! SIAM Journal on Numerical Analysis, ! Volume 17, Number 2, April 1980, pages 238-246. ! ! Parameters: ! ! Input, real ( kind = 8 ) ARG1, ARG2, two values to check. ! ! Output, real ( kind = 8 ) PCHST, ! -1.0, if ARG1 and ARG2 are of opposite sign. ! 0.0, if either argument is zero. ! +1.0, if ARG1 and ARG2 are of the same sign. ! implicit none real ( kind = 8 ) arg1 real ( kind = 8 ) arg2 real ( kind = 8 ) pchst pchst = sign ( 1.0D+00, arg1 ) * sign ( 1.0D+00, arg2 ) if ( arg1 == 0.0D+00 .or. arg2 == 0.0D+00 ) then pchst = 0.0D+00 end if return end subroutine penta ( n, a1, a2, a3, a4, a5, b, x ) !*****************************************************************************80 ! !! PENTA solves a pentadiagonal system of linear equations. ! ! Discussion: ! ! The matrix A is pentadiagonal. It is entirely zero, except for ! the main diagaonal, and the two immediate sub- and super-diagonals. ! ! The entries of Row I are stored as: ! ! A(I,I-2) -> A1(I) ! A(I,I-1) -> A2(I) ! A(I,I) -> A3(I) ! A(I,I+1) -> A4(I) ! A(I,I-2) -> A5(I) ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 07 June 2013 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Cheney, Kincaid, ! Numerical Mathematics and Computing, ! 1985, pages 233-236. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input, real ( kind = 8 ) A1(N), A2(N), A3(N), A4(N), A5(N), the nonzero ! elements of the matrix. Note that the data in A2, A3 and A4 ! is overwritten by this routine during the solution process. ! ! Input, real ( kind = 8 ) B(N), the right hand side of the linear system. ! ! Output, real ( kind = 8 ) X(N), the solution of the linear system. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a1(n) real ( kind = 8 ) a2(n) real ( kind = 8 ) a3(n) real ( kind = 8 ) a4(n) real ( kind = 8 ) a5(n) real ( kind = 8 ) b(n) integer ( kind = 4 ) i real ( kind = 8 ) x(n) real ( kind = 8 ) xmult do i = 2, n - 1 xmult = a2(i) / a3(i-1) a3(i) = a3(i) - xmult * a4(i-1) a4(i) = a4(i) - xmult * a5(i-1) b(i) = b(i) - xmult * b(i-1) xmult = a1(i+1) / a3(i-1) a2(i+1) = a2(i+1) - xmult * a4(i-1) a3(i+1) = a3(i+1) - xmult * a5(i-1) b(i+1) = b(i+1) - xmult * b(i-1) end do xmult = a2(n) / a3(n-1) a3(n) = a3(n) - xmult * a4(n-1) x(n) = ( b(n) - xmult * b(n-1) ) / a3(n) x(n-1) = ( b(n-1) - a4(n-1) * x(n) ) / a3(n-1) do i = n - 2, 1, -1 x(i) = ( b(i) - a4(i) * x(i+1) - a5(i) * x(i+2) ) / a3(i) end do return end subroutine r8_swap ( x, y ) !*****************************************************************************80 ! !! R8_SWAP swaps two real values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 01 May 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real ( kind = 8 ) X, Y. On output, the values of X and ! Y have been interchanged. ! implicit none real ( kind = 8 ) x real ( kind = 8 ) y real ( kind = 8 ) z z = x x = y y = z return end function r8_uniform_01 ( seed ) !*****************************************************************************80 ! !! R8_UNIFORM_01 is a portable pseudorandom number generator. ! ! Discussion: ! ! This routine implements the recursion ! ! seed = 16807 * seed mod ( 2**31 - 1 ) ! unif = seed / ( 2**31 - 1 ) ! ! The integer arithmetic never requires more than 32 bits, ! including a sign bit. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 11 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Paul Bratley, Bennett Fox, Linus Schrage, ! A Guide to Simulation, ! Springer Verlag, pages 201-202, 1983. ! ! Bennett Fox, ! Algorithm 647: ! Implementation and Relative Efficiency of Quasirandom ! Sequence Generators, ! ACM Transactions on Mathematical Software, ! Volume 12, Number 4, pages 362-376, 1986. ! ! Parameters: ! ! Input/output, integer ( kind = 4 ) SEED, the "seed" value, which should ! NOT be 0. (Otherwise, the output values of SEED and UNIFORM will be zero.) ! On output, SEED has been updated. ! ! Output, real ( kind = 8 ) R8_UNIFORM_01, a new pseudorandom variate, ! strictly between 0 and 1. ! implicit none integer ( kind = 4 ) k integer ( kind = 4 ) seed real ( kind = 8 ) r8_uniform_01 k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + 2147483647 end if ! ! Although SEED can be represented exactly as a 32 bit integer, ! it generally cannot be represented exactly as a 32 bit real number! ! r8_uniform_01 = real ( seed, kind = 8 ) * 4.656612875E-10 return end subroutine r83_mxv ( n, a, x, b ) !*****************************************************************************80 ! !! R83_MXV multiplies an R83 matrix times a vector. ! ! Discussion: ! ! The R83 storage format is used for a tridiagonal matrix. ! The superdiagonal is stored in entries (1,2:N), the diagonal in ! entries (2,1:N), and the subdiagonal in (3,1:N-1). Thus, the ! original matrix is "collapsed" vertically into the array. ! ! Example: ! ! Here is how an R83 matrix of order 5 would be stored: ! ! * A12 A23 A34 A45 ! A11 A22 A33 A44 A55 ! A21 A32 A43 A54 * ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 02 November 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the linear system. ! ! Input, real ( kind = 8 ) A(3,N), the R83 matrix. ! ! Input, real ( kind = 8 ) X(N), the vector to be multiplied by A. ! ! Output, real ( kind = 8 ) B(N), the product A * x. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a(3,n) real ( kind = 8 ) b(n) real ( kind = 8 ) x(n) b(1:n) = a(2,1:n) * x(1:n) b(1:n-1) = b(1:n-1) + a(1,2:n) * x(2:n) b(2:n) = b(2:n) + a(3,1:n-1) * x(1:n-1) return end subroutine r83_np_fs ( n, a, b, x ) !*****************************************************************************80 ! !! R83_NP_FS factors and solves an R83 system. ! ! Discussion: ! ! The R83 storage format is used for a tridiagonal matrix. ! The superdiagonal is stored in entries (1,2:N), the diagonal in ! entries (2,1:N), and the subdiagonal in (3,1:N-1). Thus, the ! original matrix is "collapsed" vertically into the array. ! ! This algorithm requires that each diagonal entry be nonzero. ! It does not use pivoting, and so can fail on systems that ! are actually nonsingular. ! ! Here is how an R83 matrix of order 5 would be stored: ! ! * A12 A23 A34 A45 ! A11 A22 A33 A44 A55 ! A21 A32 A43 A54 * ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 02 November 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the linear system. ! ! Input/output, real ( kind = 8 ) A(3,N). ! On input, the tridiagonal matrix. ! On output, the data in these vectors has been overwritten ! by factorization information. ! ! Input, real ( kind = 8 ) B(N), the right hand side of the linear system. ! ! Output, real ( kind = 8 ) X(N), the solution of the linear system. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a(3,n) real ( kind = 8 ) b(n) integer ( kind = 4 ) i real ( kind = 8 ) x(n) real ( kind = 8 ) xmult ! ! The diagonal entries can't be zero. ! do i = 1, n if ( a(2,i) == 0.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R83_NP_FS - Fatal error!' write ( *, '(a,i8,a)' ) ' A(2,', i, ') = 0.' return end if end do x(1:n) = b(1:n) do i = 2, n xmult = a(3,i-1) / a(2,i-1) a(2,i) = a(2,i) - xmult * a(1,i) x(i) = x(i) - xmult * x(i-1) end do x(n) = x(n) / a(2,n) do i = n - 1, 1, -1 x(i) = ( x(i) - a(1,i+1) * x(i+1) ) / a(2,i) end do return end subroutine r83_uniform ( n, seed, a ) !*****************************************************************************80 ! !! R83_UNIFORM randomizes an R83 matrix. ! ! Discussion: ! ! The R83 storage format is used for a tridiagonal matrix. ! The superdiagonal is stored in entries (1,2:N), the diagonal in ! entries (2,1:N), and the subdiagonal in (3,1:N-1). Thus, the ! original matrix is "collapsed" vertically into the array. ! ! Here is how an R83 matrix of order 5 would be stored: ! ! * A12 A23 A34 A45 ! A11 A22 A33 A44 A55 ! A21 A32 A43 A54 * ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 19 September 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the linear system. ! ! Input/output, integer ( kind = 4 ) SEED, a seed for the random number ! generator. ! ! Output, real ( kind = 8 ) A(3,N), the R83 matrix. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a(3,n) integer ( kind = 4 ) seed a(1,1) = 0.0D+00 call r8vec_uniform_01 ( n-1, seed, a(1,2:n) ) call r8vec_uniform_01 ( n, seed, a(2,1:n) ) call r8vec_uniform_01 ( n-1, seed, a(3,1:n-1) ) a(3,n) = 0.0D+00 return end subroutine r85_np_fs ( n, a, b, x ) !*****************************************************************************80 ! !! R85_NP_FS factors and solves an R85 linear system. ! ! Discussion: ! ! The R85 storage format represents a pentadiagonal matrix as a 5 ! by N array, in which each row corresponds to a diagonal, and ! column locations are preserved. Thus, the original matrix is ! "collapsed" vertically into the array. ! ! The factorization algorithm requires that each diagonal entry be nonzero. ! ! No pivoting is performed, and therefore the algorithm may fail ! in simple cases where the matrix is not singular. ! ! Example: ! ! Here is how an R85 matrix of order 6 would be stored: ! ! * * A13 A24 A35 A46 ! * A12 A23 A34 A45 A56 ! A11 A22 A33 A44 A55 A66 ! A21 A32 A43 A54 A65 * ! A31 A42 A53 A64 * * ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 17 September 2003 ! ! Author: ! ! Original FORTRAN77 version by Ward Cheney, David Kincaid. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Ward Cheney, David Kincaid, ! Numerical Mathematics and Computing, ! Brooks-Cole Publishing, 2004, ! ISBN: 0534201121. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the linear system. ! ! Input/output, real ( kind = 8 ) A(5,N), ! On input, the pentadiagonal matrix. ! On output, the data has been overwritten by factorization information. ! ! Input/output, real ( kind = 8 ) B(N). ! On input, B contains the right hand side of the linear system. ! On output, B has been overwritten by factorization information. ! ! Output, real ( kind = 8 ) X(N), the solution of the linear system. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a(5,n) real ( kind = 8 ) b(n) integer ( kind = 4 ) i real ( kind = 8 ) x(n) real ( kind = 8 ) xmult do i = 1, n if ( a(3,i) == 0.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R85_NP_FS - Fatal error!' write ( *, '(a,i8,a)' ) ' A(3,', i, ') = 0.' stop 1 end if end do do i = 2, n - 1 xmult = a(2,i) / a(3,i-1) a(3,i) = a(3,i) - xmult * a(4,i-1) a(4,i) = a(4,i) - xmult * a(5,i-1) b(i) = b(i) - xmult * b(i-1) xmult = a(1,i+1) / a(3,i-1) a(2,i+1) = a(2,i+1) - xmult * a(4,i-1) a(3,i+1) = a(3,i+1) - xmult * a(5,i-1) b(i+1) = b(i+1) - xmult * b(i-1) end do call r8vec_print ( n, b, ' Bpart1: ') xmult = a(2,n) / a(3,n-1) a(3,n) = a(3,n) - xmult * a(4,n-1) x(n) = ( b(n) - xmult * b(n-1) ) / a(3,n) x(n-1) = ( b(n-1) - a(4,n-1) * x(n) ) / a(3,n-1) do i = n - 2, 1, -1 x(i) = ( b(i) - a(4,i) * x(i+1) - a(5,i) * x(i+2) ) / a(3,i) end do return end subroutine r85_print ( n, a, title ) !*****************************************************************************80 ! !! R85_PRINT prints an R85 matrix. ! ! Discussion: ! ! The R85 storage format represents a pentadiagonal matrix as a 5 ! by N array, in which each row corresponds to a diagonal, and ! column locations are preserved. Thus, the original matrix is ! "collapsed" vertically into the array. ! ! Example: ! ! Here is how an R85 matrix of order 6 would be stored: ! ! * * A13 A24 A35 A46 ! * A12 A23 A34 A45 A56 ! A11 A22 A33 A44 A55 A66 ! A21 A32 A43 A54 A65 * ! A31 A42 A53 A64 * * ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 17 September 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! N must be positive. ! ! Input, real ( kind = 8 ) A(5,N), the matrix. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a(5,n) character ( len = * ) title call r85_print_some ( n, a, 1, 1, n, n, title ) return end subroutine r85_print_some ( n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! R85_PRINT_SOME prints some of an R85 matrix. ! ! Discussion: ! ! The R85 storage format represents a pentadiagonal matrix as a 5 ! by N array, in which each row corresponds to a diagonal, and ! column locations are preserved. Thus, the original matrix is ! "collapsed" vertically into the array. ! ! Example: ! ! Here is how an R85 matrix of order 6 would be stored: ! ! * * A13 A24 A35 A46 ! * A12 A23 A34 A45 A56 ! A11 A22 A33 A44 A55 A66 ! A21 A32 A43 A54 A65 * ! A31 A42 A53 A64 * * ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 05 January 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! N must be positive. ! ! Input, real ( kind = 8 ) A(5,N), the matrix. ! ! Input, integer ( kind = 4 ) ILO, JLO, IHI, JHI, the first row and ! column, and the last row and column, to be printed. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer ( kind = 4 ), parameter :: incx = 5 integer ( kind = 4 ) n real ( kind = 8 ) a(5,n) character ( len = 14 ) ctemp(incx) integer ( kind = 4 ) i integer ( kind = 4 ) i2hi integer ( kind = 4 ) i2lo integer ( kind = 4 ) ihi integer ( kind = 4 ) ilo integer ( kind = 4 ) inc integer ( kind = 4 ) j integer ( kind = 4 ) j2 integer ( kind = 4 ) j2hi integer ( kind = 4 ) j2lo integer ( kind = 4 ) jhi integer ( kind = 4 ) jlo character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) ! ! Print the columns of the matrix, in strips of 5. ! do j2lo = jlo, jhi, incx j2hi = j2lo + incx - 1 j2hi = min ( j2hi, n ) j2hi = min ( j2hi, jhi ) inc = j2hi + 1 - j2lo write ( *, '(a)' ) ' ' do j = j2lo, j2hi j2 = j + 1 - j2lo write ( ctemp(j2), '(i7,7x)' ) j end do write ( *, '('' Col: '',5a14)' ) ( ctemp(j2), j2 = 1, inc ) write ( *, '(a)' ) ' Row' write ( *, '(a)' ) ' ---' ! ! Determine the range of the rows in this strip. ! i2lo = max ( ilo, 1 ) i2lo = max ( i2lo, j2lo - 2 ) i2hi = min ( ihi, n ) i2hi = min ( i2hi, j2hi + 2 ) do i = i2lo, i2hi ! ! Print out (up to) 5 entries in row I, that lie in the current strip. ! do j2 = 1, inc j = j2lo - 1 + j2 if ( 2 < i-j .or. 2 < j-i ) then ctemp(j2) = ' ' else if ( j == i+2 ) then write ( ctemp(j2), '(g14.6)' ) a(1,j) else if ( j == i+1 ) then write ( ctemp(j2), '(g14.6)' ) a(2,j) else if ( j == i ) then write ( ctemp(j2), '(g14.6)' ) a(3,j) else if ( j == i-1 ) then write ( ctemp(j2), '(g14.6)' ) a(4,j) else if ( j == i-2 ) then write ( ctemp(j2), '(g14.6)' ) a(5,j) end if end do write ( *, '(i5,1x,5a14)' ) i, ( ctemp(j2), j2 = 1, inc ) end do end do return end subroutine r8ge_fs ( n, a, b, info ) !*****************************************************************************80 ! !! R8GE_FS factors and solves an R8GE system. ! ! Discussion: ! ! The R8GE storage format is used for a general M by N matrix. A storage ! space is made for each entry. The two dimensional logical ! array can be thought of as a vector of M*N entries, starting with ! the M entries in the column 1, then the M entries in column 2 ! and so on. Considered as a vector, the entry A(I,J) is then stored ! in vector location I+(J-1)*M. ! ! R8GE storage is used by LINPACK and LAPACK. ! ! R8GE_FS does not save the LU factors of the matrix, and hence cannot ! be used to efficiently solve multiple linear systems, or even to ! factor A at one time, and solve a single linear system at a later time. ! ! R8GE_FS uses partial pivoting, but no pivot vector is required. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 04 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! N must be positive. ! ! Input/output, real ( kind = 8 ) A(N,N). ! On input, A is the coefficient matrix of the linear system. ! On output, A is in unit upper triangular form, and ! represents the U factor of an LU factorization of the ! original coefficient matrix. ! ! Input/output, real ( kind = 8 ) B(N). ! On input, B is the right hand side of the linear system. ! On output, B is the solution of the linear system. ! ! Output, integer ( kind = 4 ) INFO, singularity flag. ! 0, no singularity detected. ! nonzero, the factorization failed on the INFO-th step. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a(n,n) real ( kind = 8 ) b(n) integer ( kind = 4 ) i integer ( kind = 4 ) info integer ( kind = 4 ) ipiv integer ( kind = 4 ) jcol real ( kind = 8 ) piv real ( kind = 8 ) row(n) real ( kind = 8 ) temp info = 0 do jcol = 1, n ! ! Find the maximum element in column I. ! piv = abs ( a(jcol,jcol) ) ipiv = jcol do i = jcol + 1, n if ( piv < abs ( a(i,jcol) ) ) then piv = abs ( a(i,jcol) ) ipiv = i end if end do if ( piv == 0.0D+00 ) then info = jcol write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8GE_FS - Fatal error!' write ( *, '(a,i8)' ) ' Zero pivot on step ', info stop 1 end if ! ! Switch rows JCOL and IPIV, and B. ! if ( jcol /= ipiv ) then row(1:n) = a(jcol,1:n) a(jcol,1:n) = a(ipiv,1:n) a(ipiv,1:n) = row(1:n) temp = b(jcol) b(jcol) = b(ipiv) b(ipiv) = temp end if ! ! Scale the pivot row. ! a(jcol,jcol+1:n) = a(jcol,jcol+1:n) / a(jcol,jcol) b(jcol) = b(jcol) / a(jcol,jcol) a(jcol,jcol) = 1.0D+00 ! ! Use the pivot row to eliminate lower entries in that column. ! do i = jcol + 1, n if ( a(i,jcol) /= 0.0D+00 ) then temp = - a(i,jcol) a(i,jcol) = 0.0D+00 a(i,jcol+1:n) = a(i,jcol+1:n) + temp * a(jcol,jcol+1:n) b(i) = b(i) + temp * b(jcol) end if end do end do ! ! Back solve. ! do jcol = n, 2, -1 b(1:jcol-1) = b(1:jcol-1) - a(1:jcol-1,jcol) * b(jcol) end do return end subroutine r8vec_bracket ( n, x, xval, left, right ) !*****************************************************************************80 ! !! R8VEC_BRACKET searches a sorted R8VEC for successive brackets of a value. ! ! Discussion: ! ! An R8VEC is an array of double precision real values. ! ! If the values in the vector are thought of as defining intervals ! on the real line, then this routine searches for the interval ! nearest to or containing the given value. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 06 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, length of input array. ! ! Input, real ( kind = 8 ) X(N), an array sorted into ascending order. ! ! Input, real ( kind = 8 ) XVAL, a value to be bracketed. ! ! Output, integer ( kind = 4 ) LEFT, RIGHT, the results of the search. ! Either: ! XVAL < X(1), when LEFT = 1, RIGHT = 2; ! X(N) < XVAL, when LEFT = N-1, RIGHT = N; ! or ! X(LEFT) <= XVAL <= X(RIGHT). ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) i integer ( kind = 4 ) left integer ( kind = 4 ) right real ( kind = 8 ) x(n) real ( kind = 8 ) xval do i = 2, n - 1 if ( xval < x(i) ) then left = i - 1 right = i return end if end do left = n - 1 right = n return end subroutine r8vec_bracket3 ( n, t, tval, left ) !*****************************************************************************80 ! !! R8VEC_BRACKET3 finds the interval containing or nearest a given value. ! ! Discussion: ! ! An R8VEC is an array of double precision real values. ! ! The routine always returns the index LEFT of the sorted array ! T with the property that either ! * T is contained in the interval [ T(LEFT), T(LEFT+1) ], or ! * T < T(LEFT) = T(1), or ! * T > T(LEFT+1) = T(N). ! ! The routine is useful for interpolation problems, where ! the abscissa must be located within an interval of data ! abscissas for interpolation, or the "nearest" interval ! to the (extreme) abscissa must be found so that extrapolation ! can be carried out. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 05 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, length of the input array. ! ! Input, real ( kind = 8 ) T(N), an array sorted into ascending order. ! ! Input, real ( kind = 8 ) TVAL, a value to be bracketed by entries of T. ! ! Input/output, integer ( kind = 4 ) LEFT. ! ! On input, if 1 <= LEFT <= N-1, LEFT is taken as a suggestion for the ! interval [ T(LEFT), T(LEFT+1) ] in which TVAL lies. This interval ! is searched first, followed by the appropriate interval to the left ! or right. After that, a binary search is used. ! ! On output, LEFT is set so that the interval [ T(LEFT), T(LEFT+1) ] ! is the closest to TVAL; it either contains TVAL, or else TVAL ! lies outside the interval [ T(1), T(N) ]. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) high integer ( kind = 4 ) left integer ( kind = 4 ) low integer ( kind = 4 ) mid real ( kind = 8 ) t(n) real ( kind = 8 ) tval ! ! Check the input data. ! if ( n < 2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8VEC_BRACKET3 - Fatal error!' write ( *, '(a)' ) ' N must be at least 2.' stop 1 end if ! ! If LEFT is not between 1 and N-1, set it to the middle value. ! if ( left < 1 .or. n - 1 < left ) then left = ( n + 1 ) / 2 end if ! ! CASE 1: TVAL < T(LEFT): ! Search for TVAL in [T(I), T(I+1)] for intervals I = 1 to LEFT-1. ! if ( tval < t(left) ) then if ( left == 1 ) then return else if ( left == 2 ) then left = 1 return else if ( t(left-1) <= tval ) then left = left - 1 return else if ( tval <= t(2) ) then left = 1 return end if ! ! ...Binary search for TVAL in [T(I), T(I+1)] for intervals I = 2 to LEFT-2. ! low = 2 high = left - 2 do if ( low == high ) then left = low return end if mid = ( low + high + 1 ) / 2 if ( t(mid) <= tval ) then low = mid else high = mid - 1 end if end do ! ! CASE2: T(LEFT+1) < TVAL: ! Search for TVAL in {T(I),T(I+1)] for intervals I = LEFT+1 to N-1. ! else if ( t(left+1) < tval ) then if ( left == n - 1 ) then return else if ( left == n - 2 ) then left = left + 1 return else if ( tval <= t(left+2) ) then left = left + 1 return else if ( t(n-1) <= tval ) then left = n - 1 return end if ! ! ...Binary search for TVAL in [T(I), T(I+1)] for intervals I = LEFT+2 to N-2. ! low = left + 2 high = n - 2 do if ( low == high ) then left = low return end if mid = ( low + high + 1 ) / 2 if ( t(mid) <= tval ) then low = mid else high = mid - 1 end if end do ! ! CASE3: T(LEFT) <= TVAL <= T(LEFT+1): ! T is in [T(LEFT), T(LEFT+1)], as the user said it might be. ! else end if return end function r8vec_distinct ( n, x ) !*****************************************************************************80 ! !! R8VEC_DISTINCT is true if the entries in an R8VEC are distinct. ! ! Discussion: ! ! An R8VEC is an array of double precision real values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 16 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of entries in the vector. ! ! Input, real ( kind = 8 ) X(N), the vector to be checked. ! ! Output, logical R8VEC_DISTINCT is TRUE if all N elements of X ! are distinct. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) i integer ( kind = 4 ) j logical r8vec_distinct real ( kind = 8 ) x(n) r8vec_distinct = .false. do i = 2, n do j = 1, i - 1 if ( x(i) == x(j) ) then return end if end do end do r8vec_distinct = .true. return end subroutine r8vec_even ( n, alo, ahi, a ) !*****************************************************************************80 ! !! R8VEC_EVEN returns N real values, evenly spaced between ALO and AHI. ! ! Discussion: ! ! An R8VEC is an array of double precision real values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 17 February 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of values. ! ! Input, real ( kind = 8 ) ALO, AHI, the low and high values. ! ! Output, real ( kind = 8 ) A(N), N evenly spaced values. ! Normally, A(1) = ALO and A(N) = AHI. ! However, if N = 1, then A(1) = 0.5*(ALO+AHI). ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a(n) real ( kind = 8 ) ahi real ( kind = 8 ) alo integer ( kind = 4 ) i if ( n == 1 ) then a(1) = 0.5D+00 * ( alo + ahi ) else do i = 1, n a(i) = ( real ( n - i, kind = 8 ) * alo & + real ( i - 1, kind = 8 ) * ahi ) & / real ( n - 1, kind = 8 ) end do end if return end subroutine r8vec_indicator ( n, a ) !*****************************************************************************80 ! !! R8VEC_INDICATOR sets an R8VEC to the indicator vector. ! ! Discussion: ! ! An R8VEC is an array of double precision real values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 09 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of elements of A. ! ! Output, real ( kind = 8 ) A(N), the array to be initialized. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a(n) integer ( kind = 4 ) i do i = 1, n a(i) = real ( i, kind = 8 ) end do return end subroutine r8vec_order_type ( n, a, order ) !*****************************************************************************80 ! !! R8VEC_ORDER_TYPE determines the order type of an R8VEC. ! ! Discussion: ! ! An R8VEC is an array of double precision real values. ! ! We assume the array is increasing or decreasing, and we want to ! verify that. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 20 July 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of entries of the array. ! ! Input, real ( kind = 8 ) A(N), the array to be checked. ! ! Output, integer ( kind = 4 ) ORDER, order indicator: ! -1, no discernable order; ! 0, all entries are equal; ! 1, ascending order; ! 2, strictly ascending order; ! 3, descending order; ! 4, strictly descending order. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a(n) integer ( kind = 4 ) i integer ( kind = 4 ) order ! ! Search for the first value not equal to A(1). ! i = 1 do i = i + 1 if ( n < i ) then order = 0 return end if if ( a(1) < a(i) ) then if ( i == 2 ) then order = 2 else order = 1 end if exit else if ( a(i) < a(1) ) then if ( i == 2 ) then order = 4 else order = 3 end if exit end if end do ! ! Now we have a "direction". Examine subsequent entries. ! do i = i + 1 if ( n < i ) then exit end if if ( order == 1 ) then if ( a(i) < a(i-1) ) then order = -1 exit end if else if ( order == 2 ) then if ( a(i) < a(i-1) ) then order = -1 exit else if ( a(i) == a(i-1) ) then order = 1 end if else if ( order == 3 ) then if ( a(i-1) < a(i) ) then order = -1 exit end if else if ( order == 4 ) then if ( a(i-1) < a(i) ) then order = -1 exit else if ( a(i) == a(i-1) ) then order = 3 end if end if end do return end subroutine r8vec_print ( n, a, title ) !*****************************************************************************80 ! !! R8VEC_PRINT prints an R8VEC. ! ! Discussion: ! ! An R8VEC is an array of double precision real values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 16 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of components of the vector. ! ! Input, real ( kind = 8 ) A(N), the vector to be printed. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a(n) integer ( kind = 4 ) i character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' do i = 1, n write ( *, '(i8,g14.6)' ) i, a(i) end do return end subroutine r8vec_sort_bubble_a ( n, a ) !*****************************************************************************80 ! !! R8VEC_SORT_BUBBLE_A ascending sorts an R8VEC using bubble sort. ! ! Discussion: ! ! An R8VEC is an array of double precision real values. ! ! Bubble sort is simple to program, but inefficient. It should not ! be used for large arrays. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 01 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of entries in the array. ! ! Input/output, real ( kind = 8 ) A(N). ! On input, an unsorted array. ! On output, the array has been sorted. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a(n) integer ( kind = 4 ) i integer ( kind = 4 ) j do i = 1, n-1 do j = i+1, n if ( a(j) < a(i) ) then call r8_swap ( a(i), a(j) ) end if end do end do return end subroutine r8vec_uniform_01 ( n, seed, r ) !*****************************************************************************80 ! !! R8VEC_UNIFORM_01 returns a unit pseudorandom R8VEC. ! ! Discussion: ! ! An R8VEC is an array of double precision real values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 19 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Paul Bratley, Bennett Fox, Linus Schrage, ! A Guide to Simulation, ! Springer Verlag, pages 201-202, 1983. ! ! Bennett Fox, ! Algorithm 647: ! Implementation and Relative Efficiency of Quasirandom ! Sequence Generators, ! ACM Transactions on Mathematical Software, ! Volume 12, Number 4, pages 362-376, 1986. ! ! Peter Lewis, Allen Goodman, James Miller, ! A Pseudo-Random Number Generator for the System/360, ! IBM Systems Journal, ! Volume 8, pages 136-143, 1969. ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, the number of entries in the vector. ! ! Input/output, integer ( kind = 4 ) SEED, the "seed" value, which should ! NOT be 0. On output, SEED has been updated. ! ! Output, real ( kind = 8 ) R(N), the vector of pseudorandom values. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) i integer ( kind = 4 ) k integer ( kind = 4 ) seed real ( kind = 8 ) r(n) do i = 1, n k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + 2147483647 end if r(i) = real ( seed, kind = 8 ) * 4.656612875D-10 end do return end subroutine r8vec_unique_count ( n, a, tol, unique_num ) !*****************************************************************************80 ! !! R8VEC_UNIQUE_COUNT counts the unique elements in an unsorted R8VEC. ! ! Discussion: ! ! An R8VEC is an array of double precision real values. ! ! Because the array is unsorted, this algorithm is O(N^2). ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 08 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of elements of A. ! ! Input, real ( kind = 8 ) A(N), the unsorted array to examine. ! ! Input, real ( kind = 8 ) TOL, a nonnegative tolerance for equality. ! Set it to 0.0 for the strictest test. ! ! Output, integer ( kind = 4 ) UNIQUE_NUM, the number of unique elements. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a(n) integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) unique_num real ( kind = 8 ) tol unique_num = 0 do i = 1, n unique_num = unique_num + 1 do j = 1, i - 1 if ( abs ( a(i) - a(j) ) <= tol ) then unique_num = unique_num - 1 exit end if end do end do return end subroutine spline_b_val ( ndata, tdata, ydata, tval, yval ) !*****************************************************************************80 ! !! SPLINE_B_VAL evaluates a cubic B spline approximant. ! ! Discussion: ! ! The cubic B spline will approximate the data, but is not ! designed to interpolate it. ! ! In effect, two "phantom" data values are appended to the data, ! so that the spline will interpolate the first and last data values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 11 February 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Carl deBoor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663. ! ! Parameters: ! ! Input, integer ( kind = 4 ) NDATA, the number of data values. ! ! Input, real ( kind = 8 ) TDATA(NDATA), the abscissas of the data. ! ! Input, real ( kind = 8 ) YDATA(NDATA), the data values. ! ! Input, real ( kind = 8 ) TVAL, a point at which the spline is ! to be evaluated. ! ! Output, real ( kind = 8 ) YVAL, the value of the function at TVAL. ! implicit none integer ( kind = 4 ) ndata real ( kind = 8 ) bval integer ( kind = 4 ) left integer ( kind = 4 ) right real ( kind = 8 ) tdata(ndata) real ( kind = 8 ) tval real ( kind = 8 ) u real ( kind = 8 ) ydata(ndata) real ( kind = 8 ) yval ! ! Find the nearest interval [ TDATA(LEFT), TDATA(RIGHT) ] to TVAL. ! call r8vec_bracket ( ndata, tdata, tval, left, right ) ! ! Evaluate the 5 nonzero B spline basis functions in the interval, ! weighted by their corresponding data values. ! u = ( tval - tdata(left) ) / ( tdata(right) - tdata(left) ) yval = 0.0D+00 ! ! B function associated with node LEFT - 1, (or "phantom node"), ! evaluated in its 4th interval. ! bval = ( ( ( - 1.0D+00 & * u + 3.0D+00 ) & * u - 3.0D+00 ) & * u + 1.0D+00 ) / 6.0D+00 if ( 0 < left-1 ) then yval = yval + ydata(left-1) * bval else yval = yval + ( 2.0D+00 * ydata(1) - ydata(2) ) * bval end if ! ! B function associated with node LEFT, ! evaluated in its third interval. ! bval = ( ( ( 3.0D+00 & * u - 6.0D+00 ) & * u + 0.0D+00 ) & * u + 4.0D+00 ) / 6.0D+00 yval = yval + ydata(left) * bval ! ! B function associated with node RIGHT, ! evaluated in its second interval. ! bval = ( ( ( - 3.0D+00 & * u + 3.0D+00 ) & * u + 3.0D+00 ) & * u + 1.0D+00 ) / 6.0D+00 yval = yval + ydata(right) * bval ! ! B function associated with node RIGHT+1, (or "phantom node"), ! evaluated in its first interval. ! bval = u**3 / 6.0D+00 if ( right+1 <= ndata ) then yval = yval + ydata(right+1) * bval else yval = yval + ( 2.0D+00 * ydata(ndata) - ydata(ndata-1) ) * bval end if return end subroutine spline_beta_val ( beta1, beta2, ndata, tdata, ydata, tval, yval ) !*****************************************************************************80 ! !! SPLINE_BETA_VAL evaluates a cubic beta spline approximant. ! ! Discussion: ! ! The cubic beta spline will approximate the data, but is not ! designed to interpolate it. ! ! If BETA1 = 1 and BETA2 = 0, the cubic beta spline will be the ! same as the cubic B spline approximant. ! ! With BETA1 = 1 and BETA2 large, the beta spline becomes more like ! a linear spline. ! ! In effect, two "phantom" data values are appended to the data, ! so that the spline will interpolate the first and last data values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 11 February 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) BETA1, the skew or bias parameter. ! BETA1 = 1 for no skew or bias. ! ! Input, real ( kind = 8 ) BETA2, the tension parameter. ! BETA2 = 0 for no tension. ! ! Input, integer ( kind = 4 ) NDATA, the number of data values. ! ! Input, real ( kind = 8 ) TDATA(NDATA), the abscissas of the data. ! ! Input, real ( kind = 8 ) YDATA(NDATA), the data values. ! ! Input, real ( kind = 8 ) TVAL, a point at which the spline is ! to be evaluated. ! ! Output, real ( kind = 8 ) YVAL, the value of the function at TVAL. ! implicit none integer ( kind = 4 ) ndata real ( kind = 8 ) a real ( kind = 8 ) b real ( kind = 8 ) beta1 real ( kind = 8 ) beta2 real ( kind = 8 ) bval real ( kind = 8 ) c real ( kind = 8 ) d real ( kind = 8 ) delta integer ( kind = 4 ) left integer ( kind = 4 ) right real ( kind = 8 ) tdata(ndata) real ( kind = 8 ) tval real ( kind = 8 ) u real ( kind = 8 ) ydata(ndata) real ( kind = 8 ) yval ! ! Find the nearest interval [ TDATA(LEFT), TDATA(RIGHT) ] to TVAL. ! call r8vec_bracket ( ndata, tdata, tval, left, right ) ! ! Evaluate the 5 nonzero beta spline basis functions in the interval, ! weighted by their corresponding data values. ! u = ( tval - tdata(left) ) / ( tdata(right) - tdata(left) ) delta = ( ( 2.0D+00 & * beta1 + 4.0D+00 ) & * beta1 + 4.0D+00 ) & * beta1 + 2.0D+00 + beta2 yval = 0.0D+00 ! ! Beta function associated with node LEFT - 1, (or "phantom node"), ! evaluated in its 4th interval. ! bval = 2.0D+00 * ( beta1 * ( 1.0D+00 - u ) )**3 / delta if ( 0 < left-1 ) then yval = yval + ydata(left-1) * bval else yval = yval + ( 2.0D+00 * ydata(1) - ydata(2) ) * bval end if ! ! Beta function associated with node LEFT, ! evaluated in its third interval. ! a = beta2 + ( 4.0D+00 + 4.0D+00 * beta1 ) * beta1 b = - 6.0D+00 * beta1 * ( 1.0D+00 - beta1 ) * ( 1.0D+00 + beta1 ) c = ( ( - 6.0D+00 & * beta1 - 6.0D+00 ) & * beta1 + 0.0D+00 ) & * beta1 - 3.0D+00 * beta2 d = ( ( + 2.0D+00 & * beta1 + 2.0D+00 ) & * beta1 + 2.0D+00 ) & * beta1 + 2.0D+00 * beta2 bval = ( a + u * ( b + u * ( c + u * d ) ) ) / delta yval = yval + ydata(left) * bval ! ! Beta function associated with node RIGHT, ! evaluated in its second interval. ! a = 2.0D+00 b = 6.0D+00 * beta1 c = 3.0D+00 * beta2 + 6.0D+00 * beta1 * beta1 d = - 2.0D+00 * ( 1.0D+00 + beta2 + beta1 + beta1 * beta1 ) bval = ( a + u * ( b + u * ( c + u * d ) ) ) / delta yval = yval + ydata(right) * bval ! ! Beta function associated with node RIGHT+1, (or "phantom node"), ! evaluated in its first interval. ! bval = 2.0D+00 * u**3 / delta if ( right + 1 <= ndata ) then yval = yval + ydata(right+1) * bval else yval = yval + ( 2.0D+00 * ydata(ndata) - ydata(ndata-1) ) * bval end if return end subroutine spline_bezier_val ( dim_num, interval_num, data_val, point_num, & point_t, point_val ) !*****************************************************************************80 ! !! SPLINE_BEZIER_VAL evaluates a cubic Bezier spline. ! ! Discussion: ! ! The cubic Bezier spline of N parts is defined by choosing ! 3*N+1 equally spaced T-abscissa values in the interval [0,N]. ! This defines N subintervals, each of length 1, and each containing ! 4 successives T abscissa values. ! ! At each abscissa value, a DIM_NUM-dimensional Bezier control ! value is assigned. Over each interval, a Bezier cubic function ! is used to define the value of the Bezier spline. To the left of ! the first interval, or to the right of the last interval, ! extrapolation may be used to extend the spline definition to ! the entire real line. ! ! Note that the Bezier spline will pass through the 1st, 4th, ! and in general 3*I+1 control values exactly. The other control ! values are not interpolating points. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 17 June 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) DIM_NUM, the spatial dimension. ! ! Input, integer ( kind = 4 ) INTERVAL_NUM, the number of intervals. ! ! Input, real ( kind = 8 ) DATA_VAL(DIM_NUM,3*INTERVAL_NUM+1), the control ! values. ! ! Input, integer ( kind = 4 ) POINT_NUM, the number of sample points at ! which the Bezier cubic spline is to be evaluated. ! ! Input, real ( kind = 8 ) POINT_T(POINT_NUM), the "T" values associated ! with the points. A value of T between 0 and 1, for instance, ! is associated with the first interval, and a value of T between ! INTERVAL_NUM-1 and INTERVAL_NUM is in the last interval. ! ! Output, real ( kind = 8 ) POINT_VAL(DIM_NUM,POINT_NUM), the value ! of the Bezier cubic spline at the sample points. ! implicit none integer ( kind = 4 ), parameter :: cubic = 3 integer ( kind = 4 ) interval_num integer ( kind = 4 ) dim_num integer ( kind = 4 ) point_num real ( kind = 8 ) bernstein_val(0:cubic) real ( kind = 8 ) data_val(dim_num,cubic*interval_num+1) integer ( kind = 4 ) dim integer ( kind = 4 ) interval integer ( kind = 4 ) offset integer ( kind = 4 ) point real ( kind = 8 ) point_t(point_num) real ( kind = 8 ) point_val(dim_num,point_num) real ( kind = 8 ) t real ( kind = 8 ) t_01 do point = 1, point_num t = point_t(point) interval = int ( t + 1 ) interval = max ( interval, 1 ) interval = min ( interval, interval_num ) offset = 1 + ( interval - 1 ) * cubic t_01 = t - real ( interval - 1, kind = 8 ) call bp01 ( cubic, t_01, bernstein_val ) do dim = 1, dim_num point_val(dim,point) = dot_product ( & data_val(dim,offset:offset+cubic), bernstein_val(0:cubic) ) end do end do return end subroutine spline_constant_val ( ndata, tdata, ydata, tval, yval ) !*****************************************************************************80 ! !! SPLINE_CONSTANT_VAL evaluates a piecewise constant spline at a point. ! ! Discussion: ! ! NDATA-1 points TDATA define NDATA intervals, with the first ! and last being semi-infinite. ! ! The value of the spline anywhere in interval I is YDATA(I). ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 16 November 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NDATA, the number of data points defining ! the spline. NDATA must be at least 1. ! ! Input, real ( kind = 8 ) TDATA(NDATA-1), the breakpoints. The values ! of TDATA should be distinct and increasing. ! ! Input, real ( kind = 8 ) YDATA(NDATA), the values of the spline in ! the intervals defined by the breakpoints. ! ! Input, real ( kind = 8 ) TVAL, the point at which the spline is ! to be evaluated. ! ! Output, real ( kind = 8 ) YVAL, the value of the spline at TVAL. ! implicit none integer ( kind = 4 ) ndata integer ( kind = 4 ) i real ( kind = 8 ) tdata(ndata-1) real ( kind = 8 ) tval real ( kind = 8 ) ydata(ndata) real ( kind = 8 ) yval ! ! Check NDATA. ! if ( ndata < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_CONSTANT_VAL - Fatal error!' write ( *, '(a)' ) ' NDATA < 1.' stop 1 end if do i = 1, ndata-1 if ( tval <= tdata(i) ) then yval = ydata(i) return end if end do yval = ydata(ndata) return end subroutine spline_cubic_set_full ( n, t, y, ibcbeg, ybcbeg, ibcend, ybcend, & ypp ) !*****************************************************************************80 ! !! SPLINE_CUBIC_SET_FULL computes second derivatives of piecewise cubic spline. ! ! Discussion: ! ! DELETE THIS FUNCTION WHEN YOU ARE SATISFIED WITH THE CURRENT ! VERSION OF SPLINE_CUBIC_SET. ! ! For data interpolation, the user must call SPLINE_CUBIC_SET to ! determine the second derivative data, passing in the data to be ! interpolated, and the desired boundary conditions. ! ! The data to be interpolated, plus the SPLINE_CUBIC_SET output, ! defines the spline. The user may then call SPLINE_CUBIC_VAL to ! evaluate the spline at any point. ! ! The cubic spline is a piecewise cubic polynomial. The intervals ! are determined by the "knots" or abscissas of the data to be ! interpolated. The cubic spline has continous first and second ! derivatives over the entire interval of interpolation. ! ! For any point T in the interval T(IVAL), T(IVAL+1), the form of ! the spline is ! ! SPL(T) = A(IVAL) ! + B(IVAL) * ( T - T(IVAL) ) ! + C(IVAL) * ( T - T(IVAL) )^2 ! + D(IVAL) * ( T - T(IVAL) )^3 ! ! If we assume that we know the values Y(*) and YPP(*), which represent ! the values and second derivatives of the spline at each knot, then ! the coefficients can be computed as: ! ! A(IVAL) = Y(IVAL) ! B(IVAL) = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) ) ! - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6 ! C(IVAL) = YPP(IVAL) / 2 ! D(IVAL) = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) ) ! ! Since the first derivative of the spline is ! ! SPL'(T) = B(IVAL) ! + 2 * C(IVAL) * ( T - T(IVAL) ) ! + 3 * D(IVAL) * ( T - T(IVAL) )^2, ! ! the requirement that the first derivative be continuous at interior ! knot I results in a total of N-2 equations, of the form: ! ! B(IVAL-1) + 2 C(IVAL-1) * (T(IVAL)-T(IVAL-1)) ! + 3 * D(IVAL-1) * (T(IVAL) - T(IVAL-1))^2 = B(IVAL) ! ! or, setting H(IVAL) = T(IVAL+1) - T(IVAL) ! ! ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1) ! - ( YPP(IVAL) + 2 * YPP(IVAL-1) ) * H(IVAL-1) / 6 ! + YPP(IVAL-1) * H(IVAL-1) ! + ( YPP(IVAL) - YPP(IVAL-1) ) * H(IVAL-1) / 2 ! = ! ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL) ! - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * H(IVAL) / 6 ! ! or ! ! YPP(IVAL-1) * H(IVAL-1) + 2 * YPP(IVAL) * ( H(IVAL-1) + H(IVAL) ) ! + YPP(IVAL) * H(IVAL) ! = ! 6 * ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL) ! - 6 * ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1) ! ! Boundary conditions must be applied at the first and last knots. ! The resulting tridiagonal system can be solved for the YPP values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 06 June 2013 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Carl deBoor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of data points; N must be ! at least 2. ! ! Input, real ( kind = 8 ) T(N), the points where data is specified. ! The values should be distinct, and increasing. ! ! Input, real ( kind = 8 ) Y(N), the data values to be interpolated. ! ! Input, integer ( kind = 4 ) IBCBEG, the left boundary condition flag: ! 0: the spline should be a quadratic over the first interval; ! 1: the first derivative at the left endpoint should be YBCBEG; ! 2: the second derivative at the left endpoint should be YBCBEG; ! 3: Not-a-knot: the third derivative is continuous at T(2). ! ! Input, real ( kind = 8 ) YBCBEG, the left boundary value, if needed. ! ! Input, integer ( kind = 4 ) IBCEND, the right boundary condition flag: ! 0: the spline should be a quadratic over the last interval; ! 1: the first derivative at the right endpoint should be YBCEND; ! 2: the second derivative at the right endpoint should be YBCEND; ! 3: Not-a-knot: the third derivative is continuous at T(N-1). ! ! Input, real ( kind = 8 ) YBCEND, the right boundary value, if needed. ! ! Output, real ( kind = 8 ) YPP(N), the second derivatives of ! the cubic spline. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a(n,n) real ( kind = 8 ) b(n) integer ( kind = 4 ) i integer ( kind = 4 ) ibcbeg integer ( kind = 4 ) ibcend integer ( kind = 4 ) info real ( kind = 8 ) t(n) real ( kind = 8 ) y(n) real ( kind = 8 ) ybcbeg real ( kind = 8 ) ybcend real ( kind = 8 ) ypp(n) ! ! Check. ! if ( n <= 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_CUBIC_SET - Fatal error!' write ( *, '(a)' ) ' The number of knots must be at least 2.' write ( *, '(a,i8)' ) ' The input value of N = ', n stop 1 end if do i = 1, n - 1 if ( t(i+1) <= t(i) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_CUBIC_SET - Fatal error!' write ( *, '(a)' ) ' The knots must be strictly increasing, but' write ( *, '(a,i8,a,g14.6)' ) ' T(', i,') = ', t(i) write ( *, '(a,i8,a,g14.6)' ) ' T(',i+1,') = ', t(i+1) stop 1 end if end do ! ! Zero out the matrix. ! a(1:n,1:n) = 0.0D+00 ! ! Set the first equation. ! if ( ibcbeg == 0 ) then b(1) = 0.0D+00 a(1,1) = 1.0D+00 a(1,2) = -1.0D+00 else if ( ibcbeg == 1 ) then b(1) = ( y(2) - y(1) ) / ( t(2) - t(1) ) - ybcbeg a(1,1) = ( t(2) - t(1) ) / 3.0D+00 a(1,2) = ( t(2) - t(1) ) / 6.0D+00 else if ( ibcbeg == 2 ) then b(1) = ybcbeg a(1,1) = 1.0D+00 a(1,2) = 0.0D+00 else if ( ibcbeg == 3 ) then b(1) = 0.0D+00 a(1,1) = - ( t(3) - t(2) ) a(1,2) = ( t(3) - t(1) ) a(1,3) = - ( t(2) - t(1) ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_CUBIC_SET - Fatal error!' write ( *, '(a)' ) ' The boundary flag IBCBEG must be 0, 1, 2 or 3.' write ( *, '(a,i8)' ) ' The input value is IBCBEG = ', ibcbeg stop 1 end if ! ! Set the intermediate equations. ! do i = 2, n - 1 b(i) = ( y(i+1) - y(i) ) / ( t(i+1) - t(i) ) & - ( y(i) - y(i-1) ) / ( t(i) - t(i-1) ) a(i,i-1) = ( t(i) - t(i-1) ) / 6.0D+00 a(i,i) = ( t(i+1) - t(i-1) ) / 3.0D+00 a(i,i+1) = ( t(i+1) - t(i) ) / 6.0D+00 end do ! ! Set the last equation. ! if ( ibcend == 0 ) then b(n) = 0.0D+00 a(n,n-1) = -1.0D+00 a(n,n) = 1.0D+00 else if ( ibcend == 1 ) then b(n) = ybcend - ( y(n) - y(n-1) ) / ( t(n) - t(n-1) ) a(n,n-1) = ( t(n) - t(n-1) ) / 6.0D+00 a(n,n) = ( t(n) - t(n-1) ) / 3.0D+00 else if ( ibcend == 2 ) then b(n) = ybcend a(n,n-1) = 0.0D+00 a(n,n) = 1.0D+00 else if ( ibcend == 3 ) then b(n) = 0.0D+00 a(n,n-2) = - ( t(n) - t(n-1) ) a(n,n-1) = ( t(n) - t(n-2) ) a(n,n) = - ( t(n-1) - t(n-2) ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_CUBIC_SET - Fatal error!' write ( *, '(a)' ) ' The boundary flag IBCEND must be 0, 1, 2 or 3.' write ( *, '(a,i8)' ) ' The input value is IBCEND = ', ibcend stop 1 end if ! ! Special case: ! N = 2, IBCBEG = IBCEND = 0. ! if ( n == 2 .and. ibcbeg == 0 .and. ibcend == 0 ) then ypp(1) = 0.0D+00 ypp(2) = 0.0D+00 ! ! Solve the linear system. ! else call r8ge_fs ( n, a, b, info ) ypp = b end if return end subroutine spline_cubic_set ( n, t, y, ibcbeg, ybcbeg, ibcend, ybcend, ypp ) !*****************************************************************************80 ! !! SPLINE_CUBIC_SET computes the second derivatives of a piecewise cubic spline. ! ! Discussion: ! ! For data interpolation, the user must call SPLINE_CUBIC_SET to ! determine the second derivative data, passing in the data to be ! interpolated, and the desired boundary conditions. ! ! The data to be interpolated, plus the SPLINE_CUBIC_SET output, ! defines the spline. The user may then call SPLINE_CUBIC_VAL to ! evaluate the spline at any point. ! ! The cubic spline is a piecewise cubic polynomial. The intervals ! are determined by the "knots" or abscissas of the data to be ! interpolated. The cubic spline has continous first and second ! derivatives over the entire interval of interpolation. ! ! For any point T in the interval T(IVAL), T(IVAL+1), the form of ! the spline is ! ! SPL(T) = A(IVAL) ! + B(IVAL) * ( T - T(IVAL) ) ! + C(IVAL) * ( T - T(IVAL) )^2 ! + D(IVAL) * ( T - T(IVAL) )^3 ! ! If we assume that we know the values Y(*) and YPP(*), which represent ! the values and second derivatives of the spline at each knot, then ! the coefficients can be computed as: ! ! A(IVAL) = Y(IVAL) ! B(IVAL) = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) ) ! - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6 ! C(IVAL) = YPP(IVAL) / 2 ! D(IVAL) = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) ) ! ! Since the first derivative of the spline is ! ! SPL'(T) = B(IVAL) ! + 2 * C(IVAL) * ( T - T(IVAL) ) ! + 3 * D(IVAL) * ( T - T(IVAL) )^2, ! ! the requirement that the first derivative be continuous at interior ! knot I results in a total of N-2 equations, of the form: ! ! B(IVAL-1) + 2 C(IVAL-1) * (T(IVAL)-T(IVAL-1)) ! + 3 * D(IVAL-1) * (T(IVAL) - T(IVAL-1))^2 = B(IVAL) ! ! or, setting H(IVAL) = T(IVAL+1) - T(IVAL) ! ! ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1) ! - ( YPP(IVAL) + 2 * YPP(IVAL-1) ) * H(IVAL-1) / 6 ! + YPP(IVAL-1) * H(IVAL-1) ! + ( YPP(IVAL) - YPP(IVAL-1) ) * H(IVAL-1) / 2 ! = ! ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL) ! - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * H(IVAL) / 6 ! ! or ! ! YPP(IVAL-1) * H(IVAL-1) + 2 * YPP(IVAL) * ( H(IVAL-1) + H(IVAL) ) ! + YPP(IVAL) * H(IVAL) ! = ! 6 * ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL) ! - 6 * ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1) ! ! Boundary conditions must be applied at the first and last knots. ! The resulting tridiagonal system can be solved for the YPP values. ! ! Anton Reinhard corrected the assignments: ! a2(i) = ( t(i+1) - t(i) ) / 6.0D+00 ! a3(i) = ( t(i+1) - t(i-1) ) / 3.0D+00 ! a4(i) = ( t(i) - t(i-1) ) / 6.0D+00 ! to: ! a2(i) = ( t(i) - t(i-1) ) / 6.0D+00 ! a3(i) = ( t(i+1) - t(i-1) ) / 3.0D+00 ! a4(i) = ( t(i+1) - t(i) ) / 6.0D+00 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 17 August 2020 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Carl deBoor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663. ! ! Input: ! ! integer ( kind = 4 ) N, the number of data points; N must be ! at least 2. ! ! real ( kind = 8 ) T(N), the points where data is specified. ! The values should be distinct, and increasing. ! ! real ( kind = 8 ) Y(N), the data values to be interpolated. ! ! integer ( kind = 4 ) IBCBEG, the left boundary condition flag: ! 0: the spline should be a quadratic over the first interval; ! 1: the first derivative at the left endpoint should be YBCBEG; ! 2: the second derivative at the left endpoint should be YBCBEG; ! 3: Not-a-knot: the third derivative is continuous at T(2). ! ! real ( kind = 8 ) YBCBEG, the left boundary value, if needed. ! ! integer ( kind = 4 ) IBCEND, the right boundary condition flag: ! 0: the spline should be a quadratic over the last interval; ! 1: the first derivative at the right endpoint should be YBCEND; ! 2: the second derivative at the right endpoint should be YBCEND; ! 3: Not-a-knot: the third derivative is continuous at T(N-1). ! ! real ( kind = 8 ) YBCEND, the right boundary value, if needed. ! ! Output: ! ! real ( kind = 8 ) YPP(N), the second derivatives of ! the cubic spline. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a1(n) real ( kind = 8 ) a2(n) real ( kind = 8 ) a3(n) real ( kind = 8 ) a4(n) real ( kind = 8 ) a5(n) real ( kind = 8 ) b(n) integer ( kind = 4 ) i integer ( kind = 4 ) ibcbeg integer ( kind = 4 ) ibcend real ( kind = 8 ) t(n) real ( kind = 8 ) y(n) real ( kind = 8 ) ybcbeg real ( kind = 8 ) ybcend real ( kind = 8 ) ypp(n) ! ! Check. ! if ( n <= 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_CUBIC_SET - Fatal error!' write ( *, '(a)' ) ' The number of knots must be at least 2.' write ( *, '(a,i8)' ) ' The input value of N = ', n stop 1 end if do i = 1, n - 1 if ( t(i+1) <= t(i) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_CUBIC_SET - Fatal error!' write ( *, '(a)' ) ' The knots must be strictly increasing, but' write ( *, '(a,i8,a,g14.6)' ) ' T(', i,') = ', t(i) write ( *, '(a,i8,a,g14.6)' ) ' T(',i+1,') = ', t(i+1) stop 1 end if end do ! ! Zero out the matrix. ! a1(1:n) = 0.0D+00 a2(1:n) = 0.0D+00 a3(1:n) = 0.0D+00 a4(1:n) = 0.0D+00 a5(1:n) = 0.0D+00 ! ! Set the first equation. ! if ( ibcbeg == 0 ) then b(1) = 0.0D+00 a3(1) = 1.0D+00 a4(1) = -1.0D+00 else if ( ibcbeg == 1 ) then b(1) = ( y(2) - y(1) ) / ( t(2) - t(1) ) - ybcbeg a3(1) = ( t(2) - t(1) ) / 3.0D+00 a4(1) = ( t(2) - t(1) ) / 6.0D+00 else if ( ibcbeg == 2 ) then b(1) = ybcbeg a3(1) = 1.0D+00 a4(1) = 0.0D+00 else if ( ibcbeg == 3 ) then b(1) = 0.0D+00 a3(1) = - ( t(3) - t(2) ) a4(1) = ( t(3) - t(1) ) a5(1) = - ( t(2) - t(1) ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_CUBIC_SET - Fatal error!' write ( *, '(a)' ) ' The boundary flag IBCBEG must be 0, 1, 2 or 3.' write ( *, '(a,i8)' ) ' The input value is IBCBEG = ', ibcbeg stop 1 end if ! ! Set the intermediate equations. ! Note that these lines have been corrected as suggested ! by Anton Reinhard, 17 August 2020. ! do i = 2, n - 1 b(i) = ( y(i+1) - y(i) ) / ( t(i+1) - t(i) ) & - ( y(i) - y(i-1) ) / ( t(i) - t(i-1) ) a2(i) = ( t(i) - t(i-1) ) / 6.0D+00 a3(i) = ( t(i+1) - t(i-1) ) / 3.0D+00 a4(i) = ( t(i+1) - t(i) ) / 6.0D+00 end do ! ! Set the last equation. ! if ( ibcend == 0 ) then b(n) = 0.0D+00 a2(n) = -1.0D+00 a3(n) = 1.0D+00 else if ( ibcend == 1 ) then b(n) = ybcend - ( y(n) - y(n-1) ) / ( t(n) - t(n-1) ) a2(n) = ( t(n) - t(n-1) ) / 6.0D+00 a3(n) = ( t(n) - t(n-1) ) / 3.0D+00 else if ( ibcend == 2 ) then b(n) = ybcend a2(n) = 0.0D+00 a3(n) = 1.0D+00 else if ( ibcend == 3 ) then b(n) = 0.0D+00 a1(n) = - ( t(n) - t(n-1) ) a2(n) = ( t(n) - t(n-2) ) a3(n) = - ( t(n-1) - t(n-2) ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_CUBIC_SET - Fatal error!' write ( *, '(a)' ) ' The boundary flag IBCEND must be 0, 1, 2 or 3.' write ( *, '(a,i8)' ) ' The input value is IBCEND = ', ibcend stop 1 end if ! ! Special case: ! N = 2, IBCBEG = IBCEND = 0. ! if ( n == 2 .and. ibcbeg == 0 .and. ibcend == 0 ) then ypp(1) = 0.0D+00 ypp(2) = 0.0D+00 ! ! Solve the linear system. ! else call penta ( n, a1, a2, a3, a4, a5, b, ypp ) end if return end subroutine spline_cubic_val ( n, t, y, ypp, tval, yval, ypval, yppval ) !*****************************************************************************80 ! !! SPLINE_CUBIC_VAL evaluates a piecewise cubic spline at a point. ! ! Discussion: ! ! SPLINE_CUBIC_SET must have already been called to define the ! values of YPP. ! ! For any point T in the interval T(IVAL), T(IVAL+1), the form of ! the spline is ! ! SPL(T) = A ! + B * ( T - T(IVAL) ) ! + C * ( T - T(IVAL) )^2 ! + D * ( T - T(IVAL) )^3 ! ! Here: ! A = Y(IVAL) ! B = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) ) ! - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6 ! C = YPP(IVAL) / 2 ! D = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) ) ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 20 November 2000 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Carl deBoor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663. ! ! Input: ! ! integer ( kind = 4 ) N, the number of data values. ! ! real ( kind = 8 ) T(N), the knot values. ! ! real ( kind = 8 ) Y(N), the data values at the knots. ! ! real ( kind = 8 ) YPP(N), the second derivatives of the ! spline at the knots. ! ! real ( kind = 8 ) TVAL, a point, typically between T(1) and ! T(N), at which the spline is to be evalulated. If TVAL lies outside ! this range, extrapolation is used. ! ! Output: ! ! real ( kind = 8 ) YVAL, YPVAL, YPPVAL, the value of the spline, and ! its first two derivatives at TVAL. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) dt real ( kind = 8 ) h integer ( kind = 4 ) left integer ( kind = 4 ) right real ( kind = 8 ) t(n) real ( kind = 8 ) tval real ( kind = 8 ) y(n) real ( kind = 8 ) ypp(n) real ( kind = 8 ) yppval real ( kind = 8 ) ypval real ( kind = 8 ) yval ! ! Determine the interval [T(LEFT), T(RIGHT)] that contains TVAL. ! Values below T(1) or above T(N) use extrapolation. ! call r8vec_bracket ( n, t, tval, left, right ) ! ! Evaluate the polynomial. ! dt = tval - t(left) h = t(right) - t(left) yval = y(left) & + dt * ( ( y(right) - y(left) ) / h & - ( ypp(right) / 6.0D+00 + ypp(left) / 3.0D+00 ) * h & + dt * ( 0.5D+00 * ypp(left) & + dt * ( ( ypp(right) - ypp(left) ) / ( 6.0D+00 * h ) ) ) ) ypval = ( y(right) - y(left) ) / h & - ( ypp(right) / 6.0D+00 + ypp(left) / 3.0D+00 ) * h & + dt * ( ypp(left) & + dt * ( 0.5D+00 * ( ypp(right) - ypp(left) ) / h ) ) yppval = ypp(left) + dt * ( ypp(right) - ypp(left) ) / h return end subroutine spline_cubic_val2 ( n, t, y, ypp, left, tval, yval, ypval, yppval ) !*****************************************************************************80 ! !! SPLINE_CUBIC_VAL2 evaluates a piecewise cubic spline at a point. ! ! Discussion: ! ! This routine is a modification of SPLINE_CUBIC_VAL; it allows the ! user to speed up the code by suggesting the appropriate T interval ! to search first. ! ! SPLINE_CUBIC_SET must have already been called to define the ! values of YPP. ! ! In the LEFT interval, let RIGHT = LEFT+1. The form of the spline is ! ! SPL(T) = ! A ! + B * ( T - T(LEFT) ) ! + C * ( T - T(LEFT) )^2 ! + D * ( T - T(LEFT) )^3 ! ! Here: ! A = Y(LEFT) ! B = ( Y(RIGHT) - Y(LEFT) ) / ( T(RIGHT) - T(LEFT) ) ! - ( YPP(RIGHT) + 2 * YPP(LEFT) ) * ( T(RIGHT) - T(LEFT) ) / 6 ! C = YPP(LEFT) / 2 ! D = ( YPP(RIGHT) - YPP(LEFT) ) / ( 6 * ( T(RIGHT) - T(LEFT) ) ) ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 24 February 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Carl deBoor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of knots. ! ! Input, real ( kind = 8 ) T(N), the knot values. ! ! Input, real ( kind = 8 ) Y(N), the data values at the knots. ! ! Input, real ( kind = 8 ) YPP(N), the second derivatives of the spline at ! the knots. ! ! Input/output, integer ( kind = 4 ) LEFT, the suggested T interval to ! search. LEFT should be between 1 and N-1. If LEFT is not in this range, ! then its value will be ignored. On output, LEFT is set to the ! actual interval in which TVAL lies. ! ! Input, real ( kind = 8 ) TVAL, a point, typically between T(1) and T(N), at ! which the spline is to be evalulated. If TVAL lies outside ! this range, extrapolation is used. ! ! Output, real ( kind = 8 ) YVAL, YPVAL, YPPVAL, the value of the spline, and ! its first two derivatives at TVAL. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) dt real ( kind = 8 ) h integer ( kind = 4 ) left integer ( kind = 4 ) right real ( kind = 8 ) t(n) real ( kind = 8 ) tval real ( kind = 8 ) y(n) real ( kind = 8 ) ypp(n) real ( kind = 8 ) yppval real ( kind = 8 ) ypval real ( kind = 8 ) yval ! ! Determine the interval [T(LEFT), T(RIGHT)] that contains TVAL. ! ! What you want from R8VEC_BRACKET3 is that TVAL is to be computed ! by the data in interval [T(LEFT), T(RIGHT)]. ! call r8vec_bracket3 ( n, t, tval, left ) right = left + 1 ! ! In the interval LEFT, the polynomial is in terms of a normalized ! coordinate ( DT / H ) between 0 and 1. ! dt = tval - t(left) h = t(right) - t(left) yval = y(left) + dt * ( ( y(right) - y(left) ) / h & - ( ypp(right) / 6.0D+00 + ypp(left) / 3.0D+00 ) * h & + dt * ( 0.5D+00 * ypp(left) & + dt * ( ( ypp(right) - ypp(left) ) / ( 6.0D+00 * h ) ) ) ) ypval = ( y(right) - y(left) ) / h & - ( ypp(right) / 6.0D+00 + ypp(left) / 3.0D+00 ) * h & + dt * ( ypp(left) & + dt * ( 0.5D+00 * ( ypp(right) - ypp(left) ) / h ) ) yppval = ypp(left) + dt * ( ypp(right) - ypp(left) ) / h return end subroutine spline_hermite_set ( ndata, tdata, ydata, ypdata, c ) !*****************************************************************************80 ! !! SPLINE_HERMITE_SET sets up a piecewise cubic Hermite interpolant. ! ! Discussion: ! ! Once the array C is computed, then in the interval ! (TDATA(I), TDATA(I+1)), the interpolating Hermite polynomial ! is given by ! ! SVAL(TVAL) = C(1,I) ! + ( TVAL - TDATA(I) ) * ( C(2,I) ! + ( TVAL - TDATA(I) ) * ( C(3,I) ! + ( TVAL - TDATA(I) ) * C(4,I) ) ) ! ! This is algorithm CALCCF of Conte and deBoor. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 11 February 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Samuel Conte, Carl deBoor, ! Elementary Numerical Analysis, ! Second Edition, ! McGraw Hill, 1972, ! ISBN: 07-012446-4, ! LC: QA297.C65. ! ! Parameters: ! ! Input, integer ( kind = 4 ) NDATA, the number of data points. ! NDATA must be at least 2. ! ! Input, real ( kind = 8 ) TDATA(NDATA), the abscissas of the data points. ! The entries of TDATA are assumed to be strictly increasing. ! ! Input, real ( kind = 8 ) Y(NDATA), YP(NDATA), the value of the ! function and its derivative at TDATA(1:NDATA). ! ! Output, real ( kind = 8 ) C(4,NDATA), the coefficients of the ! Hermite polynomial. ! C(1,1:NDATA) = Y(1:NDATA) and C(2,1:NDATA) = YP(1:NDATA). ! C(3,1:NDATA-1) and C(4,1:NDATA-1) are the quadratic and cubic ! coefficients. ! implicit none integer ( kind = 4 ) ndata real ( kind = 8 ) c(4,ndata) real ( kind = 8 ) divdif1 real ( kind = 8 ) divdif3 real ( kind = 8 ) dt integer ( kind = 4 ) i real ( kind = 8 ) tdata(ndata) real ( kind = 8 ) ydata(ndata) real ( kind = 8 ) ypdata(ndata) ! ! Check NDATA. ! if ( ndata < 2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_HERMITE_SET - Fatal error!' write ( *, '(a)' ) ' NDATA < 2.' stop 1 end if c(1,1:ndata) = ydata(1:ndata) c(2,1:ndata) = ypdata(1:ndata) do i = 1, ndata - 1 dt = tdata(i+1) - tdata(i) divdif1 = ( c(1,i+1) - c(1,i) ) / dt divdif3 = c(2,i) + c(2,i+1) - 2.0D+00 * divdif1 c(3,i) = ( divdif1 - c(2,i) - divdif3 ) / dt c(4,i) = divdif3 / ( dt * dt ) end do c(3,ndata) = 0.0D+00 c(4,ndata) = 0.0D+00 return end subroutine spline_hermite_val ( ndata, tdata, c, tval, sval, spval ) !*****************************************************************************80 ! !! SPLINE_HERMITE_VAL evaluates a piecewise cubic Hermite interpolant. ! ! Discussion: ! ! SPLINE_HERMITE_SET must be called first, to set up the ! spline data from the raw function and derivative data. ! ! In the interval (TDATA(I), TDATA(I+1)), the interpolating ! Hermite polynomial is given by ! ! SVAL(TVAL) = C(1,I) ! + ( TVAL - TDATA(I) ) * ( C(2,I) ! + ( TVAL - TDATA(I) ) * ( C(3,I) ! + ( TVAL - TDATA(I) ) * C(4,I) ) ) ! ! and ! ! SVAL'(TVAL) = C(2,I) ! + ( TVAL - TDATA(I) ) * ( 2 * C(3,I) ! + ( TVAL - TDATA(I) ) * 3 * C(4,I) ) ! ! This is algorithm PCUBIC of Conte and deBoor. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 11 February 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Samuel Conte, Carl deBoor, ! Elementary Numerical Analysis, ! Second Edition, ! McGraw Hill, 1972, ! ISBN: 07-012446-4, ! LC: QA297.C65. ! ! Parameters: ! ! Input, integer ( kind = 4 ) NDATA, the number of data points. ! NDATA must be at least 2. ! ! Input, real ( kind = 8 ) TDATA(NDATA), the abscissas of the data points. ! The entries of TDATA are assumed to be strictly increasing. ! ! Input, real ( kind = 8 ) C(4,NDATA), the coefficient data computed by ! SPLINE_HERMITE_SET. ! ! Input, real ( kind = 8 ) TVAL, the point where the interpolant is to ! be evaluated. ! ! Output, real ( kind = 8 ) SVAL, SPVAL, the value of the interpolant ! and its derivative at TVAL. ! implicit none integer ( kind = 4 ) ndata real ( kind = 8 ) c(4,ndata) real ( kind = 8 ) dt integer ( kind = 4 ) left integer ( kind = 4 ) right real ( kind = 8 ) spval real ( kind = 8 ) sval real ( kind = 8 ) tdata(ndata) real ( kind = 8 ) tval ! ! Check NDATA. ! if ( ndata < 2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_HERMITE_VAL - Fatal error!' write ( *, '(a)' ) ' NDATA < 2.' stop 1 end if ! ! Find the interval [ TDATA(LEFT), TDATA(RIGHT) ] that contains ! or is nearest to TVAL. ! call r8vec_bracket ( ndata, tdata, tval, left, right ) ! ! Evaluate the cubic polynomial. ! dt = tval - tdata(left) sval = c(1,left) + dt * ( c(2,left) + dt * ( c(3,left) + dt * c(4,left) ) ) spval = c(2,left) + dt * ( 2.0D+00 * c(3,left) + dt * 3.0D+00 * c(4,left) ) return end subroutine spline_linear_int ( ndata, tdata, ydata, a, b, int_val ) !*****************************************************************************80 ! !! SPLINE_LINEAR_INT evaluates the integral of a piecewise linear spline. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 01 November 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NDATA, the number of data points defining ! the spline. NDATA must be at least 2. ! ! Input, real ( kind = 8 ) TDATA(NDATA), YDATA(NDATA), the values of ! the independent and dependent variables at the data points. The ! values of TDATA should be distinct and increasing. ! ! Input, real ( kind = 8 ) A, B, the interval over which the ! integral is desired. ! ! Output, real ( kind = 8 ) INT_VAL, the value of the integral. ! implicit none integer ( kind = 4 ) ndata real ( kind = 8 ) a real ( kind = 8 ) a_copy integer ( kind = 4 ) a_left integer ( kind = 4 ) a_right real ( kind = 8 ) b real ( kind = 8 ) b_copy integer ( kind = 4 ) b_left integer ( kind = 4 ) b_right integer ( kind = 4 ) i_left real ( kind = 8 ) int_val real ( kind = 8 ) tdata(ndata) real ( kind = 8 ) tval real ( kind = 8 ) ydata(ndata) real ( kind = 8 ) yp real ( kind = 8 ) yval ! ! Check NDATA. ! if ( ndata < 2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_LINEAR_INT - Fatal error!' write ( *, '(a)' ) ' NDATA < 2.' stop 1 end if int_val = 0.0D+00 if ( a == b ) then return end if a_copy = min ( a, b ) b_copy = max ( a, b ) ! ! Find the interval [ TDATA(A_LEFT), TDATA(A_RIGHT) ] that contains, or is ! nearest to, A. ! call r8vec_bracket ( ndata, tdata, a_copy, a_left, a_right ) ! ! Find the interval [ TDATA(B_LEFT), TDATA(B_RIGHT) ] that contains, or is ! nearest to, B. ! call r8vec_bracket ( ndata, tdata, b_copy, b_left, b_right ) ! ! If A and B are in the same interval... ! if ( a_left == b_left ) then tval = ( a_copy + b_copy ) / 2.0D+00 yp = ( ydata(a_right) - ydata(a_left) ) / & ( tdata(a_right) - tdata(a_left) ) yval = ydata(a_left) + ( tval - tdata(a_left) ) * yp int_val = yval * ( b_copy - a_copy ) return end if ! ! Otherwise, integrate from: ! ! A to TDATA(A_RIGHT), ! TDATA(A_RIGHT) to TDATA(A_RIGHT+1),... ! TDATA(B_LEFT-1) to TDATA(B_LEFT), ! TDATA(B_LEFT) to B. ! ! Use the fact that the integral of a linear function is the ! value of the function at the midpoint times the width of the interval. ! tval = ( a_copy + tdata(a_right) ) / 2.0D+00 yp = ( ydata(a_right) - ydata(a_left) ) / & ( tdata(a_right) - tdata(a_left) ) yval = ydata(a_left) + ( tval - tdata(a_left) ) * yp int_val = int_val + yval * ( tdata(a_right) - a_copy ) do i_left = a_right, b_left - 1 tval = ( tdata(i_left+1) + tdata(i_left) ) / 2.0D+00 yp = ( ydata(i_left+1) - ydata(i_left) ) / & ( tdata(i_left+1) - tdata(i_left) ) yval = ydata(i_left) + ( tval - tdata(i_left) ) * yp int_val = int_val + yval * ( tdata(i_left + 1) - tdata(i_left) ) end do tval = ( tdata(b_left) + b_copy ) / 2.0D+00 yp = ( ydata(b_right) - ydata(b_left) ) / & ( tdata(b_right) - tdata(b_left) ) yval = ydata(b_left) + ( tval - tdata(b_left) ) * yp int_val = int_val + yval * ( b_copy - tdata(b_left) ) if ( b < a ) then int_val = - int_val end if return end subroutine spline_linear_intset ( n, int_x, int_v, data_x, data_y ) !*****************************************************************************80 ! !! SPLINE_LINEAR_INTSET: piecewise linear spline with given integral properties. ! ! Discussion: ! ! The user has in mind an interval, divided by INT_N+1 points into ! INT_N intervals. A linear spline is to be constructed, ! with breakpoints at the centers of each interval, and extending ! continuously to the left of the first and right of the last ! breakpoints. The constraint on the linear spline is that it is ! required that it have a given integral value over each interval. ! ! A tridiagonal linear system of equations is solved for the ! values of the spline at the breakpoints. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 27 January 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of intervals. ! ! Input, real ( kind = 8 ) INT_X(N+1), the points that define the intervals. ! Interval I lies between INT_X(I) and INT_X(I+1). ! ! Input, real ( kind = 8 ) INT_V(N), the desired value of the integral of the ! linear spline over each interval. ! ! Output, real ( kind = 8 ) DATA_X(N), DATA_Y(N), the values of the ! independent and dependent variables at the data points. The values ! of DATA_X are the interval midpoints. The values of DATA_Y are ! determined in such a way that the exact integral of the linear ! spline over interval I is equal to INT_V(I). ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a(3,n) real ( kind = 8 ) data_x(n) real ( kind = 8 ) data_y(n) real ( kind = 8 ) int_v(n) real ( kind = 8 ) int_x(n+1) ! ! Set up the easy stuff. ! data_x(1:n) = 0.5D+00 * ( int_x(1:n) + int_x(2:n+1) ) ! ! Set up C, D, E, the coefficients of the linear system. ! a(3,1:n-2) = 1.0D+00 & - ( 0.5D+00 * ( data_x(2:n-1) + int_x(2:n-1) ) - data_x(1:n-2) ) & / ( data_x(2:n-1) - data_x(1:n-2) ) a(3,n-1) = 0.0D+00 a(3,n) = 0.0D+00 a(2,1) = int_x(2) - int_x(1) a(2,2:n-1) = 1.0D+00 & + ( 0.5D+00 * ( data_x(2:n-1) + int_x(2:n-1) ) & - data_x(1:n-2) ) & / ( data_x(2:n-1) - data_x(1:n-2) ) & - ( 0.5D+00 * ( data_x(2:n-1) + int_x(3:n) ) - data_x(2:n-1) ) & / ( data_x(3:n) - data_x(2:n-1) ) a(2,n) = int_x(n+1) - int_x(n) a(1,1) = 0.0D+00 a(1,2) = 0.0D+00 a(1,3:n) = ( 0.5D+00 * ( data_x(2:n-1) + int_x(3:n) ) & - data_x(2:n-1) ) / ( data_x(3:n) - data_x(2:n-1) ) ! ! Set up DATA_Y, which begins as the right hand side of the linear system. ! data_y(1) = int_v(1) data_y(2:n-1) = 2.0D+00 * int_v(2:n-1) / ( int_x(3:n) - int_x(2:n-1) ) data_y(n) = int_v(n) ! ! Solve the linear system. ! call r83_np_fs ( n, a, data_y, data_y ) return end subroutine spline_linear_val ( ndata, tdata, ydata, tval, yval, ypval ) !*****************************************************************************80 ! !! SPLINE_LINEAR_VAL evaluates a piecewise linear spline at a point. ! ! Discussion: ! ! Because of the extremely simple form of the linear spline, ! the raw data points ( TDATA(I), YDATA(I)) can be used directly to ! evaluate the spline at any point. No processing of the data ! is required. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 06 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NDATA, the number of data points defining ! the spline. NDATA must be at least 2. ! ! Input, real ( kind = 8 ) TDATA(NDATA), YDATA(NDATA), the values of ! the independent and dependent variables at the data points. The ! values of TDATA should be distinct and increasing. ! ! Input, real ( kind = 8 ) TVAL, the point at which the spline is ! to be evaluated. ! ! Output, real ( kind = 8 ) YVAL, YPVAL, the value of the spline and ! its first derivative dYdT at TVAL. YPVAL is not reliable if TVAL ! is exactly equal to TDATA(I) for some I. ! implicit none integer ( kind = 4 ) ndata integer ( kind = 4 ) left integer ( kind = 4 ) right real ( kind = 8 ) tdata(ndata) real ( kind = 8 ) tval real ( kind = 8 ) ydata(ndata) real ( kind = 8 ) ypval real ( kind = 8 ) yval ! ! Check NDATA. ! if ( ndata < 2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_LINEAR_VAL - Fatal error!' write ( *, '(a)' ) ' NDATA < 2.' stop 1 end if ! ! Find the interval [ TDATA(LEFT), TDATA(RIGHT) ] that contains, or is ! nearest to, TVAL. ! call r8vec_bracket ( ndata, tdata, tval, left, right ) ! ! Now evaluate the piecewise linear function. ! ypval = ( ydata(right) - ydata(left) ) / ( tdata(right) - tdata(left) ) yval = ydata(left) + ( tval - tdata(left) ) * ypval return end subroutine spline_overhauser_nonuni_val ( ndata, tdata, ydata, tval, yval ) !*****************************************************************************80 ! !! SPLINE_OVERHAUSER_NONUNI_VAL evaluates the nonuniform Overhauser spline. ! ! Discussion: ! ! The nonuniformity refers to the fact that the abscissa values ! need not be uniformly spaced. ! ! Thanks to Doug Fortune for pointing out that the point distances ! used to define ALPHA and BETA should be the Euclidean distances ! between the points. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 08 May 2007 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NDATA, the number of data points. ! 3 <= NDATA is required. ! ! Input, real ( kind = 8 ) TDATA(NDATA), the abscissas of the data points. ! The values of TDATA are assumed to be distinct and increasing. ! ! Input, real ( kind = 8 ) YDATA(NDATA), the data values. ! ! Input, real ( kind = 8 ) TVAL, the value where the spline is to ! be evaluated. ! ! Output, real ( kind = 8 ) YVAL, the value of the spline at TVAL. ! implicit none integer ( kind = 4 ) ndata real ( kind = 8 ) alpha real ( kind = 8 ) beta real ( kind = 8 ) d21 real ( kind = 8 ) d32 real ( kind = 8 ) d43 integer ( kind = 4 ) left real ( kind = 8 ) mbasis(4,4) real ( kind = 8 ) mbasis_l(3,3) real ( kind = 8 ) mbasis_r(3,3) integer ( kind = 4 ) right real ( kind = 8 ) tdata(ndata) real ( kind = 8 ) tval real ( kind = 8 ) ydata(ndata) real ( kind = 8 ) yval ! ! Check NDATA. ! if ( ndata < 3 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_OVERHAUSER_NONUNI_VAL - Fatal error!' write ( *, '(a)' ) ' NDATA < 3.' stop 1 end if ! ! Find the nearest interval [ TDATA(LEFT), TDATA(RIGHT) ] to TVAL. ! call r8vec_bracket ( ndata, tdata, tval, left, right ) ! ! Evaluate the spline in the given interval. ! if ( left == 1 ) then d21 = sqrt ( ( tdata(2) - tdata(1) )**2 & + ( ydata(2) - ydata(1) )**2 ) d32 = sqrt ( ( tdata(3) - tdata(2) )**2 & + ( ydata(3) - ydata(2) )**2 ) alpha = d21 / ( d32 + d21 ) call basis_matrix_overhauser_nul ( alpha, mbasis_l ) call basis_matrix_tmp ( left, 3, mbasis_l, ndata, tdata, ydata, tval, yval ) else if ( left < ndata - 1 ) then d21 = sqrt ( ( tdata(left) - tdata(left-1) )**2 & + ( ydata(left) - ydata(left-1) )**2 ) d32 = sqrt ( ( tdata(left+1) - tdata(left) )**2 & + ( ydata(left+1) - ydata(left) )**2 ) d43 = sqrt ( ( tdata(left+2) - tdata(left+1) )**2 & + ( ydata(left+2) - ydata(left+1) )**2 ) alpha = d21 / ( d32 + d21 ) beta = d32 / ( d43 + d32 ) call basis_matrix_overhauser_nonuni ( alpha, beta, mbasis ) call basis_matrix_tmp ( left, 4, mbasis, ndata, tdata, ydata, tval, yval ) else if ( left == ndata - 1 ) then d32 = sqrt ( ( tdata(ndata-1) - tdata(ndata-2) )**2 & + ( ydata(ndata-1) - ydata(ndata-2) )**2 ) d43 = sqrt ( ( tdata(ndata) - tdata(ndata-1) )**2 & + ( ydata(ndata) - ydata(ndata-1) )**2 ) beta = d32 / ( d43 + d32 ) call basis_matrix_overhauser_nur ( beta, mbasis_r ) call basis_matrix_tmp ( left, 3, mbasis_r, ndata, tdata, ydata, tval, yval ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_OVERHAUSER_NONUNI_VAL - Fatal error!' write ( *, '(a,i8)' ) ' Nonsensical value of LEFT = ', left write ( *, '(a,i8)' ) ' but 0 < LEFT < NDATA = ', ndata write ( *, '(a)' ) ' is required.' stop 1 end if return end subroutine spline_overhauser_uni_val ( ndata, tdata, ydata, tval, yval ) !*****************************************************************************80 ! !! SPLINE_OVERHAUSER_UNI_VAL evaluates the uniform Overhauser spline. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 06 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NDATA, the number of data points. ! NDATA must be at least 3. ! ! Input, real ( kind = 8 ) TDATA(NDATA), the abscissas of the data points. ! The values of TDATA are assumed to be distinct and increasing. ! This routine also assumes that the values of TDATA are uniformly ! spaced; for instance, TDATA(1) = 10, TDATA(2) = 11, TDATA(3) = 12... ! ! Input, real ( kind = 8 ) YDATA(NDATA), the data values. ! ! Input, real ( kind = 8 ) TVAL, the value where the spline is to ! be evaluated. ! ! Output, real ( kind = 8 ) YVAL, the value of the spline at TVAL. ! implicit none integer ( kind = 4 ) ndata integer ( kind = 4 ) left real ( kind = 8 ) mbasis(4,4) real ( kind = 8 ) mbasis_l(3,3) real ( kind = 8 ) mbasis_r(3,3) integer ( kind = 4 ) right real ( kind = 8 ) tdata(ndata) real ( kind = 8 ) tval real ( kind = 8 ) ydata(ndata) real ( kind = 8 ) yval ! ! Check NDATA. ! if ( ndata < 3 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_OVERHAUSER_UNI_VAL - Fatal error!' write ( *, '(a)' ) ' NDATA < 3.' stop 1 end if ! ! Find the nearest interval [ TDATA(LEFT), TDATA(RIGHT) ] to TVAL. ! call r8vec_bracket ( ndata, tdata, tval, left, right ) ! ! Evaluate the spline in the given interval. ! if ( left == 1 ) then call basis_matrix_overhauser_uni_l ( mbasis_l ) call basis_matrix_tmp ( left, 3, mbasis_l, ndata, tdata, ydata, tval, yval ) else if ( left < ndata - 1 ) then call basis_matrix_overhauser_uni ( mbasis ) call basis_matrix_tmp ( left, 4, mbasis, ndata, tdata, ydata, tval, yval ) else if ( left == ndata - 1 ) then call basis_matrix_overhauser_uni_r ( mbasis_r ) call basis_matrix_tmp ( left, 3, mbasis_r, ndata, tdata, ydata, tval, yval ) end if return end subroutine spline_overhauser_val ( dim_num, ndata, tdata, ydata, tval, yval ) !*****************************************************************************80 ! !! SPLINE_OVERHAUSER_VAL evaluates an Overhauser spline. ! ! Discussion: ! ! Over the first and last intervals, the Overhauser spline is a ! quadratic. In the intermediate intervals, it is a piecewise cubic. ! The Overhauser spline is also known as the Catmull-Rom spline. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 08 April 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! JA Brewer, DC Anderson, ! Visual Interaction with Overhauser Curves and Surfaces, ! SIGGRAPH 77, ! in Proceedings of the 4th Annual Conference on Computer Graphics ! and Interactive Techniques, ! ASME, July 1977, pages 132-137. ! ! Edwin Catmull, Raphael Rom, ! A Class of Local Interpolating Splines, ! in Computer Aided Geometric Design, ! edited by Robert Barnhill, Richard Reisenfeld, ! Academic Press, 1974, pages 317-326, ! ISBN: 0120790505. ! ! David Rogers, Alan Adams, ! Mathematical Elements of Computer Graphics, ! Second Edition, ! McGraw Hill, 1989, ! ISBN: 0070535299. ! ! Parameters: ! ! Input, integer ( kind = 4 ) DIM_NUM, the dimension of a single data point. ! DIM_NUM must be at least 1. ! ! Input, integer ( kind = 4 ) NDATA, the number of data points. ! NDATA must be at least 3. ! ! Input, real ( kind = 8 ) TDATA(NDATA), the abscissas of the data ! points. The values in TDATA must be in strictly ascending order. ! ! Input, real ( kind = 8 ) YDATA(DIM_NUM,NDATA), the data points ! corresponding to the abscissas. ! ! Input, real ( kind = 8 ) TVAL, the abscissa value at which the spline ! is to be evaluated. Normally, TDATA(1) <= TVAL <= T(NDATA), and ! the data will be interpolated. For TVAL outside this range, ! extrapolation will be used. ! ! Output, real ( kind = 8 ) YVAL(DIM_NUM), the value of the spline at TVAL. ! implicit none integer ( kind = 4 ) ndata integer ( kind = 4 ) dim_num integer ( kind = 4 ) left integer ( kind = 4 ) order integer ( kind = 4 ) right real ( kind = 8 ) tdata(ndata) real ( kind = 8 ) tval real ( kind = 8 ) ydata(dim_num,ndata) real ( kind = 8 ) yl(dim_num) real ( kind = 8 ) yr(dim_num) real ( kind = 8 ) yval(dim_num) ! ! Check. ! call r8vec_order_type ( ndata, tdata, order ) if ( order /= 2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_OVERHAUSER_VAL - Fatal error!' write ( *, '(a)' ) ' The data abscissas are not strictly ascending.' stop 1 end if if ( ndata < 3 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_OVERHAUSER_VAL - Fatal error!' write ( *, '(a)' ) ' NDATA < 3.' stop 1 end if if ( dim_num < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_OVERHAUSER_VAL - Fatal error!' write ( *, '(a)' ) ' DIM_NUM < 1.' stop 1 end if ! ! Locate the abscissa interval T(LEFT), T(LEFT+1) nearest to or ! containing TVAL. ! call r8vec_bracket ( ndata, tdata, tval, left, right ) ! ! Evaluate the "left hand" quadratic defined at T(LEFT-1), T(LEFT), T(RIGHT). ! if ( 0 < left-1 ) then call parabola_val2 ( dim_num, ndata, tdata, ydata, left-1, tval, yl ) end if ! ! Evaluate the "right hand" quadratic defined at T(LEFT), T(RIGHT), T(RIGHT+1). ! if ( right+1 <= ndata ) then call parabola_val2 ( dim_num, ndata, tdata, ydata, left, tval, yr ) end if ! ! Average the quadratics. ! if ( left == 1 ) then yval(1:dim_num) = yr(1:dim_num) else if ( right < ndata ) then yval(1:dim_num) = & ( ( tdata(right) - tval ) * yl(1:dim_num) & + ( tval - tdata(left) ) * yr(1:dim_num) ) & / ( tdata(right) - tdata(left) ) else yval(1:dim_num) = yl(1:dim_num) end if return end subroutine spline_pchip_set ( n, x, f, d ) !*****************************************************************************80 ! !! SPLINE_PCHIP_SET sets derivatives for a piecewise cubic Hermite interpolant. ! ! Discussion: ! ! This routine computes what would normally be called a Hermite ! interpolant. However, the user is only required to supply function ! values, not derivative values as well. This routine computes ! "suitable" derivative values, so that the resulting Hermite interpolant ! has desirable shape and monotonicity properties. ! ! The interpolant will have an extremum at each point where ! monotonicity switches direction. ! ! The resulting piecewise cubic Hermite function may be evaluated ! by SPLINE_PCHIP_VAL. ! ! This routine was originally named "PCHIM". ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 15 December 2008 ! ! Author: ! ! Original FORTRAN77 version by Fred Fritsch. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Fred Fritsch, Ralph Carlson, ! Monotone Piecewise Cubic Interpolation, ! SIAM Journal on Numerical Analysis, ! Volume 17, Number 2, April 1980, pages 238-246. ! ! Fred Fritsch, Judy Butland, ! A Method for Constructing Local Monotone Piecewise Cubic Interpolants, ! SIAM Journal on Scientific and Statistical Computing, ! Volume 5, Number 2, 1984, pages 300-304. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of data points. N must be ! at least 2. ! ! Input, real ( kind = 8 ) X(N), the strictly increasing independent ! variable values. ! ! Input, real ( kind = 8 ) F(N), dependent variable values to be ! interpolated. F(I) is the value corresponding to X(I). ! This routine is designed for monotonic data, but it will work for any ! F array. It will force extrema at points where monotonicity switches ! direction. ! ! Output, real ( kind = 8 ) D(N), the derivative values at the ! data points. If the data are monotonic, these values will determine ! a monotone cubic Hermite function. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) d(n) real ( kind = 8 ) del1 real ( kind = 8 ) del2 real ( kind = 8 ) dmax real ( kind = 8 ) dmin real ( kind = 8 ) drat1 real ( kind = 8 ) drat2 real ( kind = 8 ) dsave real ( kind = 8 ) f(n) real ( kind = 8 ) h1 real ( kind = 8 ) h2 real ( kind = 8 ) hsum real ( kind = 8 ) hsumt3 integer ( kind = 4 ) i integer ( kind = 4 ) ierr integer ( kind = 4 ) nless1 ! real ( kind = 8 ) pchst real ( kind = 8 ) temp real ( kind = 8 ) w1 real ( kind = 8 ) w2 real ( kind = 8 ) x(n) ! ! Check the arguments. ! if ( n < 2 ) then ierr = -1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_PCHIP_SET - Fatal error!' write ( *, '(a)' ) ' Number of data points less than 2.' stop 1 end if do i = 2, n if ( x(i) <= x(i-1) ) then ierr = -3 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_PCHIP_SET - Fatal error!' write ( *, '(a)' ) ' X array not strictly increasing.' stop 1 end if end do ierr = 0 nless1 = n - 1 h1 = x(2) - x(1) del1 = ( f(2) - f(1) ) / h1 dsave = del1 ! ! Special case N=2, use linear interpolation. ! if ( n == 2 ) then d(1) = del1 d(n) = del1 return end if ! ! Normal case, 3 <= N. ! h2 = x(3) - x(2) del2 = ( f(3) - f(2) ) / h2 ! ! Set D(1) via non-centered three point formula, adjusted to be ! shape preserving. ! hsum = h1 + h2 w1 = ( h1 + hsum ) / hsum w2 = -h1 / hsum d(1) = w1 * del1 + w2 * del2 if ( pchst ( d(1), del1 ) <= 0.0D+00 ) then d(1) = 0.0D+00 ! ! Need do this check only if monotonicity switches. ! else if ( pchst ( del1, del2 ) < 0.0D+00 ) then dmax = 3.0D+00 * del1 if ( abs ( dmax ) < abs ( d(1) ) ) then d(1) = dmax end if end if ! ! Loop through interior points. ! do i = 2, nless1 if ( 2 < i ) then h1 = h2 h2 = x(i+1) - x(i) hsum = h1 + h2 del1 = del2 del2 = ( f(i+1) - f(i) ) / h2 end if ! ! Set D(I)=0 unless data are strictly monotonic. ! d(i) = 0.0D+00 temp = pchst ( del1, del2 ) if ( temp < 0.0D+00 ) then ierr = ierr + 1 dsave = del2 ! ! Count number of changes in direction of monotonicity. ! else if ( temp == 0.0D+00 ) then if ( del2 /= 0.0D+00 ) then if ( pchst ( dsave, del2 ) < 0.0D+00 ) then ierr = ierr + 1 end if dsave = del2 end if ! ! Use Brodlie modification of Butland formula. ! else hsumt3 = 3.0D+00 * hsum w1 = ( hsum + h1 ) / hsumt3 w2 = ( hsum + h2 ) / hsumt3 dmax = max ( abs ( del1 ), abs ( del2 ) ) dmin = min ( abs ( del1 ), abs ( del2 ) ) drat1 = del1 / dmax drat2 = del2 / dmax d(i) = dmin / ( w1 * drat1 + w2 * drat2 ) end if end do ! ! Set D(N) via non-centered three point formula, adjusted to be ! shape preserving. ! w1 = -h2 / hsum w2 = ( h2 + hsum ) / hsum d(n) = w1 * del1 + w2 * del2 if ( pchst ( d(n), del2 ) <= 0.0D+00 ) then d(n) = 0.0D+00 else if ( pchst ( del1, del2 ) < 0.0D+00 ) then ! ! Need do this check only if monotonicity switches. ! dmax = 3.0D+00 * del2 if ( abs ( dmax ) < abs ( d(n) ) ) then d(n) = dmax end if end if return end subroutine spline_pchip_val ( n, x, f, d, ne, xe, fe ) !*****************************************************************************80 ! !! SPLINE_PCHIP_VAL evaluates a piecewise cubic Hermite function. ! ! Description: ! ! This routine may be used by itself for Hermite interpolation, or as an ! evaluator for SPLINE_PCHIP_SET. ! ! This routine evaluates the cubic Hermite function at the points XE. ! ! Most of the coding between the call to CHFEV and the end of ! the IR loop could be eliminated if it were permissible to ! assume that XE is ordered relative to X. ! ! CHFEV does not assume that X1 is less than X2. Thus, it would ! be possible to write a version of SPLINE_PCHIP_VAL that assumes a strictly ! decreasing X array by simply running the IR loop backwards ! and reversing the order of appropriate tests. ! ! The present code has a minor bug, which I have decided is not ! worth the effort that would be required to fix it. ! If XE contains points in [X(N-1),X(N)], followed by points less than ! X(N-1), followed by points greater than X(N), the extrapolation points ! will be counted (at least) twice in the total returned in IERR. ! ! The evaluation will be most efficient if the elements of XE are ! increasing relative to X; that is, for all J <= K, ! X(I) <= XE(J) ! implies ! X(I) <= XE(K). ! ! If any of the XE are outside the interval [X(1),X(N)], ! values are extrapolated from the nearest extreme cubic, ! and a warning error is returned. ! ! This routine was originally called "PCHFE". ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 14 August 2005 ! ! Author: ! ! Original FORTRAN77 version by Fred Fritsch. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Fred Fritsch, Ralph Carlson, ! Monotone Piecewise Cubic Interpolation, ! SIAM Journal on Numerical Analysis, ! Volume 17, Number 2, April 1980, pages 238-246. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of data points. N must be ! at least 2. ! ! Input, real ( kind = 8 ) X(N), the strictly increasing independent ! variable values. ! ! Input, real ( kind = 8 ) F(N), the function values. ! ! Input, real ( kind = 8 ) D(N), the derivative values. ! ! Input, integer ( kind = 4 ) NE, the number of evaluation points. ! ! Input, real ( kind = 8 ) XE(NE), points at which the function is to ! be evaluated. ! ! Output, real ( kind = 8 ) FE(NE), the values of the cubic Hermite ! function at XE. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) ne real ( kind = 8 ) d(n) real ( kind = 8 ) f(n) real ( kind = 8 ) fe(ne) integer ( kind = 4 ) i integer ( kind = 4 ) ierc integer ( kind = 4 ) ierr integer ( kind = 4 ) ir integer ( kind = 4 ) j integer ( kind = 4 ) j_first integer ( kind = 4 ) j_new integer ( kind = 4 ) j_save integer ( kind = 4 ) next(2) integer ( kind = 4 ) nj real ( kind = 8 ) x(n) real ( kind = 8 ) xe(ne) ! ! Check arguments. ! if ( n < 2 ) then ierr = -1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_PCHIP_VAL - Fatal error!' write ( *, '(a)' ) ' Number of data points less than 2.' stop 1 end if do i = 2, n if ( x(i) <= x(i-1) ) then ierr = -3 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_PCHIP_VAL - Fatal error!' write ( *, '(a)' ) ' X array not strictly increasing.' stop 1 end if end do if ( ne < 1 ) then ierr = -4 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_PCHIP_VAL - Fatal error!' write ( *, '(a)' ) ' Number of evaluation points less than 1.' return end if ierr = 0 ! ! Loop over intervals. ! The interval index is IL = IR-1. ! The interval is X(IL) <= X < X(IR). ! j_first = 1 ir = 2 do ! ! Skip out of the loop if have processed all evaluation points. ! if ( ne < j_first ) then exit end if ! ! Locate all points in the interval. ! j_save = ne + 1 do j = j_first, ne if ( x(ir) <= xe(j) ) then j_save = j if ( ir == n ) then j_save = ne + 1 end if exit end if end do ! ! Have located first point beyond interval. ! j = j_save nj = j - j_first ! ! Skip evaluation if no points in interval. ! if ( nj /= 0 ) then ! ! Evaluate cubic at XE(J_FIRST:J-1). ! call chfev ( x(ir-1), x(ir), f(ir-1), f(ir), d(ir-1), d(ir), & nj, xe(j_first:j-1), fe(j_first:j-1), next, ierc ) if ( ierc < 0 ) then ierr = -5 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_PCHIP_VAL - Fatal error!' write ( *, '(a)' ) ' Error return from CHFEV.' stop 1 end if ! ! In the current set of XE points, there are NEXT(2) to the right of X(IR). ! if ( next(2) /= 0 ) then if ( ir < n ) then ierr = -5 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_PCHIP_VAL - Fatal error!' write ( *, '(a)' ) ' IR < N.' stop 1 end if ! ! These are actually extrapolation points. ! ierr = ierr + next(2) end if ! ! In the current set of XE points, there are NEXT(1) to the left of X(IR-1). ! if ( next(1) /= 0 ) then ! ! These are actually extrapolation points. ! if ( ir <= 2 ) then ierr = ierr + next(1) else j_new = -1 do i = j_first, j - 1 if ( xe(i) < x(ir-1) ) then j_new = i exit end if end do if ( j_new == -1 ) then ierr = -5 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_PCHIP_VAL - Fatal error!' write ( *, '(a)' ) ' Could not bracket the data point.' stop 1 end if ! ! Reset J. This will be the new J_FIRST. ! j = j_new ! ! Now find out how far to back up in the X array. ! do i = 1, ir-1 if ( xe(j) < x(i) ) then exit end if end do ! ! At this point, either XE(J) < X(1) or X(i-1) <= XE(J) < X(I) . ! ! Reset IR, recognizing that it will be incremented before cycling. ! ir = max ( 1, i-1 ) end if end if j_first = j end if ir = ir + 1 if ( n < ir ) then exit end if end do return end subroutine spline_quadratic_val ( ndata, tdata, ydata, tval, yval, ypval ) !*****************************************************************************80 ! !! SPLINE_QUADRATIC_VAL evaluates a piecewise quadratic spline at a point. ! ! Discussion: ! ! Because of the simple form of a piecewise quadratic spline, ! the raw data points ( TDATA(I), YDATA(I)) can be used directly to ! evaluate the spline at any point. No processing of the data ! is required. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 24 October 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NDATA, the number of data points defining ! the spline. NDATA should be odd and at least 3. ! ! Input, real ( kind = 8 ) TDATA(NDATA), YDATA(NDATA), the values of ! the independent and dependent variables at the data points. The ! values of TDATA should be distinct and increasing. ! ! Input, real ( kind = 8 ) TVAL, the point at which the spline is to ! be evaluated. ! ! Output, real ( kind = 8 ) YVAL, YPVAL, the value of the spline and ! its first derivative dYdT at TVAL. YPVAL is not reliable if TVAL ! is exactly equal to TDATA(I) for some I. ! implicit none integer ( kind = 4 ) ndata real ( kind = 8 ) dif1 real ( kind = 8 ) dif2 integer ( kind = 4 ) left integer ( kind = 4 ) right real ( kind = 8 ) t1 real ( kind = 8 ) t2 real ( kind = 8 ) t3 real ( kind = 8 ) tdata(ndata) real ( kind = 8 ) tval real ( kind = 8 ) y1 real ( kind = 8 ) y2 real ( kind = 8 ) y3 real ( kind = 8 ) ydata(ndata) real ( kind = 8 ) ypval real ( kind = 8 ) yval if ( ndata < 3 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_QUADRATIC_VAL - Fatal error!' write ( *, '(a)' ) ' NDATA < 3.' stop 1 end if if ( mod ( ndata, 2 ) == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_QUADRATIC_VAL - Fatal error!' write ( *, '(a)' ) ' NDATA must be odd.' stop 1 end if ! ! Find the interval [ TDATA(LEFT), TDATA(RIGHT) ] that contains, or is ! nearest to, TVAL. ! call r8vec_bracket ( ndata, tdata, tval, left, right ) ! ! Force LEFT to be odd. ! if ( mod ( left, 2 ) == 0 ) then left = left - 1 end if ! ! Copy out the three abscissas. ! t1 = tdata(left) t2 = tdata(left+1) t3 = tdata(left+2) if ( t2 <= t1 .or. t3 <= t2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINE_QUADRATIC_VAL - Fatal error!' write ( *, '(a)' ) ' T2 <= T1 or T3 <= T2.' stop 1 end if ! ! Construct and evaluate a parabolic interpolant for the data ! in each dimension. ! y1 = ydata(left) y2 = ydata(left+1) y3 = ydata(left+2) dif1 = ( y2 - y1 ) / ( t2 - t1 ) dif2 = ( ( y3 - y1 ) / ( t3 - t1 ) & - ( y2 - y1 ) / ( t2 - t1 ) ) / ( t3 - t2 ) yval = y1 + ( tval - t1 ) * ( dif1 + ( tval - t2 ) * dif2 ) ypval = dif1 + dif2 * ( 2.0D+00 * tval - t1 - t2 ) return end !*****************************************************************************80 end module m_spline