1# bessel.tcl --
2#    Evaluate the most common Bessel functions
3#
4# TODO:
5#    Yn - finding decent approximations seems tough
6#    Jnu - for arbitrary values of the parameter
7#    J'n - first derivative (from recurrence relation)
8#    Kn - forward application of recurrence relation?
9#
10
11# namespace special
12#    Create a convenient namespace for the "special" mathematical functions
13#
14namespace eval ::math::special {
15    #
16    # Define a number of common mathematical constants
17    #
18    ::math::constants::constants pi
19
20    #
21    # Export the functions
22    #
23    namespace export J0 J1 Jn J1/2 J-1/2 I_n
24}
25
26# J0 --
27#    Zeroth-order Bessel function
28#
29# Arguments:
30#    x         Value of the x-coordinate
31# Result:
32#    Value of J0(x)
33#
34proc ::math::special::J0 {x} {
35    Jn 0 $x
36}
37
38# J1 --
39#    First-order Bessel function
40#
41# Arguments:
42#    x         Value of the x-coordinate
43# Result:
44#    Value of J1(x)
45#
46proc ::math::special::J1 {x} {
47    Jn 1 $x
48}
49
50# Jn --
51#    Compute the Bessel function of the first kind of order n
52# Arguments:
53#    n       Order of the function (must be integer)
54#    x       Value of the argument
55# Result:
56#    Jn(x)
57# Note:
58#    This relies on the integral representation for
59#    the Bessel functions of integer order:
60#             1     I pi
61#    Jn(x) = --     I   cos(x sin t - nt) dt
62#            pi   0 I
63#
64#    For this kind of integrands the trapezoidal rule is
65#    very efficient according to Davis and Rabinowitz
66#    (Methods of numerical integration, 1984).
67#
68proc ::math::special::Jn {n x} {
69    variable pi
70
71    if { ![string is integer -strict $n] } {
72         return -code error "Order argument must be integer"
73    }
74
75    #
76    # Integrate over the interval [0,pi] using a small
77    # enough step - 40 points should do a good job
78    # with |x| < 20, n < 20 (an accuracy of 1.0e-8
79    # is reported by Davis and Rabinowitz)
80    #
81    set number 40
82    set step   [expr {$pi/double($number)}]
83    set result 0.0
84
85    for { set i 0 } { $i <= $number } { incr i } {
86        set t [expr {double($i)*$step}]
87        set f [expr {cos($x * sin($t) - $n * $t)}]
88        if { $i == 0 || $i == $number } {
89            set f [expr {$f/2.0}]
90        }
91        set result [expr {$result+$f}]
92    }
93
94    expr {$result*$step/$pi}
95}
96
97# J1/2 --
98#    Half-order Bessel function
99#
100# Arguments:
101#    x         Value of the x-coordinate
102# Result:
103#    Value of J1/2(x)
104#
105proc ::math::special::J1/2 {x} {
106    variable pi
107    #
108    # This Bessel function can be expressed in terms of elementary
109    # functions. Therefore use the explicit formula
110    #
111    if { $x != 0.0 } {
112        expr {sqrt(2.0/$pi/$x)*sin($x)}
113    } else {
114        return 0.0
115    }
116}
117
118# J-1/2 --
119#    Compute the Bessel function of the first kind of order -1/2
120# Arguments:
121#    x       Value of the argument (!= 0.0)
122# Result:
123#    J-1/2(x)
124#
125proc ::math::special::J-1/2 {x} {
126    variable pi
127    if { $x == 0.0 } {
128        return -code error "Argument must not be zero (singularity)"
129    } else {
130        return [expr {-cos($x)/sqrt($pi*$x/2.0)}]
131    }
132}
133
134# I_n --
135#    Compute the modified Bessel function of the first kind
136#
137# Arguments:
138#    n            Order of the function (must be positive integer or zero)
139#    x            Abscissa at which to compute it
140# Result:
141#    Value of In(x)
142# Note:
143#    This relies on Miller's algorithm for finding minimal solutions
144#
145namespace eval ::math::special {}
146
147proc ::math::special::I_n {n x} {
148    if { ! [string is integer $n] || $n < 0 } {
149        error "Wrong order: must be positive integer or zero"
150    }
151
152    set n2 [expr {$n+8}]  ;# Note: just a guess that this will be enough
153
154    set ynp1 0.0
155    set yn   1.0
156    set sum  1.0
157
158    while { $n2 > 0 } {
159        set ynm1 [expr {$ynp1+2.0*$n2*$yn/$x}]
160        set sum  [expr {$sum+$ynm1}]
161        if { $n2 == $n+1 } {
162           set result $ynm1
163        }
164        set ynp1 $yn
165        set yn   $ynm1
166        incr n2 -1
167    }
168
169    set quotient [expr {(2.0*$sum-$ynm1)/exp($x)}]
170
171    expr {$result/$quotient}
172}
173
174#
175# some tests --
176#
177if { 0 } {
178set tcl_precision 17
179foreach x {0.0 2.0 4.4 6.0 10.0 11.0 12.0 13.0 14.0} {
180    puts "J0($x) = [::math::special::J0 $x] - J1($x) = [::math::special::J1 $x] \
181- J1/2($x) = [::math::special::J1/2 $x]"
182}
183foreach n {0 1 2 3 4 5} {
184    puts [::math::special::I_n $n 1.0]
185}
186}
187