1 /* hrmesp.f -- translated by f2c (version 19980913).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include "f2c.h"
7
8 /* $Procedure HRMESP ( Hermite polynomial interpolation, equal spacing ) */
hrmesp_(integer * n,doublereal * first,doublereal * step,doublereal * yvals,doublereal * x,doublereal * work,doublereal * f,doublereal * df)9 /* Subroutine */ int hrmesp_(integer *n, doublereal *first, doublereal *step,
10 doublereal *yvals, doublereal *x, doublereal *work, doublereal *f,
11 doublereal *df)
12 {
13 /* System generated locals */
14 integer yvals_dim1, work_dim1, work_offset, i__1, i__2, i__3, i__4, i__5,
15 i__6, i__7;
16
17 /* Builtin functions */
18 integer s_rnge(char *, integer, char *, integer);
19
20 /* Local variables */
21 doublereal temp;
22 integer this__, prev, next;
23 doublereal newx;
24 integer i__, j;
25 extern /* Subroutine */ int chkin_(char *, ftnlen);
26 doublereal denom, c1, c2, xi;
27 extern /* Subroutine */ int sigerr_(char *, ftnlen), chkout_(char *,
28 ftnlen), setmsg_(char *, ftnlen), errint_(char *, integer *,
29 ftnlen);
30 extern logical return_(void);
31 doublereal xij;
32
33 /* $ Abstract */
34
35 /* Evaluate, at a specified point, an Hermite interpolating */
36 /* polynomial for a specified set of coordinate pairs whose */
37 /* abscissas are equally spaced. */
38
39 /* $ Disclaimer */
40
41 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
42 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
43 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
44 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
45 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
46 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
47 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
48 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
49 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
50 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
51
52 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
53 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
54 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
55 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
56 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
57 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
58
59 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
60 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
61 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
62 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
63
64 /* $ Required_Reading */
65
66 /* None. */
67
68 /* $ Keywords */
69
70 /* INTERPOLATION */
71 /* POLYNOMIAL */
72
73 /* $ Declarations */
74 /* $ Brief_I/O */
75
76 /* Variable I/O Description */
77 /* -------- --- -------------------------------------------------- */
78 /* N I Number of points defining the polynomial. */
79 /* FIRST I First abscissa value. */
80 /* STEP I Step size. */
81 /* YVALS I Ordinate and derivative values. */
82 /* X I Point at which to interpolate the polynomial. */
83 /* WORK I-O Work space array. */
84 /* F O Interpolated function value at X. */
85 /* DF O Interpolated function's derivative at X. */
86
87 /* $ Detailed_Input */
88
89 /* N is the number of points defining the polynomial. */
90 /* The array YVALS contains 2*N elements. */
91
92 /* FIRST, */
93 /* STEP are, respectively, a starting abscissa value and a */
94 /* step size that define the set of abscissa values */
95
96 /* FIRST + (I-1) * STEP, I = 1, ..., N */
97
98 /* STEP must be non-zero. */
99
100
101 /* YVALS is an array of length 2*N containing ordinate and */
102 /* derivative values for each point in the domain */
103 /* defined by FIRST, STEP, and N. The elements */
104
105 /* YVALS( 2*I - 1 ) */
106 /* YVALS( 2*I ) */
107
108 /* give the value and first derivative of the output */
109 /* polynomial at the abscissa value */
110
111 /* FIRST + I * STEP */
112
113 /* where I ranges from 1 to N. */
114
115
116 /* WORK is a work space array. It is used by this routine */
117 /* as a scratch area to hold intermediate results. */
118
119
120 /* X is the abscissa value at which the interpolating */
121 /* polynomial and its derivative are to be evaluated. */
122
123 /* $ Detailed_Output */
124
125 /* F, */
126 /* DF are the value and derivative at X of the unique */
127 /* polynomial of degree 2N-1 that fits the points and */
128 /* derivatives defined by FIRST, STEP, and YVALS. */
129
130 /* $ Parameters */
131
132 /* None. */
133
134 /* $ Exceptions */
135
136 /* 1) If STEP is zero, the error SPICE(INVALIDSTEPSIZE) will be */
137 /* signaled. */
138
139 /* 2) If N is less than 1, the error SPICE(INVALIDSIZE) is */
140 /* signaled. */
141
142 /* 3) This routine does not attempt to ward off or diagnose */
143 /* arithmetic overflows. */
144
145 /* $ Files */
146
147 /* None. */
148
149 /* $ Particulars */
150
151 /* Users of this routine must choose the number of points to use */
152 /* in their interpolation method. The authors of Reference [1] have */
153 /* this to say on the topic: */
154
155 /* Unless there is solid evidence that the interpolating function */
156 /* is close in form to the true function f, it is a good idea to */
157 /* be cautious about high-order interpolation. We */
158 /* enthusiastically endorse interpolations with 3 or 4 points, we */
159 /* are perhaps tolerant of 5 or 6; but we rarely go higher than */
160 /* that unless there is quite rigorous monitoring of estimated */
161 /* errors. */
162
163 /* The same authors offer this warning on the use of the */
164 /* interpolating function for extrapolation: */
165
166 /* ...the dangers of extrapolation cannot be overemphasized: */
167 /* An interpolating function, which is perforce an extrapolating */
168 /* function, will typically go berserk when the argument x is */
169 /* outside the range of tabulated values by more than the typical */
170 /* spacing of tabulated points. */
171
172 /* $ Examples */
173
174
175 /* 1) Fit a 7th degree polynomial through the points ( x, y, y' ) */
176
177 /* ( -1, 6, 3 ) */
178 /* ( 1, 8, 11 ) */
179 /* ( 3, 2210, 5115 ) */
180 /* ( 5, 78180, 109395 ) */
181
182 /* and evaluate this polynomial at x = 2. */
183
184
185 /* PROGRAM TEST_HRMINT */
186
187 /* DOUBLE PRECISION ANSWER */
188 /* DOUBLE PRECISION DERIV */
189 /* DOUBLE PRECISION FIRST */
190 /* DOUBLE PRECISION STEP */
191 /* DOUBLE PRECISION YVALS (8) */
192 /* DOUBLE PRECISION WORK (8,2) */
193 /* INTEGER N */
194
195
196 /* N = 4 */
197
198 /* YVALS(1) = 6.D0 */
199 /* YVALS(2) = 3.D0 */
200 /* YVALS(3) = 8.D0 */
201 /* YVALS(4) = 11.D0 */
202 /* YVALS(5) = 2210.D0 */
203 /* YVALS(6) = 5115.D0 */
204 /* YVALS(7) = 78180.D0 */
205 /* YVALS(8) = 109395.D0 */
206
207 /* FIRST = -1.D0 */
208 /* STEP = 2.D0 */
209
210 /* CALL HRMESP ( N, FIRST, STEP, YVALS, */
211 /* . 2.D0, WORK, ANSWER, DERIV ) */
212
213 /* WRITE (*,*) 'ANSWER = ', ANSWER */
214 /* WRITE (*,*) 'DERIV = ', DERIV */
215 /* END */
216
217
218 /* The returned value of ANSWER should be 141.D0, and the returned */
219 /* derivative value should be 456.D0, since the unique 7th degree */
220 /* polynomial that fits these constraints is */
221
222 /* 7 2 */
223 /* f(x) = x + 2x + 5 */
224
225 /* $ Restrictions */
226
227 /* None. */
228
229 /* $ Literature_References */
230
231 /* [1] "Numerical Recipes---The Art of Scientific Computing" by */
232 /* William H. Press, Brian P. Flannery, Saul A. Teukolsky, */
233 /* William T. Vetterling (see sections 3.0 and 3.1). */
234
235 /* [2] "Elementary Numerical Analysis---An Algorithmic Approach" */
236 /* by S. D. Conte and Carl de Boor. See p. 64. */
237
238 /* $ Author_and_Institution */
239
240 /* N.J. Bachman (JPL) */
241
242 /* $ Version */
243
244 /* - SPICELIB Version 1.2.1, 28-JAN-2014 (NJB) */
245
246 /* Fixed a few comment typos. */
247
248 /* - SPICELIB Version 1.2.0, 31-JAN-2002 (EDW) */
249
250 /* Added the use of DBLE to convert integer values */
251 /* used in DOUBLE PRECISION calculations. */
252
253 /* - SPICELIB Version 1.1.0, 28-DEC-2001 (NJB) */
254
255 /* Blanks following final newline were truncated to */
256 /* suppress compilation warnings on the SGI-N32 platform. */
257
258 /* - SPICELIB Version 1.0.0, 01-MAR-2000 (NJB) */
259
260 /* -& */
261 /* $ Index_Entries */
262
263 /* interpolate function using Hermite polynomial */
264 /* Hermite interpolation */
265
266 /* -& */
267
268 /* SPICELIB functions */
269
270
271 /* Local variables */
272
273
274 /* Check in only if an error is detected. */
275
276 /* Parameter adjustments */
277 work_dim1 = *n << 1;
278 work_offset = work_dim1 + 1;
279 yvals_dim1 = *n << 1;
280
281 /* Function Body */
282 if (return_()) {
283 return 0;
284 }
285
286 /* No data, no interpolation. */
287
288 if (*n < 1) {
289 chkin_("HRMESP", (ftnlen)6);
290 setmsg_("Array size must be positive; was #.", (ftnlen)35);
291 errint_("#", n, (ftnlen)1);
292 sigerr_("SPICE(INVALIDSIZE)", (ftnlen)18);
293 chkout_("HRMESP", (ftnlen)6);
294 return 0;
295 }
296
297 /* The step size must be non-zero. */
298
299 if (*step == 0.) {
300 chkin_("HRMESP", (ftnlen)6);
301 setmsg_("Step size was zero.", (ftnlen)19);
302 sigerr_("SPICE(INVALIDSTEPSIZE)", (ftnlen)22);
303 chkout_("HRMESP", (ftnlen)6);
304 return 0;
305 }
306
307 /* We can simplify the interpolation problem by shifting */
308 /* and scaling the abscissa values so that they start at 1 */
309 /* and are separated by a unit step. All we need to do here is */
310 /* shift and scale X. */
311
312 newx = (*x - *first) / *step + 1.;
313
314 /* For consistency with our scaled horizontal axis, we'll have */
315 /* scale our local derivative values by STEP, and scale our final */
316 /* computed derivative by 1/STEP. */
317
318 /* Copy the input array into WORK. Scale the derivatives at this */
319 /* step. After this, the first column of WORK represents the first */
320 /* column of our triangular interpolation table. */
321
322 i__1 = (*n << 1) - 1;
323 for (i__ = 1; i__ <= i__1; i__ += 2) {
324 work[(i__2 = i__ + work_dim1 - work_offset) < work_dim1 << 1 && 0 <=
325 i__2 ? i__2 : s_rnge("work", i__2, "hrmesp_", (ftnlen)331)] =
326 yvals[(i__3 = i__ - 1) < yvals_dim1 && 0 <= i__3 ? i__3 :
327 s_rnge("yvals", i__3, "hrmesp_", (ftnlen)331)];
328 }
329 i__1 = *n << 1;
330 for (i__ = 2; i__ <= i__1; i__ += 2) {
331 work[(i__2 = i__ + work_dim1 - work_offset) < work_dim1 << 1 && 0 <=
332 i__2 ? i__2 : s_rnge("work", i__2, "hrmesp_", (ftnlen)335)] =
333 yvals[(i__3 = i__ - 1) < yvals_dim1 && 0 <= i__3 ? i__3 :
334 s_rnge("yvals", i__3, "hrmesp_", (ftnlen)335)] * *step;
335 }
336
337 /* Compute the second column of the interpolation table: this */
338 /* consists of the N-1 values obtained by evaluating the */
339 /* first-degree interpolants at NEWX. We'll also evaluate the */
340 /* derivatives of these interpolants at NEWX and save the results in */
341 /* the second column of WORK. Because the derivative computations */
342 /* depend on the function computations from the previous column in */
343 /* the interpolation table, and because the function interpolation */
344 /* overwrites the previous column of interpolated function values, */
345 /* we must evaluate the derivatives first. */
346
347 i__1 = *n - 1;
348 for (i__ = 1; i__ <= i__1; ++i__) {
349 c1 = (doublereal) (i__ + 1) - newx;
350 c2 = newx - (doublereal) i__;
351
352 /* The second column of WORK contains interpolated derivative */
353 /* values. */
354
355 /* The odd-indexed interpolated derivatives are simply the input */
356 /* derivatives, after scaling. */
357
358 prev = (i__ << 1) - 1;
359 this__ = prev + 1;
360 next = this__ + 1;
361 work[(i__2 = prev + (work_dim1 << 1) - work_offset) < work_dim1 << 1
362 && 0 <= i__2 ? i__2 : s_rnge("work", i__2, "hrmesp_", (ftnlen)
363 366)] = work[(i__3 = this__ + work_dim1 - work_offset) <
364 work_dim1 << 1 && 0 <= i__3 ? i__3 : s_rnge("work", i__3,
365 "hrmesp_", (ftnlen)366)];
366
367 /* The even-indexed interpolated derivatives are the slopes of */
368 /* the linear interpolating polynomials for adjacent input */
369 /* abscissa/ordinate pairs. No scaling is needed here. */
370
371 work[(i__2 = this__ + (work_dim1 << 1) - work_offset) < work_dim1 <<
372 1 && 0 <= i__2 ? i__2 : s_rnge("work", i__2, "hrmesp_", (
373 ftnlen)373)] = work[(i__3 = next + work_dim1 - work_offset) <
374 work_dim1 << 1 && 0 <= i__3 ? i__3 : s_rnge("work", i__3,
375 "hrmesp_", (ftnlen)373)] - work[(i__4 = prev + work_dim1 -
376 work_offset) < work_dim1 << 1 && 0 <= i__4 ? i__4 : s_rnge(
377 "work", i__4, "hrmesp_", (ftnlen)373)];
378
379 /* The first column of WORK contains interpolated function values. */
380 /* The odd-indexed entries are the linear Taylor polynomials, */
381 /* each input abscissa value, evaluated at NEWX. */
382
383 temp = work[(i__2 = this__ + work_dim1 - work_offset) < work_dim1 <<
384 1 && 0 <= i__2 ? i__2 : s_rnge("work", i__2, "hrmesp_", (
385 ftnlen)380)] * (newx - (doublereal) i__) + work[(i__3 = prev
386 + work_dim1 - work_offset) < work_dim1 << 1 && 0 <= i__3 ?
387 i__3 : s_rnge("work", i__3, "hrmesp_", (ftnlen)380)];
388 work[(i__2 = this__ + work_dim1 - work_offset) < work_dim1 << 1 && 0
389 <= i__2 ? i__2 : s_rnge("work", i__2, "hrmesp_", (ftnlen)383)]
390 = c1 * work[(i__3 = prev + work_dim1 - work_offset) <
391 work_dim1 << 1 && 0 <= i__3 ? i__3 : s_rnge("work", i__3,
392 "hrmesp_", (ftnlen)383)] + c2 * work[(i__4 = next + work_dim1
393 - work_offset) < work_dim1 << 1 && 0 <= i__4 ? i__4 : s_rnge(
394 "work", i__4, "hrmesp_", (ftnlen)383)];
395 work[(i__2 = prev + work_dim1 - work_offset) < work_dim1 << 1 && 0 <=
396 i__2 ? i__2 : s_rnge("work", i__2, "hrmesp_", (ftnlen)386)] =
397 temp;
398 }
399
400 /* The last column entries were not computed by the preceding loop; */
401 /* compute them now. */
402
403 work[(i__1 = (*n << 1) - 1 + (work_dim1 << 1) - work_offset) < work_dim1
404 << 1 && 0 <= i__1 ? i__1 : s_rnge("work", i__1, "hrmesp_", (
405 ftnlen)394)] = work[(i__2 = (*n << 1) + work_dim1 - work_offset) <
406 work_dim1 << 1 && 0 <= i__2 ? i__2 : s_rnge("work", i__2, "hrme"
407 "sp_", (ftnlen)394)];
408 work[(i__1 = (*n << 1) - 1 + work_dim1 - work_offset) < work_dim1 << 1 &&
409 0 <= i__1 ? i__1 : s_rnge("work", i__1, "hrmesp_", (ftnlen)395)] =
410 work[(i__2 = (*n << 1) + work_dim1 - work_offset) < work_dim1 <<
411 1 && 0 <= i__2 ? i__2 : s_rnge("work", i__2, "hrmesp_", (ftnlen)
412 395)] * (newx - *n) + work[(i__3 = (*n << 1) - 1 + work_dim1 -
413 work_offset) < work_dim1 << 1 && 0 <= i__3 ? i__3 : s_rnge("work",
414 i__3, "hrmesp_", (ftnlen)395)];
415
416 /* Compute columns 3 through 2*N of the table. */
417
418 i__1 = (*n << 1) - 1;
419 for (j = 2; j <= i__1; ++j) {
420 i__2 = (*n << 1) - j;
421 for (i__ = 1; i__ <= i__2; ++i__) {
422
423 /* In the theoretical construction of the interpolation table, */
424 /* there are 2*N abscissa values, since each input abcissa */
425 /* value occurs with multiplicity two. In this theoretical */
426 /* construction, the Jth column of the interpolation table */
427 /* contains results of evaluating interpolants that span J+1 */
428 /* consecutive abscissa values. The indices XI and XIJ below */
429 /* are used to pick the correct abscissa values out of this */
430 /* sequence of 2*N values. */
431
432 xi = (doublereal) ((i__ + 1) / 2);
433 xij = (doublereal) ((i__ + j + 1) / 2);
434 c1 = xij - newx;
435 c2 = newx - xi;
436 denom = xij - xi;
437
438 /* Compute the interpolated derivative at NEWX for the Ith */
439 /* interpolant. This is the derivative with respect to NEWX of */
440 /* the expression for the interpolated function value, which is */
441 /* the second expression below. This derivative computation */
442 /* is done first because it relies on the interpolated function */
443 /* values from the previous column of the interpolation table. */
444
445 /* The derivative expression here corresponds to equation */
446 /* 2.35 on page 64 in reference [2]. */
447
448 work[(i__3 = i__ + (work_dim1 << 1) - work_offset) < work_dim1 <<
449 1 && 0 <= i__3 ? i__3 : s_rnge("work", i__3, "hrmesp_", (
450 ftnlen)433)] = (c1 * work[(i__4 = i__ + (work_dim1 << 1)
451 - work_offset) < work_dim1 << 1 && 0 <= i__4 ? i__4 :
452 s_rnge("work", i__4, "hrmesp_", (ftnlen)433)] + c2 * work[
453 (i__5 = i__ + 1 + (work_dim1 << 1) - work_offset) <
454 work_dim1 << 1 && 0 <= i__5 ? i__5 : s_rnge("work", i__5,
455 "hrmesp_", (ftnlen)433)] + (work[(i__6 = i__ + 1 +
456 work_dim1 - work_offset) < work_dim1 << 1 && 0 <= i__6 ?
457 i__6 : s_rnge("work", i__6, "hrmesp_", (ftnlen)433)] -
458 work[(i__7 = i__ + work_dim1 - work_offset) < work_dim1 <<
459 1 && 0 <= i__7 ? i__7 : s_rnge("work", i__7, "hrmesp_", (
460 ftnlen)433)])) / denom;
461
462 /* Compute the interpolated function value at NEWX for the Ith */
463 /* interpolant. */
464
465 work[(i__3 = i__ + work_dim1 - work_offset) < work_dim1 << 1 && 0
466 <= i__3 ? i__3 : s_rnge("work", i__3, "hrmesp_", (ftnlen)
467 440)] = (c1 * work[(i__4 = i__ + work_dim1 - work_offset)
468 < work_dim1 << 1 && 0 <= i__4 ? i__4 : s_rnge("work",
469 i__4, "hrmesp_", (ftnlen)440)] + c2 * work[(i__5 = i__ +
470 1 + work_dim1 - work_offset) < work_dim1 << 1 && 0 <=
471 i__5 ? i__5 : s_rnge("work", i__5, "hrmesp_", (ftnlen)440)
472 ]) / denom;
473 }
474 }
475
476 /* Our interpolated function value is sitting in WORK(1,1) at this */
477 /* point. The interpolated derivative is located in WORK(1,2). */
478 /* We must undo the scaling of the derivative. We've already */
479 /* checked that STEP is non-zero. */
480
481 *f = work[(i__1 = work_dim1 + 1 - work_offset) < work_dim1 << 1 && 0 <=
482 i__1 ? i__1 : s_rnge("work", i__1, "hrmesp_", (ftnlen)452)];
483 *df = work[(i__1 = (work_dim1 << 1) + 1 - work_offset) < work_dim1 << 1 &&
484 0 <= i__1 ? i__1 : s_rnge("work", i__1, "hrmesp_", (ftnlen)453)]
485 / *step;
486 return 0;
487 } /* hrmesp_ */
488
489