# HG changeset patch # User Mike Miller # Date 1367189451 14400 # Node ID f423873d32755e17a98e93cddfe3bc708bb72d8b # Parent 1a3bfb14b5da2c0932e0cf9465cab7e58ca82fcb Style fixes for ellipj.cc, ellipke.m, and expint.m * ellipj.cc, ellipke.m, expint.m: Style fixes. diff -r 1a3bfb14b5da -r f423873d3275 libinterp/corefcn/ellipj.cc --- a/libinterp/corefcn/ellipj.cc Sun Apr 28 17:08:11 2013 -0400 +++ b/libinterp/corefcn/ellipj.cc Sun Apr 28 18:50:51 2013 -0400 @@ -1,40 +1,23 @@ /* - Copyright (C) 2001 Leopoldo Cerbaro - - 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 3 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. +Copyright (C) 2001 Leopoldo Cerbaro - You should have received a copy of the GNU General Public License - along with this program; If not, see . - - Compute the Jacobi elliptic functions sn(u|m), cn(u|m) and dn(u|m) for - argument u (real or complex) and parameter m. +This file is part of Octave. - usage: [sn,cn,dn] = ellipj(u,m[,tol]) - - u and can be complex. - m is restricted to 0 <= m <= 1. - They can be scalars, matrix and scalar, scalar and matrix, - column and row, conformant matrices. +Octave 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 3 of the License, or (at your +option) any later version. - modified so u can be complex. Leopoldo Cerbaro redbliss@libero.it - - Ref: Abramowitz, Milton and Stegun, Irene A - Handbook of Mathematical Functions, Dover, 1965 - Chapter 16 (Sections 16.4, 16.13 and 16.15) +Octave 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. - Based upon ellipj.m made by David Billinghurst - and besselj.cc +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +. - Author: Leopoldo Cerbaro - Created: 15 December 2001 */ #ifdef HAVE_CONFIG_H @@ -46,115 +29,114 @@ #include "lo-ieee.h" static void -gripe_ellipj_arg ( const char *arg) +gripe_ellipj_arg (const char *arg) { error ("ellipj: expecting scalar or matrix as %s argument", arg); } -const double eps = 2.220446049e-16; -const int Nmax = 16; +static void +sncndn (double u, double m, double& sn, double& cn, double& dn, double& err) +{ + static const int Nmax = 16; + double m1, t=0, si_u, co_u, se_u, ta_u, b, c[Nmax], a[Nmax], phi; + int n, Nn, ii; -static void -sncndn ( double u, double m, double& sn, double& cn, double& dn, double& err) { -/* real */ -double sqrt_eps, m1, t=0., si_u, co_u, se_u, ta_u, b, c[Nmax], a[Nmax], phi; -int n, Nn, ii; + if (m < 0 || m > 1) + { + warning ("ellipj: expecting 0 <= m <= 1"); /* -lc- */ + sn = cn = dn = lo_ieee_nan_value (); + return; + } - if (m < 0. || m > 1.) { - warning ("ellipj: expecting 0. <= m <= 1."); /* -lc- */ - sn = cn = dn = lo_ieee_nan_value (); - return; - } - sqrt_eps = sqrt(eps); - if (m < sqrt_eps) { - /* # For small m, ( Abramowitz and Stegun, Section 16.13 ) */ - /*{{{*/ - si_u = sin(u); - co_u = cos(u); - t = 0.25*m*(u-si_u*co_u); - sn = si_u - t * co_u; - cn = co_u + t * si_u; - dn = 1.0 - 0.5*m*si_u*si_u; -/*}}}*/ - } else if ( (1.0 - m) < sqrt_eps ) { - /* For m1 = (1-m) small ( Abramowitz and Stegun, Section 16.15 ) */ - /*{{{*/ - m1 = 1.0-m; - si_u = sinh(u); - co_u = cosh(u); - ta_u = tanh(u); - se_u = 1.0/co_u; - sn = ta_u + 0.25*m1*(si_u*co_u-u)*se_u*se_u; - cn = se_u - 0.25*m1*(si_u*co_u-u)*ta_u*se_u; - dn = se_u + 0.25*m1*(si_u*co_u+u)*ta_u*se_u; -/*}}}*/ - } else { - /*{{{*/ - /* - // Arithmetic-Geometric Mean (AGM) algorithm - // ( Abramowitz and Stegun, Section 16.4 ) - */ - - a[0] = 1.0; - b = sqrt(1.0-m); - c[0] = sqrt(m); - for (n = 1; n::epsilon ()); + if (m < sqrt_eps) + { + /* # For small m, ( Abramowitz and Stegun, Section 16.13 ) */ + si_u = sin (u); + co_u = cos (u); + t = 0.25*m*(u - si_u*co_u); + sn = si_u - t * co_u; + cn = co_u + t * si_u; + dn = 1 - 0.5*m*si_u*si_u; + } + else if ((1 - m) < sqrt_eps) + { + /* For m1 = (1-m) small ( Abramowitz and Stegun, Section 16.15 ) */ + m1 = 1 - m; + si_u = sinh (u); + co_u = cosh (u); + ta_u = tanh (u); + se_u = 1/co_u; + sn = ta_u + 0.25*m1*(si_u*co_u - u)*se_u*se_u; + cn = se_u - 0.25*m1*(si_u*co_u - u)*ta_u*se_u; + dn = se_u + 0.25*m1*(si_u*co_u + u)*ta_u*se_u; + } + else + { + /* + // Arithmetic-Geometric Mean (AGM) algorithm + // ( Abramowitz and Stegun, Section 16.4 ) + */ + + a[0] = 1; + b = sqrt (1 - m); + c[0] = sqrt (m); + for (n = 1; n < Nmax; ++n) + { + a[n] = (a[n - 1] + b)/2; + c[n] = (a[n - 1] - b)/2; + b = sqrt (a[n - 1]*b); + if (c[n]/a[n] < std::numeric_limits::epsilon ()) break; } - if ( n >= Nmax-1) { - // fprintf(stderr, "Not enough workspace\n"); - err = 1.; - return; + if (n >= Nmax - 1) + { + err = 1; + return; } - Nn = n; - for ( ii = 1; n>0; ii = ii*2, --n) ; // pow(2, Nn) - phi = ii*a[Nn]*u; - for ( n = Nn; n > 0; --n) { + Nn = n; + for (ii = 1; n > 0; ii = ii*2, --n) ; // ii = pow(2,Nn) + phi = ii*a[Nn]*u; + for (n = Nn; n > 0; --n) + { t = phi; - phi = (asin((c[n]/a[n])* sin(phi))+phi)/2.; + phi = (asin ((c[n]/a[n])* sin (phi)) + phi)/2; } - sn = sin(phi); - cn = cos(phi); - dn = cn/cos(t-phi); -/*}}}*/ - } - return; + sn = sin (phi); + cn = cos (phi); + dn = cn/cos (t - phi); + } } static void -sncndn ( Complex& u, double m, - Complex& sn, Complex& cn, Complex& dn, double& err) { -double m1 = 1.-m, ss1, cc1, dd1; +sncndn (Complex& u, double m, Complex& sn, Complex& cn, Complex& dn, + double& err) +{ + double m1 = 1 - m, ss1, cc1, dd1; - sncndn( imag(u), m1, ss1, cc1, dd1, err); - if ( real(u) == 0.) { - /* u is pure imag: Jacoby imag. transf. */ - /*{{{*/ - sn = Complex (0. , ss1/cc1); - cn = 1/cc1; // cn.imag = 0.; - dn = dd1/cc1; // dn.imag = 0.; - /*}}}*/ - } else { - /* u is generic complex */ - /*{{{*/ - double ss, cc, dd, ddd; + sncndn (imag (u), m1, ss1, cc1, dd1, err); + if (real (u) == 0) + { + /* u is pure imag: Jacoby imag. transf. */ + sn = Complex (0, ss1/cc1); + cn = 1/cc1; // cn.imag = 0; + dn = dd1/cc1; // dn.imag = 0; + } + else + { + /* u is generic complex */ + double ss, cc, dd, ddd; - sncndn( real(u), m, ss, cc, dd, err); + sncndn (real (u), m, ss, cc, dd, err); ddd = cc1*cc1 + m*ss*ss*ss1*ss1; - sn = Complex (ss*dd1/ddd, cc*dd*ss1*cc1/ddd); + sn = Complex (ss*dd1/ddd, cc*dd*ss1*cc1/ddd); cn = Complex (cc*cc1/ddd, -ss*dd*ss1*dd1/ddd); dn = Complex (dd*cc1*dd1/ddd, -m*ss*cc*ss1/ddd); - /*}}}*/ - } - return; + } } DEFUN (ellipj, args, nargout, "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {[@var{sn}, @var{cn}, @var{dn}] =} \ +@deftypefn {Built-in Function} {[@var{sn}, @var{cn}, @var{dn}] =} \ ellipj (@var{u}, @var{m}, @var{err})\n\ Compute the Jacobi elliptic functions sn, cn, dn of complex argument and real parameter.\n\ \n\ @@ -166,7 +148,7 @@ @var{m} must conform and the results will be the same size.\n\ \n\ The value of @var{u} may be complex.\n\ -The value of @var{m} must be 0 <= m <= 1. .\n\ +The value of @var{m} must be 0 <= m <= 1.\n\ \n\ If requested, @var{err} contains the following status information\n\ and is the same size as the result.\n\ @@ -178,218 +160,258 @@ Error---no computation, algorithm termination condition not met,\n\ return @code{NaN}.\n\ @end enumerate\n\ + Ref: Abramowitz, Milton and Stegun, Irene A\n\ + Handbook of Mathematical Functions, Dover, 1965\n\ + Chapter 16 (Sections 16.4, 16.13 and 16.15)\n\ @end deftypefn") { octave_value_list retval; int nargin = args.length (); - if (nargin == 2 ) { - octave_value u_arg = args(0); - octave_value m_arg = args(1); + if (nargin != 2 ) + { + print_usage (); + return retval; + } - if (m_arg.is_scalar_type ()) { // m is scalar - double m = args(1).double_value (); + octave_value u_arg = args(0); + octave_value m_arg = args(1); - if (! error_state) { + if (m_arg.is_scalar_type ()) + { + double m = args(1).double_value (); - if (u_arg.is_scalar_type ()) { /* u scalar */ - /*{{{*/ - if (u_arg.is_real_type ()) { // u real - double u = args(0).double_value (); + if (error_state) + { + gripe_ellipj_arg ("second"); + return retval; + } - if (! error_state) { - double sn, cn, dn; - double err=0; - octave_value result; + if (u_arg.is_scalar_type ()) + { + if (u_arg.is_real_type ()) + { // u real + double u = args(0).double_value (); - sncndn(u, m, sn, cn, dn, err); - retval (0) = sn; - retval (1) = cn; - retval (2) = dn; - if (nargout > 3) - retval(3) = err; - } else - gripe_ellipj_arg ( "first"); + if (error_state) + { + gripe_ellipj_arg ("first"); + return retval; + } + double sn, cn, dn; + double err = 0; - } else { // u complex + sncndn (u, m, sn, cn, dn, err); + retval (0) = sn; + retval (1) = cn; + retval (2) = dn; + if (nargout > 3) retval(3) = err; + } + else + { // u complex Complex u = u_arg.complex_value (); - if (! error_state) { - Complex sn, cn, dn; - double err; - octave_value result; + if (error_state) + { + gripe_ellipj_arg ("second"); + return retval; + } + + Complex sn, cn, dn; + double err; + + sncndn (u, m, sn, cn, dn, err); - sncndn( u, m, sn, cn, dn, err); + retval (0) = sn; + retval (1) = cn; + retval (2) = dn; + if (nargout > 3) retval(3) = err; + } + } + else + { /* u is matrix ( m is scalar ) */ + ComplexMatrix u = u_arg.complex_matrix_value (); + + if (error_state) + { + gripe_ellipj_arg ("first"); + return retval; + } + + octave_idx_type nr = u.rows (); + octave_idx_type nc = u.cols (); + + ComplexMatrix sn (nr, nc), cn (nr, nc), dn (nr, nc); + Matrix err (nr, nc); + + for (octave_idx_type j = 0; j < nc; j++) + for (octave_idx_type i = 0; i < nr; i++) + sncndn (u(i,j), m, sn(i,j), cn(i,j), dn(i,j), err(i,j)); - retval (0) = sn; - retval (1) = cn; - retval (2) = dn; - if (nargout > 3) retval(3) = err; - } else - gripe_ellipj_arg ( "second"); + retval (0) = sn; + retval (1) = cn; + retval (2) = dn; + if (nargout > 3) retval(3) = err; + } + } + else + { + Matrix m = args(1).matrix_value (); + + if (error_state) + { + gripe_ellipj_arg ("second"); + return retval; + } + + octave_idx_type mr = m.rows (); + octave_idx_type mc = m.cols (); + + if (u_arg.is_scalar_type ()) + { /* u is scalar */ + octave_idx_type nr = m.rows (); + octave_idx_type nc = m.cols (); + Matrix err (nr, nc); + + if (u_arg.is_real_type ()) + { + double u = u_arg.double_value (); + + if (error_state) + { + gripe_ellipj_arg ("first"); + return retval; + } + + Matrix sn (nr, nc), cn (nr, nc), dn (nr, nc); + for (octave_idx_type j = 0; j < nc; j++) + for (octave_idx_type i = 0; i < nr; i++) + sncndn (u, m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j)); + + retval (0) = sn; + retval (1) = cn; + retval (2) = dn; + if (nargout > 3) retval(3) = err; } - /*}}}*/ - } else { /* u is matrix ( m is scalar ) */ - /*{{{*/ - ComplexMatrix u = u_arg.complex_matrix_value (); - - if (! error_state) { - octave_value result; - int nr = u.rows (); - int nc = u.cols (); + else + { + Complex u = u_arg.complex_value (); + if (error_state) + { + gripe_ellipj_arg ("first"); + return retval; + } ComplexMatrix sn (nr, nc), cn (nr, nc), dn (nr, nc); - Matrix err (nr, nc); - - for (int j = 0; j < nc; j++) - for (int i = 0; i < nr; i++) - sncndn (u(i,j), m, sn(i,j), cn(i,j), dn(i,j), err(i,j)); - - retval (0) = sn; - retval (1) = cn; - retval (2) = dn; - if (nargout > 3) retval(3) = err; - } else - gripe_ellipj_arg ( "first"); - /*}}}*/ - } - } else - gripe_ellipj_arg ( "second"); - } else { // m is matrix - Matrix m = args(1).matrix_value (); - - if (! error_state) { - int mr = m.rows (); - int mc = m.cols (); - - if (u_arg.is_scalar_type ()) { /* u is scalar */ - /*{{{*/ - octave_value result; - int nr = m.rows (); - int nc = m.cols (); - Matrix err (nr, nc); - - if (u_arg.is_real_type ()) { - double u = u_arg.double_value (); - Matrix sn (nr, nc), cn (nr, nc), dn (nr, nc); - if (! error_state) { - for (int j = 0; j < nc; j++) - for (int i = 0; i < nr; i++) - sncndn (u, m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j)); - - retval (0) = sn; - retval (1) = cn; - retval (2) = dn; - if (nargout > 3) retval(3) = err; - } else - gripe_ellipj_arg ( "first"); - } else { - Complex u = u_arg.complex_value (); - ComplexMatrix sn (nr, nc), cn (nr, nc), dn (nr, nc); - if (! error_state) { - for (int j = 0; j < nc; j++) - for (int i = 0; i < nr; i++) - sncndn (u, m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j)); - retval (0) = sn; - retval (1) = cn; - retval (2) = dn; - if (nargout > 3) retval(3) = err; - } else - gripe_ellipj_arg ( "first"); - } - /*}}}*/ - } else { // u is matrix (m is matrix) - /*{{{*/ - if (u_arg.is_real_type ()) { // u real matrix + for (octave_idx_type j = 0; j < nc; j++) + for (octave_idx_type i = 0; i < nr; i++) + sncndn (u, m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j)); + retval (0) = sn; + retval (1) = cn; + retval (2) = dn; + if (nargout > 3) retval(3) = err; + } + } + else + { // u is matrix (m is matrix) + if (u_arg.is_real_type ()) + { // u real matrix Matrix u = u_arg.matrix_value (); - if (! error_state) { - int ur = u.rows (); - int uc = u.cols (); - - if (mr == 1 && uc == 1) { // u column, m row - RowVector rm = m.row ((octave_idx_type)0); - ColumnVector cu = u.column ((octave_idx_type)0); + if (error_state) + { + gripe_ellipj_arg ("first "); + return retval; + } - Matrix sn (ur, mc), cn (ur, mc), dn (ur, mc); - Matrix err(ur,mc); -// octave_value result; + octave_idx_type ur = u.rows (); + octave_idx_type uc = u.cols (); - for (int j = 0; j < mc; j++) - for (int i = 0; i < ur; i++) - sncndn (cu(i), rm(j), sn(i,j), cn(i,j), dn(i,j), err(i,j)); + if (mr == 1 && uc == 1) + { // u column, m row + RowVector rm = m.row (0); + ColumnVector cu = u.column (0); + + Matrix sn (ur, mc), cn (ur, mc), dn (ur, mc); + Matrix err (ur,mc); - retval (0) = sn; - retval (1) = cn; - retval (2) = dn; - if (nargout > 3) retval(3) = err; - } else if (ur == mr && uc == mc) { - Matrix sn (ur, mc), cn (ur, mc), dn (ur, mc); - Matrix err(ur,mc); -// octave_value result; - - for (int j = 0; j < uc; j++) - for (int i = 0; i < ur; i++) - sncndn (u(i,j), m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j)); + for (octave_idx_type j = 0; j < mc; j++) + for (octave_idx_type i = 0; i < ur; i++) + sncndn (cu(i), rm(j), sn(i,j), cn(i,j), dn(i,j), err(i,j)); - retval (0) = sn; - retval (1) = cn; - retval (2) = dn; - if (nargout > 3) retval(3) = err; - } else - error("u m invalid"); - } else - gripe_ellipj_arg ( "first "); - } else { // u complex matrix - ComplexMatrix u = u_arg.complex_matrix_value (); - if (! error_state) { - int ur = u.rows (); - int uc = u.cols (); + retval (0) = sn; + retval (1) = cn; + retval (2) = dn; + if (nargout > 3) retval(3) = err; + } + else if (ur == mr && uc == mc) + { + Matrix sn (ur, mc), cn (ur, mc), dn (ur, mc); + Matrix err (ur,mc); + + for (octave_idx_type j = 0; j < uc; j++) + for (octave_idx_type i = 0; i < ur; i++) + sncndn (u(i,j), m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j)); - if (mr == 1 && uc == 1) { - RowVector rm = m.row ((octave_idx_type)0); - ComplexColumnVector cu = u.column ((octave_idx_type)0); - - ComplexMatrix sn (ur, mc), cn (ur, mc), dn (ur, mc); - Matrix err(ur,mc); -// octave_value result; + retval (0) = sn; + retval (1) = cn; + retval (2) = dn; + if (nargout > 3) retval(3) = err; + } + else + error ("u m invalid"); + } + else + { // u complex matrix + ComplexMatrix u = u_arg.complex_matrix_value (); + if (error_state) + { + gripe_ellipj_arg ("second"); + return retval; + } - for (int j = 0; j < mc; j++) - for (int i = 0; i < ur; i++) - sncndn (cu(i), rm(j), sn(i,j), cn(i,j), dn(i,j), err(i,j)); + octave_idx_type ur = u.rows (); + octave_idx_type uc = u.cols (); - retval (0) = sn; - retval (1) = cn; - retval (2) = dn; - if (nargout > 3) retval(3) = err; - } else if (ur == mr && uc == mc) { + if (mr == 1 && uc == 1) + { + RowVector rm = m.row (0); + ComplexColumnVector cu = u.column (0); - ComplexMatrix sn (ur, mc), cn (ur, mc), dn (ur, mc); - Matrix err(ur,mc); -// octave_value result; + ComplexMatrix sn (ur, mc), cn (ur, mc), dn (ur, mc); + Matrix err (ur,mc); - for (int j = 0; j < uc; j++) - for (int i = 0; i < ur; i++) - sncndn (u(i,j), m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j)); + for (octave_idx_type j = 0; j < mc; j++) + for (octave_idx_type i = 0; i < ur; i++) + sncndn (cu(i), rm(j), sn(i,j), cn(i,j), dn(i,j), err(i,j)); + + retval (0) = sn; + retval (1) = cn; + retval (2) = dn; + if (nargout > 3) retval(3) = err; + } + else if (ur == mr && uc == mc) + { + ComplexMatrix sn (ur, mc), cn (ur, mc), dn (ur, mc); + Matrix err (ur,mc); - retval (0) = sn; - retval (1) = cn; - retval (2) = dn; - if (nargout > 3) retval(3) = err; - } else - error("u m invalid"); - } else - gripe_ellipj_arg ( "second"); + for (octave_idx_type j = 0; j < uc; j++) + for (octave_idx_type i = 0; i < ur; i++) + sncndn (u(i,j), m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j)); + + retval (0) = sn; + retval (1) = cn; + retval (2) = dn; + if (nargout > 3) retval(3) = err; + } + else + error ("u m invalid"); } - /*}}}*/ - } - } else - gripe_ellipj_arg ( "second"); - } // m matrix - } else // wrong n. of argin - print_usage (); - return retval; + } + } // m matrix + + return retval; } /* @@ -444,7 +466,7 @@ %!test %! k = (tan(pi/8.))^2; m = k*k; -%! SN = [ +%! SN = [ %! -1. + I * 0. , -0.8392965923 + 0. * I %! -1. + I * 0.2 , -0.8559363407 + 0.108250955 * I %! -1. + I * 0.4 , -0.906529758 + 0.2204040232 * I diff -r 1a3bfb14b5da -r f423873d3275 scripts/specfun/ellipke.m --- a/scripts/specfun/ellipke.m Sun Apr 28 17:08:11 2013 -0400 +++ b/scripts/specfun/ellipke.m Sun Apr 28 18:50:51 2013 -0400 @@ -2,25 +2,30 @@ ## Copyright (C) 2001 Paul Kienzle ## Copyright (C) 2003 Jaakko Ruohio ## -## 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 3 of the License, or (at your option) any later -## version. +## This file is part of Octave. +## +## Octave 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 3 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. +## Octave 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 . +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## . ## -*- texinfo -*- -## @deftypefn {Function File} {[@var{k}, @var{e}] =} ellipke (@var{m}[,@var{tol}]) +## @deftypefn {Function File} {} ellipke (@var{m}) +## @deftypefnx {Function File} {} ellipke (@var{m}, @var{tol}) +## @deftypefnx {Function File} {[@var{k}, @var{e}] =} ellipke (@dots{}) ## Compute complete elliptic integral of first K(@var{m}) and second E(@var{m}). ## ## @var{m} is either real array or scalar with 0 <= m <= 1 -## +## ## @var{tol} will be ignored (@sc{Matlab} uses this to allow faster, less ## accurate approximation) ## @@ -29,91 +34,89 @@ ## @seealso{ellipj} ## @end deftypefn -function [k,e] = ellipke( m ) +function [k, e] = ellipke (m) if (nargin < 1 || nargin > 2) - print_usage; + print_usage (); endif - k = e = zeros(size(m)); + k = e = zeros (size (m)); m = m(:); - if any(~isreal(m)) - error("ellipke must have real m"); + if any (!isreal (m)) + error ("ellipke must have real m"); endif - if any(m>1) - error("ellipke must have m <= 1"); + if any (m > 1) + error ("ellipke must have m <= 1"); endif Nmax = 16; - idx = find(m == 1); - if (!isempty(idx)) + idx = find (m == 1); + if (!isempty (idx)) k(idx) = Inf; - e(idx) = 1.0; + e(idx) = 1; endif - - idx = find(m == -Inf); - if (!isempty(idx)) - k(idx) = 0.0; + + idx = find (m == -Inf); + if (!isempty (idx)) + k(idx) = 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 = sqrt(1.0-m(idx)); - c = sqrt(m(idx)); + 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 = sqrt (1 - m(idx)); + c = 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); + 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 + if (all (c./a < eps)) break; endif endfor - if n >= Nmax, error("ellipke: not enough workspace"); endif + if (n >= Nmax) error ("ellipke: not enough workspace"); endif k(idx) = 0.5*pi./a; - e(idx) = 0.5*pi.*(1.0-sum)./a; + e(idx) = 0.5*pi.*(1 - sum)./a; k(idx_neg) = mult_k.*k(idx_neg); e(idx_neg) = mult_e.*e(idx_neg); endif endfunction -## Test complete elliptic functions of first and second kind -## against "exact" solution from Mathematica 3.0 +%% Test complete elliptic functions of first and second kind +%% against "exact" solution from Mathematica 3.0 %!test %! m = [0.0; 0.01; 0.1; 0.5; 0.9; 0.99; 1.0 ]; -%! [k,e] = ellipke(m); +%! [k,e] = ellipke (m); %! %! # K(1.0) is really infinity - see below -%! k_exp = [ -%! 1.5707963267948966192; -%! 1.5747455615173559527; -%! 1.6124413487202193982; -%! 1.8540746773013719184; -%! 2.5780921133481731882; -%! 3.6956373629898746778; -%! 0.0 ]; -%! e_exp = [ -%! 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_exp,k,8*eps); -%! assert(e,e_exp,8*eps); +%! k_exp = [1.5707963267948966192; +%! 1.5747455615173559527; +%! 1.6124413487202193982; +%! 1.8540746773013719184; +%! 2.5780921133481731882; +%! 3.6956373629898746778; +%! 0.0 ]; +%! e_exp = [1.5707963267948966192; +%! 1.5668619420216682912; +%! 1.5307576368977632025; +%! 1.3506438810476755025; +%! 1.1047747327040733261; +%! 1.0159935450252239356; +%! 1.0 ]; +%! if k(7)==Inf, k(7)=0; endif; +%! assert (k, k_exp, 8*eps); +%! assert (e, e_exp, 8*eps); %% Test against A&S Table 17.1 %!test diff -r 1a3bfb14b5da -r f423873d3275 scripts/specfun/expint.m --- a/scripts/specfun/expint.m Sun Apr 28 17:08:11 2013 -0400 +++ b/scripts/specfun/expint.m Sun Apr 28 18:50:51 2013 -0400 @@ -1,20 +1,23 @@ ## Copyright (C) 2006 Sylvain Pelissier ## -## 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 3 of the License, or (at your option) any later -## version. +## This file is part of Octave. +## +## Octave 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 3 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. +## Octave 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 . +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## . ## -*- texinfo -*- -## @deftypefn {Function File} {@var{y} =} expint (@var{x}) +## @deftypefn {Function File} {} expint (@var{x}) ## Compute the exponential integral, ## @verbatim ## infinity @@ -23,14 +26,16 @@ ## / ## x ## @end verbatim -## @seealso{expint_E1, expint_Ei} ## @end deftypefn -function y = expint(x) +function y = expint (x) + if (nargin != 1) - print_usage; + print_usage (); endif - y = expint_E1(x); + + y = expint_E1 (x); + endfunction ## -*- texinfo -*- @@ -43,18 +48,28 @@ ## / ## x ## @end verbatim -## @seealso{expint, expint_Ei} ## @end deftypefn -function y = expint_E1(x) +function y = expint_E1 (x) + if (nargin != 1) - print_usage; + print_usage (); endif + y = x; - y(imag(x) > 0 & imag(x) != 0) = -expint_Ei(-y(imag(x) > 0 & imag(x) != 0)) -i.*pi; - y(imag(x) < 0 & imag(x) != 0) = -expint_Ei(-y(imag(x) < 0 & imag(x) != 0)) +i.*pi; - y(real(x) >= 0 & imag(x)==0) = -expint_Ei(-y(real(x) >= 0 & imag(x)==0)); - y(real(x) < 0 & imag(x)==0) = -expint_Ei(-y(real(x) < 0 & imag(x)==0)) -i.*pi; + + idx = (imag (x) > 0 & imag (x) != 0); + y(idx) = -expint_Ei (-y(idx)) - i.*pi; + + idx = (imag (x) < 0 & imag (x) != 0); + y(idx) = -expint_Ei (-y(idx)) + i.*pi; + + idx = (real (x) >= 0 & imag (x) == 0); + y(idx) = -expint_Ei (-y(idx)); + + idx = (real (x) < 0 & imag (x) == 0); + y(idx) = -expint_Ei (-y(idx)) - i.*pi; + endfunction ## -*- texinfo -*- @@ -67,39 +82,43 @@ ## / ## -x ## @end verbatim -## @seealso{expint, expint_E1} ## @end deftypefn -function y = expint_Ei(x) +function y = expint_Ei (x) + if (nargin != 1) - print_usage; + print_usage (); endif - y = zeros(size(x)); - F = @(x) exp(-x)./x; - s = prod(size(x)); + + y = zeros (size (x)); + F = @(x) exp (-x)./x; + s = prod (size (x)); + for t = 1:s; - if(x(t)<0 && imag(x(t)) == 0) - y(t) = -quad(F,-x(t),Inf); + if (x(t) < 0 && imag (x(t)) == 0) + y(t) = -quad (F, -x(t), Inf); else - if(abs(x(t)) > 2 && imag(x(t)) == 0) - y(t) = expint_Ei(2) - quad(F,-x(t),-2); + if (abs (x(t)) > 2 && imag (x(t)) == 0) + y(t) = expint_Ei (2) - quad (F, -x(t), -2); else - if(abs(x(t)) >= 10) - if(imag(x(t)) <= 0) + if (abs (x(t)) >= 10) + if (imag (x(t)) <= 0) a1 = 4.03640; a2 = 1.15198; b1 = 5.03637; b2 = 4.19160; - y(t) = -(x(t).^2 - a1.*x(t) + a2)./((x(t).^2-b1.*x(t)+b2).*(-x(t)).*exp(-x(t)))-i.*pi; + y(t) = -(x(t).^2 - a1.*x(t) + a2) ... + ./ ((x(t).^2 - b1.*x(t) + b2) .* (-x(t)) .* exp (-x(t))) ... + - i.*pi; else - y(t) = conj(expint_Ei(conj(x(t)))); + y(t) = conj (expint_Ei (conj (x(t)))); endif; ## Serie Expansion else for k = 1:100; - y(t) = y(t) + x(t).^k./(k.*factorial(k)); + y(t) = y(t) + x(t).^k ./ (k.*factorial (k)); endfor - y(t) = 0.577215664901532860606512090082402431 + log(x(t)) + y(t); + y(t) = 0.577215664901532860606512090082402431 + log (x(t)) + y(t); endif endif endif