1# elliptic.tcl -- 2# Compute elliptic functions and integrals 3# 4# Computation of elliptic functions cn, dn and sn 5# adapted from: 6# Michael W. Pashea 7# Numerical computation of elliptic functions 8# Doctor Dobbs' Journal, May 2005 9# 10 11# namespace ::math::special 12# 13namespace eval ::math::special { 14 namespace export cn sn dn 15 16 ::math::constants::constants pi 17 18 variable halfpi [expr {$pi/2.0}] 19 variable tol 20 21 set tol 1.0e-10 22} 23 24# elliptic_K -- 25# Compute the complete elliptic integral of the first kind 26# 27# Arguments: 28# k Parameter of the integral 29# Result: 30# Value of K(k) 31# Note: 32# This relies on the arithmetic-geometric mean 33# 34proc ::math::special::elliptic_K {k} { 35 variable halfpi 36 if { $k < 0.0 || $k >= 1.0 } { 37 error "Domain error: must be between 0 (inclusive) and 1 (not inclusive)" 38 } 39 40 if { $k == 0.0 } { 41 return $halfpi 42 } 43 44 set a 1.0 45 set b [expr {sqrt(1.0-$k*$k)}] 46 47 for {set i 0} {$i < 10} {incr i} { 48 set anew [expr {($a+$b)/2.0}] 49 set bnew [expr {sqrt($a*$b)}] 50 51 set a $anew 52 set b $bnew 53 #puts "$a $b" 54 } 55 56 return [expr {$halfpi/$a}] 57} 58 59# elliptic_E -- 60# Compute the complete elliptic integral of the second kind 61# 62# Arguments: 63# k Parameter of the integral 64# Result: 65# Value of E(k) 66# Note: 67# This relies on the arithmetic-geometric mean 68# 69proc ::math::special::elliptic_E {k} { 70 variable halfpi 71 if { $k < 0.0 || $k >= 1.0 } { 72 error "Domain error: must be between 0 (inclusive) and 1 (not inclusive)" 73 } 74 75 if { $k == 0.0 } { 76 return $halfpi 77 } 78 if { $k == 1.0 } { 79 return 1.0 80 } 81 82 set a 1.0 83 set b [expr {sqrt(1.0-$k*$k)}] 84 set sumc [expr {$k*$k/2.0}] 85 set factor 0.25 86 87 for {set i 0} {$i < 10} {incr i} { 88 set anew [expr {($a+$b)/2.0}] 89 set bnew [expr {sqrt($a*$b)}] 90 set sumc [expr {$sumc+$factor*($a-$b)*($a-$b)}] 91 set factor [expr {$factor*2.0}] 92 93 set a $anew 94 set b $bnew 95 #puts "$a $b" 96 } 97 98 set Kk [expr {$halfpi/$a}] 99 return [expr {(1.0-$sumc)*$Kk}] 100} 101 102namespace eval ::math::special { 103} 104 105# Nextk -- 106# Auxiliary function for computing next value of k 107# 108# Arguments: 109# k Parameter 110# Return value: 111# Next value to be used 112# 113proc ::math::special::Nextk { k } { 114 set ksq [expr {sqrt(1.0-$k*$k)}] 115 return [expr {(1.0-$ksq)/(1+$ksq)}] 116} 117 118# IterateUK -- 119# Auxiliary function to compute the raw value (phi) 120# 121# Arguments: 122# u Independent variable 123# k Parameter 124# Return value: 125# phi 126# 127proc ::math::special::IterateUK { u k } { 128 variable tol 129 set kvalues {} 130 set nmax 1 131 while { $k > $tol } { 132 set k [Nextk $k] 133 set kvalues [concat $k $kvalues] 134 set u [expr {$u*2.0/(1.0+$k)}] 135 incr nmax 136 #puts "$nmax -$u - $k" 137 } 138 foreach k $kvalues { 139 set u [expr {( $u + asin($k*sin($u)) )/2.0}] 140 } 141 return $u 142} 143 144# cn -- 145# Compute the elliptic function cn 146# 147# Arguments: 148# u Independent variable 149# k Parameter 150# Return value: 151# cn(u,k) 152# Note: 153# If k == 1, then the iteration does not stop 154# 155proc ::math::special::cn { u k } { 156 if { $k > 1.0 } { 157 return -code error "Parameter out of range - must be <= 1.0" 158 } 159 if { $k == 1.0 } { 160 return [expr {1.0/cosh($u)}] 161 } else { 162 set u [IterateUK $u $k] 163 return [expr {cos($u)}] 164 } 165} 166 167# sn -- 168# Compute the elliptic function sn 169# 170# Arguments: 171# u Independent variable 172# k Parameter 173# Return value: 174# sn(u,k) 175# Note: 176# If k == 1, then the iteration does not stop 177# 178proc ::math::special::sn { u k } { 179 if { $k > 1.0 } { 180 return -code error "Parameter out of range - must be <= 1.0" 181 } 182 if { $k == 1.0 } { 183 return [expr {tanh($u)}] 184 } else { 185 set u [IterateUK $u $k] 186 return [expr {sin($u)}] 187 } 188} 189 190# dn -- 191# Compute the elliptic function sn 192# 193# Arguments: 194# u Independent variable 195# k Parameter 196# Return value: 197# dn(u,k) 198# Note: 199# If k == 1, then the iteration does not stop 200# 201proc ::math::special::sn { u k } { 202 if { $k > 1.0 } { 203 return -code error "Parameter out of range - must be <= 1.0" 204 } 205 if { $k == 1.0 } { 206 return [expr {1.0/cosh($u)}] 207 } else { 208 set u [IterateUK $u $k] 209 return [expr {sqrt(1.0-$k*$k*sin($u))}] 210 } 211} 212 213 214# main -- 215# Simple tests 216# 217if { 0 } { 218puts "Special cases:" 219puts "cos(1): [::math::special::cn 1.0 0.0] -- [expr {cos(1.0)}]" 220puts "1/cosh(1): [::math::special::cn 1.0 0.999] -- [expr {1.0/cosh(1.0)}]" 221} 222 223# some tests -- 224# 225if { 0 } { 226set tcl_precision 17 227#foreach k {0.0 0.1 0.2 0.4 0.6 0.8 0.9} { 228# puts "$k: [::math::special::elliptic_K $k]" 229#} 230foreach k2 {0.0 0.1 0.2 0.4 0.6 0.8 0.9} { 231 set k [expr {sqrt($k2)}] 232 puts "$k2: [::math::special::elliptic_K $k] \ 233[::math::special::elliptic_E $k]" 234} 235} 236 237