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