--- a/scripts/ChangeLog	Sat Feb 14 19:50:43 2009 +0100
+++ b/scripts/ChangeLog	Sun Feb 15 19:07:05 2009 -0500
@@ -1,3 +1,7 @@
+2009-02-15  Ben Abbott <bpabbott@mac.com>
+
+	* plot/subplot.m: Compatible placement of subplots.
+
 2009-02-13 Ben Abbott <bpabott@mac.com>
 
 	* plot/__go_draw_axes__.m: Respect axes ticklength property.

--- a/scripts/plot/subplot.m	Sat Feb 14 19:50:43 2009 +0100
+++ b/scripts/plot/subplot.m	Sun Feb 15 19:07:05 2009 -0500
@@ -82,8 +82,10 @@
     tmp = (tmp - columns) / 10;
     rows = rem (tmp, 10);
 
-  elseif (! (isscalar (columns) && isscalar (rows) && isscalar (index)))
-    error ("subplot: columns, rows, and index have to be scalars");
+  elseif (! (isscalar (columns) && isscalar (rows)))
+    error ("subplot: columns, and rows must be scalars");
+  elseif (any (index < 1) || any (index > rows*columns))
+    error ("subplot: index value must be greater than 1 and less than rows*columns")
   endif
 
   columns = round (columns);
@@ -98,70 +100,163 @@
     error ("subplot: columns,rows,index must be be positive");
   endif
 
-  xsize = 1 / columns;
-  ysize = 1 / rows;
-
-  yp = fix ((index-1)/columns);
-  xp = index - yp*columns - 1;
+  units = get (0, "defaultaxesunits");
+  unwind_protect
+    set (0, "defaultaxesunits", "normalized")
+    pos = subplot_position (rows, columns, index, "outerposition", units);
 
-  x0 = xp * xsize;
-  y0 = (rows - yp - 1) * ysize;
-
-  pos = [x0, y0, xsize, ysize];
+    cf = gcf ();
 
-  x1 = x0 + xsize;
-  y1 = y0 + ysize;
-
-  cf = gcf ();
+    set (cf, "nextplot", "add");
 
-  set (cf, "nextplot", "add");
-
-  found = false;
-  kids = get (cf, "children");
-  for child = reshape (kids, 1, numel (kids))
-    ## Check whether this child is still valid; this might not be the
-    ## case anymore due to the deletion of previous children (due to
-    ## "deletefcn" callback or for legends/colorbars that are deleted
-    ## with their corresponding axes).
-    if (! ishandle (child))
-      continue;
-    endif
-    if (strcmp (get (child, "type"), "axes"))
-      ## Skip legend and colorbar objects.
-      if (strcmp (get (child, "tag"), "legend") || 
-	  strcmp (get (child, "tag"), "colorbar"))
+    found = false;
+    kids = get (cf, "children");
+    for child = reshape (kids, 1, numel (kids))
+      ## Check whether this child is still valid; this might not be the
+      ## case anymore due to the deletion of previous children (due to
+      ## "deletefcn" callback or for legends/colorbars that are deleted
+      ## with their corresponding axes).
+      if (! ishandle (child))
         continue;
       endif
-      objpos = get (child, "outerposition");
-      if (objpos == pos)
-	## If the new axes are in exactly the same position as an
-	## existing axes object, use the existing axes.
-	found = true;
-	tmp = child;
-      else
-	## If the new axes overlap an old axes object, delete the old
-	## axes.
-	objx0 = objpos(1);
-	objx1 = objx0 + objpos(3);
-	objy0 = objpos(2);
-	objy1 = objy0 + objpos(4);
-	if (! (x0 >= objx1 || x1 <= objx0 || y0 >= objy1 || y1 <= objy0))
-	  delete (child);
-	endif
+      if (strcmp (get (child, "type"), "axes"))
+        ## Skip legend and colorbar objects.
+        if (strcmp (get (child, "tag"), "legend") || 
+	    strcmp (get (child, "tag"), "colorbar"))
+          continue;
+        endif
+        objpos = get (child, "outerposition");
+        if (all (objpos == pos))
+	  ## If the new axes are in exactly the same position as an
+	  ## existing axes object, use the existing axes.
+	  found = true;
+	  tmp = child;
+        else
+	  ## If the new axes overlap an old axes object, delete the old
+	  ## axes.
+	  x0 = pos(1);
+	  x1 = x0 + pos(3);
+	  y0 = pos(2);
+	  y1 = y0 + pos(4);	
+	  objx0 = objpos(1);
+	  objx1 = objx0 + objpos(3);
+	  objy0 = objpos(2);
+	  objy1 = objy0 + objpos(4);
+	  if (! (x0 >= objx1 || x1 <= objx0 || y0 >= objy1 || y1 <= objy0))
+	    delete (child);
+	  endif
+        endif
       endif
