root/gliderproc/trunk/MATLAB/opnml/FCAST_1.2/contourfill3d.m

function [ho,hb]=contourfill3d(Q,Z,mesh,qlev,slev)

% [H,HB]=contour3d(Q,Z,MESH,QLEV,SLEV)

nnv=size(Z,1);
xint=slev(1);
yint=slev(2);
zint=slev(3);
h=[];
hb=[];

hold on

% CONSTANT Z

zrange= zint*ceil(min(Z(1,:))/zint+0.5):...
        abs(zint):0;

if ~isnan(zint)
for zin=zrange
        if zin<0
                l=horzslicefem(Z',Q',zin);
                c=contourctri([mesh.X mesh.Y],l,mesh.E,qlev,1);
        else
                c=contourctri([mesh.X mesh.Y],Q(nnv,:)',mesh.E,qlev,1);
        end
        while ~isempty(c)
                np=c(2,1);
                z=ones(1,np)*zin;
                x=c(1,2:np+1);
                y=c(2,2:np+1);
                h=[h;patch(x,y,z,'w',...
                        'facecolor','none',...
                        'edgecolor','c')];
                c=c(:,np+2:size(c,2));
        end
end
end

% BATHYMETRY

c=contoursurf([mesh.X mesh.Y],mesh.E',-mesh.H,zrange);
while ~isempty(c)
        np=c(2,1);
        z=ones(1,np)*c(1,1);
        x=c(1,2:np+1);
        y=c(2,2:np+1);
        hb=[hb;plot3(x,y,z,'-',...
                'color',[0.3 0.3 0.3])];
        c=c(:,np+2:size(c,2));
end

% CONSTANT X

xrange=xint*ceil(min(mesh.X)/xint):xint:xint*floor(max(mesh.X)/xint);

if ~isnan(xint)
for xin=xrange
        [l,z,d,xl,yl]=vslice(Q,Z,xin,[],60,mesh);
        hb=[hb;plot3(xl(1,:),yl(1,:),z(1,:),'-',...
                'color',[0.3 0.3 0.3])];
        
        [c,ht]=contourf(yl,z,l,[qlev qlev]);
        h=[h;ht];
        for in=ht(:)'
                z=get(in,'ydata');
                y=get(in,'xdata');
                x=ones(size(y))*xin;
                set(in,'xdata',x,...
                        'ydata',y,...
                        'zdata',z,...
                        'facecolor','none',...
                        'edgecolor','c')
        end
end
end

% CONSTANT Y

yrange=yint*ceil(min(mesh.Y)/yint):yint:yint*floor(max(mesh.Y)/yint);

if ~isnan(yint)
for yin=yrange
        [l,z,d,xl,yl]=vslice(Q,Z,[],yin,100,mesh);
        hb=[hb;plot3(xl(1,:),yl(1,:),z(1,:),'-',...
                'color',[0.3 0.3 0.3])];

        [c,ht]=contourf(xl,z,l,[qlev qlev]);
        h=[h;ht];
        for in=ht(:)'
                z=get(in,'ydata');
                x=get(in,'xdata');
                y=ones(size(x))*yin;
                set(in,'xdata',x,...
                        'ydata',y,...
                        'zdata',z,...
                        'facecolor','none',...
                        'edgecolor','c')
        end
end
end

if nargout
        ho=h;
end