Raymond E. Rogers wrote: > Either you or I could add a ellint_[KE]comp comment to the ellipke > code. I could enter the k^2/m directions in the ellint_[KE]comp > documents; or perhaps that would go in both cases?
I think there is often confusion regarding whether elliptic integrals are defined with their input argument squared or not, ie. whether it's E(k) = \int_0^\pi/2 dt \sqrt((1 - k^2 \sin^2(t))) or E(k) = \int_0^\pi/2 dt \sqrt((1 - k \sin^2(t))) With a little effort, people using either ellint_[KE]comp or ellipke will find how each function is defined, but I think we should decide whether it's worth explicitly documenting the definitions used within each function. To me the arguments in favor of adding this documentation are: (i) it's a common point of confusion, (ii) that among Octave-Forge functions different definitions are used. An argument against adding this documentation might be that giving references to Abramowitz and Stegun and the GSL documentation is sufficient. > I went back and looked at the ellipke code and your changes. For > some reason the changes don't jibe with my source-forge copy. It's > like they are already half done even though my source copy was svn > on 4/28 like your diff header. Perhaps somebody inserted the fixes? > In any case they don't seem to work; although the only failure > seems to be [x,y,z]. I am obviously out of sync somewhere. Here > are my test cases: > [a,b]=ellipke([[.3, .5,.6];[.4,.6,.6];[.4,.5,.6]]) > [a,b]=ellipke([[.3, .5,.6];[.4,.6,.6]]) > [a,b]=ellipke([[.3, .5];[.4,.6];[.4,.5]]) > [a,b]=ellipke([.3;.5]) > [a,b]=ellipke([.3,.5]) With the ellipke.m below, all these cases work, whereas with the SVN version, the last case fails. Kristjan ## Copyright (C) 2001 David Billinghurst ## ## This program is free software; you can redistribute it and/or modify ## it under the terms of the GNU General Public License as published by ## the Free Software Foundation; either version 2 of the License, or ## (at your option) any later version. ## ## This program is distributed in the hope that it will be useful, ## but WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ## GNU General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with this program; If not, see <http://www.gnu.org/licenses/>. ## Compute: ## complete elliptic integral of first K(m) ## complete elliptic integral of second E(m) ## ## usage: [k,e] = ellipke(m[,tol]) ## ## m is either real array or scalar with 0 <= m <= 1 ## ## tol Ignored. ## (Matlab uses this to allow faster, less accurate approximation) ## ## Ref: Abramowitz, Milton and Stegun, Irene A ## Handbook of Mathematical Functions, Dover, 1965 ## Chapter 17 ## ## See also: ellipj ## Author: David Billinghurst <[EMAIL PROTECTED]> ## Created: 31 January 2001 ## 2001-02-01 Paul Kienzle ## * vectorized ## * included function name in error messages ## 2003-1-18 Jaakko Ruohio ## * extended for m < 0 function [k,e] = ellipke( m ) if (nargin < 1 || nargin > 2) usage("[k, e] = ellipke (m)"); endif k = e = zeros(size(m)); if any(~isreal(m)) error("ellipke must have real m"); endif if any(m>1) error("ellipke must have m <= 1"); endif Nmax = 16; idx = find(m == 1); if (!isempty(idx)) k(idx) = Inf; e(idx) = 1.0; endif idx = find(m == -Inf); if (!isempty(idx)) k(idx) = 0.0; e(idx) = Inf; endif ## Arithmetic-Geometric Mean (AGM) algorithm ## ( Abramowitz and Stegun, Section 17.6 ) idx = find(m != 1 & m != -Inf); if (!isempty(idx)) idx_neg = find(m < 0 & m != -Inf); mult_k = 1./sqrt(1-m(idx_neg)); mult_e = sqrt(1-m(idx_neg)); m(idx_neg) = -m(idx_neg)./(1-m(idx_neg)); a = ones(length(idx),1); b = NaN(length(idx),1); b(:,1) = sqrt(1.0-m(idx)); c = NaN(length(idx),1); c(:,1) = sqrt(m(idx)); f = 0.5; sum = f*c.*c; for n = 2:Nmax t = (a+b)/2; c = (a-b)/2; b = sqrt(a.*b); a = t; f = f * 2; sum = sum + f*c.*c; if all(c./a < eps), break; endif endfor if n >= Nmax, error("ellipke: not enough workspace"); endif k(idx) = 0.5*pi./a; e(idx) = 0.5*pi.*(1.0-sum)./a; k(idx_neg) = mult_k.*k(idx_neg); e(idx_neg) = mult_e.*e(idx_neg); endif endfunction %!test %! ## Test complete elliptic functions of first and second kind %! ## against "exact" solution from Mathematica 3.0 %! ## %! ## David Billinghurst <[EMAIL PROTECTED]> %! ## 1 February 2001 %! m = [0.0; 0.01; 0.1; 0.5; 0.9; 0.99; 1.0 ]; %! [k,e] = ellipke(m); %! %! # K(1.0) is really infinity - see below %! K = [ %! 1.5707963267948966192; %! 1.5747455615173559527; %! 1.6124413487202193982; %! 1.8540746773013719184; %! 2.5780921133481731882; %! 3.6956373629898746778; %! 0.0 ]; %! E = [ %! 1.5707963267948966192; %! 1.5668619420216682912; %! 1.5307576368977632025; %! 1.3506438810476755025; %! 1.1047747327040733261; %! 1.0159935450252239356; %! 1.0 ]; %! if k(7)==Inf, k(7)=0.0; endif; %! assert(K,k,8*eps); %! assert(E,e,8*eps); ------------------------------------------------------------------------- This SF.net email is sponsored by the 2008 JavaOne(SM) Conference Don't miss this year's exciting event. There's still time to save $100. Use priority code J8TL2D2. http://ad.doubleclick.net/clk;198757673;13503038;p?http://java.sun.com/javaone _______________________________________________ Octave-dev mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/octave-dev
