plot_contour

PURPOSE ^

PLOT_CONTOUR - plot contours of a given field

SYNOPSIS ^

function plot_contour(md,data,options_structure);

DESCRIPTION ^

PLOT_CONTOUR - plot contours of a given field

   Usage:
      plot_contour(md,data,options_structure);

   See also: PLOTMODEL

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function plot_contour(md,data,options_structure);
0002 %PLOT_CONTOUR - plot contours of a given field
0003 %
0004 %   Usage:
0005 %      plot_contour(md,data,options_structure);
0006 %
0007 %   See also: PLOTMODEL
0008 
0009 %process data and model
0010 [x y z index is2d]=processmesh(md,options_structure);
0011 [data isongrid]=processdata(md,data,options_structure);
0012 
0013 %check is2d
0014 if ~is2d,
0015     error('plot_contour error message: contour not supported for 3d meshes, project on a layer');
0016 end
0017 
0018 %first, process data: must be on grids
0019 if ~isongrid
0020     data=averaging(md,data,0);
0021 end
0022 
0023 %prepare colors
0024 if isnan(options_structure.contouronly)
0025     %contour color
0026     if isnan(options_structure.contourcolor)
0027         options_structure.contourcolor='y';
0028     end
0029 else
0030     %remove the previous plots
0031     cla
0032 end
0033 color=options_structure.contourcolor;
0034 
0035 %get contours levels
0036 if isnumeric(options_structure.contourlevels),
0037     levels=round_ice(linspace(max(data),min(data),options_structure.contourlevels),2);
0038 else
0039     levels=[];
0040     for i=1:length(options_structure.contourlevels)
0041         levels(end+1)=options_structure.contourlevels{i};
0042     end
0043     levels=sort(unique(levels),'descend');
0044 end
0045 numlevels=length(levels);
0046 
0047 %initialization of some variables
0048 numberofelements=md.numberofelements;
0049 elementslist=1:numberofelements;
0050 c=[];
0051 h=[];
0052 
0053 %get unique edges in mesh
0054 %1: list of edges
0055 edges=[index(:,[1,2]); index(:,[2,3]); index(:,[3,1])];
0056 %2: find unique edges
0057 [edges,I,J]=unique(sort(edges,2),'rows');
0058 %3: unique edge numbers
0059 vec=J;
0060 %4: unique edges numbers in each triangle (2 triangles sharing the same edge will have
0061 %   the same edge number)
0062 edges_tria=[vec(elementslist), vec(elementslist+numberofelements), vec(elementslist+2*numberofelements)];
0063 
0064 %segments [grids1 gruids2]
0065 Seg1=index(:,[1 2]);
0066 Seg2=index(:,[2 3]);
0067 Seg3=index(:,[3 1]);
0068 
0069 %segment numbers [1;4;6;...]
0070 Seg1_num=edges_tria(:,1);
0071 Seg2_num=edges_tria(:,2);
0072 Seg3_num=edges_tria(:,3);
0073 
0074 %value of data on each tips of the segments
0075 Data1=data(Seg1);
0076 Data2=data(Seg2);
0077 Data3=data(Seg3);
0078 
0079 %get the ranges for each segment
0080 Range1=sort(Data1,2);
0081 Range2=sort(Data2,2);
0082 Range3=sort(Data3,2);
0083 
0084 for i=1:numlevels
0085 
0086     level=levels(i);
0087 
0088     %find the segments that contain this value
0089     pos1=(Range1(:,1)<level & Range1(:,2)>level);
0090     pos2=(Range2(:,1)<level & Range2(:,2)>level);
0091     pos3=(Range3(:,1)<level & Range3(:,2)>level);
0092 
0093     %get elements
0094     poselem12=(pos1 & pos2);
0095     poselem13=(pos1 & pos3);
0096     poselem23=(pos2 & pos3);
0097     poselem=find(poselem12 | poselem13 | poselem23);
0098     numelems=length(poselem);
0099 
0100     %if no element has been flagged, skip to the next level
0101     if numelems==0,
0102         continue,
0103     end
0104 
0105     %go through the elements and build the coordinates for each segment (1 by element)
0106     x1=zeros(numelems,1);
0107     x2=zeros(numelems,1);
0108     y1=zeros(numelems,1);
0109     y2=zeros(numelems,1);
0110     edge_l=zeros(numelems,2);
0111 
0112     for j=1:numelems,
0113 
0114         weight1=(level-Data1(poselem(j),1))/(Data1(poselem(j),2)-Data1(poselem(j),1));
0115         weight2=(level-Data2(poselem(j),1))/(Data2(poselem(j),2)-Data2(poselem(j),1));
0116         weight3=(level-Data3(poselem(j),1))/(Data3(poselem(j),2)-Data3(poselem(j),1));
0117 
0118         if poselem12(poselem(j));
0119 
0120             x1(j)=x(Seg1(poselem(j),1))+weight1*(x(Seg1(poselem(j),2))-x(Seg1(poselem(j),1)));
0121             x2(j)=x(Seg2(poselem(j),1))+weight2*(x(Seg2(poselem(j),2))-x(Seg2(poselem(j),1)));
0122             y1(j)=y(Seg1(poselem(j),1))+weight1*(y(Seg1(poselem(j),2))-y(Seg1(poselem(j),1)));
0123             y2(j)=y(Seg2(poselem(j),1))+weight2*(y(Seg2(poselem(j),2))-y(Seg2(poselem(j),1)));
0124             edge_l(j,1)=Seg1_num(poselem(j));
0125             edge_l(j,2)=Seg2_num(poselem(j));
0126 
0127         elseif poselem13(poselem(j)),
0128 
0129             x1(j)=x(Seg1(poselem(j),1))+weight1*(x(Seg1(poselem(j),2))-x(Seg1(poselem(j),1)));
0130             x2(j)=x(Seg3(poselem(j),1))+weight3*(x(Seg3(poselem(j),2))-x(Seg3(poselem(j),1)));
0131             y1(j)=y(Seg1(poselem(j),1))+weight1*(y(Seg1(poselem(j),2))-y(Seg1(poselem(j),1)));
0132             y2(j)=y(Seg3(poselem(j),1))+weight3*(y(Seg3(poselem(j),2))-y(Seg3(poselem(j),1)));
0133             edge_l(j,1)=Seg1_num(poselem(j));
0134             edge_l(j,2)=Seg3_num(poselem(j));
0135 
0136         elseif poselem23(poselem(j)),
0137 
0138             x1(j)=x(Seg2(poselem(j),1))+weight2*(x(Seg2(poselem(j),2))-x(Seg2(poselem(j),1)));
0139             x2(j)=x(Seg3(poselem(j),1))+weight3*(x(Seg3(poselem(j),2))-x(Seg3(poselem(j),1)));
0140             y1(j)=y(Seg2(poselem(j),1))+weight2*(y(Seg2(poselem(j),2))-y(Seg2(poselem(j),1)));
0141             y2(j)=y(Seg3(poselem(j),1))+weight3*(y(Seg3(poselem(j),2))-y(Seg3(poselem(j),1)));
0142             edge_l(j,1)=Seg2_num(poselem(j));
0143             edge_l(j,2)=Seg3_num(poselem(j));
0144         else
0145             %it shoud not go here
0146         end
0147     end
0148 
0149     %now that we have the segments, we must try to connect them...
0150 
0151     %loop over the subcontours
0152     while ~isempty(edge_l),
0153 
0154         %take the right edge of the second segment and connect it to the next segments if any
0155         e1=edge_l(1,1);   e2=edge_l(1,2);
0156         xc=[x1(1);x2(1)]; yc=[y1(1);y2(1)];
0157 
0158         %erase the lines corresponding to this edge
0159         edge_l(1,:)=[];
0160         x1(1)=[]; x2(1)=[];
0161         y1(1)=[]; y2(1)=[];
0162 
0163         [ro1,co1]=find(edge_l==e1);
0164 
0165         while ~isempty(ro1)
0166 
0167             if co1==1,
0168                 xc=[x2(ro1);xc]; yc=[y2(ro1);yc];
0169 
0170                 %next edge:
0171                 e1=edge_l(ro1,2);
0172 
0173             else
0174                 xc=[x1(ro1);xc]; yc=[y1(ro1);yc];
0175 
0176                 %next edge:
0177                 e1=edge_l(ro1,1);
0178             end
0179 
0180             %erase the lines of this
0181             edge_l(ro1,:)=[];
0182             x1(ro1)=[]; x2(ro1)=[];
0183             y1(ro1)=[]; y2(ro1)=[];
0184 
0185             %next connection
0186             [ro1,co1]=find(edge_l==e1);
0187         end
0188 
0189         %same thing the other way (to the right)
0190         [ro2,co2]=find(edge_l==e2);
0191 
0192         while ~isempty(ro2)
0193 
0194             if co2==1,
0195                 xc=[xc;x2(ro2)]; yc=[yc;y2(ro2)];
0196 
0197                 %next edge:
0198                 e2=edge_l(ro2,2);
0199             else
0200                 xc=[xc;x1(ro2)]; yc=[yc;y1(ro2)];
0201 
0202                 %next edge:
0203                 e2=edge_l(ro2,1);
0204             end
0205 
0206             %erase the lines of this
0207             edge_l(ro2,:)=[];
0208             x1(ro2)=[]; x2(ro2)=[];
0209             y1(ro2)=[]; y2(ro2)=[];
0210 
0211             %next connection
0212             [ro2,co2]=find(edge_l==e2);
0213         end
0214 
0215         %we now have one subcontour ready to be plotted
0216         zc=level*ones(length(xc)+1,1);
0217         if isnan(color),
0218             h=[h;patch('Xdata',[xc;NaN],'Ydata',[yc;NaN],'Zdata',zc,'Cdata',zc,'facecolor','none','edgecolor','flat')];
0219             hold on      
0220         else
0221             h=[h;patch('Xdata',[xc;NaN],'Ydata',[yc;NaN],'facecolor','none','edgecolor',color)];
0222             hold on
0223         end
0224 
0225         % Update the CS data structure as per "contours.m"
0226         % so that clabel works
0227         c = horzcat(c,[level, xc'; length(xc), yc']);
0228 
0229     end
0230 end
0231 
0232 %labels?
0233 if ~strcmpi(options_structure.contourticks,'off')
0234     if ~isnan(options_structure.contourcolor)
0235         clabel(c,h,'color',color);
0236     else
0237         clabel(c,h);
0238     end
0239 end

Generated on Sun 29-Mar-2009 20:22:55 by m2html © 2003