CoastalME (Coastal Modelling Environment)
Simulates the long-term behaviour of complex coastlines
Loading...
Searching...
No Matches
hermite_cubic.cpp
Go to the documentation of this file.
1
12
13/* ===============================================================================================================================
14
15 This file is part of CoastalME, the Coastal Modelling Environment.
16
17 CoastalME is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.
18
19 This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
20
21 You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
22
23===============================================================================================================================*/
24#include <iostream>
25using std::cerr;
26using std::endl;
27using std::ios;
28
29#include "hermite_cubic.h"
30
31//===============================================================================================================================
32//
33// Purpose:
34//
35// R8VEC_BRACKET3 finds the interval containing or nearest a given value.
36//
37// Discussion:
38//
39// An R8VEC is a vector of R8's.
40//
41// The routine always returns the index LEFT of the sorted array
42// T with the property that either
43// * T is contained in the interval [ T[LEFT], T[LEFT+1] ], or
44// * T < T[LEFT] = T[0], or
45// * T > T[LEFT+1] = T[N-1].
46//
47// The routine is useful for interpolation problems, where
48// the abscissa must be located within an interval of data
49// abscissas for interpolation, or the "nearest" interval
50// to the (extreme) abscissa must be found so that extrapolation
51// can be carried out.
52//
53// This version of the function has been revised so that the value of
54// LEFT that is returned uses the 0-based indexing natural to C++.
55//
56// Licensing:
57//
58// This code is distributed under the GNU LGPL license.
59//
60// Modified:
61//
62// 30 April 2009
63//
64// Author:
65//
66// John Burkardt
67//
68// Parameters:
69//
70// Input, int N, length of the input array.
71//
72// Input, double T[N], an array that has been sorted into ascending order.
73//
74// Input, double TVAL, a value to be bracketed by entries of T.
75//
76// Input/output, int *LEFT.
77// On input, if 0 <= LEFT <= N-2, LEFT is taken as a suggestion for the
78// interval [ T[LEFT-1] T[LEFT] ] in which TVAL lies. This interval
79// is searched first, followed by the appropriate interval to the left
80// or right. After that, a binary search is used.
81// On output, LEFT is set so that the interval [ T[LEFT], T[LEFT+1] ]
82// is the closest to TVAL; it either contains TVAL, or else TVAL
83// lies outside the interval [ T[0], T[N-1] ].
84//
85//===============================================================================================================================
87//===============================================================================================================================
88void r8vec_bracket3(int const n, double const* t, double const tval, int* left)
89{
90 int high;
91 int low;
92 int mid;
93
94 // Check the input data
95 if (n < 2)
96 {
97 cerr << endl
98 << "R8VEC_BRACKET3 - Fatal error! N must be at least 2." << endl;
99 return;
100 }
101
102 // If *left is not between 0 and n-2, set it to the middle value
103 if ((*left < 0) || (n - 2 < *left))
104 {
105 *left = (n - 1) / 2;
106 }
107
108 // CASE 1: TVAL < T[*LEFT]: search for TVAL in (T[I],T[I+1]), for I = 0 to *LEFT-1.
109 if (tval < t[*left])
110 {
111 if (*left == 0)
112 {
113 return;
114 }
115
116 else if (*left == 1)
117 {
118 *left = 0;
119 return;
120 }
121
122 else if (t[*left - 1] <= tval)
123 {
124 *left = *left - 1;
125 return;
126 }
127
128 else if (tval <= t[1])
129 {
130 *left = 0;
131 return;
132 }
133
134 //
135 // ...Binary search for TVAL in (T[I],T[I+1]), for I = 1 to *LEFT-2.
136 //
137 low = 1;
138 high = *left - 2;
139
140 for (;;)
141 {
142 if (low == high)
143 {
144 *left = low;
145 return;
146 }
147
148 mid = (low + high + 1) / 2;
149
150 if (t[mid] <= tval)
151 {
152 low = mid;
153 }
154
155 else
156 {
157 high = mid - 1;
158 }
159 }
160 }
161
162 //
163 // CASE 2: T[*LEFT+1] < TVAL:
164 // Search for TVAL in (T[I],T[I+1]) for intervals I = *LEFT+1 to N-2.
165 //
166 else if (t[*left + 1] < tval)
167 {
168 if (*left == n - 2)
169 {
170 return;
171 }
172
173 else if (*left == n - 3)
174 {
175 *left = *left + 1;
176 return;
177 }
178
179 else if (tval <= t[*left + 2])
180 {
181 *left = *left + 1;
182 return;
183 }
184
185 else if (t[n - 2] <= tval)
186 {
187 *left = n - 2;
188 return;
189 }
190
191 // ...Binary search for TVAL in (T[I],T[I+1]) for intervals I = *LEFT+2 to N-3.
192 low = *left + 2;
193 high = n - 3;
194
195 for (;;)
196 {
197 if (low == high)
198 {
199 *left = low;
200 return;
201 }
202
203 mid = (low + high + 1) / 2;
204
205 if (t[mid] <= tval)
206 {
207 low = mid;
208 }
209
210 else
211 {
212 high = mid - 1;
213 }
214 }
215 }
216
217 // CASE 3: T[*LEFT] <= TVAL <= T[*LEFT+1]: T is just where the user said it might be.
218 else
219 {
220 }
221
222 return;
223}
224
225//===============================================================================================================================
226//
227// Purpose:
228//
229// HERMITE_CUBIC_VALUE evaluates a Hermite cubic polynomial.
230//
231// Discussion:
232//
233// The input arguments can be vectors.
234//
235// Licensing:
236//
237// This code is distributed under the GNU LGPL license.
238//
239// Modified:
240//
241// 13 February 2011
242//
243// Author:
244//
245// John Burkardt
246//
247// Reference:
248//
249// Fred Fritsch, Ralph Carlson,
250// Monotone Piecewise Cubic Interpolation,
251// SIAM Journal on Numerical Analysis,
252// Volume 17, Number 2, April 1980, pages 238-246.
253//
254// Parameters:
255//
256// Input, double X1, F1, D1, the left endpoint, function value
257// and derivative.
258//
259// Input, double X2, F2, D2, the right endpoint, function value
260// and derivative.
261//
262// Input, int N, the number of evaluation points.
263//
264// Input, double X[N], the points at which the Hermite cubic
265// is to be evaluated.
266//
267// Output, double F[N], D[N], S[N], T[N], the value and first
268// three derivatives of the Hermite cubic at X.
269//
270//===============================================================================================================================
272void hermite_cubic_value(double const x1, double const f1, double const d1, double const x2, double const f2, double const d2, int const n, double const* x, double* f, double* d, double* s, double* t)
273{
274 double c2;
275 double c3;
276 double df;
277 double h;
278 int i;
279
280 h = x2 - x1;
281 df = (f2 - f1) / h;
282
283 c2 = -(2.0 * d1 - 3.0 * df + d2) / h;
284 c3 = (d1 - 2.0 * df + d2) / h / h;
285
286 for (i = 0; i < n; i++)
287 {
288 f[i] = f1 + (x[i] - x1) * (d1 + (x[i] - x1) * (c2 + (x[i] - x1) * c3));
289 d[i] = d1 + (x[i] - x1) * (2.0 * c2 + (x[i] - x1) * 3.0 * c3);
290 s[i] = 2.0 * c2 + (x[i] - x1) * 6.0 * c3;
291 t[i] = 6.0 * c3;
292 }
293
294 return;
295}
296
297//===============================================================================================================================
298//
299// Purpose:
300//
301// HERMITE_CUBIC_SPLINE_VALUE evaluates a Hermite cubic spline.
302//
303// Licensing:
304//
305// This code is distributed under the GNU LGPL license.
306//
307// Modified:
308//
309// 13 February 2011
310//
311// Author:
312//
313// John Burkardt
314//
315// Reference:
316//
317// Fred Fritsch, Ralph Carlson,
318// Monotone Piecewise Cubic Interpolation,
319// SIAM Journal on Numerical Analysis,
320// Volume 17, Number 2, April 1980, pages 238-246.
321//
322// Parameters:
323//
324// Input, int NN, the number of data points.
325//
326// Input, double XN[NN], the coordinates of the data points.
327// The entries in XN must be in strictly ascending order.
328//
329// Input, double FN[NN], the function values.
330//
331// Input, double DN[NN], the derivative values.
332//
333// Input, int N, the number of sample points.
334//
335// Input, double X[N], the coordinates of the sample points.
336//
337// Output, double F[N], the function value at the sample points.
338//
339// Output, double D[N], the derivative value at the sample points.
340//
341// Output, double S[N], the second derivative value at the
342// sample points.
343//
344// Output, double T[N], the third derivative value at the
345// sample points.
346//
347//===============================================================================================================================
349//===============================================================================================================================
350void hermite_cubic_spline_value(int const nn, double* const xn, double* const fn, double* const dn, int const n, double* const x, double* f, double* d, double* s, double* t)
351{
352 int left = n / 2;
353
354 for (int i = 0; i < n; i++)
355 {
356 r8vec_bracket3(nn, xn, x[i], &left);
357
358 hermite_cubic_value(xn[left], fn[left], dn[left], xn[left + 1], fn[left + 1], dn[left + 1], 1, x + i, f + i, d + i, s + i, t + i);
359 }
360
361 return;
362}
void hermite_cubic_value(double const x1, double const f1, double const d1, double const x2, double const f2, double const d2, int const n, double const *x, double *f, double *d, double *s, double *t)
This is part of a C++ library from http://people.sc.fsu.edu/~jburkardt/cpp_src/hermite_cubic/hermite_...
void r8vec_bracket3(int const n, double const *t, double const tval, int *left)
This is part of a C++ library from http://people.sc.fsu.edu/~jburkardt/cpp_src/hermite_cubic/hermite_...
void hermite_cubic_spline_value(int const nn, double *const xn, double *const fn, double *const dn, int const n, double *const x, double *f, double *d, double *s, double *t)
This is part of a C++ library from http://people.sc.fsu.edu/~jburkardt/cpp_src/hermite_cubic/hermite_...
Definitions of some routines from the hermite_cubic library.