NGSolve
4.9
|
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 }