Mercurial > octave
changeset 19099:55b613e5183d
legendre.m: Overhaul function.
* legendre.m: Rewrite docstring. Add input validation for real value of N.
Add input validation for X. Shorten error() messages. Use in-place operators
for performance. Add %!error tests for input validation.
author | Rik <rik@octave.org> |
---|---|
date | Fri, 19 Sep 2014 14:30:37 -0700 |
parents | 4591a1238ee0 |
children | c5a0b995b996 |
files | scripts/specfun/legendre.m |
diffstat | 1 files changed, 42 insertions(+), 38 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/specfun/legendre.m Fri Sep 19 14:25:21 2014 -0700 +++ b/scripts/specfun/legendre.m Fri Sep 19 14:30:37 2014 -0700 @@ -21,15 +21,19 @@ ## @deftypefn {Function File} {@var{l} =} legendre (@var{n}, @var{x}) ## @deftypefnx {Function File} {@var{l} =} legendre (@var{n}, @var{x}, @var{normalization}) ## Compute the Legendre function of degree @var{n} and order -## @var{m} = 0 @dots{} N@. The optional argument, @var{normalization}, -## may be one of @qcode{"unnorm"}, @qcode{"sch"}, or @qcode{"norm"}. -## The default is @qcode{"unnorm"}. The value of @var{n} must be a -## non-negative scalar integer. +## @var{m} = 0 @dots{} @var{n}. +## +## The value @var{n} must be a real non-negative integer. @var{x} is a +## vector with real-valued elements in the range [-1, 1]. ## -## If the optional argument @var{normalization} is missing or is -## @qcode{"unnorm"}, compute the Legendre function of degree @var{n} and -## order @var{m} and return all values for @var{m} = 0 @dots{} @var{n}. -## The return value has one dimension more than @var{x}. +## The optional argument @var{normalization} may be one of @qcode{"unnorm"}, +## @qcode{"sch"}, or @qcode{"norm"}. The default if no normalization is given +## is @qcode{"unnorm"}. +## +## When the optional argument @var{normalization} is @qcode{"unnorm"}, compute +## the Legendre function of degree @var{n} and order @var{m} and return all +## values for @var{m} = 0 @dots{} @var{n}. The return value has one dimension +## more than @var{x}. ## ## The Legendre Function of degree @var{n} and order @var{m}: ## @@ -42,7 +46,7 @@ ## ## @example ## @group -## m m 2 m/2 d^m +## m m 2 m/2 d^m ## P(x) = (-1) * (1-x ) * ---- P(x) ## n dx^m n ## @end group @@ -84,12 +88,12 @@ ## @end group ## @end example ## -## If the optional argument @code{normalization} is @qcode{"sch"}, +## When the optional argument @code{normalization} is @qcode{"sch"}, ## compute the Schmidt semi-normalized associated Legendre function. ## The Schmidt semi-normalized associated Legendre function is related ## to the unnormalized Legendre functions by the following: ## -## For Legendre functions of degree n and order 0: +## For Legendre functions of degree @var{n} and order 0: ## ## @tex ## $$ @@ -127,7 +131,7 @@ ## ## @end ifnottex ## -## If the optional argument @var{normalization} is @qcode{"norm"}, +## When the optional argument @var{normalization} is @qcode{"norm"}, ## compute the fully normalized associated Legendre function. ## The fully normalized associated Legendre function is related ## to the unnormalized Legendre functions by the following: @@ -162,11 +166,9 @@ print_usage (); endif - if (!isscalar (n) || n < 0 || n != fix (n)) - error ("legendre: N must be a non-negative scalar integer"); - endif - - if (!isreal (x) || any (x(:) < -1 | x(:) > 1)) + if (! isreal (n) || ! isscalar (n) || n < 0 || n != fix (n)) + error ("legendre: N must be a real non-negative integer"); + elseif (! isreal (x) || any (x(:) < -1 | x(:) > 1)) error ("legendre: X must be real-valued vector in the range -1 <= X <= 1"); endif @@ -176,15 +178,17 @@ normalization = "unnorm"; endif + unnorm = false; switch (normalization) + case "unnorm" + scale = 1; + unnorm = true; case "norm" scale = sqrt (n+0.5); case "sch" scale = sqrt (2); - case "unnorm" - scale = 1; otherwise - error ('legendre: expecting NORMALIZATION option to be "norm", "sch", or "unnorm"'); + error ('legendre: NORMALIZATION option must be "unnorm", "norm", or "sch"'); endswitch scale = scale * ones (size (x)); @@ -204,11 +208,11 @@ for k = m+1:n lpm3a = (2*k-1) .* x .* lpm2; lpm3b = (k+m-2) .* lpm1; - lpm3 = (lpm3a - lpm3b)/(k-m+1); + lpm3 = (lpm3a - lpm3b) / (k-m+1); lpm1 = lpm2; lpm2 = lpm3; if (! warned_overflow) - if (any (abs (lpm3a) > realmax) + if ( any (abs (lpm3a) > realmax) || any (abs (lpm3b) > realmax) || any (abs (lpm3) > realmax)) overflow = true; @@ -216,24 +220,23 @@ endif endfor retval(m,:) = lpm3(:); - if (strcmp (normalization, "unnorm")) - scale = -scale * (2*m-1); - else - ## normalization == "sch" or normalization == "norm" - scale = scale / sqrt ((n-m+1)*(n+m))*(2*m-1); + if (unnorm) + scale *= -(2*m-1); + else # normalization = "sch" or "norm" + scale *= (2*m-1) / sqrt ((n-m+1)*(n+m)); endif - scale = scale .* sqrt (1-x.^2); + scale .*= sqrt (1-x.^2); endfor retval(n+1,:) = scale(:); if (isvector (x)) - ## vector case is special + ## vector case is special retval = reshape (retval, n + 1, length (x)); endif if (strcmp (normalization, "sch")) - retval(1,:) = retval(1,:) / sqrt (2); + retval(1,:) ./= sqrt (2); endif if (overflow && ! warned_overflow) @@ -289,8 +292,8 @@ %! assert (result, full (ones (1,11))); %!test +%! ## Test matrix input %! result = legendre (3, [-1,0,1;1,0,-1]); -%! ## Test matrix input %! expected(:,:,1) = [-1,1;0,0;0,0;0,0]; %! expected(:,:,2) = [0,0;1.5,1.5;0,0;-15,-15]; %! expected(:,:,3) = [1,-1;0,0;0,0;0,0]; @@ -306,11 +309,12 @@ %!error legendre () %!error legendre (1) %!error legendre (1,2,3,4) -%!error legendre ([1, 2], [-1, 0, 1]) -%!error legendre (-1, [-1, 0, 1]) -%!error legendre (1.1, [-1, 0, 1]) -%!error legendre (1, [-1+i, 0, 1]) -%!error legendre (1, [-2, 0, 1]) -%!error legendre (1, [-1, 0, 2]) -%!error legendre (1, [-1, 0, 1], "badnorm") +%!error <must be a real non-negative integer> legendre (i, [-1, 0, 1]) +%!error <must be a real non-negative integer> legendre ([1, 2], [-1, 0, 1]) +%!error <must be a real non-negative integer> legendre (-1, [-1, 0, 1]) +%!error <must be a real non-negative integer> legendre (1.1, [-1, 0, 1]) +%!error <must be real-valued vector> legendre (1, [-1+i, 0, 1]) +%!error <in the range -1 .= X .= 1> legendre (1, [-2, 0, 1]) +%!error <in the range -1 .= X .= 1> legendre (1, [-1, 0, 2]) +%!error <NORMALIZATION option must be> legendre (1, [-1, 0, 1], "badnorm")