NGSolve  5.3
ngstd/bessel.hpp
00001 /*                                                      j0.c
00002  *
00003  *      Bessel function of order zero
00004  *
00005  *      thanks to Thorsten Hohage
00006  *
00007  *
00008  * SYNOPSIS:
00009  *
00010  * double x, y, j0();
00011  *
00012  * y = j0( x );
00013  *
00014  *
00015  *
00016  * DESCRIPTION:
00017  *
00018  * Returns Bessel function of order zero of the argument.
00019  *
00020  * The domain is divided into the intervals [0, 5] and
00021  * (5, infinity). In the first interval the following rational
00022  * approximation is used:
00023  *
00024  *
00025  *        2         2
00026  * (w - r  ) (w - r  ) P (w) / Q (w)
00027  *       1         2    3       8
00028  *
00029  *            2
00030  * where w = x  and the two r's are zeros of the function.
00031  *
00032  * In the second interval, the Hankel asymptotic expansion
00033  * is employed with two rational functions of degree 6/6
00034  * and 7/7.
00035  *
00036  *
00037  *
00038  * ACCURACY:
00039  *
00040  *                      Absolute error:
00041  * arithmetic   domain     # trials      peak         rms
00042  *    DEC       0, 30       10000       4.4e-17     6.3e-18
00043  *    IEEE      0, 30       60000       4.2e-16     1.1e-16
00044  *
00045  */
00046 /*                                                      y0.c
00047  *
00048  *      Bessel function of the second kind, order zero
00049  *
00050  *
00051  *
00052  * SYNOPSIS:
00053  *
00054  * double x, y, y0();
00055  *
00056  * y = y0( x );
00057  *
00058  *
00059  *
00060  * DESCRIPTION:
00061  *
00062  * Returns Bessel function of the second kind, of order
00063  * zero, of the argument.
00064  *
00065  * The domain is divided into the intervals [0, 5] and
00066  * (5, infinity). In the first interval a rational approximation
00067  * R(x) is employed to compute
00068  *   y0(x)  = R(x)  +   2 * log(x) * j0(x) / PI.
00069  * Thus a call to j0() is required.
00070  *
00071  * In the second interval, the Hankel asymptotic expansion
00072  * is employed with two rational functions of degree 6/6
00073  * and 7/7.
00074  *
00075  *
00076  *
00077  * ACCURACY:
00078  *
00079  *  Absolute error, when y0(x) < 1; else relative error:
00080  *
00081  * arithmetic   domain     # trials      peak         rms
00082  *    DEC       0, 30        9400       7.0e-17     7.9e-18
00083  *    IEEE      0, 30       30000       1.3e-15     1.6e-16
00084  *
00085  */
00086 /*
00087 Cephes Math Library Release 2.1:  January, 1989
00088 Copyright 1984, 1987, 1989 by Stephen L. Moshier
00089 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
00090 */
00091 
00092 /* Note: all coefficients satisfy the relative error criterion
00093  * except YP, YQ which are designed for absolute error. */
00094 
00095 //#include "mconf.h"
00096 #define UNK
00097 #ifdef UNK
00098 static double PP[7] = {
00099   7.96936729297347051624E-4,
00100   8.28352392107440799803E-2,
00101   1.23953371646414299388E0,
00102   5.44725003058768775090E0,
00103   8.74716500199817011941E0,
00104   5.30324038235394892183E0,
00105   9.99999999999999997821E-1,
00106 };
00107 static double PQ[7] = {
00108   9.24408810558863637013E-4,
00109   8.56288474354474431428E-2,
00110   1.25352743901058953537E0,
00111   5.47097740330417105182E0,
00112   8.76190883237069594232E0,
00113   5.30605288235394617618E0,
00114   1.00000000000000000218E0,
00115 };
00116 #endif
00117 #ifdef DEC
00118 static unsigned short PP[28] = {
00119 0035520,0164604,0140733,0054470,
00120 0037251,0122605,0115356,0107170,
00121 0040236,0124412,0071500,0056303,
00122 0040656,0047737,0045720,0045263,
00123 0041013,0172143,0045004,0142103,
00124 0040651,0132045,0026241,0026406,
00125 0040200,0000000,0000000,0000000,
00126 };
00127 static unsigned short PQ[28] = {
00128 0035562,0052006,0070034,0134666,
00129 0037257,0057055,0055242,0123424,
00130 0040240,0071626,0046630,0032371,
00131 0040657,0011077,0032013,0012731,
00132 0041014,0030307,0050331,0006414,
00133 0040651,0145457,0065021,0150304,
00134 0040200,0000000,0000000,0000000,
00135 };
00136 #endif
00137 #ifdef IBMPC
00138 static unsigned short PP[28] = {
00139 0x6b27,0x983b,0x1d30,0x3f4a,
00140 0xd1cf,0xb35d,0x34b0,0x3fb5,
00141 0x0b98,0x4e68,0xd521,0x3ff3,
00142 0x0956,0xe97a,0xc9fb,0x4015,
00143 0x9888,0x6940,0x7e8c,0x4021,
00144 0x25a1,0xa594,0x3684,0x4015,
00145 0x0000,0x0000,0x0000,0x3ff0,
00146 };
00147 static unsigned short PQ[28] = {
00148 0x9737,0xce03,0x4a80,0x3f4e,
00149 0x54e3,0xab54,0xebc5,0x3fb5,
00150 0x069f,0xc9b3,0x0e72,0x3ff4,
00151 0x62bb,0xe681,0xe247,0x4015,
00152 0x21a1,0xea1b,0x8618,0x4021,
00153 0x3a19,0xed42,0x3965,0x4015,
00154 0x0000,0x0000,0x0000,0x3ff0,
00155 };
00156 #endif
00157 #ifdef MIEEE
00158 static unsigned short PP[28] = {
00159 0x3f4a,0x1d30,0x983b,0x6b27,
00160 0x3fb5,0x34b0,0xb35d,0xd1cf,
00161 0x3ff3,0xd521,0x4e68,0x0b98,
00162 0x4015,0xc9fb,0xe97a,0x0956,
00163 0x4021,0x7e8c,0x6940,0x9888,
00164 0x4015,0x3684,0xa594,0x25a1,
00165 0x3ff0,0x0000,0x0000,0x0000,
00166 };
00167 static unsigned short PQ[28] = {
00168 0x3f4e,0x4a80,0xce03,0x9737,
00169 0x3fb5,0xebc5,0xab54,0x54e3,
00170 0x3ff4,0x0e72,0xc9b3,0x069f,
00171 0x4015,0xe247,0xe681,0x62bb,
00172 0x4021,0x8618,0xea1b,0x21a1,
00173 0x4015,0x3965,0xed42,0x3a19,
00174 0x3ff0,0x0000,0x0000,0x0000,
00175 };
00176 #endif
00177 
00178 #ifdef UNK
00179 static double QP[8] = {
00180 -1.13663838898469149931E-2,
00181 -1.28252718670509318512E0,
00182 -1.95539544257735972385E1,
00183 -9.32060152123768231369E1,
00184 -1.77681167980488050595E2,
00185 -1.47077505154951170175E2,
00186 -5.14105326766599330220E1,
00187 -6.05014350600728481186E0,
00188 };
00189 static double QQ[7] = {
00190 /*  1.00000000000000000000E0,*/
00191   6.43178256118178023184E1,
00192   8.56430025976980587198E2,
00193   3.88240183605401609683E3,
00194   7.24046774195652478189E3,
00195   5.93072701187316984827E3,
00196   2.06209331660327847417E3,
00197   2.42005740240291393179E2,
00198 };
00199 #endif
00200 #ifdef DEC
00201 static unsigned short QP[32] = {
00202 0136472,0035021,0142451,0141115,
00203 0140244,0024731,0150620,0105642,
00204 0141234,0067177,0124161,0060141,
00205 0141672,0064572,0151557,0043036,
00206 0142061,0127141,0003127,0043517,
00207 0142023,0011727,0060271,0144544,
00208 0141515,0122142,0126620,0143150,
00209 0140701,0115306,0106715,0007344,
00210 };
00211 static unsigned short QQ[28] = {
00212 /*0040200,0000000,0000000,0000000,*/
00213 0041600,0121272,0004741,0026544,
00214 0042526,0015605,0105654,0161771,
00215 0043162,0123155,0165644,0062645,
00216 0043342,0041675,0167576,0130756,
00217 0043271,0052720,0165631,0154214,
00218 0043000,0160576,0034614,0172024,
00219 0042162,0000570,0030500,0051235,
00220 };
00221 #endif
00222 #ifdef IBMPC
00223 static unsigned short QP[32] = {
00224 0x384a,0x38a5,0x4742,0xbf87,
00225 0x1174,0x3a32,0x853b,0xbff4,
00226 0x2c0c,0xf50e,0x8dcf,0xc033,
00227 0xe8c4,0x5a6d,0x4d2f,0xc057,
00228 0xe8ea,0x20ca,0x35cc,0xc066,
00229 0x392d,0xec17,0x627a,0xc062,
00230 0x18cd,0x55b2,0xb48c,0xc049,
00231 0xa1dd,0xd1b9,0x3358,0xc018,
00232 };
00233 static unsigned short QQ[28] = {
00234 /*0x0000,0x0000,0x0000,0x3ff0,*/
00235 0x25ac,0x413c,0x1457,0x4050,
00236 0x9c7f,0xb175,0xc370,0x408a,
00237 0x8cb5,0xbd74,0x54cd,0x40ae,
00238 0xd63e,0xbdef,0x4877,0x40bc,
00239 0x3b11,0x1d73,0x2aba,0x40b7,
00240 0x9e82,0xc731,0x1c2f,0x40a0,
00241 0x0a54,0x0628,0x402f,0x406e,
00242 };
00243 #endif
00244 #ifdef MIEEE
00245 static unsigned short QP[32] = {
00246 0xbf87,0x4742,0x38a5,0x384a,
00247 0xbff4,0x853b,0x3a32,0x1174,
00248 0xc033,0x8dcf,0xf50e,0x2c0c,
00249 0xc057,0x4d2f,0x5a6d,0xe8c4,
00250 0xc066,0x35cc,0x20ca,0xe8ea,
00251 0xc062,0x627a,0xec17,0x392d,
00252 0xc049,0xb48c,0x55b2,0x18cd,
00253 0xc018,0x3358,0xd1b9,0xa1dd,
00254 };
00255 static unsigned short QQ[28] = {
00256 /*0x3ff0,0x0000,0x0000,0x0000,*/
00257 0x4050,0x1457,0x413c,0x25ac,
00258 0x408a,0xc370,0xb175,0x9c7f,
00259 0x40ae,0x54cd,0xbd74,0x8cb5,
00260 0x40bc,0x4877,0xbdef,0xd63e,
00261 0x40b7,0x2aba,0x1d73,0x3b11,
00262 0x40a0,0x1c2f,0xc731,0x9e82,
00263 0x406e,0x402f,0x0628,0x0a54,
00264 };
00265 #endif
00266 
00267 
00268 #ifdef UNK
00269 static double YP[8] = {
00270  1.55924367855235737965E4,
00271 -1.46639295903971606143E7,
00272  5.43526477051876500413E9,
00273 -9.82136065717911466409E11,
00274  8.75906394395366999549E13,
00275 -3.46628303384729719441E15,
00276  4.42733268572569800351E16,
00277 -1.84950800436986690637E16,
00278 };
00279 static double YQ[7] = {
00280 /* 1.00000000000000000000E0,*/
00281  1.04128353664259848412E3,
00282  6.26107330137134956842E5,
00283  2.68919633393814121987E8,
00284  8.64002487103935000337E10,
00285  2.02979612750105546709E13,
00286  3.17157752842975028269E15,
00287  2.50596256172653059228E17,
00288 };
00289 #endif
00290 #ifdef DEC
00291 static unsigned short YP[32] = {
00292 0043563,0120677,0042264,0046166,
00293 0146137,0140371,0113444,0042260,
00294 0050241,0175707,0100502,0063344,
00295 0152144,0125737,0007265,0164526,
00296 0053637,0051621,0163035,0060546,
00297 0155105,0004416,0107306,0060023,
00298 0056035,0045133,0030132,0000024,
00299 0155603,0065132,0144061,0131732,
00300 };
00301 static unsigned short YQ[28] = {
00302 /*0040200,0000000,0000000,0000000,*/
00303 0042602,0024422,0135557,0162663,
00304 0045030,0155665,0044075,0160135,
00305 0047200,0035432,0105446,0104005,
00306 0051240,0167331,0056063,0022743,
00307 0053223,0127746,0025764,0012160,
00308 0055064,0044206,0177532,0145545,
00309 0056536,0111375,0163715,0127201,
00310 };
00311 #endif
00312 #ifdef IBMPC
00313 static unsigned short YP[32] = {
00314 0x898f,0xe896,0x7437,0x40ce,
00315 0x8896,0x32e4,0xf81f,0xc16b,
00316 0x4cdd,0xf028,0x3f78,0x41f4,
00317 0xbd2b,0xe1d6,0x957b,0xc26c,
00318 0xac2d,0x3cc3,0xea72,0x42d3,
00319 0xcc02,0xd1d8,0xa121,0xc328,
00320 0x4003,0x660b,0xa94b,0x4363,
00321 0x367b,0x5906,0x6d4b,0xc350,
00322 };
00323 static unsigned short YQ[28] = {
00324 /*0x0000,0x0000,0x0000,0x3ff0,*/
00325 0xfcb6,0x576d,0x4522,0x4090,
00326 0xbc0c,0xa907,0x1b76,0x4123,
00327 0xd101,0x5164,0x0763,0x41b0,
00328 0x64bc,0x2b86,0x1ddb,0x4234,
00329 0x828e,0xc57e,0x75fc,0x42b2,
00330 0x596d,0xdfeb,0x8910,0x4326,
00331 0xb5d0,0xbcf9,0xd25f,0x438b,
00332 };
00333 #endif
00334 #ifdef MIEEE
00335 static unsigned short YP[32] = {
00336 0x40ce,0x7437,0xe896,0x898f,
00337 0xc16b,0xf81f,0x32e4,0x8896,
00338 0x41f4,0x3f78,0xf028,0x4cdd,
00339 0xc26c,0x957b,0xe1d6,0xbd2b,
00340 0x42d3,0xea72,0x3cc3,0xac2d,
00341 0xc328,0xa121,0xd1d8,0xcc02,
00342 0x4363,0xa94b,0x660b,0x4003,
00343 0xc350,0x6d4b,0x5906,0x367b,
00344 };
00345 static unsigned short YQ[28] = {
00346 /*0x3ff0,0x0000,0x0000,0x0000,*/
00347 0x4090,0x4522,0x576d,0xfcb6,
00348 0x4123,0x1b76,0xa907,0xbc0c,
00349 0x41b0,0x0763,0x5164,0xd101,
00350 0x4234,0x1ddb,0x2b86,0x64bc,
00351 0x42b2,0x75fc,0xc57e,0x828e,
00352 0x4326,0x8910,0xdfeb,0x596d,
00353 0x438b,0xd25f,0xbcf9,0xb5d0,
00354 };
00355 #endif
00356 
00357 #ifdef UNK
00358 /*  5.783185962946784521175995758455807035071 */
00359 static double DR1 = 5.78318596294678452118E0;
00360 /* 30.47126234366208639907816317502275584842 */
00361 static double DR2 = 3.04712623436620863991E1;
00362 #endif
00363 
00364 #ifdef DEC
00365 static unsigned short R1[] = {0040671,0007734,0001061,0056734};
00366 #define DR1 *(double *)R1
00367 static unsigned short R2[] = {0041363,0142445,0030416,0165567};
00368 #define DR2 *(double *)R2
00369 #endif
00370 
00371 #ifdef IBMPC
00372 static unsigned short R1[] = {0x2bbb,0x8046,0x21fb,0x4017};
00373 #define DR1 *(double *)R1
00374 static unsigned short R2[] = {0xdd6f,0xa621,0x78a4,0x403e};
00375 #define DR2 *(double *)R2
00376 #endif
00377 
00378 #ifdef MIEEE
00379 static unsigned short R1[] = {0x4017,0x21fb,0x8046,0x2bbb};
00380 #define DR1 *(double *)R1
00381 static unsigned short R2[] = {0x403e,0x78a4,0xa621,0xdd6f};
00382 #define DR2 *(double *)R2
00383 #endif
00384 
00385 #ifdef UNK
00386 static double RP[4] = {
00387 -4.79443220978201773821E9,
00388  1.95617491946556577543E12,
00389 -2.49248344360967716204E14,
00390  9.70862251047306323952E15,
00391 };
00392 static double RQ[8] = {
00393 /* 1.00000000000000000000E0,*/
00394  4.99563147152651017219E2,
00395  1.73785401676374683123E5,
00396  4.84409658339962045305E7,
00397  1.11855537045356834862E10,
00398  2.11277520115489217587E12,
00399  3.10518229857422583814E14,
00400  3.18121955943204943306E16,
00401  1.71086294081043136091E18,
00402 };
00403 #endif
00404 #ifdef DEC
00405 static unsigned short RP[16] = {
00406 0150216,0161235,0064344,0014450,
00407 0052343,0135216,0035624,0144153,
00408 0154142,0130247,0003310,0003667,
00409 0055411,0173703,0047772,0176635,
00410 };
00411 static unsigned short RQ[32] = {
00412 /*0040200,0000000,0000000,0000000,*/
00413 0042371,0144025,0032265,0136137,
00414 0044451,0133131,0132420,0151466,
00415 0046470,0144641,0072540,0030636,
00416 0050446,0126600,0045042,0044243,
00417 0052365,0172633,0110301,0071063,
00418 0054215,0032424,0062272,0043513,
00419 0055742,0005013,0171731,0072335,
00420 0057275,0170646,0036663,0013134,
00421 };
00422 #endif
00423 #ifdef IBMPC
00424 static unsigned short RP[16] = {
00425 0x8325,0xad1c,0xdc53,0xc1f1,
00426 0x990d,0xc772,0x7751,0x427c,
00427 0x00f7,0xe0d9,0x5614,0xc2ec,
00428 0x5fb4,0x69ff,0x3ef8,0x4341,
00429 };
00430 static unsigned short RQ[32] = {
00431 /*0x0000,0x0000,0x0000,0x3ff0,*/
00432 0xb78c,0xa696,0x3902,0x407f,
00433 0x1a67,0x36a2,0x36cb,0x4105,
00434 0x0634,0x2eac,0x1934,0x4187,
00435 0x4914,0x0944,0xd5b0,0x4204,
00436 0x2e46,0x7218,0xbeb3,0x427e,
00437 0x48e9,0x8c97,0xa6a2,0x42f1,
00438 0x2e9c,0x7e7b,0x4141,0x435c,
00439 0x62cc,0xc7b6,0xbe34,0x43b7,
00440 };
00441 #endif
00442 #ifdef MIEEE
00443 static unsigned short RP[16] = {
00444 0xc1f1,0xdc53,0xad1c,0x8325,
00445 0x427c,0x7751,0xc772,0x990d,
00446 0xc2ec,0x5614,0xe0d9,0x00f7,
00447 0x4341,0x3ef8,0x69ff,0x5fb4,
00448 };
00449 static unsigned short RQ[32] = {
00450 /*0x3ff0,0x0000,0x0000,0x0000,*/
00451 0x407f,0x3902,0xa696,0xb78c,
00452 0x4105,0x36cb,0x36a2,0x1a67,
00453 0x4187,0x1934,0x2eac,0x0634,
00454 0x4204,0xd5b0,0x0944,0x4914,
00455 0x427e,0xbeb3,0x7218,0x2e46,
00456 0x42f1,0xa6a2,0x8c97,0x48e9,
00457 0x435c,0x4141,0x7e7b,0x2e9c,
00458 0x43b7,0xbe34,0xc7b6,0x62cc,
00459 };
00460 #endif
00461 
00462 #define ANSIPROT
00463 #ifndef ANSIPROT
00464 double bessj0(), polevl(), p1evl(), log(), sin(), cos(), sqrt();
00465 #endif
00466 
00467 double polevl(double x,double coef[],int N)
00468 {
00469 double ans;
00470 int i;
00471 double *p;
00472 
00473 p = coef;
00474 ans = *p++;
00475 i = N;
00476 
00477 do
00478         ans = ans * x  +  *p++;
00479 while( --i );
00480 
00481 return( ans );
00482 }
00483 
00484 /*                                                      p1evl() */
00485 /*                                          N
00486  * Evaluate polynomial when coefficient of x  is 1.0.
00487  * Otherwise same as polevl.
00488  */
00489 
00490 double p1evl(double x, double coef[],int N)
00491 {
00492 double ans;
00493 double *p;
00494 int i;
00495 
00496 p = coef;
00497 ans = x + *p++;
00498 i = N-1;
00499 
00500 do
00501         ans = ans * x  + *p++;
00502 while( --i );
00503 
00504 return( ans );
00505 }
00506 
00507 double TWOOPI =  6.36619772367581343075535E-1; /* 2/pi */
00508 double THPIO4 =  2.35619449019234492885;       /* 3*pi/4 */
00509 double SQ2OPI =  7.9788456080286535587989E-1;  /* sqrt( 2/pi ) */
00510 double PIO4   =  7.85398163397448309616E-1;    /* pi/4 */
00511 
00512 double bessj0(double x)
00513 {
00514 double w, z, p, q, xn;
00515 
00516 if( x < 0 )
00517         x = -x;
00518 
00519 if( x <= 5.0 )
00520         {
00521         z = x * x;
00522         if( x < 1.0e-5 )
00523                 return( 1.0 - z/4.0 );
00524 
00525         p = (z - DR1) * (z - DR2);
00526         p = p * polevl( z, RP, 3)/p1evl( z, RQ, 8 );
00527         return( p );
00528         }
00529 
00530 w = 5.0/x;
00531 q = 25.0/(x*x);
00532 p = polevl( q, PP, 6)/polevl( q, PQ, 6 );
00533 q = polevl( q, QP, 7)/p1evl( q, QQ, 7 );
00534 xn = x - PIO4;
00535 p = p * cos(xn) - w * q * sin(xn);
00536 return( p * SQ2OPI / sqrt(x) );
00537 }
00538 /*                                                      bessy0() 2      */
00539 /* Bessel function of second kind, order zero   */
00540 
00541 /* Rational approximation coefficients YP[], YQ[] are used here.
00542  * The function computed is  bessy0(x)  -  2 * log(x) * bessj0(x) / PI,
00543  * whose value at x = 0 is  2 * ( log(0.5) + EUL ) / PI
00544  * = 0.073804295108687225.
00545  */
00546 
00547 
00548 double bessy0(double x)
00549 {
00550 double w, z, p, q, xn;
00551 
00552 if( x <= 5.0 )
00553         {
00554         if( x <= 0.0 )
00555           throw Exception ("arg<=0 in bessy0");
00556         z = x * x;
00557         w = polevl( z, YP, 7) / p1evl( z, YQ, 7 );
00558         w += TWOOPI * log(x) * bessj0(x);
00559         return( w );
00560         }
00561 
00562 w = 5.0/x;
00563 z = 25.0 / (x * x);
00564 p = polevl( z, PP, 6)/polevl( z, PQ, 6 );
00565 q = polevl( z, QP, 7)/p1evl( z, QQ, 7 );
00566 xn = x - PIO4;
00567 p = p * sin(xn) + w * q * cos(xn);
00568 return( p * SQ2OPI / sqrt(x) );
00569 }
00570 
00571 /*                                                      j1.c
00572  *
00573  *      Bessel function of order one
00574  *
00575  *
00576  *
00577  * SYNOPSIS:
00578  *
00579  * double x, y, j1();
00580  *
00581  * y = j1( x );
00582  *
00583  *
00584  *
00585  * DESCRIPTION:
00586  *
00587  * Returns Bessel function of order one of the argument.
00588  *
00589  * The domain is divided into the intervals [0, 8] and
00590  * (8, infinity). In the first interval a 24 term Chebyshev
00591  * expansion is used. In the second, the asymptotic
00592  * trigonometric representation is employed using two
00593  * rational functions of degree 5/5.
00594  *
00595  *
00596  *
00597  * ACCURACY:
00598  *
00599  *                      Absolute error:
00600  * arithmetic   domain      # trials      peak         rms
00601  *    DEC       0, 30       10000       4.0e-17     1.1e-17
00602  *    IEEE      0, 30       30000       2.6e-16     1.1e-16
00603  *
00604  *
00605  */
00606 /*                                                      y1.c
00607  *
00608  *      Bessel function of second kind of order one
00609  *
00610  *
00611  *
00612  * SYNOPSIS:
00613  *
00614  * double x, y, y1();
00615  *
00616  * y = y1( x );
00617  *
00618  *
00619  *
00620  * DESCRIPTION:
00621  *
00622  * Returns Bessel function of the second kind of order one
00623  * of the argument.
00624  *
00625  * The domain is divided into the intervals [0, 8] and
00626  * (8, infinity). In the first interval a 25 term Chebyshev
00627  * expansion is used, and a call to j1() is required.
00628  * In the second, the asymptotic trigonometric representation
00629  * is employed using two rational functions of degree 5/5.
00630  *
00631  *
00632  *
00633  * ACCURACY:
00634  *
00635  *                      Absolute error:
00636  * arithmetic   domain      # trials      peak         rms
00637  *    DEC       0, 30       10000       8.6e-17     1.3e-17
00638  *    IEEE      0, 30       30000       1.0e-15     1.3e-16
00639  *
00640  * (error criterion relative when |y1| > 1).
00641  *
00642  */
00643 /*
00644 Cephes Math Library Release 2.1:  January, 1989
00645 Copyright 1984, 1987, 1989 by Stephen L. Moshier
00646 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
00647 */
00648 
00649 #define PIO4 .78539816339744830962
00650 #define THPIO4 2.35619449019234492885
00651 #define SQ2OPI .79788456080286535588
00652 
00653 // #include "mconf.h"
00654 
00655 #ifdef UNK
00656 static double RP1[4] = {
00657 -8.99971225705559398224E8,
00658  4.52228297998194034323E11,
00659 -7.27494245221818276015E13,
00660  3.68295732863852883286E15,
00661 };
00662 static double RQ1[8] = {
00663 /* 1.00000000000000000000E0,*/
00664  6.20836478118054335476E2,
00665  2.56987256757748830383E5,
00666  8.35146791431949253037E7,
00667  2.21511595479792499675E10,
00668  4.74914122079991414898E12,
00669  7.84369607876235854894E14,
00670  8.95222336184627338078E16,
00671  5.32278620332680085395E18,
00672 };
00673 #endif
00674 #ifdef DEC
00675 static unsigned short RP1[16] = {
00676 0147526,0110742,0063322,0077052,
00677 0051722,0112720,0065034,0061530,
00678 0153604,0052227,0033147,0105650,
00679 0055121,0055025,0032276,0022015,
00680 };
00681 static unsigned short RQ1[32] = {
00682 /*0040200,0000000,0000000,0000000,*/
00683 0042433,0032610,0155604,0033473,
00684 0044572,0173320,0067270,0006616,
00685 0046637,0045246,0162225,0006606,
00686 0050645,0004773,0157577,0053004,
00687 0052612,0033734,0001667,0176501,
00688 0054462,0054121,0173147,0121367,
00689 0056237,0002777,0121451,0176007,
00690 0057623,0136253,0131601,0044710,
00691 };
00692 #endif
00693 #ifdef IBMPC
00694 static unsigned short RP1[16] = {
00695 0x4fc5,0x4cda,0xd23c,0xc1ca,
00696 0x8c6b,0x0d43,0x52ba,0x425a,
00697 0xf175,0xe6cc,0x8a92,0xc2d0,
00698 0xc482,0xa697,0x2b42,0x432a,
00699 };
00700 static unsigned short RQ1[32] = {
00701 /*0x0000,0x0000,0x0000,0x3ff0,*/
00702 0x86e7,0x1b70,0x66b1,0x4083,
00703 0x01b2,0x0dd7,0x5eda,0x410f,
00704 0xa1b1,0xdc92,0xe954,0x4193,
00705 0xeac1,0x7bef,0xa13f,0x4214,
00706 0xffa8,0x8076,0x46fb,0x4291,
00707 0xf45f,0x3ecc,0x4b0a,0x4306,
00708 0x3f81,0xf465,0xe0bf,0x4373,
00709 0x2939,0x7670,0x7795,0x43d2,
00710 };
00711 #endif
00712 #ifdef MIEEE
00713 static unsigned short RP1[16] = {
00714 0xc1ca,0xd23c,0x4cda,0x4fc5,
00715 0x425a,0x52ba,0x0d43,0x8c6b,
00716 0xc2d0,0x8a92,0xe6cc,0xf175,
00717 0x432a,0x2b42,0xa697,0xc482,
00718 };
00719 static unsigned short RQ1[32] = {
00720 /*0x3ff0,0x0000,0x0000,0x0000,*/
00721 0x4083,0x66b1,0x1b70,0x86e7,
00722 0x410f,0x5eda,0x0dd7,0x01b2,
00723 0x4193,0xe954,0xdc92,0xa1b1,
00724 0x4214,0xa13f,0x7bef,0xeac1,
00725 0x4291,0x46fb,0x8076,0xffa8,
00726 0x4306,0x4b0a,0x3ecc,0xf45f,
00727 0x4373,0xe0bf,0xf465,0x3f81,
00728 0x43d2,0x7795,0x7670,0x2939,
00729 };
00730 #endif
00731 
00732 #ifdef UNK
00733 static double PP1[7] = {
00734  7.62125616208173112003E-4,
00735  7.31397056940917570436E-2,
00736  1.12719608129684925192E0,
00737  5.11207951146807644818E0,
00738  8.42404590141772420927E0,
00739  5.21451598682361504063E0,
00740  1.00000000000000000254E0,
00741 };
00742 static double PQ1[7] = {
00743  5.71323128072548699714E-4,
00744  6.88455908754495404082E-2,
00745  1.10514232634061696926E0,
00746  5.07386386128601488557E0,
00747  8.39985554327604159757E0,
00748  5.20982848682361821619E0,
00749  9.99999999999999997461E-1,
00750 };
00751 #endif
00752 #ifdef DEC
00753 static unsigned short PP1[28] = {
00754 0035507,0144542,0061543,0024326,
00755 0037225,0145105,0017766,0022661,
00756 0040220,0043766,0010254,0133255,
00757 0040643,0113047,0142611,0151521,
00758 0041006,0144344,0055351,0074261,
00759 0040646,0156520,0120574,0006416,
00760 0040200,0000000,0000000,0000000,
00761 };
00762 static unsigned short PQ1[28] = {
00763 0035425,0142330,0115041,0165514,
00764 0037214,0177352,0145105,0052026,
00765 0040215,0072515,0141207,0073255,
00766 0040642,0056427,0137222,0106405,
00767 0041006,0062716,0166427,0165450,
00768 0040646,0133352,0035425,0123304,
00769 0040200,0000000,0000000,0000000,
00770 };
00771 #endif
00772 #ifdef IBMPC
00773 static unsigned short PP1[28] = {
00774 0x651b,0x4c6c,0xf92c,0x3f48,
00775 0xc4b6,0xa3fe,0xb948,0x3fb2,
00776 0x96d6,0xc215,0x08fe,0x3ff2,
00777 0x3a6a,0xf8b1,0x72c4,0x4014,
00778 0x2f16,0x8b5d,0xd91c,0x4020,
00779 0x81a2,0x142f,0xdbaa,0x4014,
00780 0x0000,0x0000,0x0000,0x3ff0,
00781 };
00782 static unsigned short PQ1[28] = {
00783 0x3d69,0x1344,0xb89b,0x3f42,
00784 0xaa83,0x5948,0x9fdd,0x3fb1,
00785 0xeed6,0xb850,0xaea9,0x3ff1,
00786 0x51a1,0xf7d2,0x4ba2,0x4014,
00787 0xfd65,0xdda2,0xccb9,0x4020,
00788 0xb4d9,0x4762,0xd6dd,0x4014,
00789 0x0000,0x0000,0x0000,0x3ff0,
00790 };
00791 #endif
00792 #ifdef MIEEE
00793 static unsigned short PP1[28] = {
00794 0x3f48,0xf92c,0x4c6c,0x651b,
00795 0x3fb2,0xb948,0xa3fe,0xc4b6,
00796 0x3ff2,0x08fe,0xc215,0x96d6,
00797 0x4014,0x72c4,0xf8b1,0x3a6a,
00798 0x4020,0xd91c,0x8b5d,0x2f16,
00799 0x4014,0xdbaa,0x142f,0x81a2,
00800 0x3ff0,0x0000,0x0000,0x0000,
00801 };
00802 static unsigned short PQ1[28] = {
00803 0x3f42,0xb89b,0x1344,0x3d69,
00804 0x3fb1,0x9fdd,0x5948,0xaa83,
00805 0x3ff1,0xaea9,0xb850,0xeed6,
00806 0x4014,0x4ba2,0xf7d2,0x51a1,
00807 0x4020,0xccb9,0xdda2,0xfd65,
00808 0x4014,0xd6dd,0x4762,0xb4d9,
00809 0x3ff0,0x0000,0x0000,0x0000,
00810 };
00811 #endif
00812 
00813 #ifdef UNK
00814 static double QP1[8] = {
00815  5.10862594750176621635E-2,
00816  4.98213872951233449420E0,
00817  7.58238284132545283818E1,
00818  3.66779609360150777800E2,
00819  7.10856304998926107277E2,
00820  5.97489612400613639965E2,
00821  2.11688757100572135698E2,
00822  2.52070205858023719784E1,
00823 };
00824 static double QQ1[7] = {
00825 /* 1.00000000000000000000E0,*/
00826  7.42373277035675149943E1,
00827  1.05644886038262816351E3,
00828  4.98641058337653607651E3,
00829  9.56231892404756170795E3,
00830  7.99704160447350683650E3,
00831  2.82619278517639096600E3,
00832  3.36093607810698293419E2,
00833 };
00834 #endif
00835 #ifdef DEC
00836 static unsigned short QP1[32] = {
00837 0037121,0037723,0055605,0151004,
00838 0040637,0066656,0031554,0077264,
00839 0041627,0122714,0153170,0161466,
00840 0042267,0061712,0036520,0140145,
00841 0042461,0133315,0131573,0071176,
00842 0042425,0057525,0147500,0013201,
00843 0042123,0130122,0061245,0154131,
00844 0041311,0123772,0064254,0172650,
00845 };
00846 static unsigned short QQ1[28] = {
00847 /*0040200,0000000,0000000,0000000,*/
00848 0041624,0074603,0002112,0101670,
00849 0042604,0007135,0010162,0175565,
00850 0043233,0151510,0157757,0172010,
00851 0043425,0064506,0112006,0104276,
00852 0043371,0164125,0032271,0164242,
00853 0043060,0121425,0122750,0136013,
00854 0042250,0005773,0053472,0146267,
00855 };
00856 #endif
00857 #ifdef IBMPC
00858 static unsigned short QP1[32] = {
00859 0xba40,0x6b70,0x27fa,0x3faa,
00860 0x8fd6,0xc66d,0xedb5,0x4013,
00861 0x1c67,0x9acf,0xf4b9,0x4052,
00862 0x180d,0x47aa,0xec79,0x4076,
00863 0x6e50,0xb66f,0x36d9,0x4086,
00864 0x02d0,0xb9e8,0xabea,0x4082,
00865 0xbb0b,0x4c54,0x760a,0x406a,
00866 0x9eb5,0x4d15,0x34ff,0x4039,
00867 };
00868 static unsigned short QQ1[28] = {
00869 /*0x0000,0x0000,0x0000,0x3ff0,*/
00870 0x5077,0x6089,0x8f30,0x4052,
00871 0x5f6f,0xa20e,0x81cb,0x4090,
00872 0xfe81,0x1bfd,0x7a69,0x40b3,
00873 0xd118,0xd280,0xad28,0x40c2,
00874 0x3d14,0xa697,0x3d0a,0x40bf,
00875 0x1781,0xb4bd,0x1462,0x40a6,
00876 0x5997,0x6ae7,0x017f,0x4075,
00877 };
00878 #endif
00879 #ifdef MIEEE
00880 static unsigned short QP1[32] = {
00881 0x3faa,0x27fa,0x6b70,0xba40,
00882 0x4013,0xedb5,0xc66d,0x8fd6,
00883 0x4052,0xf4b9,0x9acf,0x1c67,
00884 0x4076,0xec79,0x47aa,0x180d,
00885 0x4086,0x36d9,0xb66f,0x6e50,
00886 0x4082,0xabea,0xb9e8,0x02d0,
00887 0x406a,0x760a,0x4c54,0xbb0b,
00888 0x4039,0x34ff,0x4d15,0x9eb5,
00889 };
00890 static unsigned short QQ1[28] = {
00891 /*0x3ff0,0x0000,0x0000,0x0000,*/
00892 0x4052,0x8f30,0x6089,0x5077,
00893 0x4090,0x81cb,0xa20e,0x5f6f,
00894 0x40b3,0x7a69,0x1bfd,0xfe81,
00895 0x40c2,0xad28,0xd280,0xd118,
00896 0x40bf,0x3d0a,0xa697,0x3d14,
00897 0x40a6,0x1462,0xb4bd,0x1781,
00898 0x4075,0x017f,0x6ae7,0x5997,
00899 };
00900 #endif
00901 
00902 #ifdef UNK
00903 static double YP1[6] = {
00904  1.26320474790178026440E9,
00905 -6.47355876379160291031E11,
00906  1.14509511541823727583E14,
00907 -8.12770255501325109621E15,
00908  2.02439475713594898196E17,
00909 -7.78877196265950026825E17,
00910 };
00911 static double YQ1[8] = {
00912 /* 1.00000000000000000000E0,*/
00913  5.94301592346128195359E2,
00914  2.35564092943068577943E5,
00915  7.34811944459721705660E7,
00916  1.87601316108706159478E10,
00917  3.88231277496238566008E12,
00918  6.20557727146953693363E14,
00919  6.87141087355300489866E16,
00920  3.97270608116560655612E18,
00921 };
00922 #endif
00923 #ifdef DEC
00924 static unsigned short YP1[24] = {
00925 0047626,0112763,0013715,0133045,
00926 0152026,0134552,0142033,0024411,
00927 0053720,0045245,0102210,0077565,
00928 0155347,0000321,0136415,0102031,
00929 0056463,0146550,0055633,0032605,
00930 0157054,0171012,0167361,0054265,
00931 };
00932 static unsigned short YQ1[32] = {
00933 /*0040200,0000000,0000000,0000000,*/
00934 0042424,0111515,0044773,0153014,
00935 0044546,0005405,0171307,0075774,
00936 0046614,0023575,0047105,0063556,
00937 0050613,0143034,0101533,0156026,
00938 0052541,0175367,0166514,0114257,
00939 0054415,0014466,0134350,0171154,
00940 0056164,0017436,0025075,0022101,
00941 0057534,0103614,0103663,0121772,
00942 };
00943 #endif
00944 #ifdef IBMPC
00945 static unsigned short YP1[24] = {
00946 0xb6c5,0x62f9,0xd2be,0x41d2,
00947 0x6521,0x5883,0xd72d,0xc262,
00948 0x0fef,0xb091,0x0954,0x42da,
00949 0xb083,0x37a1,0xe01a,0xc33c,
00950 0x66b1,0x0b73,0x79ad,0x4386,
00951 0x2b17,0x5dde,0x9e41,0xc3a5,
00952 };
00953 static unsigned short YQ1[32] = {
00954 /*0x0000,0x0000,0x0000,0x3ff0,*/
00955 0x7ac2,0xa93f,0x9269,0x4082,
00956 0xef7f,0xbe58,0xc160,0x410c,
00957 0xacee,0xa9c8,0x84ef,0x4191,
00958 0x7b83,0x906b,0x78c3,0x4211,
00959 0x9316,0xfda9,0x3f5e,0x428c,
00960 0x1e4e,0xd71d,0xa326,0x4301,
00961 0xa488,0xc547,0x83e3,0x436e,
00962 0x747f,0x90f6,0x90f1,0x43cb,
00963 };
00964 #endif
00965 #ifdef MIEEE
00966 static unsigned short YP1[24] = {
00967 0x41d2,0xd2be,0x62f9,0xb6c5,
00968 0xc262,0xd72d,0x5883,0x6521,
00969 0x42da,0x0954,0xb091,0x0fef,
00970 0xc33c,0xe01a,0x37a1,0xb083,
00971 0x4386,0x79ad,0x0b73,0x66b1,
00972 0xc3a5,0x9e41,0x5dde,0x2b17,
00973 };
00974 static unsigned short YQ1[32] = {
00975 /*0x3ff0,0x0000,0x0000,0x0000,*/
00976 0x4082,0x9269,0xa93f,0x7ac2,
00977 0x410c,0xc160,0xbe58,0xef7f,
00978 0x4191,0x84ef,0xa9c8,0xacee,
00979 0x4211,0x78c3,0x906b,0x7b83,
00980 0x428c,0x3f5e,0xfda9,0x9316,
00981 0x4301,0xa326,0xd71d,0x1e4e,
00982 0x436e,0x83e3,0xc547,0xa488,
00983 0x43cb,0x90f1,0x90f6,0x747f,
00984 };
00985 #endif
00986 
00987 
00988 #ifdef UNK
00989 static double Z1 = 1.46819706421238932572E1;
00990 static double Z2 = 4.92184563216946036703E1;
00991 #endif
00992 
00993 #ifdef DEC
00994 static unsigned short DZ1[] = {0041152,0164532,0006114,0010540};
00995 static unsigned short DZ2[] = {0041504,0157663,0001625,0020621};
00996 #define Z1 (*(double *)DZ1)
00997 #define Z2 (*(double *)DZ2)
00998 #endif
00999 
01000 #ifdef IBMPC
01001 static unsigned short DZ1[] = {0x822c,0x4189,0x5d2b,0x402d};
01002 static unsigned short DZ2[] = {0xa432,0x6072,0x9bf6,0x4048};
01003 #define Z1 (*(double *)DZ1)
01004 #define Z2 (*(double *)DZ2)
01005 #endif
01006 
01007 #ifdef MIEEE
01008 static unsigned short DZ1[] = {0x402d,0x5d2b,0x4189,0x822c};
01009 static unsigned short DZ2[] = {0x4048,0x9bf6,0x6072,0xa432};
01010 #define Z1 (*(double *)DZ1)
01011 #define Z2 (*(double *)DZ2)
01012 #endif
01013 
01014 #ifndef ANSIPROT
01015 double bessj1(), polevl(), p1evl(), log(), sin(), cos(), sqrt();
01016 #endif
01017 
01018 double bessj1(double x)
01019 {
01020 double w, z, p, q, xn;
01021 
01022 w = x;
01023 if( x < 0 )
01024         w = -x;
01025 
01026 if( w <= 5.0 )
01027         {
01028         z = x * x;      
01029         w = polevl( z, RP1, 3 ) / p1evl( z, RQ1, 8 );
01030         w = w * x * (z - Z1) * (z - Z2);
01031         return( w );
01032         }
01033 
01034 w = 5.0/x;
01035 z = w * w;
01036 p = polevl( z, PP1, 6)/polevl( z, PQ1, 6 );
01037 q = polevl( z, QP1, 7)/p1evl( z, QQ1, 7 );
01038 xn = x - THPIO4;
01039 p = p * cos(xn) - w * q * sin(xn);
01040 return( p * SQ2OPI / sqrt(x) );
01041 }
01042 
01043 double bessy1(double x)
01044 {
01045 double w, z, p, q, xn;
01046 
01047 if( x <= 5.0 )
01048         {
01049         if( x <= 0.0 )
01050           throw Exception("arg<=0 in bessy1");
01051         z = x * x;
01052         w = x * (polevl( z, YP1, 5 ) / p1evl( z, YQ1, 8 ));
01053         w += TWOOPI * ( bessj1(x) * log(x)  -  1.0/x );
01054         return( w );
01055         }
01056 
01057 w = 5.0/x;
01058 z = w * w;
01059 p = polevl( z, PP1, 6)/polevl( z, PQ1, 6 );
01060 q = polevl( z, QP1, 7)/p1evl( z, QQ1, 7 );
01061 xn = x - THPIO4;
01062 p = p * sin(xn) + w * q * cos(xn);
01063 return( p * SQ2OPI / sqrt(x) );
01064 }