1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188
|
# special.tcl --
# Provide well-known special mathematical functions
#
# This file contains a collection of tests for one or more of the Tcllib
# procedures. Sourcing this file into Tcl runs the tests and
# generates output for errors. No output means no errors were found.
#
# Copyright (c) 2004 by Arjen Markus. All rights reserved.
#
# RCS: @(#) $Id: special.tcl,v 1.9 2005/10/06 05:16:37 andreas_kupries Exp $
#
package require math
package require math::constants
package require math::statistics
# 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
variable halfpi [expr {$pi/2.0}]
#
# Functions defined in other math submodules
#
if { [info commands Beta] == {} } {
namespace import ::math::Beta
}
if { [info commands Beta] == {} } {
namespace import ::math::ln_Gamma
}
#
# Export the various functions
#
namespace export Beta ln_Gamma Gamma erf erfc fresnel_C fresnel_S sinc
}
# Gamma --
# The Gamma function - synonym for "factorial"
#
proc ::math::special::Gamma {x} {
::math::factorial [expr { $x + 1 }]
}
# erf --
# The error function
# Arguments:
# x The value for which the function must be evaluated
# Result:
# erf(x)
#
proc ::math::special::erf {x} {
set x2 $x
if { $x > 0.0 } {
set x2 [expr {-$x}]
}
if { $x2 != 0.0 } {
set r [::math::statistics::cdf-normal 0.0 [expr {sqrt(0.5)}] $x2]
if { $x > 0.0 } {
return [expr {1.0-2.0*$r}]
} else {
return [expr {2.0*$r-1.0}]
}
} else {
return 0.0
}
}
# erfc --
# The complement of the error function
# Arguments:
# x The value for which the function must be evaluated
# Result:
# erfc(x) = 1.0-erf(x)
#
proc ::math::special::erfc {x} {
set x2 $x
if { $x > 0.0 } {
set x2 [expr {-$x}]
}
if { $x2 != 0.0 } {
set r [::math::statistics::cdf-normal 0.0 [expr {sqrt(0.5)}] $x2]
if { $x > 0.0 } {
return [expr {2.0*$r}]
} else {
return [expr {2.0-2.0*$r}]
}
} else {
return 1.0
}
}
# ComputeFG --
# Compute the auxiliary functions f and g
#
# Arguments:
# x Parameter of the integral (x>=0)
# Result:
# Approximate values for f and g
# Note:
# See Abramowitz and Stegun. The accuracy is 2.0e-3.
#
proc ::math::special::ComputeFG {x} {
list [expr {(1.0+0.926*$x)/(2.0+1.792*$x+3.104*$x*$x)}] \
[expr {1.0/(2.0+4.142*$x+3.492*$x*$x+6.670*$x*$x*$x)}]
}
# fresnel_C --
# Compute the Fresnel cosine integral
#
# Arguments:
# x Parameter of the integral (x>=0)
# Result:
# Value of C(x) = integral from 0 to x of cos(0.5*pi*x^2)
# Note:
# This relies on a rational approximation of the two auxiliary functions f and g
#
proc ::math::special::fresnel_C {x} {
variable halfpi
if { $x < 0.0 } {
error "Domain error: x must be non-negative"
}
if { $x == 0.0 } {
return 0.0
}
foreach {f g} [ComputeFG $x] {break}
set xarg [expr {$halfpi*$x*$x}]
return [expr {0.5+$f*sin($xarg)-$g*cos($xarg)}]
}
# fresnel_S --
# Compute the Fresnel sine integral
#
# Arguments:
# x Parameter of the integral (x>=0)
# Result:
# Value of S(x) = integral from 0 to x of sin(0.5*pi*x^2)
# Note:
# This relies on a rational approximation of the two auxiliary functions f and g
#
proc ::math::special::fresnel_S {x} {
variable halfpi
if { $x < 0.0 } {
error "Domain error: x must be non-negative"
}
if { $x == 0.0 } {
return 0.0
}
foreach {f g} [ComputeFG $x] {break}
set xarg [expr {$halfpi*$x*$x}]
return [expr {0.5-$f*cos($xarg)-$g*sin($xarg)}]
}
# sinc --
# Compute the sinc function
# Arguments:
# x Value of the argument
# Result:
# sin(x)/x
#
proc ::math::special::sinc {x} {
if { $x == 0.0 } {
return 1.0
} else {
return [expr {sin($x)/$x}]
}
}
# Bessel functions and elliptic integrals --
#
source [file join [file dirname [info script]] "bessel.tcl"]
source [file join [file dirname [info script]] "classic_polyns.tcl"]
source [file join [file dirname [info script]] "elliptic.tcl"]
source [file join [file dirname [info script]] "exponential.tcl"]
package provide math::special 0.2
|