source: issm/trunk/src/m/model/plot/plot_contour.m@ 5901

function plot_contour(md,datain,options);
%PLOT_CONTOUR - plot contours of a given field
%
%   Usage:
%      plot_contour(md,data,options);
%
%   See also: PLOTMODEL

%process data and model
[x y z index is2d]=processmesh(md,[],options);
[data datatype]=processdata(md,datain,options);

%check is2d
if ~is2d,
        error('plot_contour error message: contour not supported for 3d meshes, project on a layer');
end

%first, process data: must be on grids
if datatype==1,
        %elements -> take average
        data=averaging(md,data,0);
elseif datatype==2,
        %nodes -> do nothing
elseif datatype==3,
        %quiver -> take norm
        data=sqrt(sum(datain.*datain,2));
else
        error('datatype not supported yet');
end

%prepare colors
if exist(options,'contouronly')
        %remove the previous plots
        cla
end
color=getfieldvalue(options,'contourcolor','y');

%get contours levels
contourlevels=getfieldvalue(options,'contourlevels');
if isnumeric(contourlevels),
        levels=round_ice(linspace(max(data),min(data),contourlevels),2);
else
        levels=[];
        for i=1:length(contourlevels)
                levels(end+1)=contourlevels{i};
        end
        levels=sort(unique(levels),'descend');
end
numlevels=length(levels);

%initialization of some variables
numberofelements=size(index,1);
elementslist=1:numberofelements;
c=[];
h=[];

%get unique edges in mesh
%1: list of edges
edges=[index(:,[1,2]); index(:,[2,3]); index(:,[3,1])];
%2: find unique edges
[edges,I,J]=unique(sort(edges,2),'rows');
%3: unique edge numbers
vec=J;
%4: unique edges numbers in each triangle (2 triangles sharing the same edge will have
%   the same edge number)
edges_tria=[vec(elementslist), vec(elementslist+numberofelements), vec(elementslist+2*numberofelements)];

%segments [grids1 gruids2]
Seg1=index(:,[1 2]);
Seg2=index(:,[2 3]);
Seg3=index(:,[3 1]);

%segment numbers [1;4;6;...]
Seg1_num=edges_tria(:,1);
Seg2_num=edges_tria(:,2);
Seg3_num=edges_tria(:,3);

%value of data on each tips of the segments
Data1=data(Seg1);
Data2=data(Seg2);
Data3=data(Seg3);

%get the ranges for each segment
Range1=sort(Data1,2);
Range2=sort(Data2,2);
Range3=sort(Data3,2);

