new file mode 100644
--- /dev/null
+++ b/scripts/plot/surfl.m
@@ -0,0 +1,178 @@
+## Copyright (C) 2009 Kai Habel
+##
+## 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.
+##
+## 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 Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} surfl (@var{x}, @var{y}, @var{z})
+## @deftypefnx {Function File} {} surfl (@var{z})
+## @deftypefnx {Function File} {} surfl (@var{x}, @var{y}, @var{z}, @var{L})
+## @deftypefnx {Function File} {} surfl (@var{x}, @var{y}, @var{z}, @var{L}, @var{P})
+## @deftypefnx {Function File} {} surfl (...,"light")
+## Plot a lighted surface given matrices @var{x}, and @var{y} from @code{meshgrid} and
+## a matrix @var{z} corresponding to the @var{x} and @var{y} coordinates of
+## the mesh.  If @var{x} and @var{y} are vectors, then a typical vertex
+## is (@var{x}(j), @var{y}(i), @var{z}(i,j)).  Thus, columns of @var{z}
+## correspond to different @var{x} values and rows of @var{z} correspond
+## to different @var{y} values.
+##
+## The light direction can be specified using @var{L}. It can be
+## given as 2-element vector [azimuth, elevation] in degrees or as 3-element vector [lx, ly, lz].
+## The default value is rotated 45° counter-clockwise from the current view.
+##
+## The material properties of the surface can specified using a 4-element vector
+## @var{P} = [@var{AM} @var{D} @var{SP} @var{exp}] which defaults to
+## @var{p} = [0.55 0.6 0.4 10]. 
+## @table @code
+## @item "AM" strength of ambient light
+## @item "D" strength of diffuse reflection
+## @item "SP" strength of specular reflection
+## @item "EXP" specular exponent
+## @end table
+## 
+## The default lighting mode "cdata", changes the cdata property to give the impression
+## of a lighted surface. Please note: the alternative "light" mode, which creates a light
+## object to iluminate the the surface is not implemented (yet).
+##
+## Example:
+##
+## @example
+## colormap(bone);
+## surfl(peaks);
+## shading interp;
+## @end example
+## @seealso{surf, diffuse, specular, surface}
+## @end deftypefn
+
+## Author: Kai Habel <kai.habel@gmx.de>
+
+function retval = surfl (varargin)
+
+  [h, varargin] = __plt_get_axis_arg__ ("surfl", varargin{:});
+
+  oldh = gca ();
+  unwind_protect
+    axes (h);
+    newplot ();
+
+    ## check for lighting type
+    use_cdata = true;
+    if (ischar (varargin{end}))
+      lstr = varargin{end};
+      if strncmp (tolower (lstr), "light", 5)
+        warning ("light method not supported (yet), using cdata method instead");
+        # this can be implemented when light objects are being
+        # supported.
+        use_cdata = false;
+      elseif (strncmp (tolower (lstr), "cdata", 5))
+        use_cdata = true;
+      else
+        usage ("unknown lighting method");
+      endif
+      varargin(end) = [];
+    endif
+
+    ## check for reflection properties argument
+    ## r = [ambient light strength,
+    ##      diffuse reflection strength,
+    ##      specular reflection strength,
+    ##      specular shine] 
+    if ((length (varargin{end}) == 4) && isnumeric (varargin{end}))
+      r = varargin{end};
+      varargin(end) = [];
+    else
+      ## default values
+      r = [0.55, 0.6, 0.4, 10];
+    endif
+    
+
+    ## check for light vector (lv) argument
+    have_lv = false;
+    if (isnumeric (varargin{end}))
+      len = numel (varargin{end});
+      lastarg = varargin{end};
+      if (len == 3)
+        lv = lastarg;
+        varargin(end) = [];
+        have_lv = true;
+      elseif (len == 2)
+        [lv(1), lv(2), lv(3)] = sph2cart ((lastarg(1) - 90) * pi/180, lastarg(2) * pi/180, 1.0);
+        varargin(end) = [];
+        have_lv = true;
+      endif
+    endif
+    
+    tmp = surface (varargin{:});
+    if (! ishold ())
+      set (h, "view", [-37.5, 30],
+	   "xgrid", "on", "ygrid", "on", "zgrid", "on", "clim", [0 1]);
+    endif
+
+    ## get view vector (vv)
+    a = axis;
+    [az, el] = view;
+    [vv(1), vv(2), vv(3)] = sph2cart ((az - 90) * pi/180.0, el * pi/180.0, 1.0);
+    vv /= norm (vv);
+
+    if (!have_lv)
+      ## calculate light vector (lv) from view vector
+      Phi = 45.0 / 180.0 * pi;
+      R = [cos(Phi), -sin(Phi), 0;\
+           sin(Phi),  cos(Phi), 0;\
+           0,          0,         1];
+      lv = (R * vv.').';
+    endif
+
+    vn = get (tmp, "vertexnormals");
+    dar = get (h, "dataaspectratio");
+    vn(:, :, 1) *= dar(1);
+    vn(:, :, 2) *= dar(2);
+    vn(:, :, 3) *= dar(3);
+
+    ## normalize vn
+    vn = vn ./ repmat (sqrt (sumsq (vn, 3)), [1, 1, 3]);
+    [nr, nc] = size(get(tmp, "zdata"));
+
+    ## ambient, diffuse, and specular term
+    cdata = r(1) * ones (nr, nc) \
+          + r(2) * diffuse  (vn(:, :, 1), vn(:, :, 2), vn(:, :, 3), lv) \
+          + r(3) * specular (vn(:, :, 1), vn(:, :, 2), vn(:, :, 3), lv, vv, r(4));
+
+    set (tmp, "cdata", cdata ./ sum (r(1:3)));
+    
+  unwind_protect_cleanup
+    axes (oldh);
+  end_unwind_protect
+
+  if (nargout > 0)
+    retval = tmp;
+  endif
+
+endfunction
+
+%!demo
+%! [X,Y,Z]=sombrero;
+%! colormap(copper);
+%! surfl(X,Y,Z);
+%! shading interp;
+
+%!demo
+%! [X,Y,Z]=sombrero;
+%! colormap(copper);
+%! [az, el] = view;
+%! surfl(X,Y,Z,[az+225,el],[0.2 0.6 0.4 25]);
+%! shading interp;
+