1 /* chbder.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 CHBDER ( Derivatives of a Chebyshev expansion ) */
chbder_(doublereal * cp,integer * degp,doublereal * x2s,doublereal * x,integer * nderiv,doublereal * partdp,doublereal * dpdxs)9 /* Subroutine */ int chbder_(doublereal *cp, integer *degp, doublereal *x2s,
10 doublereal *x, integer *nderiv, doublereal *partdp, doublereal *dpdxs)
11 {
12 /* System generated locals */
13 integer i__1;
14
15 /* Local variables */
16 integer i__, j;
17 doublereal s, scale, s2;
18
19 /* $ Abstract */
20
21 /* Given the coefficients for the Chebyshev expansion of a */
22 /* polynomial, this returns the value of the polynomial and its */
23 /* first NDERIV derivatives evaluated at the input X. */
24
25 /* $ Disclaimer */
26
27 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
28 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
29 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
30 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
31 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
32 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
33 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
34 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
35 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
36 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
37
38 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
39 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
40 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
41 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
42 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
43 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
44
45 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
46 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
47 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
48 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
49
50 /* $ Required_Reading */
51
52 /* None. */
53
54 /* $ Keywords */
55
56 /* INTERPOLATION, MATH, POLYNOMIAL */
57
58 /* $ Declarations */
59 /* $ Brief_I/O */
60
61 /* VARIABLE I/O DESCRIPTION */
62 /* -------- --- -------------------------------------------------- */
63 /* CP I NDEG+1 Chebyshev polynomial coefficients. */
64 /* DEGP I Degree of polynomial. */
65 /* X2S I Transformation parameters of polynomial. */
66 /* X I Value for which the polynomial is to be evaluated */
67 /* NDERIV I The number of derivatives to compute */
68 /* PARTDP - Workspace provided for computing derivatives */
69 /* DPDXS(I) O Value of the I'th derivative of the polynomial */
70
71 /* $ Detailed_Input */
72
73 /* CP is an array of coefficients a polynomial with respect */
74 /* to the Chebyshev basis. The polynomial to be */
75 /* evaluated is assumed to be of the form: */
76
77 /* CP(DEGP+1)*T(DEGP,S) + CP(DEGP)*T(DEGP-1,S) + ... */
78
79 /* ... + CP(2)*T(1,S) + CP(1)*T(0,S) */
80
81 /* where T(I,S) is the I'th Chebyshev polynomial */
82 /* evaluated at a number S whose double precision */
83 /* value lies between -1 and 1. The value of S is */
84 /* computed from the input variables P(1), P(2) and X. */
85
86 /* DEGP is the degree of the Chebyshev polynomial to be */
87 /* evaluated. */
88
89 /* X2S is an array of two parameters. These parameters are */
90 /* used to transform the domain of the input variable X */
91 /* into the standard domain of the Chebyshev polynomial. */
92 /* X2S(1) should be a reference point in the domain of */
93 /* X; X2S(2) should be the radius by which points are */
94 /* allowed to deviate from the reference point and while */
95 /* remaining within the domain of X. The value of */
96 /* X is transformed into the value S given by */
97
98 /* S = ( X - X2S(1) ) / X2S(2) */
99
100 /* Typically X2S(1) is the midpoint of the interval over */
101 /* which X is allowed to vary and X2S(2) is the radius */
102 /* of the interval. */
103
104 /* The main reason for doing this is that a Chebyshev */
105 /* expansion is usually fit to data over a span */
106 /* from A to B where A and B are not -1 and 1 */
107 /* respectively. Thus to get the "best fit" the */
108 /* data was transformed to the interval [-1,1] and */
109 /* coefficients generated. These coefficients are */
110 /* not rescaled to the interval of the data so that */
111 /* the numerical "robustness" of the Chebyshev fit will */
112 /* not be lost. Consequently, when the "best fitting" */
113 /* polynomial needs to be evaluated at an intermediate */
114 /* point, the point of evaluation must be transformed */
115 /* in the same way that the generating points were */
116 /* transformed. */
117
118 /* X Value for which the polynomial is to be evaluated. */
119
120 /* NDERIV is the number of derivatives to be computed by the */
121 /* routine. NDERIV should be non-negative. */
122
123 /* PARTDP Is a work space used by the program to compute */
124 /* all of the desired derivatives. It should be declared */
125 /* in the calling program as */
126
127 /* DOUBLE PRECISION PARTDP(3, 0:NDERIV) */
128
129 /* $ Detailed_Output */
130
131 /* DPDXS(0) The value of the polynomial to be evaluated. It */
132 /* is given by */
133
134 /* CP(DEGP+1)*T(DEGP,S) + CP(DEGP)*T(DEGP-1,S) + ... */
135
136 /* ... + CP(2)*T(1,S) + CP(1)*T(0,S) */
137
138 /* where T(I,S) is the I'th Chebyshev polynomial */
139 /* evaluated at a number S = ( X - P(1) )/P(2) */
140
141 /* DPDXS(I) The value of the I'th derivative of the polynomial at */
142 /* X. (I ranges from 1 to NDERIV) It is given by */
143
144 /* [i] */
145 /* (1/P(2)**I) ( CP(DEGP+1)*T (DEGP,S) */
146
147 /* [i] */
148 /* + CP(DEGP)*T (DEGP-1,S) + ... */
149
150 /* . */
151 /* . */
152 /* . */
153 /* [i] */
154 /* ... + CP(2)*T (1,S) */
155
156 /* [i] */
157 /* + CP(1)*T (0,S) ) */
158
159 /* [i] */
160 /* where T(k,S) and T (I,S) are the k'th Chebyshev */
161 /* polynomial and its i'th derivative respectively, */
162 /* evaluated at the number S = ( X - X2S(1) )/X2S(2). */
163
164 /* $ Parameters */
165
166 /* None. */
167
168 /* $ Particulars */
169
170 /* This routine computes the value of a Chebyshev polynomial */
171 /* expansion and the derivatives of the expansion with respect to X. */
172 /* The polynomial is given by */
173
174 /* CP(DEGP+1)*T(DEGP,S) + CP(DEGP)*T(DEGP-1,S) + ... */
175
176 /* ... + CP(2)*T(1,S) + CP(1)*T(0,S) */
177
178 /* where */
179
180 /* S = ( X - X2S(1) ) / X2S(2) */
181
182 /* and */
183
184 /* T(i,S) is the i'th Chebyshev polynomial of the first kind */
185 /* evaluated at S. */
186
187
188 /* $ Examples */
189
190 /* Depending upon the user's needs, there are 3 routines available */
191 /* for evaluating Chebyshev polynomials. */
192
193 /* CHBVAL for evaluating a Chebyshev polynomial when no */
194 /* derivatives are desired. */
195
196 /* CHBINT for evaluating a Chebyshev polynomial and its */
197 /* first derivative. */
198
199 /* CHBDER for evaluating a Chebyshev polynomial and a user */
200 /* or application dependent number of derivatives. */
201
202 /* Of these 3 the one most commonly employed by NAIF software */
203 /* is CHBINT as it is used to interpolate ephemeris state */
204 /* vectors which requires the evaluation of a polynomial */
205 /* and its derivative. When no derivatives are desired one */
206 /* should use CHBVAL, or when more than one or an unknown */
207 /* number of derivatives are desired one should use CHBDER. */
208
209 /* The code fragment below illustrates how this routine might */
210 /* be used to obtain points for plotting a polynomial */
211 /* and its derivatives. */
212
213 /* fetch the pieces needed for describing the polynomial */
214 /* to be evaluated. */
215
216 /* READ (*,*) DEGP, ( CP(I), I = 1, DEG+1 ), NDERIV, BEG, END */
217
218 /* check to see that BEG is actually less than END */
219
220 /* IF ( BEG .GE. END ) THEN */
221
222 /* take some appropriate action */
223
224 /* ELSE */
225
226 /* X2S(1) = ( END + BEG ) / 2.0D0 */
227 /* X2S(2) = ( END - BEG ) / 2.0D0 */
228
229 /* END IF */
230
231 /* STEP = END - BEG / <number of points used for plotting> */
232 /* X = BEG */
233
234 /* DO WHILE ( X .LE. END ) */
235
236 /* CALL CHBDER ( CP, DEGP, X2S , X, NDERIV, PARTDP, DPDXS ) */
237
238 /* do something with the pairs ( X, DPDXS(0)),(X,DPDXS(1)), */
239 /* (X,DPDXS(2)) ... (X,DPDXS(NDERIV)) */
240
241 /* X = X + STEP */
242
243 /* END DO */
244
245 /* $ Restrictions */
246
247 /* The user must be sure that the provided workspace is declared */
248 /* properly in the calling routine. The proper declaration is: */
249
250 /* INTEGER NDERIV */
251 /* PARAMETER ( NDERIV = the desired number of derivatives ) */
252 /* DOUBLE PRECISION PARTDP (3, 0:NDERIV) */
253
254 /* If for some reason a parameter is not passed to this routine in */
255 /* NDERIV, the user should make sure that the value of NDERIV is not */
256 /* so large that the work space provided is inadequate. */
257
258 /* One needs to be careful that the value (X-X2S(1)) / X2S(2) lies */
259 /* between -1 and 1. Otherwise, the routine may fail spectacularly */
260 /* (for example with a floating point overflow). */
261
262 /* While this routine will compute derivatives of the input */
263 /* polynomial, the user should consider how accurately the */
264 /* derivatives of the Chebyshev fit, match the derivatives of the */
265 /* function it approximates. */
266
267
268
269 /* $ Exceptions */
270
271 /* Error free */
272
273 /* No tests are performed for exceptional values ( NDERIV negative, */
274 /* DEGP negative, etc.) This routine is expected to be used at a low */
275 /* level in ephemeris evaluations. For that reason it has been */
276 /* elected as a routine that will not participate in error handling. */
277
278 /* $ Files */
279
280 /* None. */
281
282 /* $ Author_and_Institution */
283
284 /* N.J. Bachman (JPL) */
285 /* W.L. Taber (JPL) */
286
287 /* $ Literature_References */
288
289 /* "Numerical Recipes -- The Art of Scientific Computing" by */
290 /* William H. Press, Brian P. Flannery, Saul A. Teukolsky, */
291 /* Willam T. Vetterling. (See Clenshaw's Recurrence Formula) */
292
293 /* "The Chebyshev Polynomials" by Theodore J. Rivlin */
294
295 /* "CRC Handbook of Tables for Mathematics" */
296
297 /* $ Version */
298
299 /* - SPICELIB Version 1.0.1, 10-MAR-1992 (WLT) */
300
301 /* Comment section for permuted index source lines was added */
302 /* following the header. */
303
304 /* - SPICELIB Version 1.0.0, 31-JAN-1990 (WLT) */
305
306 /* -& */
307 /* $ Index_Entries */
308
309 /* derivatives of a chebyshev expansion */
310
311 /* -& */
312 /* $ Revisions */
313
314 /* - Beta Version 1.0.1, 16-FEB-1988 (WLT) (NJB) */
315
316 /* The Error free specification was added to the routine as */
317 /* well as an explanation for this designation. Examples added. */
318 /* Declaration of unused variable RECIP removed. */
319 /* -& */
320
321 /* Local variables */
322
323
324 /* Transform X to S and initialize temporary variables. */
325
326 s = (*x - x2s[0]) / x2s[1];
327 s2 = s * 2.;
328 j = *degp + 1;
329 i__1 = *nderiv;
330 for (i__ = 0; i__ <= i__1; ++i__) {
331 partdp[i__ * 3] = 0.;
332 partdp[i__ * 3 + 1] = 0.;
333 }
334
335 /* Evaluate the polynomial ... */
336
337 while(j > 1) {
338 partdp[2] = partdp[1];
339 partdp[1] = partdp[0];
340 partdp[0] = cp[j - 1] + (s2 * partdp[1] - partdp[2]);
341
342 /* ... and its derivatives using recursion. */
343
344 scale = 2.;
345 i__1 = *nderiv;
346 for (i__ = 1; i__ <= i__1; ++i__) {
347 partdp[i__ * 3 + 2] = partdp[i__ * 3 + 1];
348 partdp[i__ * 3 + 1] = partdp[i__ * 3];
349 partdp[i__ * 3] = partdp[(i__ - 1) * 3 + 1] * scale + partdp[i__ *
350 3 + 1] * s2 - partdp[i__ * 3 + 2];
351 scale += 2.;
352 }
353 --j;
354 }
355 dpdxs[0] = cp[0] + (s * partdp[0] - partdp[1]);
356 scale = 1.;
357 i__1 = *nderiv;
358 for (i__ = 1; i__ <= i__1; ++i__) {
359 dpdxs[i__] = partdp[(i__ - 1) * 3] * scale + partdp[i__ * 3] * s -
360 partdp[i__ * 3 + 1];
361 scale += 1;
362 }
363
364 /* Scale the k'th derivative w.r.t S by (1/X2S(2)**k) so that we have */
365 /* the derivatives */
366
367 /* 2 3 4 5 */
368 /* d P(S) d P(S) d P(S) d P(S) d P(S) */
369 /* ------ ------ ------ ------ ------ */
370 /* 2 3 4 5 */
371 /* dX dX dX dX dX */
372
373
374 /* NOTE: In the loop that follows we perform division instead of */
375 /* multiplying by reciprocals so that the algorithm matches */
376 /* CHBINT. If multiplication by reciprocals is performed */
377 /* CHBINT and CHBDER (although mathematically equivalent) will */
378 /* not produce identical results for the first derivative. */
379
380
381 scale = x2s[1];
382 i__1 = *nderiv;
383 for (i__ = 1; i__ <= i__1; ++i__) {
384 dpdxs[i__] /= scale;
385 scale = x2s[1] * scale;
386 }
387 return 0;
388 } /* chbder_ */
389
390