www.gusucode.com > elfun18工具箱matlab源码程序 > elfun18/elfun18v1_3/elfun18v1_2/elfun18/el1.m
function result = el1( x, kc ) %EL1 Evaluates the Bulirsch's elliptic integral of 1st kind % % x % | | % | dt % el1(x,kc) = | ------------------------------ % | sqrt((1 + t^2)*(1 + kc^2*t^2)) % | | % 0 % % Result: % el1(x,kc) -- real scalar or NaN if either argument is invalid or % convergence failes. % % Arguments: % x -- real scalar % kc -- real scalar, complementary modulus % % Functions called: % cel1 % % Matlab functions called: % abs, asinh, atan, isinf, isnan, pi, sqrt % % Reference: % [1] R.Bulirsch, "Numerical calculation of elliptic integrals and % elliptic functions", Numerische Mathematik 7, 1965, pp 78-90 % % Check input if isnan(x) || isnan(kc) result = NaN; return end % Special cases if isinf(kc) result = 0; return end if isinf(x) %|| abs(x) > 1e-16*max(1,1/abs(kc)) result = sign(x)*cel1(kc); return end if x == 0 result = 0; return end if kc == 0 result = asinh(x); return end kc = abs(kc); if kc == 1 result = atan(x); return end % General case is translation on Algol procedure el1 from [1]. % CA = 10^(-D/2) and CB = 10^(-(D+2)) where D is number of % significant digits. CA = 0.5e-8; % by tests. CB = 1e-19; y = abs(1/x); m = 1; n = 0; while true e = m*kc; g = m; m = m + kc; y = y - e/y; if y == 0 y = sqrt(e)*CB; end if abs(g - kc) <= CA*g break end kc = sqrt(e)*2; if isinf(kc) || kc == 0 result = NaN; return end n = n + n; if y < 0 n = n + 1; end end if y < 0 n = n + 1; end e = (atan(m/y) + pi*n)/m; if x < 0 y = -e; else y = e; end result = y; end