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