GNU Octave  4.2.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
besselj.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1997-2017 John W. Eaton
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #if defined (HAVE_CONFIG_H)
24 # include "config.h"
25 #endif
26 
27 #include "lo-specfun.h"
28 #include "quit.h"
29 
30 #include "defun.h"
31 #include "error.h"
32 #include "ovl.h"
33 #include "utils.h"
34 
35 enum bessel_type
36 {
37  BESSEL_J,
38  BESSEL_Y,
39  BESSEL_I,
40  BESSEL_K,
41  BESSEL_H1,
42  BESSEL_H2
43 };
44 
45 #define DO_BESSEL(type, alpha, x, scaled, ierr, result) \
46  do \
47  { \
48  switch (type) \
49  { \
50  case BESSEL_J: \
51  result = octave::math::besselj (alpha, x, scaled, ierr); \
52  break; \
53  \
54  case BESSEL_Y: \
55  result = octave::math::bessely (alpha, x, scaled, ierr); \
56  break; \
57  \
58  case BESSEL_I: \
59  result = octave::math::besseli (alpha, x, scaled, ierr); \
60  break; \
61  \
62  case BESSEL_K: \
63  result = octave::math::besselk (alpha, x, scaled, ierr); \
64  break; \
65  \
66  case BESSEL_H1: \
67  result = octave::math::besselh1 (alpha, x, scaled, ierr); \
68  break; \
69  \
70  case BESSEL_H2: \
71  result = octave::math::besselh2 (alpha, x, scaled, ierr); \
72  break; \
73  \
74  default: \
75  break; \
76  } \
77  } \
78  while (0)
79 
80 octave_value_list
81 do_bessel (enum bessel_type type, const char *fn,
82  const octave_value_list& args, int nargout)
83 {
84  int nargin = args.length ();
85 
86  if (nargin < 2 || nargin > 3)
87  print_usage ();
88 
89  octave_value_list retval (nargout > 1 ? 2 : 1);
90 
91  bool scaled = false;
92  if (nargin == 3)
93  {
94  octave_value opt_arg = args(2);
95  bool rpt_error = false;
96 
97  if (! opt_arg.is_scalar_type ())
98  rpt_error = true;
99  else if (opt_arg.is_numeric_type ())
100  {
101  double opt_val = opt_arg.double_value ();
102  if (opt_val != 0.0 && opt_val != 1.0)
103  rpt_error = true;
104  scaled = (opt_val == 1.0);
105  }
106  else if (opt_arg.is_bool_type ())
107  scaled = opt_arg.bool_value ();
108 
109  if (rpt_error)
110  error ("%s: OPT must be 0 (or false) or 1 (or true)", fn);
111  }
112 
113  octave_value alpha_arg = args(0);
114  octave_value x_arg = args(1);
115 
116  if (alpha_arg.is_single_type () || x_arg.is_single_type ())
117  {
118  if (alpha_arg.is_scalar_type ())
119  {
120  float alpha = args(0).xfloat_value ("%s: ALPHA must be a scalar or matrix", fn);
121 
122  if (x_arg.is_scalar_type ())
123  {
124  FloatComplex x = x_arg.xfloat_complex_value ("%s: X must be a scalar or matrix", fn);
125 
126  octave_idx_type ierr;
127  octave_value result;
128 
129  DO_BESSEL (type, alpha, x, scaled, ierr, result);
130 
131  retval(0) = result;
132  if (nargout > 1)
133  retval(1) = static_cast<float> (ierr);
134  }
135  else
136  {
137  FloatComplexNDArray x = x_arg.xfloat_complex_array_value ("%s: X must be a scalar or matrix", fn);
138 
139  Array<octave_idx_type> ierr;
140  octave_value result;
141 
142  DO_BESSEL (type, alpha, x, scaled, ierr, result);
143 
144  retval(0) = result;
145  if (nargout > 1)
146  retval(1) = NDArray (ierr);
147  }
148  }
149  else
150  {
151  dim_vector dv0 = args(0).dims ();
152  dim_vector dv1 = args(1).dims ();
153 
154  bool args0_is_row_vector = (dv0(1) == dv0.numel ());
155  bool args1_is_col_vector = (dv1(0) == dv1.numel ());
156 
157  if (args0_is_row_vector && args1_is_col_vector)
158  {
159  FloatRowVector ralpha = args(0).xfloat_row_vector_value ("%s: ALPHA must be a scalar or matrix", fn);
160 
161  FloatComplexColumnVector cx =
162  x_arg.xfloat_complex_column_vector_value ("%s: X must be a scalar or matrix", fn);
163 
164  Array<octave_idx_type> ierr;
165  octave_value result;
166 
167  DO_BESSEL (type, ralpha, cx, scaled, ierr, result);
168 
169  retval(0) = result;
170  if (nargout > 1)
171  retval(1) = NDArray (ierr);
172  }
173  else
174  {
175  FloatNDArray alpha = args(0).xfloat_array_value ("%s: ALPHA must be a scalar or matrix", fn);
176 
177  if (x_arg.is_scalar_type ())
178  {
179  FloatComplex x = x_arg.xfloat_complex_value ("%s: X must be a scalar or matrix", fn);
180 
181  Array<octave_idx_type> ierr;
182  octave_value result;
183 
184  DO_BESSEL (type, alpha, x, scaled, ierr, result);
185 
186  retval(0) = result;
187  if (nargout > 1)
188  retval(1) = NDArray (ierr);
189  }
190  else
191  {
192  FloatComplexNDArray x = x_arg.xfloat_complex_array_value ("%s: X must be a scalar or matrix", fn);
193 
194  Array<octave_idx_type> ierr;
195  octave_value result;
196 
197  DO_BESSEL (type, alpha, x, scaled, ierr, result);
198 
199  retval(0) = result;
200  if (nargout > 1)
201  retval(1) = NDArray (ierr);
202  }
203  }
204  }
205  }
206  else
207  {
208  if (alpha_arg.is_scalar_type ())
209  {
210  double alpha = args(0).xdouble_value ("%s: ALPHA must be a scalar or matrix", fn);
211 
212  if (x_arg.is_scalar_type ())
213  {
214  Complex x = x_arg.xcomplex_value ("%s: X must be a scalar or matrix", fn);
215 
216  octave_idx_type ierr;
217  octave_value result;
218 
219  DO_BESSEL (type, alpha, x, scaled, ierr, result);
220 
221  retval(0) = result;
222  if (nargout > 1)
223  retval(1) = static_cast<double> (ierr);
224  }
225  else
226  {
227  ComplexNDArray x = x_arg.xcomplex_array_value ("%s: X must be a scalar or matrix", fn);
228 
229  Array<octave_idx_type> ierr;
230  octave_value result;
231 
232  DO_BESSEL (type, alpha, x, scaled, ierr, result);
233 
234  retval(0) = result;
235  if (nargout > 1)
236  retval(1) = NDArray (ierr);
237  }
238  }
239  else
240  {
241  dim_vector dv0 = args(0).dims ();
242  dim_vector dv1 = args(1).dims ();
243 
244  bool args0_is_row_vector = (dv0(1) == dv0.numel ());
245  bool args1_is_col_vector = (dv1(0) == dv1.numel ());
246 
247  if (args0_is_row_vector && args1_is_col_vector)
248  {
249  RowVector ralpha = args(0).xrow_vector_value ("%s: ALPHA must be a scalar or matrix", fn);
250 
251  ComplexColumnVector cx =
252  x_arg.xcomplex_column_vector_value ("%s: X must be a scalar or matrix", fn);
253 
254  Array<octave_idx_type> ierr;
255  octave_value result;
256 
257  DO_BESSEL (type, ralpha, cx, scaled, ierr, result);
258 
259  retval(0) = result;
260  if (nargout > 1)
261  retval(1) = NDArray (ierr);
262  }
263  else
264  {
265  NDArray alpha = args(0).xarray_value ("%s: ALPHA must be a scalar or matrix", fn);
266 
267  if (x_arg.is_scalar_type ())
268  {
269  Complex x = x_arg.xcomplex_value ("%s: X must be a scalar or matrix", fn);
270 
271  Array<octave_idx_type> ierr;
272  octave_value result;
273 
274  DO_BESSEL (type, alpha, x, scaled, ierr, result);
275 
276  retval(0) = result;
277  if (nargout > 1)
278  retval(1) = NDArray (ierr);
279  }
280  else
281  {
282  ComplexNDArray x = x_arg.xcomplex_array_value ("%s: X must be a scalar or matrix", fn);
283 
284  Array<octave_idx_type> ierr;
285  octave_value result;
286 
287  DO_BESSEL (type, alpha, x, scaled, ierr, result);
288 
289  retval(0) = result;
290  if (nargout > 1)
291  retval(1) = NDArray (ierr);
292  }
293  }
294  }
295  }
296 
297  return retval;
298 }
299 
300 DEFUN (besselj, args, nargout,
301  doc: /* -*- texinfo -*-
302 @deftypefn {} {[@var{j}, @var{ierr}] =} besselj (@var{alpha}, @var{x}, @var{opt})
303 @deftypefnx {} {[@var{y}, @var{ierr}] =} bessely (@var{alpha}, @var{x}, @var{opt})
304 @deftypefnx {} {[@var{i}, @var{ierr}] =} besseli (@var{alpha}, @var{x}, @var{opt})
305 @deftypefnx {} {[@var{k}, @var{ierr}] =} besselk (@var{alpha}, @var{x}, @var{opt})
306 @deftypefnx {} {[@var{h}, @var{ierr}] =} besselh (@var{alpha}, @var{k}, @var{x}, @var{opt})
307 Compute Bessel or Hankel functions of various kinds:
308 
309 @table @code
310 @item besselj
311 Bessel functions of the first kind. If the argument @var{opt} is 1 or true,
312 the result is multiplied by @w{@code{exp (-abs (imag (@var{x})))}}.
313 
314 @item bessely
315 Bessel functions of the second kind. If the argument @var{opt} is 1 or
316 true, the result is multiplied by @code{exp (-abs (imag (@var{x})))}.
317 
318 @item besseli
319 
320 Modified Bessel functions of the first kind. If the argument @var{opt} is 1
321 or true, the result is multiplied by @code{exp (-abs (real (@var{x})))}.
322 
323 @item besselk
324 
325 Modified Bessel functions of the second kind. If the argument @var{opt} is
326 1 or true, the result is multiplied by @code{exp (@var{x})}.
327 
328 @item besselh
329 Compute Hankel functions of the first (@var{k} = 1) or second (@var{k}
330 = 2) kind. If the argument @var{opt} is 1 or true, the result is multiplied
331 by @code{exp (-I*@var{x})} for @var{k} = 1 or @code{exp (I*@var{x})} for
332 @var{k} = 2.
333 @end table
334 
335 If @var{alpha} is a scalar, the result is the same size as @var{x}.
336 If @var{x} is a scalar, the result is the same size as @var{alpha}.
337 If @var{alpha} is a row vector and @var{x} is a column vector, the
338 result is a matrix with @code{length (@var{x})} rows and
339 @code{length (@var{alpha})} columns. Otherwise, @var{alpha} and
340 @var{x} must conform and the result will be the same size.
341 
342 The value of @var{alpha} must be real. The value of @var{x} may be
343 complex.
344 
345 If requested, @var{ierr} contains the following status information
346 and is the same size as the result.
347 
348 @enumerate 0
349 @item
350 Normal return.
351 
352 @item
353 Input error, return @code{NaN}.
354 
355 @item
356 Overflow, return @code{Inf}.
357 
358 @item
359 Loss of significance by argument reduction results in less than
360 half of machine accuracy.
361 
362 @item
363 Complete loss of significance by argument reduction, return @code{NaN}.
364 
365 @item
366 Error---no computation, algorithm termination condition not met,
367 return @code{NaN}.
368 @end enumerate
369 @end deftypefn */)
370 {
371  return do_bessel (BESSEL_J, "besselj", args, nargout);
372 }
373 
374 /*
375 %!# Function besselj is tested along with other bessels at the end of this file
376 */
377 
378 DEFUN (bessely, args, nargout,
379  doc: /* -*- texinfo -*-
380 @deftypefn {} {[@var{y}, @var{ierr}] =} bessely (@var{alpha}, @var{x}, @var{opt})
381 See besselj.
382 @end deftypefn */)
383 {
384  return do_bessel (BESSEL_Y, "bessely", args, nargout);
385 }
386 
387 /*
388 %!# Function bessely is tested along with other bessels at the end of this file
389 */
390 
391 DEFUN (besseli, args, nargout,
392  doc: /* -*- texinfo -*-
393 @deftypefn {} {[@var{i}, @var{ierr}] =} besseli (@var{alpha}, @var{x}, @var{opt})
394 See besselj.
395 @end deftypefn */)
396 {
397  return do_bessel (BESSEL_I, "besseli", args, nargout);
398 }
399 
400 /*
401 %!# Function besseli is tested along with other bessels at the end of this file
402 */
403 
404 DEFUN (besselk, args, nargout,
405  doc: /* -*- texinfo -*-
406 @deftypefn {} {[@var{k}, @var{ierr}] =} besselk (@var{alpha}, @var{x}, @var{opt})
407 See besselj.
408 @end deftypefn */)
409 {
410  return do_bessel (BESSEL_K, "besselk", args, nargout);
411 }
412 
413 /*
414 %!# Function besselk is tested along with other bessels at the end of this file
415 */
416 
417 DEFUN (besselh, args, nargout,
418  doc: /* -*- texinfo -*-
419 @deftypefn {} {[@var{h}, @var{ierr}] =} besselh (@var{alpha}, @var{k}, @var{x}, @var{opt})
420 See besselj.
421 @end deftypefn */)
422 {
423  int nargin = args.length ();
424 
425  if (nargin < 2 || nargin > 4)
426  print_usage ();
427 
428  octave_value_list retval;
429 
430  if (nargin == 2)
431  {
432  retval = do_bessel (BESSEL_H1, "besselh", args, nargout);
433  }
434  else
435  {
436  octave_idx_type kind = args(1).xint_value ("besselh: invalid value of K");
437 
438  octave_value_list tmp_args;
439 
440  if (nargin == 4)
441  tmp_args(2) = args(3);
442 
443  tmp_args(1) = args(2);
444  tmp_args(0) = args(0);
445 
446  if (kind == 1)
447  retval = do_bessel (BESSEL_H1, "besselh", tmp_args, nargout);
448  else if (kind == 2)
449  retval = do_bessel (BESSEL_H2, "besselh", tmp_args, nargout);
450  else
451  error ("besselh: K must be 1 or 2");
452  }
453 
454  return retval;
455 }
456 
457 /*
458 %!# Function besselh is tested along with other bessels at the end of this file
459 */
460 
461 DEFUN (airy, args, nargout,
462  doc: /* -*- texinfo -*-
463 @deftypefn {} {[@var{a}, @var{ierr}] =} airy (@var{k}, @var{z}, @var{opt})
464 Compute Airy functions of the first and second kind, and their derivatives.
465 
466 @example
467 @group
468  K Function Scale factor (if "opt" is supplied)
469 --- -------- ---------------------------------------
470  0 Ai (Z) exp ((2/3) * Z * sqrt (Z))
471  1 dAi(Z)/dZ exp ((2/3) * Z * sqrt (Z))
472  2 Bi (Z) exp (-abs (real ((2/3) * Z * sqrt (Z))))
473  3 dBi(Z)/dZ exp (-abs (real ((2/3) * Z * sqrt (Z))))
474 @end group
475 @end example
476 
477 The function call @code{airy (@var{z})} is equivalent to
478 @code{airy (0, @var{z})}.
479 
480 The result is the same size as @var{z}.
481 
482 If requested, @var{ierr} contains the following status information and
483 is the same size as the result.
484 
485 @enumerate 0
486 @item
487 Normal return.
488 
489 @item
490 Input error, return @code{NaN}.
491 
492 @item
493 Overflow, return @code{Inf}.
494 
495 @item
496 Loss of significance by argument reduction results in less than half
497  of machine accuracy.
498 
499 @item
500 Complete loss of significance by argument reduction, return @code{NaN}.
501 
502 @item
503 Error---no computation, algorithm termination condition not met,
504 return @code{NaN}.
505 @end enumerate
506 @end deftypefn */)
507 {
508  int nargin = args.length ();
509 
510  if (nargin < 1 || nargin > 3)
511  print_usage ();
512 
513  octave_value_list retval (nargout > 1 ? 2 : 1);
514 
515  int kind = 0;
516  if (nargin > 1)
517  {
518  kind = args(0).xint_value ("airy: K must be an integer value");
519 
520  if (kind < 0 || kind > 3)
521  error ("airy: K must be 0, 1, 2, or 3");
522  }
523 
524  bool scale = (nargin == 3);
525 
526  int idx = (nargin == 1 ? 0 : 1);
527 
528  if (args(idx).is_single_type ())
529  {
530  FloatComplexNDArray z = args(idx).xfloat_complex_array_value ("airy: Z must be a complex matrix");
531 
532  Array<octave_idx_type> ierr;
533  octave_value result;
534 
535  if (kind > 1)
536  result = octave::math::biry (z, kind == 3, scale, ierr);
537  else
538  result = octave::math::airy (z, kind == 1, scale, ierr);
539 
540  retval(0) = result;
541  if (nargout > 1)
542  retval(1) = NDArray (ierr);
543  }
544  else
545  {
546  ComplexNDArray z = args(idx).xcomplex_array_value ("airy: Z must be a complex matrix");
547 
548  Array<octave_idx_type> ierr;
549  octave_value result;
550 
551  if (kind > 1)
552  result = octave::math::biry (z, kind == 3, scale, ierr);
553  else
554  result = octave::math::airy (z, kind == 1, scale, ierr);
555 
556  retval(0) = result;
557  if (nargout > 1)
558  retval(1) = NDArray (ierr);
559  }
560 
561  return retval;
562 }
563 
564 /*
565 FIXME: Function airy does not yet have BIST tests
566 */
567 
568 /*
569 ## Test values computed with GP/PARI version 2.3.3
570 %!shared alpha, x, jx, yx, ix, kx, nix
571 %!
572 %! ## Bessel functions, even order, positive and negative x
573 %! alpha = 2; x = 1.25;
574 %! jx = 0.1710911312405234823613091417;
575 %! yx = -1.193199310178553861283790424;
576 %! ix = 0.2220184483766341752692212604;
577 %! kx = 0.9410016167388185767085460540;
578 %!
579 %!assert (besselj (alpha,x), jx, 100*eps)
580 %!assert (bessely (alpha,x), yx, 100*eps)
581 %!assert (besseli (alpha,x), ix, 100*eps)
582 %!assert (besselk (alpha,x), kx, 100*eps)
583 %!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
584 %!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
585 %!
586 %!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
587 %!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
588 %!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
589 %!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
590 %!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
591 %!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
592 %!
593 %!assert (besselj (-alpha,x), jx, 100*eps)
594 %!assert (bessely (-alpha,x), yx, 100*eps)
595 %!assert (besseli (-alpha,x), ix, 100*eps)
596 %!assert (besselk (-alpha,x), kx, 100*eps)
597 %!assert (besselh (-alpha,1,x), jx + I*yx, 100*eps)
598 %!assert (besselh (-alpha,2,x), jx - I*yx, 100*eps)
599 %!
600 %!assert (besselj (-alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
601 %!assert (bessely (-alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
602 %!assert (besseli (-alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
603 %!assert (besselk (-alpha,x,1), kx*exp(x), 100*eps)
604 %!assert (besselh (-alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
605 %!assert (besselh (-alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
606 %!
607 %! x *= -1;
608 %! yx = -1.193199310178553861283790424 + 0.3421822624810469647226182835*I;
609 %! kx = 0.9410016167388185767085460540 - 0.6974915263814386815610060884*I;
610 %!
611 %!assert (besselj (alpha,x), jx, 100*eps)
612 %!assert (bessely (alpha,x), yx, 100*eps)
613 %!assert (besseli (alpha,x), ix, 100*eps)
614 %!assert (besselk (alpha,x), kx, 100*eps)
615 %!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
616 %!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
617 %!
618 %!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
619 %!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
620 %!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
621 %!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
622 %!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
623 %!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
624 %!
625 %! ## Bessel functions, odd order, positive and negative x
626 %! alpha = 3; x = 2.5;
627 %! jx = 0.2166003910391135247666890035;
628 %! yx = -0.7560554967536709968379029772;
629 %! ix = 0.4743704087780355895548240179;
630 %! kx = 0.2682271463934492027663765197;
631 %!
632 %!assert (besselj (alpha,x), jx, 100*eps)
633 %!assert (bessely (alpha,x), yx, 100*eps)
634 %!assert (besseli (alpha,x), ix, 100*eps)
635 %!assert (besselk (alpha,x), kx, 100*eps)
636 %!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
637 %!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
638 %!
639 %!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
640 %!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
641 %!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
642 %!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
643 %!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
644 %!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
645 %!
646 %!assert (besselj (-alpha,x), -jx, 100*eps)
647 %!assert (bessely (-alpha,x), -yx, 100*eps)
648 %!assert (besseli (-alpha,x), ix, 100*eps)
649 %!assert (besselk (-alpha,x), kx, 100*eps)
650 %!assert (besselh (-alpha,1,x), -(jx + I*yx), 100*eps)
651 %!assert (besselh (-alpha,2,x), -(jx - I*yx), 100*eps)
652 %!
653 %!assert (besselj (-alpha,x,1), -jx*exp(-abs(imag(x))), 100*eps)
654 %!assert (bessely (-alpha,x,1), -yx*exp(-abs(imag(x))), 100*eps)
655 %!assert (besseli (-alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
656 %!assert (besselk (-alpha,x,1), kx*exp(x), 100*eps)
657 %!assert (besselh (-alpha,1,x,1), -(jx + I*yx)*exp(-I*x), 100*eps)
658 %!assert (besselh (-alpha,2,x,1), -(jx - I*yx)*exp(I*x), 100*eps)
659 %!
660 %! x *= -1;
661 %! jx = -jx;
662 %! yx = 0.7560554967536709968379029772 - 0.4332007820782270495333780070*I;
663 %! ix = -ix;
664 %! kx = -0.2682271463934492027663765197 - 1.490278591297463775542004240*I;
665 %!
666 %!assert (besselj (alpha,x), jx, 100*eps)
667 %!assert (bessely (alpha,x), yx, 100*eps)
668 %!assert (besseli (alpha,x), ix, 100*eps)
669 %!assert (besselk (alpha,x), kx, 100*eps)
670 %!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
671 %!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
672 %!
673 %!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
674 %!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
675 %!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
676 %!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
677 %!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
678 %!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
679 %!
680 %! ## Bessel functions, fractional order, positive and negative x
681 %!
682 %! alpha = 3.5; x = 2.75;
683 %! jx = 0.1691636439842384154644784389;
684 %! yx = -0.8301381935499356070267953387;
685 %! ix = 0.3930540878794826310979363668;
686 %! kx = 0.2844099013460621170288192503;
687 %!
688 %!assert (besselj (alpha,x), jx, 100*eps)
689 %!assert (bessely (alpha,x), yx, 100*eps)
690 %!assert (besseli (alpha,x), ix, 100*eps)
691 %!assert (besselk (alpha,x), kx, 100*eps)
692 %!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
693 %!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
694 %!
695 %!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
696 %!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
697 %!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
698 %!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
699 %!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
700 %!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
701 %!
702 %! nix = 0.2119931212254662995364461998;
703 %!
704 %!assert (besselj (-alpha,x), yx, 100*eps)
705 %!assert (bessely (-alpha,x), -jx, 100*eps)
706 %!assert (besseli (-alpha,x), nix, 100*eps)
707 %!assert (besselk (-alpha,x), kx, 100*eps)
708 %!assert (besselh (-alpha,1,x), -I*(jx + I*yx), 100*eps)
709 %!assert (besselh (-alpha,2,x), I*(jx - I*yx), 100*eps)
710 %!
711 %!assert (besselj (-alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
712 %!assert (bessely (-alpha,x,1), -jx*exp(-abs(imag(x))), 100*eps)
713 %!assert (besseli (-alpha,x,1), nix*exp(-abs(real(x))), 100*eps)
714 %!assert (besselk (-alpha,x,1), kx*exp(x), 100*eps)
715 %!assert (besselh (-alpha,1,x,1), -I*(jx + I*yx)*exp(-I*x), 100*eps)
716 %!assert (besselh (-alpha,2,x,1), I*(jx - I*yx)*exp(I*x), 100*eps)
717 %!
718 %! x *= -1;
719 %! jx *= -I;
720 %! yx = -0.8301381935499356070267953387*I;
721 %! ix *= -I;
722 %! kx = -0.9504059335995575096509874508*I;
723 %!
724 %!assert (besselj (alpha,x), jx, 100*eps)
725 %!assert (bessely (alpha,x), yx, 100*eps)
726 %!assert (besseli (alpha,x), ix, 100*eps)
727 %!assert (besselk (alpha,x), kx, 100*eps)
728 %!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
729 %!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
730 %!
731 %!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
732 %!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
733 %!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
734 %!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
735 %!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
736 %!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
737 %!
738 %! ## Bessel functions, even order, complex x
739 %!
740 %! alpha = 2; x = 1.25 + 3.625 * I;
741 %! jx = -1.299533366810794494030065917 + 4.370833116012278943267479589*I;
742 %! yx = -4.370357232383223896393056727 - 1.283083391453582032688834041*I;
743 %! ix = -0.6717801680341515541002273932 - 0.2314623443930774099910228553*I;
744 %! kx = -0.01108009888623253515463783379 + 0.2245218229358191588208084197*I;
745 %!
746 %!assert (besselj (alpha,x), jx, 100*eps)
747 %!assert (bessely (alpha,x), yx, 100*eps)
748 %!assert (besseli (alpha,x), ix, 100*eps)
749 %!assert (besselk (alpha,x), kx, 100*eps)
750 %!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
751 %!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
752 %!
753 %!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
754 %!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
755 %!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
756 %!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
757 %!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
758 %!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
759 %!
760 %!assert (besselj (-alpha,x), jx, 100*eps)
761 %!assert (bessely (-alpha,x), yx, 100*eps)
762 %!assert (besseli (-alpha,x), ix, 100*eps)
763 %!assert (besselk (-alpha,x), kx, 100*eps)
764 %!assert (besselh (-alpha,1,x), jx + I*yx, 100*eps)
765 %!assert (besselh (-alpha,2,x), jx - I*yx, 100*eps)
766 %!
767 %!assert (besselj (-alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
768 %!assert (bessely (-alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
769 %!assert (besseli (-alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
770 %!assert (besselk (-alpha,x,1), kx*exp(x), 100*eps)
771 %!assert (besselh (-alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
772 %!assert (besselh (-alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
773 %!
774 %! ## Bessel functions, odd order, complex x
775 %!
776 %! alpha = 3; x = 2.5 + 1.875 * I;
777 %! jx = 0.1330721523048277493333458596 + 0.5386295217249660078754395597*I;
778 %! yx = -0.6485072392105829901122401551 + 0.2608129289785456797046996987*I;
779 %! ix = -0.6182064685486998097516365709 + 0.4677561094683470065767989920*I;
780 %! kx = -0.1568585587733540007867882337 - 0.05185853709490846050505141321*I;
781 %!
782 %!assert (besselj (alpha,x), jx, 100*eps)
783 %!assert (bessely (alpha,x), yx, 100*eps)
784 %!assert (besseli (alpha,x), ix, 100*eps)
785 %!assert (besselk (alpha,x), kx, 100*eps)
786 %!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
787 %!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
788 %!
789 %!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
790 %!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
791 %!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
792 %!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
793 %!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
794 %!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
795 %!
796 %!assert (besselj (-alpha,x), -jx, 100*eps)
797 %!assert (bessely (-alpha,x), -yx, 100*eps)
798 %!assert (besseli (-alpha,x), ix, 100*eps)
799 %!assert (besselk (-alpha,x), kx, 100*eps)
800 %!assert (besselh (-alpha,1,x), -(jx + I*yx), 100*eps)
801 %!assert (besselh (-alpha,2,x), -(jx - I*yx), 100*eps)
802 %!
803 %!assert (besselj (-alpha,x,1), -jx*exp(-abs(imag(x))), 100*eps)
804 %!assert (bessely (-alpha,x,1), -yx*exp(-abs(imag(x))), 100*eps)
805 %!assert (besseli (-alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
806 %!assert (besselk (-alpha,x,1), kx*exp(x), 100*eps)
807 %!assert (besselh (-alpha,1,x,1), -(jx + I*yx)*exp(-I*x), 100*eps)
808 %!assert (besselh (-alpha,2,x,1), -(jx - I*yx)*exp(I*x), 100*eps)
809 %!
810 %! ## Bessel functions, fractional order, complex x
811 %!
812 %! alpha = 3.5; x = 1.75 + 4.125 * I;
813 %! jx = -3.018566131370455929707009100 - 0.7585648436793900607704057611*I;
814 %! yx = 0.7772278839106298215614791107 - 3.018518722313849782683792010*I;
815 %! ix = 0.2100873577220057189038160913 - 0.6551765604618246531254970926*I;
816 %! kx = 0.1757147290513239935341488069 + 0.08772348296883849205562558311*I;
817 %!
818 %!assert (besselj (alpha,x), jx, 100*eps)
819 %!assert (bessely (alpha,x), yx, 100*eps)
820 %!assert (besseli (alpha,x), ix, 100*eps)
821 %!assert (besselk (alpha,x), kx, 100*eps)
822 %!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
823 %!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
824 %!
825 %!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
826 %!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
827 %!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
828 %!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
829 %!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
830 %!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
831 %!
832 %! nix = 0.09822388691172060573913739253 - 0.7110230642207380127317227407*I;
833 %!
834 %!assert (besselj (-alpha,x), yx, 100*eps)
835 %!assert (bessely (-alpha,x), -jx, 100*eps)
836 %!assert (besseli (-alpha,x), nix, 100*eps)
837 %!assert (besselk (-alpha,x), kx, 100*eps)
838 %!assert (besselh (-alpha,1,x), -I*(jx + I*yx), 100*eps)
839 %!assert (besselh (-alpha,2,x), I*(jx - I*yx), 100*eps)
840 %!
841 %!assert (besselj (-alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
842 %!assert (bessely (-alpha,x,1), -jx*exp(-abs(imag(x))), 100*eps)
843 %!assert (besseli (-alpha,x,1), nix*exp(-abs(real(x))), 100*eps)
844 %!assert (besselk (-alpha,x,1), kx*exp(x), 100*eps)
845 %!assert (besselh (-alpha,1,x,1), -I*(jx + I*yx)*exp(-I*x), 100*eps)
846 %!assert (besselh (-alpha,2,x,1), I*(jx - I*yx)*exp(I*x), 100*eps)
847 
848 Tests contributed by Robert T. Short.
849 Tests are based on the properties and tables in A&S:
850  Abramowitz and Stegun, "Handbook of Mathematical Functions",
851  1972.
852 
853 For regular Bessel functions, there are 3 tests. These compare octave
854 results against Tables 9.1, 9.2, and 9.4 in A&S. Tables 9.1 and 9.2
855 are good to only a few decimal places, so any failures should be
856 considered a broken implementation. Table 9.4 is an extended table
857 for larger orders and arguments. There are some differences between
858 Octave and Table 9.4, mostly in the last decimal place but in a very
859 few instances the errors are in the last two places. The comparison
860 tolerance has been changed to reflect this.
861 
862 Similarly for modified Bessel functions, there are 3 tests. These
863 compare octave results against Tables 9.8, 9.9, and 9.11 in A&S.
864 Tables 9.8 and 9.9 are good to only a few decimal places, so any
865 failures should be considered a broken implementation. Table 9.11 is
866 an extended table for larger orders and arguments. There are some
867 differences between octave and Table 9.11, mostly in the last decimal
868 place but in a very few instances the errors are in the last two
869 places. The comparison tolerance has been changed to reflect this.
870 
871 For spherical Bessel functions, there are also three tests, comparing
872 octave results to Tables 10.1, 10.2, and 10.4 in A&S. Very similar
873 comments may be made here as in the previous lines. At this time,
874 modified spherical Bessel function tests are not included.
875 
876 % Table 9.1 - J and Y for integer orders 0, 1, 2.
877 % Compare against excerpts of Table 9.1, Abramowitz and Stegun.
878 %!test
879 %! n = 0:2;
880 %! z = (0:2.5:17.5)';
881 %!
882 %! Jt = [[ 1.000000000000000, 0.0000000000, 0.0000000000];
883 %! [-0.048383776468198, 0.4970941025, 0.4460590584];
884 %! [-0.177596771314338, -0.3275791376, 0.0465651163];
885 %! [ 0.266339657880378, 0.1352484276, -0.2302734105];
886 %! [-0.245935764451348, 0.0434727462, 0.2546303137];
887 %! [ 0.146884054700421, -0.1654838046, -0.1733614634];
888 %! [-0.014224472826781, 0.2051040386, 0.0415716780];
889 %! [-0.103110398228686, -0.1634199694, 0.0844338303]];
890 %!
891 %! Yt = [[-Inf, -Inf, -Inf ];
892 %! [ 0.4980703596, 0.1459181380, -0.38133585 ];
893 %! [-0.3085176252, 0.1478631434, 0.36766288 ];
894 %! [ 0.1173132861, -0.2591285105, -0.18641422 ];
895 %! [ 0.0556711673, 0.2490154242, -0.00586808 ];
896 %! [-0.1712143068, -0.1538382565, 0.14660019 ];
897 %! [ 0.2054642960, 0.0210736280, -0.20265448 ];
898 %! [-0.1604111925, 0.0985727987, 0.17167666 ]];
899 %!
900 %! J = besselj (n,z);
901 %! Y = bessely (n,z);
902 %! assert (Jt(:,1), J(:,1), 0.5e-10);
903 %! assert (Yt(:,1), Y(:,1), 0.5e-10);
904 %! assert (Jt(:,2:3), J(:,2:3), 0.5e-10);
905 
906 Table 9.2 - J and Y for integer orders 3-9.
907 
908 %!test
909 %! n = (3:9);
910 %! z = (0:2:20).';
911 %!
912 %! Jt = [[ 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00];
913 %! [ 1.2894e-01, 3.3996e-02, 7.0396e-03, 1.2024e-03, 1.7494e-04, 2.2180e-05, 2.4923e-06];
914 %! [ 4.3017e-01, 2.8113e-01, 1.3209e-01, 4.9088e-02, 1.5176e-02, 4.0287e-03, 9.3860e-04];
915 %! [ 1.1477e-01, 3.5764e-01, 3.6209e-01, 2.4584e-01, 1.2959e-01, 5.6532e-02, 2.1165e-02];
916 %! [-2.9113e-01,-1.0536e-01, 1.8577e-01, 3.3758e-01, 3.2059e-01, 2.2345e-01, 1.2632e-01];
917 %! [ 5.8379e-02,-2.1960e-01,-2.3406e-01,-1.4459e-02, 2.1671e-01, 3.1785e-01, 2.9186e-01];
918 %! [ 1.9514e-01, 1.8250e-01,-7.3471e-02,-2.4372e-01,-1.7025e-01, 4.5095e-02, 2.3038e-01];
919 %! [-1.7681e-01, 7.6244e-02, 2.2038e-01, 8.1168e-02,-1.5080e-01,-2.3197e-01,-1.1431e-01];
920 %! [-4.3847e-02,-2.0264e-01,-5.7473e-02, 1.6672e-01, 1.8251e-01,-7.0211e-03,-1.8953e-01];
921 %! [ 1.8632e-01, 6.9640e-02,-1.5537e-01,-1.5596e-01, 5.1399e-02, 1.9593e-01, 1.2276e-01];
922 %! [-9.8901e-02, 1.3067e-01, 1.5117e-01,-5.5086e-02,-1.8422e-01,-7.3869e-02, 1.2513e-01]];
923 %!
924 %! Yt = [[ -Inf, -Inf, -Inf, -Inf, -Inf, -Inf, -Inf];
925 %! [-1.1278e+00,-2.7659e+00,-9.9360e+00,-4.6914e+01,-2.7155e+02,-1.8539e+03,-1.4560e+04];
926 %! [-1.8202e-01,-4.8894e-01,-7.9585e-01,-1.5007e+00,-3.7062e+00,-1.1471e+01,-4.2178e+01];
927 %! [ 3.2825e-01, 9.8391e-02,-1.9706e-01,-4.2683e-01,-6.5659e-01,-1.1052e+00,-2.2907e+00];
928 %! [ 2.6542e-02, 2.8294e-01, 2.5640e-01, 3.7558e-02,-2.0006e-01,-3.8767e-01,-5.7528e-01];
929 %! [-2.5136e-01,-1.4495e-01, 1.3540e-01, 2.8035e-01, 2.0102e-01, 1.0755e-03,-1.9930e-01];
930 %! [ 1.2901e-01,-1.5122e-01,-2.2982e-01,-4.0297e-02, 1.8952e-01, 2.6140e-01, 1.5902e-01];
931 %! [ 1.2350e-01, 2.0393e-01,-6.9717e-03,-2.0891e-01,-1.7209e-01, 3.6816e-02, 2.1417e-01];
932 %! [-1.9637e-01,-7.3222e-05, 1.9633e-01, 1.2278e-01,-1.0425e-01,-2.1399e-01,-1.0975e-01];
933 %! [ 3.3724e-02,-1.7722e-01,-1.1249e-01, 1.1472e-01, 1.8897e-01, 3.2253e-02,-1.6030e-01];
934 %! [ 1.4967e-01, 1.2409e-01,-1.0004e-01,-1.7411e-01,-4.4312e-03, 1.7101e-01, 1.4124e-01]];
935 %!
936 %! n = (3:9);
937 %! z = (0:2:20).';
938 %! J = besselj (n,z);
939 %! Y = bessely (n,z);
940 %!
941 %! assert (J(1,:), zeros (1, columns (J)));
942 %! assert (J(2:end,:), Jt(2:end,:), -5e-5);
943 %! assert (Yt(1,:), Y(1,:));
944 %! assert (Y(2:end,:), Yt(2:end,:), -5e-5);
945 
946 Table 9.4 - J and Y for various integer orders and arguments.
947 
948 %!test
949 %! Jt = [[ 7.651976866e-01, 2.238907791e-01, -1.775967713e-01, -2.459357645e-01, 5.581232767e-02, 1.998585030e-02];
950 %! [ 2.497577302e-04, 7.039629756e-03, 2.611405461e-01, -2.340615282e-01, -8.140024770e-02, -7.419573696e-02];
951 %! [ 2.630615124e-10, 2.515386283e-07, 1.467802647e-03, 2.074861066e-01, -1.138478491e-01, -5.473217694e-02];
952 %! [ 2.297531532e-17, 7.183016356e-13, 4.796743278e-07, 4.507973144e-03, -1.082255990e-01, 1.519812122e-02];
953 %! [ 3.873503009e-25, 3.918972805e-19, 2.770330052e-11, 1.151336925e-05, -1.167043528e-01, 6.221745850e-02];
954 %! [ 3.482869794e-42, 3.650256266e-33, 2.671177278e-21, 1.551096078e-12, 4.843425725e-02, 8.146012958e-02];
955 %! [ 1.107915851e-60, 1.196077458e-48, 8.702241617e-33, 6.030895312e-21, -1.381762812e-01, 7.270175482e-02];
956 %! [ 2.906004948e-80, 3.224095839e-65, 2.294247616e-45, 1.784513608e-30, 1.214090219e-01, -3.869833973e-02];
957 %! [ 8.431828790e-189, 1.060953112e-158, 6.267789396e-119, 6.597316064e-89, 1.115927368e-21, 9.636667330e-02]];
958 %!
959 %! Yt = [[ 8.825696420e-02, 5.103756726e-01, -3.085176252e-01, 5.567116730e-02, -9.806499547e-02, -7.724431337e-02]
960 %! [-2.604058666e+02, -9.935989128e+00, -4.536948225e-01, 1.354030477e-01, -7.854841391e-02, -2.948019628e-02]
961 %! [-1.216180143e+08, -1.291845422e+05, -2.512911010e+01, -3.598141522e-01, 5.723897182e-03, 5.833157424e-02]
962 %! [-9.256973276e+14, -2.981023646e+10, -4.694049564e+04, -6.364745877e+00, 4.041280205e-02, 7.879068695e-02]
963 %! [-4.113970315e+22, -4.081651389e+16, -5.933965297e+08, -1.597483848e+03, 1.644263395e-02, 5.124797308e-02]
964 %! [-3.048128783e+39, -2.913223848e+30, -4.028568418e+18, -7.256142316e+09, -1.164572349e-01, 6.138839212e-03]
965 %! [-7.184874797e+57, -6.661541235e+45, -9.216816571e+29, -1.362803297e+18, -4.530801120e-02, 4.074685217e-02]
966 %! [-2.191142813e+77, -1.976150576e+62, -2.788837017e+42, -3.641066502e+27, -2.103165546e-01, 7.650526394e-02]
967 %! [-3.775287810e+185, -3.000826049e+155, -5.084863915e+115, -4.849148271e+85, -3.293800188e+18, -1.669214114e-01]];
968 %!
969 %! n = [(0:5:20).';30;40;50;100];
970 %! z = [1,2,5,10,50,100];
971 %! J = besselj (n.', z.').';
972 %! Y = bessely (n.', z.').';
973 %! assert (J, Jt, -1e-9);
974 %! assert (Y, Yt, -1e-9);
975 
976 Table 9.8 - I and K for integer orders 0, 1, 2.
977 
978 %!test
979 %! n = 0:2;
980 %! z1 = [0.1;2.5;5.0];
981 %! z2 = [7.5;10.0;15.0;20.0];
982 %! rtbl = [[ 0.9071009258 0.0452984468 0.1251041992 2.6823261023 10.890182683 1.995039646 ];
983 %! [ 0.2700464416 0.2065846495 0.2042345837 0.7595486903 0.9001744239 0.759126289 ];
984 %! [ 0.1835408126 0.1639722669 0.7002245988 0.5478075643 0.6002738588 0.132723593 ];
985 %! [ 0.1483158301 0.1380412115 0.111504840 0.4505236991 0.4796689336 0.57843541 ];
986 %! [ 0.1278333372 0.1212626814 0.103580801 0.3916319344 0.4107665704 0.47378525 ];
987 %! [ 0.1038995314 0.1003741751 0.090516308 0.3210023535 0.3315348950 0.36520701 ];
988 %! [ 0.0897803119 0.0875062222 0.081029690 0.2785448768 0.2854254970 0.30708743 ]];
989 %!
990 %! tbl = [besseli(n,z1,1), besselk(n,z1,1)];
991 %! tbl(:,3) = tbl(:,3) .* (exp (z1) .* z1.^(-2));
992 %! tbl(:,6) = tbl(:,6) .* (exp (-z1) .* z1.^(2));
993 %! tbl = [tbl;[besseli(n,z2,1),besselk(n,z2,1)]];
994 %!
995 %! assert (tbl, rtbl, -2e-8);
996 
997 Table 9.9 - I and K for orders 3-9.
998 
999 %!test
1000 %! It = [[ 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00];
1001 %! [ 2.8791e-02 6.8654e-03 1.3298e-03 2.1656e-04 3.0402e-05 3.7487e-06 4.1199e-07];
1002 %! [ 6.1124e-02 2.5940e-02 9.2443e-03 2.8291e-03 7.5698e-04 1.7968e-04 3.8284e-05];
1003 %! [ 7.4736e-02 4.1238e-02 1.9752e-02 8.3181e-03 3.1156e-03 1.0484e-03 3.1978e-04];
1004 %! [ 7.9194e-02 5.0500e-02 2.8694e-02 1.4633e-02 6.7449e-03 2.8292e-03 1.0866e-03];
1005 %! [ 7.9830e-02 5.5683e-02 3.5284e-02 2.0398e-02 1.0806e-02 5.2694e-03 2.3753e-03];
1006 %! [ 7.8848e-02 5.8425e-02 3.9898e-02 2.5176e-02 1.4722e-02 8.0010e-03 4.0537e-03];
1007 %! [ 7.7183e-02 5.9723e-02 4.3056e-02 2.8969e-02 1.8225e-02 1.0744e-02 5.9469e-03];
1008 %! [ 7.5256e-02 6.0155e-02 4.5179e-02 3.1918e-02 2.1240e-02 1.3333e-02 7.9071e-03];
1009 %! [ 7.3263e-02 6.0059e-02 4.6571e-02 3.4186e-02 2.3780e-02 1.5691e-02 9.8324e-03];
1010 %! [ 7.1300e-02 5.9640e-02 4.7444e-02 3.5917e-02 2.5894e-02 1.7792e-02 1.1661e-02]];
1011 %!
1012 %! Kt = [[ Inf Inf Inf Inf Inf Inf Inf];
1013 %! [ 4.7836e+00 1.6226e+01 6.9687e+01 3.6466e+02 2.2576e+03 1.6168e+04 1.3160e+05];
1014 %! [ 1.6317e+00 3.3976e+00 8.4268e+00 2.4465e+01 8.1821e+01 3.1084e+02 1.3252e+03];
1015 %! [ 9.9723e-01 1.6798e+00 3.2370e+00 7.0748e+00 1.7387e+01 4.7644e+01 1.4444e+02];
1016 %! [ 7.3935e-01 1.1069e+00 1.8463e+00 3.4148e+00 6.9684e+00 1.5610e+01 3.8188e+01];
1017 %! [ 6.0028e-01 8.3395e-01 1.2674e+00 2.1014e+00 3.7891e+00 7.4062e+00 1.5639e+01];
1018 %! [ 5.1294e-01 6.7680e-01 9.6415e-01 1.4803e+00 2.4444e+00 4.3321e+00 8.2205e+00];
1019 %! [ 4.5266e-01 5.7519e-01 7.8133e-01 1.1333e+00 1.7527e+00 2.8860e+00 5.0510e+00];
1020 %! [ 4.0829e-01 5.0414e-01 6.6036e-01 9.1686e-01 1.3480e+00 2.0964e+00 3.4444e+00];
1021 %! [ 3.7411e-01 4.5162e-01 5.7483e-01 7.7097e-01 1.0888e+00 1.6178e+00 2.5269e+00];
1022 %! [ 3.4684e-01 4.1114e-01 5.1130e-01 6.6679e-01 9.1137e-01 1.3048e+00 1.9552e+00]];
1023 %!
1024 %! n = (3:9);
1025 %! z = (0:2:20).';
1026 %! I = besseli (n,z,1);
1027 %! K = besselk (n,z,1);
1028 %!
1029 %! assert (abs (I(1,:)), zeros (1, columns (I)));
1030 %! assert (I(2:end,:), It(2:end,:), -5e-5);
1031 %! assert (Kt(1,:), K(1,:));
1032 %! assert (K(2:end,:), Kt(2:end,:), -5e-5);
1033 
1034 Table 9.11 - I and K for various integer orders and arguments.
1035 
1036 %!test
1037 %! It = [[ 1.266065878e+00 2.279585302e+00 2.723987182e+01 2.815716628e+03 2.93255378e+20 1.07375171e+42 ];
1038 %! [ 2.714631560e-04 9.825679323e-03 2.157974547e+00 7.771882864e+02 2.27854831e+20 9.47009387e+41 ];
1039 %! [ 2.752948040e-10 3.016963879e-07 4.580044419e-03 2.189170616e+01 1.07159716e+20 6.49897552e+41 ];
1040 %! [ 2.370463051e-17 8.139432531e-13 1.047977675e-06 1.043714907e-01 3.07376455e+19 3.47368638e+41 ];
1041 %! [ 3.966835986e-25 4.310560576e-19 5.024239358e-11 1.250799736e-04 5.44200840e+18 1.44834613e+41 ];
1042 %! [ 3.539500588e-42 3.893519664e-33 3.997844971e-21 7.787569783e-12 4.27499365e+16 1.20615487e+40 ];
1043 %! [ 1.121509741e-60 1.255869192e-48 1.180426980e-32 2.042123274e-20 6.00717897e+13 3.84170550e+38 ];
1044 %! [ 2.934635309e-80 3.353042830e-65 2.931469647e-45 4.756894561e-30 1.76508024e+10 4.82195809e+36 ];
1045 %! [ 8.473674008e-189 1.082171475e-158 7.093551489e-119 1.082344202e-88 2.72788795e-16 4.64153494e+21 ]];
1046 %!
1047 %! Kt = [[ 4.210244382e-01 1.138938727e-01 3.691098334e-03 1.778006232e-05 3.41016774e-23 4.65662823e-45 ];
1048 %! [ 3.609605896e+02 9.431049101e+00 3.270627371e-02 5.754184999e-05 4.36718224e-23 5.27325611e-45 ];
1049 %! [ 1.807132899e+08 1.624824040e+05 9.758562829e+00 1.614255300e-03 9.15098819e-23 7.65542797e-45 ];
1050 %! [ 1.403066801e+15 4.059213332e+10 3.016976630e+04 2.656563849e-01 3.11621117e-22 1.42348325e-44 ];
1051 %! [ 6.294369360e+22 5.770856853e+16 4.827000521e+08 1.787442782e+02 1.70614838e-21 3.38520541e-44 ];
1052 %! [ 4.706145527e+39 4.271125755e+30 4.112132063e+18 2.030247813e+09 2.00581681e-19 3.97060205e-43 ];
1053 %! [ 1.114220651e+58 9.940839886e+45 1.050756722e+30 5.938224681e+17 1.29986971e-16 1.20842080e-41 ];
1054 %! [ 3.406896854e+77 2.979981740e+62 3.394322243e+42 2.061373775e+27 4.00601349e-13 9.27452265e-40 ];
1055 %! [ 5.900333184e+185 4.619415978e+155 7.039860193e+115 4.596674084e+85 1.63940352e+13 7.61712963e-25 ]];
1056 %!
1057 %! n = [(0:5:20).';30;40;50;100];
1058 %! z = [1,2,5,10,50,100];
1059 %! I = besseli (n.', z.').';
1060 %! K = besselk (n.', z.').';
1061 %! assert (I, It, -5e-9);
1062 %! assert (K, Kt, -5e-9);
1063 
1064 The next section checks that negative integer orders and positive
1065 integer orders are appropriately related.
1066 
1067 %!test
1068 %! n = (0:2:20);
1069 %! assert (besselj (n,1), besselj (-n,1), 1e-8);
1070 %! assert (-besselj (n+1,1), besselj (-n-1,1), 1e-8);
1071 
1072 besseli (n,z) = besseli (-n,z);
1073 
1074 %!test
1075 %! n = (0:2:20);
1076 %! assert (besseli (n,1), besseli (-n,1), 1e-8);
1077 
1078 Table 10.1 - j and y for integer orders 0, 1, 2.
1079 Compare against excerpts of Table 10.1, Abramowitz and Stegun.
1080 
1081 %!test
1082 %! n = (0:2);
1083 %! z = [0.1;(2.5:2.5:10.0).'];
1084 %!
1085 %! jt = [[ 9.9833417e-01 3.33000119e-02 6.6619061e-04 ];
1086 %! [ 2.3938886e-01 4.16212989e-01 2.6006673e-01 ];
1087 %! [-1.9178485e-01 -9.50894081e-02 1.3473121e-01 ];
1088 %! [ 1.2507e-01 -2.9542e-02 -1.3688e-01 ];
1089 %! [ -5.4402e-02 7.8467e-02 7.7942e-02 ]];
1090 %!
1091 %! yt = [[-9.9500417e+00 -1.0049875e+02 -3.0050125e+03 ];
1092 %! [ 3.2045745e-01 -1.1120588e-01 -4.5390450e-01 ];
1093 %! [-5.6732437e-02 1.8043837e-01 1.6499546e-01 ];
1094 %! [ -4.6218e-02 -1.3123e-01 -6.2736e-03 ];
1095 %! [ 8.3907e-02 6.2793e-02 -6.5069e-02 ]];
1096 %!
1097 %! j = sqrt ((pi/2)./z) .* besselj (n+1/2,z);
1098 %! y = sqrt ((pi/2)./z) .* bessely (n+1/2,z);
1099 %! assert (jt, j, -5e-5);
1100 %! assert (yt, y, -5e-5);
1101 
1102 Table 10.2 - j and y for orders 3-8.
1103 Compare against excerpts of Table 10.2, Abramowitzh and Stegun.
1104 
1105  Important note: In A&S, y_4(0.1) = -1.0507e+7, but Octave returns
1106  y_4(0.1) = -1.0508e+07 (-10507503.75). If I compute the same term using
1107  a series, the difference is in the eighth significant digit so I left
1108  the Octave results in place.
1109 
1110 %!test
1111 %! n = (3:8);
1112 %! z = (0:2.5:10).'; z(1) = 0.1;
1113 %!
1114 %! jt = [[ 9.5185e-06 1.0577e-07 9.6163e-10 7.3975e-12 4.9319e-14 2.9012e-16];
1115 %! [ 1.0392e-01 3.0911e-02 7.3576e-03 1.4630e-03 2.5009e-04 3.7516e-05];
1116 %! [ 2.2982e-01 1.8702e-01 1.0681e-01 4.7967e-02 1.7903e-02 5.7414e-03];
1117 %! [-6.1713e-02 7.9285e-02 1.5685e-01 1.5077e-01 1.0448e-01 5.8188e-02];
1118 %! [-3.9496e-02 -1.0559e-01 -5.5535e-02 4.4501e-02 1.1339e-01 1.2558e-01]];
1119 %!
1120 %! yt = [[-1.5015e+05 -1.0508e+07 -9.4553e+08 -1.0400e+11 -1.3519e+13 -2.0277e+15];
1121 %! [-7.9660e-01 -1.7766e+00 -5.5991e+00 -2.2859e+01 -1.1327e+02 -6.5676e+02];
1122 %! [-1.5443e-02 -1.8662e-01 -3.2047e-01 -5.1841e-01 -1.0274e+00 -2.5638e+00];
1123 %! [ 1.2705e-01 1.2485e-01 2.2774e-02 -9.1449e-02 -1.8129e-01 -2.7112e-01];
1124 %! [-9.5327e-02 -1.6599e-03 9.3834e-02 1.0488e-01 4.2506e-02 -4.1117e-02]];
1125 %!
1126 %! j = sqrt ((pi/2)./z) .* besselj (n+1/2,z);
1127 %! y = sqrt ((pi/2)./z) .* bessely (n+1/2,z);
1128 %!
1129 %! assert (jt, j, -5e-5);
1130 %! assert (yt, y, -5e-5);
1131 
1132 Table 10.4 - j and y for various integer orders and arguments.
1133 
1134 %!test
1135 %! jt = [[ 8.414709848e-01 4.546487134e-01 -1.917848549e-01 -5.440211109e-02 -5.247497074e-03 -5.063656411e-03];
1136 %! [ 9.256115861e-05 2.635169770e-03 1.068111615e-01 -5.553451162e-02 -2.004830056e-02 -9.290148935e-03];
1137 %! [ 7.116552640e-11 6.825300865e-08 4.073442442e-04 6.460515449e-02 -1.503922146e-02 -1.956578597e-04];
1138 %! [ 5.132686115e-18 1.606982166e-13 1.084280182e-07 1.063542715e-03 -1.129084539e-02 7.877261748e-03];
1139 %! [ 7.537795722e-26 7.632641101e-20 5.427726761e-12 2.308371961e-06 -1.578502990e-02 1.010767128e-02];
1140 %! [ 5.566831267e-43 5.836617888e-34 4.282730217e-22 2.512057385e-13 -1.494673454e-03 8.700628514e-03];
1141 %! [ 1.538210374e-61 1.660978779e-49 1.210347583e-33 8.435671634e-22 -2.606336952e-02 1.043410851e-02];
1142 %! [ 3.615274717e-81 4.011575290e-66 2.857479350e-46 2.230696023e-31 1.882910737e-02 5.797140882e-04];
1143 %! [7.444727742e-190 9.367832591e-160 5.535650303e-120 5.832040182e-90 1.019012263e-22 1.088047701e-02]];
1144 %!
1145 %! yt = [[ -5.403023059e-01 2.080734183e-01 -5.673243709e-02 8.390715291e-02 -1.929932057e-02 -8.623188723e-03]
1146 %! [ -9.994403434e+02 -1.859144531e+01 -3.204650467e-01 9.383354168e-02 -6.971131965e-04 3.720678486e-03]
1147 %! [ -6.722150083e+08 -3.554147201e+05 -2.665611441e+01 -1.724536721e-01 1.352468751e-02 1.002577737e-02]
1148 %! [ -6.298007233e+15 -1.012182944e+11 -6.288146513e+04 -3.992071745e+00 1.712319725e-02 6.258641510e-03]
1149 %! [ -3.239592219e+23 -1.605436493e+17 -9.267951403e+08 -1.211210605e+03 1.375953130e-02 5.631729379e-05]
1150 %! [ -2.946428547e+40 -1.407393871e+31 -7.760717570e+18 -6.908318646e+09 -2.241226812e-02 -5.412929349e-03]
1151 %! [ -8.028450851e+58 -3.720929322e+46 -2.055758716e+30 -1.510304919e+18 4.978797221e-05 -7.048420407e-04]
1152 %! [ -2.739192285e+78 -1.235021944e+63 -6.964109188e+42 -4.528227272e+27 -4.190000150e-02 1.074782297e-02]
1153 %! [-6.683079463e+186 -2.655955830e+156 -1.799713983e+116 -8.573226309e+85 -1.125692891e+18 -2.298385049e-02]];
1154 %!
1155 %! n = [(0:5:20).';30;40;50;100];
1156 %! z = [1,2,5,10,50,100];
1157 %! j = sqrt ((pi/2)./z) .* besselj ((n+1/2).', z.').';
1158 %! y = sqrt ((pi/2)./z) .* bessely ((n+1/2).', z.').';
1159 %! assert (j, jt, -1e-9);
1160 %! assert (y, yt, -1e-9);
1161 */
ierr
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE const F77_INT F77_INT * ierr
Definition: lo-slatec-proto.h:109
octave_value
Definition: ov.h:70
octave_value::xfloat_complex_column_vector_value
FloatComplexColumnVector xfloat_complex_column_vector_value(const char *fmt,...) const
Definition: ov.cc:2140
print_usage
OCTINTERP_API void print_usage(void)
Definition: defun.cc:52
octave_value_list::length
octave_idx_type length(void) const
Definition: ovl.h:96
DO_BESSEL
#define DO_BESSEL(type, alpha, x, scaled, ierr, result)
Definition: besselj.cc:45
octave::math::besselk
Complex besselk(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
Definition: lo-specfun.cc:1337
octave_value::is_scalar_type
bool is_scalar_type(void) const
Definition: ov.h:673
octave_value::is_numeric_type
bool is_numeric_type(void) const
Definition: ov.h:679
DEFUN
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:46
error
void error(const char *fmt,...)
Definition: error.cc:570
octave_value::xcomplex_array_value
ComplexNDArray xcomplex_array_value(const char *fmt,...) const
Definition: ov.cc:2081
octave_value::xcomplex_value
Complex xcomplex_value(const char *fmt,...) const
Definition: ov.cc:2075
octave_value_list
Definition: ovl.h:39
dim_vector::numel
octave_idx_type numel(int n=0) const
Number of elements that a matrix with this dimensions would have.
Definition: dim-vector.h:389
args
JNIEnv void * args
Definition: ov-java.cc:67
octave_value::xfloat_complex_array_value
FloatComplexNDArray xfloat_complex_array_value(const char *fmt,...) const
Definition: ov.cc:2082
nargout
OCTAVE_EXPORT octave_value_list return the number of command line arguments passed to Octave If called with the optional argument the function xample nargout(@histc)
Definition: ov-usr-fcn.cc:935
octave::math::besseli
Complex besseli(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
Definition: lo-specfun.cc:1336
octave_value::is_bool_type
bool is_bool_type(void) const
Definition: ov.h:661
do_bessel
octave_value_list do_bessel(enum bessel_type type, const char *fn, const octave_value_list &args, int nargout)
Definition: besselj.cc:81
nargin
int nargin
Definition: graphics.cc:10115
retval
octave_value retval
Definition: data.cc:6294
type
idx type
Definition: ov.cc:3129
result
With real return the complex result
Definition: data.cc:3375
octave_value::bool_value
bool bool_value(bool warn=false) const
Definition: ov.h:819
octave_value::xfloat_complex_value
FloatComplex xfloat_complex_value(const char *fmt,...) const
Definition: ov.cc:2076
scale
void scale(Matrix &m, double x, double y, double z)
Definition: graphics.cc:5150
octave_idx_type
FloatComplex
std::complex< float > FloatComplex
Definition: oct-cmplx.h:32
octave::math::airy
Complex airy(const Complex &z, bool deriv, bool scaled, octave_idx_type &ierr)
Definition: lo-specfun.cc:1971
Complex
std::complex< double > Complex
Definition: oct-cmplx.h:31
bessel_type
bessel_type
Definition: besselj.cc:35
octave_value::is_single_type
bool is_single_type(void) const
Definition: ov.h:627
octave_value::xcomplex_column_vector_value
ComplexColumnVector xcomplex_column_vector_value(const char *fmt,...) const
Definition: ov.cc:2134
octave_value::double_value
double double_value(bool frc_str_conv=false) const
Definition: ov.h:775
dim_vector
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
x
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE * x
Definition: lo-slatec-proto.h:109
octave::math::besselj
Complex besselj(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
Definition: lo-specfun.cc:1334
octave::math::bessely
Complex bessely(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
Definition: lo-specfun.cc:1335
octave::math::biry
Complex biry(const Complex &z, bool deriv, bool scaled, octave_idx_type &ierr)
Definition: lo-specfun.cc:2005
Generated on Wed May 10 2017 15:42:44 for GNU Octave by   doxygen 1.8.8