--- a/src/data.cc
+++ b/src/data.cc
@@ -74,6 +74,148 @@
   return retval;
 }
 
+// These mapping functions may also be useful in other places, eh?
+
+typedef double (*d_dd_fcn) (double, double);
+
+static Matrix
+map (d_dd_fcn f, double x, const Matrix& y)
+{
+  int nr = y.rows ();
+  int nc = y.columns ();
+
+  Matrix retval (nr, nc);
+
+  for (int j = 0; j < nc; j++)
+    for (int i = 0; i < nr; i++)
+      retval.elem (i, j) = f (x, y.elem (i, j));
+
+  return retval;
+}
+
+static Matrix
+map (d_dd_fcn f, const Matrix& x, double y)
+{
+  int nr = x.rows ();
+  int nc = x.columns ();
+
+  Matrix retval (nr, nc);
+
+  for (int j = 0; j < nc; j++)
+    for (int i = 0; i < nr; i++)
+      retval.elem (i, j) = f (x.elem (i, j), y);
+
+  return retval;
+}
+
+static Matrix
+map (d_dd_fcn f, const Matrix& x, const Matrix& y)
+{
+  int x_nr = x.rows ();
+  int x_nc = x.columns ();
+
+  int y_nr = y.rows ();
+  int y_nc = y.columns ();
+
+  assert (x_nr == x_nc && y_nr == y_nc);
+
+  Matrix retval (x_nr, x_nc);
+
+  for (int j = 0; j < x_nc; j++)
+    for (int i = 0; i < x_nr; i++)
+      retval.elem (i, j) = f (x.elem (i, j), y.elem (i, j));
+
+  return retval;
+}
+
+DEFUN ("atan2", Fatan2, Satan2, 3, 1,
+  "atan2 (Y, X): atan (Y / X) in range -pi to pi")
+{
+  Octave_object retval;
+
+  if (args.length () != 3)
+    print_usage ("atan2");
+  else
+    {
+      tree_constant arg_y = args(1);
+      tree_constant arg_x = args(2);
+
+      int y_nr = arg_y.rows ();
+      int y_nc = arg_y.columns ();
+
+      int x_nr = arg_x.rows ();
+      int x_nc = arg_x.columns ();
+
+      int arg_y_empty = empty_arg ("atan2", y_nr, y_nc);
+      int arg_x_empty = empty_arg ("atan2", x_nr, x_nc);
+
+      if (arg_y_empty < 0 || arg_x_empty < 0)
+	return retval;
+
+      if (arg_y_empty || arg_x_empty)
+	{
+	  retval = Matrix ();
+	  return retval;
+	}
+
+      int y_is_scalar = (y_nr == 1 && y_nc == 1);
+      int x_is_scalar = (x_nr == 1 && x_nc == 1);
+
+      if (y_is_scalar && x_is_scalar)
+	{
+	  double y = arg_y.double_value ();
+
+	  if (! error_state)
+	    {
+	      double x = arg_x.double_value ();
+
+	      if (! error_state)
+		retval = atan2 (y, x);
+	    }
+	}
+      else if (y_is_scalar)
+	{
+	  double y = arg_y.double_value ();
+
+	  if (! error_state)
+	    {
+	      Matrix x = arg_x.matrix_value ();
+
+	      if (! error_state)
+		retval = map (atan2, y, x);
+	    }
+	}
+      else if (x_is_scalar)
+	{
+	  Matrix y = arg_y.matrix_value ();
+
+	  if (! error_state)
+	    {
+	      double x = arg_x.double_value ();
+
+	      if (! error_state)
+		retval = map (atan2, y, x);
+	    }
+	}
+      else if (y_nr == x_nr && y_nc == x_nc)
+	{
+	  Matrix y = arg_y.matrix_value ();
+
+	  if (! error_state)
+	    {
+	      Matrix x = arg_x.matrix_value ();
+
+	      if (! error_state)
+		retval = map (atan2, y, x);
+	    }
+	}
+      else
+	error ("atan2: nonconformant matrices");
+    }
+
+  return retval;
+}
+
 DEFUN ("cumprod", Fcumprod, Scumprod, 2, 1,
   "cumprod (X): cumulative products")
 {