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