--- a/main/specfun/NEWS	Wed Dec 05 01:51:51 2012 +0000
+++ b/main/specfun/NEWS	Wed Dec 05 16:38:22 2012 +0000
@@ -3,7 +3,17 @@
 
  ** The following functions are new:
 
-      big_factorial  lfactorial
+      big_factorial
+      lfactorial
+      zernike_cartesian
+      zernike_name
+      zernike_noll_to_mn
+      zernike_polar
+      zernike_R_poly
+
+ ** A group functions to deal with Zernike polynomials have been implemented.
+    These are a sequence of polynomials orthogonal on the unit disk. They
+    are used, for example, to characterize distortions in optical systems.
 
  ** The following functions have removed since they now are part of
     Octave core:

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/specfun/inst/zernike_R_poly.m	Wed Dec 05 16:38:22 2012 +0000
@@ -0,0 +1,57 @@
+## Copyright (C) 2012 Andreas Weber <andy.weber.aw@gmail.com>
+##
+## 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.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{R} =} zernike_R_poly (@var{m}, @var{n})
+## Return the first part of the radial zernike polynom R^m_n.
+##
+## The polynom returned has a length of N+1.
+## @seealso{zernike_cartesian, zernike_name, zernike_noll_to_nm, zernike_polar}
+## @end deftypefn
+
+function ret = zernike_R_poly (m, n)
+  if (nargin != 2)
+    print_usage ();
+  elseif (! isscalar (m) || m < 0 || m != fix (m) ||
+          ! isscalar (n) || n < 0 || n != fix (n))
+    error ("zernike_R_poly: M and N must all be non-negative integers");
+  endif
+  ret = zeros (1, n+1);
+  if (! mod (n-m, 2))
+    #TODO: try to omit for-loop
+    for k = 0:(n-m)/2
+      ret(2*k+1) = (-1)^k * bincoeff (n-k, k) * bincoeff (n-2*k, (n-m)/2-k);
+    endfor
+  endif
+endfunction
+
+## see http://en.wikipedia.org/wiki/Zernike_polynomials#Radial_polynomials
+## added all examples up to order 6
+%!assert(zernike_R_poly(0,0),[1])
+%!assert(zernike_R_poly(1,1),[1 0])
+%!assert(zernike_R_poly(0,2),[2 0 -1])
+%!assert(zernike_R_poly(2,2),[1 0 0])
+%!assert(zernike_R_poly(1,3),[3 0 -2 0])
+%!assert(zernike_R_poly(3,3),[1 0 0 0])
+%!assert(zernike_R_poly(0,4),[6 0 -6 0 1])
+%!assert(zernike_R_poly(2,4),[4 0 -3 0 0])
+%!assert(zernike_R_poly(4,4),[1 0 0 0 0])
+%!assert(zernike_R_poly(1,5),[10 0 -12 0 3 0])
+%!assert(zernike_R_poly(3,5),[5 0 -4 0 0 0])
+%!assert(zernike_R_poly(5,5),[1 0  0 0 0 0])
+%!assert(zernike_R_poly(0,6),[20 0 -30 0 12 0 -1])
+%!assert(zernike_R_poly(2,6),[15 0 -20 0  6 0  0])
+%!assert(zernike_R_poly(4,6),[ 6 0  -5 0  0 0  0])
+%!assert(zernike_R_poly(6,6),[ 1 0   0 0  0 0  0])

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/specfun/inst/zernike_cartesian.m	Wed Dec 05 16:38:22 2012 +0000
@@ -0,0 +1,61 @@
+## Copyright (C) 2012 Andreas Weber <andy.weber.aw@gmail.com>
+##
+## 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.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {@var{Z} =} zernike_cartesian (@var{x}, @var{y}, @var{n})
+## @deftypefnx {Function File} {@var{Z} =} zernike_cartesian (@var{x}, @var{y}, @var{n}, @var{limit_r})
+## Return the cartesian zernikes up to order n (as noll's index).
+##
+## If @var{limit_r} is false (default true), the polynoms for r>1 are @emph{not} set to NaN
+## because strictly, the polynoms are only defined for 0 <= r <= 1.
+## 
+## Size of @var{x} must be equal size of @var{y}.
+##
+## Demo: type "demo zernike_cartesian"
+##
+## @seealso{zernike_name, zernike_noll_to_nm, zernike_polar, zernike_R_poly}
+## @end deftypefn
+
+## TODO: The cartesian zernikes can be calculated quicker for example using a method
+## described in "Hedser van Brug: Efficient Cartesian representation of Zernike polynomials."
+## Until then the cartesians get mapped to the polar ones.
+
+function Z = zernike_cartesian (x, y, n, limit_r = true)
+  if (nargin < 3 || nargin > 4)
+    print_usage ();
+  elseif (! isscalar (n) || n < 1 || n != fix (n))
+    error ("zernike_cartesian: N must be a integer >=1");
+  elseif (ndims (x) != ndims (y) || any (size (x) != size (y)))
+    error ("zernike_cartesian: X and Y must have the same size")
+  endif
+  r   = sqrt (x.*x + y.*y);
+  phi = atan2 (y, x);
+  Z   = zernike_polar (r, phi ,n ,limit_r);
+endfunction
+
+%!demo
+%! t = linspace (-1, 1, 150);
+%! [x, y] = meshgrid (t, t);
+%! max_order = 16;
+%! Z = zernike_cartesian (x, y, max_order);
+%! for k = 1:max_order
+%!   subplot (4, 4, k)
+%!   factors = zeros (max_order, 1);
+%!   factors(k) = 1;
+%!   z = reshape (Z*factors, size (x));
+%!   imagesc (z)
+%!   axis ("off", "equal")
+%!   title (zernike_name (k))
+%! endfor

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/specfun/inst/zernike_name.m	Wed Dec 05 16:38:22 2012 +0000
@@ -0,0 +1,70 @@
+## Copyright (C) 2012 Andreas Weber <andy.weber.aw@gmail.com>
+##
+## 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.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{name} =} zernike_name (@var{n})
+## Return the classic name for noll's index @var{n} or
+## "-" (no name defined) without warning if @var{n} > 21.
+##
+## Examples:
+## @example
+## @group
+## zernike_name(4)
+##     @result{} defocus
+## zernike_name(21)
+##     @result{} vertical pentafoil
+## @end group
+## @end example
+##
+## @seealso{zernike_cartesian, zernike_noll_to_nm, zernike_polar, zernike_R_poly}
+## @end deftypefn
+
+function name = zernike_name (n)
+  ## taken from http://www.telescope-optics.net/zernike_coefficients.htm
+  persistent classical_names = ...
+  {
+    "piston";
+    "horizontal tilt";
+    "vertical tilt";
+    "defocus";
+    "oblique primary astigmatism";
+    "vertical primary astigmatism";
+    "vertical coma";
+    "horizontal coma";
+    "vertical trefoil";
+    "oblique trefoil";
+    "primary spherical";
+    "vertical secondary astigmatism";
+    "oblique secondary astigmatism";
+    "vertical quadrafoil";
+    "oblique quadrafoil";
+    "horizontal secondary coma";
+    "vertical secondary coma";
+    "oblique secondary trefoil";
+    "vertical secondary trefoil";
+    "oblique pentafoil";
+    "vertical pentafoil";
+  };
+  if (nargin != 1)
+    print_usage ();
+  elseif (! isscalar (n) || n < 1 || n != fix (n))
+    error ("zernike_name: n must be a integer >=1");
+  endif
+  if (n > numel (classical_names))
+    name = "-";
+  else
+    name = classical_names{n};
+  endif
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/specfun/inst/zernike_noll_to_mn.m	Wed Dec 05 16:38:22 2012 +0000
@@ -0,0 +1,77 @@
+## Copyright (C) 2012 Andreas Weber <andy.weber.aw@gmail.com>
+##
+## 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.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{m}, @var{n}] =} zernike_noll_to_mn (@var{j})
+## Convert Noll's index @var{j} to @var{m} (Azimuthal degree) and @var{n} (Radial degree).
+##
+## See sequence A176988 in OEIS (http://oeis.org/A176988)
+## @seealso{zernike_cartesian, zernike_name, zernike_polar, zernike_R_poly}
+## @end deftypefn
+
+function [m, n] = zernike_noll_to_mn (j)
+  if (nargin != 1 || nargout != 2)
+    print_usage ();
+  elseif (any (j(:) < 1 | j(:) != fix (j(:))) || ! isvector (j))
+    error ("zernike_noll_to_mn: j has to be a vector with integers >=1");
+  endif
+  n  = fix (sqrt (2*j-1) + 0.5) - 1;
+  s  = mod (n, 2);
+  me = 2 * fix ((2*j + 1 - n.*(n+1)) / 4);        #the even ones
+  mo = 2 * fix ((2*(j+1) - n.*(n+1)) / 4) - 1;    #the odd ones
+  m  = (mo.*s + me.*(1-s)).*(1 - 2*mod(j,2));
+endfunction
+
+## see http://oeis.org/A176988
+##
+%!test
+%! [m,n]=zernike_noll_to_mn(1);
+%! assert([m n],[0 0])
+%!test
+%! [m,n]=zernike_noll_to_mn(2);
+%! assert([m n],[1 1]) 
+%!test
+%! [m,n]=zernike_noll_to_mn(3);
+%! assert([m n],[-1 1]) 
+%!test
+%! [m,n]=zernike_noll_to_mn(4);
+%! assert([m n],[0 2]) 
+%!test
+%! [m,n]=zernike_noll_to_mn(5);
+%! assert([m n],[-2 2]) 
+
+## skipp noll indices 6..19 TODO: or should we add all to this test?
+
+%!test
+%! [m,n]=zernike_noll_to_mn(20);
+%! assert([m n],[5 5]) 
+
+## skipp noll indices 21..33 TODO: or should we add all to this test?
+
+%!test
+%! [m,n]=zernike_noll_to_mn(34);
+%! assert([m n],[5 7]) 
+%!test
+%! [m,n]=zernike_noll_to_mn(35);
+%! assert([m n],[-7 7]) 
+%!test
+%! [m,n]=zernike_noll_to_mn(36);
+%! assert([m n],[7 7]) 
+
+## vector test
+%!test
+%! [m,n]=zernike_noll_to_mn([2 5 8 35]);
+%! assert(m,[1 -2 1 -7])
+%! assert(n,[1 2 3 7])

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/specfun/inst/zernike_polar.m	Wed Dec 05 16:38:22 2012 +0000
@@ -0,0 +1,66 @@
+## Copyright (C) 2012 Andreas Weber <andy.weber.aw@gmail.com>
+##
+## 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.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {@var{Z} =} zernike_polar (@var{r}, @var{phi}, @var{n})
+## @deftypefnx {Function File} {@var{Z} =} zernike_polar (@var{r}, @var{phi}, @var{n}, @var{limit_r})
+## Return the polar zernikes up to order n (as noll's index).
+##
+## If @var{limit_r} is false (default true), the polynoms for r>1 are @emph{not} set to NaN 
+## because strictly, the polynoms are only defined for 0 <= r <= 1.
+## 
+## The first argument @var{r} is a matrix containing the radial distance,
+## the second argument @var{phi} a matrix with the angles.
+##
+## Size of @var{r} must be equal size of @var{phi}.
+##
+## This file hasn't a demo yet but have a look on "demo zernike_cartesian"
+## @seealso{zernike_cartesian, zernike_name, zernike_noll_to_nm, zernike_R_poly}
+## @end deftypefn
+
+function Z = zernike_polar (r, phi, n, limit_r = true)
+  if (nargin < 3 || nargin > 4)
+    print_usage ();
+  elseif (! isscalar (n) || n < 1 || n != fix (n))
+    error ("zernike_polar: n must be a integer >=1");
+  elseif (any (size (r) != size (phi)))
+    error ("zernike_polar: R and PHI must have the same size")
+  endif
+  Z = zeros (numel (r), n);
+  for k = 1:n
+    [m, n] = zernike_noll_to_mn (k);
+    P      = zernike_R_poly (abs (m), n);
+    r_eval = polyval (P, r(:));
+    if (m < 0)
+      r_eval = r_eval .* sin (abs (m) * phi(:));
+    else
+      r_eval = r_eval .* cos (m * phi(:));
+    endif
+    Z(:, k) = r_eval;
+  endfor
+  if (limit_r)
+    Z(r(:) > 1, :) = NaN;
+  endif
+endfunction
+
+%!test
+%! r=0.2; phi=1.23; n=4;
+%! ret=zernike_polar(r,phi,n);
+%! assert(ret,[1 r*cos(phi) r*sin(phi) 2*r^2-1],5*eps)
+
+%!test
+%! r=[0.5 0.8]; phi=[pi/4 pi/4]; n=4;
+%! ret=zernike_polar(r,phi,n);
+%! assert(ret,[1 r(1)*cos(phi(1)) r(1)*sin(phi(1)) 2*r(1)^2-1; 1 r(2)*cos(phi(2)) r(2)*sin(phi(2)) 2*r(2)^2-1],5*eps)