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

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

% [H,HB]=slice3d(Q,Z,MESH,SLEV)
%
% Generates a wireframe of the isosurface of Q at value QLEV, with the
% interval between frames being described by SLEV
%
%     Q = Scalar field (nnv x nn)
%     Z = Z locations of scalar field (nnv x nn)
%     mesh = Structure loaded by MESHSTR (with E,X,Y fields)
%     slev = intervals between x, y and z wireframes (3 x 1)

% OPNML COMPLIANT 3-18-99

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

if ~isnan(xint)
        xrange=xint*ceil(min(mesh.x)/xint):...
                xint:...
                xint*floor(max(mesh.x)/xint);
else
        xrange=[];
end
if ~isnan(zint)
        zrange= zint*ceil(min(Z(1,:))/zint+0.5):...
                abs(zint):-5;
else
        zrange=[];
end
if ~isnan(yint)
        yrange=yint*ceil(min(mesh.y)/yint):...
                yint:...
                yint*floor(max(mesh.y)/yint);
else
        yrange=[];
end
hold on

% CONSTANT Z

for zin=zrange
        x=mesh.x;
        y=mesh.y;
        e=mesh.e;
        if zin<0
                l=horzslicefem(Z,Q,zin);
                ht=colormesh2d(mesh,l);
        else
                ht=colormesh2d(mesh,Q(nnv,:)');
        end
        set(ht,'zdata',get(ht,'xdata')*0+zin)
        h=[h;ht];
end

% CONSTANT X

for xin=xrange
        [l,z,d,xl,yl]=vslice(Q,Z,xin,[],100,mesh);
        ht=surf(xl,yl,z,l);
        shading('interp')
        h=[h;ht];
        plot3(xl(:,1),yl(:,1),z(:,1),'k')
        plot3(xl(:,size(xl,2)),yl(:,size(xl,2)),z(:,size(xl,2)),'k')
        plot3(xl(size(xl,1),:),yl(size(xl,1),:),z(size(xl,1),:),'k')
        plot3(xl(1,:),yl(1,:),z(1,:),'k')
end

% CONSTANT Y

for yin=yrange
        [l,z,d,xl,yl]=vslice(Q,Z,[],yin,150,mesh);
        ht=surf(xl,yl,z,l);
        shading('interp')
        h=[h;ht];
        plot3(xl(:,1),yl(:,1),z(:,1),'k')
        plot3(xl(:,size(xl,2)),yl(:,size(xl,2)),z(:,size(xl,2)),'k')
        plot3(xl(size(xl,1),:),yl(size(xl,1),:),z(size(xl,1),:),'k')
        plot3(xl(1,:),yl(1,:),z(1,:),'k')
end

if nargout
        ho=h;
end

for xin=xrange
for yin=yrange
        plot3(  xin*ones(2,1),...
                yin*ones(2,1),...
                [min(Z(:)) max(Z(:))],...
                'k-');
end
end

axis([-inf inf -inf inf min(Z(:)) max(Z(:))]);