subroutine basis_function_b_val ( tdata, tval, yval )
!
!*******************************************************************************
!
!! 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 / at TDATA(), / at TDATA() and TDATA();
! * it is for TVAL <= TDATA() or TDATA() <= TVAL;
! * it is strictly increasing from TDATA() to TDATA(),
! and strictly decreasing from TDATA() to TDATA();
! * the function and its first two derivatives are continuous
! at each node TDATA(I).
!
! Reference:
!
! Alan Davies and Philip Samuels,
! An Introduction to Computational Geometry for Curves and Surfaces,
! Clarendon Press, .
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, real TDATA(), the nodes associated with the basis function.
! The entries of TDATA are assumed to be distinct and increasing.
!
! Input, real TVAL, a point at which the B spline basis function is
! to be evaluated.
!
! Output, real YVAL, the value of the function at TVAL.
!
implicit none
!
integer, parameter :: ndata =
!
integer left
integer right
real tdata(ndata)
real tval
real u
real yval
!
if ( tval <= tdata() .or. tval >= tdata(ndata) ) then
yval = 0.0E+00
return
end if
!
! Find the interval [ TDATA(LEFT), TDATA(RIGHT) ] containing TVAL.
!
call rvec_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() ) then
yval = u** / 6.0E+00
else if ( tval < tdata() ) then
yval = ( 1.0E+00 + 3.0E+00 * u + 3.0E+00 * u** - 3.0E+00 * u** ) / 6.0E+00
else if ( tval < tdata() ) then
yval = ( 4.0E+00 - 6.0E+00 * u** + 3.0E+00 * u** ) / 6.0E+00
else if ( tval < tdata() ) then
yval = ( 1.0E+00 - u )** / 6.0E+00
end if return
end
subroutine basis_function_beta_val ( beta1, beta2, tdata, tval, yval )
!
!*******************************************************************************
!
!! BASIS_FUNCTION_BETA_VAL evaluates the beta spline basis function.
!
!
! Discussion:
!
! With BETA1 = and BETA2 = , the beta spline basis function
! equals the B spline basis function.
!
! With BETA1 large, and BETA2 = , the beta spline basis function
! skews to the right, that is, its maximum increases, and occurs
! to the right of the center.
!
! With BETA1 = 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 at the second node, at the
! third, and at the fourth.
!
! Reference:
!
! Alan Davies and Philip Samuels,
! An Introduction to Computational Geometry for Curves and Surfaces,
! Clarendon Press, , page .
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, real BETA1, the skew or bias parameter.
! BETA1 = for no skew or bias.
!
! Input, real BETA2, the tension parameter.
! BETA2 = for no tension.
!
! Input, real TDATA(), the nodes associated with the basis function.
! The entries of TDATA are assumed to be distinct and increasing.
!
! Input, real TVAL, a point at which the B spline basis function is
! to be evaluated.
!
! Output, real YVAL, the value of the function at TVAL.
!
implicit none
!
integer, parameter :: ndata =
!
real a
real b
real beta1
real beta2
real c
real d
integer left
integer right
real tdata(ndata)
real tval
real u
real yval
!
if ( tval <= tdata() .or. tval >= tdata(ndata) ) then
yval = 0.0E+00
return
end if
!
! Find the interval [ TDATA(LEFT), TDATA(RIGHT) ] containing TVAL.
!
call rvec_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() ) then yval = 2.0E+00 * u** else if ( tval < tdata() ) then a = beta2 + 4.0E+00 * beta1 + 4.0E+00 * beta1** &
+ 6.0E+00 * ( 1.0E+00 - beta1** ) &
- 3.0E+00 * ( 2.0E+00 + beta2 + 2.0E+00 * beta1 ) &
+ 2.0E+00 * ( 1.0E+00 + beta2 + beta1 + beta1** ) b = - 6.0E+00 * ( 1.0E+00 - beta1** ) &
+ 6.0E+00 * ( 2.0E+00 + beta2 + 2.0E+00 * beta1 ) &
- 6.0E+00 * ( 1.0E+00 + beta2 + beta1 + beta1** ) c = - 3.0E+00 * ( 2.0E+00 + beta2 + 2.0E+00 * beta1 ) &
+ 6.0E+00 * ( 1.0E+00 + beta2 + beta1 + beta1** ) d = - 2.0E+00 * ( 1.0E+00 + beta2 + beta1 + beta1** ) yval = a + b * u + c * u** + d * u** else if ( tval < tdata() ) then a = beta2 + 4.0E+00 * beta1 + 4.0E+00 * beta1** b = - 6.0E+00 * ( beta1 - beta1** ) c = - 3.0E+00 * ( beta2 + 2.0E+00 * beta1** + 2.0E+00 * beta1** ) d = 2.0E+00 * ( beta2 + beta1 + beta1** + beta1** ) yval = a + b * u + c * u** + d * u** else if ( tval < tdata() ) then yval = 2.0E+00 * beta1** * ( 1.0E+00 - u )** end if yval = yval / ( 2.0E+00 + beta2 + 4.0E+00 * beta1 + 4.0E+00 * beta1** &
+ 2.0E+00 * beta1** ) return
end
subroutine basis_matrix_b_uni ( mbasis )
!
!*******************************************************************************
!
!! BASIS_MATRIX_B_UNI sets up the uniform B spline basis matrix.
!
!
! Reference:
!
! Foley, van Dam, Feiner, Hughes,
! Computer Graphics: Principles and Practice,
! page .
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Output, real MBASIS(,), the basis matrix.
!
implicit none
!
real mbasis(,)
!
mbasis(,) = - 1.0E+00 / 6.0E+00
mbasis(,) = 3.0E+00 / 6.0E+00
mbasis(,) = - 3.0E+00 / 6.0E+00
mbasis(,) = 1.0E+00 / 6.0E+00 mbasis(,) = 3.0E+00 / 6.0E+00
mbasis(,) = - 6.0E+00 / 6.0E+00
mbasis(,) = 3.0E+00 / 6.0E+00
mbasis(,) = 0.0E+00 mbasis(,) = - 3.0E+00 / 6.0E+00
mbasis(,) = 0.0E+00
mbasis(,) = 3.0E+00 / 6.0E+00
mbasis(,) = 0.0E+00 mbasis(,) = 1.0E+00 / 6.0E+00
mbasis(,) = 4.0E+00 / 6.0E+00
mbasis(,) = 1.0E+00 / 6.0E+00
mbasis(,) = 0.0E+00 return
end
subroutine basis_matrix_beta_uni ( beta1, beta2, mbasis )
!
!*******************************************************************************
!
!! BASIS_MATRIX_BETA_UNI sets up the uniform beta spline basis matrix.
!
!
! Discussion:
!
! If BETA1 = and BETA2 = , then the beta spline reduces to
! the B spline.
!
! Reference:
!
! Foley, van Dam, Feiner, Hughes,
! Computer Graphics: Principles and Practice,
! page .
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, real BETA1, the skew or bias parameter.
! BETA1 = for no skew or bias.
!
! Input, real BETA2, the tension parameter.
! BETA2 = for no tension.
!
! Output, real MBASIS(,), the basis matrix.
!
implicit none
!
real beta1
real beta2
real delta
integer i
integer j
real mbasis(,)
!
mbasis(,) = - 2.0E+00 * beta1**
mbasis(,) = 2.0E+00 * ( beta2 + beta1** + beta1** + beta1 )
mbasis(,) = - 2.0E+00 * ( beta2 + beta1** + beta1 + 1.0E+00 )
mbasis(,) = 2.0E+00 mbasis(,) = 6.0E+00 * beta1**
mbasis(,) = - 3.0E+00 * ( beta2 + 2.0E+00 * beta1** + 2.0E+00 * beta1** )
mbasis(,) = 3.0E+00 * ( beta2 + 2.0E+00 * beta1** )
mbasis(,) = 0.0E+00 mbasis(,) = - 6.0E+00 * beta1**
mbasis(,) = 6.0E+00 * beta1 * ( beta1 - 1.0E+00 ) * ( beta1 + 1.0E+00 )
mbasis(,) = 6.0E+00 * beta1
mbasis(,) = 0.0E+00 mbasis(,) = 2.0E+00 * beta1**
mbasis(,) = 4.0E+00 * beta1 * ( beta1 + 1.0E+00 ) + beta2
mbasis(,) = 2.0E+00
mbasis(,) = 0.0E+00 delta = beta2 + 2.0E+00 * beta1** + 4.0E+00 * beta1** &
+ 4.0E+00 * beta1 + 2.0E+00 mbasis(:,:) = mbasis(:,:) / delta return
end
subroutine basis_matrix_bezier ( mbasis )
!
!*******************************************************************************
!
!! BASIS_MATRIX_BEZIER_UNI 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 = , while
! P2 is used to approximate the derivative at T = by
! dP/dt = * ( P2 - P1 ). Similarly, P4 is the function value
! at T = , and P3 is used to approximate the derivative at T =
! by dP/dT = * ( P4 - P3 ).
!
! Reference:
!
! Foley, van Dam, Feiner, Hughes,
! Computer Graphics: Principles and Practice,
! page .
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Output, real MBASIS(,), the basis matrix.
!
implicit none
!
real mbasis(,)
!
mbasis(,) = - 1.0E+00
mbasis(,) = 3.0E+00
mbasis(,) = - 3.0E+00
mbasis(,) = 1.0E+00 mbasis(,) = 3.0E+00
mbasis(,) = - 6.0E+00
mbasis(,) = 3.0E+00
mbasis(,) = 0.0E+00 mbasis(,) = - 3.0E+00
mbasis(,) = 3.0E+00
mbasis(,) = 0.0E+00
mbasis(,) = 0.0E+00 mbasis(,) = 1.0E+00
mbasis(,) = 0.0E+00
mbasis(,) = 0.0E+00
mbasis(,) = 0.0E+00 return
end
subroutine basis_matrix_hermite ( mbasis )
!
!*******************************************************************************
!
!! 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 = , while P2 and P2' apply at T = 1.
!
! Reference:
!
! Foley, van Dam, Feiner, Hughes,
! Computer Graphics: Principles and Practice,
! page .
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Output, real MBASIS(,), the basis matrix.
!
implicit none
!
real mbasis(,)
!
mbasis(,) = 2.0E+00
mbasis(,) = - 2.0E+00
mbasis(,) = 1.0E+00
mbasis(,) = 1.0E+00 mbasis(,) = - 3.0E+00
mbasis(,) = 3.0E+00
mbasis(,) = - 2.0E+00
mbasis(,) = - 1.0E+00 mbasis(,) = 0.0E+00
mbasis(,) = 0.0E+00
mbasis(,) = 1.0E+00
mbasis(,) = 0.0E+00 mbasis(,) = 1.0E+00
mbasis(,) = 0.0E+00
mbasis(,) = 0.0E+00
mbasis(,) = 0.0E+00 return
end
subroutine basis_matrix_overhauser_nonuni ( alpha, beta, mbasis )
!
!*******************************************************************************
!
!! BASIS_MATRIX_OVERHAUSER_NONUNI sets up the 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 = ,
! and P3 to T = .
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, real ALPHA, BETA.
! ALPHA = | P2 - P1 | / ( | P3 - P2 | + | P2 - P1 | )
! BETA = | P3 - P2 | / ( | P4 - P3 | + | P3 - P2 | ).
!
! Output, real MBASIS(,), the basis matrix.
!
implicit none
!
real alpha
real beta
real mbasis(,)
!
mbasis(,) = - ( 1.0E+00 - alpha )** / alpha
mbasis(,) = beta + ( 1.0E+00 - alpha ) / alpha
mbasis(,) = alpha - 1.0E+00 / ( 1.0E+00 - beta )
mbasis(,) = beta** / ( 1.0E+00 - beta ) mbasis(,) = 2.0E+00 * ( 1.0E+00 - alpha )** / alpha
mbasis(,) = ( - 2.0E+00 * ( 1.0E+00 - alpha ) - alpha * beta ) / alpha
mbasis(,) = ( 2.0E+00 * ( 1.0E+00 - alpha ) &
- beta * ( 1.0E+00 - 2.0E+00 * alpha ) ) / ( 1.0E+00 - beta )
mbasis(,) = - beta** / ( 1.0E+00 - beta ) mbasis(,) = - ( 1.0E+00 - alpha )** / alpha
mbasis(,) = ( 1.0E+00 - 2.0E+00 * alpha ) / alpha
mbasis(,) = alpha
mbasis(,) = 0.0E+00 mbasis(,) = 0.0E+00
mbasis(,) = 1.0E+00
mbasis(,) = 0.0E+00
mbasis(,) = 0.0E+00 return
end
subroutine basis_matrix_overhauser_nul ( alpha, mbasis )
!
!*******************************************************************************
!
!! BASIS_MATRIX_OVERHAUSER_NUL sets up the left nonuniform Overhauser spline 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 = ,
! and P2 to T = . (???)
!
! Modified:
!
! August
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, real ALPHA.
! ALPHA = | P2 - P1 | / ( | P3 - P2 | + | P2 - P1 | )
!
! Output, real MBASIS(,), the basis matrix.
!
implicit none
!
real alpha
real mbasis(,)
!
mbasis(,) = 1.0E+00 / alpha
mbasis(,) = - 1.0E+00 / ( alpha * ( 1.0E+00 - alpha ) )
mbasis(,) = 1.0E+00 / ( 1.0E+00 - alpha ) mbasis(,) = - ( 1.0E+00 + alpha ) / alpha
mbasis(,) = 1.0E+00 / ( alpha * ( 1.0E+00 - alpha ) )
mbasis(,) = - alpha / ( 1.0E+00 - alpha ) mbasis(,) = 1.0E+00
mbasis(,) = 0.0E+00
mbasis(,) = 0.0E+00 return
end
subroutine basis_matrix_overhauser_nur ( beta, mbasis )
!
!*******************************************************************************
!
!! BASIS_MATRIX_OVERHAUSER_NUR sets up the right nonuniform Overhauser spline basis matrix.
!
!
! Discussion:
!
! This basis matrix assumes that the data points PN-, PN-, and
! PN are not uniformly spaced in T, and that PN- corresponds to T = ,
! and PN to T = . (???)
!
! Modified:
!
! August
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, real BETA.
! BETA = | P(N) - P(N-) | / ( | P(N) - P(N-) | + | P(N-) - P(N-) | )
!
! Output, real MBASIS(,), the basis matrix.
!
implicit none
!
real beta
real mbasis(,)
!
mbasis(,) = 1.0E+00 / beta
mbasis(,) = - 1.0E+00 / ( beta * ( 1.0E+00 - beta ) )
mbasis(,) = 1.0E+00 / ( 1.0E+00 - beta ) mbasis(,) = - ( 1.0E+00 + beta ) / beta
mbasis(,) = 1.0E+00 / ( beta * ( 1.0E+00 - beta ) )
mbasis(,) = - beta / ( 1.0E+00 - beta ) mbasis(,) = 1.0E+00
mbasis(,) = 0.0E+00
mbasis(,) = 0.0E+00 return
end
subroutine basis_matrix_overhauser_uni ( mbasis )
!
!*******************************************************************************
!
!! BASIS_MATRIX_OVERHAUSER_UNI sets up 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 = ,
! and P3 to T = .
!
! Reference:
!
! Foley, van Dam, Feiner, Hughes,
! Computer Graphics: Principles and Practice,
! page .
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Output, real MBASIS(,), the basis matrix.
!
implicit none
!
real mbasis(,)
!
mbasis(,) = - 1.0E+00 / 2.0E+00
mbasis(,) = 3.0E+00 / 2.0E+00
mbasis(,) = - 3.0E+00 / 2.0E+00
mbasis(,) = 1.0E+00 / 2.0E+00 mbasis(,) = 2.0E+00 / 2.0E+00
mbasis(,) = - 5.0E+00 / 2.0E+00
mbasis(,) = 4.0E+00 / 2.0E+00
mbasis(,) = - 1.0E+00 / 2.0E+00 mbasis(,) = - 1.0E+00 / 2.0E+00
mbasis(,) = 0.0E+00
mbasis(,) = 1.0E+00 / 2.0E+00
mbasis(,) = 0.0E+00 mbasis(,) = 0.0E+00
mbasis(,) = 2.0E+00 / 2.0E+00
mbasis(,) = 0.0E+00
mbasis(,) = 0.0E+00 return
end
subroutine basis_matrix_overhauser_uni_l ( mbasis )
!
!*******************************************************************************
!
!! BASIS_MATRIX_OVERHAUSER_UNI_L sets up the left uniform Overhauser spline 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 = ,
! and P2 to T = .
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Output, real MBASIS(,), the basis matrix.
!
implicit none
!
real mbasis(,)
!
mbasis(,) = 2.0E+00
mbasis(,) = - 4.0E+00
mbasis(,) = 2.0E+00 mbasis(,) = - 3.0E+00
mbasis(,) = 4.0E+00
mbasis(,) = - 1.0E+00 mbasis(,) = 1.0E+00
mbasis(,) = 0.0E+00
mbasis(,) = 0.0E+00 return
end
subroutine basis_matrix_overhauser_uni_r ( mbasis )
!
!*******************************************************************************
!
!! BASIS_MATRIX_OVERHAUSER_UNI_R sets up the right uniform Overhauser spline basis matrix.
!
!
! Discussion:
!
! This basis matrix assumes that the data points P(N-), P(N-),
! and P(N) are uniformly spaced in T, and that P(N-) corresponds to
! T = , and P(N) to T = .
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Output, real MBASIS(,), the basis matrix.
!
implicit none
!
real mbasis(,)
!
mbasis(,) = 2.0E+00
mbasis(,) = - 4.0E+00
mbasis(,) = 2.0E+00 mbasis(,) = - 3.0E+00
mbasis(,) = 4.0E+00
mbasis(,) = - 1.0E+00 mbasis(,) = 1.0E+00
mbasis(,) = 0.0E+00
mbasis(,) = 0.0E+00 return
end
subroutine basis_matrix_tmp ( left, n, mbasis, ndata, tdata, ydata, tval, yval )
!
!*******************************************************************************
!
!! 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.
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer LEFT, indicats that TVAL is in the interval
! [ TDATA(LEFT), TDATA(LEFT+) ], or that this is the "nearest"
! interval to TVAL.
! For TVAL < TDATA(), use LEFT = .
! For TVAL > TDATA(NDATA), use LEFT = NDATA - .
!
! Input, integer N, the order of the basis matrix.
!
! Input, real MBASIS(N,N), the basis matrix.
!
! Input, integer NDATA, the dimension of the vectors TDATA and YDATA.
!
! Input, real TDATA(NDATA), the abscissa values. This routine
! assumes that the TDATA values are uniformly spaced, with an
! increment of 1.0.
!
! Input, real YDATA(NDATA), the data values to be interpolated or
! approximated.
!
! Input, real TVAL, the value of T at which the spline is to be
! evaluated.
!
! Output, real YVAL, the value of the spline at TVAL.
!
implicit none
!
integer, parameter :: maxn =
integer n
integer ndata
!
real arg
integer first
integer i
integer j
integer left
real mbasis(n,n)
real tdata(ndata)
real temp
real tval
real tvec(maxn)
real ydata(ndata)
real yval
!
if ( left == ) then
arg = 0.5 * ( tval - tdata(left) )
first = left
else if ( left < ndata - ) then
arg = tval - tdata(left)
first = left -
else if ( left == ndata - ) then
arg = 0.5E+00 * ( 1.0E+00 + tval - tdata(left) )
first = left -
end if do i = , n -
tvec(i) = arg**( n - i )
end do
tvec(n) = 1.0E+00 yval = 0.0E+00
do j = , n
yval = yval + dot_product ( tvec(:n), mbasis(:n,j) ) &
* ydata(first - + j)
end do return
end
subroutine bc_val ( n, t, xcon, ycon, xval, yval )
!
!*******************************************************************************
!
!! BC_VAL evaluates a parameterized Bezier curve.
!
!
! 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 (i = to N) XCON(I) * BERN(I,N)(T)
! Y(T) = Sum (i = to N) YCON(I) * BERN(I,N)(T)
!
! where BERN(I,N)(T) is the I-th Bernstein polynomial of order N
! defined on the interval [,], and where XCON(*) and YCON(*)
! are the coordinates of N+ "control points".
!
! Modified:
!
! March
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer N, the order of the Bezier curve, which
! must be at least .
!
! Input, real T, the point at which the Bezier curve should
! be evaluated. The best results are obtained within the interval
! [,] but T may be anywhere.
!
! Input, real XCON(:N), YCON(:N), the X and Y coordinates
! of the control points. The Bezier curve will pass through
! the points ( XCON(), YCON() ) and ( XCON(N), YCON(N) ), but
! generally NOT through the other control points.
!
! Output, real XVAL, YVAL, the X and Y coordinates of the point
! on the Bezier curve corresponding to the given T value.
!
implicit none
!
integer n
!
real bval(:n)
integer i
real t
real xcon(:n)
real xval
real ycon(:n)
real yval
!
call bp01 ( n, bval, t ) xval = dot_product ( xcon(:n), bval(:n) )
yval = dot_product ( ycon(:n), bval(:n) ) return
end
function bez_val ( n, x, a, b, y )
!
!*******************************************************************************
!
!! BEZ_VAL evaluates a Bezier function at a point.
!
!
! Discussion:
!
! The Bezier function has the form:
!
! BEZ(X) = Sum (i = to N) Y(I) * BERN(N,I)( (X-A)/(B-A) )
!
! where BERN(N,I)(X) is the I-th Bernstein polynomial of order N
! defined on the interval [,], and Y(*) is a set of
! coefficients.
!
! If we define the N+ equally spaced
!
! X(I) = ( (N-I)*A + I*B)/ N,
!
! for I = to N, then the pairs ( X(I), Y(I) ) can be regarded as
! "control points".
!
! Modified:
!
! March
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer N, the order of the Bezier function, which
! must be at least .
!
! Input, real 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 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() 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 Y(:N), a set of data defining the Y coordinates
! of the control points.
!
! Output, real BEZ_VAL, the value of the Bezier function at X.
!
implicit none
!
integer n
!
real a
real b
real bez_val
real bval(:n)
integer i
integer ncopy
real x
real x01
real y(:n)
!
if ( b - a == 0.0E+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'BEZ_VAL - Fatal error!'
write ( *, '(a,g14.6)' ) ' Null interval, A = B = ', a
stop
end if x01 = ( x - a ) / ( b - a ) ncopy = n
call bp01 ( ncopy, bval, x01 ) bez_val = dot_product ( y(:n), bval(:n) ) return
end
subroutine bp01 ( n, bern, x )
!
!*******************************************************************************
!
!! BP01 evaluates the Bernstein basis polynomials for [,] at a point X.
!
!
! Formula:
!
! BERN(N,I,X) = [N!/(I!*(N-I)!)] * (-X)**(N-I) * X**I
!
! First values:
!
! B(,,X) =
!
! B(,,X) = -X
! B(,,X) = X
!
! B(,,X) = (-X)**
! B(,,X) = * (-X) * X
! B(,,X) = X**
!
! B(,,X) = (-X)**
! B(,,X) = * (-X)** * X
! B(,,X) = * (-X) * X**
! B(,,X) = X**
!
! B(,,X) = (-X)**
! B(,,X) = * (-X)** * X
! B(,,X) = * (-X)** * X**
! B(,,X) = * (-X) * X**
! B(,,X) = X**
!
! Special values:
!
! B(N,I,/) = C(N,K) / **N
!
! Modified:
!
! January
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer N, the degree of the Bernstein basis polynomials.
! For any N greater than or equal to , there is a set of N+ Bernstein
! basis polynomials, each of degree N, which form a basis for
! all polynomials on [,].
!
! Output, real BERN(:N), the values of the N+ Bernstein basis
! polynomials at X.
!
! Input, real X, the point at which the polynomials are to be
! evaluated.
!
implicit none
!
integer n
!
real bern(:n)
integer i
integer j
real x
!
if ( n == ) then bern() = 1.0E+00 else if ( n > ) then bern() = 1.0E+00 - x
bern() = x do i = , n
bern(i) = x * bern(i-)
do j = i-, , -
bern(j) = x * bern(j-) + ( 1.0E+00 - x ) * bern(j)
end do
bern() = ( 1.0E+00 - x ) * bern()
end do end if return
end
subroutine bp_approx ( n, a, b, ydata, xval, yval )
!
!*******************************************************************************
!
!! BP_APPROX evaluates the Bernstein polynomial for F(X) on [A,B].
!
!
! Formula:
!
! BERN(F)(X) = SUM ( <= 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.
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer N, the degree of the Bernstein polynomial to be used.
!
! Input, real A, B, the endpoints of the interval on which the
! approximant is based. A and B should not be equal.
!
! Input, real YDATA(:N), the data values at N+ equally spaced points
! in [A,B]. If N = , then the evaluation point should be 0.5 * ( A + B).
! Otherwise, evaluation point I should be ( (N-I)*A + I*B ) / N ).
!
! Input, real 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 YVAL, the value of the Bernstein polynomial approximant
! for F, based in [A,B], evaluated at XVAL.
!
implicit none
!
integer n
!
real a
real b
real bvec(:n)
integer i
real xval
real ydata(:n)
real yval
!
! Evaluate the Bernstein basis polynomials at XVAL.
!
call bpab ( n, bvec, xval, a, b )
!
! Now compute the sum of YDATA(I) * BVEC(I).
!
yval = dot_product ( ydata(:n), bvec(:n) ) return
end
subroutine bpab ( n, bern, x, a, b )
!
!*******************************************************************************
!
!! BPAB evaluates the Bernstein basis polynomials for [A,B] at a point X.
!
!
! Formula:
!
! BERN(N,I,X) = [N!/(I!*(N-I)!)] * (B-X)**(N-I) * (X-A)**I / (B-A)**N
!
! First values:
!
! B(,,X) =
!
! B(,,X) = ( B-X ) / (B-A)
! B(,,X) = ( X-A ) / (B-A)
!
! B(,,X) = ( (B-X)** ) / (B-A)**
! B(,,X) = ( * (B-X) * (X-A) ) / (B-A)**
! B(,,X) = ( (X-A)** ) / (B-A)**
!
! B(,,X) = ( (B-X)** ) / (B-A)**
! B(,,X) = ( * (B-X)** * (X-A) ) / (B-A)**
! B(,,X) = ( * (B-X) * (X-A)** ) / (B-A)**
! B(,,X) = ( (X-A)** ) / (B-A)**
!
! B(,,X) = ( (B-X)** ) / (B-A)**
! B(,,X) = ( * (B-X)** * (X-A) ) / (B-A)**
! B(,,X) = ( * (B-X)** * (X-A)** ) / (B-A)**
! B(,,X) = ( * (B-X) * (X-A)** ) / (B-A)**
! B(,,X) = ( (X-A)** ) / (B-A)**
!
! Modified:
!
! January
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer N, the degree of the Bernstein basis polynomials.
! For any N greater than or equal to , there is a set of N+
! Bernstein basis polynomials, each of degree N, which form a basis
! for polynomials on [A,B].
!
! Output, real BERN(:N), the values of the N+ Bernstein basis
! polynomials at X.
!
! Input, real X, the point at which the polynomials are to be
! evaluated. X need not lie in the interval [A,B].
!
! Input, real A, B, the endpoints of the interval on which the
! polynomials are to be based. A and B should not be equal.
!
implicit none
!
integer n
!
real a
real b
real bern(:n)
integer i
integer j
real x
!
if ( b == a ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'BPAB - Fatal error!'
write ( *, '(a,g14.6)' ) ' A = B = ', a
stop
end if if ( n == ) then bern() = 1.0E+00 else if ( n > ) then bern() = ( b - x ) / ( b - a )
bern() = ( x - a ) / ( b - a ) do i = , n
bern(i) = ( x - a ) * bern(i-) / ( b - a )
do j = i-, , -
bern(j) = ( ( b - x ) * bern(j) + ( x - a ) * bern(j-) ) / ( b - a )
end do
bern() = ( b - x ) * bern() / ( b - a )
end do end if return
end
subroutine data_to_dif ( diftab, ntab, xtab, ytab )
!
!*******************************************************************************
!
!! 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.
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Output, real DIFTAB(NTAB), the divided difference coefficients
! corresponding to the input (XTAB,YTAB).
!
! Input, integer 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 XTAB(NTAB), the X values at which data was taken.
! These values must be distinct.
!
! Input, real YTAB(NTAB), the corresponding Y values.
!
implicit none
!
integer ntab
!
real diftab(ntab)
integer i
integer j
logical rvec_distinct
real xtab(ntab)
real ytab(ntab)
!
if ( .not. rvec_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(:ntab) = ytab(:ntab)
!
! Compute the divided differences.
!
do i = , ntab
do j = ntab, i, - diftab(j) = ( diftab(j) - diftab(j-) ) / ( xtab(j) - xtab(j+-i) ) end do
end do return
end
subroutine dif_val ( diftab, ntab, xtab, xval, yval )
!
!*******************************************************************************
!
!! DIF_VAL evaluates a divided difference polynomial at a point.
!
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, real DIFTAB(NTAB), the divided difference polynomial coefficients.
!
! Input, integer NTAB, the number of divided difference
! coefficients in DIFTAB, and the number of points XTAB.
!
! Input, real XTAB(NTAB), the X values upon which the
! divided difference polynomial is based.
!
! Input, real XVAL, a value of X at which the polynomial
! is to be evaluated.
!
! Output, real YVAL, the value of the polynomial at XVAL.
!
implicit none
!
integer ntab
!
real diftab(ntab)
integer i
real xtab(ntab)
real xval
real yval
!
yval = diftab(ntab)
do i = , ntab-
yval = diftab(ntab-i) + ( xval - xtab(ntab-i) ) * yval
end do return
end
subroutine least_set ( ntab, xtab, ytab, ndeg, ptab, array, eps, ierror )
!
!*******************************************************************************
!
!! LEAST_SET constructs the least squares polynomial approximation to data.
!
!
! Discussion:
!
! The routine LEAST_EVAL must be used to evaluate the approximation at a
! point.
!
! Modified:
!
! November
!
! Parameters:
!
! Input, integer NTAB, the number of data points.
!
! Input, real XTAB(NTAB), the X data. The values in XTAB
! should be distinct, and in increasing order.
!
! Input, real YTAB(NTAB), the Y data values corresponding
! to the X data in XTAB.
!
! Input, integer NDEG, the degree of the polynomial which the
! program is to use. NDEG must be at least , and less than or
! equal to NTAB-.
!
! Output, real PTAB(NTAB). PTAB(I) is the value of the
! least squares polynomial at the point XTAB(I).
!
! Output, real ARRAY(*NTAB+*NDEG), an array containing data about
! the polynomial.
!
! Output, real EPS, the root-mean-square discrepancy of the
! polynomial fit.
!
! Output, integer IERROR, error flag.
! zero, no error occurred;
! nonzero, an error occurred, and the polynomial could not be computed.
!
implicit none
!
integer ndeg
integer ntab
!
real array(*ntab+*ndeg)
real eps
real error
integer i
integer i0l1
integer i1l1
integer ierror
integer it
integer k
integer mdeg
real ptab(ntab)
real rn0
real rn1
real s
real sum2
real xtab(ntab)
real y_sum
real ytab(ntab)
!
ierror =
!
! Check NDEG.
!
if ( ndeg < ) then
ierror =
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'LEAST_SET - Fatal error!'
write ( *, '(a)' ) ' NDEG < 1.'
return
else if ( ndeg >= ntab ) then
ierror =
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'LEAST_SET - Fatal error!'
write ( *, '(a)' ) ' NDEG >= NTAB.'
return
end if
!
! Check that the abscissas are strictly increasing.
!
do i = , ntab-
if ( xtab(i) >= xtab(i+) ) then
ierror =
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'LEAST_SET - Fatal error!'
write ( *, '(a)' ) ' XTAB must be strictly increasing, but'
write ( *, '(a,i6,a,g14.6)' ) ' XTAB(', i, ') = ', xtab(i)
write ( *, '(a,i6,a,g14.6)' ) ' XTAB(', i+, ') = ', xtab(i+)
return
end if
end do i0l1 = * ndeg
i1l1 = * ndeg + ntab y_sum = sum ( ytab )
rn0 = ntab
array(*ndeg) = y_sum / real ( ntab ) ptab(:ntab) = y_sum / real ( ntab ) error = 0.0E+00
do i = , ntab
error = error + ( y_sum / real ( ntab ) - ytab(i) )**
end do if ( ndeg == ) then
eps = sqrt ( error / real ( ntab ) )
return
end if array() = sum ( xtab ) / real ( ntab ) s = 0.0E+00
sum2 = 0.0E+00
do i = , ntab
array(i1l1+i) = xtab(i) - array()
s = s + array(i1l1+i)**
sum2 = sum2 + array(i1l1+i) * ( ytab(i) - ptab(i) )
end do rn1 = s
array(*ndeg+) = sum2 / s do i = , ntab
ptab(i) = ptab(i) + sum2 * array(i1l1+i) / s
end do error = sum ( ( ptab(:ntab) - ytab(:ntab) )** ) if ( ndeg == ) then
eps = sqrt ( error / real ( ntab ) )
return
end if do i = , ntab
array(*ndeg+i) = 1.0E+00
end do mdeg =
k = do array(ndeg-+k) = rn1 / rn0 sum2 = 0.0E+00
do i = , ntab
sum2 = sum2 + xtab(i) * array(i1l1+i)**
end do array(k) = sum2 / rn1 s = 0.0E+00
sum2 = 0.0E+00
do i = , ntab
array(i0l1+i) = ( xtab(i) - array(k) ) * array(i1l1+i) &
- array(ndeg-+k) * array(i0l1+i)
s = s + array(i0l1+i)**
sum2 = sum2 + array(i0l1+i) * ( ytab(i) - ptab(i) )
end do rn0 = rn1
rn1 = s
it = i0l1
i0l1 = i1l1
i1l1 = it
array(*ndeg+k) = sum2 / rn1 do i = , ntab
ptab(i) = ptab(i) + sum2 * array(i1l1+i) / rn1
end do error = sum ( ( ptab(:ntab) - ytab(:ntab) )** ) if ( mdeg >= ndeg ) then
exit
end if mdeg = mdeg +
k = k + end do eps = sqrt ( error / real ( ntab ) ) return
end
subroutine least_val ( x, ndeg, array, value )
!
!*******************************************************************************
!
!! LEAST_VAL evaluates a least squares polynomial defined by LEAST_SET.
!
!
! Modified:
!
! March
!
! Parameters:
!
! Input, real X, the point at which the polynomial is to be evaluated.
!
! Input, integer NDEG, the degree of the polynomial fit used.
! This is the value of NDEG as returned from LEAST_SET.
!
! Input, real ARRAY(*), an array of a certain dimension.
! See LEAST_SET for details on the size of ARRAY.
! ARRAY contains information about the polynomial, as set up by LEAST_SET.
!
! Output, real VALUE, the value of the polynomial at X.
!
implicit none
!
real array(*)
real dk
real dkp1
real dkp2
integer k
integer l
integer ndeg
real value
real x
!
if ( ndeg <= ) then value = array(*ndeg) else if ( ndeg == ) then value = array(*ndeg) + array(*ndeg+) * ( x - array() ) else dkp2 = array(*ndeg)
dkp1 = array(*ndeg-) + ( x - array(ndeg) ) * dkp2 do l = , ndeg- k = ndeg - - l dk = array(*ndeg+k) + ( x - array(k+) ) * dkp1 - array(ndeg++k) * dkp2 dkp2 = dkp1 dkp1 = dk end do value = array(*ndeg) + ( x - array() ) * dkp1 - array(ndeg+) * dkp2 end if return
end
subroutine parabola_val2 ( ndim, ndata, tdata, ydata, left, tval, yval )
!
!*******************************************************************************
!
!! 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
! consecutive entries of a table and evaluates this interpolant
! at a given abscissa value.
!
! Modified:
!
! March
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer NDIM, the dimension of a single data point.
! NDIM must be at least .
!
! Input, integer NDATA, the number of data points.
! NDATA must be at least .
!
! Input, real TDATA(NDATA), the abscissas of the data points. The
! values in TDATA must be in strictly ascending order.
!
! Input, real YDATA(NDIM,NDATA), the data points corresponding to
! the abscissas.
!
! Input, integer LEFT, the location of the first of the three
! consecutive data points through which the parabolic interpolant
! must pass. <= LEFT <= NDATA - .
!
! Input, real TVAL, the value of T at which the parabolic interpolant
! is to be evaluated. Normally, TDATA() <= TVAL <= T(NDATA), and
! the data will be interpolated. For TVAL outside this range,
! extrapolation will be used.
!
! Output, real YVAL(NDIM), the value of the parabolic interpolant at TVAL.
!
implicit none
!
integer ndata
integer ndim
!
real dif1
real dif2
integer i
integer left
real t1
real t2
real t3
real tval
real tdata(ndata)
real ydata(ndim,ndata)
real y1
real y2
real y3
real yval(ndim)
!
! Check.
!
if ( left < .or. left > ndata- ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'PARABOLA_VAL2 - Fatal error!'
write ( *, '(a)' ) ' LEFT < 1 or LEFT > NDATA-2.'
stop
end if if ( ndim < ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'PARABOLA_VAL2 - Fatal error!'
write ( *, '(a)' ) ' NDIM < 1.'
stop
end if
!
! Copy out the three abscissas.
!
t1 = tdata(left)
t2 = tdata(left+)
t3 = tdata(left+) if ( t1 >= t2 .or. t2 >= t3 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'PARABOLA_VAL2 - Fatal error!'
write ( *, '(a)' ) ' T1 >= T2 or T2 >= T3.'
stop
end if
!
! Construct and evaluate a parabolic interpolant for the data
! in each dimension.
!
do i = , ndim y1 = ydata(i,left)
y2 = ydata(i,left+)
y3 = ydata(i,left+) 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
subroutine r_random ( rlo, rhi, r )
!
!*******************************************************************************
!
!! R_RANDOM returns a random real in a given range.
!
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, real RLO, RHI, the minimum and maximum values.
!
! Output, real R, the randomly chosen value.
!
implicit none
!
real r
real rhi
real rlo
real t
!
! Pick T, a random number in (,).
!
call random_number ( harvest = t )
!
! Set R in ( RLO, RHI ).
!
r = ( 1.0E+00 - t ) * rlo + t * rhi return
end
subroutine r_swap ( x, y )
!
!*******************************************************************************
!
!! R_SWAP swaps two real values.
!
!
! Modified:
!
! May
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input/output, real X, Y. On output, the values of X and
! Y have been interchanged.
!
implicit none
!
real x
real y
real z
!
z = x
x = y
y = z return
end
subroutine rvec_bracket ( n, x, xval, left, right )
!
!*******************************************************************************
!
!! RVEC_BRACKET searches a sorted array for successive brackets of a value.
!
!
! Discussion:
!
! 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.
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer N, length of input array.
!
! Input, real X(N), an array sorted into ascending order.
!
! Input, real XVAL, a value to be bracketed.
!
! Output, integer LEFT, RIGHT, the results of the search.
! Either:
! XVAL < X(), when LEFT = , RIGHT = ;
! XVAL > X(N), when LEFT = N-, RIGHT = N;
! or
! X(LEFT) <= XVAL <= X(RIGHT).
!
implicit none
!
integer n
!
integer i
integer left
integer right
real x(n)
real xval
!
do i = , n - if ( xval < x(i) ) then
left = i -
right = i
return
end if end do left = n -
right = n return
end
subroutine rvec_bracket3 ( n, t, tval, left )
!
!*******************************************************************************
!
!! RVEC_BRACKET3 finds the interval containing or nearest a given value.
!
!
! Discussion:
!
! 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+) ], or
! * T < T(LEFT) = T(), or
! * T > T(LEFT+) = 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.
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer N, length of the input array.
!
! Input, real T(N), an array sorted into ascending order.
!
! Input, real TVAL, a value to be bracketed by entries of T.
!
! Input/output, integer LEFT.
!
! On input, if <= LEFT <= N-, LEFT is taken as a suggestion for the
! interval [ T(LEFT), T(LEFT+) ] 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+) ]
! is the closest to TVAL; it either contains TVAL, or else TVAL
! lies outside the interval [ T(), T(N) ].
!
implicit none
!
integer n
!
integer high
integer left
integer low
integer mid
real t(n)
real tval
!
! Check the input data.
!
if ( n < ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'RVEC_BRACKET3 - Fatal error!'
write ( *, '(a)' ) ' N must be at least 2.'
stop
end if
!
! If LEFT is not between and N-, set it to the middle value.
!
if ( left < .or. left > n - ) then
left = ( n + ) /
end if
!
! CASE : TVAL < T(LEFT):
! Search for TVAL in [T(I), T(I+)] for intervals I = to LEFT-.
!
if ( tval < t(left) ) then if ( left == ) then
return
else if ( left == ) then
left =
return
else if ( tval >= t(left-) ) then
left = left -
return
else if ( tval <= t() ) then
left =
return
end if
!
! ...Binary search for TVAL in [T(I), T(I+)] for intervals I = to LEFT-.
!
low =
high = left - do if ( low == high ) then
left = low
return
end if mid = ( low + high + ) / if ( tval >= t(mid) ) then
low = mid
else
high = mid -
end if end do
!
! CASE2: T(LEFT+) < TVAL:
! Search for TVAL in {T(I),T(I+)] for intervals I = LEFT+ to N-.
!
else if ( tval > t(left+) ) then if ( left == n - ) then
return
else if ( left == n - ) then
left = left +
return
else if ( tval <= t(left+) ) then
left = left +
return
else if ( tval >= t(n-) ) then
left = n -
return
end if
!
! ...Binary search for TVAL in [T(I), T(I+)] for intervals I = LEFT+ to N-.
!
low = left +
high = n - do if ( low == high ) then
left = low
return
end if mid = ( low + high + ) / if ( tval >= t(mid) ) then
low = mid
else
high = mid -
end if end do
!
! CASE3: T(LEFT) <= TVAL <= T(LEFT+):
! T is in [T(LEFT), T(LEFT+)], as the user said it might be.
!
else end if return
end
function rvec_distinct ( n, x )
!
!*******************************************************************************
!
!! RVEC_DISTINCT is true if the entries in a real vector are distinct.
!
!
! Modified:
!
! September
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer N, the number of entries in the vector.
!
! Input, real X(N), the vector to be checked.
!
! Output, logical RVEC_DISTINCT is .TRUE. if all N elements of X
! are distinct.
!
implicit none
!
integer n
!
integer i
integer j
logical rvec_distinct
real x(n)
!
rvec_distinct = .false. do i = , n
do j = , i -
if ( x(i) == x(j) ) then
return
end if
end do
end do rvec_distinct = .true. return
end
subroutine rvec_even ( alo, ahi, n, a )
!
!*******************************************************************************
!
!! RVEC_EVEN returns N real values, evenly spaced between ALO and AHI.
!
!
! Modified:
!
! October
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, real ALO, AHI, the low and high values.
!
! Input, integer N, the number of values.
!
! Output, real A(N), N evenly spaced values.
! Normally, A() = ALO and A(N) = AHI.
! However, if N = , then A() = 0.5*(ALO+AHI).
!
implicit none
!
integer n
!
real a(n)
real ahi
real alo
integer i
!
if ( n == ) then a() = 0.5E+00 * ( alo + ahi ) else do i = , n
a(i) = ( real ( n - i ) * alo + real ( i - ) * ahi ) / real ( n - )
end do end if return
end
subroutine rvec_order_type ( n, a, order )
!
!*******************************************************************************
!
!! RVEC_ORDER_TYPE determines if a real array is (non)strictly ascending/descending.
!
!
! Modified:
!
! July
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer N, the number of entries of the array.
!
! Input, real A(N), the array to be checked.
!
! Output, integer ORDER, order indicator:
! -, no discernable order;
! , all entries are equal;
! , ascending order;
! , strictly ascending order;
! , descending order;
! , strictly descending order.
!
implicit none
!
integer n
!
real a(n)
integer i
integer order
!
! Search for the first value not equal to A().
!
i = do i = i + if ( i > n ) then
order =
return
end if if ( a(i) > a() ) then if ( i == ) then
order =
else
order =
end if exit else if ( a(i) < a() ) then if ( i == ) then
order =
else
order =
end if exit end if end do
!
! Now we have a "direction". Examine subsequent entries.
!
do i = i +
if ( i > n ) then
exit
end if if ( order == ) then if ( a(i) < a(i-) ) then
order = -
exit
end if else if ( order == ) then if ( a(i) < a(i-) ) then
order = -
exit
else if ( a(i) == a(i-) ) then
order =
end if else if ( order == ) then if ( a(i) > a(i-) ) then
order = -
exit
end if else if ( order == ) then if ( a(i) > a(i-) ) then
order = -
exit
else if ( a(i) == a(i-) ) then
order =
end if end if end do return
end
subroutine rvec_print ( n, a, title )
!
!*******************************************************************************
!
!! RVEC_PRINT prints a real vector.
!
!
! Modified:
!
! December
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer N, the number of components of the vector.
!
! Input, real A(N), the vector to be printed.
!
! Input, character ( len = * ) TITLE, a title to be printed first.
! TITLE may be blank.
!
implicit none
!
integer n
!
real a(n)
integer i
character ( len = * ) title
!
if ( title /= ' ' ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) trim ( title )
end if write ( *, '(a)' ) ' '
do i = , n
write ( *, '(i6,g14.6)' ) i, a(i)
end do return
end
subroutine rvec_random ( alo, ahi, n, a )
!
!*******************************************************************************
!
!! RVEC_RANDOM returns a random real vector in a given range.
!
!
! Modified:
!
! February
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, real ALO, AHI, the range allowed for the entries.
!
! Input, integer N, the number of entries in the vector.
!
! Output, real A(N), the vector of randomly chosen values.
!
implicit none
!
integer n
!
real a(n)
real ahi
real alo
integer i
!
do i = , n
call r_random ( alo, ahi, a(i) )
end do return
end
subroutine rvec_sort_bubble_a ( n, a )
!
!*******************************************************************************
!
!! RVEC_SORT_BUBBLE_A ascending sorts a real array using bubble sort.
!
!
! Discussion:
!
! Bubble sort is simple to program, but inefficient. It should not
! be used for large arrays.
!
! Modified:
!
! February
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer N, the number of entries in the array.
!
! Input/output, real A(N).
! On input, an unsorted array.
! On output, the array has been sorted.
!
implicit none
!
integer n
!
real a(n)
integer i
integer j
!
do i = , n-
do j = i+, n
if ( a(i) > a(j) ) then
call r_swap ( a(i), a(j) )
end if
end do
end do return
end
subroutine s3_fs ( a1, a2, a3, n, b, x )
!
!*******************************************************************************
!
!! S3_FS factors and solves a tridiagonal linear system.
!
!
! Note:
!
! This algorithm requires that each diagonal entry be nonzero.
!
! Modified:
!
! December
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input/output, real A1(:N), A2(:N), A3(:N-).
! On input, the nonzero diagonals of the linear system.
! On output, the data in these vectors has been overwritten
! by factorization information.
!
! Input, integer N, the order of the linear system.
!
! Input/output, real 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 X(N), the solution of the linear system.
!
implicit none
!
integer n
!
real a1(:n)
real a2(:n)
real a3(:n-)
real b(n)
integer i
real x(n)
real xmult
!
! The diagonal entries can't be zero.
!
do i = , n
if ( a2(i) == 0.0E+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'S3_FS - Fatal error!'
write ( *, '(a,i6,a)' ) ' A2(', i, ') = 0.'
return
end if
end do do i = , n- xmult = a1(i) / a2(i-)
a2(i) = a2(i) - xmult * a3(i-) b(i) = b(i) - xmult * b(i-) end do xmult = a1(n) / a2(n-)
a2(n) = a2(n) - xmult * a3(n-) x(n) = ( b(n) - xmult * b(n-) ) / a2(n)
do i = n-, , -
x(i) = ( b(i) - a3(i) * x(i+) ) / a2(i)
end do return
end
subroutine sgtsl ( n, c, d, e, b, info )
!
!*******************************************************************************
!
!! SGTSL solves a general tridiagonal linear system.
!
!
! Reference:
!
! Dongarra, Moler, Bunch and Stewart,
! LINPACK User's Guide,
! SIAM, (Society for Industrial and Applied Mathematics),
! University City Science Center,
! Philadelphia, PA, -.
! ISBN ---X
!
! Modified:
!
! October
!
! Parameters:
!
! Input, integer N, the order of the tridiagonal matrix.
!
! Input/output, real C(N), contains the subdiagonal of the tridiagonal
! matrix in entries C(:N). On output, C is destroyed.
!
! Input/output, real D(N). On input, the diagonal of the matrix.
! On output, D is destroyed.
!
! Input/output, real E(N), contains the superdiagonal of the tridiagonal
! matrix in entries E(:N-). On output E is destroyed.
!
! Input/output, real B(N). On input, the right hand side. On output,
! the solution.
!
! Output, integer INFO, error flag.
! , normal value.
! K, the K-th element of the diagonal becomes exactly zero. The
! subroutine returns if this error condition is detected.
!
implicit none
!
integer n
!
real b(n)
real c(n)
real d(n)
real e(n)
integer info
integer k
real t
!
info =
c() = d() if ( n >= ) then d() = e()
e() = 0.0E+00
e(n) = 0.0E+00 do k = , n -
!
! Find the larger of the two rows, and interchange if necessary.
!
if ( abs ( c(k+) ) >= abs ( c(k) ) ) then
call r_swap ( c(k), c(k+) )
call r_swap ( d(k), d(k+) )
call r_swap ( e(k), e(k+) )
call r_swap ( b(k), b(k+) )
end if
!
! Fail if no nonzero pivot could be found.
!
if ( c(k) == 0.0E+00 ) then
info = k
return
end if
!
! Zero elements.
!
t = -c(k+) / c(k)
c(k+) = d(k+) + t * d(k)
d(k+) = e(k+) + t * e(k)
e(k+) = 0.0E+00
b(k+) = b(k+) + t * b(k) end do end if if ( c(n) == 0.0E+00 ) then
info = n
return
end if
!
! Back solve.
!
b(n) = b(n) / c(n) if ( n > ) then b(n-) = ( b(n-) - d(n-) * b(n) ) / c(n-) do k = n-, , -
b(k) = ( b(k) - d(k) * b(k+) - e(k) * b(k+) ) / c(k)
end do end if return
end
subroutine spline_b_val ( ndata, tdata, ydata, tval, yval )
!
!*******************************************************************************
!
!! 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.
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer NDATA, the number of data values.
!
! Input, real TDATA(NDATA), the abscissas of the data.
!
! Input, real YDATA(NDATA), the data values.
!
! Input, real TVAL, a point at which the spline is to be evaluated.
!
! Output, real YVAL, the value of the function at TVAL.
!
implicit none
!
integer ndata
!
real bval
integer left
integer right
real tdata(ndata)
real tval
real u
real ydata(ndata)
real yval
!
! Find the nearest interval [ TDATA(LEFT), TDATA(RIGHT) ] to TVAL.
!
call rvec_bracket ( ndata, tdata, tval, left, right )
!
! Evaluate the nonzero B spline basis functions in the interval,
! weighted by their corresponding data values.
!
u = ( tval - tdata(left) ) / ( tdata(right) - tdata(left) )
yval = 0.0E+00
!
! B function associated with node LEFT - , (or "phantom node"),
! evaluated in its 4th interval.
!
bval = ( 1.0E+00 - 3.0E+00 * u + 3.0E+00 * u** - u** ) / 6.0E+00
if ( left- > ) then
yval = yval + ydata(left-) * bval
else
yval = yval + ( 2.0E+00 * ydata() - ydata() ) * bval
end if
!
! B function associated with node LEFT,
! evaluated in its third interval.
!
bval = ( 4.0E+00 - 6.0E+00 * u** + 3.0E+00 * u** ) / 6.0E+00
yval = yval + ydata(left) * bval
!
! B function associated with node RIGHT,
! evaluated in its second interval.
!
bval = ( 1.0E+00 + 3.0E+00 * u + 3.0E+00 * u** - 3.0E+00 * u** ) / 6.0E+00
yval = yval + ydata(right) * bval
!
! B function associated with node RIGHT+, (or "phantom node"),
! evaluated in its first interval.
!
bval = u** / 6.0E+00
if ( right+ <= ndata ) then
yval = yval + ydata(right+) * bval
else
yval = yval + ( 2.0E+00 * ydata(ndata) - ydata(ndata-) ) * bval
end if return
end
subroutine spline_beta_val ( beta1, beta2, ndata, tdata, ydata, tval, yval )
!
!*******************************************************************************
!
!! 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 = and BETA2 = , the cubic beta spline will be the
! same as the cubic B spline approximant.
!
! With BETA1 = 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.
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, real BETA1, the skew or bias parameter.
! BETA1 = for no skew or bias.
!
! Input, real BETA2, the tension parameter.
! BETA2 = for no tension.
!
! Input, integer NDATA, the number of data values.
!
! Input, real TDATA(NDATA), the abscissas of the data.
!
! Input, real YDATA(NDATA), the data values.
!
! Input, real TVAL, a point at which the spline is to be evaluated.
!
! Output, real YVAL, the value of the function at TVAL.
!
implicit none
!
integer ndata
!
real a
real b
real beta1
real beta2
real bval
real c
real d
real delta
integer left
integer right
real tdata(ndata)
real tval
real u
real ydata(ndata)
real yval
!
! Find the nearest interval [ TDATA(LEFT), TDATA(RIGHT) ] to TVAL.
!
call rvec_bracket ( ndata, tdata, tval, left, right )
!
! Evaluate the nonzero beta spline basis functions in the interval,
! weighted by their corresponding data values.
!
u = ( tval - tdata(left) ) / ( tdata(right) - tdata(left) ) delta = 2.0E+00 + beta2 + 4.0E+00 * beta1 + 4.0E+00 * beta1** &
+ 2.0E+00 * beta1** yval = 0.0E+00
!
! Beta function associated with node LEFT - , (or "phantom node"),
! evaluated in its 4th interval.
!
bval = ( 2.0E+00 * beta1** * ( 1.0E+00 - u )** ) / delta if ( left- > ) then
yval = yval + ydata(left-) * bval
else
yval = yval + ( 2.0E+00 * ydata() - ydata() ) * bval
end if
!
! Beta function associated with node LEFT,
! evaluated in its third interval.
!
a = beta2 + 4.0E+00 * beta1 + 4.0E+00 * beta1** b = - 6.0E+00 * beta1 * ( 1.0E+00 - beta1 ) * ( 1.0E+00 + beta1 ) c = - 3.0E+00 * ( beta2 + 2.0E+00 * beta1** + 2.0E+00 * beta1** ) d = 2.0E+00 * ( beta2 + beta1 + beta1** + beta1** ) 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.0E+00 b = 6.0E+00 * beta1 c = 3.0E+00 * beta2 + 6.0E+00 * beta1** d = - 2.0E+00 * ( 1.0E+00 + beta2 + beta1 + beta1** ) bval = ( a + u * ( b + u * ( c + u * d ) ) ) / delta yval = yval + ydata(right) * bval
!
! Beta function associated with node RIGHT+, (or "phantom node"),
! evaluated in its first interval.
!
bval = 2.0E+00 * u** / delta if ( right+ <= ndata ) then
yval = yval + ydata(right+) * bval
else
yval = yval + ( 2.0E+00 * ydata(ndata) - ydata(ndata-) ) * bval
end if return
end
subroutine spline_constant_val ( ndata, tdata, ydata, tval, yval )
!
!*******************************************************************************
!
!! SPLINE_CONSTANT_VAL evaluates a piecewise constant spline at a point.
!
!
! Discussion:
!
! NDATA- 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).
!
! Modified:
!
! November
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer NDATA, the number of data points defining the spline.
!
! Input, real TDATA(NDATA-), the breakpoints. The values of TDATA should
! be distinct and increasing.
!
! Input, YDATA(NDATA), the values of the spline in the intervals
! defined by the breakpoints.
!
! Input, real TVAL, the point at which the spline is to be evaluated.
!
! Output, real YVAL, the value of the spline at TVAL.
!
implicit none
!
integer ndata
!
integer i
real tdata(ndata-)
real tval
real ydata(ndata)
real yval
!
do i = , ndata-
if ( tval <= tdata(i) ) then
yval = ydata(i)
return
end if
end do yval = ydata(ndata) return
end
subroutine spline_cubic_set ( n, t, y, ibcbeg, ybcbeg, ibcend, ybcend, ypp )
!
!*******************************************************************************
!
!! SPLINE_CUBIC_SET computes the second derivatives of a 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+), the form of
! the spline is
!
! SPL(T) = A(IVAL)
! + B(IVAL) * ( T - T(IVAL) )
! + C(IVAL) * ( T - T(IVAL) )**
! + D(IVAL) * ( T - T(IVAL) )**
!
! 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+) - Y(IVAL) ) / ( T(IVAL+) - T(IVAL) )
! - ( YPP(IVAL+) + * YPP(IVAL) ) * ( T(IVAL+) - T(IVAL) ) /
! C(IVAL) = YPP(IVAL) /
! D(IVAL) = ( YPP(IVAL+) - YPP(IVAL) ) / ( * ( T(IVAL+) - T(IVAL) ) )
!
! Since the first derivative of the spline is
!
! SPL'(T) = B(IVAL)
! + * C(IVAL) * ( T - T(IVAL) )
! + * D(IVAL) * ( T - T(IVAL) )**,
!
! the requirement that the first derivative be continuous at interior
! knot I results in a total of N- equations, of the form:
!
! B(IVAL-) + C(IVAL-) * (T(IVAL)-T(IVAL-))
! + * D(IVAL-) * (T(IVAL) - T(IVAL-))** = B(IVAL)
!
! or, setting H(IVAL) = T(IVAL+) - T(IVAL)
!
! ( Y(IVAL) - Y(IVAL-) ) / H(IVAL-)
! - ( YPP(IVAL) + * YPP(IVAL-) ) * H(IVAL-) /
! + YPP(IVAL-) * H(IVAL-)
! + ( YPP(IVAL) - YPP(IVAL-) ) * H(IVAL-) /
! =
! ( Y(IVAL+) - Y(IVAL) ) / H(IVAL)
! - ( YPP(IVAL+) + * YPP(IVAL) ) * H(IVAL) /
!
! or
!
! YPP(IVAL-) * H(IVAL-) + * YPP(IVAL) * ( H(IVAL-) + H(IVAL) )
! + YPP(IVAL) * H(IVAL)
! =
! * ( Y(IVAL+) - Y(IVAL) ) / H(IVAL)
! - * ( Y(IVAL) - Y(IVAL-) ) / H(IVAL-)
!
! Boundary conditions must be applied at the first and last knots.
! The resulting tridiagonal system can be solved for the YPP values.
!
! Modified:
!
! November
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer N, the number of data points; N must be at least .
!
! Input, real T(N), the points where data is specified.
! The values should be distinct, and increasing.
!
! Input, real Y(N), the data values to be interpolated.
!
! Input, integer IBCBEG, the left boundary condition flag:
!
! : the spline should be a quadratic over the first interval;
! : the first derivative at the left endpoint should be YBCBEG;
! : the second derivative at the left endpoint should be YBCBEG.
!
! Input, real YBCBEG, the left boundary value, if needed.
!
! Input, integer IBCEND, the right boundary condition flag:
!
! : the spline should be a quadratic over the last interval;
! : the first derivative at the right endpoint should be YBCEND;
! : the second derivative at the right endpoint should be YBCEND.
!
! Input, real YBCEND, the right boundary value, if needed.
!
! Output, real YPP(N), the second derivatives of the cubic spline.
!
implicit none
!
integer n
!
real diag(n)
integer i
integer ibcbeg
integer ibcend
real sub(:n)
real sup(:n-)
real t(n)
real y(n)
real ybcbeg
real ybcend
real ypp(n)
!
! Check.
!
if ( n <= ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'SPLINE_CUBIC_SET - Fatal error!'
write ( *, '(a)' ) ' The number of knots must be at least 2.'
write ( *, '(a,i6)' ) ' The input value of N = ', n
stop
end if do i = , n-
if ( t(i) >= t(i+) ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'SPLINE_CUBIC_SET - Fatal error!'
write ( *, '(a)' ) ' The knots must be strictly increasing, but'
write ( *, '(a,i6,a,g14.6)' ) ' T(', i,') = ', t(i)
write ( *, '(a,i6,a,g14.6)' ) ' T(',i+,') = ', t(i+)
stop
end if
end do
!
! Set the first equation.
!
if ( ibcbeg == ) then
ypp() = 0.0E+00
diag() = 1.0E+00
sup() = -1.0E+00
else if ( ibcbeg == ) then
ypp() = ( y() - y() ) / ( t() - t() ) - ybcbeg
diag() = ( t() - t() ) / 3.0E+00
sup() = ( t() - t() ) / 6.0E+00
else if ( ibcbeg == ) then
ypp() = ybcbeg
diag() = 1.0E+00
sup() = 0.0E+00
else
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'SPLINE_CUBIC_SET - Fatal error!'
write ( *, '(a)' ) ' The boundary flag IBCBEG must be 0, 1 or 2.'
write ( *, '(a,i6)' ) ' The input value is IBCBEG = ', ibcbeg
stop
end if
!
! Set the intermediate equations.
!
do i = , n-
ypp(i) = ( y(i+) - y(i) ) / ( t(i+) - t(i) ) &
- ( y(i) - y(i-) ) / ( t(i) - t(i-) )
sub(i) = ( t(i) - t(i-) ) / 6.0E+00
diag(i) = ( t(i+) - t(i-) ) / 3.0E+00
sup(i) = ( t(i+) - t(i) ) / 6.0E+00
end do
!
! Set the last equation.
!
if ( ibcend == ) then
ypp(n) = 0.0E+00
sub(n) = -1.0E+00
diag(n) = 1.0E+00
else if ( ibcend == ) then
ypp(n) = ybcend - ( y(n) - y(n-) ) / ( t(n) - t(n-) )
sub(n) = ( t(n) - t(n-) ) / 6.0E+00
diag(n) = ( t(n) - t(n-) ) / 3.0E+00
else if ( ibcend == ) then
ypp(n) = ybcend
sub(n) = 0.0E+00
diag(n) = 1.0E+00
else
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'SPLINE_CUBIC_SET - Fatal error!'
write ( *, '(a)' ) ' The boundary flag IBCEND must be 0, 1 or 2.'
write ( *, '(a,i6)' ) ' The input value is IBCEND = ', ibcend
stop
end if
!
! Special case:
! N = , IBCBEG = IBCEND = .
!
if ( n == .and. ibcbeg == .and. ibcend == ) then ypp() = 0.0E+00
ypp() = 0.0E+00
!
! Solve the linear system.
!
else call s3_fs ( sub, diag, sup, n, ypp, ypp ) end if return
end
subroutine spline_cubic_val ( n, t, y, ypp, tval, yval, ypval, yppval )
!
!*******************************************************************************
!
!! SPLINE_CUBIC_VAL evaluates a cubic spline at a specific 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+), the form of
! the spline is
!
! SPL(T) = A
! + B * ( T - T(IVAL) )
! + C * ( T - T(IVAL) )**
! + D * ( T - T(IVAL) )**
!
! Here:
! A = Y(IVAL)
! B = ( Y(IVAL+) - Y(IVAL) ) / ( T(IVAL+) - T(IVAL) )
! - ( YPP(IVAL+) + * YPP(IVAL) ) * ( T(IVAL+) - T(IVAL) ) /
! C = YPP(IVAL) /
! D = ( YPP(IVAL+) - YPP(IVAL) ) / ( * ( T(IVAL+) - T(IVAL) ) )
!
! Modified:
!
! November
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer N, the number of data values.
!
! Input, real T(N), the knot values.
!
! Input, real Y(N), the data values at the knots.
!
! Input, real YPP(N), the second derivatives of the spline at the knots.
!
! Input, real TVAL, a point, typically between T() and T(N), at
! which the spline is to be evalulated. If TVAL lies outside
! this range, extrapolation is used.
!
! Output, real YVAL, YPVAL, YPPVAL, the value of the spline, and
! its first two derivatives at TVAL.
!
implicit none
!
integer n
!
real dt
real h
integer left
integer right
real t(n)
real tval
real y(n)
real ypp(n)
real yppval
real ypval
real yval
!
! Determine the interval [T(LEFT), T(RIGHT)] that contains TVAL.
! Values below T() or above T(N) use extrapolation.
!
call rvec_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.0E+00 + ypp(left) / 3.0E+00 ) * h &
+ dt * ( 0.5E+00 * ypp(left) &
+ dt * ( ( ypp(right) - ypp(left) ) / ( 6.0E+00 * h ) ) ) ) ypval = ( y(right) - y(left) ) / h &
- ( ypp(right) / 6.0E+00 + ypp(left) / 3.0E+00 ) * h &
+ dt * ( ypp(left) &
+ dt * ( 0.5E+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 )
!
!*******************************************************************************
!
!! SPLINE_CUBIC_VAL2 evaluates a cubic spline at a specific 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+. The form of the spline is
!
! SPL(T) =
! A
! + B * ( T - T(LEFT) )
! + C * ( T - T(LEFT) )**
! + D * ( T - T(LEFT) )**
!
! Here:
! A = Y(LEFT)
! B = ( Y(RIGHT) - Y(LEFT) ) / ( T(RIGHT) - T(LEFT) )
! - ( YPP(RIGHT) + * YPP(LEFT) ) * ( T(RIGHT) - T(LEFT) ) /
! C = YPP(LEFT) /
! D = ( YPP(RIGHT) - YPP(LEFT) ) / ( * ( T(RIGHT) - T(LEFT) ) )
!
! Modified:
!
! November
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer N, the number of knots.
!
! Input, real T(N), the knot values.
!
! Input, real Y(N), the data values at the knots.
!
! Input, real YPP(N), the second derivatives of the spline at
! the knots.
!
! Input/output, integer LEFT, the suggested T interval to search.
! LEFT should be between and N-. 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 TVAL, a point, typically between T() and T(N), at
! which the spline is to be evalulated. If TVAL lies outside
! this range, extrapolation is used.
!
! Output, real YVAL, YPVAL, YPPVAL, the value of the spline, and
! its first two derivatives at TVAL.
!
implicit none
!
integer n
!
real dt
real h
integer left
integer right
real t(n)
real tval
real y(n)
real ypp(n)
real yppval
real ypval
real yval
!
! Determine the interval [T(LEFT), T(RIGHT)] that contains TVAL.
!
! What you want from RVEC_BRACKET3 is that TVAL is to be computed
! by the data in interval {T(LEFT), T(RIGHT)].
!
left =
call rvec_bracket3 ( n, t, tval, left )
right = left +
!
! In the interval LEFT, the polynomial is in terms of a normalized
! coordinate ( DT / H ) between and .
!
dt = tval - t(left)
h = t(right) - t(left) yval = y(left) + dt * ( ( y(right) - y(left) ) / h &
- ( ypp(right) / 6.0E+00 + ypp(left) / 3.0E+00 ) * h &
+ dt * ( 0.5E+00 * ypp(left) &
+ dt * ( ( ypp(right) - ypp(left) ) / ( 6.0E+00 * h ) ) ) ) ypval = ( y(right) - y(left) ) / h &
- ( ypp(right) / 6.0E+00 + ypp(left) / 3.0E+00 ) * h &
+ dt * ( ypp(left) &
+ dt * ( 0.5E+00 * ( ypp(right) - ypp(left) ) / h ) ) yppval = ypp(left) + dt * ( ypp(right) - ypp(left) ) / h return
end
subroutine spline_hermite_set ( ndata, tdata, c )
!
!*************************************************************************
!
!! SPLINE_HERMITE_SET sets up a piecewise cubic Hermite interpolant.
!
!
! Reference:
!
! Conte and de Boor,
! Algorithm CALCCF,
! Elementary Numerical Analysis,
! , page .
!
! Modified:
!
! April
!
! Parameters:
!
! Input, integer NDATA, the number of data points.
! NDATA must be at least .
!
! Input, real TDATA(NDATA), the abscissas of the data points.
! The entries of TDATA are assumed to be strictly increasing.
!
! Input/output, real C(,NDATA).
!
! On input, C(,I) and C(,I) should contain the value of the
! function and its derivative at TDATA(I), for I = to NDATA.
! These values will not be changed by this routine.
!
! On output, C(,I) and C(,I) contain the quadratic
! and cubic coefficients of the Hermite polynomial
! in the interval (TDATA(I), TDATA(I+)), for I= to NDATA-.
! C(,NDATA) and C(,NDATA) are set to .
!
! In the interval (TDATA(I), TDATA(I+)), the interpolating Hermite
! polynomial is given by
!
! SVAL(TVAL) = C(,I)
! + ( TVAL - TDATA(I) ) * ( C(,I)
! + ( TVAL - TDATA(I) ) * ( C(,I)
! + ( TVAL - TDATA(I) ) * C(,I) ) )
!
implicit none
!
integer ndata
!
real c(,ndata)
real divdif1
real divdif3
real dt
integer i
real tdata(ndata)
!
do i = , ndata-
dt = tdata(i+) - tdata(i)
divdif1 = ( c(,i+) - c(,i) ) / dt
divdif3 = c(,i) + c(,i+) - 2.0E+00 * divdif1
c(,i) = ( divdif1 - c(,i) - divdif3 ) / dt
c(,i) = divdif3 / dt**
end do c(,ndata) = 0.0E+00
c(,ndata) = 0.0E+00 return
end
subroutine spline_hermite_val ( ndata, tdata, c, tval, sval )
!
!*************************************************************************
!
!! 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.
!
! Reference:
!
! Conte and de Boor,
! Algorithm PCUBIC,
! Elementary Numerical Analysis,
! , page .
!
! Modified:
!
! April
!
! Parameters:
!
! Input, integer NDATA, the number of data points.
! NDATA is assumed to be at least .
!
! Input, real TDATA(NDATA), the abscissas of the data points.
! The entries of TDATA are assumed to be strictly increasing.
!
! Input, real C(,NDATA), contains the data computed by
! SPLINE_HERMITE_SET.
!
! Input, real TVAL, the point where the interpolant is to
! be evaluated.
!
! Output, real SVAL, the value of the interpolant at TVAL.
!
implicit none
!
integer ndata
!
real c(,ndata)
real dt
integer left
integer right
real sval
real tdata(ndata)
real tval
!
! Find the interval [ TDATA(LEFT), TDATA(RIGHT) ] that contains
! or is nearest to TVAL.
!
call rvec_bracket ( ndata, tdata, tval, left, right )
!
! Evaluate the cubic polynomial.
!
dt = tval - tdata(left) sval = c(,left) + dt * ( c(,left) + dt * ( c(,left) + dt * c(,left) ) ) return
end
subroutine spline_linear_int ( ndata, tdata, ydata, a, b, int_val )
!
!*******************************************************************************
!
!! SPLINE_LINEAR_INT evaluates the integral of a linear spline.
!
!
! Modified:
!
! November
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer NDATA, the number of data points defining the spline.
!
! Input, real 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 A, B, the interval over which the integral is desired.
!
! Output, real INT_VAL, the value of the integral.
!
implicit none
!
integer ndata
!
real a
real a_copy
integer a_left
integer a_right
real b
real b_copy
integer b_left
integer b_right
integer i_left
real int_val
real tdata(ndata)
real tval
real ydata(ndata)
real yp
real yval
!
int_val = 0.0E+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 rvec_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 rvec_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.0E+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+),...
! TDATA(B_LEFT-) 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.0E+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 - tval = ( tdata(i_left+) + tdata(i_left) ) / 2.0E+00 yp = ( ydata(i_left+) - ydata(i_left) ) / &
( tdata(i_left+) - tdata(i_left) ) yval = ydata(i_left) + ( tval - tdata(i_left) ) * yp int_val = int_val + yval * ( tdata(i_left + ) - tdata(i_left) ) end do tval = ( tdata(b_left) + b_copy ) / 2.0E+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 ( int_n, int_x, int_v, data_n, data_x, data_y )
!
!*******************************************************************************
!
!! SPLINE_LINEAR_INTSET sets a linear spline with given integral properties.
!
!
! Discussion:
!
! The user has in mind an interval, divided by INT_N+ 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.
!
! Modified:
!
! November
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer INT_N, the number of intervals.
!
! Input, real INT_X(INT_N+), the points that define the intervals.
! Interval I lies between INT_X(I) and INT_X(I+).
!
! Input, real INT_V(INT_N), the desired value of the integral of the
! linear spline over each interval.
!
! Output, integer DATA_N, the number of data points defining the spline.
! (This is the same as INT_N).
!
! Output, real DATA_X(DATA_N), DATA_Y(DATA_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 int_n
!
real c(int_n)
real d(int_n)
integer data_n
real data_x(int_n)
real data_y(int_n)
real e(int_n)
integer info
real int_v(int_n)
real int_x(int_n+)
!
! Set up the easy stuff.
!
data_n = int_n
data_x(:data_n) = 0.5E+00 * ( int_x(:data_n) + int_x(:data_n+) )
!
! Set up C, D, E, the coefficients of the linear system.
!
c() = 0.0E+00
c(:data_n-) = 1.0 &
- ( 0.5 * ( data_x(:data_n-) + int_x(:data_n-) ) &
- data_x(:data_n-) ) &
/ ( data_x(:data_n-) - data_x(:data_n-) )
c(data_n) = 0.0E+00 d() = int_x() - int_x() d(:data_n-) = 1.0 &
+ ( 0.5 * ( data_x(:data_n-) + int_x(:data_n-) ) &
- data_x(:data_n-) ) &
/ ( data_x(:data_n-) - data_x(:data_n-) ) &
- ( 0.5 * ( data_x(:data_n-) + int_x(:data_n) ) - data_x(:data_n-) ) &
/ ( data_x(:data_n) - data_x(:data_n-) ) d(data_n) = int_x(data_n+) - int_x(data_n) e() = 0.0E+00 e(:data_n-) = ( 0.5 * ( data_x(:data_n-) + int_x(:data_n) ) &
- data_x(:data_n-) ) / ( data_x(:data_n) - data_x(:data_n-) ) e(data_n) = 0.0E+00
!
! Set up DATA_Y, which begins as the right hand side of the linear system.
!
data_y() = int_v()
data_y(:data_n-) = 2.0E+00 * int_v(:data_n-) &
/ ( int_x(:int_n) - int_x(:int_n-) )
data_y(data_n) = int_v(data_n)
!
! Solve the linear system.
!
call sgtsl ( data_n, c, d, e, data_y, info ) if ( info /= ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'SPLINE_LINEAR_INTSET - Fatal error!'
write ( *, '(a)' ) ' The linear system is singular.'
stop
end if return
end
subroutine spline_linear_val ( ndata, tdata, ydata, tval, yval, ypval )
!
!*******************************************************************************
!
!! SPLINE_LINEAR_VAL evaluates a linear spline at a specific 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.
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer NDATA, the number of data points defining the spline.
!
! Input, real 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 TVAL, the point at which the spline is to be evaluated.
!
! Output, real 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 ndata
!
integer left
integer right
real tdata(ndata)
real tval
real ydata(ndata)
real ypval
real yval
!
! Find the interval [ TDATA(LEFT), TDATA(RIGHT) ] that contains, or is
! nearest to, TVAL.
!
call rvec_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 )
!
!*******************************************************************************
!
!! SPLINE_OVERHAUSER_NONUNI_VAL evaluates the nonuniform Overhauser spline.
!
!
! Discussion:
!
! The nonuniformity refers to the fact that the abscissas values
! need not be uniformly spaced.
!
! Diagnostics:
!
! The values of ALPHA and BETA have to be properly assigned.
! The basis matrices for the first and last interval have to
! be computed.
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer NDATA, the number of data points.
!
! Input, real TDATA(NDATA), the abscissas of the data points.
! The values of TDATA are assumed to be distinct and increasing.
!
! Input, real YDATA(NDATA), the data values.
!
! Input, real TVAL, the value where the spline is to
! be evaluated.
!
! Output, real YVAL, the value of the spline at TVAL.
!
implicit none
!
integer ndata
!
real alpha
real beta
integer left
real mbasis(,)
real mbasis_l(,)
real mbasis_r(,)
integer right
real tdata(ndata)
real tval
real ydata(ndata)
real yval
!
! Find the nearest interval [ TDATA(LEFT), TDATA(RIGHT) ] to TVAL.
!
call rvec_bracket ( ndata, tdata, tval, left, right )
!
! Evaluate the spline in the given interval.
!
if ( left == ) then alpha = 1.0E+00
call basis_matrix_overhauser_nul ( alpha, mbasis_l ) call basis_matrix_tmp ( , , mbasis_l, ndata, tdata, ydata, tval, yval ) else if ( left < ndata- ) then alpha = 1.0E+00
beta = 1.0E+00
call basis_matrix_overhauser_nonuni ( alpha, beta, mbasis ) call basis_matrix_tmp ( left, , mbasis, ndata, tdata, ydata, tval, yval ) else if ( left == ndata- ) then beta = 1.0E+00
call basis_matrix_overhauser_nur ( beta, mbasis_r ) call basis_matrix_tmp ( left, , mbasis_r, ndata, tdata, ydata, tval, yval ) end if return
end
subroutine spline_overhauser_uni_val ( ndata, tdata, ydata, tval, yval )
!
!*******************************************************************************
!
!! SPLINE_OVERHAUSER_UNI_VAL evaluates the uniform Overhauser spline.
!
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer NDATA, the number of data points.
!
! Input, real 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() = , TDATA() = , TDATA() = ...
!
! Input, real YDATA(NDATA), the data values.
!
! Input, real TVAL, the value where the spline is to
! be evaluated.
!
! Output, real YVAL, the value of the spline at TVAL.
!
implicit none
!
integer ndata
!
integer left
real mbasis(,)
real mbasis_l(,)
real mbasis_r(,)
integer right
real tdata(ndata)
real tval
real ydata(ndata)
real yval
!
! Find the nearest interval [ TDATA(LEFT), TDATA(RIGHT) ] to TVAL.
!
call rvec_bracket ( ndata, tdata, tval, left, right )
!
! Evaluate the spline in the given interval.
!
if ( left == ) then call basis_matrix_overhauser_uni_l ( mbasis_l ) call basis_matrix_tmp ( , , mbasis_l, ndata, tdata, ydata, tval, yval ) else if ( left < ndata- ) then call basis_matrix_overhauser_uni ( mbasis ) call basis_matrix_tmp ( left, , mbasis, ndata, tdata, ydata, tval, yval ) else if ( left == ndata- ) then call basis_matrix_overhauser_uni_r ( mbasis_r ) call basis_matrix_tmp ( left, , mbasis_r, ndata, tdata, ydata, tval, yval ) end if return
end
subroutine spline_overhauser_val ( ndim, ndata, tdata, ydata, tval, yval )
!
!*******************************************************************************
!
!! 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.
!
! Reference:
!
! H Brewer and D Anderson,
! Visual Interaction with Overhauser Curves and Surfaces,
! SIGGRAPH , pages -.
!
! E Catmull and R Rom,
! A Class of Local Interpolating Splines,
! in Computer Aided Geometric Design,
! edited by R Barnhill and R Reisenfeld,
! Academic Press, , pages -.
!
! David Rogers and Alan Adams,
! Mathematical Elements of Computer Graphics,
! McGraw Hill, , Second Edition, pages -.
!
! Modified:
!
! April
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer NDIM, the dimension of a single data point.
! NDIM must be at least . There is an internal limit on NDIM,
! called MAXDIM, which is presently set to .
!
! Input, integer NDATA, the number of data points.
! NDATA must be at least .
!
! Input, real TDATA(NDATA), the abscissas of the data points. The
! values in TDATA must be in strictly ascending order.
!
! Input, real YDATA(NDIM,NDATA), the data points corresponding to
! the abscissas.
!
! Input, real TVAL, the abscissa value at which the spline
! is to be evaluated. Normally, TDATA() <= TVAL <= T(NDATA), and
! the data will be interpolated. For TVAL outside this range,
! extrapolation will be used.
!
! Output, real YVAL(NDIM), the value of the spline at TVAL.
!
implicit none
!
integer, parameter :: MAXDIM =
integer ndata
integer ndim
!
integer i
integer left
integer order
integer right
real tdata(ndata)
real tval
real ydata(ndim,ndata)
real yl(MAXDIM)
real yr(MAXDIM)
real yval(ndim)
!
! Check.
!
call rvec_order_type ( ndata, tdata, order ) if ( order /= ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'SPLINE_OVERHAUSER_VAL - Fatal error!'
write ( *, '(a)' ) ' The data abscissas are not strictly ascending.'
stop
end if if ( ndata < ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'SPLINE_OVERHAUSER_VAL - Fatal error!'
write ( *, '(a)' ) ' NDATA < 3.'
stop
end if if ( ndim < ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'SPLINE_OVERHAUSER_VAL - Fatal error!'
write ( *, '(a)' ) ' NDIM < 1.'
stop
end if if ( ndim > maxdim ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'SPLINE_OVERHAUSER_VAL - Fatal error!'
write ( *, '(a)' ) ' NDIM > MAXDIM.'
stop
end if
!
! Locate the abscissa interval T(LEFT), T(LEFT+) nearest to or
! containing TVAL.
!
call rvec_bracket ( ndata, tdata, tval, left, right )
!
! Evaluate the "left hand" quadratic defined at T(LEFT-), T(LEFT), T(RIGHT).
!
if ( left- > ) then
call parabola_val2 ( ndim, ndata, tdata, ydata, left-, tval, yl )
end if
!
! Evaluate the "right hand" quadratic defined at T(LEFT), T(RIGHT), T(RIGHT+).
!
if ( right+ <= ndata ) then
call parabola_val2 ( ndim, ndata, tdata, ydata, left, tval, yr )
end if
!
! Average the quadratics.
!
if ( left == ) then yval(:ndim) = yr(:ndim) else if ( right < ndata ) then yval(:ndim) = ( ( tdata(right) - tval ) * yl(:ndim) &
+ ( tval - tdata(left) ) * yr(:ndim) ) / ( tdata(right) - tdata(left) ) else yval(:ndim) = yl(:ndim) end if return
end
subroutine spline_quadratic_val ( ndata, tdata, ydata, tval, yval, ypval )
!
!*******************************************************************************
!
!! SPLINE_QUADRATIC_VAL evaluates a quadratic spline at a specific 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.
!
! Modified:
!
! October
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer NDATA, the number of data points defining the spline.
! NDATA should be odd.
!
! Input, real 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 TVAL, the point at which the spline is to be evaluated.
!
! Output, real 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 ndata
!
real dif1
real dif2
integer left
integer right
real t1
real t2
real t3
real tdata(ndata)
real tval
real y1
real y2
real y3
real ydata(ndata)
real ypval
real yval
!
if ( mod ( ndata, ) == ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'SPLINE_QUADRATIC_VAL - Fatal error!'
write ( *, '(a)' ) ' NDATA must be odd.'
stop
end if
!
! Find the interval [ TDATA(LEFT), TDATA(RIGHT) ] that contains, or is
! nearest to, TVAL.
!
call rvec_bracket ( ndata, tdata, tval, left, right )
!
! Force LEFT to be odd.
!
if ( mod ( left, ) == ) then
left = left -
end if
!
! Copy out the three abscissas.
!
t1 = tdata(left)
t2 = tdata(left+)
t3 = tdata(left+) if ( t1 >= t2 .or. t2 >= t3 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'SPLINE_QUADRATIC_VAL - Fatal error!'
write ( *, '(a)' ) ' T1 >= T2 or T2 >= T3.'
stop
end if
!
! Construct and evaluate a parabolic interpolant for the data
! in each dimension.
!
y1 = ydata(left)
y2 = ydata(left+)
y3 = ydata(left+) 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.0E+00 * tval - t1 - t2 ) return
end
subroutine timestamp ( )
!
!*******************************************************************************
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!
! Example:
!
! May ::54.872 AM
!
! Modified:
!
! May
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! None
!
implicit none
!
character ( len = ) ampm
integer d
character ( len = ) date
integer h
integer m
integer mm
character ( len = ), parameter, dimension() :: month = (/ &
'January ', 'February ', 'March ', 'April ', &
'May ', 'June ', 'July ', 'August ', &
'September', 'October ', 'November ', 'December ' /)
integer n
integer s
character ( len = ) time
integer values()
integer y
character ( len = ) zone
!
call date_and_time ( date, time, zone, values ) y = values()
m = values()
d = values()
h = values()
n = values()
s = values()
mm = values() if ( h < ) then
ampm = 'AM'
else if ( h == ) then
if ( n == .and. s == ) then
ampm = 'Noon'
else
ampm = 'PM'
end if
else
h = h -
if ( h < ) then
ampm = 'PM'
else if ( h == ) then
if ( n == .and. s == ) then
ampm = 'Midnight'
else
ampm = 'AM'
end if
end if
end if write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return
end
样条曲线的Fortran程序的更多相关文章
-
求解轨道力学二体意义下的Lambert方程（兰伯特方程）的Fortran程序
轨道力学中二体问题下求解兰伯特方程. 老外写的Matlab程序,我把它转成了Fortran程序. !************************************************** ...
-
获取硬盘序列号的Fortran程序
以前写了个获取硬盘序列号的fortran程序,但未经实证 program FortranDemo Use Kernel32 Implicit None Interface SUBROUTINE Get ...
-
c++调用fortran程序中遇到的问题
一.C++动态调用Fortran DLL (1)创建FORTRAN DLL工程,生成forsubs.dll文件供调用. ! forsubs.f90 ! ! FUNCTIONS/SUBROUTINES ...
-
C++调用Fortran程序----动态链接方式
参考http://yxbwuhee.blog.sohu.com/143577510.html 一.C++动态调用Fortran DLL (1)创建FORTRAN DLL工程,生成forsubs.dll ...
-
Fortran程序调试中的“吐核”错误
在CentOS7上安装了PGI编译器,但是调试过程中遇到的“段错误(吐核)”一直让人很头疼. 通常采用在程序中增加屏幕输出代码的方式来追踪和定位出错的变量,比如下面这个样例程序就在第16行和第18行增 ...
-
visual studio 2005 编fortran程序，运行后dos窗口显示问题
比如程序: program main implicit none write(*,*) "AAAAAAAAAAAAAAAAAAAAAAAA" stop end 虽然可以看见DOS窗 ...
-
利用ctypes调用Fortran程序
本来python下面调用fortran最傻瓜方便的办法就是f2py,但是若fortran和C混合编程的代码,分别指定gfortran和gcc为编译器,在windows下面f2py直接报错 那么ctyp ...
-
Linux 下编译 有多个子程序文件的Fortran程序
第一种方法 ifort -o outprogram Source1.f90 Source2.f90 第二种 在主程序中include 'Source2.f90' program main call p ...
-
[转载：]C#与Fortran混合编程之本地调用Fortran动态链接库
前言 C#发展到现在,已是一门相当完善的语言,他基于C语言风格,演化于C++.并依靠强大的.NET底层框架.C#可以用来快速构建桌面及Web应用.然而在我们的实际工作中,尽管C#已经非常完善,但还是不 ...
随机推荐
-
JavaBean学习总结（上）
一.何为JavaBean: 遵循特定规则的Java类,必须为共有类: 1. 需要对成员属性私有化: 2. 需要无参的构造函数: 3. 需要通过public方法将私有属性暴露给其他程序,且方法遵循一定命 ...
-
高性能JavaScript笔记二（算法和流程控制、快速响应用户界面、Ajax）
循环 在javaScript中的四种循环中(for.for-in.while.do-while),只有for-in循环比其它几种明显要慢,另外三种速度区别不大 有一点需要注意的是,javascript ...
-
Windows Phone零距离开发(Composite Thread组合线程)
简洁流畅,快速响应是Windows Phone的特点也是他的买点,我们在开发App时候,也要在手机有限的硬件性能上最大限度做到UI快速响应,微软在优化手机快速响应这块做了很多底层优化工作,其中有一个就 ...
-
Echarts自适应浏览器大小
var myChart = echarts.init(document.getElementById('sitesChar')); var option = { title : { text: 'No ...
-
python---os模块使用详解
os模块调用操作系统接口的模块 相关方法或属性: getcwd() --- 获取当前的操作目录,等同于linux中的pwd命令. 调用:os.getcwd() chdir() --- 改变python ...
-
bzoj 2752: [HAOI2012]高速公路(road)
Description Y901高速公路是一条重要的交通纽带,*部门建设初期的投入以及使用期间的养护费用都不低,因此*在这条高速公路上设立了许多收费站.Y901高速公路是一条由N-1段路以及N个收 ...
-
c# http请求ajax页面
我们在用Http请求的时候,某些页面是ajax加载的,所以请求过来的页面数据不完整.也就是说ajax局部加载数据的地方,我们请求不到,这时候该怎么办呢? WebDriver+phantomjs 这两个 ...
-
docker安装及配置
docker下载安装(官方) 卸载旧版本 sudo yum remove docker docker-client docker-client-latest docker-common docker- ...
-
Day14 集合（一）
集合总体介绍 Java集合是java提供的工具包,包含了常用的数据结构:集合.链表.队列.栈.数组.映射等.Java集合工具包位置是java.util.*Java集合主要可以划分为4个部分:List列 ...
-
PC网页js调用本地应用程序
最近要现实一个在PC网页中实现点击按钮调用本地应用程序的功能 其实实现原理也非常简单, 首先注册一个本地注册表文件,指向本地应用程序路径 其次在网页中用js指向这个注册表文件,就可以实现网页调用本地应 ...