-    endif
-  endfor
+    endfor
 
-  if (found)
-    set (cf, "currentaxes", tmp);
-  else
-    border = [0.130, 0.110, 0.225, 0.185] .* [xsize, ysize, xsize, ysize];
-    pos2 = [pos(1:2) + border(1:2), pos(3:4) - border(1:2) - border(3:4)];
-    tmp = axes ("outerposition", pos, "position", pos2);
-  endif
+    if (found)
+      set (cf, "currentaxes", tmp);
+    else
+      pos2 = subplot_position (rows, columns, index, "position", units);
+      tmp = axes ("outerposition", pos, "position", pos2);
+    endif
+
+  unwind_protect_cleanup
+    set (0, "defaultaxesunits", units);
+  end_unwind_protect
 
   if (nargout > 0)
     h = tmp;
   endif
 
 endfunction
+
+function pos = subplot_position (rows, columns, index, position_property, units)
+
+  ## For 1 row and 1 column return the usual default.
+  if (rows == 1 && columns == 1)
+    if (strcmpi (position_property, "position"))
+      pos = get (0, "defaultaxesposition");
+    else
+      pos = get (0, "defaultaxesouterposition");
+    endif
+    return
+  endif
+
+  ## This produces compatible behavior for the "position" property.
+  margins.left   = 0.130;
+  margins.right  = 0.095;
+  margins.top    = 0.075;
+  margins.bottom = 0.110;
+  pc = 1 ./ [0.1860, (margins.left + margins.right - 1)];
+  margins.column = 1 ./ polyval (pc , columns);
+  pr = 1 ./ [0.2282, (margins.top + margins.bottom - 1)];
+  margins.row    = 1 ./ polyval (pr , rows);
+
+  ## Calculate the width/height of the subplot axes.
+  width = 1 - margins.left - margins.right - (columns-1)*margins.column;
+  width = width / columns;
+  height = 1 - margins.top - margins.bottom - (rows-1)*margins.row;
+  height = height / rows;
+
+  if (strcmp (position_property, "outerposition") )
+    ## Calculate the outerposition/position inset
+    if (rows > 1)
+      inset.top    = 8/420;
+      inset.bottom = max (polyval ([0.1382,-0.0026], width), 16/420);
+    else
+      inset.bottom = margins.bottom;
+      inset.top = margins.top;
+    endif
+    if (columns > 1)
+      if (strcmpi (units, "normalized"))
+        inset.right = max (polyval ([0.1200,-0.0014], width), 5/560);
+      else
+        inset.right = max (polyval ([0.1252,-0.0023], width), 5/560);
+      endif
+      inset.left   = 22/560;
+    else
+      inset.left  = margins.left;
+      inset.right = margins.right;
+    endif
+    ## Apply the inset to the geometries for the "position" property.
+    margins.column = margins.column - inset.right - inset.left;
+    margins.row = margins.row - inset.top - inset.bottom;
+    width = width + inset.right + inset.left;
+    height = height + inset.top + inset.bottom;
+  endif
+
+  yp = fix ((index(:)-1)/columns);
+  xp = index(:) - yp*columns - 1;
+  yp = (rows - 1) - yp;
+
+  x0 = xp .* (width + margins.column) + margins.left;
+  y0 = yp .* (height + margins.row) + margins.bottom;
+
+  if (numel(x0) > 1)
+    x1 = max (x0) + width;
+    y1 = max (y0) + height;
+    x0 = min (x0);
+    y0 = min (y0);
+    pos = [x0, y0, x1-x0, y1-y0];
+  else
+    pos = [x0, y0, width, height];
+  endif
+
+endfunction
+
+%!demo
+%! clf
+%! r = 3;
+%! c = 3;
+%! fmt = {'horizontalalignment', 'center', 'verticalalignment', 'middle'};
+%! for n = 1:(r*c)
+%!   subplot (r, c, n)
+%!   xlabel (sprintf ("xlabel #%d", n))
+%!   ylabel (sprintf ("ylabel #%d", n))
+%!   title (sprintf ("title #%d", n))
+%!   text (0.5, 0.5, sprintf('subplot(%d,%d,%d)', r, c, n), fmt{:})
+%!   axis ([0 1 0 1])
+%! endfor
+%! subplot (r, c, 1:3)
+%! xlabel (sprintf ("xlabel #%d:%d", 1, 3))
+%! ylabel (sprintf ("ylabel #%d:%d", 1, 3))
+%! title (sprintf ("title #%d:%d", 1, 3))
+%! text (0.5, 0.5, sprintf('subplot(%d,%d,%d:%d)', r, c, 1, 3), fmt{:})
+%! axis ([0 1 0 1])
+