1 /* 2 3 -Procedure hrmint_c ( Hermite polynomial interpolation ) 4 5 -Abstract 6 7 Evaluate a Hermite interpolating polynomial at a specified 8 abscissa value. 9 10 -Disclaimer 11 12 THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE 13 CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. 14 GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE 15 ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE 16 PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" 17 TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY 18 WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A 19 PARTICULAR USE OR PURPOSE (AS set_c FORTH IN UNITED STATES UCC 20 SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE 21 SOFTWARE AND RELATED MATERIALS, HOWEVER USED. 22 23 IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA 24 BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT 25 LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, 26 INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, 27 REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE 28 REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. 29 30 RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF 31 THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY 32 CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE 33 ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. 34 35 -Required_Reading 36 37 None. 38 39 -Keywords 40 41 INTERPOLATION 42 MATH 43 POLYNOMIAL 44 45 */ 46 47 #include "SpiceUsr.h" 48 #include "SpiceZfc.h" 49 #undef hrmint_c 50 hrmint_c(SpiceInt n,ConstSpiceDouble * xvals,ConstSpiceDouble * yvals,SpiceDouble x,SpiceDouble * work,SpiceDouble * f,SpiceDouble * df)51 void hrmint_c ( SpiceInt n, 52 ConstSpiceDouble * xvals, 53 ConstSpiceDouble * yvals, 54 SpiceDouble x, 55 SpiceDouble * work, 56 SpiceDouble * f, 57 SpiceDouble * df ) 58 59 /* 60 61 -Brief_I/O 62 63 Variable I/O Description 64 -------- --- -------------------------------------------------- 65 n I Number of points defining the polynomial. 66 xvals I Abscissa values. 67 yvals I Ordinate and derivative values. 68 x I Point at which to interpolate the polynomial. 69 work I-O Work space array. 70 f O Interpolated function value at x. 71 df O Interpolated function's derivative at x. 72 73 -Detailed_Input 74 75 n is the number of points defining the polynomial. 76 The arrays xvals and yvals contain n and 2*n 77 elements respectively. 78 79 xvals is an array of length n containing abscissa values. 80 81 yvals is an array of length 2*n containing ordinate and 82 derivative values for each point in the domain 83 defined by xvals and n. The elements 84 85 yvals[ 2*i ] 86 yvals[ 2*i +1 ] 87 88 give the value and first derivative of the output 89 polynomial at the abscissa value 90 91 xvals[i] 92 93 where i ranges from 0 to n - 1. 94 95 work is a work space array. It is used by this routine 96 as a scratch area to hold intermediate results. 97 Generally sized at number of elements in yvals 98 time two. 99 100 x is the abscissa value at which the interpolating 101 polynomial and its derivative are to be evaluated. 102 103 -Detailed_Output 104 105 f, 106 df are the value and derivative at x of the unique 107 polynomial of degree 2n-1 that fits the points and 108 derivatives defined by xvals and yvals. 109 110 -Parameters 111 112 None. 113 114 -Exceptions 115 116 1) The error SPICE(DIVIDEBYZERO) signals from a routine 117 in the call tree if two input abscissas are equal, 118 119 2) The error SPICE(INVALIDSIZE) signals from a routine 120 in the call tree if n is less than 1. 121 122 3) This routine does not attempt to ward off or diagnose 123 arithmetic overflows. 124 125 -Files 126 127 None. 128 129 -Particulars 130 131 Users of this routine must choose the number of points to use 132 in their interpolation method. The authors of Reference [1] have 133 this to say on the topic: 134 135 Unless there is solid evidence that the interpolating function 136 is close in form to the true function f, it is a good idea to 137 be cautious about high-order interpolation. We 138 enthusiastically endorse interpolations with 3 or 4 points, we 139 are perhaps tolerant of 5 or 6; but we rarely go higher than 140 that unless there is quite rigorous monitoring of estimated 141 errors. 142 143 The same authors offer this warning on the use of the 144 interpolating function for extrapolation: 145 146 ...the dangers of extrapolation cannot be overemphasized: 147 An interpolating function, which is perforce an extrapolating 148 function, will typically go berserk when the argument x is 149 outside the range of tabulated values by more than the typical 150 spacing of tabulated points. 151 152 -Examples 153 154 Example: 155 156 Fit a 7th degree polynomial through the points ( x, y, y' ) 157 158 ( -1, 6, 3 ) 159 ( 0, 5, 0 ) 160 ( 3, 2210, 5115 ) 161 ( 5, 78180, 109395 ) 162 163 and evaluate this polynomial at x = 2. 164 165 #include <stdio.h> 166 #include "SpiceUsr.h" 167 168 int main() 169 { 170 171 /. 172 Local variables. 173 ./ 174 SpiceDouble answer; 175 SpiceDouble deriv; 176 SpiceDouble xvals [] = {-1., 0., 3., 5.}; 177 SpiceDouble yvals [] = {6., 3., 5., 0., 178 2210., 5115., 78180., 109395.}; 179 SpiceDouble work [2*8]; 180 SpiceDouble x = 2; 181 SpiceInt n = 4; 182 183 hrmint_c ( n, xvals, yvals, x, work, &answer, &deriv ); 184 185 printf( "answer = %lf\n", answer ); 186 printf( "deriv = %lf\n", deriv ); 187 188 return(0); 189 } 190 191 The returned value of 'answer' should be 141., and the returned 192 value of 'deriv' should be 456., since the unique 7th degree 193 polynomial that fits these constraints is 194 195 7 2 196 f(x) = x + 2x + 5 197 198 -Restrictions 199 200 None. 201 202 -Literature_References 203 204 [1] "Numerical Recipes---The Art of Scientific Computing" by 205 William H. Press, Brian P. Flannery, Saul A. Teukolsky, 206 William T. Vetterling (see sections 3.0 and 3.1). 207 208 [2] "Elementary Numerical Analysis---An Algorithmic Approach" 209 by S. D. Conte and Carl de Boor. See p. 64. 210 211 -Author_and_Institution 212 213 N.J. Bachman (JPL) 214 E.D. Wright (JPL) 215 216 -Version 217 218 -CSPICE Version 1.0.0, 24-AUG-2015 (EDW) 219 220 -Index_Entries 221 222 interpolate function using Hermite polynomial 223 Hermite interpolation 224 225 -& 226 */ 227 228 { /* Begin hrmint_c */ 229 230 /* 231 Local constants 232 */ 233 234 235 /* 236 Local macros 237 */ 238 239 240 /* 241 Local variables 242 */ 243 244 245 /* 246 Static variables 247 */ 248 249 250 /* 251 Participate in error tracing. 252 */ 253 254 chkin_c ( "hrmint_c" ); 255 256 /* 257 The f2c'd routine does the work. 258 */ 259 hrmint_( (integer * ) &n, 260 (doublereal *) xvals, 261 (doublereal *) yvals, 262 (doublereal *) &x, 263 (doublereal *) work, 264 (doublereal *) f, 265 (doublereal *) df); 266 267 chkout_c ( "hrmint_c" ); 268 269 } /* End hrmint_c */ 270