for i=1:numlevels

        level=levels(i);

        %find the segments that contain this value
        pos1=(Range1(:,1)<level & Range1(:,2)>level);
        pos2=(Range2(:,1)<level & Range2(:,2)>level);
        pos3=(Range3(:,1)<level & Range3(:,2)>level);

        %get elements
        poselem12=(pos1 & pos2);
        poselem13=(pos1 & pos3);
        poselem23=(pos2 & pos3);
        poselem=find(poselem12 | poselem13 | poselem23);
        numelems=length(poselem);

        %if no element has been flagged, skip to the next level
        if numelems==0,
                continue,
        end

        %go through the elements and build the coordinates for each segment (1 by element)
        x1=zeros(numelems,1);
        x2=zeros(numelems,1);
        y1=zeros(numelems,1);
        y2=zeros(numelems,1);
        edge_l=zeros(numelems,2);

        for j=1:numelems,

                weight1=(level-Data1(poselem(j),1))/(Data1(poselem(j),2)-Data1(poselem(j),1));
                weight2=(level-Data2(poselem(j),1))/(Data2(poselem(j),2)-Data2(poselem(j),1));
                weight3=(level-Data3(poselem(j),1))/(Data3(poselem(j),2)-Data3(poselem(j),1));

                if poselem12(poselem(j));

                        x1(j)=x(Seg1(poselem(j),1))+weight1*(x(Seg1(poselem(j),2))-x(Seg1(poselem(j),1)));
                        x2(j)=x(Seg2(poselem(j),1))+weight2*(x(Seg2(poselem(j),2))-x(Seg2(poselem(j),1)));
                        y1(j)=y(Seg1(poselem(j),1))+weight1*(y(Seg1(poselem(j),2))-y(Seg1(poselem(j),1)));
                        y2(j)=y(Seg2(poselem(j),1))+weight2*(y(Seg2(poselem(j),2))-y(Seg2(poselem(j),1)));
                        edge_l(j,1)=Seg1_num(poselem(j));
                        edge_l(j,2)=Seg2_num(poselem(j));

                elseif poselem13(poselem(j)),

                        x1(j)=x(Seg1(poselem(j),1))+weight1*(x(Seg1(poselem(j),2))-x(Seg1(poselem(j),1)));
                        x2(j)=x(Seg3(poselem(j),1))+weight3*(x(Seg3(poselem(j),2))-x(Seg3(poselem(j),1)));
                        y1(j)=y(Seg1(poselem(j),1))+weight1*(y(Seg1(poselem(j),2))-y(Seg1(poselem(j),1)));
                        y2(j)=y(Seg3(poselem(j),1))+weight3*(y(Seg3(poselem(j),2))-y(Seg3(poselem(j),1)));
                        edge_l(j,1)=Seg1_num(poselem(j));
                        edge_l(j,2)=Seg3_num(poselem(j));

                elseif poselem23(poselem(j)),

                        x1(j)=x(Seg2(poselem(j),1))+weight2*(x(Seg2(poselem(j),2))-x(Seg2(poselem(j),1)));
                        x2(j)=x(Seg3(poselem(j),1))+weight3*(x(Seg3(poselem(j),2))-x(Seg3(poselem(j),1)));
                        y1(j)=y(Seg2(poselem(j),1))+weight2*(y(Seg2(poselem(j),2))-y(Seg2(poselem(j),1)));
                        y2(j)=y(Seg3(poselem(j),1))+weight3*(y(Seg3(poselem(j),2))-y(Seg3(poselem(j),1)));
                        edge_l(j,1)=Seg2_num(poselem(j));
                        edge_l(j,2)=Seg3_num(poselem(j));
                else
                        %it shoud not go here
                end
        end

        %now that we have the segments, we must try to connect them...

        %loop over the subcontours
        while ~isempty(edge_l),

                %take the right edge of the second segment and connect it to the next segments if any
                e1=edge_l(1,1);   e2=edge_l(1,2);
                xc=[x1(1);x2(1)]; yc=[y1(1);y2(1)];

                %erase the lines corresponding to this edge
                edge_l(1,:)=[];
                x1(1)=[]; x2(1)=[];
                y1(1)=[]; y2(1)=[];

                [ro1,co1]=find(edge_l==e1);

                while ~isempty(ro1)

                        if co1==1,
                                xc=[x2(ro1);xc]; yc=[y2(ro1);yc];

                                %next edge:
                                e1=edge_l(ro1,2);

                        else
                                xc=[x1(ro1);xc]; yc=[y1(ro1);yc];

                                %next edge:
                                e1=edge_l(ro1,1);
                        end

                        %erase the lines of this
                        edge_l(ro1,:)=[];
                        x1(ro1)=[]; x2(ro1)=[];
                        y1(ro1)=[]; y2(ro1)=[];

                        %next connection
                        [ro1,co1]=find(edge_l==e1);
                end

                %same thing the other way (to the right)
                [ro2,co2]=find(edge_l==e2);

                while ~isempty(ro2)

                        if co2==1,
                                xc=[xc;x2(ro2)]; yc=[yc;y2(ro2)];

                                %next edge:
                                e2=edge_l(ro2,2);
                        else
                                xc=[xc;x1(ro2)]; yc=[yc;y1(ro2)];

                                %next edge:
                                e2=edge_l(ro2,1);
                        end

                        %erase the lines of this
                        edge_l(ro2,:)=[];
                        x1(ro2)=[]; x2(ro2)=[];
                        y1(ro2)=[]; y2(ro2)=[];

                        %next connection
                        [ro2,co2]=find(edge_l==e2);
                end

                %we now have one subcontour ready to be plotted
                zc=level*ones(length(xc)+1,1);
                if getfieldvalue(options,'contouronly',0),
                        h=[h;patch('Xdata',[xc;NaN],'Ydata',[yc;NaN],'Zdata',zc,'Cdata',zc,'facecolor','none','edgecolor','flat')];
                        hold on      
                else
                        h=[h;patch('Xdata',[xc;NaN],'Ydata',[yc;NaN],'facecolor','none','edgecolor',color)];
                        hold on
                end

                % Update the CS data structure as per "contours.m"
                % so that clabel works
                c = horzcat(c,[level, xc'; length(xc), yc']);

        end
end

%labels?
if (~strcmpi(getfieldvalue(options,'contourticks','on'),'off') & ~isempty(c) & ~isempty(h))
        if exist(options,'contouronly')
                clabel(c,h);
        else
                clabel(c,h,'color',color,'FontSize',14);
        end
end