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