function [rgb_index, rgb_map]=shade_relief_fn(z,zmin,zmax,g,gmin,gmax, huevec,grayvec)
% [rgb_index, rgb_map]=shade_relief_fn(z,zmin,zmax,g,gmin,gmax huevec,grayvec);
%                                                                 *      * 
% z is input depth array;                                         optional
% g is input gray-shade array.
% zmin, zmax are limits, matching z, to map to hues. e.g., 3500,4500
% gmin, gmax are limits, matching g, to map to gray scales, e.g., -0.8, 0.8
%
% Optional huevec and grayvec control color map.
% e.g., huevec = 0:15:300 (  0=red;   60=yellow;  120=green; 180=cyan; 
%                          240=blue; 300=magenta; 360=red)
%       grayvec = 0.2:0.1:0.9  (0=black; 0.5=full color; 1.0=white)
%
% Your 2 basic plotting choices are (caxis call is significant!):
%image(rgb_index); colormap(rgb_map)
%pcolor(rgb_index); caxis([1 size(rgb_map,1)]); shading flat; colormap(rgb_map)
%
% Peter Shaw, WHOI, 11/93

% Peter Shaw
% Woods Hole Oceanographic Institution
% prs@cook.whoi.edu
% 508 457-2000 x2473

if nargin <8, 
  Ngrays=8;    Imin=0.2; Imax=1.0;
  grayvec = Imin + (Imax-Imin)*[0:Ngrays-1]'/(Ngrays-1);
end

if nargin<7,
  Nhuesub=5; %5 max These are subdivisions of 60 degrees in hue.
  % Nhues will be simply 6+5*Nhuesub.
  dhue=60/Nhuesub;
  huevec=(0:dhue:300)';
end

Nhues=length(huevec);
Ngrays=length(grayvec);

disp([Nhues Ngrays Nhues*Ngrays])

g_index = round( 1 + (g-gmin)*(Ngrays-1) /(gmax-gmin) );
g_index(g_index<1) = ones(size( g_index(g_index<1) ));
g_index(g_index>Ngrays) = Ngrays*ones(size( g_index(g_index>Ngrays) ));

z_index = round( 1 + (z-zmin)*(Nhues-1) /(zmax-zmin) );
z_index(z_index<1) = ones(size( z_index(z_index<1) ));
z_index(z_index>Nhues) = Nhues*ones(size( z_index(z_index>Nhues) ));

rgb_index = Nhues*(g_index-1) + z_index; % :-)

rgb_map = shaded_rgbmake(huevec,grayvec);

