1 /* stmp03.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 STMP03 ( Stumpff functions 0 through 3 ) */
stmp03_(doublereal * x,doublereal * c0,doublereal * c1,doublereal * c2,doublereal * c3)9 /* Subroutine */ int stmp03_(doublereal *x, doublereal *c0, doublereal *c1,
10 doublereal *c2, doublereal *c3)
11 {
12 /* Initialized data */
13
14 static logical first = TRUE_;
15
16 /* System generated locals */
17 integer i__1;
18
19 /* Builtin functions */
20 integer s_rnge(char *, integer, char *, integer);
21 double log(doublereal), sqrt(doublereal), cosh(doublereal), sinh(
22 doublereal), cos(doublereal), sin(doublereal);
23
24 /* Local variables */
25 integer i__;
26 doublereal y, z__;
27 extern /* Subroutine */ int chkin_(char *, ftnlen);
28 extern doublereal dpmax_(void);
29 extern /* Subroutine */ int errdp_(char *, doublereal *, ftnlen);
30 static doublereal pairs[20], lbound;
31 extern /* Subroutine */ int sigerr_(char *, ftnlen), chkout_(char *,
32 ftnlen), setmsg_(char *, ftnlen);
33
34 /* $ Abstract */
35
36 /* Compute the values of the Stumpff functions C_0 through C_3 at */
37 /* a specified point. */
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 /* CONIC */
71 /* MATH */
72 /* UTILITY */
73
74 /* $ Declarations */
75 /* $ Brief_I/O */
76
77 /* Variable I/O Description */
78 /* -------- --- -------------------------------------------------- */
79 /* X I Argument to each Stumpff function C_0 to C_3. */
80 /* C0 O Value of C_0(X) */
81 /* C1 O Value of C_1(X) */
82 /* C2 O Value of C_2(X) */
83 /* C3 O Value of C_3(X) */
84 /* TRUNC P Number of terms needed in Maclaurin series for C_3. */
85
86 /* $ Detailed_Input */
87
88 /* X is the argument to use in each of the Stumpff functions */
89 /* C_0, C_1, C_2, and C_3. */
90
91 /* $ Detailed_Output */
92
93 /* C0 are the values of the Stumpff functions */
94 /* C1 C_0(X), C_1(X), C_2(X), and C_3(X). */
95 /* C2 */
96 /* C3 */
97
98 /* $ Parameters */
99
100 /* TRUNC The Maclaurin series for C_3 and C_2 respectively are: */
101
102 /* 2 3 k */
103 /* 1 X X X (-X) */
104 /* C_3(X) = --- - --- + --- - --- + . . . + ----------. . . */
105 /* 3! 5! 7! 9! (3 + 2*K)! */
106
107 /* and */
108
109 /* 2 3 k */
110 /* 1 X X X (-X) */
111 /* C_2(X) = --- - --- + --- - --- + . . . + ----------. . . */
112 /* 2! 4! 6! 8! (2 + 2*K)! */
113
114 /* These series are used in the evaluation of C_3 and C_2. */
115 /* Thus, it is necessary to make a decision about where to */
116 /* truncate the series in our evaluation of C_3 and C_2. */
117
118 /* TRUNC is used to tell this routine where to truncate */
119 /* the Maclaurin series for C_3 and C_2. */
120
121 /* The value of TRUNC for your machine is the smallest */
122 /* integer such that */
123
124 /* 1 */
125 /* 1.0D0 + ---------- = 1.0D0 */
126 /* (2*TRUNC)! */
127
128 /* The following program will (if compiled and linked) */
129 /* will produce the values of TRUNC for your machine. */
130
131 /* INTEGER TRUNC */
132
133 /* DOUBLE PRECISION DENOM */
134 /* DOUBLE PRECISION FACTR */
135
136 /* DOUBLE PRECISION X */
137
138 /* DENOM = 2.0D0 */
139 /* FACTR = 2.0D0 */
140 /* TRUNC = 1 */
141
142 /* X = 1.0D0 / DENOM */
143
144 /* DO WHILE ( 1.0D0 + X .GT. 1.0D0 ) */
145 /* DENOM = DENOM * (2.0D0+FACTR) * (1.0D0+FACTR) */
146 /* FACTR = FACTR + 2.0D0 */
147 /* TRUNC = TRUNC + 1 */
148 /* X = 1.0D0 / DENOM */
149 /* END DO */
150
151 /* WRITE (*,*) 'The value of TRUNC is: ', TRUNC */
152
153 /* END */
154
155 /* $ Exceptions */
156
157 /* 1) If the input value of X is not in the domain of values */
158 /* for which the Stumpff functions can be computed, the error */
159 /* SPICE(VALUEOUTOFRANGE) is signalled. */
160
161 /* The range of valid inputs is from -[ln(2) + ln(DPMAX)]**2 */
162 /* to DPMAX. */
163
164 /* $ Files */
165
166 /* None. */
167
168 /* $ Particulars */
169
170 /* This routine computes the values of the Stumpff functions C_0, */
171 /* C_1, C_2, and C_3 at the input X. */
172
173 /* The Stumpff function C_k(X) for k = 0, 1, ... is given by the */
174 /* series: */
175
176 /* 2 3 m */
177 /* 1 X X X (-X) */
178 /* C_k(X) = --- - ------ + ------ - ------ + . . . + ------- + . . . */
179 /* k! (k+2)! (k+4)! (k+6)! (k+2m)! */
180
181
182 /* These series converge for all real values of X. */
183
184
185 /* $ Examples */
186
187 /* For positive X, */
188
189 /* C_0(X) = COS ( DSQRT(X) ) */
190
191
192 /* SIN ( DSQRT(X) ) */
193 /* C_1(X) = --------------- */
194 /* DSQRT(X) */
195
196
197 /* 1 - COS ( DSQRT(X) ) */
198 /* C_2(X) = --------------------- */
199 /* X */
200
201
202
203 /* 1 - SIN ( DSQRT(X) ) / DSQRT(X) */
204 /* C_3(X) = ---------------------------------- */
205 /* X */
206
207 /* Thus the following block of code can be used to check this */
208 /* routine for reasonableness: */
209
210 /* INTEGER I */
211
212 /* DOUBLE PRECISION X */
213 /* DOUBLE PRECISION ROOTX */
214
215 /* DOUBLE PRECISION TC0 */
216 /* DOUBLE PRECISION TC1 */
217 /* DOUBLE PRECISION TC2 */
218 /* DOUBLE PRECISION TC3 */
219
220 /* DOUBLE PRECISION C0 */
221 /* DOUBLE PRECISION C1 */
222 /* DOUBLE PRECISION C2 */
223 /* DOUBLE PRECISION C3 */
224
225 /* DO I = 1, 10 */
226
227 /* X = DBLE (I) */
228 /* ROOTX = DSQRT(X) */
229
230 /* TC0 = COS ( ROOTX ) */
231 /* TC1 = SIN ( ROOTX ) / ROOTX */
232
233 /* TC2 = ( 1.0D0 - COS( ROOTX ) ) / X */
234 /* TC3 = ( 1.0D0 - SIN( ROOTX ) / ROOTX ) / X */
235
236 /* CALL STMP03 ( X, C0, C1, C2, C3 ) */
237
238 /* WRITE (*,*) */
239 /* WRITE (*,*) 'Expected - Computed for X = ', X */
240 /* WRITE (*,*) */
241 /* WRITE (*,*) 'Delta C0 :', TC0 - C0 */
242 /* WRITE (*,*) 'Delta C1 :', TC1 - C1 */
243 /* WRITE (*,*) 'Delta C2 :', TC2 - C2 */
244 /* WRITE (*,*) 'Delta C3 :', TC3 - C3 */
245
246 /* END DO */
247
248 /* END */
249
250 /* You should expect all of the differences to be on the order of */
251 /* the precision of the machine on which this program is executed. */
252
253 /* $ Restrictions */
254
255 /* None. */
256
257 /* $ Literature_References */
258
259 /* [1] `Fundamentals of Celestial Mechanics', Second Edition */
260 /* by J.M.A. Danby; Willman-Bell, Inc., P.O. Box 35025 */
261 /* Richmond Virginia; pp 168-180 */
262
263 /* $ Author_and_Institution */
264
265 /* N.J. Bachman (JPL) */
266 /* H.A. Neilan (JPL) */
267 /* W.L. Taber (JPL) */
268 /* B.V. Semenov (JPL) */
269
270 /* $ Version */
271
272 /* - SPICELIB Version 3.27.0, 08-APR-2014 (NJB) */
273
274 /* Updated in-line documentation and cleaned up */
275 /* code following changes made in version 3.21.0. */
276
277 /* - SPICELIB Version 3.26.0, 10-MAR-2014 (BVS) */
278
279 /* Updated for SUN-SOLARIS-64BIT-INTEL. */
280
281 /* - SPICELIB Version 3.25.0, 10-MAR-2014 (BVS) */
282
283 /* Updated for PC-LINUX-64BIT-IFORT. */
284
285 /* - SPICELIB Version 3.24.0, 10-MAR-2014 (BVS) */
286
287 /* Updated for PC-CYGWIN-GFORTRAN. */
288
289 /* - SPICELIB Version 3.23.0, 10-MAR-2014 (BVS) */
290
291 /* Updated for PC-CYGWIN-64BIT-GFORTRAN. */
292
293 /* - SPICELIB Version 3.22.0, 10-MAR-2014 (BVS) */
294
295 /* Updated for PC-CYGWIN-64BIT-GCC_C. */
296
297 /* - SPICELIB Version 3.21.0, 09-APR-2012 (WLT) */
298
299 /* Code was updated to correct execessive round-off */
300 /* errors in the case where |X| > 1. */
301
302 /* - SPICELIB Version 3.20.0, 13-MAY-2010 (BVS) */
303
304 /* Updated for SUN-SOLARIS-INTEL. */
305
306 /* - SPICELIB Version 3.19.0, 13-MAY-2010 (BVS) */
307
308 /* Updated for SUN-SOLARIS-INTEL-CC_C. */
309
310 /* - SPICELIB Version 3.18.0, 13-MAY-2010 (BVS) */
311
312 /* Updated for SUN-SOLARIS-INTEL-64BIT-CC_C. */
313
314 /* - SPICELIB Version 3.17.0, 13-MAY-2010 (BVS) */
315
316 /* Updated for SUN-SOLARIS-64BIT-NATIVE_C. */
317
318 /* - SPICELIB Version 3.16.0, 13-MAY-2010 (BVS) */
319
320 /* Updated for PC-WINDOWS-64BIT-IFORT. */
321
322 /* - SPICELIB Version 3.15.0, 13-MAY-2010 (BVS) */
323
324 /* Updated for PC-LINUX-64BIT-GFORTRAN. */
325
326 /* - SPICELIB Version 3.14.0, 13-MAY-2010 (BVS) */
327
328 /* Updated for PC-64BIT-MS_C. */
329
330 /* - SPICELIB Version 3.13.0, 13-MAY-2010 (BVS) */
331
332 /* Updated for MAC-OSX-64BIT-INTEL_C. */
333
334 /* - SPICELIB Version 3.12.0, 13-MAY-2010 (BVS) */
335
336 /* Updated for MAC-OSX-64BIT-IFORT. */
337
338 /* - SPICELIB Version 3.11.0, 13-MAY-2010 (BVS) */
339
340 /* Updated for MAC-OSX-64BIT-GFORTRAN. */
341
342 /* - SPICELIB Version 3.10.0, 18-MAR-2009 (BVS) */
343
344 /* Updated for PC-LINUX-GFORTRAN. */
345
346 /* - SPICELIB Version 3.9.0, 18-MAR-2009 (BVS) */
347
348 /* Updated for MAC-OSX-GFORTRAN. */
349
350 /* - SPICELIB Version 3.8.0, 19-FEB-2008 (BVS) */
351
352 /* Updated for PC-LINUX-IFORT. */
353
354 /* - SPICELIB Version 3.7.0, 14-NOV-2006 (BVS) */
355
356 /* Updated for PC-LINUX-64BIT-GCC_C. */
357
358 /* - SPICELIB Version 3.6.0, 14-NOV-2006 (BVS) */
359
360 /* Updated for MAC-OSX-INTEL_C. */
361
362 /* - SPICELIB Version 3.5.0, 14-NOV-2006 (BVS) */
363
364 /* Updated for MAC-OSX-IFORT. */
365
366 /* - SPICELIB Version 3.4.0, 14-NOV-2006 (BVS) */
367
368 /* Updated for PC-WINDOWS-IFORT. */
369
370 /* - SPICELIB Version 3.3.0, 26-OCT-2005 (BVS) */
371
372 /* Updated for SUN-SOLARIS-64BIT-GCC_C. */
373
374 /* - SPICELIB Version 3.2.0, 03-JAN-2005 (BVS) */
375
376 /* Updated for PC-CYGWIN_C. */
377
378 /* - SPICELIB Version 3.1.0, 03-JAN-2005 (BVS) */
379
380 /* Updated for PC-CYGWIN. */
381
382 /* - SPICELIB Version 3.0.5, 17-JUL-2002 (BVS) */
383
384 /* Added MAC-OSX environments. */
385
386 /* - SPICELIB Version 3.0.4, 08-OCT-1999 (WLT) */
387
388 /* The environment lines were expanded so that the supported */
389 /* environments are now explicitely given. New */
390 /* environments are WIN-NT */
391
392 /* - SPICELIB Version 3.0.3, 24-SEP-1999 (NJB) */
393
394 /* CSPICE environments were added. Some typos were corrected. */
395
396 /* - SPICELIB Version 3.0.2, 28-JUL-1999 (WLT) */
397
398 /* The environment lines were expanded so that the supported */
399 /* environments are now explicitly given. New */
400 /* environments are PC-DIGITAL, SGI-O32 and SGI-N32. */
401
402 /* - SPICELIB Version 3.0.1, 18-MAR-1999 (WLT) */
403
404 /* The environment lines were expanded so that the supported */
405 /* environments are now explicitly given. Previously, */
406 /* environments such as SUN-SUNOS and SUN-SOLARIS were implied */
407 /* by the environment label SUN. */
408
409 /* - SPICELIB Version 3.0.0, 08-APR-1998 (NJB) */
410
411 /* Module was updated for the PC-LINUX platform. */
412
413 /* - SPICELIB Version 2.0.0, 11-NOV-1993 (HAN) */
414
415 /* The file was modified to include values for other platforms. */
416 /* Also, the file was formatted for use by the program that */
417 /* creates the environment specific source files. */
418
419 /* - SPICELIB Version 1.0.0, 17-FEB-1992 (WLT) */
420
421 /* -& */
422 /* $ Index_Entries */
423
424 /* Evaluate the first four Stumpff functions */
425
426 /* -& */
427 /* $ Revisions */
428
429 /* - SPICELIB Version 3.27.0, 08-APR-2014 (NJB) (WLT) */
430
431 /* In version 3.21.0, the routine was re-written to use the */
432 /* standard trigonometric and hyperbolic trigonometric formulas */
433 /* for the Stumpff functions for input arguments having magnitude */
434 /* greater than or equal to 1. This was done to prevent loss of */
435 /* accuracy for some input values. */
436
437 /* In version 3.27.0, the code was cleaned up: unreachable code */
438 /* was deleted, and comments were changed to match the updated */
439 /* code. */
440
441 /* The derivation of the argument mapping formulas has been */
442 /* retained as an appendix in a comment section at the end of */
443 /* this source file. These formulas may be used a future revision */
444 /* of this routine. */
445
446 /* - SPICELIB Version 3.0.0, 08-APR-1998 (NJB) */
447
448 /* Module was updated for the PC-LINUX platform. */
449
450 /* - SPICELIB Version 2.0.0, 11-NOV-1993 (HAN) */
451
452 /* The file was modified to include values for other platforms. */
453 /* Also, the file was formatted for use by the program that */
454 /* creates the environment specific source files. */
455
456 /* -& */
457
458 /* SPICELIB Functions */
459
460
461 /* Local Parameters */
462
463
464 /* The integers NPAIRS, LPAIR2, and LPAIR3 are used to declare */
465 /* space for Maclaurin series coefficients and for determining how */
466 /* many terms of these series to use in the computation of */
467 /* C_2 and C_3. */
468
469 /* Here's what is supposed to be true. */
470
471 /* 1/(TRUNC*2)! + 1.0D0 = 1.0D0 */
472
473 /* using this machine's double precision arithmetic. */
474
475 /* We will map the input X to a value y between -1 and 1 and then */
476 /* construct the values of the functions at X from their values at y. */
477 /* Since we will only evaluate the series expansion for C_2 and C_3 */
478 /* for values of y between -1 and 1, its easy to show that we don't */
479 /* need to consider terms in the series whose coefficients have */
480 /* magnitudes less than or equal 1/(2*TRUNC)! . */
481
482 /* If the value of TRUNC is 10, then the series expansions for */
483 /* C_2(y) and C_3(y) are can be truncated as shown here: */
484
485 /* 2 7 8 */
486 /* . 1 y y y y */
487 /* C_3(y) = --- - --- + --- + ... - --- + --- */
488 /* 3! 5! 7! 17! 19! */
489
490
491 /* 1 y y y y */
492 /* = ---( 1 - --- ( 1 - --- (...( 1- ----- ( 1 - ----- )...) */
493 /* 2*3 4*5 6*7 16*17 18*19 */
494
495
496
497
498 /* . 1 y y y y */
499 /* C_2(y) = --- - --- + --- + ... + --- - --- */
500 /* 2! 4! 6! 16! 18! */
501
502
503 /* 1 y y y y */
504 /* = ---( 1 - --- ( 1 - --- (...( 1- ----- ( 1 - ----- )...) */
505 /* 1*2 3*4 5*6 15*16 17*18 */
506
507 /* As is evident from the above, we are going to need the */
508 /* "reciprocal pairs" */
509
510 /* 1/(1*2), 1/(2*3), 1/(3*4), 1/(4*5), ... */
511
512 /* The number of such fractions be computed directly from */
513 /* TRUNC. LPAIR3 and LPAIR2 indicate which of these pairs */
514 /* (counting 1/(1*2) as the first) will be the last one needed in */
515 /* the evaluation of C_2 and C_3. */
516
517
518 /* Local variables */
519
520
521 /* Saved variables */
522
523
524 /* Initial values */
525
526
527 /* We are going to need the numbers */
528
529 /* 1/(2*3), 1/(3*4), 1/(4*5), ... */
530
531 /* but we don't want to compute them every time this routine is */
532 /* called. So the first time this routine is called we compute */
533 /* them and put them in the array PAIRS for use on subsequent */
534 /* calls. (This could be done via parameters, but computing them */
535 /* at run time seems to have a better chance of being */
536 /* easily maintained.) */
537
538 /* In addition we will need to compute the lower bound for which */
539 /* C_0,...,C_3 can be computed. This lower bound is computed by */
540 /* noting that C_0 has the largest magnitude of all the Stumpff */
541 /* functions over the domain from -infinity to -1. Moreover, in this */
542 /* range */
543
544 /* C_0(X) = Cosh( SQRT(-X) ) */
545
546 /* Thus the range of X for which the Stumpff functions can be */
547 /* computed is bounded below by the value of X for which */
548
549 /* Cosh ( SQRT(-X) ) = DPMAX */
550
551 /* Which implies the lower bound for valid inputs is at */
552
553 /* X = - ( DLOG ( 2.0 ) + DLOG( DPMAX ) ) ** 2 */
554
555 /* = - ( DLOG ( 2*N ) + DLOG ( DPMAX/N ) ) ** 2 */
556
557 /* We point out the second formulation of the bound just in case */
558 /* your compiler can't handle the computation of DLOG ( DPMAX ). */
559 /* If this unfortunate situation should arise, complain to the */
560 /* company that produces your compiler and in the code below */
561 /* compute LBOUND using the second form above with N equal to */
562 /* some large power of 2 (say 2**20). */
563
564 if (first) {
565 first = FALSE_;
566 for (i__ = 1; i__ <= 20; ++i__) {
567 pairs[(i__1 = i__ - 1) < 20 && 0 <= i__1 ? i__1 : s_rnge("pairs",
568 i__1, "stmp03_", (ftnlen)589)] = 1. / ((doublereal) i__ *
569 (doublereal) (i__ + 1));
570 }
571 y = log(2.) + log(dpmax_());
572 lbound = -y * y;
573 }
574
575 /* First we make sure that the input value of X is within the */
576 /* range that we are confident we can use to compute the Stumpff */
577 /* functions. */
578
579 if (*x <= lbound) {
580 chkin_("STMP03", (ftnlen)6);
581 setmsg_("The input value of X must be greater than #. The input val"
582 "ue was #", (ftnlen)67);
583 errdp_("#", &lbound, (ftnlen)1);
584 errdp_("#", x, (ftnlen)1);
585 sigerr_("SPICE(VALUEOUTOFRANGE)", (ftnlen)22);
586 chkout_("STMP03", (ftnlen)6);
587 return 0;
588 }
589
590 /* From the definition of the Stumpff functions it can be seen that */
591 /* C_0(X), C_1(X) are given by */
592
593 /* COS ( DSQRT(X) ) and SIN ( DSQRT(X) ) / DSQRT(X) */
594
595 /* for positive X. Moreover, the series used to define them converges */
596 /* for all real X. */
597
598 /* These functions have a number of simple relationships that make */
599 /* their computations practical. Among these are: */
600
601 /* 1 */
602 /* x*C_k+2(x) = --- - C_k(x) */
603 /* k! */
604
605
606
607 /* 2 */
608 /* C_0(4x) = 2*[ C_0(x) ] - 1 */
609
610
611
612
613 /* C_1(4x) = C_1(x)*C_0(x) */
614
615
616
617 /* 2 */
618 /* C_2(4x) = [C_1(x)] / 2 */
619
620
621
622
623 /* C_3(4x) = [ C_2(x) + C_0(x)*C_3(x) ] / 4 */
624
625 /* These can be used to derive formulae for C_0(16x) ... C_3(16x) */
626 /* that involve only C_0(x) ... C_3(x). If we let */
627
628 /* 2 */
629 /* Z = C_0(x) - 0.5 */
630
631 /* and */
632
633 /* W = 2*C_0(x)*C_1(x) */
634
635 /* then */
636
637 /* 2 */
638 /* C_0(16x) = 8*Z - 1 */
639
640
641 /* C_1(16x) = W*Z */
642
643
644 /* 2 */
645 /* C_2(16x) = W / 8 */
646
647
648 /* 2 */
649 /* C_1(x) + Z*[C_2(x) + C_0(x)*C_3(x)] */
650 /* C_3(16x) = ---------------------------------- */
651 /* 8 */
652
653
654 if (*x < -1.) {
655 z__ = sqrt(-(*x));
656 *c0 = cosh(z__);
657 *c1 = sinh(z__) / z__;
658 *c2 = (1 - *c0) / *x;
659 *c3 = (1 - *c1) / *x;
660 return 0;
661 }
662 if (*x > 1.) {
663 z__ = sqrt(*x);
664 *c0 = cos(z__);
665 *c1 = sin(z__) / z__;
666 *c2 = (1 - *c0) / *x;
667 *c3 = (1 - *c1) / *x;
668 return 0;
669 }
670
671 /* If the magnitude of X is less than or equal to 1, we compute */
672 /* the function values directly from their power series */
673 /* representations. */
674
675
676 /* Compute C_3 of x : */
677
678 /* 2 7 8 */
679 /* . 1 x x x x */
680 /* C_3(x) = --- - --- + --- + ... - --- + --- */
681 /* 3! 5! 7! 17! 19! */
682
683
684 /* 1 x x x x */
685 /* = ---( 1 - --- ( 1 - --- (...( 1- ----- ( 1 - ----- )...) */
686 /* 2*3 4*5 6*7 16*17 18*19 */
687
688 /* ^ ^ ^ ^ ^ */
689 /* | | | | | */
690 /* | | | | | */
691 /* PAIR(2) PAIR(4) PAIR(6) ... PAIR(16) PAIR(18) */
692
693 /* Assuming that we don't need to go beyond the term with 1/19!, */
694 /* LPAIR3 will be 18. */
695
696 *c3 = 1.;
697 for (i__ = 20; i__ >= 4; i__ += -2) {
698 *c3 = 1. - *x * pairs[(i__1 = i__ - 1) < 20 && 0 <= i__1 ? i__1 :
699 s_rnge("pairs", i__1, "stmp03_", (ftnlen)733)] * *c3;
700 }
701 *c3 = pairs[1] * *c3;
702
703 /* Compute C_2 of x : */
704
705 /* Here's how we do it. */
706 /* 2 7 8 */
707 /* . 1 x x x x */
708 /* C_2(x) = --- - --- + --- + ... + --- - --- */
709 /* 2! 4! 6! 16! 18! */
710
711
712 /* 1 x x x x */
713 /* = ---( 1 - --- ( 1 - --- (...( 1- ----- ( 1 - ----- )...) */
714 /* 1*2 3*4 5*6 15*16 17*18 */
715
716 /* ^ ^ ^ ^ ^ */
717 /* | | | | | */
718 /* | | | | | */
719 /* PAIR(1) PAIR(3) PAIR(5) ... PAIR(15) PAIR(17) */
720
721 /* Assuming that we don't need to go beyond the term with 1/18!, */
722 /* LPAIR2 will be 17. */
723
724 *c2 = 1.;
725 for (i__ = 19; i__ >= 3; i__ += -2) {
726 *c2 = 1. - *x * pairs[(i__1 = i__ - 1) < 20 && 0 <= i__1 ? i__1 :
727 s_rnge("pairs", i__1, "stmp03_", (ftnlen)764)] * *c2;
728 }
729 *c2 = pairs[0] * *c2;
730
731 /* Get C1 and C0 via the recursion formula: */
732
733 /* 1 */
734 /* x*C_k+2(y) = --- - C_k(x) */
735 /* k! */
736
737 *c1 = 1. - *x * *c3;
738 *c0 = 1. - *x * *c2;
739 return 0;
740 } /* stmp03_ */
741
742