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