bessel.tcl

Go to the documentation of this file.

/*  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 ::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)*/
/* */
ret  ::math::special::J0 (type x) {
    Jn 0 $x
}

/*  J1 --*/
/*     First-order Bessel function*/
/* */
/*  Arguments:*/
/*     x         Value of the x-coordinate*/
/*  Result:*/
/*     Value of J1(x)*/
/* */
ret  ::math::special::J1 (type 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).*/
/* */
ret  ::math::special::Jn (type n , type x) {
    variable pi

    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)*/
/* */
ret  ::math::special::J1/2 (type 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)*/
/* */
ret  ::math::special::J-1/2 (type 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 ::math::special {}

ret  ::math::special::I_n (type n , type x) {
    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 } {
 tcl = _precision 17
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]
}
}