00001
00002
00003 #include "cddefines.h"
00004 #include "thirdparty.h"
00005 #include "physconst.h"
00006
00007 inline double polevl(double x, const double coef[], int N);
00008 inline double p1evl(double x, const double coef[], int N);
00009 inline double chbevl(double, const double[], int);
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045 bool linfit(
00046 long n,
00047 double x[],
00048 double y[],
00049 double &a,
00050 double &siga,
00051 double &b,
00052 double &sigb
00053 )
00054 {
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092 DEBUG_ENTRY( "linfit()" );
00093
00094
00095 a = 0.0;
00096 siga = 0.0;
00097 b = 0.0;
00098 sigb = 0.0;
00099
00100
00101 double s1 = 0.0;
00102 double s2 = 0.0;
00103 for( long i=0; i < n; i++ )
00104 {
00105 s1 += x[i];
00106 s2 += y[i];
00107 }
00108 double rn = (double)n;
00109 double xavg = s1/rn;
00110 double yavg = s2/rn;
00111 double sxx = 0.0;
00112 double sxy = 0.0;
00113 for( long i=0; i < n; i++ )
00114 {
00115 x[i] -= xavg;
00116 y[i] -= yavg;
00117 sxx += pow2(x[i]);
00118 sxy += x[i]*y[i];
00119 }
00120
00121 if( sxy == 0.0 )
00122 {
00123 DEBUG_EXIT( "linfit()" );
00124 return true;
00125 }
00126
00127
00128 b = sxy/sxx;
00129
00130
00131 a = yavg - b*xavg;
00132
00133
00134 double sum1 = 0.0;
00135 for( long i=0; i < n; i++ )
00136 sum1 += pow2(x[i]*(y[i] - b*x[i]));
00137
00138
00139 sigb = sum1/pow2(sxx);
00140
00141
00142 for( long i=0; i < n; i++ )
00143 siga += pow2((y[i] - b*x[i])*(1.0 - rn*xavg*x[i]/sxx));
00144
00145
00146 sigb = sqrt(sigb);
00147 siga = sqrt(siga)/rn;
00148
00149
00150 for( long i=0; i < n; i++ )
00151 {
00152 x[i] += xavg;
00153 y[i] += yavg;
00154 }
00155
00156 DEBUG_EXIT( "linfit()" );
00157 return false;
00158 }
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168 static const double pre_factorial[NPRE_FACTORIAL] =
00169 {
00170 1.00000000000000000000e+00,
00171 1.00000000000000000000e+00,
00172 2.00000000000000000000e+00,
00173 6.00000000000000000000e+00,
00174 2.40000000000000000000e+01,
00175 1.20000000000000000000e+02,
00176 7.20000000000000000000e+02,
00177 5.04000000000000000000e+03,
00178 4.03200000000000000000e+04,
00179 3.62880000000000000000e+05,
00180 3.62880000000000000000e+06,
00181 3.99168000000000000000e+07,
00182 4.79001600000000000000e+08,
00183 6.22702080000000000000e+09,
00184 8.71782912000000000000e+10,
00185 1.30767436800000000000e+12,
00186 2.09227898880000000000e+13,
00187 3.55687428096000000000e+14,
00188 6.40237370572800000000e+15,
00189 1.21645100408832000000e+17,
00190 2.43290200817664000000e+18,
00191 5.10909421717094400000e+19,
00192 1.12400072777760768000e+21,
00193 2.58520167388849766400e+22,
00194 6.20448401733239439360e+23,
00195 1.55112100433309859840e+25,
00196 4.03291461126605635592e+26,
00197 1.08888694504183521614e+28,
00198 3.04888344611713860511e+29,
00199 8.84176199373970195470e+30,
00200 2.65252859812191058647e+32,
00201 8.22283865417792281807e+33,
00202 2.63130836933693530178e+35,
00203 8.68331761881188649615e+36,
00204 2.95232799039604140861e+38,
00205 1.03331479663861449300e+40,
00206 3.71993326789901217463e+41,
00207 1.37637530912263450457e+43,
00208 5.23022617466601111726e+44,
00209 2.03978820811974433568e+46,
00210 8.15915283247897734264e+47,
00211 3.34525266131638071044e+49,
00212 1.40500611775287989839e+51,
00213 6.04152630633738356321e+52,
00214 2.65827157478844876773e+54,
00215 1.19622220865480194551e+56,
00216 5.50262215981208894940e+57,
00217 2.58623241511168180614e+59,
00218 1.24139155925360726691e+61,
00219 6.08281864034267560801e+62,
00220 3.04140932017133780398e+64,
00221 1.55111875328738227999e+66,
00222 8.06581751709438785585e+67,
00223 4.27488328406002556374e+69,
00224 2.30843697339241380441e+71,
00225 1.26964033536582759243e+73,
00226 7.10998587804863451749e+74,
00227 4.05269195048772167487e+76,
00228 2.35056133128287857145e+78,
00229 1.38683118545689835713e+80,
00230 8.32098711274139014271e+81,
00231 5.07580213877224798711e+83,
00232 3.14699732603879375200e+85,
00233 1.98260831540444006372e+87,
00234 1.26886932185884164078e+89,
00235 8.24765059208247066512e+90,
00236 5.44344939077443063905e+92,
00237 3.64711109181886852801e+94,
00238 2.48003554243683059915e+96,
00239 1.71122452428141311337e+98,
00240 1.19785716699698917933e+100,
00241 8.50478588567862317347e+101,
00242 6.12344583768860868500e+103,
00243 4.47011546151268434004e+105,
00244 3.30788544151938641157e+107,
00245 2.48091408113953980872e+109,
00246 1.88549470166605025458e+111,
00247 1.45183092028285869606e+113,
00248 1.13242811782062978295e+115,
00249 8.94618213078297528506e+116,
00250 7.15694570462638022794e+118,
00251 5.79712602074736798470e+120,
00252 4.75364333701284174746e+122,
00253 3.94552396972065865030e+124,
00254 3.31424013456535326627e+126,
00255 2.81710411438055027626e+128,
00256 2.42270953836727323750e+130,
00257 2.10775729837952771662e+132,
00258 1.85482642257398439069e+134,
00259 1.65079551609084610774e+136,
00260 1.48571596448176149700e+138,
00261 1.35200152767840296226e+140,
00262 1.24384140546413072522e+142,
00263 1.15677250708164157442e+144,
00264 1.08736615665674307994e+146,
00265 1.03299784882390592592e+148,
00266 9.91677934870949688836e+149,
00267 9.61927596824821198159e+151,
00268 9.42689044888324774164e+153,
00269 9.33262154439441526381e+155,
00270 9.33262154439441526388e+157,
00271 9.42594775983835941673e+159,
00272 9.61446671503512660515e+161,
00273 9.90290071648618040340e+163,
00274 1.02990167451456276198e+166,
00275 1.08139675824029090008e+168,
00276 1.14628056373470835406e+170,
00277 1.22652020319613793888e+172,
00278 1.32464181945182897396e+174,
00279 1.44385958320249358163e+176,
00280 1.58824554152274293982e+178,
00281 1.76295255109024466316e+180,
00282 1.97450685722107402277e+182,
00283 2.23119274865981364576e+184,
00284 2.54355973347218755612e+186,
00285 2.92509369349301568964e+188,
00286 3.39310868445189820004e+190,
00287 3.96993716080872089396e+192,
00288 4.68452584975429065488e+194,
00289 5.57458576120760587943e+196,
00290 6.68950291344912705515e+198,
00291 8.09429852527344373681e+200,
00292 9.87504420083360135884e+202,
00293 1.21463043670253296712e+205,
00294 1.50614174151114087918e+207,
00295 1.88267717688892609901e+209,
00296 2.37217324288004688470e+211,
00297 3.01266001845765954361e+213,
00298 3.85620482362580421582e+215,
00299 4.97450422247728743840e+217,
00300 6.46685548922047366972e+219,
00301 8.47158069087882050755e+221,
00302 1.11824865119600430699e+224,
00303 1.48727070609068572828e+226,
00304 1.99294274616151887582e+228,
00305 2.69047270731805048244e+230,
00306 3.65904288195254865604e+232,
00307 5.01288874827499165889e+234,
00308 6.91778647261948848943e+236,
00309 9.61572319694108900019e+238,
00310 1.34620124757175246000e+241,
00311 1.89814375907617096864e+243,
00312 2.69536413788816277557e+245,
00313 3.85437071718007276916e+247,
00314 5.55029383273930478744e+249,
00315 8.04792605747199194159e+251,
00316 1.17499720439091082343e+254,
00317 1.72724589045463891049e+256,
00318 2.55632391787286558753e+258,
00319 3.80892263763056972532e+260,
00320 5.71338395644585458806e+262,
00321 8.62720977423324042775e+264,
00322 1.31133588568345254503e+267,
00323 2.00634390509568239384e+269,
00324 3.08976961384735088657e+271,
00325 4.78914290146339387432e+273,
00326 7.47106292628289444390e+275,
00327 1.17295687942641442768e+278,
00328 1.85327186949373479574e+280,
00329 2.94670227249503832518e+282,
00330 4.71472363599206132029e+284,
00331 7.59070505394721872577e+286,
00332 1.22969421873944943358e+289,
00333 2.00440157654530257674e+291,
00334 3.28721858553429622598e+293,
00335 5.42391066613158877297e+295,
00336 9.00369170577843736335e+297,
00337 1.50361651486499903974e+300,
00338 2.52607574497319838672e+302,
00339 4.26906800900470527345e+304,
00340 7.25741561530799896496e+306
00341 };
00342
00343 double factorial( long n )
00344 {
00345 if( n < 0 || n >= NPRE_FACTORIAL )
00346 {
00347 fprintf( ioQQQ, "factorial: domain error\n" );
00348 puts( "[Stop in factorial]" );
00349 cdEXIT(EXIT_FAILURE);
00350 }
00351
00352 return pre_factorial[n];
00353 }
00354
00355
00356
00357 class t_lfact : public Singleton<t_lfact>
00358 {
00359 friend class Singleton<t_lfact>;
00360 protected:
00361 t_lfact()
00362 {
00363 lf.reserve( 512 );
00364 lf.push_back( 0. );
00365 lf.push_back( 0. );
00366 }
00367 private:
00368 vector<double> lf;
00369 public:
00370 double get_lfact( unsigned long n )
00371 {
00372 if( n < lf.size() )
00373 {
00374 return lf[n];
00375 }
00376 else
00377 {
00378 for( unsigned long i=(unsigned long)lf.size(); i <= n; i++ )
00379 lf.push_back( lf[i-1] + log10(static_cast<double>(i)) );
00380 return lf[n];
00381 }
00382 }
00383 };
00384
00385 double lfactorial( long n )
00386 {
00387
00388
00389
00390
00391
00392
00393
00394
00395 if( n < 0 )
00396 {
00397 fprintf( ioQQQ, "lfactorial: domain error\n" );
00398 puts( "[Stop in lfactorial]" );
00399 cdEXIT(EXIT_FAILURE);
00400 }
00401
00402 return t_lfact::Inst().get_lfact( static_cast<unsigned long>(n) );
00403 }
00404
00405
00406
00407
00408
00409
00410
00411
00412 complex<double> cdgamma(complex<double> x)
00413 {
00414 double xr, xi, wr, wi, ur, ui, vr, vi, yr, yi, t;
00415
00416 DEBUG_ENTRY( "cdgamma()" );
00417
00418 xr = x.real();
00419 xi = x.imag();
00420 if(xr < 0)
00421 {
00422 wr = 1. - xr;
00423 wi = -xi;
00424 }
00425 else
00426 {
00427 wr = xr;
00428 wi = xi;
00429 }
00430 ur = wr + 6.00009857740312429;
00431 vr = ur * (wr + 4.99999857982434025) - wi * wi;
00432 vi = wi * (wr + 4.99999857982434025) + ur * wi;
00433 yr = ur * 13.2280130755055088 + vr * 66.2756400966213521 +
00434 0.293729529320536228;
00435 yi = wi * 13.2280130755055088 + vi * 66.2756400966213521;
00436 ur = vr * (wr + 4.00000003016801681) - vi * wi;
00437 ui = vi * (wr + 4.00000003016801681) + vr * wi;
00438 vr = ur * (wr + 2.99999999944915534) - ui * wi;
00439 vi = ui * (wr + 2.99999999944915534) + ur * wi;
00440 yr += ur * 91.1395751189899762 + vr * 47.3821439163096063;
00441 yi += ui * 91.1395751189899762 + vi * 47.3821439163096063;
00442 ur = vr * (wr + 2.00000000000603851) - vi * wi;
00443 ui = vi * (wr + 2.00000000000603851) + vr * wi;
00444 vr = ur * (wr + 0.999999999999975753) - ui * wi;
00445 vi = ui * (wr + 0.999999999999975753) + ur * wi;
00446 yr += ur * 10.5400280458730808 + vr;
00447 yi += ui * 10.5400280458730808 + vi;
00448 ur = vr * wr - vi * wi;
00449 ui = vi * wr + vr * wi;
00450 t = ur * ur + ui * ui;
00451 vr = yr * ur + yi * ui + t * 0.0327673720261526849;
00452 vi = yi * ur - yr * ui;
00453 yr = wr + 7.31790632447016203;
00454 ur = log(yr * yr + wi * wi) * 0.5 - 1;
00455 ui = atan2(wi, yr);
00456 yr = exp(ur * (wr - 0.5) - ui * wi - 3.48064577727581257) / t;
00457 yi = ui * (wr - 0.5) + ur * wi;
00458 ur = yr * cos(yi);
00459 ui = yr * sin(yi);
00460 yr = ur * vr - ui * vi;
00461 yi = ui * vr + ur * vi;
00462 if(xr < 0)
00463 {
00464 wr = xr * 3.14159265358979324;
00465 wi = exp(xi * 3.14159265358979324);
00466 vi = 1 / wi;
00467 ur = (vi + wi) * sin(wr);
00468 ui = (vi - wi) * cos(wr);
00469 vr = ur * yr + ui * yi;
00470 vi = ui * yr - ur * yi;
00471 ur = 6.2831853071795862 / (vr * vr + vi * vi);
00472 yr = ur * vr;
00473 yi = ur * vi;
00474 }
00475
00476 DEBUG_EXIT( "cdgamma()" );
00477 return complex<double>( yr, yi );
00478 }
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594 static const double b0_PP[7] = {
00595 7.96936729297347051624e-4,
00596 8.28352392107440799803e-2,
00597 1.23953371646414299388e0,
00598 5.44725003058768775090e0,
00599 8.74716500199817011941e0,
00600 5.30324038235394892183e0,
00601 9.99999999999999997821e-1,
00602 };
00603
00604 static const double b0_PQ[7] = {
00605 9.24408810558863637013e-4,
00606 8.56288474354474431428e-2,
00607 1.25352743901058953537e0,
00608 5.47097740330417105182e0,
00609 8.76190883237069594232e0,
00610 5.30605288235394617618e0,
00611 1.00000000000000000218e0,
00612 };
00613
00614 static const double b0_QP[8] = {
00615 -1.13663838898469149931e-2,
00616 -1.28252718670509318512e0,
00617 -1.95539544257735972385e1,
00618 -9.32060152123768231369e1,
00619 -1.77681167980488050595e2,
00620 -1.47077505154951170175e2,
00621 -5.14105326766599330220e1,
00622 -6.05014350600728481186e0,
00623 };
00624
00625 static const double b0_QQ[7] = {
00626
00627 6.43178256118178023184e1,
00628 8.56430025976980587198e2,
00629 3.88240183605401609683e3,
00630 7.24046774195652478189e3,
00631 5.93072701187316984827e3,
00632 2.06209331660327847417e3,
00633 2.42005740240291393179e2,
00634 };
00635
00636 static const double b0_YP[8] = {
00637 1.55924367855235737965e4,
00638 -1.46639295903971606143e7,
00639 5.43526477051876500413e9,
00640 -9.82136065717911466409e11,
00641 8.75906394395366999549e13,
00642 -3.46628303384729719441e15,
00643 4.42733268572569800351e16,
00644 -1.84950800436986690637e16,
00645 };
00646
00647 static const double b0_YQ[7] = {
00648
00649 1.04128353664259848412e3,
00650 6.26107330137134956842e5,
00651 2.68919633393814121987e8,
00652 8.64002487103935000337e10,
00653 2.02979612750105546709e13,
00654 3.17157752842975028269e15,
00655 2.50596256172653059228e17,
00656 };
00657
00658
00659 static const double DR1 = 5.78318596294678452118e0;
00660
00661 static const double DR2 = 3.04712623436620863991e1;
00662
00663 static double b0_RP[4] = {
00664 -4.79443220978201773821e9,
00665 1.95617491946556577543e12,
00666 -2.49248344360967716204e14,
00667 9.70862251047306323952e15,
00668 };
00669
00670 static double b0_RQ[8] = {
00671
00672 4.99563147152651017219e2,
00673 1.73785401676374683123e5,
00674 4.84409658339962045305e7,
00675 1.11855537045356834862e10,
00676 2.11277520115489217587e12,
00677 3.10518229857422583814e14,
00678 3.18121955943204943306e16,
00679 1.71086294081043136091e18,
00680 };
00681
00682 static const double TWOOPI = 2./PI;
00683 static const double SQ2OPI = sqrt(2./PI);
00684 static const double PIO4 = PI/4.;
00685
00686 double bessel_j0(double x)
00687 {
00688 double w, z, p, q, xn;
00689
00690 DEBUG_ENTRY( "bessel_j0()" );
00691
00692 if( x < 0 )
00693 x = -x;
00694
00695 if( x <= 5.0 )
00696 {
00697 z = x * x;
00698 if( x < 1.0e-5 )
00699 return 1.0 - z/4.0;
00700
00701 p = (z - DR1) * (z - DR2);
00702 p = p * polevl( z, b0_RP, 3)/p1evl( z, b0_RQ, 8 );
00703
00704 DEBUG_EXIT( "bessel_j0()" );
00705 return p;
00706 }
00707
00708 w = 5.0/x;
00709 q = 25.0/(x*x);
00710 p = polevl( q, b0_PP, 6)/polevl( q, b0_PQ, 6 );
00711 q = polevl( q, b0_QP, 7)/p1evl( q, b0_QQ, 7 );
00712 xn = x - PIO4;
00713 p = p * cos(xn) - w * q * sin(xn);
00714
00715 DEBUG_EXIT( "bessel_j0()" );
00716 return p * SQ2OPI / sqrt(x);
00717 }
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728 double bessel_y0(double x)
00729 {
00730 double w, z, p, q, xn;
00731
00732 DEBUG_ENTRY( "bessel_y0()" );
00733
00734 if( x <= 5.0 )
00735 {
00736 if( x <= 0.0 )
00737 {
00738 fprintf( ioQQQ, "bessel_y0: domain error\n" );
00739 puts( "[Stop in bessel_y0]" );
00740 cdEXIT(EXIT_FAILURE);
00741 }
00742 z = x * x;
00743 w = polevl( z, b0_YP, 7 ) / p1evl( z, b0_YQ, 7 );
00744 w += TWOOPI * log(x) * bessel_j0(x);
00745
00746 DEBUG_EXIT( "bessel_y0()" );
00747 return w;
00748 }
00749
00750 w = 5.0/x;
00751 z = 25.0 / (x * x);
00752 p = polevl( z, b0_PP, 6)/polevl( z, b0_PQ, 6 );
00753 q = polevl( z, b0_QP, 7)/p1evl( z, b0_QQ, 7 );
00754 xn = x - PIO4;
00755 p = p * sin(xn) + w * q * cos(xn);
00756
00757 DEBUG_EXIT( "bessel_y0()" );
00758 return p * SQ2OPI / sqrt(x);
00759 }
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839 static const double b1_RP[4] = {
00840 -8.99971225705559398224e8,
00841 4.52228297998194034323e11,
00842 -7.27494245221818276015e13,
00843 3.68295732863852883286e15,
00844 };
00845
00846 static const double b1_RQ[8] = {
00847
00848 6.20836478118054335476e2,
00849 2.56987256757748830383e5,
00850 8.35146791431949253037e7,
00851 2.21511595479792499675e10,
00852 4.74914122079991414898e12,
00853 7.84369607876235854894e14,
00854 8.95222336184627338078e16,
00855 5.32278620332680085395e18,
00856 };
00857
00858 static const double b1_PP[7] = {
00859 7.62125616208173112003e-4,
00860 7.31397056940917570436e-2,
00861 1.12719608129684925192e0,
00862 5.11207951146807644818e0,
00863 8.42404590141772420927e0,
00864 5.21451598682361504063e0,
00865 1.00000000000000000254e0,
00866 };
00867
00868 static const double b1_PQ[7] = {
00869 5.71323128072548699714e-4,
00870 6.88455908754495404082e-2,
00871 1.10514232634061696926e0,
00872 5.07386386128601488557e0,
00873 8.39985554327604159757e0,
00874 5.20982848682361821619e0,
00875 9.99999999999999997461e-1,
00876 };
00877
00878 static const double b1_QP[8] = {
00879 5.10862594750176621635e-2,
00880 4.98213872951233449420e0,
00881 7.58238284132545283818e1,
00882 3.66779609360150777800e2,
00883 7.10856304998926107277e2,
00884 5.97489612400613639965e2,
00885 2.11688757100572135698e2,
00886 2.52070205858023719784e1,
00887 };
00888
00889 static const double b1_QQ[7] = {
00890
00891 7.42373277035675149943e1,
00892 1.05644886038262816351e3,
00893 4.98641058337653607651e3,
00894 9.56231892404756170795e3,
00895 7.99704160447350683650e3,
00896 2.82619278517639096600e3,
00897 3.36093607810698293419e2,
00898 };
00899
00900 static const double b1_YP[6] = {
00901 1.26320474790178026440e9,
00902 -6.47355876379160291031e11,
00903 1.14509511541823727583e14,
00904 -8.12770255501325109621e15,
00905 2.02439475713594898196e17,
00906 -7.78877196265950026825e17,
00907 };
00908
00909 static const double b1_YQ[8] = {
00910
00911 5.94301592346128195359E2,
00912 2.35564092943068577943E5,
00913 7.34811944459721705660E7,
00914 1.87601316108706159478E10,
00915 3.88231277496238566008E12,
00916 6.20557727146953693363E14,
00917 6.87141087355300489866E16,
00918 3.97270608116560655612E18,
00919 };
00920
00921 static const double Z1 = 1.46819706421238932572E1;
00922 static const double Z2 = 4.92184563216946036703E1;
00923
00924 static const double THPIO4 = 3.*PI/4.;
00925
00926 double bessel_j1(double x)
00927 {
00928 double w, z, p, q, xn;
00929
00930 DEBUG_ENTRY( "bessel_j1()" );
00931
00932 w = x;
00933 if( x < 0 )
00934 w = -x;
00935
00936 if( w <= 5.0 )
00937 {
00938 z = x * x;
00939 w = polevl( z, b1_RP, 3 ) / p1evl( z, b1_RQ, 8 );
00940 w = w * x * (z - Z1) * (z - Z2);
00941
00942 DEBUG_EXIT( "bessel_j1()" );
00943 return w;
00944 }
00945
00946 w = 5.0/x;
00947 z = w * w;
00948 p = polevl( z, b1_PP, 6)/polevl( z, b1_PQ, 6 );
00949 q = polevl( z, b1_QP, 7)/p1evl( z, b1_QQ, 7 );
00950 xn = x - THPIO4;
00951 p = p * cos(xn) - w * q * sin(xn);
00952
00953 DEBUG_EXIT( "bessel_j1()" );
00954 return p * SQ2OPI / sqrt(x);
00955 }
00956
00957 double bessel_y1(double x)
00958 {
00959 double w, z, p, q, xn;
00960
00961 DEBUG_ENTRY( "bessel_y1()" );
00962
00963 if( x <= 5.0 )
00964 {
00965 if( x <= 0.0 )
00966 {
00967 fprintf( ioQQQ, "bessel_y1: domain error\n" );
00968 puts( "[Stop in bessel_y1]" );
00969 cdEXIT(EXIT_FAILURE);
00970 }
00971 z = x * x;
00972 w = x * (polevl( z, b1_YP, 5 ) / p1evl( z, b1_YQ, 8 ));
00973 w += TWOOPI * ( bessel_j1(x) * log(x) - 1.0/x );
00974
00975 DEBUG_EXIT( "bessel_y1()" );
00976 return w;
00977 }
00978
00979 w = 5.0/x;
00980 z = w * w;
00981 p = polevl( z, b1_PP, 6 )/polevl( z, b1_PQ, 6 );
00982 q = polevl( z, b1_QP, 7 )/p1evl( z, b1_QQ, 7 );
00983 xn = x - THPIO4;
00984 p = p * sin(xn) + w * q * cos(xn);
00985
00986 DEBUG_EXIT( "bessel_y1()" );
00987 return p * SQ2OPI / sqrt(x);
00988 }
00989
00990
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038 double bessel_jn(int n, double x)
01039 {
01040 double pkm2, pkm1, pk, xk, r, ans;
01041 int k, sign;
01042
01043 DEBUG_ENTRY( "bessel_jn()" );
01044
01045 if( n < 0 )
01046 {
01047 n = -n;
01048 if( (n & 1) == 0 )
01049 sign = 1;
01050 else
01051 sign = -1;
01052 }
01053 else
01054 sign = 1;
01055
01056 if( x < 0.0 )
01057 {
01058 if( n & 1 )
01059 sign = -sign;
01060 x = -x;
01061 }
01062
01063 if( n == 0 )
01064 {
01065 DEBUG_EXIT( "bessel_jn()" );
01066 return sign * bessel_j0(x);
01067 }
01068 if( n == 1 )
01069 {
01070 DEBUG_EXIT( "bessel_jn()" );
01071 return sign * bessel_j1(x);
01072 }
01073 if( n == 2 )
01074 {
01075 DEBUG_EXIT( "bessel_jn()" );
01076 return sign * (2.0 * bessel_j1(x) / x - bessel_j0(x));
01077 }
01078
01079 if( x < DBL_EPSILON )
01080 {
01081 DEBUG_EXIT( "bessel_jn()" );
01082 return 0.0;
01083 }
01084
01085
01086 k = 52;
01087
01088 pk = 2 * (n + k);
01089 ans = pk;
01090 xk = x * x;
01091
01092 do
01093 {
01094 pk -= 2.0;
01095 ans = pk - (xk/ans);
01096 }
01097 while( --k > 0 );
01098 ans = x/ans;
01099
01100
01101 pk = 1.0;
01102 pkm1 = 1.0/ans;
01103 k = n-1;
01104 r = 2 * k;
01105
01106 do
01107 {
01108 pkm2 = (pkm1 * r - pk * x) / x;
01109 pk = pkm1;
01110 pkm1 = pkm2;
01111 r -= 2.0;
01112 }
01113 while( --k > 0 );
01114
01115 if( fabs(pk) > fabs(pkm1) )
01116 ans = bessel_j1(x)/pk;
01117 else
01118 ans = bessel_j0(x)/pkm1;
01119
01120 DEBUG_EXIT( "bessel_jn()" );
01121 return sign * ans;
01122 }
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178 double bessel_yn(int n, double x)
01179 {
01180 double an, anm1, anm2, r;
01181 int k, sign;
01182
01183 DEBUG_ENTRY( "bessel_yn()" );
01184
01185 if( n < 0 )
01186 {
01187 n = -n;
01188 if( (n & 1) == 0 )
01189 sign = 1;
01190 else
01191 sign = -1;
01192 }
01193 else
01194 sign = 1;
01195
01196 if( n == 0 )
01197 {
01198 DEBUG_EXIT( "bessel_yn()" );
01199 return sign * bessel_y0(x);
01200 }
01201 if( n == 1 )
01202 {
01203 DEBUG_EXIT( "bessel_yn()" );
01204 return sign * bessel_y1(x);
01205 }
01206
01207
01208 if( x <= 0.0 )
01209 {
01210 fprintf( ioQQQ, "bessel_yn: domain error\n" );
01211 puts( "[Stop in bessel_yn]" );
01212 cdEXIT(EXIT_FAILURE);
01213 }
01214
01215
01216 anm2 = bessel_y0(x);
01217 anm1 = bessel_y1(x);
01218 k = 1;
01219 r = 2 * k;
01220 do
01221 {
01222 an = r * anm1 / x - anm2;
01223 anm2 = anm1;
01224 anm1 = an;
01225 r += 2.0;
01226 ++k;
01227 }
01228 while( k < n );
01229
01230 DEBUG_EXIT( "bessel_yn()" );
01231 return sign * an;
01232 }
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317 static const double k0_A[] =
01318 {
01319 1.37446543561352307156e-16,
01320 4.25981614279661018399e-14,
01321 1.03496952576338420167e-11,
01322 1.90451637722020886025e-9,
01323 2.53479107902614945675e-7,
01324 2.28621210311945178607e-5,
01325 1.26461541144692592338e-3,
01326 3.59799365153615016266e-2,
01327 3.44289899924628486886e-1,
01328 -5.35327393233902768720e-1
01329 };
01330
01331
01332
01333
01334
01335
01336
01337 static const double k0_B[] = {
01338 5.30043377268626276149e-18,
01339 -1.64758043015242134646e-17,
01340 5.21039150503902756861e-17,
01341 -1.67823109680541210385e-16,
01342 5.51205597852431940784e-16,
01343 -1.84859337734377901440e-15,
01344 6.34007647740507060557e-15,
01345 -2.22751332699166985548e-14,
01346 8.03289077536357521100e-14,
01347 -2.98009692317273043925e-13,
01348 1.14034058820847496303e-12,
01349 -4.51459788337394416547e-12,
01350 1.85594911495471785253e-11,
01351 -7.95748924447710747776e-11,
01352 3.57739728140030116597e-10,
01353 -1.69753450938905987466e-9,
01354 8.57403401741422608519e-9,
01355 -4.66048989768794782956e-8,
01356 2.76681363944501510342e-7,
01357 -1.83175552271911948767e-6,
01358 1.39498137188764993662e-5,
01359 -1.28495495816278026384e-4,
01360 1.56988388573005337491e-3,
01361 -3.14481013119645005427e-2,
01362 2.44030308206595545468e0
01363 };
01364
01365 double bessel_k0(double x)
01366 {
01367 double y, z;
01368
01369 DEBUG_ENTRY( "bessel_k0()" );
01370
01371 if( x <= 0.0 )
01372 {
01373 fprintf( ioQQQ, "bessel_k0: domain error\n" );
01374 puts( "[Stop in bessel_k0]" );
01375 cdEXIT(EXIT_FAILURE);
01376 }
01377
01378 if( x <= 2.0 )
01379 {
01380 y = x * x - 2.0;
01381 y = chbevl( y, k0_A, 10 ) - log( 0.5 * x ) * bessel_i0(x);
01382
01383 DEBUG_EXIT( "bessel_k0()" );
01384 return y;
01385 }
01386 z = 8.0/x - 2.0;
01387 y = exp(-x) * chbevl( z, k0_B, 25 ) / sqrt(x);
01388
01389 DEBUG_EXIT( "bessel_k0()" );
01390 return y;
01391 }
01392
01393 double bessel_k0_scaled(double x)
01394 {
01395 double y;
01396
01397 DEBUG_ENTRY( "bessel_k0_scaled()" );
01398
01399 if( x <= 0.0 )
01400 {
01401 fprintf( ioQQQ, "bessel_k0_scaled: domain error\n" );
01402 puts( "[Stop in bessel_k0_scaled]" );
01403 cdEXIT(EXIT_FAILURE);
01404 }
01405
01406 if( x <= 2.0 )
01407 {
01408 y = x * x - 2.0;
01409 y = chbevl( y, k0_A, 10 ) - log( 0.5 * x ) * bessel_i0(x);
01410
01411 DEBUG_EXIT( "bessel_k0_scaled()" );
01412 return y * exp(x);
01413 }
01414
01415 DEBUG_EXIT( "bessel_k0_scaled()" );
01416 return chbevl( 8.0/x - 2.0, k0_B, 25 ) / sqrt(x);
01417 }
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465
01466
01467
01468
01469
01470
01471
01472
01473
01474
01475
01476
01477
01478
01479
01480
01481
01482
01483
01484
01485
01486
01487
01488
01489
01490
01491
01492
01493
01494
01495
01496
01497
01498
01499
01500
01501 static const double k1_A[] =
01502 {
01503 -7.02386347938628759343e-18,
01504 -2.42744985051936593393e-15,
01505 -6.66690169419932900609e-13,
01506 -1.41148839263352776110e-10,
01507 -2.21338763073472585583e-8,
01508 -2.43340614156596823496e-6,
01509 -1.73028895751305206302e-4,
01510 -6.97572385963986435018e-3,
01511 -1.22611180822657148235e-1,
01512 -3.53155960776544875667e-1,
01513 1.52530022733894777053e0
01514 };
01515
01516
01517
01518
01519
01520
01521
01522 static const double k1_B[] =
01523 {
01524 -5.75674448366501715755e-18,
01525 1.79405087314755922667e-17,
01526 -5.68946255844285935196e-17,
01527 1.83809354436663880070e-16,
01528 -6.05704724837331885336e-16,
01529 2.03870316562433424052e-15,
01530 -7.01983709041831346144e-15,
01531 2.47715442448130437068e-14,
01532 -8.97670518232499435011e-14,
01533 3.34841966607842919884e-13,
01534 -1.28917396095102890680e-12,
01535 5.13963967348173025100e-12,
01536 -2.12996783842756842877e-11,
01537 9.21831518760500529508e-11,
01538 -4.19035475934189648750e-10,
01539 2.01504975519703286596e-9,
01540 -1.03457624656780970260e-8,
01541 5.74108412545004946722e-8,
01542 -3.50196060308781257119e-7,
01543 2.40648494783721712015e-6,
01544 -1.93619797416608296024e-5,
01545 1.95215518471351631108e-4,
01546 -2.85781685962277938680e-3,
01547 1.03923736576817238437e-1,
01548 2.72062619048444266945e0
01549 };
01550
01551 double bessel_k1(double x)
01552 {
01553 double y, z;
01554
01555 DEBUG_ENTRY( "bessel_k1()" );
01556
01557 z = 0.5 * x;
01558 if( z <= 0.0 )
01559 {
01560 fprintf( ioQQQ, "bessel_k1: domain error\n" );
01561 puts( "[Stop in bessel_k1]" );
01562 cdEXIT(EXIT_FAILURE);
01563 }
01564
01565 if( x <= 2.0 )
01566 {
01567 y = x * x - 2.0;
01568 y = log(z) * bessel_i1(x) + chbevl( y, k1_A, 11 ) / x;
01569
01570 DEBUG_EXIT( "bessel_k1()" );
01571 return y;
01572 }
01573
01574 DEBUG_EXIT( "bessel_k1()" );
01575 return exp(-x) * chbevl( 8.0/x - 2.0, k1_B, 25 ) / sqrt(x);
01576 }
01577
01578 double bessel_k1_scaled(double x)
01579 {
01580 double y;
01581
01582 DEBUG_ENTRY( "bessel_k1_scaled()" );
01583
01584 if( x <= 0.0 )
01585 {
01586 fprintf( ioQQQ, "bessel_k1_scaled: domain error\n" );
01587 puts( "[Stop in bessel_k1_scaled]" );
01588 cdEXIT(EXIT_FAILURE);
01589 }
01590
01591 if( x <= 2.0 )
01592 {
01593 y = x * x - 2.0;
01594 y = log( 0.5 * x ) * bessel_i1(x) + chbevl( y, k1_A, 11 ) / x;
01595
01596 DEBUG_EXIT( "bessel_k1_scaled()" );
01597 return y * exp(x);
01598 }
01599
01600 DEBUG_EXIT( "bessel_k1_scaled()" );
01601 return chbevl( 8.0/x - 2.0, k1_B, 25 ) / sqrt(x);
01602 }
01603
01604
01605
01606
01607
01608
01609
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622
01623
01624
01625
01626
01627
01628
01629
01630
01631
01632
01633
01634
01635
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656
01657
01658
01659
01660
01661
01662
01663
01664
01665
01666
01667
01668
01669
01670
01671
01672
01673
01674
01675
01676
01677
01678
01679
01680
01681
01682
01683 static const double i0_A[] =
01684 {
01685 -4.41534164647933937950e-18,
01686 3.33079451882223809783e-17,
01687 -2.43127984654795469359e-16,
01688 1.71539128555513303061e-15,
01689 -1.16853328779934516808e-14,
01690 7.67618549860493561688e-14,
01691 -4.85644678311192946090e-13,
01692 2.95505266312963983461e-12,
01693 -1.72682629144155570723e-11,
01694 9.67580903537323691224e-11,
01695 -5.18979560163526290666e-10,
01696 2.65982372468238665035e-9,
01697 -1.30002500998624804212e-8,
01698 6.04699502254191894932e-8,
01699 -2.67079385394061173391e-7,
01700 1.11738753912010371815e-6,
01701 -4.41673835845875056359e-6,
01702 1.64484480707288970893e-5,
01703 -5.75419501008210370398e-5,
01704 1.88502885095841655729e-4,
01705 -5.76375574538582365885e-4,
01706 1.63947561694133579842e-3,
01707 -4.32430999505057594430e-3,
01708 1.05464603945949983183e-2,
01709 -2.37374148058994688156e-2,
01710 4.93052842396707084878e-2,
01711 -9.49010970480476444210e-2,
01712 1.71620901522208775349e-1,
01713 -3.04682672343198398683e-1,
01714 6.76795274409476084995e-1
01715 };
01716
01717
01718
01719
01720
01721
01722
01723 static const double i0_B[] =
01724 {
01725 -7.23318048787475395456e-18,
01726 -4.83050448594418207126e-18,
01727 4.46562142029675999901e-17,
01728 3.46122286769746109310e-17,
01729 -2.82762398051658348494e-16,
01730 -3.42548561967721913462e-16,
01731 1.77256013305652638360e-15,
01732 3.81168066935262242075e-15,
01733 -9.55484669882830764870e-15,
01734 -4.15056934728722208663e-14,
01735 1.54008621752140982691e-14,
01736 3.85277838274214270114e-13,
01737 7.18012445138366623367e-13,
01738 -1.79417853150680611778e-12,
01739 -1.32158118404477131188e-11,
01740 -3.14991652796324136454e-11,
01741 1.18891471078464383424e-11,
01742 4.94060238822496958910e-10,
01743 3.39623202570838634515e-9,
01744 2.26666899049817806459e-8,
01745 2.04891858946906374183e-7,
01746 2.89137052083475648297e-6,
01747 6.88975834691682398426e-5,
01748 3.36911647825569408990e-3,
01749 8.04490411014108831608e-1
01750 };
01751
01752 double bessel_i0(double x)
01753 {
01754 double y;
01755
01756 DEBUG_ENTRY( "bessel_i0()" );
01757
01758 if( x < 0 )
01759 x = -x;
01760
01761 if( x <= 8.0 )
01762 {
01763 y = (x/2.0) - 2.0;
01764
01765 DEBUG_EXIT( "bessel_i0()" );
01766 return exp(x) * chbevl( y, i0_A, 30 );
01767 }
01768
01769 DEBUG_EXIT( "bessel_i0()" );
01770 return exp(x) * chbevl( 32.0/x - 2.0, i0_B, 25 ) / sqrt(x);
01771 }
01772
01773 double bessel_i0_scaled(double x)
01774 {
01775 double y;
01776
01777 DEBUG_ENTRY( "bessel_i0_scaled()" );
01778
01779 if( x < 0 )
01780 x = -x;
01781
01782 if( x <= 8.0 )
01783 {
01784 y = (x/2.0) - 2.0;
01785
01786 DEBUG_EXIT( "bessel_i0_scaled()" );
01787 return chbevl( y, i0_A, 30 );
01788 }
01789
01790 DEBUG_EXIT( "bessel_i0_scaled()" );
01791 return chbevl( 32.0/x - 2.0, i0_B, 25 ) / sqrt(x);
01792 }
01793
01794
01795
01796
01797
01798
01799
01800
01801
01802
01803
01804
01805
01806
01807
01808
01809
01810
01811
01812
01813
01814
01815
01816
01817
01818
01819
01820
01821
01822
01823
01824
01825
01826
01827
01828
01829
01830
01831
01832
01833
01834
01835
01836
01837
01838
01839
01840
01841
01842
01843
01844
01845
01846
01847
01848
01849
01850
01851
01852
01853
01854
01855
01856
01857
01858
01859
01860
01861
01862
01863
01864
01865
01866
01867
01868
01869
01870
01871
01872
01873
01874 static double i1_A[] =
01875 {
01876 2.77791411276104639959e-18,
01877 -2.11142121435816608115e-17,
01878 1.55363195773620046921e-16,
01879 -1.10559694773538630805e-15,
01880 7.60068429473540693410e-15,
01881 -5.04218550472791168711e-14,
01882 3.22379336594557470981e-13,
01883 -1.98397439776494371520e-12,
01884 1.17361862988909016308e-11,
01885 -6.66348972350202774223e-11,
01886 3.62559028155211703701e-10,
01887 -1.88724975172282928790e-9,
01888 9.38153738649577178388e-9,
01889 -4.44505912879632808065e-8,
01890 2.00329475355213526229e-7,
01891 -8.56872026469545474066e-7,
01892 3.47025130813767847674e-6,
01893 -1.32731636560394358279e-5,
01894 4.78156510755005422638e-5,
01895 -1.61760815825896745588e-4,
01896 5.12285956168575772895e-4,
01897 -1.51357245063125314899e-3,
01898 4.15642294431288815669e-3,
01899 -1.05640848946261981558e-2,
01900 2.47264490306265168283e-2,
01901 -5.29459812080949914269e-2,
01902 1.02643658689847095384e-1,
01903 -1.76416518357834055153e-1,
01904 2.52587186443633654823e-1
01905 };
01906
01907
01908
01909
01910
01911
01912
01913 static double i1_B[] =
01914 {
01915 7.51729631084210481353e-18,
01916 4.41434832307170791151e-18,
01917 -4.65030536848935832153e-17,
01918 -3.20952592199342395980e-17,
01919 2.96262899764595013876e-16,
01920 3.30820231092092828324e-16,
01921 -1.88035477551078244854e-15,
01922 -3.81440307243700780478e-15,
01923 1.04202769841288027642e-14,
01924 4.27244001671195135429e-14,
01925 -2.10154184277266431302e-14,
01926 -4.08355111109219731823e-13,
01927 -7.19855177624590851209e-13,
01928 2.03562854414708950722e-12,
01929 1.41258074366137813316e-11,
01930 3.25260358301548823856e-11,
01931 -1.89749581235054123450e-11,
01932 -5.58974346219658380687e-10,
01933 -3.83538038596423702205e-9,
01934 -2.63146884688951950684e-8,
01935 -2.51223623787020892529e-7,
01936 -3.88256480887769039346e-6,
01937 -1.10588938762623716291e-4,
01938 -9.76109749136146840777e-3,
01939 7.78576235018280120474e-1
01940 };
01941
01942 double bessel_i1(double x)
01943 {
01944 double y, z;
01945
01946 DEBUG_ENTRY( "bessel_i1()" );
01947
01948 z = fabs(x);
01949 if( z <= 8.0 )
01950 {
01951 y = (z/2.0) - 2.0;
01952 z = chbevl( y, i1_A, 29 ) * z * exp(z);
01953 }
01954 else
01955 {
01956 z = exp(z) * chbevl( 32.0/z - 2.0, i1_B, 25 ) / sqrt(z);
01957 }
01958 if( x < 0.0 )
01959 z = -z;
01960
01961 DEBUG_EXIT( "bessel_i1()" );
01962 return z;
01963 }
01964
01965 double bessel_i1_scaled(double x)
01966 {
01967 double y, z;
01968
01969 DEBUG_ENTRY( "bessel_i1_scaled()" );
01970
01971 z = fabs(x);
01972 if( z <= 8.0 )
01973 {
01974 y = (z/2.0) - 2.0;
01975 z = chbevl( y, i1_A, 29 ) * z;
01976 }
01977 else
01978 {
01979 z = chbevl( 32.0/z - 2.0, i1_B, 25 ) / sqrt(z);
01980 }
01981 if( x < 0.0 )
01982 z = -z;
01983
01984 DEBUG_EXIT( "bessel_i1_scaled()" );
01985 return z;
01986 }
01987
01988
01989
01990
01991
01992
01993
01994
01995
01996
01997
01998
01999
02000
02001
02002
02003
02004
02005
02006
02007
02008
02009
02010
02011
02012
02013
02014
02015
02016
02017
02018
02019
02020
02021
02022
02023
02024
02025
02026
02027
02028
02029
02030
02031
02032
02033
02034
02035
02036
02037
02038
02039
02040
02041
02042
02043
02044
02045
02046
02047 static const double elk_P[] =
02048 {
02049 1.37982864606273237150e-4,
02050 2.28025724005875567385e-3,
02051 7.97404013220415179367e-3,
02052 9.85821379021226008714e-3,
02053 6.87489687449949877925e-3,
02054 6.18901033637687613229e-3,
02055 8.79078273952743772254e-3,
02056 1.49380448916805252718e-2,
02057 3.08851465246711995998e-2,
02058 9.65735902811690126535e-2,
02059 1.38629436111989062502e0
02060 };
02061
02062 static const double elk_Q[] =
02063 {
02064 2.94078955048598507511e-5,
02065 9.14184723865917226571e-4,
02066 5.94058303753167793257e-3,
02067 1.54850516649762399335e-2,
02068 2.39089602715924892727e-2,
02069 3.01204715227604046988e-2,
02070 3.73774314173823228969e-2,
02071 4.88280347570998239232e-2,
02072 7.03124996963957469739e-2,
02073 1.24999999999870820058e-1,
02074 4.99999999999999999821e-1
02075 };
02076
02077 static const double C1 = 1.3862943611198906188e0;
02078
02079 double ellpk(double x)
02080 {
02081 if( x < 0.0 || x > 1.0 )
02082 {
02083 fprintf( ioQQQ, "ellpk: domain error\n" );
02084 puts( "[Stop in ellpk]" );
02085 cdEXIT(EXIT_FAILURE);
02086 }
02087
02088 if( x > DBL_EPSILON )
02089 {
02090 return polevl(x,elk_P,10) - log(x) * polevl(x,elk_Q,10);
02091 }
02092 else
02093 {
02094 if( x == 0.0 )
02095 {
02096 fprintf( ioQQQ, "ellpk: domain error\n" );
02097 puts( "[Stop in ellpk]" );
02098 cdEXIT(EXIT_FAILURE);
02099 }
02100 else
02101 {
02102 return C1 - 0.5 * log(x);
02103 }
02104 }
02105 }
02106
02107
02108
02109
02110
02111
02112
02113
02114
02115
02116
02117
02118
02119
02120
02121
02122
02123
02124
02125
02126
02127
02128
02129
02130
02131
02132
02133
02134
02135
02136
02137
02138
02139
02140
02141
02142
02143
02144
02145
02146
02147
02148
02149
02150
02151
02152
02153
02154
02155 static const double MAXLOG = log(DBL_MAX);
02156 static const double BIG = 1.44115188075855872E+17;
02157
02158 double expn(int n, double x)
02159 {
02160 double ans, r, t, yk, xk;
02161 double pk, pkm1, pkm2, qk, qkm1, qkm2;
02162 double psi, z;
02163 int i, k;
02164
02165 DEBUG_ENTRY( "expn()" );
02166
02167 if( n < 0 || x < 0. )
02168 {
02169 fprintf( ioQQQ, "expn: domain error\n" );
02170 puts( "[Stop in expn]" );
02171 cdEXIT(EXIT_FAILURE);
02172 }
02173
02174 if( x > MAXLOG )
02175 {
02176 DEBUG_EXIT( "expn()" );
02177 return 0.0;
02178 }
02179
02180 if( x == 0.0 )
02181 {
02182 if( n < 2 )
02183 {
02184 fprintf( ioQQQ, "expn: domain error\n" );
02185 puts( "[Stop in expn]" );
02186 cdEXIT(EXIT_FAILURE);
02187 }
02188 else
02189 {
02190 DEBUG_EXIT( "expn()" );
02191 return 1.0/((double)n-1.0);
02192 }
02193 }
02194
02195 if( n == 0 )
02196 {
02197 DEBUG_EXIT( "expn()" );
02198 return exp(-x)/x;
02199 }
02200
02201
02202 if( n > 5000 )
02203 {
02204 xk = x + n;
02205 yk = 1.0 / (xk * xk);
02206 t = n;
02207 ans = yk * t * (6.0 * x * x - 8.0 * t * x + t * t);
02208 ans = yk * (ans + t * (t - 2.0 * x));
02209 ans = yk * (ans + t);
02210 ans = (ans + 1.0) * exp( -x ) / xk;
02211
02212 DEBUG_EXIT( "expn()" );
02213 return ans;
02214 }
02215
02216 if( x <= 1.0 )
02217 {
02218
02219 psi = -EULER - log(x);
02220 for( i=1; i < n; i++ )
02221 psi = psi + 1.0/i;
02222
02223 z = -x;
02224 xk = 0.0;
02225 yk = 1.0;
02226 pk = 1.0 - n;
02227 if( n == 1 )
02228 ans = 0.0;
02229 else
02230 ans = 1.0/pk;
02231 do
02232 {
02233 xk += 1.0;
02234 yk *= z/xk;
02235 pk += 1.0;
02236 if( pk != 0.0 )
02237 {
02238 ans += yk/pk;
02239 }
02240 if( ans != 0.0 )
02241 t = fabs(yk/ans);
02242 else
02243 t = 1.0;
02244 }
02245 while( t > DBL_EPSILON );
02246 ans = powi(z,n-1)*psi/factorial(n-1) - ans;
02247
02248 DEBUG_EXIT( "expn()" );
02249 return ans;
02250 }
02251 else
02252 {
02253
02254 k = 1;
02255 pkm2 = 1.0;
02256 qkm2 = x;
02257 pkm1 = 1.0;
02258 qkm1 = x + n;
02259 ans = pkm1/qkm1;
02260
02261 do
02262 {
02263 k += 1;
02264 if( is_odd(k) )
02265 {
02266 yk = 1.0;
02267 xk = static_cast<double>(n + (k-1)/2);
02268 }
02269 else
02270 {
02271 yk = x;
02272 xk = static_cast<double>(k/2);
02273 }
02274 pk = pkm1 * yk + pkm2 * xk;
02275 qk = qkm1 * yk + qkm2 * xk;
02276 if( qk != 0 )
02277 {
02278 r = pk/qk;
02279 t = fabs( (ans - r)/r );
02280 ans = r;
02281 }
02282 else
02283 t = 1.0;
02284 pkm2 = pkm1;
02285 pkm1 = pk;
02286 qkm2 = qkm1;
02287 qkm1 = qk;
02288 if( fabs(pk) > BIG )
02289 {
02290 pkm2 /= BIG;
02291 pkm1 /= BIG;
02292 qkm2 /= BIG;
02293 qkm1 /= BIG;
02294 }
02295 }
02296 while( t > DBL_EPSILON );
02297
02298 ans *= exp( -x );
02299
02300 DEBUG_EXIT( "expn()" );
02301 return ans;
02302 }
02303 }
02304
02305
02306
02307
02308
02309
02310
02311
02312
02313
02314
02315
02316
02317
02318
02319
02320
02321
02322
02323
02324
02325
02326
02327
02328
02329
02330
02331
02332
02333
02334
02335
02336
02337
02338
02339
02340
02341
02342
02343
02344
02345
02346
02347
02348
02349
02350
02351
02352
02353
02354
02355 inline double polevl(double x, const double coef[], int N)
02356 {
02357 double ans;
02358 int i;
02359 const double *p = coef;
02360
02361 ans = *p++;
02362 i = N;
02363
02364 do
02365 ans = ans * x + *p++;
02366 while( --i );
02367
02368 return ans;
02369 }
02370
02371
02372
02373
02374
02375
02376
02377 inline double p1evl(double x, const double coef[], int N)
02378 {
02379 double ans;
02380 const double *p = coef;
02381 int i;
02382
02383 ans = x + *p++;
02384 i = N-1;
02385
02386 do
02387 ans = ans * x + *p++;
02388 while( --i );
02389
02390 return ans;
02391 }
02392
02393
02394
02395
02396
02397
02398
02399
02400
02401
02402
02403
02404
02405
02406
02407
02408
02409
02410
02411
02412
02413
02414
02415
02416
02417
02418
02419
02420
02421
02422
02423
02424
02425
02426
02427
02428
02429
02430
02431
02432
02433
02434
02435
02436
02437
02438
02439
02440
02441
02442
02443
02444
02445
02446
02447
02448
02449
02450
02451 inline double chbevl(double x, const double array[], int n)
02452 {
02453 double b0, b1, b2;
02454 const double *p = array;
02455 int i;
02456
02457 b0 = *p++;
02458 b1 = 0.0;
02459 i = n - 1;
02460
02461 do
02462 {
02463 b2 = b1;
02464 b1 = b0;
02465 b0 = x * b1 - b2 + *p++;
02466 }
02467 while( --i );
02468
02469 return 0.5*(b0-b2);
02470 }