NGSolve 5.3
bessel.hpp
1/* j0.c
2 *
3 * Bessel function of order zero
4 *
5 * thanks to Thorsten Hohage
6 *
7 *
8 * SYNOPSIS:
9 *
10 * double x, y, j0();
11 *
12 * y = j0( x );
13 *
14 *
15 *
16 * DESCRIPTION:
17 *
18 * Returns Bessel function of order zero of the argument.
19 *
20 * The domain is divided into the intervals [0, 5] and
21 * (5, infinity). In the first interval the following rational
22 * approximation is used:
23 *
24 *
25 * 2 2
26 * (w - r ) (w - r ) P (w) / Q (w)
27 * 1 2 3 8
28 *
29 * 2
30 * where w = x and the two r's are zeros of the function.
31 *
32 * In the second interval, the Hankel asymptotic expansion
33 * is employed with two rational functions of degree 6/6
34 * and 7/7.
35 *
36 *
37 *
38 * ACCURACY:
39 *
40 * Absolute error:
41 * arithmetic domain # trials peak rms
42 * DEC 0, 30 10000 4.4e-17 6.3e-18
43 * IEEE 0, 30 60000 4.2e-16 1.1e-16
44 *
45 */
46/* y0.c
47 *
48 * Bessel function of the second kind, order zero
49 *
50 *
51 *
52 * SYNOPSIS:
53 *
54 * double x, y, y0();
55 *
56 * y = y0( x );
57 *
58 *
59 *
60 * DESCRIPTION:
61 *
62 * Returns Bessel function of the second kind, of order
63 * zero, of the argument.
64 *
65 * The domain is divided into the intervals [0, 5] and
66 * (5, infinity). In the first interval a rational approximation
67 * R(x) is employed to compute
68 * y0(x) = R(x) + 2 * log(x) * j0(x) / PI.
69 * Thus a call to j0() is required.
70 *
71 * In the second interval, the Hankel asymptotic expansion
72 * is employed with two rational functions of degree 6/6
73 * and 7/7.
74 *
75 *
76 *
77 * ACCURACY:
78 *
79 * Absolute error, when y0(x) < 1; else relative error:
80 *
81 * arithmetic domain # trials peak rms
82 * DEC 0, 30 9400 7.0e-17 7.9e-18
83 * IEEE 0, 30 30000 1.3e-15 1.6e-16
84 *
85 */
86/*
87Cephes Math Library Release 2.1: January, 1989
88Copyright 1984, 1987, 1989 by Stephen L. Moshier
89Direct inquiries to 30 Frost Street, Cambridge, MA 02140
90*/
91
92/* Note: all coefficients satisfy the relative error criterion
93 * except YP, YQ which are designed for absolute error. */
94
95//#include "mconf.h"
96#define UNK
97#ifdef UNK
98static double PP[7] = {
99 7.96936729297347051624E-4,
100 8.28352392107440799803E-2,
101 1.23953371646414299388E0,
102 5.44725003058768775090E0,
103 8.74716500199817011941E0,
104 5.30324038235394892183E0,
105 9.99999999999999997821E-1,
106};
107static double PQ[7] = {
108 9.24408810558863637013E-4,
109 8.56288474354474431428E-2,
110 1.25352743901058953537E0,
111 5.47097740330417105182E0,
112 8.76190883237069594232E0,
113 5.30605288235394617618E0,
114 1.00000000000000000218E0,
115};
116#endif
117#ifdef DEC
118static unsigned short PP[28] = {
1190035520,0164604,0140733,0054470,
1200037251,0122605,0115356,0107170,
1210040236,0124412,0071500,0056303,
1220040656,0047737,0045720,0045263,
1230041013,0172143,0045004,0142103,
1240040651,0132045,0026241,0026406,
1250040200,0000000,0000000,0000000,
126};
127static unsigned short PQ[28] = {
1280035562,0052006,0070034,0134666,
1290037257,0057055,0055242,0123424,
1300040240,0071626,0046630,0032371,
1310040657,0011077,0032013,0012731,
1320041014,0030307,0050331,0006414,
1330040651,0145457,0065021,0150304,
1340040200,0000000,0000000,0000000,
135};
136#endif
137#ifdef IBMPC
138static unsigned short PP[28] = {
1390x6b27,0x983b,0x1d30,0x3f4a,
1400xd1cf,0xb35d,0x34b0,0x3fb5,
1410x0b98,0x4e68,0xd521,0x3ff3,
1420x0956,0xe97a,0xc9fb,0x4015,
1430x9888,0x6940,0x7e8c,0x4021,
1440x25a1,0xa594,0x3684,0x4015,
1450x0000,0x0000,0x0000,0x3ff0,
146};
147static unsigned short PQ[28] = {
1480x9737,0xce03,0x4a80,0x3f4e,
1490x54e3,0xab54,0xebc5,0x3fb5,
1500x069f,0xc9b3,0x0e72,0x3ff4,
1510x62bb,0xe681,0xe247,0x4015,
1520x21a1,0xea1b,0x8618,0x4021,
1530x3a19,0xed42,0x3965,0x4015,
1540x0000,0x0000,0x0000,0x3ff0,
155};
156#endif
157#ifdef MIEEE
158static unsigned short PP[28] = {
1590x3f4a,0x1d30,0x983b,0x6b27,
1600x3fb5,0x34b0,0xb35d,0xd1cf,
1610x3ff3,0xd521,0x4e68,0x0b98,
1620x4015,0xc9fb,0xe97a,0x0956,
1630x4021,0x7e8c,0x6940,0x9888,
1640x4015,0x3684,0xa594,0x25a1,
1650x3ff0,0x0000,0x0000,0x0000,
166};
167static unsigned short PQ[28] = {
1680x3f4e,0x4a80,0xce03,0x9737,
1690x3fb5,0xebc5,0xab54,0x54e3,
1700x3ff4,0x0e72,0xc9b3,0x069f,
1710x4015,0xe247,0xe681,0x62bb,
1720x4021,0x8618,0xea1b,0x21a1,
1730x4015,0x3965,0xed42,0x3a19,
1740x3ff0,0x0000,0x0000,0x0000,
175};
176#endif
177
178#ifdef UNK
179static double QP[8] = {
180-1.13663838898469149931E-2,
181-1.28252718670509318512E0,
182-1.95539544257735972385E1,
183-9.32060152123768231369E1,
184-1.77681167980488050595E2,
185-1.47077505154951170175E2,
186-5.14105326766599330220E1,
187-6.05014350600728481186E0,
188};
189static double QQ[7] = {
190/* 1.00000000000000000000E0,*/
191 6.43178256118178023184E1,
192 8.56430025976980587198E2,
193 3.88240183605401609683E3,
194 7.24046774195652478189E3,
195 5.93072701187316984827E3,
196 2.06209331660327847417E3,
197 2.42005740240291393179E2,
198};
199#endif
200#ifdef DEC
201static unsigned short QP[32] = {
2020136472,0035021,0142451,0141115,
2030140244,0024731,0150620,0105642,
2040141234,0067177,0124161,0060141,
2050141672,0064572,0151557,0043036,
2060142061,0127141,0003127,0043517,
2070142023,0011727,0060271,0144544,
2080141515,0122142,0126620,0143150,
2090140701,0115306,0106715,0007344,
210};
211static unsigned short QQ[28] = {
212/*0040200,0000000,0000000,0000000,*/
2130041600,0121272,0004741,0026544,
2140042526,0015605,0105654,0161771,
2150043162,0123155,0165644,0062645,
2160043342,0041675,0167576,0130756,
2170043271,0052720,0165631,0154214,
2180043000,0160576,0034614,0172024,
2190042162,0000570,0030500,0051235,
220};
221#endif
222#ifdef IBMPC
223static unsigned short QP[32] = {
2240x384a,0x38a5,0x4742,0xbf87,
2250x1174,0x3a32,0x853b,0xbff4,
2260x2c0c,0xf50e,0x8dcf,0xc033,
2270xe8c4,0x5a6d,0x4d2f,0xc057,
2280xe8ea,0x20ca,0x35cc,0xc066,
2290x392d,0xec17,0x627a,0xc062,
2300x18cd,0x55b2,0xb48c,0xc049,
2310xa1dd,0xd1b9,0x3358,0xc018,
232};
233static unsigned short QQ[28] = {
234/*0x0000,0x0000,0x0000,0x3ff0,*/
2350x25ac,0x413c,0x1457,0x4050,
2360x9c7f,0xb175,0xc370,0x408a,
2370x8cb5,0xbd74,0x54cd,0x40ae,
2380xd63e,0xbdef,0x4877,0x40bc,
2390x3b11,0x1d73,0x2aba,0x40b7,
2400x9e82,0xc731,0x1c2f,0x40a0,
2410x0a54,0x0628,0x402f,0x406e,
242};
243#endif
244#ifdef MIEEE
245static unsigned short QP[32] = {
2460xbf87,0x4742,0x38a5,0x384a,
2470xbff4,0x853b,0x3a32,0x1174,
2480xc033,0x8dcf,0xf50e,0x2c0c,
2490xc057,0x4d2f,0x5a6d,0xe8c4,
2500xc066,0x35cc,0x20ca,0xe8ea,
2510xc062,0x627a,0xec17,0x392d,
2520xc049,0xb48c,0x55b2,0x18cd,
2530xc018,0x3358,0xd1b9,0xa1dd,
254};
255static unsigned short QQ[28] = {
256/*0x3ff0,0x0000,0x0000,0x0000,*/
2570x4050,0x1457,0x413c,0x25ac,
2580x408a,0xc370,0xb175,0x9c7f,
2590x40ae,0x54cd,0xbd74,0x8cb5,
2600x40bc,0x4877,0xbdef,0xd63e,
2610x40b7,0x2aba,0x1d73,0x3b11,
2620x40a0,0x1c2f,0xc731,0x9e82,
2630x406e,0x402f,0x0628,0x0a54,
264};
265#endif
266
267
268#ifdef UNK
269static double YP[8] = {
270 1.55924367855235737965E4,
271-1.46639295903971606143E7,
272 5.43526477051876500413E9,
273-9.82136065717911466409E11,
274 8.75906394395366999549E13,
275-3.46628303384729719441E15,
276 4.42733268572569800351E16,
277-1.84950800436986690637E16,
278};
279static double YQ[7] = {
280/* 1.00000000000000000000E0,*/
281 1.04128353664259848412E3,
282 6.26107330137134956842E5,
283 2.68919633393814121987E8,
284 8.64002487103935000337E10,
285 2.02979612750105546709E13,
286 3.17157752842975028269E15,
287 2.50596256172653059228E17,
288};
289#endif
290#ifdef DEC
291static unsigned short YP[32] = {
2920043563,0120677,0042264,0046166,
2930146137,0140371,0113444,0042260,
2940050241,0175707,0100502,0063344,
2950152144,0125737,0007265,0164526,
2960053637,0051621,0163035,0060546,
2970155105,0004416,0107306,0060023,
2980056035,0045133,0030132,0000024,
2990155603,0065132,0144061,0131732,
300};
301static unsigned short YQ[28] = {
302/*0040200,0000000,0000000,0000000,*/
3030042602,0024422,0135557,0162663,
3040045030,0155665,0044075,0160135,
3050047200,0035432,0105446,0104005,
3060051240,0167331,0056063,0022743,
3070053223,0127746,0025764,0012160,
3080055064,0044206,0177532,0145545,
3090056536,0111375,0163715,0127201,
310};
311#endif
312#ifdef IBMPC
313static unsigned short YP[32] = {
3140x898f,0xe896,0x7437,0x40ce,
3150x8896,0x32e4,0xf81f,0xc16b,
3160x4cdd,0xf028,0x3f78,0x41f4,
3170xbd2b,0xe1d6,0x957b,0xc26c,
3180xac2d,0x3cc3,0xea72,0x42d3,
3190xcc02,0xd1d8,0xa121,0xc328,
3200x4003,0x660b,0xa94b,0x4363,
3210x367b,0x5906,0x6d4b,0xc350,
322};
323static unsigned short YQ[28] = {
324/*0x0000,0x0000,0x0000,0x3ff0,*/
3250xfcb6,0x576d,0x4522,0x4090,
3260xbc0c,0xa907,0x1b76,0x4123,
3270xd101,0x5164,0x0763,0x41b0,
3280x64bc,0x2b86,0x1ddb,0x4234,
3290x828e,0xc57e,0x75fc,0x42b2,
3300x596d,0xdfeb,0x8910,0x4326,
3310xb5d0,0xbcf9,0xd25f,0x438b,
332};
333#endif
334#ifdef MIEEE
335static unsigned short YP[32] = {
3360x40ce,0x7437,0xe896,0x898f,
3370xc16b,0xf81f,0x32e4,0x8896,
3380x41f4,0x3f78,0xf028,0x4cdd,
3390xc26c,0x957b,0xe1d6,0xbd2b,
3400x42d3,0xea72,0x3cc3,0xac2d,
3410xc328,0xa121,0xd1d8,0xcc02,
3420x4363,0xa94b,0x660b,0x4003,
3430xc350,0x6d4b,0x5906,0x367b,
344};
345static unsigned short YQ[28] = {
346/*0x3ff0,0x0000,0x0000,0x0000,*/
3470x4090,0x4522,0x576d,0xfcb6,
3480x4123,0x1b76,0xa907,0xbc0c,
3490x41b0,0x0763,0x5164,0xd101,
3500x4234,0x1ddb,0x2b86,0x64bc,
3510x42b2,0x75fc,0xc57e,0x828e,
3520x4326,0x8910,0xdfeb,0x596d,
3530x438b,0xd25f,0xbcf9,0xb5d0,
354};
355#endif
356
357#ifdef UNK
358/* 5.783185962946784521175995758455807035071 */
359static double DR1 = 5.78318596294678452118E0;
360/* 30.47126234366208639907816317502275584842 */
361static double DR2 = 3.04712623436620863991E1;
362#endif
363
364#ifdef DEC
365static unsigned short R1[] = {0040671,0007734,0001061,0056734};
366#define DR1 *(double *)R1
367static unsigned short R2[] = {0041363,0142445,0030416,0165567};
368#define DR2 *(double *)R2
369#endif
370
371#ifdef IBMPC
372static unsigned short R1[] = {0x2bbb,0x8046,0x21fb,0x4017};
373#define DR1 *(double *)R1
374static unsigned short R2[] = {0xdd6f,0xa621,0x78a4,0x403e};
375#define DR2 *(double *)R2
376#endif
377
378#ifdef MIEEE
379static unsigned short R1[] = {0x4017,0x21fb,0x8046,0x2bbb};
380#define DR1 *(double *)R1
381static unsigned short R2[] = {0x403e,0x78a4,0xa621,0xdd6f};
382#define DR2 *(double *)R2
383#endif
384
385#ifdef UNK
386static double RP[4] = {
387-4.79443220978201773821E9,
388 1.95617491946556577543E12,
389-2.49248344360967716204E14,
390 9.70862251047306323952E15,
391};
392static double RQ[8] = {
393/* 1.00000000000000000000E0,*/
394 4.99563147152651017219E2,
395 1.73785401676374683123E5,
396 4.84409658339962045305E7,
397 1.11855537045356834862E10,
398 2.11277520115489217587E12,
399 3.10518229857422583814E14,
400 3.18121955943204943306E16,
401 1.71086294081043136091E18,
402};
403#endif
404#ifdef DEC
405static unsigned short RP[16] = {
4060150216,0161235,0064344,0014450,
4070052343,0135216,0035624,0144153,
4080154142,0130247,0003310,0003667,
4090055411,0173703,0047772,0176635,
410};
411static unsigned short RQ[32] = {
412/*0040200,0000000,0000000,0000000,*/
4130042371,0144025,0032265,0136137,
4140044451,0133131,0132420,0151466,
4150046470,0144641,0072540,0030636,
4160050446,0126600,0045042,0044243,
4170052365,0172633,0110301,0071063,
4180054215,0032424,0062272,0043513,
4190055742,0005013,0171731,0072335,
4200057275,0170646,0036663,0013134,
421};
422#endif
423#ifdef IBMPC
424static unsigned short RP[16] = {
4250x8325,0xad1c,0xdc53,0xc1f1,
4260x990d,0xc772,0x7751,0x427c,
4270x00f7,0xe0d9,0x5614,0xc2ec,
4280x5fb4,0x69ff,0x3ef8,0x4341,
429};
430static unsigned short RQ[32] = {
431/*0x0000,0x0000,0x0000,0x3ff0,*/
4320xb78c,0xa696,0x3902,0x407f,
4330x1a67,0x36a2,0x36cb,0x4105,
4340x0634,0x2eac,0x1934,0x4187,
4350x4914,0x0944,0xd5b0,0x4204,
4360x2e46,0x7218,0xbeb3,0x427e,
4370x48e9,0x8c97,0xa6a2,0x42f1,
4380x2e9c,0x7e7b,0x4141,0x435c,
4390x62cc,0xc7b6,0xbe34,0x43b7,
440};
441#endif
442#ifdef MIEEE
443static unsigned short RP[16] = {
4440xc1f1,0xdc53,0xad1c,0x8325,
4450x427c,0x7751,0xc772,0x990d,
4460xc2ec,0x5614,0xe0d9,0x00f7,
4470x4341,0x3ef8,0x69ff,0x5fb4,
448};
449static unsigned short RQ[32] = {
450/*0x3ff0,0x0000,0x0000,0x0000,*/
4510x407f,0x3902,0xa696,0xb78c,
4520x4105,0x36cb,0x36a2,0x1a67,
4530x4187,0x1934,0x2eac,0x0634,
4540x4204,0xd5b0,0x0944,0x4914,
4550x427e,0xbeb3,0x7218,0x2e46,
4560x42f1,0xa6a2,0x8c97,0x48e9,
4570x435c,0x4141,0x7e7b,0x2e9c,
4580x43b7,0xbe34,0xc7b6,0x62cc,
459};
460#endif
461
462#define ANSIPROT
463#ifndef ANSIPROT
464double bessj0(), polevl(), p1evl(), log(), sin(), cos(), sqrt();
465#endif
466
467double polevl(double x,double coef[],int N)
468{
469double ans;
470int i;
471double *p;
472
473p = coef;
474ans = *p++;
475i = N;
476
477do
478 ans = ans * x + *p++;
479while( --i );
480
481return( ans );
482}
483
484/* p1evl() */
485/* N
486 * Evaluate polynomial when coefficient of x is 1.0.
487 * Otherwise same as polevl.
488 */
489
490double p1evl(double x, double coef[],int N)
491{
492double ans;
493double *p;
494int i;
495
496p = coef;
497ans = x + *p++;
498i = N-1;
499
500do
501 ans = ans * x + *p++;
502while( --i );
503
504return( ans );
505}
506
507double TWOOPI = 6.36619772367581343075535E-1; /* 2/pi */
508double THPIO4 = 2.35619449019234492885; /* 3*pi/4 */
509double SQ2OPI = 7.9788456080286535587989E-1; /* sqrt( 2/pi ) */
510double PIO4 = 7.85398163397448309616E-1; /* pi/4 */
511
512double bessj0(double x)
513{
514double w, z, p, q, xn;
515
516if( x < 0 )
517 x = -x;
518
519if( x <= 5.0 )
520 {
521 z = x * x;
522 if( x < 1.0e-5 )
523 return( 1.0 - z/4.0 );
524
525 p = (z - DR1) * (z - DR2);
526 p = p * polevl( z, RP, 3)/p1evl( z, RQ, 8 );
527 return( p );
528 }
529
530w = 5.0/x;
531q = 25.0/(x*x);
532p = polevl( q, PP, 6)/polevl( q, PQ, 6 );
533q = polevl( q, QP, 7)/p1evl( q, QQ, 7 );
534xn = x - PIO4;
535p = p * cos(xn) - w * q * sin(xn);
536return( p * SQ2OPI / sqrt(x) );
537}
538/* bessy0() 2 */
539/* Bessel function of second kind, order zero */
540
541/* Rational approximation coefficients YP[], YQ[] are used here.
542 * The function computed is bessy0(x) - 2 * log(x) * bessj0(x) / PI,
543 * whose value at x = 0 is 2 * ( log(0.5) + EUL ) / PI
544 * = 0.073804295108687225.
545 */
546
547
548double bessy0(double x)
549{
550double w, z, p, q, xn;
551
552if( x <= 5.0 )
553 {
554 if( x <= 0.0 )
555 throw Exception ("arg<=0 in bessy0");
556 z = x * x;
557 w = polevl( z, YP, 7) / p1evl( z, YQ, 7 );
558 w += TWOOPI * log(x) * bessj0(x);
559 return( w );
560 }
561
562w = 5.0/x;
563z = 25.0 / (x * x);
564p = polevl( z, PP, 6)/polevl( z, PQ, 6 );
565q = polevl( z, QP, 7)/p1evl( z, QQ, 7 );
566xn = x - PIO4;
567p = p * sin(xn) + w * q * cos(xn);
568return( p * SQ2OPI / sqrt(x) );
569}
570
571/* j1.c
572 *
573 * Bessel function of order one
574 *
575 *
576 *
577 * SYNOPSIS:
578 *
579 * double x, y, j1();
580 *
581 * y = j1( x );
582 *
583 *
584 *
585 * DESCRIPTION:
586 *
587 * Returns Bessel function of order one of the argument.
588 *
589 * The domain is divided into the intervals [0, 8] and
590 * (8, infinity). In the first interval a 24 term Chebyshev
591 * expansion is used. In the second, the asymptotic
592 * trigonometric representation is employed using two
593 * rational functions of degree 5/5.
594 *
595 *
596 *
597 * ACCURACY:
598 *
599 * Absolute error:
600 * arithmetic domain # trials peak rms
601 * DEC 0, 30 10000 4.0e-17 1.1e-17
602 * IEEE 0, 30 30000 2.6e-16 1.1e-16
603 *
604 *
605 */
606/* y1.c
607 *
608 * Bessel function of second kind of order one
609 *
610 *
611 *
612 * SYNOPSIS:
613 *
614 * double x, y, y1();
615 *
616 * y = y1( x );
617 *
618 *
619 *
620 * DESCRIPTION:
621 *
622 * Returns Bessel function of the second kind of order one
623 * of the argument.
624 *
625 * The domain is divided into the intervals [0, 8] and
626 * (8, infinity). In the first interval a 25 term Chebyshev
627 * expansion is used, and a call to j1() is required.
628 * In the second, the asymptotic trigonometric representation
629 * is employed using two rational functions of degree 5/5.
630 *
631 *
632 *
633 * ACCURACY:
634 *
635 * Absolute error:
636 * arithmetic domain # trials peak rms
637 * DEC 0, 30 10000 8.6e-17 1.3e-17
638 * IEEE 0, 30 30000 1.0e-15 1.3e-16
639 *
640 * (error criterion relative when |y1| > 1).
641 *
642 */
643/*
644Cephes Math Library Release 2.1: January, 1989
645Copyright 1984, 1987, 1989 by Stephen L. Moshier
646Direct inquiries to 30 Frost Street, Cambridge, MA 02140
647*/
648
649#define PIO4 .78539816339744830962
650#define THPIO4 2.35619449019234492885
651#define SQ2OPI .79788456080286535588
652
653// #include "mconf.h"
654
655#ifdef UNK
656static double RP1[4] = {
657-8.99971225705559398224E8,
658 4.52228297998194034323E11,
659-7.27494245221818276015E13,
660 3.68295732863852883286E15,
661};
662static double RQ1[8] = {
663/* 1.00000000000000000000E0,*/
664 6.20836478118054335476E2,
665 2.56987256757748830383E5,
666 8.35146791431949253037E7,
667 2.21511595479792499675E10,
668 4.74914122079991414898E12,
669 7.84369607876235854894E14,
670 8.95222336184627338078E16,
671 5.32278620332680085395E18,
672};
673#endif
674#ifdef DEC
675static unsigned short RP1[16] = {
6760147526,0110742,0063322,0077052,
6770051722,0112720,0065034,0061530,
6780153604,0052227,0033147,0105650,
6790055121,0055025,0032276,0022015,
680};
681static unsigned short RQ1[32] = {
682/*0040200,0000000,0000000,0000000,*/
6830042433,0032610,0155604,0033473,
6840044572,0173320,0067270,0006616,
6850046637,0045246,0162225,0006606,
6860050645,0004773,0157577,0053004,
6870052612,0033734,0001667,0176501,
6880054462,0054121,0173147,0121367,
6890056237,0002777,0121451,0176007,
6900057623,0136253,0131601,0044710,
691};
692#endif
693#ifdef IBMPC
694static unsigned short RP1[16] = {
6950x4fc5,0x4cda,0xd23c,0xc1ca,
6960x8c6b,0x0d43,0x52ba,0x425a,
6970xf175,0xe6cc,0x8a92,0xc2d0,
6980xc482,0xa697,0x2b42,0x432a,
699};
700static unsigned short RQ1[32] = {
701/*0x0000,0x0000,0x0000,0x3ff0,*/
7020x86e7,0x1b70,0x66b1,0x4083,
7030x01b2,0x0dd7,0x5eda,0x410f,
7040xa1b1,0xdc92,0xe954,0x4193,
7050xeac1,0x7bef,0xa13f,0x4214,
7060xffa8,0x8076,0x46fb,0x4291,
7070xf45f,0x3ecc,0x4b0a,0x4306,
7080x3f81,0xf465,0xe0bf,0x4373,
7090x2939,0x7670,0x7795,0x43d2,
710};
711#endif
712#ifdef MIEEE
713static unsigned short RP1[16] = {
7140xc1ca,0xd23c,0x4cda,0x4fc5,
7150x425a,0x52ba,0x0d43,0x8c6b,
7160xc2d0,0x8a92,0xe6cc,0xf175,
7170x432a,0x2b42,0xa697,0xc482,
718};
719static unsigned short RQ1[32] = {
720/*0x3ff0,0x0000,0x0000,0x0000,*/
7210x4083,0x66b1,0x1b70,0x86e7,
7220x410f,0x5eda,0x0dd7,0x01b2,
7230x4193,0xe954,0xdc92,0xa1b1,
7240x4214,0xa13f,0x7bef,0xeac1,
7250x4291,0x46fb,0x8076,0xffa8,
7260x4306,0x4b0a,0x3ecc,0xf45f,
7270x4373,0xe0bf,0xf465,0x3f81,
7280x43d2,0x7795,0x7670,0x2939,
729};
730#endif
731
732#ifdef UNK
733static double PP1[7] = {
734 7.62125616208173112003E-4,
735 7.31397056940917570436E-2,
736 1.12719608129684925192E0,
737 5.11207951146807644818E0,
738 8.42404590141772420927E0,
739 5.21451598682361504063E0,
740 1.00000000000000000254E0,
741};
742static double PQ1[7] = {
743 5.71323128072548699714E-4,
744 6.88455908754495404082E-2,
745 1.10514232634061696926E0,
746 5.07386386128601488557E0,
747 8.39985554327604159757E0,
748 5.20982848682361821619E0,
749 9.99999999999999997461E-1,
750};
751#endif
752#ifdef DEC
753static unsigned short PP1[28] = {
7540035507,0144542,0061543,0024326,
7550037225,0145105,0017766,0022661,
7560040220,0043766,0010254,0133255,
7570040643,0113047,0142611,0151521,
7580041006,0144344,0055351,0074261,
7590040646,0156520,0120574,0006416,
7600040200,0000000,0000000,0000000,
761};
762static unsigned short PQ1[28] = {
7630035425,0142330,0115041,0165514,
7640037214,0177352,0145105,0052026,
7650040215,0072515,0141207,0073255,
7660040642,0056427,0137222,0106405,
7670041006,0062716,0166427,0165450,
7680040646,0133352,0035425,0123304,
7690040200,0000000,0000000,0000000,
770};
771#endif
772#ifdef IBMPC
773static unsigned short PP1[28] = {
7740x651b,0x4c6c,0xf92c,0x3f48,
7750xc4b6,0xa3fe,0xb948,0x3fb2,
7760x96d6,0xc215,0x08fe,0x3ff2,
7770x3a6a,0xf8b1,0x72c4,0x4014,
7780x2f16,0x8b5d,0xd91c,0x4020,
7790x81a2,0x142f,0xdbaa,0x4014,
7800x0000,0x0000,0x0000,0x3ff0,
781};
782static unsigned short PQ1[28] = {
7830x3d69,0x1344,0xb89b,0x3f42,
7840xaa83,0x5948,0x9fdd,0x3fb1,
7850xeed6,0xb850,0xaea9,0x3ff1,
7860x51a1,0xf7d2,0x4ba2,0x4014,
7870xfd65,0xdda2,0xccb9,0x4020,
7880xb4d9,0x4762,0xd6dd,0x4014,
7890x0000,0x0000,0x0000,0x3ff0,
790};
791#endif
792#ifdef MIEEE
793static unsigned short PP1[28] = {
7940x3f48,0xf92c,0x4c6c,0x651b,
7950x3fb2,0xb948,0xa3fe,0xc4b6,
7960x3ff2,0x08fe,0xc215,0x96d6,
7970x4014,0x72c4,0xf8b1,0x3a6a,
7980x4020,0xd91c,0x8b5d,0x2f16,
7990x4014,0xdbaa,0x142f,0x81a2,
8000x3ff0,0x0000,0x0000,0x0000,
801};
802static unsigned short PQ1[28] = {
8030x3f42,0xb89b,0x1344,0x3d69,
8040x3fb1,0x9fdd,0x5948,0xaa83,
8050x3ff1,0xaea9,0xb850,0xeed6,
8060x4014,0x4ba2,0xf7d2,0x51a1,
8070x4020,0xccb9,0xdda2,0xfd65,
8080x4014,0xd6dd,0x4762,0xb4d9,
8090x3ff0,0x0000,0x0000,0x0000,
810};
811#endif
812
813#ifdef UNK
814static double QP1[8] = {
815 5.10862594750176621635E-2,
816 4.98213872951233449420E0,
817 7.58238284132545283818E1,
818 3.66779609360150777800E2,
819 7.10856304998926107277E2,
820 5.97489612400613639965E2,
821 2.11688757100572135698E2,
822 2.52070205858023719784E1,
823};
824static double QQ1[7] = {
825/* 1.00000000000000000000E0,*/
826 7.42373277035675149943E1,
827 1.05644886038262816351E3,
828 4.98641058337653607651E3,
829 9.56231892404756170795E3,
830 7.99704160447350683650E3,
831 2.82619278517639096600E3,
832 3.36093607810698293419E2,
833};
834#endif
835#ifdef DEC
836static unsigned short QP1[32] = {
8370037121,0037723,0055605,0151004,
8380040637,0066656,0031554,0077264,
8390041627,0122714,0153170,0161466,
8400042267,0061712,0036520,0140145,
8410042461,0133315,0131573,0071176,
8420042425,0057525,0147500,0013201,
8430042123,0130122,0061245,0154131,
8440041311,0123772,0064254,0172650,
845};
846static unsigned short QQ1[28] = {
847/*0040200,0000000,0000000,0000000,*/
8480041624,0074603,0002112,0101670,
8490042604,0007135,0010162,0175565,
8500043233,0151510,0157757,0172010,
8510043425,0064506,0112006,0104276,
8520043371,0164125,0032271,0164242,
8530043060,0121425,0122750,0136013,
8540042250,0005773,0053472,0146267,
855};
856#endif
857#ifdef IBMPC
858static unsigned short QP1[32] = {
8590xba40,0x6b70,0x27fa,0x3faa,
8600x8fd6,0xc66d,0xedb5,0x4013,
8610x1c67,0x9acf,0xf4b9,0x4052,
8620x180d,0x47aa,0xec79,0x4076,
8630x6e50,0xb66f,0x36d9,0x4086,
8640x02d0,0xb9e8,0xabea,0x4082,
8650xbb0b,0x4c54,0x760a,0x406a,
8660x9eb5,0x4d15,0x34ff,0x4039,
867};
868static unsigned short QQ1[28] = {
869/*0x0000,0x0000,0x0000,0x3ff0,*/
8700x5077,0x6089,0x8f30,0x4052,
8710x5f6f,0xa20e,0x81cb,0x4090,
8720xfe81,0x1bfd,0x7a69,0x40b3,
8730xd118,0xd280,0xad28,0x40c2,
8740x3d14,0xa697,0x3d0a,0x40bf,
8750x1781,0xb4bd,0x1462,0x40a6,
8760x5997,0x6ae7,0x017f,0x4075,
877};
878#endif
879#ifdef MIEEE
880static unsigned short QP1[32] = {
8810x3faa,0x27fa,0x6b70,0xba40,
8820x4013,0xedb5,0xc66d,0x8fd6,
8830x4052,0xf4b9,0x9acf,0x1c67,
8840x4076,0xec79,0x47aa,0x180d,
8850x4086,0x36d9,0xb66f,0x6e50,
8860x4082,0xabea,0xb9e8,0x02d0,
8870x406a,0x760a,0x4c54,0xbb0b,
8880x4039,0x34ff,0x4d15,0x9eb5,
889};
890static unsigned short QQ1[28] = {
891/*0x3ff0,0x0000,0x0000,0x0000,*/
8920x4052,0x8f30,0x6089,0x5077,
8930x4090,0x81cb,0xa20e,0x5f6f,
8940x40b3,0x7a69,0x1bfd,0xfe81,
8950x40c2,0xad28,0xd280,0xd118,
8960x40bf,0x3d0a,0xa697,0x3d14,
8970x40a6,0x1462,0xb4bd,0x1781,
8980x4075,0x017f,0x6ae7,0x5997,
899};
900#endif
901
902#ifdef UNK
903static double YP1[6] = {
904 1.26320474790178026440E9,
905-6.47355876379160291031E11,
906 1.14509511541823727583E14,
907-8.12770255501325109621E15,
908 2.02439475713594898196E17,
909-7.78877196265950026825E17,
910};
911static double YQ1[8] = {
912/* 1.00000000000000000000E0,*/
913 5.94301592346128195359E2,
914 2.35564092943068577943E5,
915 7.34811944459721705660E7,
916 1.87601316108706159478E10,
917 3.88231277496238566008E12,
918 6.20557727146953693363E14,
919 6.87141087355300489866E16,
920 3.97270608116560655612E18,
921};
922#endif
923#ifdef DEC
924static unsigned short YP1[24] = {
9250047626,0112763,0013715,0133045,
9260152026,0134552,0142033,0024411,
9270053720,0045245,0102210,0077565,
9280155347,0000321,0136415,0102031,
9290056463,0146550,0055633,0032605,
9300157054,0171012,0167361,0054265,
931};
932static unsigned short YQ1[32] = {
933/*0040200,0000000,0000000,0000000,*/
9340042424,0111515,0044773,0153014,
9350044546,0005405,0171307,0075774,
9360046614,0023575,0047105,0063556,
9370050613,0143034,0101533,0156026,
9380052541,0175367,0166514,0114257,
9390054415,0014466,0134350,0171154,
9400056164,0017436,0025075,0022101,
9410057534,0103614,0103663,0121772,
942};
943#endif
944#ifdef IBMPC
945static unsigned short YP1[24] = {
9460xb6c5,0x62f9,0xd2be,0x41d2,
9470x6521,0x5883,0xd72d,0xc262,
9480x0fef,0xb091,0x0954,0x42da,
9490xb083,0x37a1,0xe01a,0xc33c,
9500x66b1,0x0b73,0x79ad,0x4386,
9510x2b17,0x5dde,0x9e41,0xc3a5,
952};
953static unsigned short YQ1[32] = {
954/*0x0000,0x0000,0x0000,0x3ff0,*/
9550x7ac2,0xa93f,0x9269,0x4082,
9560xef7f,0xbe58,0xc160,0x410c,
9570xacee,0xa9c8,0x84ef,0x4191,
9580x7b83,0x906b,0x78c3,0x4211,
9590x9316,0xfda9,0x3f5e,0x428c,
9600x1e4e,0xd71d,0xa326,0x4301,
9610xa488,0xc547,0x83e3,0x436e,
9620x747f,0x90f6,0x90f1,0x43cb,
963};
964#endif
965#ifdef MIEEE
966static unsigned short YP1[24] = {
9670x41d2,0xd2be,0x62f9,0xb6c5,
9680xc262,0xd72d,0x5883,0x6521,
9690x42da,0x0954,0xb091,0x0fef,
9700xc33c,0xe01a,0x37a1,0xb083,
9710x4386,0x79ad,0x0b73,0x66b1,
9720xc3a5,0x9e41,0x5dde,0x2b17,
973};
974static unsigned short YQ1[32] = {
975/*0x3ff0,0x0000,0x0000,0x0000,*/
9760x4082,0x9269,0xa93f,0x7ac2,
9770x410c,0xc160,0xbe58,0xef7f,
9780x4191,0x84ef,0xa9c8,0xacee,
9790x4211,0x78c3,0x906b,0x7b83,
9800x428c,0x3f5e,0xfda9,0x9316,
9810x4301,0xa326,0xd71d,0x1e4e,
9820x436e,0x83e3,0xc547,0xa488,
9830x43cb,0x90f1,0x90f6,0x747f,
984};
985#endif
986
987
988#ifdef UNK
989static double Z1 = 1.46819706421238932572E1;
990static double Z2 = 4.92184563216946036703E1;
991#endif
992
993#ifdef DEC
994static unsigned short DZ1[] = {0041152,0164532,0006114,0010540};
995static unsigned short DZ2[] = {0041504,0157663,0001625,0020621};
996#define Z1 (*(double *)DZ1)
997#define Z2 (*(double *)DZ2)
998#endif
999
1000#ifdef IBMPC
1001static unsigned short DZ1[] = {0x822c,0x4189,0x5d2b,0x402d};
1002static unsigned short DZ2[] = {0xa432,0x6072,0x9bf6,0x4048};
1003#define Z1 (*(double *)DZ1)
1004#define Z2 (*(double *)DZ2)
1005#endif
1006
1007#ifdef MIEEE
1008static unsigned short DZ1[] = {0x402d,0x5d2b,0x4189,0x822c};
1009static unsigned short DZ2[] = {0x4048,0x9bf6,0x6072,0xa432};
1010#define Z1 (*(double *)DZ1)
1011#define Z2 (*(double *)DZ2)
1012#endif
1013
1014#ifndef ANSIPROT
1015double bessj1(), polevl(), p1evl(), log(), sin(), cos(), sqrt();
1016#endif
1017
1018double bessj1(double x)
1019{
1020double w, z, p, q, xn;
1021
1022w = x;
1023if( x < 0 )
1024 w = -x;
1025
1026if( w <= 5.0 )
1027 {
1028 z = x * x;
1029 w = polevl( z, RP1, 3 ) / p1evl( z, RQ1, 8 );
1030 w = w * x * (z - Z1) * (z - Z2);
1031 return( w );
1032 }
1033
1034w = 5.0/x;
1035z = w * w;
1036p = polevl( z, PP1, 6)/polevl( z, PQ1, 6 );
1037q = polevl( z, QP1, 7)/p1evl( z, QQ1, 7 );
1038xn = x - THPIO4;
1039p = p * cos(xn) - w * q * sin(xn);
1040return( p * SQ2OPI / sqrt(x) );
1041}
1042
1043double bessy1(double x)
1044{
1045double w, z, p, q, xn;
1046
1047if( x <= 5.0 )
1048 {
1049 if( x <= 0.0 )
1050 throw Exception("arg<=0 in bessy1");
1051 z = x * x;
1052 w = x * (polevl( z, YP1, 5 ) / p1evl( z, YQ1, 8 ));
1053 w += TWOOPI * ( bessj1(x) * log(x) - 1.0/x );
1054 return( w );
1055 }
1056
1057w = 5.0/x;
1058z = w * w;
1059p = polevl( z, PP1, 6)/polevl( z, PQ1, 6 );
1060q = polevl( z, QP1, 7)/p1evl( z, QQ1, 7 );
1061xn = x - THPIO4;
1062p = p * sin(xn) + w * q * cos(xn);
1063return( p * SQ2OPI / sqrt(x) );
1064}