gusucode.com > DSISoft是由加拿大地质调查局发布的用于垂直地震剖面(VSP)数据处理的免费软件包 > dsisoftv3/cdpt3d/plotbound.m
function plotbound(bound,bins) %plotbound(bound,bins) % %this function plots boundary layers on current figure (3D) %bound is a cell array describing one or more layers %each layer can be either a 1x5 matrix of x,y,z coordinates of a point %on the layer followed by strike and dip in degrees %or else as an n x 3 matrix of x,y,z points on the surface %bins describes the number of grid points to be plotted % %written by K.S. Beaty %last modified November 5, 1997 w=pi./180; %convert strike and dip to radians colormap(cool); limx = get(gca, 'XLim'); limy = get(gca, 'YLim'); xmax = limx(1,2); ymax = limy(1,2); xmin = limx(1,1); ymin = limy(1,1); dy = ymax-ymin; dx = xmax-xmin; %creating a meshgrid for boundaries x = [xmin:(xmax-xmin)/bins:xmax]; y = [ymin:(ymax-ymin)/bins:ymax]; [X,Y] = meshgrid(x,y); C = ones(bins+1,bins+1); %start loop to plot all boundaries boundnum=size(bound); layers=boundnum(1,2); for i=1:layers col=size(bound{i}); if col(1,2)==5 strike = bound{i}(1,4).*pi./180; dip = bound{i}(1,5).*pi./180; zb = bound{i}(1,3); yb = bound{i}(1,2); xb = bound{i}(1,1); z0 = zb - ((xmax-xb)*cos(strike)-(yb-ymin)*sin(strike))*tan(dip); z1 = z0 - dy*sin(strike)*tan(dip); z2 = z0 + (dx-dy*tan(strike))*cos(strike)*tan(dip); z3 = z0 + dx*cos(strike)*tan(dip); corners(1, :) = [xmin, ymin, z3]; corners(2, :) = [xmin, ymax, z2]; corners(3, :) = [xmax, ymax, z1]; corners(4, :) = [xmax, ymin, z0]; %corners(5, :) = corners(1,:); points=corners; else points=bound{i}(:,1:3); end; %if/else %plot3(points(:,1), points(:,2), points(:,3), 'kx'); Z=griddata(points(:,1), points(:,2), points(:,3), X, Y, 'linear'); %mesh(x,y,Z,C); H=surf(x,y,Z,C); set(H,'facecolor', 'none'); set(H,'edgecolor','flat'); C = C+1; %increment colour matrix for next plane planegrid{i}=Z; end %plotting boundaries loop hidden off;