source: trunk/src/plot_field.m @ 688

Last change on this file since 688 was 688, checked in by sommeria, 11 years ago

bug correction on xml time reading
some cleaning of the GUI uvmat
improvement of option 'contours' in plot_field

File size: 58.1 KB
RevLine 
[201]1%'plot_field': plot any field with the structure defined in the uvmat package
2%------------------------------------------------------------------------
3%
4%  This function is used by uvmat to plot fields. It automatically chooses the representation
5% appropriate to the input field structure:
[292]6%     2D vector fields are represented by arrows, 2D scalar fields by grey scale images or contour plots, 1D fields are represented by usual plot with (abscissa, ordinate).
[201]7%  The input field structure is first tested by check_field_structure.m,
[512]8%  then split into blocks of related variables  by find_field_cells.m.
[292]9%  The dimensionality of each block is obtained  by this function
[201]10%  considering the presence of variables with the attribute .Role='coord_x'
11%  and/or coord_y and/or coord_z (case of unstructured coordinates), or
12%  dimension variables (case of matrices).
13%
[625]14% function [PlotType,PlotParamOut,haxes]= plot_field(Data,haxes,PlotParam,PosColorbar)
[201]15%
16% OUPUT:
17% PlotType: type of plot: 'text','line'(curve plot),'plane':2D view,'volume'
18% PlotParamOut: structure, representing the updated  plotting parameters, in case of automatic scaling
19% haxes: handle of the plotting axis, when a new figure is created.
20%
21%INPUT
22%    Data:   structure describing the field to plot
23%         (optional) .ListGlobalAttribute: cell listing the names of the global attributes
24%                    .Att_1,Att_2... : values of the global attributes
25%         (requested)  .ListVarName: list of variable names to select (cell array of  char strings {'VarName1', 'VarName2',...} )
26%         (requested)  .VarDimName: list of dimension names for each element of .ListVarName (cell array of string cells)
27%                      .VarAttribute: cell of attributes for each element of .ListVarName (cell array of structures of the form VarAtt.key=value)
28%         (requested) .Var1, .Var2....: variables (Matlab arrays) with names listed in .ListVarName
29
30%            Variable attribute .Role :
31%    The only variable attribute used for plotting purpose is .Role which can take
32%    the values
33%       Role = 'scalar':  (default) represents a scalar field
34%            = 'coord_x', 'coord_y',  'coord_z': represents a separate set of
35%                        unstructured coordinate x, y  or z
36%            = 'vector': represents a vector field whose number of components
37%                is given by the last dimension (called 'nb_dim')
38%            = 'vector_x', 'vector_y', 'vector_z'  :represents the x, y or z  component of a vector 
39%            = 'warnflag' : provides a warning flag about the quality of data in a 'Field', default=0, no warning
40%            = 'errorflag': provides an error flag marking false data,
41%                   default=0, no error. Different non zero values can represent different criteria of elimination.
42%
43%   haxes: handle of the plotting axes to update with the new plot. If this input is absent or not a valid axes handle, a new figure is created.
44%
[294]45%   PlotParam: structure containing the parameters for plotting, as read on the uvmat or view_field GUI (by function 'read_GUI.m').
46%      Contains three substructures:
47%     .Coordinates: coordinate parameters:
48%           .CheckFixLimits:=0 (default) adjust axes limit to the X,Y data, =1: preserves the previous axes limits
[428]49%     .Coordinates.CheckFixAspectRatio: =0 (default):automatic adjustment of the graph, keep 1 to 1 aspect ratio for x and y scales.
[429]50%     .Coordinates.AspectRatio: imposed aspect ratio y/x of axis unit plots
[201]51%            --scalars--
52%    .Scalar.MaxA: upper bound (saturation color) for the scalar representation, max(field) by default
53%    .Scalar.MinA: lower bound (saturation) for the scalar representation, min(field) by default
[292]54%    .Scalar.CheckFixScal: =0 (default) lower and upper bounds of the scalar representation set to the min and max of the field
[201]55%               =1 lower and upper bound imposed by .AMax and .MinA
[421]56%    .Scalar.CheckBW= 1: black and white representation imposed, =0 color imposed (color scale or rgb),
57%                   =[]: automatic (B/W for integer positive scalars, color  else)
[292]58%    .Scalar.CheckContours= 1: represent scalars by contour plots (Matlab function 'contour'); =0 by default
[201]59%    .IncrA : contour interval
60%            -- vectors--
61%    .Vectors.VecScale: scale for the vector representation
[292]62%    .Vectors.CheckFixVec: =0 (default) automatic length for vector representation, =1: length set by .VecScale
63%    .Vectors.CheckHideFalse= 0 (default) false vectors represented in magenta, =1: false vectors not represented;
64%    .Vectors.CheckHideWarning= 0 (default) vectors marked by warnflag~=0 marked in black, 1: no warning representation;
65%    .Vectors.CheckDecimate4 = 0 (default) all vectors reprtesented, =1: half of  the vectors represented along each coordinate
[201]66%         -- vector color--
67%    .Vectors.ColorCode= 'black','white': imposed color  (default ='blue')
68%                        'rgb', : three colors red, blue, green depending
69%                        on thresholds .colcode1 and .colcode2 on the input  scalar value (C)
70%                        'brg': like rgb but reversed color order (blue, green, red)
71%                        '64 colors': continuous color from blue to red (multijet)
72%    .Vectors.colcode1 : first threshold for rgb, first value for'continuous'
73%    .Vectors.colcode2 : second threshold for rgb, last value (saturation) for 'continuous'
[292]74%    .Vectors.CheckFixedCbounds;  =0 (default): the bounds on C representation are min and max, =1: they are fixed by .Minc and .MaxC
[201]75%    .Vectors.MinC = imposed minimum of the scalar field used for vector color;
76%    .Vectors.MaxC = imposed maximum of the scalar field used for vector color;
77%
[546]78% PosColorbar: % if absent, no action on colorbar
79%              % if empty, suppress any existing colorbar
80%              % if not empty, display a colorbar for B&W images at position PosColorbar
81%                expressed in figure relative unit (ex [0.821 0.471 0.019 0.445])
[201]82
83%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
84%  Copyright Joel Sommeria, 2008, LEGI / CNRS-UJF-INPG, sommeria@coriolis-legi.org.
85%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
86%     This file is part of the toolbox UVMAT.
87%
88%     UVMAT is free software; you can redistribute it and/or modify
89%     it under the terms of the GNU General Public License as published by
90%     the Free Software Foundation; either version 2 of the License, or
91%     (at your option) any later version.
92%
93%     UVMAT is distributed in the hope that it will be useful,
94%     but WITHOUT ANY WARRANTY; without even the implied warranty of
95%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
96%     GNU General Public License (file UVMAT/COPYING.txt) for more details.
97%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
98
99function [PlotType,PlotParamOut,haxes]= plot_field(Data,haxes,PlotParam,PosColorbar)
[247]100
[388]101%% default input and output
[201]102if ~exist('PlotParam','var'),PlotParam=[];end;
103PlotType='text'; %default
[411]104if ~isfield(PlotParam,'Coordinates')
105    PlotParam.Coordinates=[];
[582]106    if isfield(Data,'CoordUnit')
107        PlotParam.Coordinates.CheckFixAspectRatio=1;
108        PlotParam.Coordinates.AspectRatio=1; %set axes equal by default if CoordUnit is defined
109    end
[292]110end
[569]111PlotParamOut=PlotParam;%default
[201]112
[388]113%% check input structure
114index_2D=[];
115index_1D=[];
116index_0D=[];
117% check the cells of fields :
[530]118[CellInfo,NbDimArray,errormsg]=find_field_cells(Data);
119%[CellVarIndex,NbDim,CoordType,VarRole,errormsg]=find_field_cells(Data);
[388]120if ~isempty(errormsg)
[512]121    msgbox_uvmat('ERROR',['input of plot_field/find_field_cells: ' errormsg]);
[388]122    return
123end
[530]124
125index_3D=find(NbDimArray>2,1);
[388]126if ~isempty(index_3D)
[546]127    msgbox_uvmat('ERROR','volume plot not implemented yet');
128    return
[388]129end
[530]130index_2D=find(NbDimArray==2);%find 2D fields
131index_1D=find(NbDimArray==1);
132index_0D=find(NbDimArray==0);
[294]133
[388]134%% test axes and figure
[429]135testnewfig=1;%test to create a new figure (default)
136testzoomaxes=0;%test for the existence of a zoom secondary figure attached to the plotting axes
137if exist('haxes','var')
138    if ishandle(haxes)
139        if isequal(get(haxes,'Type'),'axes')
140            testnewfig=0;
141            AxeData=get(haxes,'UserData');
142            if isfield(AxeData,'ZoomAxes')&& ishandle(AxeData.ZoomAxes)
143                if isequal(get(AxeData.ZoomAxes,'Type'),'axes')
144                    testzoomaxes=1;
145                    zoomaxes=AxeData.ZoomAxes;
[388]146                end
147            end
148        end
[201]149    end
[429]150end
151
152%% create a new figure and axes if the plotting axes does not exist
153if testnewfig
154    hfig=figure;
155    set(hfig,'Units','normalized')
156    haxes=axes;
157    set(haxes,'position',[0.13,0.2,0.775,0.73])
[569]158    PlotParamOut.NextPlot='add'; %parameter for plot_profile
[429]159else
160    hfig=get(haxes,'parent');
161    set(0,'CurrentFigure',hfig)% the parent of haxes becomes the current figure
162    set(hfig,'CurrentAxes',haxes)%  haxes becomes the current axes of the parent figure
163end
164
165%% set axes properties
[569]166if isfield(PlotParamOut.Coordinates,'CheckFixLimits') && isequal(PlotParamOut.Coordinates.CheckFixLimits,1)  %adjust the graph limits
[429]167    set(haxes,'XLimMode', 'manual')
168    set(haxes,'YLimMode', 'manual')
169else
170    set(haxes,'XLimMode', 'auto')
171    set(haxes,'YLimMode', 'auto')
172end
173% if ~isfield(PlotParam.Coordinates,'CheckFixAspectRatio')&& isfield(Data,'CoordUnit')
174%     PlotParam.Coordinates.CheckFixAspectRatio=1;% if CoordUnit is defined, the two coordiantes should be plotted with equal scale by default
175% end
176errormsg='';
[569]177%PlotParamOut.Coordinates=[]; %default output
[429]178AxeData=get(haxes,'UserData');
179
180%% 2D plots
181if isempty(index_2D)
[546]182    plot_plane([],[],haxes,[],[]);%removes images or vector plots in the absence of 2D field plot
[429]183else  %plot 2D field
[546]184    if ~exist('PosColorbar','var'),PosColorbar=[];end;
[569]185    [tild,PlotParamOut,PlotType,errormsg]=plot_plane(Data,CellInfo(index_2D),haxes,PlotParamOut,PosColorbar);
[429]186    AxeData.NbDim=2;
187    if testzoomaxes && isempty(errormsg)
[569]188        [zoomaxes,PlotParamOut,tild,errormsg]=plot_plane(Data,CellInfo(index_2D),zoomaxes,PlotParamOut,PosColorbar);
[429]189        AxeData.ZoomAxes=zoomaxes;
[201]190    end
[429]191end
192
193%% 1D plot (usual graph y vs x)
194if isempty(index_1D)
195    if ~isempty(haxes)
[530]196        plot_profile([],[],haxes);%removes usual praphs y vs x in the absence of 1D field plot
[201]197    end
[429]198else %plot 1D field (usual graph y vs x)
[644]199    CheckHold=0;
200    if isfield(PlotParam,'CheckHold')
201        CheckHold= PlotParam.CheckHold;
202    end       
203    PlotParamOut.Coordinates=plot_profile(Data,CellInfo(index_1D),haxes,PlotParamOut.Coordinates,CheckHold);%
[429]204    if testzoomaxes
[644]205        [zoomaxes,PlotParamOut.Coordinates]=plot_profile(Data,CellInfo(index_1D),zoomaxes,PlotParamOut.Coordinates,CheckHold);
[429]206        AxeData.ZoomAxes=zoomaxes;
[201]207    end
[429]208    PlotType='line';
209end
210
211%% text display
212if isempty(index_2D) && isempty(index_1D)%text display alone
213    htext=findobj(hfig,'Tag','TableDisplay');
214else  %text display added to plot
[388]215    htext=findobj(hfig,'Tag','text_display');
[429]216end
217if ~isempty(htext)
218    if isempty(index_0D)
[434]219        if strcmp(get(htext,'Type'),'uitable')
220            set(htext,'Data',{})
221        else
222            set(htext,'String',{''})
223        end
[429]224    else
[530]225        [errormsg]=plot_text(Data,CellInfo(index_0D),htext);
[294]226    end
[201]227end
228
[388]229%% display error message
[201]230if ~isempty(errormsg)
231    msgbox_uvmat('ERROR', errormsg)
232end
233
234%% update the parameters stored in AxeData
[429]235if ishandle(haxes)&&( ~isempty(index_2D)|| ~isempty(index_1D))
236%     AxeData=[];
[388]237    if isfield(PlotParamOut,'MinX')
[569]238        AxeData.RangeX=[PlotParamOut.MinX PlotParamOut.MaxX];
239        AxeData.RangeY=[PlotParamOut.MinY PlotParamOut.MaxY];
[388]240    end
241    set(haxes,'UserData',AxeData)
242end
[201]243
[206]244%% update the plotted field stored in parent figure
[429]245if ~isempty(index_2D)|| ~isempty(index_1D)
246    FigData=get(hfig,'UserData');
[511]247    if strcmp(get(hfig,'tag'),'view_field')||strcmp(get(hfig,'tag'),'uvmat')
248        FigData.(get(haxes,'tag'))=Data;
[429]249        set(hfig,'UserData',FigData)
250    end
[388]251end
[292]252
[201]253%-------------------------------------------------------------------
[530]254function errormsg=plot_text(FieldData,CellInfo,htext)
[201]255%-------------------------------------------------------------------
[581]256errormsg='';
[201]257txt_cell={};
[429]258Data={};
[530]259for icell=1:length(CellInfo)
[535]260   
261    % select types of  variables to be projected
262    ListProj={'VarIndex_scalar','VarIndex_image','VarIndex_color','VarIndex_vector_x','VarIndex_vector_y'};
263    check_proj=false(size(FieldData.ListVarName));
264    for ilist=1:numel(ListProj)
265        if isfield(CellInfo{icell},ListProj{ilist})
266            check_proj(CellInfo{icell}.(ListProj{ilist}))=1;
267        end
268    end
269    VarIndex=find(check_proj);
270    %
271    %     VarIndex=CellInfo{icell}.VarIndex;%  indices of the selected variables in the list data.ListVarName
272    %     for ivar=1:length(VarIndex)
273    %         checkancillary=0;
274    %         if length(FieldData.VarAttribute)>=VarIndex(ivar)
275    %             VarAttribute=FieldData.VarAttribute{VarIndex(ivar)};
276    %             if isfield(VarAttribute,'Role')&&(strcmp(VarAttribute.Role,'ancillary')||strcmp(VarAttribute.Role,'coord_tps')...
277    %                     ||strcmp(VarAttribute.Role,'vector_x_tps')||strcmp(VarAttribute.Role,'vector_y_tps'))
278    %                 checkancillary=1;
279    %             end
280    %         end
281    %         if ~checkancillary% does not display variables with attribute '.Role=ancillary'
[201]282    for ivar=1:length(VarIndex)
[535]283    VarName=FieldData.ListVarName{VarIndex(ivar)};
284    VarValue=FieldData.(VarName);
[609]285    if isvector(VarValue')
[535]286        VarValue=VarValue';% put the different values on a line
287    end
288    if numel(VarValue)>1 && numel(VarValue)<10
289        for ind=1:numel(VarValue)
290            VarNameCell{1,ind}=[VarName '_' num2str(ilist)];
[389]291        end
[535]292    else
293        VarNameCell={VarName};
[201]294    end
[535]295    if numel(VarValue)<10
[537]296        if isempty(VarValue)
297            VarValueCell={'[]'};
298        else
299            VarValueCell=num2cell(VarValue);
300        end
301        if isempty(Data)
302            Data =[VarNameCell; VarValueCell];
303        else
304            Data =[Data [VarNameCell; VarValueCell]];
305        end
[535]306    else
[537]307        if isempty(Data)
308            Data =[VarNameCell; num2cell(VarValue)];
309        else
[535]310        Data =[Data [VarNameCell; {['size ' num2str(size(VarValue))]}]];
[537]311        end
[535]312    end
313    if size(VarValue,1)==1
314        txt=[VarName '=' num2str(VarValue)];
315        txt_cell=[txt_cell;{txt}];
316    end
317    end
[201]318end
[429]319if strcmp(get(htext,'Type'),'uitable')
320    get(htext,'ColumnName')
321    set(htext,'ColumnName',Data(1,:))
322    set(htext,'Data',Data(2:end,:))
323else
324    set(htext,'String',txt_cell)
325    set(htext,'UserData',txt_cell)% for storage during mouse display
326end
[201]327
[429]328
[201]329%-------------------------------------------------------------------
[644]330function CoordinatesOut=plot_profile(data,CellInfo,haxes,Coordinates,CheckHold)
[201]331%-------------------------------------------------------------------
332
[644]333%% initialization
[292]334if ~exist('Coordinates','var')
335    Coordinates=[];
[201]336end
[292]337CoordinatesOut=Coordinates; %default
[201]338hfig=get(haxes,'parent');
[644]339legend_str={};
340
341%% suppress existing plot if empty data
[201]342if isempty(data)
343    hplot=findobj(haxes,'tag','plot_line');
344    if ~isempty(hplot)
345        delete(hplot)
346    end
347    hlegend=findobj(hfig,'tag','legend');
348    if ~isempty(hlegend)
349        delete(hlegend)
350    end
351    return
352end
353
[644]354%% set the colors of the successive plots (designed to produce rgb for the three components of color images)
[651]355ColorOrder=[1 0 0;0 1 0;0 0 1;0 0.75 0.75;0.75 0 0.75;0.75 0.75 0;0.25 0.25 0.25];
356set(hfig,'DefaultAxesColorOrder',ColorOrder)
[644]357% if isfield(Coordinates,'NextPlot')
358%     set(haxes,'NextPlot',Coordinates.NextPlot)
359% end
360if CheckHold
361     set(haxes,'NextPlot','add')
362else
363    set(haxes,'NextPlot','replace')
[201]364end
365
366%% prepare the string for plot command
367plotstr='hhh=plot(';
[644]368% coord_x_index=[];
[201]369xtitle='';
370ytitle='';
[644]371test_newplot=~CheckHold;
[569]372MinX=[];
373MaxX=[];
374MinY_cell=[];
375MaxY_cell=[];
[201]376%loop on input  fields
[530]377for icell=1:numel(CellInfo)
378    VarIndex=CellInfo{icell}.VarIndex;%  indices of the selected variables in the list data.ListVarName
379    coord_x_index=CellInfo{icell}.CoordIndex;
[201]380    testplot=ones(size(data.ListVarName));%default test for plotted variables
[426]381    coord_x_name{icell}=data.ListVarName{coord_x_index};
382    coord_x{icell}=data.(data.ListVarName{coord_x_index});%coordinate variable set as coord_x
[569]383    if isempty(find(strcmp(coord_x_name{icell},coord_x_name(1:end-1)), 1)) %xtitle not already selected
[426]384        xtitle=[xtitle coord_x_name{icell}];
385        if isfield(data,'VarAttribute')&& numel(data.VarAttribute)>=coord_x_index && isfield(data.VarAttribute{coord_x_index},'units')
386            xtitle=[xtitle '(' data.VarAttribute{coord_x_index}.units '), '];
387        else
388            xtitle=[xtitle ', '];
389        end
[201]390    end
[569]391    MinX(icell)=min(coord_x{icell});
392    MaxX(icell)=max(coord_x{icell});
[201]393    testplot(coord_x_index)=0;
[530]394    if isfield(CellInfo{icell},'VarIndex_ancillary')
395        testplot(CellInfo{icell}.VarIndex_ancillary)=0;
[201]396    end
[530]397    if isfield(CellInfo{icell},'VarIndex_warnflag')
398        testplot(CellInfo{icell}.VarIndex_warnflag)=0;
[201]399    end
400    if isfield(data,'VarAttribute')
401        VarAttribute=data.VarAttribute;
402        for ivar=1:length(VarIndex)
403            if length(VarAttribute)>=VarIndex(ivar) && isfield(VarAttribute{VarIndex(ivar)},'long_name')
404                plotname{VarIndex(ivar)}=VarAttribute{VarIndex(ivar)}.long_name;
405            else
406                plotname{VarIndex(ivar)}=data.ListVarName{VarIndex(ivar)};%name for display in plot A METTRE
407            end
408        end
409    end
[530]410    if isfield(CellInfo{icell},'VarIndex_discrete')
[201]411        charplot_0='''+''';
412    else
413        charplot_0='''-''';
414    end
[569]415    MinY=[];
416    MaxY=[];%default
417   
418    nbplot=0;
[201]419    for ivar=1:length(VarIndex)
420        if testplot(VarIndex(ivar))
421            VarName=data.ListVarName{VarIndex(ivar)};
[569]422            nbplot=nbplot+1;
[201]423            ytitle=[ytitle VarName];
424            if isfield(data,'VarAttribute')&& numel(data.VarAttribute)>=VarIndex(ivar) && isfield(data.VarAttribute{VarIndex(ivar)},'units')
425                ytitle=[ytitle '(' data.VarAttribute{VarIndex(ivar)}.units '), '];
426            else
427                ytitle=[ytitle ', '];
428            end
429            eval(['data.' VarName '=squeeze(data.' VarName ');'])
[569]430            MinY(ivar)=min(min(data.(VarName)));
431            MaxY(ivar)=max(max(data.(VarName)));
[201]432            plotstr=[plotstr 'coord_x{' num2str(icell) '},data.' VarName ',' charplot_0 ','];
433            eval(['nbcomponent2=size(data.' VarName ',2);']);
434            eval(['nbcomponent1=size(data.' VarName ',1);']);
435            if numel(coord_x{icell})==2
436                coord_x{icell}=linspace(coord_x{icell}(1),coord_x{icell}(2),nbcomponent1);
437            end
438            if nbcomponent1==1|| nbcomponent2==1
439                legend_str=[legend_str {VarName}]; %variable with one component
440            else  %variable with severals  components
441                for ic=1:min(nbcomponent1,nbcomponent2)
442                    legend_str=[legend_str [VarName '_' num2str(ic)]]; %variable with severals  components
443                end                                                   % labeled by their index (e.g. color component)
444            end
445        end
446    end
[569]447    if ~isempty(MinY)
448    MinY_cell(icell)=min(MinY);
449    MaxY_cell(icell)=max(MaxY);
450    end
[201]451end
452
453%% activate the plot
[644]454if  ~isequal(plotstr,'hhh=plot(') 
[201]455    set(hfig,'CurrentAxes',haxes)
456    tag=get(haxes,'tag');   
457    %%%
458    plotstr=[plotstr '''tag'',''plot_line'');'];   
459    eval(plotstr)                  %execute plot (instruction  plotstr)
460    %%%
[651]461    set(haxes,'tag',tag)% restitute the axes tag (removed by the command plot)
462    set(haxes,'ColorOrder',ColorOrder)% restitute the plot color order (to get red green blue for histograms or cuts of color images)
[201]463    grid(haxes, 'on')
464    hxlabel=xlabel(xtitle(1:end-2));% xlabel (removes ', ' at the end)
465    set(hxlabel,'Interpreter','none')% desable tex interpreter
466    if length(legend_str)>=1
467        hylabel=ylabel(ytitle(1:end-2));% ylabel (removes ', ' at the end)
468        set(hylabel,'Interpreter','none')% desable tex interpreter
469    end
470    if ~isempty(legend_str)
471        hlegend=findobj(hfig,'Tag','legend');
472        if isempty(hlegend)
473            hlegend=legend(legend_str);
474            txt=ver('MATLAB');
475            Release=txt.Release;
476            relnumb=str2double(Release(3:4));% should be changed to Version for better compatibility
477            if relnumb >= 14
478                set(hlegend,'Interpreter','none')% desable tex interpreter
479            end
480        else
481            legend_old=get(hlegend,'String');
482            if isequal(size(legend_old,1),size(legend_str,1))&&~isequal(legend_old,legend_str)
483                set(hlegend,'String',[legend_old legend_str]);
484            end
485        end
486    end
487    title_str='';
488    if isfield(data,'filename')
489       [Path, title_str, ext]=fileparts(data.filename);
490       title_str=[title_str ext];
491    end
[477]492    if isfield(data,'Action')&&isfield(data.Action,'ActionName')
[201]493        if ~isequal(title_str,'')
494            title_str=[title_str ', '];
495        end
[477]496        title_str=[title_str data.Action.ActionName];
[201]497    end
498    htitle=title(title_str);
[428]499    set(htitle,'Interpreter','none')% desable tex interpreter
[201]500end
501
502%% determine axes bounds
[294]503fix_lim=isfield(Coordinates,'CheckFixLimits') && Coordinates.CheckFixLimits;
[569]504check_lim=isfield(Coordinates,'MinX')&&isfield(Coordinates,'MaxX')&&isfield(Coordinates,'MinY')&&isfield(Coordinates,'MaxY');
[201]505if fix_lim
[569]506    if ~check_lim
507        fix_lim=0; %free limits if limits are not set,
[221]508    end
[201]509end
[294]510if fix_lim
511    set(haxes,'XLim',[Coordinates.MinX Coordinates.MaxX])
512    set(haxes,'YLim',[Coordinates.MinY Coordinates.MaxY])
[569]513else
514    if ~isempty(MinX)
515        if check_lim
516            CoordinatesOut.MinX=min(min(MinX),CoordinatesOut.MinX);
517            CoordinatesOut.MaxX=max(max(MaxX),CoordinatesOut.MaxX);
518        else
519            CoordinatesOut.MinX=min(MinX);
520            CoordinatesOut.MaxX=max(MaxX);
521        end
522    end
523    if ~isempty(MinY_cell)
524        if check_lim
525            CoordinatesOut.MinY=min(min(MinY_cell),CoordinatesOut.MinY);
526            CoordinatesOut.MaxY=max(max(MaxY_cell),CoordinatesOut.MaxY);
527        else
528            CoordinatesOut.MinY=min(MinY_cell);
529            CoordinatesOut.MaxY=max(MaxY_cell);
530        end
531    end
[201]532end
533
[428]534%% determine plot aspect ratio
[429]535if isfield(Coordinates,'CheckFixAspectRatio') && isequal(Coordinates.CheckFixAspectRatio,1)&&isfield(Coordinates,'AspectRatio')
[428]536    set(haxes,'DataAspectRatioMode','manual')
537    set(haxes,'DataAspectRatio',[Coordinates.AspectRatio 1 1])
538else
539    set(haxes,'DataAspectRatioMode','auto')%automatic aspect ratio
[429]540    AspectRatio=get(haxes,'DataAspectRatio');
[428]541    CoordinatesOut.AspectRatio=AspectRatio(1)/AspectRatio(2);
542end
543
[201]544%-------------------------------------------------------------------
[530]545function [haxes,PlotParamOut,PlotType,errormsg]=plot_plane(Data,CellInfo,haxes,PlotParam,PosColorbar)
[201]546%-------------------------------------------------------------------
[674]547PlotType='plane';
548grid(haxes, 'off')% remove grid (possibly remaining from other graphs)
[411]549
[201]550%default plotting parameters
551if ~isfield(PlotParam,'Scalar')
552    PlotParam.Scalar=[];
553end
554if ~isfield(PlotParam,'Vectors')
555    PlotParam.Vectors=[];
556end
[674]557PlotParamOut=PlotParam;%default
558errormsg='';%default
[201]559
[674]560hfig=get(haxes,'parent');%handle of the figure containing the plot axes
[201]561hcol=findobj(hfig,'Tag','Colorbar'); %look for colorbar axes
562hima=findobj(haxes,'Tag','ima');% search existing image in the current axes
563test_ima=0; %default: test for image or map plot
564test_vec=0; %default: test for vector plots
565test_black=0;
566test_false=0;
567test_C=0;
568XName='';
569x_units='';
570YName='';
571y_units='';
[674]572
573% loop on the input field cells
574for icell=1:numel(CellInfo)
[530]575    if strcmp(CellInfo{icell}.CoordType,'tps') %do not plot directly tps data (used for projection only)
[402]576        continue
577    end
[530]578    ivar_X=CellInfo{icell}.CoordIndex(end); % defines (unique) index for the variable representing unstructured x coordinate (default =[])
579    ivar_Y=CellInfo{icell}.CoordIndex(end-1); % defines (unique)index for the variable representing unstructured y coordinate (default =[])
580    ivar_C=[];
581    if isfield(CellInfo{icell},'VarIndex_scalar')
582        ivar_C=[ivar_C CellInfo{icell}.VarIndex_scalar];
583    end
584    if isfield(CellInfo{icell},'VarIndex_image')
585        ivar_C=[ivar_C CellInfo{icell}.VarIndex_image];
586    end
587    if isfield(CellInfo{icell},'VarIndex_color')
588        ivar_C=[ivar_C CellInfo{icell}.VarIndex_color];
589    end
590    if isfield(CellInfo{icell},'VarIndex_ancillary')
591        ivar_C=[ivar_C CellInfo{icell}.VarIndex_ancillary];
592    end
[201]593    if numel(ivar_C)>1
594        errormsg= 'error in plot_field: too many scalar inputs';
595        return
596    end
[530]597    ivar_F=[];
598    if isfield(CellInfo{icell},'VarIndex_warnflag')
[674]599        ivar_F=CellInfo{icell}.VarIndex_warnflag; %defines index (unique) for warning flag variable
[201]600    end
[674]601    ivar_FF_vec=[];
[530]602    if isfield(CellInfo{icell},'VarIndex_vector_x')&&isfield(CellInfo{icell},'VarIndex_vector_y') % vector components detected
[674]603        if test_vec% a vector field has been already detected
[515]604            errormsg='error in plot_field: attempt to plot two vector fields: to get the difference project on a plane with mode interp';
[201]605            return
606        else
607            test_vec=1;
[674]608            if isfield(CellInfo{icell},'VarIndex_errorflag')
609                ivar_FF_vec=CellInfo{icell}.VarIndex_errorflag; %defines index (unique) for error flag variable
610            end
611            vec_U=Data.(Data.ListVarName{CellInfo{icell}.VarIndex_vector_x});
[530]612            vec_V=Data.(Data.ListVarName{CellInfo{icell}.VarIndex_vector_y});
[674]613            if strcmp(CellInfo{icell}.CoordType,'scattered')%2D field with unstructured coordinates
[530]614                XName=Data.ListVarName{CellInfo{icell}.CoordIndex(end)};
615                YName=Data.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
[517]616                vec_X=reshape(Data.(XName),[],1); %transform vectors in column matlab vectors
617                vec_Y=reshape(Data.(YName),[],1);
[674]618            elseif strcmp(CellInfo{icell}.CoordType,'grid')%2D field with structured coordinates
[530]619                y=Data.(Data.ListVarName{CellInfo{icell}.CoordIndex(end-1)});
620                x=Data.(Data.ListVarName{CellInfo{icell}.CoordIndex(end)});
[201]621                if numel(y)==2 % y defined by first and last values on aregular mesh
622                    y=linspace(y(1),y(2),size(vec_U,1));
623                end
624                if numel(x)==2 % y defined by first and last values on aregular mesh
625                    x=linspace(x(1),x(2),size(vec_U,2));
626                end
[674]627                [vec_X,vec_Y]=meshgrid(x,y);
[201]628            end
[517]629            if isfield(PlotParam.Vectors,'ColorScalar') && ~isempty(PlotParam.Vectors.ColorScalar)
630                [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp([],Data,PlotParam.Vectors.ColorScalar);
[530]631                if ~isempty(VarVal)
632                    vec_C=reshape(VarVal{1},1,numel(VarVal{1}));
633                    test_C=1;
634                end
[201]635            end
[674]636            if ~isempty(ivar_F)%~(isfield(PlotParam.Vectors,'HideWarning')&& isequal(PlotParam.Vectors.HideWarning,1))
637%                 if test_vec
[313]638                    vec_F=Data.(Data.ListVarName{ivar_F}); % warning flags for  dubious vectors
[674]639                    if  ~(isfield(PlotParam.Vectors,'CheckHideWarning') && isequal(PlotParam.Vectors.CheckHideWarning,1))
[201]640                        test_black=1;
641                    end
[674]642%                 end
[201]643            end
[674]644            if ~isempty(ivar_FF_vec) %&& ~test_false
645%                 if test_vec% TODO: deal with FF for structured coordinates
646                    vec_FF=Data.(Data.ListVarName{ivar_FF_vec}); % flags for false vectors
647%                 end
[201]648            end
649        end
650    elseif ~isempty(ivar_C) %scalar or image
651        if test_ima
[674]652            errormsg='attempt to plot two scalar fields or images';
[201]653            return
654        end
[530]655        A=squeeze(Data.(Data.ListVarName{ivar_C}));% scalar represented as color image
[201]656        test_ima=1;
[674]657        if strcmp(CellInfo{icell}.CoordType,'scattered')%2D field with unstructured coordinates
[227]658            A=reshape(A,1,[]);
[201]659            XName=Data.ListVarName{ivar_X};
660            YName=Data.ListVarName{ivar_Y};
[674]661            eval(['AX=reshape(Data.' XName ',1,[]);'])
[227]662            eval(['AY=reshape(Data.' YName ',1,[]);'])
[674]663            [A,AX,AY]=proj_grid(AX',AY',A',[],[],'np>256');  % interpolate on a grid
[201]664            if isfield(Data,'VarAttribute')
665                if numel(Data.VarAttribute)>=ivar_X && isfield(Data.VarAttribute{ivar_X},'units')
666                    x_units=[' (' Data.VarAttribute{ivar_X}.units ')'];
667                end
668                if numel(Data.VarAttribute)>=ivar_Y && isfield(Data.VarAttribute{ivar_Y},'units')
669                    y_units=[' (' Data.VarAttribute{ivar_Y}.units ')'];
670                end
[674]671            end
672        elseif strcmp(CellInfo{icell}.CoordType,'grid')%2D field with structured coordinates
[530]673            YName=Data.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
[674]674            AY=Data.(YName);
[530]675            AX=Data.(Data.ListVarName{CellInfo{icell}.CoordIndex(end)});
[201]676            test_interp_X=0; %default, regularly meshed X coordinate
677            test_interp_Y=0; %default, regularly meshed Y coordinate
678            if isfield(Data,'VarAttribute')
[530]679                if numel(Data.VarAttribute)>=CellInfo{icell}.CoordIndex(end) && isfield(Data.VarAttribute{CellInfo{icell}.CoordIndex(end)},'units')
680                    x_units=Data.VarAttribute{CellInfo{icell}.CoordIndex(end)}.units;
[201]681                end
[530]682                if numel(Data.VarAttribute)>=CellInfo{icell}.CoordIndex(end-1) && isfield(Data.VarAttribute{CellInfo{icell}.CoordIndex(end-1)},'units')
683                    y_units=Data.VarAttribute{CellInfo{icell}.CoordIndex(end-1)}.units;
[201]684                end
[674]685            end
[201]686            if numel(AY)>2
687                DAY=diff(AY);
688                DAY_min=min(DAY);
689                DAY_max=max(DAY);
690                if sign(DAY_min)~=sign(DAY_max);% =1 for increasing values, 0 otherwise
[674]691                    errormsg=['errror in plot_field.m: non monotonic dimension variable ' Data.ListVarName{VarRole.coord(1)} ];
692                    return
693                end
[201]694                test_interp_Y=(DAY_max-DAY_min)> 0.0001*abs(DAY_max);
695            end
696            if numel(AX)>2
697                DAX=diff(AX);
698                DAX_min=min(DAX);
699                DAX_max=max(DAX);
700                if sign(DAX_min)~=sign(DAX_max);% =1 for increasing values, 0 otherwise
[674]701                    errormsg=['errror in plot_field.m: non monotonic dimension variable ' Data.ListVarName{VarRole.coord(2)} ];
702                    return
703                end
[201]704                test_interp_X=(DAX_max-DAX_min)> 0.0001*abs(DAX_max);
[674]705            end
706            if test_interp_Y
[201]707                npxy(1)=max([256 floor((AY(end)-AY(1))/DAY_min) floor((AY(end)-AY(1))/DAY_max)]);
708                yI=linspace(AY(1),AY(end),npxy(1));
709                if ~test_interp_X
[674]710                    xI=linspace(AX(1),AX(end),size(A,2));%default
[201]711                    AX=xI;
712                end
713            end
[674]714            if test_interp_X
[201]715                npxy(2)=max([256 floor((AX(end)-AX(1))/DAX_min) floor((AX(end)-AX(1))/DAX_max)]);
[674]716                xI=linspace(AX(1),AX(end),npxy(2));
[201]717                if ~test_interp_Y
[674]718                    yI=linspace(AY(1),AY(end),size(A,1));
719                    AY=yI;
[201]720                end
721            end
[674]722            if test_interp_X || test_interp_Y
[201]723                [AX2D,AY2D]=meshgrid(AX,AY);
724                A=interp2(AX2D,AY2D,double(A),xI,yI');
725            end
[674]726            AX=[AX(1) AX(end)];% keep only the lower and upper bounds for image represnetation
[201]727            AY=[AY(1) AY(end)];
728        end
729    end
730    %define coordinates as CoordUnits, if not defined as attribute for each variable
731    if isfield(Data,'CoordUnit')
732        if isempty(x_units)
733            x_units=Data.CoordUnit;
734        end
735        if isempty(y_units)
736            y_units=Data.CoordUnit;
737        end
[674]738    end   
739end
[688]740PlotParamOut=PlotParam; % output plot parameters equal to input by default
[201]741
742%%   image or scalar plot %%%%%%%%%%%%%%%%%%%%%%%%%%
743if test_ima
[688]744   
[201]745    % distinguish B/W and color images
746    np=size(A);%size of image
747    siz=numel(np);
748    if siz>3
[428]749        errormsg=['unrecognized scalar type: ' num2str(siz) ' dimensions'];
750        return
[201]751    end
752    if siz==3
753        if np(3)==1
754            siz=2;%B W image
755        elseif np(3)==3
756            siz=3;%color image
757        else
[206]758            errormsg=['unrecognized scalar type in plot_field: considered as 2D field with ' num2str(np(3)) ' color components'];
[201]759            return
760        end
761    end
762   
763    %set the color map
[421]764    if isfield(PlotParam.Scalar,'CheckBW') && ~isempty(PlotParam.Scalar.CheckBW)
765        BW=PlotParam.Scalar.CheckBW; %BW=0 color imposed, else gray scale imposed.
766    else % BW imposed automatically chosen
[201]767        BW=(siz==2) && (isa(A,'uint8')|| isa(A,'uint16'));% non color images represented in gray scale by default
[421]768        PlotParamOut.Scalar.CheckBW=BW;
[428]769    end
[688]770   
771            % determine the plot option 'image' or 'contours'
772    if isfield(PlotParam.Scalar,'ListContour')
773        CheckContour=strcmp(PlotParam.Scalar.ListContour,'contours');% =1 for contour plot option
774    else
775        CheckContour=0; %default
776    end
777   
[201]778    %case of grey level images or contour plot
[428]779    if siz==2
[294]780        if ~isfield(PlotParam.Scalar,'CheckFixScalar')
[688]781            PlotParam.Scalar.CheckFixScalar=0;% free scalar threshold value scale (from min to max) by default
[201]782        end
783        if ~isfield(PlotParam.Scalar,'MinA')
[688]784            PlotParam.Scalar.MinA=[];%no min scalar threshold value set
[201]785        end
786        if ~isfield(PlotParam.Scalar,'MaxA')
[688]787            PlotParam.Scalar.MaxA=[];%no max scalar threshold value set
[201]788        end
[688]789       
790        % determine the min scalar value
791        if PlotParam.Scalar.CheckFixScalar && ~isempty(PlotParam.Scalar.MinA) && isnumeric(PlotParam.Scalar.MinA) 
792            MinA=double(PlotParam.Scalar.MinA); % min value set as input
793        else
794            MinA=double(nanmin(nanmin(A))); % min value set as min of non NaN scalar values
795        end
796       
797        % error if the input scalar is NaN everywhere
798        if isnan(MinA)
[428]799                errormsg='NaN input scalar or image in plot_field';
[206]800                return
[688]801        end
802       
803        % determine the max scalar value
804        if PlotParam.Scalar.CheckFixScalar && ~isempty(PlotParam.Scalar.MaxA) && isnumeric(PlotParam.Scalar.MaxA) 
805            MaxA=double(PlotParam.Scalar.MaxA); % max value set as input
[201]806        else
[688]807            MaxA=double(nanmax(nanmax(A))); % max value set as min of non NaN scalar values
808        end
809       
[201]810        PlotParamOut.Scalar.MinA=MinA;
811        PlotParamOut.Scalar.MaxA=MaxA;
[569]812        PlotParamOut.Scalar.Npx=size(A,2);
813        PlotParamOut.Scalar.Npy=size(A,1);
[688]814       
[201]815        % case of contour plot
[294]816        if CheckContour
[201]817            if ~isempty(hima) && ishandle(hima)
[688]818                delete(hima) % delete existing image
[201]819            end
[688]820           
821            % set the contour values
[201]822            if ~isfield(PlotParam.Scalar,'IncrA')
[688]823                PlotParam.Scalar.IncrA=[];% automatic contour interval
[201]824            end
[688]825            if ~isempty(PlotParam.Scalar.IncrA) && isnumeric(PlotParam.Scalar.IncrA)
826                interval=PlotParam.Scalar.IncrA;
827            else % automatic contour interval
[201]828                cont=colbartick(MinA,MaxA);
[688]829                interval=cont(2)-cont(1);%default
830                PlotParamOut.Scalar.IncrA=interval;% set the interval as output for display on the GUI
[201]831            end
[688]832            %B=A;
833            abscontmin=interval*floor(MinA/interval);
834            abscontmax=interval*ceil(MaxA/interval);
835            contmin=interval*floor(min(min(A))/interval);
836            contmax=interval*ceil(max(max(A))/interval);
837            cont_pos_plus=0:interval:contmax;% zero and positive contour values (plotted as solid lines)
838            cont_pos_min=double(contmin):interval:-interval;% negative contour values (plotted as dashed lines)
839            cont_pos=[cont_pos_min cont_pos_plus];% set of all contour values
840           
[201]841            sizpx=(AX(end)-AX(1))/(np(2)-1);
842            sizpy=(AY(1)-AY(end))/(np(1)-1);
[428]843            x_cont=AX(1):sizpx:AX(end); % pixel x coordinates for image display
[201]844            y_cont=AY(1):-sizpy:AY(end); % pixel x coordinates for image display
[688]845           
846            %axes(haxes)% set the input axes handle as current axis
847
848           % colormap(map);
849           tag_axes=get(haxes,'Tag');% axes tag
850           Opacity=1;
851           if isfield(PlotParam.Scalar,'Opacity')&&~isempty(PlotParam.Scalar.Opacity)
852               Opacity=PlotParam.Scalar.Opacity;
853           end
854           % fill the space between contours if opacity is undefined or =1
855           if isequal(Opacity,1)
856               [var,hcontour]=contour(haxes,x_cont,y_cont,A,cont_pos);% determine all contours
857               set(hcontour,'Fill','on')% fill the space between contours
858               set(hcontour,'LineStyle','none')
859               hold on
860           end
861           [var_p,hcontour_p]=contour(haxes,x_cont,y_cont,A,cont_pos_plus,'k-');% draw the contours for positive values as solid lines
862           hold on
863           [var_m,hcontour_m]=contour(haxes,x_cont,y_cont,A,cont_pos_min,'--');% draw the contours for negative values as dashed lines
864           if isequal(Opacity,1)
865               set(hcontour_m,'LineColor',[1 1 1])% draw negative contours in white (better visibility in dark background)
866           end
867           set(haxes,'Tag',tag_axes);% restore axes tag (removed by the matlab fct contour !)
868           hold off
869           
870            %determine the color scale and map
871            caxis([abscontmin abscontmax])
872            if BW
873                vec=linspace(0,1,(abscontmax-abscontmin)/interval);%define a greyscale colormap with steps interval
[201]874                map=[vec' vec' vec'];
875                colormap(map);
[688]876            else
877                colormap('default'); % default matlab colormap ('jet')
[201]878            end
[688]879           
[428]880            if isfield(PlotParam.Coordinates,'CheckFixAspectRatio') && isequal(PlotParam.Coordinates.CheckFixAspectRatio,1)
[415]881                set(haxes,'DataAspectRatioMode','manual')
[428]882                if isfield(PlotParam.Coordinates,'AspectRatio')
[429]883                    set(haxes,'DataAspectRatio',[PlotParam.Coordinates.AspectRatio 1 1])
[428]884                else
885                    set(haxes,'DataAspectRatio',[1 1 1])
886                end
887            end
[688]888        else     
[201]889        % set  colormap for  image display
890            % rescale the grey levels with min and max, put a grey scale colorbar
[688]891%             B=A;
[201]892            if BW
893                vec=linspace(0,1,255);%define a linear greyscale colormap
894                map=[vec' vec' vec'];
[428]895                colormap(map);  %grey scale color map
[201]896            else
[688]897                colormap('default'); % standard false colors for div, vort , scalar fields
[201]898            end
899        end
900       
[428]901        % case of color images
902    else
[201]903        if BW
[688]904            A=uint16(sum(A,3));
[201]905        else
[688]906            A=uint8(A);
[201]907        end
908        MinA=0;
909        MaxA=255;
910    end
911   
912    % display usual image
[428]913    if ~CheckContour
[201]914        % interpolate field to increase resolution of image display
[652]915        test_interp=0;
[688]916        if size(A,3)==1 % scalar of B/W image
[652]917            test_interp=1;
918            if max(np) <= 64
919                npxy=8*np;% increase the resolution 8 times
920            elseif max(np) <= 128
921                npxy=4*np;% increase the resolution 4 times
922            elseif max(np) <= 256
923                npxy=2*np;% increase the resolution 2 times
924            else
925                npxy=np;
926                test_interp=0; % no interpolation done
927            end
[201]928        end
[652]929        if test_interp%if we interpolate
[201]930            x=linspace(AX(1),AX(2),np(2));
931            y=linspace(AY(1),AY(2),np(1));
932            [X,Y]=meshgrid(x,y);
933            xi=linspace(AX(1),AX(2),npxy(2));
934            yi=linspace(AY(1),AY(2),npxy(1));
[688]935            A = interp2(X,Y,double(A),xi,yi');
[428]936        end
[201]937        % create new image if there  no image handle is found
[428]938        if isempty(hima)
[201]939            tag=get(haxes,'Tag');
940            if MinA<MaxA
[688]941                hima=imagesc(AX,AY,A,[MinA MaxA]);
[201]942            else % to deal with uniform field
[688]943                hima=imagesc(AX,AY,A,[MaxA-1 MaxA]);
[201]944            end
[428]945            % the function imagesc reset the axes 'DataAspectRatioMode'='auto', change if .CheckFixAspectRatio is
[411]946            % requested:
947            set(hima,'Tag','ima')
948            set(hima,'HitTest','off')
[428]949            set(haxes,'Tag',tag);%preserve the axes tag (removed by image fct !!!)
[201]950            uistack(hima, 'bottom')
[428]951            % update an existing image
[201]952        else
[688]953            set(hima,'CData',A);
[201]954            if MinA<MaxA
955                set(haxes,'CLim',[MinA MaxA])
956            else
[238]957                set(haxes,'CLim',[MinA MaxA+1])
[201]958            end
959            set(hima,'XData',AX);
960            set(hima,'YData',AY);
961        end
[546]962       
[405]963        % set the transparency to 0.5 if vectors are also plotted
[428]964        if isfield(PlotParam.Scalar,'Opacity')&& ~isempty(PlotParam.Scalar.Opacity)
965            set(hima,'AlphaData',PlotParam.Scalar.Opacity)
[405]966        else
[428]967            if test_vec
968                set(hima,'AlphaData',0.5)%set opacity to 0.5 by default in the presence of vectors
969                PlotParamOut.Scalar.Opacity=0.5;
970            else
971                set(hima,'AlphaData',1)% full opacity (no transparency) by default
972            end
[405]973        end
[201]974    end
975    test_ima=1;
976   
977    %display the colorbar code for B/W images if Poscolorbar not empty
[546]978    if ~isempty(PosColorbar)
979        if siz==2 && exist('PosColorbar','var')
980            if isempty(hcol)||~ishandle(hcol)
981                hcol=colorbar;%create new colorbar
982            end
983            if length(PosColorbar)==4
984                set(hcol,'Position',PosColorbar)
985            end
986            %YTick=0;%default
987            if MaxA>MinA
988                if CheckContour
989                    colbarlim=get(hcol,'YLim');
990                    scale_bar=(colbarlim(2)-colbarlim(1))/(abscontmax-abscontmin);
991                    YTick=cont_pos(2:end-1);
992                    YTick_scaled=colbarlim(1)+scale_bar*(YTick-abscontmin);
993                    set(hcol,'YTick',YTick_scaled);
994                elseif (isfield(PlotParam.Scalar,'CheckBW') && isequal(PlotParam.Scalar.CheckBW,1))||isa(A,'uint8')|| isa(A,'uint16')%images
995                    hi=get(hcol,'children');
996                    if iscell(hi)%multiple images in colorbar
997                        hi=hi{1};
998                    end
999                    set(hi,'YData',[MinA MaxA])
1000                    set(hi,'CData',(1:256)')
1001                    set(hcol,'YLim',[MinA MaxA])
1002                    YTick=colbartick(MinA,MaxA);
1003                    set(hcol,'YTick',YTick)
1004                else
1005                    hi=get(hcol,'children');
1006                    if iscell(hi)%multiple images in colorbar
1007                        hi=hi{1};
1008                    end
1009                    set(hi,'YData',[MinA MaxA])
1010                    set(hi,'CData',(1:64)')
1011                    YTick=colbartick(MinA,MaxA);
1012                    set(hcol,'YLim',[MinA MaxA])
1013                    set(hcol,'YTick',YTick)
[201]1014                end
[546]1015                set(hcol,'Yticklabel',num2str(YTick'));
[201]1016            end
[546]1017        elseif ishandle(hcol)
1018            delete(hcol); %erase existing colorbar if not needed
[201]1019        end
1020    end
1021else%no scalar plot
[428]1022    if ~isempty(hima) && ishandle(hima)
[201]1023        delete(hima)
1024    end
[546]1025    if ~isempty(PosColorbar) && ~isempty(hcol)&& ishandle(hcol)
[428]1026        delete(hcol)
[201]1027    end
1028    PlotParamOut=rmfield(PlotParamOut,'Scalar');
1029end
1030
1031%%   vector plot %%%%%%%%%%%%%%%%%%%%%%%%%%
1032if test_vec
1033   %vector scale representation
1034    if size(vec_U,1)==numel(vec_Y) && size(vec_U,2)==numel(vec_X); % x, y  coordinate variables
1035        [vec_X,vec_Y]=meshgrid(vec_X,vec_Y);
1036    end   
1037    vec_X=reshape(vec_X,1,numel(vec_X));%reshape in matlab vectors
1038    vec_Y=reshape(vec_Y,1,numel(vec_Y));
1039    vec_U=reshape(vec_U,1,numel(vec_U));
1040    vec_V=reshape(vec_V,1,numel(vec_V));
1041     MinMaxX=max(vec_X)-min(vec_X);
[294]1042    if  isfield(PlotParam.Vectors,'CheckFixVectors') && isequal(PlotParam.Vectors.CheckFixVectors,1)&& isfield(PlotParam.Vectors,'VecScale')...
[210]1043               &&~isempty(PlotParam.Vectors.VecScale) && isa(PlotParam.Vectors.VecScale,'double') %fixed vector scale
1044        scale=PlotParam.Vectors.VecScale;  %impose the length of vector representation
1045    else
1046        if ~test_false %remove false vectors   
[201]1047            indsel=1:numel(vec_X);%
1048        end
1049        if isempty(vec_U)
1050            scale=1;
1051        else
1052            if isempty(indsel)
1053                MaxU=max(abs(vec_U));
1054                MaxV=max(abs(vec_V));
1055            else
1056                MaxU=max(abs(vec_U(indsel)));
1057                MaxV=max(abs(vec_V(indsel)));
1058            end
1059            scale=MinMaxX/(max(MaxU,MaxV)*50);
1060            PlotParam.Vectors.VecScale=scale;%update the 'scale' display
1061        end
[210]1062    end
[201]1063   
1064    %record vectors on the plotting axes
1065    if test_C==0
1066        vec_C=ones(1,numel(vec_X));
1067    end
1068   
1069    %decimate by a factor 2 in vector mesh(4 in nbre of vectors)
[581]1070    check_decimate=0;
[294]1071    if isfield(PlotParam.Vectors,'CheckDecimate4') && PlotParam.Vectors.CheckDecimate4
[581]1072        check_decimate=1;
[201]1073        diffy=diff(vec_Y); %difference dy=vec_Y(i+1)-vec_Y(i)
1074        dy_thresh=max(abs(diffy))/2;
1075        ind_jump=find(abs(diffy) > dy_thresh); %indices with diff(vec_Y)> max/2, detect change of line
1076        ind_sel=1:ind_jump(1);%select the first line
1077        for i=2:2:length(ind_jump)-1
1078            ind_sel=[ind_sel (ind_jump(i)+1:ind_jump(i+1))];% select the odd lines
1079        end
1080        nb_sel=length(ind_sel);
1081        ind_sel=ind_sel(1:2:nb_sel);% take half the points on a line
[581]1082    elseif isfield(PlotParam.Vectors,'CheckDecimate16') && PlotParam.Vectors.CheckDecimate16
1083        check_decimate=1;
1084        diffy=diff(vec_Y); %difference dy=vec_Y(i+1)-vec_Y(i)
1085        dy_thresh=max(abs(diffy))/2;
1086        ind_jump=find(abs(diffy) > dy_thresh); %indices with diff(vec_Y)> max/2, detect change of line
1087        ind_sel=1:ind_jump(1);%select the first line
1088        for i=2:4:length(ind_jump)-1
1089            ind_sel=[ind_sel (ind_jump(i)+1:ind_jump(i+1))];% select the odd lines
1090        end
1091        nb_sel=length(ind_sel);
1092        ind_sel=ind_sel(1:4:nb_sel);% take half the points on a line
1093    end
1094    if check_decimate
[201]1095        vec_X=vec_X(ind_sel);
1096        vec_Y=vec_Y(ind_sel);
1097        vec_U=vec_U(ind_sel);
1098        vec_V=vec_V(ind_sel);
1099        vec_C=vec_C(ind_sel);
1100        if ~isempty(ivar_F)
1101           vec_F=vec_F(ind_sel);
1102        end
[674]1103        if ~isempty(ivar_FF_vec)
[201]1104           vec_FF=vec_FF(ind_sel);
1105        end
1106    end
1107   
1108    %get main level color code
1109    [colorlist,col_vec,PlotParamOut.Vectors]=set_col_vec(PlotParam.Vectors,vec_C);
1110   
1111    % take flags into account: add flag colors to the list of colors
1112    sizlist=size(colorlist);
1113    nbcolor=sizlist(1);
1114    if test_black
1115       nbcolor=nbcolor+1;
1116       colorlist(nbcolor,:)=[0 0 0]; %add black to the list of colors
[674]1117       if ~isempty(ivar_FF_vec)
[210]1118            col_vec(vec_F~=1 & vec_F~=0 & vec_FF==0)=nbcolor;
[201]1119       else
[210]1120            col_vec(vec_F~=1 & vec_F~=0)=nbcolor;
[201]1121       end
1122    end
1123    nbcolor=nbcolor+1;
[674]1124    if ~isempty(ivar_FF_vec)
[294]1125        if isfield(PlotParam.Vectors,'CheckHideFalse') && PlotParam.Vectors.CheckHideFalse==1
[201]1126            colorlist(nbcolor,:)=[NaN NaN NaN];% no plot of false vectors
1127        else
1128            colorlist(nbcolor,:)=[1 0 1];% magenta color
1129        end
[210]1130        col_vec(vec_FF~=0)=nbcolor;
[201]1131    end
1132    %plot vectors:
1133    quiresetn(haxes,vec_X,vec_Y,vec_U,vec_V,scale,colorlist,col_vec);   
1134
1135else
1136    hvec=findobj(haxes,'Tag','vel');
1137    if ~isempty(hvec)
1138        delete(hvec);
1139    end
1140    PlotParamOut=rmfield(PlotParamOut,'Vectors');
1141end
[546]1142% nbvar=0;
[201]1143
1144%store the coordinate extrema occupied by the field
1145if ~isempty(Data)
[569]1146    MinX=[];
1147    MaxX=[];
1148    MinY=[];
1149    MaxY=[];
[334]1150    fix_lim=isfield(PlotParam.Coordinates,'CheckFixLimits') && PlotParam.Coordinates.CheckFixLimits;
[201]1151    if fix_lim
[334]1152        if isfield(PlotParam.Coordinates,'MinX')&&isfield(PlotParam.Coordinates,'MaxX')&&isfield(PlotParam.Coordinates,'MinY')&&isfield(PlotParam.Coordinates,'MaxY')
[569]1153            MinX=PlotParam.Coordinates.MinX;
1154            MaxX=PlotParam.Coordinates.MaxX;
1155            MinY=PlotParam.Coordinates.MinY;
1156            MaxY=PlotParam.Coordinates.MaxY;
1157        end  %else PlotParamOut.MinX =PlotParam.MinX...
[252]1158    else
1159        if test_ima %both background image and vectors coexist, take the wider bound
[569]1160            MinX=min(AX);
1161            MaxX=max(AX);
1162            MinY=min(AY);
1163            MaxY=max(AY);
[252]1164            if test_vec
[569]1165                MinX=min(MinX,min(vec_X));
1166                MaxX=max(MaxX,max(vec_X));
1167                MinY=min(MinY,min(vec_Y));
1168                MaxY=max(MaxY,max(vec_Y));
[252]1169            end
1170        elseif test_vec
[569]1171            MinX=min(vec_X);
1172            MaxX=max(vec_X);
1173            MinY=min(vec_Y);
1174            MaxY=max(vec_Y);
[201]1175        end
[252]1176    end
[569]1177    PlotParamOut.Coordinates.MinX=MinX;
1178    PlotParamOut.Coordinates.MaxX=MaxX;
1179    PlotParamOut.Coordinates.MinY=MinY;
1180    PlotParamOut.Coordinates.MaxY=MaxY;
1181    if MaxX>MinX
1182        set(haxes,'XLim',[MinX MaxX]);% set x limits of frame in axes coordinates
[546]1183    end
[569]1184    if MaxY>MinY
1185        set(haxes,'YLim',[MinY MaxY]);% set x limits of frame in axes coordinates
[546]1186    end
[201]1187    set(haxes,'YDir','normal')
1188    set(get(haxes,'XLabel'),'String',[XName ' (' x_units ')']);
1189    set(get(haxes,'YLabel'),'String',[YName ' (' y_units ')']);
[334]1190    PlotParamOut.Coordinates.x_units=x_units;
1191    PlotParamOut.Coordinates.y_units=y_units;
[201]1192end
[546]1193if isfield(PlotParam,'Coordinates') && isfield(PlotParam.Coordinates,'CheckFixAspectRatio') && isequal(PlotParam.Coordinates.CheckFixAspectRatio,1)
1194    set(haxes,'DataAspectRatioMode','manual')
1195    if isfield(PlotParam.Coordinates,'AspectRatio')
1196        set(haxes,'DataAspectRatio',[PlotParam.Coordinates.AspectRatio 1 1])
1197    end
1198else
1199    set(haxes,'DataAspectRatioMode','auto')
1200end
[201]1201%-------------------------------------------------------------------
1202% --- function for plotting vectors
1203%INPUT:
1204% haxes: handles of the plotting axes
1205% x,y,u,v: vectors coordinates and vector components to plot, arrays withb the same dimension
1206% scale: scaling factor for vector length representation
1207% colorlist(icolor,:): list of vector colors, dim (nbcolor,3), depending on color #i
1208% col_vec: matlab vector setting the color number #i for each velocity vector
1209function quiresetn(haxes,x,y,u,v,scale,colorlist,col_vec)
1210%-------------------------------------------------------------------
1211%define arrows
1212theta=0.5 ;%angle arrow
1213alpha=0.3 ;%length arrow
1214rot=alpha*[cos(theta) -sin(theta); sin(theta) cos(theta)]';
1215%find the existing lines
1216h=findobj(haxes,'Tag','vel');% search existing lines in the current axes
1217sizh=size(h);
1218set(h,'EraseMode','xor');
1219set(haxes,'NextPlot','replacechildren');
1220
1221%drawnow
1222%create lines (if no lines) or modify them
1223if ~isequal(size(col_vec),size(x))
1224    col_vec=ones(size(x));% case of error in col_vec input
1225end
1226sizlist=size(colorlist);
1227ncolor=sizlist(1);
1228
1229for icolor=1:ncolor
1230    %determine the line positions for each color icolor
1231    ind=find(col_vec==icolor);
1232    xc=x(ind);
1233    yc=y(ind);
1234    uc=u(ind)*scale;
1235    vc=v(ind)*scale;
1236    n=size(xc);
1237    xN=NaN*ones(size(xc));
1238    matx=[xc(:)-uc(:)/2 xc(:)+uc(:)/2 xN(:)]';
1239    matx=reshape(matx,1,3*n(2));
1240    maty=[yc(:)-vc(:)/2 yc(:)+vc(:)/2 xN(:)]';
1241    maty=reshape(maty,1,3*n(2));
1242   
1243    %determine arrow heads
1244    arrowplus=rot*[uc;vc];
1245    arrowmoins=rot'*[uc;vc];
1246    x1=xc+uc/2-arrowplus(1,:);
1247    x2=xc+uc/2;
1248    x3=xc+uc/2-arrowmoins(1,:);
1249    y1=yc+vc/2-arrowplus(2,:);
1250    y2=yc+vc/2;
1251    y3=yc+vc/2-arrowmoins(2,:);
1252    matxar=[x1(:) x2(:) x3(:) xN(:)]';
1253    matxar=reshape(matxar,1,4*n(2));
1254    matyar=[y1(:) y2(:) y3(:) xN(:)]';
1255    matyar=reshape(matyar,1,4*n(2));
1256    %draw the line or modify the existing ones
[421]1257%     tri=reshape(1:3*length(uc),3,[])';   
[201]1258    isn=isnan(colorlist(icolor,:));%test if color NaN
1259    if 2*icolor > sizh(1) %if icolor exceeds the number of existing ones
1260        if ~isn(1) %if the vectors are visible color not nan
1261            if n(2)>0
1262                hold on
1263                line(matx,maty,'Color',colorlist(icolor,:),'Tag','vel');% plot new lines
1264                line(matxar,matyar,'Color',colorlist(icolor,:),'Tag','vel');% plot arrows
1265            end
1266        end
1267    else
1268        if isn(1)
1269            delete(h(2*icolor-1))
1270            delete(h(2*icolor))
1271        else
1272            set(h(2*icolor-1),'Xdata',matx,'Ydata',maty);
1273            set(h(2*icolor-1),'Color',colorlist(icolor,:));
1274            set(h(2*icolor-1),'EraseMode','xor');
1275            set(h(2*icolor),'Xdata',matxar,'Ydata',matyar);
1276            set(h(2*icolor),'Color',colorlist(icolor,:));
1277            set(h(2*icolor),'EraseMode','xor');
1278        end
1279    end
1280end
1281if sizh(1) > 2*ncolor
1282    for icolor=ncolor+1 : sizh(1)/2%delete additional objects
1283        delete(h(2*icolor-1))
1284        delete(h(2*icolor))
1285    end
1286end
1287
1288%-------------------------------------------------------------------
1289% ---- determine tick positions for colorbar
1290function YTick=colbartick(MinA,MaxA)
1291%-------------------------------------------------------------------
1292%determine tick positions with "simple" values between MinA and MaxA
1293YTick=0;%default
1294maxabs=max([abs(MinA) abs(MaxA)]);
1295if maxabs>0
1296ord=10^(floor(log10(maxabs)));%order of magnitude
1297div=1;
1298siz2=1;
1299while siz2<2
1300    values=-10:div:10;
1301    ind=find((ord*values-MaxA)<0 & (ord*values-MinA)>0);%indices of 'values' such that MinA<ord*values<MaxA
1302    siz=size(ind);
1303    if siz(2)<4%if there are less than 4 selected values (4 levels)
1304        values=-9:0.5*div:9;
1305        ind=find((ord*values-MaxA)<0 & (ord*values-MinA)>0);
1306    end
1307    siz2=size(ind,2);
1308    div=div/10;
1309end
1310YTick=ord*values(ind);
1311end
1312
[356]1313% -------------------------------------------------------------------------
[389]1314% --- 'proj_grid': project  fields with unstructured coordinantes on a regular grid
[356]1315function [A,rangx,rangy]=proj_grid(vec_X,vec_Y,vec_A,rgx_in,rgy_in,npxy_in)
[389]1316% -------------------------------------------------------------------------
[356]1317if length(vec_Y)<2
1318    msgbox_uvmat('ERROR','less than 2 points in proj_grid.m');
1319    return;
1320end
1321diffy=diff(vec_Y); %difference dy=vec_Y(i+1)-vec_Y(i)
1322index=find(diffy);% find the indices of vec_Y after wich a change of horizontal line occurs(diffy non zero)
1323if isempty(index); msgbox_uvmat('ERROR','points aligned along abscissa in proj_grid.m'); return; end;%points aligned% A FAIRE: switch to line plot.
1324diff2=diff(diffy(index));% diff2 = fluctuations of the detected vertical grid mesh dy
1325if max(abs(diff2))>0.001*abs(diffy(index(1))) % if max(diff2) is larger than 1/1000 of the first mesh dy
1326    % the data are not regularly spaced and must be interpolated  on a regular grid
1327    if exist('rgx_in','var') & ~isempty (rgx_in) & isnumeric(rgx_in) & length(rgx_in)==2%  positions imposed from input
1328        rangx=rgx_in; % first and last positions
1329        rangy=rgy_in;
1330        dxy(1)=1/(npxy_in(1)-1);%grid mesh in y
1331        dxy(2)=1/(npxy_in(2)-1);%grid mesh in x
1332        dxy(1)=(rangy(2)-rangy(1))/(npxy_in(1)-1);%grid mesh in y
1333        dxy(2)=(rangx(2)-rangx(1))/(npxy_in(2)-1);%grid mesh in x
1334    else % interpolation grid automatically determined
1335        rangx(1)=min(vec_X);
1336        rangx(2)=max(vec_X);
1337        rangy(2)=min(vec_Y);
1338        rangy(1)=max(vec_Y);
1339        dxymod=sqrt((rangx(2)-rangx(1))*(rangy(1)-rangy(2))/length(vec_X));
1340        dxy=[-dxymod/4 dxymod/4];% increase the resolution 4 times
1341    end
1342    xi=[rangx(1):dxy(2):rangx(2)];
1343    yi=[rangy(1):dxy(1):rangy(2)];
1344    A=griddata_uvmat(vec_X,vec_Y,vec_A,xi,yi');
1345    A=reshape(A,length(yi),length(xi));
1346else
1347    x=vec_X(1:index(1));% the set of abscissa (obtained on the first line)
1348    indexend=index(end);% last vector index of line change
1349    ymax=vec_Y(indexend+1);% y coordinate AFTER line change
1350    ymin=vec_Y(index(1));
1351    y=vec_Y(index);
1352    y(length(y)+1)=ymax;
1353    nx=length(x);   %number of grid points in x
1354    ny=length(y);   % number of grid points in y
1355    B=(reshape(vec_A,nx,ny))'; %vec_A reshaped as a rectangular matrix
1356    [X,Y]=meshgrid(x,y);% positions X and Y also reshaped as matrix
1357
1358    %linear interpolation to improve the image resolution and/or adjust
1359    %to prescribed positions
1360    test_interp=1;
1361    if exist('rgx_in','var') & ~isempty (rgx_in) & isnumeric(rgx_in) & length(rgx_in)==2%  positions imposed from input
1362        rangx=rgx_in; % first and last positions
1363        rangy=rgy_in;
1364        npxy=npxy_in;
1365    else       
1366        rangx=[vec_X(1) vec_X(nx)];% first and last position found for x
1367          rangy=[max(ymax,ymin) min(ymax,ymin)];
1368        if max(nx,ny) <= 64 & isequal(npxy_in,'np>256')
1369            npxy=[8*ny 8*nx];% increase the resolution 8 times
1370        elseif max(nx,ny) <= 128 & isequal(npxy_in,'np>256')
1371            npxy=[4*ny 4*nx];% increase the resolution 4 times
1372        elseif max(nx,ny) <= 256 & isequal(npxy_in,'np>256')
1373            npxy=[2*ny 2*nx];% increase the resolution 2 times
1374        else
1375            npxy=[ny nx];
1376            test_interp=0; % no interpolation done
1377        end
1378    end
1379    if test_interp==1%if we interpolate
1380        xi=[rangx(1):(rangx(2)-rangx(1))/(npxy(2)-1):rangx(2)];
1381        yi=[rangy(1):(rangy(2)-rangy(1))/(npxy(1)-1):rangy(2)];
1382        [XI,YI]=meshgrid(xi,yi);
1383        A = interp2(X,Y,B,XI,YI);
1384    else %no interpolation for a resolution higher than 256
1385        A=B;
1386    end
1387end
Note: See TracBrowser for help on using the repository browser.