You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
200 lines
4.5 KiB
200 lines
4.5 KiB
# bessel.tcl -- |
|
# Evaluate the most common Bessel functions |
|
# |
|
# TODO: |
|
# Yn - finding decent approximations seems tough |
|
# Jnu - for arbitrary values of the parameter |
|
# J'n - first derivative (from recurrence relation) |
|
# Kn - forward application of recurrence relation? |
|
# |
|
|
|
# namespace special |
|
# Create a convenient namespace for the "special" mathematical functions |
|
# |
|
namespace eval ::math::special { |
|
# |
|
# Define a number of common mathematical constants |
|
# |
|
::math::constants::constants pi |
|
|
|
# |
|
# Export the functions |
|
# |
|
namespace export J0 J1 Jn J1/2 J-1/2 I_n |
|
} |
|
|
|
# J0 -- |
|
# Zeroth-order Bessel function |
|
# |
|
# Arguments: |
|
# x Value of the x-coordinate |
|
# Result: |
|
# Value of J0(x) |
|
# |
|
proc ::math::special::J0 {x} { |
|
Jn 0 $x |
|
} |
|
|
|
# J1 -- |
|
# First-order Bessel function |
|
# |
|
# Arguments: |
|
# x Value of the x-coordinate |
|
# Result: |
|
# Value of J1(x) |
|
# |
|
proc ::math::special::J1 {x} { |
|
Jn 1 $x |
|
} |
|
|
|
# Jn -- |
|
# Compute the Bessel function of the first kind of order n |
|
# Arguments: |
|
# n Order of the function (must be integer) |
|
# x Value of the argument |
|
# Result: |
|
# Jn(x) |
|
# Note: |
|
# This relies on the integral representation for |
|
# the Bessel functions of integer order: |
|
# 1 I pi |
|
# Jn(x) = -- I cos(x sin t - nt) dt |
|
# pi 0 I |
|
# |
|
# For this kind of integrands the trapezoidal rule is |
|
# very efficient according to Davis and Rabinowitz |
|
# (Methods of numerical integration, 1984). |
|
# |
|
proc ::math::special::Jn {n x} { |
|
variable pi |
|
|
|
##nagelfar ignore |
|
if { ![string is integer -strict $n] } { |
|
return -code error "Order argument must be integer" |
|
} |
|
|
|
# |
|
# Integrate over the interval [0,pi] using a small |
|
# enough step - 40 points should do a good job |
|
# with |x| < 20, n < 20 (an accuracy of 1.0e-8 |
|
# is reported by Davis and Rabinowitz) |
|
# |
|
set number 40 |
|
set step [expr {$pi/double($number)}] |
|
set result 0.0 |
|
|
|
for { set i 0 } { $i <= $number } { incr i } { |
|
set t [expr {double($i)*$step}] |
|
set f [expr {cos($x * sin($t) - $n * $t)}] |
|
if { $i == 0 || $i == $number } { |
|
set f [expr {$f/2.0}] |
|
} |
|
set result [expr {$result+$f}] |
|
} |
|
|
|
expr {$result*$step/$pi} |
|
} |
|
|
|
# J1/2 -- |
|
# Half-order Bessel function |
|
# |
|
# Arguments: |
|
# x Value of the x-coordinate |
|
# Result: |
|
# Value of J1/2(x) |
|
# |
|
proc ::math::special::J1/2 {x} { |
|
variable pi |
|
# |
|
# This Bessel function can be expressed in terms of elementary |
|
# functions. Therefore use the explicit formula |
|
# |
|
if { $x != 0.0 } { |
|
expr {sqrt(2.0/$pi/$x)*sin($x)} |
|
} else { |
|
return 0.0 |
|
} |
|
} |
|
|
|
# J-1/2 -- |
|
# Compute the Bessel function of the first kind of order -1/2 |
|
# Arguments: |
|
# x Value of the argument (!= 0.0) |
|
# Result: |
|
# J-1/2(x) |
|
# |
|
proc ::math::special::J-1/2 {x} { |
|
variable pi |
|
if { $x == 0.0 } { |
|
return -code error "Argument must not be zero (singularity)" |
|
} else { |
|
return [expr {-cos($x)/sqrt($pi*$x/2.0)}] |
|
} |
|
} |
|
|
|
# I_n -- |
|
# Compute the modified Bessel function of the first kind |
|
# |
|
# Arguments: |
|
# n Order of the function (must be positive integer or zero) |
|
# x Abscissa at which to compute it |
|
# Result: |
|
# Value of In(x) |
|
# Note: |
|
# This relies on Miller's algorithm for finding minimal solutions |
|
# |
|
namespace eval ::math::special {} |
|
|
|
proc ::math::special::I_n {n x} { |
|
##nagelfar ignore |
|
if { ! [string is integer $n] || $n < 0 } { |
|
error "Wrong order: must be positive integer or zero" |
|
} |
|
|
|
set n2 [expr {$n+8}] ;# Note: just a guess that this will be enough |
|
|
|
set ynp1 0.0 |
|
set yn 1.0 |
|
set sum 1.0 |
|
|
|
while { $n2 > 0 } { |
|
set ynm1 [expr {$ynp1+2.0*$n2*$yn/$x}] |
|
set sum [expr {$sum+$ynm1}] |
|
if { $n2 == $n+1 } { |
|
set result $ynm1 |
|
} |
|
set ynp1 $yn |
|
set yn $ynm1 |
|
incr n2 -1 |
|
} |
|
|
|
set quotient [expr {(2.0*$sum-$ynm1)/exp($x)}] |
|
|
|
expr {$result/$quotient} |
|
} |
|
|
|
# |
|
# some tests -- |
|
# |
|
if { 0 } { |
|
set prec $::tcl_precision |
|
if {![package vsatisfies [package provide Tcl] 9]} { |
|
if {![package vsatisfies [package provide Tcl] 8.5 9]} { |
|
set ::tcl_precision 17 |
|
} else { |
|
set ::tcl_precision 0 |
|
} |
|
} |
|
|
|
foreach x {0.0 2.0 4.4 6.0 10.0 11.0 12.0 13.0 14.0} { |
|
puts "J0($x) = [::math::special::J0 $x] - J1($x) = [::math::special::J1 $x] \ |
|
- J1/2($x) = [::math::special::J1/2 $x]" |
|
} |
|
foreach n {0 1 2 3 4 5} { |
|
puts [::math::special::I_n $n 1.0] |
|
} |
|
|
|
if {![package vsatisfies [package provide Tcl] 9]} { |
|
set ::tcl_precision $prec |
|
} |
|
}
|
|
|