function[contour_levels]=colorbar4(axis_handle,cols,bar);
if (nargin==0)
    disp('_____________________________________________________________________')
    disp('Usage:       colorbar4(axis_handle,cols,bar)')
    disp('axis_handle: is a handle to current subplot axes')
    disp('cols:        is a string of gray scale [0.4 0.6 0.8 ...]')
    disp('             or cols(:,3) matrix of RGBs')
    disp('bar:         is a key, <bar=1> plots color palette')
    disp('_____________________________________________________________________')
    disp('NOTE, program works inappropriately if min(contour_levels) used')
    disp('      by <contourf(array,contour_levels> is larger than min(array).')
    disp('      If min(contour_levels)<min(array) add additional color just')
    disp('      the same as for the lowest level, i.e. ')
    disp('      colorbar4(ah,[0.4 0.4 1 0.8],bar.)')
    disp('Program returns contour levels found in the contour plot.')
    disp('If the only input argument is the axis_handle')
    disp('Program returns contour levels only.')
    disp('_____________________________________________________________________')
    disp('written by Senya: senya@atmos.umd.edu ')
    disp('ver 0.2 of 11/14/01')
    disp('ver 0.3 of 11/29/01')
    disp('ver 0.4 of 12/18/01') % adjust color palette
    disp('_____________________________________________________________________')
    return
end
ap=gca;
%switch axes layer to <top>
set(axis_handle,'layer','top');

%get handles to patch objects in current axes
hh = findobj('type','patch','parent',axis_handle);


%find contour levels of the current figure
count =0;
for i = 1:length(hh);
   if (length((get(hh(i),'cdata'))) == 1 & ~isnan(get(hh(i),'cdata')))
      count = count + 1;
      h(count) = hh(i);
      cdata(count) = get(hh(i),'cdata');
   end
end
contour_levels = unique(cdata);

if (nargin==1)
    return
end


%check if the input colors are gray scale or real colors
[m,n]=size(cols);
if(m==1)
    cols=cols'*ones(1,3);
end
if(n==1)
    cols=cols*ones(1,3);
end

%assign input colors
for i = 1:length(contour_levels)
   ind = find(cdata==contour_levels(i));
   set(h(ind),'facecolor',cols(i,:));
end

%draw color palette
if (bar==1)
        i1=isequal(cols(1,:),cols(2,:));
    %positioning of color palette
pos=get(axis_handle,'position');
ax=axes;
set(ax,'gridlinestyle','-','xgrid','on','ygrid','on',...
'position',[pos(1)+pos(3)+0.015,pos(2),max(pos(3)*0.04,0.0245),pos(4)],...
'fontsize',10,'linewidth',1,'xlim',[0 1],'ylim',[i1 length(contour_levels)],...
'xtick',[0 1],'xticklabel','','YAxisLocation', 'right',...
'ytick',[i1:length(contour_levels)])
% draw color palette
for i=i1:length(contour_levels)-1
    patch([0 1 1 0 0],[i i i+1 i+1 i],cols(i+1,:));
end
% adjust y-axis labels
set(ax,'yticklabel',[contour_levels(1+i1:end) contour_levels(end)])
lb=cellstr(get(ax,'yticklabel'));
lb(1)={''}; lb(end)={''};
set(ax,'yticklabel',lb)
end
axes(ap);