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