MercuryDPM  Beta
besselFunc Namespace Reference

Namespace for evaluating the zeroth modified Bessel function of the first kind, I0(x), required in StatisticsPoint.hcc. More...

## Functions

Mdouble chebyshev (Mdouble x, const Mdouble coef[], int N)

Mdouble I0_exp (Mdouble x)

Mdouble I0 (Mdouble x)

## Detailed Description

Namespace for evaluating the zeroth modified Bessel function of the first kind, I0(x), required in StatisticsPoint.hcc.

## Function Documentation

 Mdouble besselFunc::chebyshev ( Mdouble x, const Mdouble coef[], int N )

Definition at line 159 of file ExtendedMath.cc.

Referenced by I0_exp().

160 {
161  const Mdouble *p = coef;
162  Mdouble b0 = *p++;
163  Mdouble b1 = 0, b2;
164  int i = N - 1;
165
166  assert(i > 0);
167  do
168  {
169  b2 = b1;
170  b1 = b0;
171  b0 = x * b1 - b2 + *p++;
172  }
173  while ( --i);
174
175  return (0.5 * (b0 - b2));
176 }
double Mdouble
 Mdouble besselFunc::I0 ( Mdouble x )

Definition at line 257 of file ExtendedMath.cc.

References I0_exp().

258 {
259  if (x < 0)
260  x = -x;
261  return exp(x) * I0_exp(x);
262 }
Mdouble I0_exp(Mdouble x)
 Mdouble besselFunc::I0_exp ( Mdouble x )

Definition at line 178 of file ExtendedMath.cc.

References A, and chebyshev().

Referenced by I0().

179 {
180  // Cooefficients for [0..8]
181  const Mdouble A[] =
182  {
183  -4.415341646479339379501E-18,
184  3.330794518822238097831E-17,
185  -2.431279846547954693591E-16,
186  1.715391285555133030611E-15,
187  -1.168533287799345168081E-14,
188  7.676185498604935616881E-14,
189  -4.856446783111929460901E-13,
190  2.955052663129639834611E-12,
191  -1.726826291441555707231E-11,
192  9.675809035373236912241E-11,
193  -5.189795601635262906661E-10,
194  2.659823724682386650351E-9,
195  -1.300025009986248042121E-8,
196  6.046995022541918949321E-8,
197  -2.670793853940611733911E-7,
198  1.117387539120103718151E-6,
199  -4.416738358458750563591E-6,
200  1.644844807072889708931E-5,
201  -5.754195010082103703981E-5,
202  1.885028850958416557291E-4,
203  -5.763755745385823658851E-4,
204  1.639475616941335798421E-3,
205  -4.324309995050575944301E-3,
206  1.054646039459499831831E-2,
207  -2.373741480589946881561E-2,
208  4.930528423967070848781E-2,
209  -9.490109704804764442101E-2,
210  1.716209015222087753491E-1,
211  -3.046826723431983986831E-1,
212  6.767952744094760849951E-1
213  };
214
215  // Cooefficients for [8..infinity]
216  const Mdouble B[] =
217  {
218  -7.233180487874753954561E-18,
219  -4.830504485944182071261E-18,
220  4.465621420296759999011E-17,
221  3.461222867697461093101E-17,
222  -2.827623980516583484941E-16,
223  -3.425485619677219134621E-16,
224  1.772560133056526383601E-15,
225  3.811680669352622420751E-15,
226  -9.554846698828307648701E-15,
227  -4.150569347287222086631E-14,
228  1.540086217521409826911E-14,
229  3.852778382742142701141E-13,
230  7.180124451383666233671E-13,
231  -1.794178531506806117781E-12,
232  -1.321581184044771311881E-11,
233  -3.149916527963241364541E-11,
234  1.188914710784643834241E-11,
235  4.940602388224969589101E-10,
236  3.396232025708386345151E-9,
237  2.266668990498178064591E-8,
238  2.048918589469063741831E-7,
239  2.891370520834756482971E-6,
240  6.889758346916823984261E-5,
241  3.369116478255694089901E-3,
242  8.044904110141088316081E-1
243  };
244
245  if (x < 0)
246  x = -x;
247
248  if (x <= 8.0)
249  {
250  Mdouble y = (x / 2.0) - 2.0;
251  return (chebyshev(y, A, 30));
252  }
253
254  return (chebyshev(32.0 / x - 2.0, B, 25) / sqrt(x));
255 }
Mdouble chebyshev(Mdouble x, const Mdouble coef[], int N)
double Mdouble