| Line |
Branch |
Exec |
Source |
| 1 |
|
|
module pchsp_95_m |
| 2 |
|
|
|
| 3 |
|
|
implicit none |
| 4 |
|
|
|
| 5 |
|
|
contains |
| 6 |
|
|
|
| 7 |
|
✗ |
function pchsp_95(x, f, ibeg, iend, vc_beg, vc_end) |
| 8 |
|
|
|
| 9 |
|
|
! PURPOSE: Set derivatives needed to determine the Hermite |
| 10 |
|
|
! representation of the cubic spline interpolant to given data, |
| 11 |
|
|
! with specified boundary conditions. |
| 12 |
|
|
|
| 13 |
|
|
! Part of the "pchip" package. |
| 14 |
|
|
|
| 15 |
|
|
! CATEGORY: E1A |
| 16 |
|
|
|
| 17 |
|
|
! KEYWORDS: cubic hermite interpolation, piecewise cubic |
| 18 |
|
|
! interpolation, spline interpolation |
| 19 |
|
|
|
| 20 |
|
|
! DESCRIPTION: "pchsp" stands for "Piecewise Cubic Hermite Spline" |
| 21 |
|
|
! Computes the Hermite representation of the cubic spline |
| 22 |
|
|
! interpolant to the data given in X and F satisfying the boundary |
| 23 |
|
|
! conditions specified by Ibeg, iend, vc_beg and VC_end. |
| 24 |
|
|
|
| 25 |
|
|
! The resulting piecewise cubic Hermite function may be evaluated |
| 26 |
|
|
! by "pchfe" or "pchfd". |
| 27 |
|
|
|
| 28 |
|
|
! NOTE: This is a modified version of C. de Boor's cubic spline |
| 29 |
|
|
! routine "cubspl". |
| 30 |
|
|
|
| 31 |
|
|
! REFERENCE: Carl de Boor, A Practical Guide to Splines, Springer, |
| 32 |
|
|
! 2001, pages 43-47 |
| 33 |
|
|
|
| 34 |
|
|
use assert_eq_m, only: assert_eq |
| 35 |
|
|
|
| 36 |
|
|
real, intent(in):: x(:) |
| 37 |
|
|
! independent variable values |
| 38 |
|
|
! The elements of X must be strictly increasing: |
| 39 |
|
|
! X(I-1) < X(I), I = 2...N. |
| 40 |
|
|
! (Error return if not.) |
| 41 |
|
|
! (error if size(x) < 2) |
| 42 |
|
|
|
| 43 |
|
|
real, intent(in):: f(:) |
| 44 |
|
|
! dependent variable values to be interpolated |
| 45 |
|
|
! F(I) is value corresponding to X(I). |
| 46 |
|
|
|
| 47 |
|
|
INTEGER, intent(in):: ibeg |
| 48 |
|
|
! desired boundary condition at beginning of data |
| 49 |
|
|
|
| 50 |
|
|
! IBEG = 0 to set pchsp_95(1) so that the third derivative is con- |
| 51 |
|
|
! tinuous at X(2). This is the "not a knot" condition |
| 52 |
|
|
! provided by de Boor's cubic spline routine CUBSPL. |
| 53 |
|
|
! This is the default boundary condition. |
| 54 |
|
|
! IBEG = 1 if first derivative at X(1) is given in VC_BEG. |
| 55 |
|
|
! IBEG = 2 if second derivative at X(1) is given in VC_BEG. |
| 56 |
|
|
! IBEG = 3 to use the 3-point difference formula for pchsp_95(1). |
| 57 |
|
|
! (Reverts to the default boundary condition if size(x) < 3 .) |
| 58 |
|
|
! IBEG = 4 to use the 4-point difference formula for pchsp_95(1). |
| 59 |
|
|
! (Reverts to the default boundary condition if size(x) < 4 .) |
| 60 |
|
|
|
| 61 |
|
|
! NOTES: |
| 62 |
|
|
! 1. An error return is taken if IBEG is out of range. |
| 63 |
|
|
! 2. For the "natural" boundary condition, use IBEG=2 and |
| 64 |
|
|
! VC_BEG=0. |
| 65 |
|
|
|
| 66 |
|
|
INTEGER, intent(in):: iend |
| 67 |
|
|
! IC(2) = IEND, desired condition at end of data. |
| 68 |
|
|
! IEND may take on the same values as IBEG, but applied to |
| 69 |
|
|
! derivative at X(N). In case IEND = 1 or 2, The value is given in VC_END. |
| 70 |
|
|
|
| 71 |
|
|
! NOTES: |
| 72 |
|
|
! 1. An error return is taken if IEND is out of range. |
| 73 |
|
|
! 2. For the "natural" boundary condition, use IEND=2 and |
| 74 |
|
|
! VC_END=0. |
| 75 |
|
|
|
| 76 |
|
|
REAL, intent(in), optional:: vc_beg |
| 77 |
|
|
! desired boundary value, as indicated above. |
| 78 |
|
|
! VC_BEG need be set only if IBEG = 1 or 2 . |
| 79 |
|
|
|
| 80 |
|
|
REAL, intent(in), optional:: vc_end |
| 81 |
|
|
! desired boundary value, as indicated above. |
| 82 |
|
|
! VC_END need be set only if Iend = 1 or 2 . |
| 83 |
|
|
|
| 84 |
|
|
real pchsp_95(size(x)) |
| 85 |
|
|
! derivative values at the data points |
| 86 |
|
|
! These values will determine the cubic spline interpolant |
| 87 |
|
|
! with the requested boundary conditions. |
| 88 |
|
|
! The value corresponding to X(I) is stored in |
| 89 |
|
|
! PCHSP_95(I), I=1...N. |
| 90 |
|
|
|
| 91 |
|
|
! LOCAL VARIABLES: |
| 92 |
|
✗ |
real wk(2, size(x)) ! real array of working storage |
| 93 |
|
|
INTEGER n ! number of data points |
| 94 |
|
|
integer ierr, ic(2) |
| 95 |
|
|
real vc(2) |
| 96 |
|
|
|
| 97 |
|
|
!------------------------------------------------------------------- |
| 98 |
|
|
|
| 99 |
|
✗ |
n = assert_eq(size(x), size(f), "pchsp_95 n") |
| 100 |
|
✗ |
if ((ibeg == 1 .or. ibeg == 2) .and. .not. present(vc_beg)) then |
| 101 |
|
✗ |
print *, "vc_beg required for IBEG = 1 or 2" |
| 102 |
|
✗ |
stop 1 |
| 103 |
|
|
end if |
| 104 |
|
✗ |
if ((iend == 1 .or. iend == 2) .and. .not. present(vc_end)) then |
| 105 |
|
✗ |
print *, "vc_end required for IEND = 1 or 2" |
| 106 |
|
✗ |
stop 1 |
| 107 |
|
|
end if |
| 108 |
|
✗ |
ic = (/ibeg, iend/) |
| 109 |
|
✗ |
if (present(vc_beg)) vc(1) = vc_beg |
| 110 |
|
✗ |
if (present(vc_end)) vc(2) = vc_end |
| 111 |
|
✗ |
call PCHSP(IC, VC, N, X, F, pchsp_95, 1, WK, size(WK), IERR) |
| 112 |
|
✗ |
if (ierr /= 0) stop 1 |
| 113 |
|
|
|
| 114 |
|
|
END function pchsp_95 |
| 115 |
|
|
|
| 116 |
|
|
end module pchsp_95_m |
| 117 |
|
|
|