Artifact 1cc2bdddd97b447ba98ad5f935d2ca25e4a5392f:
# elliptic.tcl --
# Compute elliptic functions and integrals
#
# Computation of elliptic functions cn, dn and sn
# adapted from:
# Michael W. Pashea
# Numerical computation of elliptic functions
# Doctor Dobbs' Journal, May 2005
#
# namespace ::math::special
#
namespace eval ::math::special {
namespace export cn sn dn
::math::constants::constants pi
variable halfpi [expr {$pi/2.0}]
variable tol
set tol 1.0e-10
}
# elliptic_K --
# Compute the complete elliptic integral of the first kind
#
# Arguments:
# k Parameter of the integral
# Result:
# Value of K(k)
# Note:
# This relies on the arithmetic-geometric mean
#
proc ::math::special::elliptic_K {k} {
variable halfpi
if { $k < 0.0 || $k >= 1.0 } {
error "Domain error: must be between 0 (inclusive) and 1 (not inclusive)"
}
if { $k == 0.0 } {
return $halfpi
}
set a 1.0
set b [expr {sqrt(1.0-$k*$k)}]
for {set i 0} {$i < 10} {incr i} {
set anew [expr {($a+$b)/2.0}]
set bnew [expr {sqrt($a*$b)}]
set a $anew
set b $bnew
#puts "$a $b"
}
return [expr {$halfpi/$a}]
}
# elliptic_E --
# Compute the complete elliptic integral of the second kind
#
# Arguments:
# k Parameter of the integral
# Result:
# Value of E(k)
# Note:
# This relies on the arithmetic-geometric mean
#
proc ::math::special::elliptic_E {k} {
variable halfpi
if { $k < 0.0 || $k >= 1.0 } {
error "Domain error: must be between 0 (inclusive) and 1 (not inclusive)"
}
if { $k == 0.0 } {
return $halfpi
}
if { $k == 1.0 } {
return 1.0
}
set a 1.0
set b [expr {sqrt(1.0-$k*$k)}]
set sumc [expr {$k*$k/2.0}]
set factor 0.25
for {set i 0} {$i < 10} {incr i} {
set anew [expr {($a+$b)/2.0}]
set bnew [expr {sqrt($a*$b)}]
set sumc [expr {$sumc+$factor*($a-$b)*($a-$b)}]
set factor [expr {$factor*2.0}]
set a $anew
set b $bnew
#puts "$a $b"
}
set Kk [expr {$halfpi/$a}]
return [expr {(1.0-$sumc)*$Kk}]
}
namespace eval ::math::special {
}
# Nextk --
# Auxiliary function for computing next value of k
#
# Arguments:
# k Parameter
# Return value:
# Next value to be used
#
proc ::math::special::Nextk { k } {
set ksq [expr {sqrt(1.0-$k*$k)}]
return [expr {(1.0-$ksq)/(1+$ksq)}]
}
# IterateUK --
# Auxiliary function to compute the raw value (phi)
#
# Arguments:
# u Independent variable
# k Parameter
# Return value:
# phi
#
proc ::math::special::IterateUK { u k } {
variable tol
set kvalues {}
set nmax 1
while { $k > $tol } {
set k [Nextk $k]
set kvalues [concat $k $kvalues]
set u [expr {$u*2.0/(1.0+$k)}]
incr nmax
#puts "$nmax -$u - $k"
}
foreach k $kvalues {
set u [expr {( $u + asin($k*sin($u)) )/2.0}]
}
return $u
}
# cn --
# Compute the elliptic function cn
#
# Arguments:
# u Independent variable
# k Parameter
# Return value:
# cn(u,k)
# Note:
# If k == 1, then the iteration does not stop
#
proc ::math::special::cn { u k } {
if { $k > 1.0 } {
return -code error "Parameter out of range - must be <= 1.0"
}
if { $k == 1.0 } {
return [expr {1.0/cosh($u)}]
} else {
set u [IterateUK $u $k]
return [expr {cos($u)}]
}
}
# sn --
# Compute the elliptic function sn
#
# Arguments:
# u Independent variable
# k Parameter
# Return value:
# sn(u,k)
# Note:
# If k == 1, then the iteration does not stop
#
proc ::math::special::sn { u k } {
if { $k > 1.0 } {
return -code error "Parameter out of range - must be <= 1.0"
}
if { $k == 1.0 } {
return [expr {tanh($u)}]
} else {
set u [IterateUK $u $k]
return [expr {sin($u)}]
}
}
# dn --
# Compute the elliptic function sn
#
# Arguments:
# u Independent variable
# k Parameter
# Return value:
# dn(u,k)
# Note:
# If k == 1, then the iteration does not stop
#
proc ::math::special::sn { u k } {
if { $k > 1.0 } {
return -code error "Parameter out of range - must be <= 1.0"
}
if { $k == 1.0 } {
return [expr {1.0/cosh($u)}]
} else {
set u [IterateUK $u $k]
return [expr {sqrt(1.0-$k*$k*sin($u))}]
}
}
# main --
# Simple tests
#
if { 0 } {
puts "Special cases:"
puts "cos(1): [::math::special::cn 1.0 0.0] -- [expr {cos(1.0)}]"
puts "1/cosh(1): [::math::special::cn 1.0 0.999] -- [expr {1.0/cosh(1.0)}]"
}
# some tests --
#
if { 0 } {
set tcl_precision 17
#foreach k {0.0 0.1 0.2 0.4 0.6 0.8 0.9} {
# puts "$k: [::math::special::elliptic_K $k]"
#}
foreach k2 {0.0 0.1 0.2 0.4 0.6 0.8 0.9} {
set k [expr {sqrt($k2)}]
puts "$k2: [::math::special::elliptic_K $k] \
[::math::special::elliptic_E $k]"
}
}