function haxes = bar3D(xfp,yfp,DI,xlab,ylab,zlab,cutoff,haxes)
%
%  Name: bar3D
% 
%  Usage: haxes = bar3D(xfp,yfp,DI,xlab,ylab,zlab,cutoff,haxes)
%
% 
%   Description:
%     Draws a 3D bar chart of the variable DI
%
%   inputs:
%	yfp(j,i) = y coordinate of the (j,i) element of DI
%	xfp(j,i) = x coordinate of the (j,i) element of DI
%	DI(j,i)  = Value of the function at x=xfp(j,i) and y=yfp(j,i)
%	haxes    = handle of the axes to use to make the chart
%	xlab, ylab, zlab = label strings for axes x, y, z
%	cutoff   = 'critical' value of DI -- bars where DI > cutoff 
%		    will be colored red (only used when cutoff > 0)

%  Version SWD970108
%
%  Author(s): P. Cornwell, cornwell@lanl.gov
%	modified by Scott W. Doebling, doebling@lanl.gov
%
%  Institution: Los Alamos National Laboratory,
%		ESA-EA
%
%  Copyright (c) 1996-1997, The Regents of the University
%     	of California. (See top level of code for
%	full copyright statement.)

if ~exist('haxes')
	haxes = axes;
elseif ~ishandle(haxes),
	haxes = axes;
end

if ~exist('xfpcount'),
	xfpcount = size(DI,1)*ones(1,size(DI,2));
end

axes(haxes);

cla

dx=abs((xfp(1,1)-xfp(1,2)))*.4;
dy=abs((yfp(1,1)-yfp(2,1)))*.4;
numdivy=size(DI,1);
numdivx=size(DI,2);

for i=1:numdivx
   for j=1:numdivy 
    	
	x1=xfp(j,i)+dx;
	x2=xfp(j,i)-dx;
	y1=yfp(j,i)-dy;
	y2=yfp(j,i)+dy;

	z=DI(j,i);
	tx11=[x1,x1,x1,x1];
	tx22=[x2,x2,x2,x2];
	ty11=[y1,y1,y1,y1];
	ty22=[y2,y2,y2,y2];
	tx12=[x1 x2 x2 x1];
	ty12=[y1 y2 y2 y1];
	ty1122=[y1 y1 y2 y2];
	tz=[z z 0 0];
	tzz=[z z z z];
	t00=[0 0 0 0];

	if cutoff > 0,
	  if abs(z)<cutoff
		color='r';
	  else,
		color='r';
	  end;
	else
      color='r';	
   end

	fill3(tx11,ty12,tz,color,tx12,ty22,tz,color,tx22,ty12,tz,color,tx12,ty11,tz,color,tx12,ty1122,tzz,color,tx12,ty1122,t00,color,'edgecolor','b');
	set(haxes,'nextplot','add')

    end
end

set(haxes,'nextplot','replace','box','on')
 
if exist('xlab'),
   xlabel(xlab),ylabel(ylab),zlabel(zlab)
end

return