Mercurial > hg > octave-lojdl
changeset 12893:72ffa81a68d4 stable
legendre.m: Allow ND-array inputs (Bug #33526).
* legendre.m: Allow ND-array inputs (Bug #33526).
author | Marco Caliari <marco.caliari@univr.it> |
---|---|
date | Wed, 27 Jul 2011 11:21:21 -0700 |
parents | 67bf9b30f3f9 |
children | ef5ebbf2a657 |
files | scripts/specfun/legendre.m |
diffstat | 1 files changed, 27 insertions(+), 9 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/specfun/legendre.m +++ b/scripts/specfun/legendre.m @@ -166,7 +166,7 @@ error ("legendre: N must be a non-negative scalar integer"); endif - if (!isvector (x) || !isreal (x) || any (x < -1 | x > 1)) + if (!isreal (x) || any (x(:) < -1 | x(:) > 1)) error ("legendre: X must be real-valued vector in the range -1 <= X <= 1"); endif @@ -187,10 +187,7 @@ error ('legendre: expecting NORMALIZATION option to be "norm", "sch", or "unnorm"'); endswitch - if (rows (x) != 1) - x = x'; - endif - scale = scale * ones (1, numel (x)); + scale = scale * ones (size (x)); ## Based on the recurrence relation below ## m m m @@ -199,6 +196,7 @@ ## http://en.wikipedia.org/wiki/Associated_Legendre_function overflow = false; + retval = zeros([n+1, size(x)]); for m = 1:n lpm1 = scale; lpm2 = (2*m-1) .* x .* scale; @@ -217,7 +215,7 @@ endif endif endfor - retval(m,:) = lpm3; + retval(m,:) = lpm3(:); if (strcmp (normalization, "unnorm")) scale = -scale * (2*m-1); else @@ -227,7 +225,12 @@ scale = scale .* sqrt(1-x.^2); endfor - retval(n+1,:) = scale; + retval(n+1,:) = scale(:); + + if (isvector (x)) + ## vector case is special + retval = reshape (retval, n + 1, length (x)); + endif if (strcmp (normalization, "sch")) retval(1,:) = retval(1,:) / sqrt (2); @@ -240,6 +243,7 @@ endfunction + %!test %! result = legendre (3, [-1.0 -0.9 -0.8]); %! expected = [ @@ -278,11 +282,25 @@ %!test %! result = legendre (150, 0); %! ## This agrees with Matlab's result. -%! assert (result(end), 3.7532741115719e+306, 0.0000000000001e+306) +%! assert (result(end), 3.7532741115719e+306, 0.0000000000001e+306); %!test %! result = legendre (0, 0:0.1:1); -%! assert (result, full(ones(1,11))) +%! assert (result, full(ones(1,11))); + +%!test +%! 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]; +%! assert (result, expected); + +%!test +%! result = legendre (3, [-1,0,1;1,0,-1]'); +%! expected(:,:,1) = [-1,0,1;0,1.5,0;0,0,0;0,-15,0]; +%! expected(:,:,2) = [1,0,-1;0,1.5,0;0,0,0;0,-15,0]; +%! assert (result, expected); %% Check correct invocation %!error legendre ();