source: trunk/src/plot_field.m @ 210

Last change on this file since 210 was 210, checked in by sommeria, 10 years ago

correction bug for vector plot (Fixscale did not work), extensions for 3D fields

File size: 50.1 KB
Line 
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:
6%     2D vector fields are represented by arrows, 2D scalar fiedlds by grey scale images or contour plots, 1D fields are represented by usual plot with (abscissa, ordinate).
7%  The input field structure is first tested by check_field_structure.m,
8%  then split into blocks of related variables  by find_field_indices.m.
9%  The dimensionality of each block is obtained  by this fuction
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%
14% function [PlotType,PlotParamOut,haxes]= plot_field(Data,haxes,PlotParam,htext,PosColorbar)
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%
44%         additional elements characterizing the projection object (should not be necessary)--
45%            Data.Style : style of projection object
46%            Data.XObject,.YObject: set of coordinates defining the object position;
47%            Data.ProjMode=type of projection ;
48%            Data.ProjAngle=angle of projection;
49%            Data.DX,.DY,.DZ=increments;
50%            Data.MaxY,MinY: min and max Y
51
52%   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.
53%
54%   PlotParam: parameters for plotting, as read on the uvmat interface (by function 'read_plot_param.m')
55%     .FixLimits:=0 (default) adjust axes limit to the X,Y data, =1: preserves the previous axes limits
56%     .FixEqual: =0 (default):automatic adjustment of the graph, keep 1 to 1 aspect ratio for x and y scales.
57%            --scalars--
58%    .Scalar.MaxA: upper bound (saturation color) for the scalar representation, max(field) by default
59%    .Scalar.MinA: lower bound (saturation) for the scalar representation, min(field) by default
60%    .Scalar.FixScal: =0 (default) lower and upper bounds of the scalar representation set to the min and max of the field
61%               =1 lower and upper bound imposed by .AMax and .MinA
62%    .Scalar.BW= 1 black and white representation imposed, =0 by default.
63%    .Scalar.Contours= 1: represent scalars by contour plots (Matlab function 'contour'); =0 by default
64%    .IncrA : contour interval
65%            -- vectors--
66%    .Vectors.VecScale: scale for the vector representation
67%    .Vectors.FixVec: =0 (default) automatic length for vector representation, =1: length set by .VecScale
68%    .Vectors.HideFalse= 0 (default) false vectors represented in magenta, =1: false vectors not represented;
69%    .Vectors.HideWarning= 0 (default) vectors marked by warnflag~=0 marked in black, 1: no warning representation;
70%    .Vectors.decimate4 = 0 (default) all vectors reprtesented, =1: half of  the vectors represented along each coordinate
71%         -- vector color--
72%    .Vectors.ColorCode= 'black','white': imposed color  (default ='blue')
73%                        'rgb', : three colors red, blue, green depending
74%                        on thresholds .colcode1 and .colcode2 on the input  scalar value (C)
75%                        'brg': like rgb but reversed color order (blue, green, red)
76%                        '64 colors': continuous color from blue to red (multijet)
77%    .Vectors.colcode1 : first threshold for rgb, first value for'continuous'
78%    .Vectors.colcode2 : second threshold for rgb, last value (saturation) for 'continuous'
79%    .Vectors.FixedCbounds;  =0 (default): the bounds on C representation are min and max, =1: they are fixed by .Minc and .MaxC
80%    .Vectors.MinC = imposed minimum of the scalar field used for vector color;
81%    .Vectors.MaxC = imposed maximum of the scalar field used for vector color;
82%
83%
84% PosColorbar: if not empty, display a colorbar for B&W images
85%               imposed position of the colorbar (ex [0.821 0.471 0.019 0.445])
86
87%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
88%  Copyright Joel Sommeria, 2008, LEGI / CNRS-UJF-INPG, sommeria@coriolis-legi.org.
89%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
90%     This file is part of the toolbox UVMAT.
91%
92%     UVMAT is free software; you can redistribute it and/or modify
93%     it under the terms of the GNU General Public License as published by
94%     the Free Software Foundation; either version 2 of the License, or
95%     (at your option) any later version.
96%
97%     UVMAT is distributed in the hope that it will be useful,
98%     but WITHOUT ANY WARRANTY; without even the implied warranty of
99%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
100%     GNU General Public License (file UVMAT/COPYING.txt) for more details.
101%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
102
103function [PlotType,PlotParamOut,haxes]= plot_field(Data,haxes,PlotParam,PosColorbar)
104% TODO:
105% use htext: handles of the text edit box (uicontrol)
106% introduce PlotParam.Hold: 'on' or 'off' (for curves)
107%default output
108if ~exist('PlotParam','var'),PlotParam=[];end;
109if ~exist('PosColorbar','var'),PosColorbar=[];end;
110PlotType='text'; %default
111PlotParamOut=PlotParam;%default
112
113%% test axes and figure
114testnewfig=1;%test to create a new figure (default)
115testzoomaxes=0;%test for the existence of a zoom secondary figure attached to the plotting axes
116if exist('haxes','var')
117    if ishandle(haxes)
118        if isequal(get(haxes,'Type'),'axes')
119            testnewfig=0;
120            AxeData=get(haxes,'UserData');
121            if isfield(AxeData,'ZoomAxes')&& ishandle(AxeData.ZoomAxes)
122                if isequal(get(AxeData.ZoomAxes,'Type'),'axes')
123                    testzoomaxes=1;
124                    zoomaxes=AxeData.ZoomAxes;
125                end
126            end
127        end
128    end
129end
130% if isfield(PlotParam,'text_display_1') && ishandle(PlotParam.text_display_1)
131%     PlotParam=read_plot_param(PlotParam);   
132% end
133% create a new figure and axes if the plotting axes does not exist
134if testnewfig
135    hfig=figure;
136    set(hfig,'Units','normalized')
137    haxes=axes;
138    set(haxes,'position',[0.13,0.2,0.775,0.73])
139     PlotParam.NextPlot='add'; %parameter for plot_profile and plot_his
140else
141    hfig=get(haxes,'parent');
142    set(0,'CurrentFigure',hfig)% the parent of haxes becomes the current figure
143    set(hfig,'CurrentAxes',haxes)%  haxes becomes the current axes of the parent figure
144end
145
146%% check input structure
147if ~isempty(Data)
148    [Data,errormsg]=check_field_structure(Data);
149    if ~isempty(errormsg)
150        msgbox_uvmat('ERROR',['input of plot_field/check_field_structure: ' errormsg])
151        display(['input of plot_field/check_field_structure:: ' errormsg])
152        return
153    end
154
155    %% check the cells of fields :
156    [CellVarIndex,NbDim,VarType,errormsg]=find_field_indices(Data);
157    if ~isempty(errormsg)
158        msgbox_uvmat('ERROR',['input of plot_field/find_field_indices: ' errormsg]);
159        return
160    end
161    index_2D=find(NbDim==2,2);%find 2D fields (at most 2)
162    index_3D=find(NbDim>2,1);
163    if ~isempty(index_3D)
164        if isfield(Data,'NbDim')&& isequal(Data.NbDim,2)
165            index_2D=[index_2D index_3D];
166        else
167        msgbox_uvmat('ERROR','volume plot not implemented yet');
168        return
169        end
170    end
171   
172    index_1D=find(NbDim==1);
173    index_0D=find(NbDim==0);
174   
175        %% set axes properties
176    if isfield(PlotParam,'FixLimits') && isequal(PlotParam.FixLimits,1)  %adjust the graph limits*
177        set(haxes,'XLimMode', 'manual')
178        set(haxes,'YLimMode', 'manual')
179    else
180        set(haxes,'XLimMode', 'auto')
181        set(haxes,'YLimMode', 'auto')
182    end
183    if ~isfield(PlotParam,'FixEqual')&& isfield(Data,'CoordUnit')
184        PlotParam.FixEqual=1;
185    end
186    if isfield(PlotParam,'FixEqual') && isequal(PlotParam.FixEqual,1)
187        set(haxes,'DataAspectRatioMode','manual')
188       set(haxes,'DataAspectRatio',[1 1 1])
189    else
190        set(haxes,'DataAspectRatioMode','auto')%automatic aspect ratio
191    end
192else
193    index_2D=[];
194    index_1D=[];
195    index_0D=[];
196end
197
198%% plot if the input field is valid
199PlotType='text';
200errormsg=[];
201AxeData=get(haxes,'UserData');
202if isempty(index_2D)
203    plot_plane([],[],[],haxes);%removes images or vector plots if any
204else
205    [xx,PlotParamOut,PlotType,errormsg]=plot_plane(Data,CellVarIndex(index_2D),VarType(index_2D),haxes,PlotParam,PosColorbar);
206    AxeData.NbDim=2;
207    if testzoomaxes && isempty(errormsg)
208        [zoomaxes,PlotParamOut,xx,errormsg]=plot_plane(Data,CellVarIndex(index_2D),VarType(index_2D),zoomaxes,PlotParam,PosColorbar);
209        AxeData.ZoomAxes=zoomaxes;
210    end
211    %remove coordinates variables from 1D plot
212    for ivar=1:length(index_1D)
213        if isequal(CellVarIndex{index_1D(ivar)},VarType{index_1D(ivar)}.coord)
214            index_1D(ivar)=0;
215        end
216    end
217    index_1D=find(index_1D);
218end
219
220if isempty(index_1D)
221     plot_profile([],[],[],haxes);%
222else
223        PlotParamOut=plot_profile(Data,CellVarIndex(index_1D),VarType(index_1D),haxes,PlotParam);%
224        if testzoomaxes
225            [zoomaxes,PlotParamOut]=plot_profile(Data,CellVarIndex(index_1D),VarType(index_1D),zoomaxes,PlotParam);
226            AxeData.ZoomAxes=zoomaxes;
227        end
228        PlotType='line';
229end
230htext=findobj(hfig,'Tag','text_display');
231if ~isempty(htext)
232    if isempty(index_0D)
233        set(htext,'String',{''})
234    else
235        [errormsg]=plot_text(Data,CellVarIndex(index_0D),htext);
236    end
237end
238
239if ~isempty(errormsg)
240    msgbox_uvmat('ERROR', errormsg)
241end
242if isfield(PlotParamOut,'MinX')
243    set(haxes,'XLim',[PlotParamOut.MinX PlotParamOut.MaxX])
244    set(haxes,'YLim',[PlotParamOut.MinY PlotParamOut.MaxY])
245    AxeData.RangeX=[PlotParamOut.MinX PlotParamOut.MaxX];
246    AxeData.RangeY=[PlotParamOut.MinY PlotParamOut.MaxY];
247end
248
249%% update the parameters stored in AxeData
250set(haxes,'UserData',AxeData)
251
252%% update the plotted field stored in parent figure
253FigData=get(hfig,'UserData');
254tagaxes=get(haxes,'tag');
255if isfield(FigData,tagaxes)
256    eval(['FigData.' tagaxes '=Data;'])
257    set(hfig,'UserData',FigData)
258end
259             
260%-------------------------------------------------------------------
261function errormsg=plot_text(FieldData,CellVarIndex,htext)
262%-------------------------------------------------------------------
263% if exist('hdisplay_in','var') && ~isempty(hdisplay_in) && ishandle(hdisplay_in) && isequal(get(hdisplay_in,'Type'),'uicontrol')
264%     hdisplay=hdisplay_in;
265% else
266%     figure;%create new figure
267%     hdisplay=uicontrol('Style','edit', 'Units','normalized','Position', [0 0 1 1],'Max',2,'FontName','monospaced');
268% end
269errormsg=[];
270txt_cell={};
271for icell=1:length(CellVarIndex)
272    VarIndex=CellVarIndex{icell};%  indices of the selected variables in the list data.ListVarName
273% ff=fields(FieldData);%list of field names
274% vv=struct2cell(FieldData);%list of field values
275%
276% for icell=1:length(vv)
277    for ivar=1:length(VarIndex)
278         VarName=FieldData.ListVarName{VarIndex(ivar)};
279         eval(['VarValue=FieldData.' VarName ';'])
280         if size(VarValue,1)~=1
281             VarValue=VarValue';
282         end
283         txt=[VarName '=' num2str(VarValue)];
284         txt_cell=[txt_cell;{txt}];
285    end
286end
287
288set(htext,'String',txt_cell)
289%     txt_cell=[txt_cell {num2str(
290%     Tabcell{icell,1}=ff{icell};
291%     ss=vv{icell};
292%     sizss=size(ss);
293%     if isnumeric(ss)
294%         if sizss(1)<=1 && length(ss)<5
295%             displ{icell}=num2str(ss);
296%         else
297%             displ{icell}=[class(ss) ', size ' num2str(size(ss))];
298%         end
299%     elseif ischar(ss)
300%         displ{icell}=ss;
301%     elseif iscell(ss)
302%         sizcell=size(ss);
303%         if sizcell(1)==1 && length(sizcell)==2 %line cell
304%            ssline='{''';
305%            for icolumn=1:sizcell(2)
306%                if isnumeric(ss{icolumn})
307%                    if size(ss{icolumn},1)<=1 && length(ss{icolumn})<5
308%                       sscolumn=num2str(ss{icolumn});%line vector
309%                    else
310%                       sscolumn=[class(ss{icolumn}) ', size ' num2str(size(ss{icolumn}))];
311%                    end
312%                elseif ischar(ss{icolumn})
313%                    sscolumn=ss{icolumn};
314%                else
315%                    sscolumn=class(ss{icolumn});
316%                end
317%                if icolumn==1
318%                    ssline=[ssline sscolumn];
319%                else
320%                    ssline=[ssline ''',''' sscolumn];
321%                end
322%            end
323%            displ{icell}=[ssline '''}'];
324%         else
325%            displ{icell}=[class(ss) ', size ' num2str(sizcell)];
326%         end
327%     else
328%         displ{icell}=class(ss);
329%     end
330%     Tabcell{icell,2}=displ{icell};
331% end
332% Tabchar=cell2tab(Tabcell,': ');
333% set(hdisplay,'String', Tabchar)
334
335
336%-------------------------------------------------------------------
337function PlotParamOut=plot_profile(data,CellVarIndex,VarType,haxes,PlotParam)
338%-------------------------------------------------------------------
339
340if ~exist('PlotParam','var')
341    PlotParam=[];
342end
343PlotParamOut=PlotParam; %default
344hfig=get(haxes,'parent');
345%suppress existing plot isf empty data
346if isempty(data)
347    hplot=findobj(haxes,'tag','plot_line');
348    if ~isempty(hplot)
349        delete(hplot)
350    end
351    hlegend=findobj(hfig,'tag','legend');
352    if ~isempty(hlegend)
353        delete(hlegend)
354    end
355    return
356end
357
358
359ColorOrder=[1 0 0;0 0.5 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];
360set(haxes,'ColorOrder',ColorOrder)
361if isfield(PlotParam,'NextPlot')
362    set(haxes,'NextPlot',PlotParam.NextPlot)
363end
364% adjust the size of the plot to include the whole field,
365
366legend_str={};
367
368%% prepare the string for plot command
369%initiate  the plot command
370plotstr='hhh=plot(';
371% textmean={};
372coord_x_index=[];
373xtitle='';
374ytitle='';
375test_newplot=1;
376
377%loop on input  fields
378for icell=1:length(CellVarIndex)
379    VarIndex=CellVarIndex{icell};%  indices of the selected variables in the list data.ListVarName
380    if ~isempty(VarType{icell}.coord_x)
381        coord_x_index=VarType{icell}.coord_x;
382    else
383        coord_x_index_cell=VarType{icell}.coord(1);
384        if isequal(coord_x_index_cell,0)
385             continue  % the cell has no abscissa, skip it
386        end
387        coord_x_index=coord_x_index_cell;
388    end
389    testplot=ones(size(data.ListVarName));%default test for plotted variables
390    xtitle=[xtitle data.ListVarName{coord_x_index}];
391    eval(['coord_x{icell}=data.' data.ListVarName{coord_x_index} ';']);%coordinate variable set as coord_x
392    if isfield(data,'VarAttribute')&& numel(data.VarAttribute)>=coord_x_index && isfield(data.VarAttribute{coord_x_index},'units')
393        xtitle=[xtitle '(' data.VarAttribute{coord_x_index}.units '), '];
394    else
395        xtitle=[xtitle ', '];
396    end
397    eval(['coord_x{icell}=data.' data.ListVarName{coord_x_index} ';']);%coordinate variable set as coord_x
398    XMin(icell)=min(coord_x{icell});
399    XMax(icell)=max(coord_x{icell});
400    testplot(coord_x_index)=0;
401    if ~isempty(VarType{icell}.ancillary')
402        testplot(VarType{icell}.ancillary)=0;
403    end
404    if ~isempty(VarType{icell}.warnflag')
405        testplot(VarType{icell}.warnflag)=0;
406    end
407    if isfield(data,'VarAttribute')
408        VarAttribute=data.VarAttribute;
409        for ivar=1:length(VarIndex)
410            if length(VarAttribute)>=VarIndex(ivar) && isfield(VarAttribute{VarIndex(ivar)},'long_name')
411                plotname{VarIndex(ivar)}=VarAttribute{VarIndex(ivar)}.long_name;
412            else
413                plotname{VarIndex(ivar)}=data.ListVarName{VarIndex(ivar)};%name for display in plot A METTRE
414            end
415        end
416    end
417    if ~isempty(VarType{icell}.discrete')
418        charplot_0='''+''';
419    else
420        charplot_0='''-''';
421    end
422    YMin=0;
423    YMax=1;%default
424    for ivar=1:length(VarIndex)
425        if testplot(VarIndex(ivar))
426            VarName=data.ListVarName{VarIndex(ivar)};
427            ytitle=[ytitle VarName];
428            if isfield(data,'VarAttribute')&& numel(data.VarAttribute)>=VarIndex(ivar) && isfield(data.VarAttribute{VarIndex(ivar)},'units')
429                ytitle=[ytitle '(' data.VarAttribute{VarIndex(ivar)}.units '), '];
430            else
431                ytitle=[ytitle ', '];
432            end
433            eval(['data.' VarName '=squeeze(data.' VarName ');'])
434            %eval(['min(data.' VarName ')'])
435            eval(['YMin(ivar)=min(min(data.' VarName '));'])
436            eval(['YMax(ivar)=max(max(data.' VarName '));'])
437            plotstr=[plotstr 'coord_x{' num2str(icell) '},data.' VarName ',' charplot_0 ','];
438            eval(['nbcomponent2=size(data.' VarName ',2);']);
439            eval(['nbcomponent1=size(data.' VarName ',1);']);
440            if numel(coord_x{icell})==2
441                coord_x{icell}=linspace(coord_x{icell}(1),coord_x{icell}(2),nbcomponent1);
442            end
443            if nbcomponent1==1|| nbcomponent2==1
444                legend_str=[legend_str {VarName}]; %variable with one component
445            else  %variable with severals  components
446                for ic=1:min(nbcomponent1,nbcomponent2)
447                    legend_str=[legend_str [VarName '_' num2str(ic)]]; %variable with severals  components
448                end                                                   % labeled by their index (e.g. color component)
449            end
450        end
451    end
452    YMin_cell(icell)=min(YMin);
453    YMax_cell(icell)=max(YMax);
454end
455
456%% activate the plot
457if test_newplot && ~isequal(plotstr,'hhh=plot(') 
458    set(hfig,'CurrentAxes',haxes)
459    tag=get(haxes,'tag');   
460    %%%
461    plotstr=[plotstr '''tag'',''plot_line'');'];   
462    eval(plotstr)                  %execute plot (instruction  plotstr)
463    %%%
464    set(haxes,'tag',tag)
465    grid(haxes, 'on')
466    hxlabel=xlabel(xtitle(1:end-2));% xlabel (removes ', ' at the end)
467    set(hxlabel,'Interpreter','none')% desable tex interpreter
468    if length(legend_str)>=1
469        hylabel=ylabel(ytitle(1:end-2));% ylabel (removes ', ' at the end)
470        set(hylabel,'Interpreter','none')% desable tex interpreter
471    end
472    if ~isempty(legend_str)
473        hlegend=findobj(hfig,'Tag','legend');
474        if isempty(hlegend)
475            hlegend=legend(legend_str);
476            txt=ver('MATLAB');
477            Release=txt.Release;
478            relnumb=str2double(Release(3:4));% should be changed to Version for better compatibility
479            if relnumb >= 14
480                set(hlegend,'Interpreter','none')% desable tex interpreter
481            end
482        else
483            legend_old=get(hlegend,'String');
484            if isequal(size(legend_old,1),size(legend_str,1))&&~isequal(legend_old,legend_str)
485                set(hlegend,'String',[legend_old legend_str]);
486            end
487        end
488    end
489    title_str='';
490    if isfield(data,'filename')
491       [Path, title_str, ext]=fileparts(data.filename);
492       title_str=[title_str ext];
493    end
494    if isfield(data,'Action')
495        if ~isequal(title_str,'')
496            title_str=[title_str ', '];
497        end
498        title_str=[title_str data.Action];
499    end
500    htitle=title(title_str);
501    txt=ver('MATLAB');
502    Release=txt.Release;
503    relnumb=str2double(Release(3:4));
504    if relnumb >= 14
505        set(htitle,'Interpreter','none')% desable tex interpreter
506    end
507end
508
509%% determine axes bounds
510 fix_lim=isfield(PlotParam,'FixLimits') && PlotParam.FixLimits;
511if fix_lim
512    if ~isfield(PlotParam,'MinX')||~isfield(PlotParam,'MaxX')||~isfield(PlotParam,'MinY')||~isfield(PlotParam,'MaxY')
513        fix_lim=0; %free limits if lits are not set,
514    end  %else PlotParamOut.XMin =PlotParam.XMin...
515end
516if ~fix_lim
517    PlotParamOut.MinX=min(XMin);
518    PlotParamOut.MaxX=max(XMax);
519    PlotParamOut.MinY=min(YMin_cell);
520    PlotParamOut.MaxY=max(YMax_cell);
521end
522
523%-------------------------------------------------------------------
524function [haxes,PlotParamOut,PlotType,errormsg]=plot_plane(Data,CellVarIndex,VarTypeCell,haxes,PlotParam,PosColorbar)
525%-------------------------------------------------------------------
526grid(haxes, 'off')
527%default plotting parameters
528PlotType='plane';%default
529if ~exist('PlotParam','var')
530    PlotParam=[];
531end
532if ~isfield(PlotParam,'Scalar')
533    PlotParam.Scalar=[];
534end
535if ~isfield(PlotParam,'Vectors')
536    PlotParam.Vectors=[];
537end
538
539PlotParamOut=PlotParam;%default
540hfig=get(haxes,'parent');
541hcol=findobj(hfig,'Tag','Colorbar'); %look for colorbar axes
542hima=findobj(haxes,'Tag','ima');% search existing image in the current axes
543errormsg='';%default
544test_ima=0; %default: test for image or map plot
545test_vec=0; %default: test for vector plots
546test_black=0;
547test_false=0;
548test_C=0;
549XName='';
550x_units='';
551YName='';
552y_units='';
553for icell=1:length(CellVarIndex) % length(CellVarIndex) =1 or 2 (from the calling function)
554    VarType=VarTypeCell{icell};
555    ivar_X=VarType.coord_x; % defines (unique) index for the variable representing unstructured x coordinate (default =[])
556    ivar_Y=VarType.coord_y; % defines (unique)index for the variable representing unstructured y coordinate (default =[])
557    ivar_U=VarType.vector_x; % defines (unique) index for the variable representing x vector component (default =[])
558    ivar_V=VarType.vector_y; % defines (unique) index for the variable representing y vector component (default =[])
559    ivar_C=[VarType.scalar VarType.image VarType.color VarType.ancillary]; %defines index (indices) for the scalar or ancillary fields
560    if numel(ivar_C)>1
561        errormsg= 'error in plot_field: too many scalar inputs';
562        return
563    end
564    ivar_F=VarType.warnflag; %defines index (unique) for warning flag variable
565    ivar_FF=VarType.errorflag; %defines index (unique) for error flag variable
566    ind_coord=find(VarType.coord);
567    if numel(ind_coord)==2
568        VarType.coord=VarType.coord(ind_coord);
569    end
570    if ~isempty(ivar_U) && ~isempty(ivar_V)% vector components detected
571        if test_vec
572            errormsg='error in plot_field: attempt to plot two vector fields';
573            return
574        else
575            test_vec=1;
576            eval(['vec_U=Data.' Data.ListVarName{ivar_U} ';'])
577            eval(['vec_V=Data.' Data.ListVarName{ivar_V} ';'])
578            if ~isempty(ivar_X) && ~isempty(ivar_Y)% 2D field (with unstructured coordinates or structured ones (then ivar_X and ivar_Y empty)
579                XName=Data.ListVarName{ivar_X};
580                YName=Data.ListVarName{ivar_Y};
581                eval(['vec_X=Data.' XName ';'])
582                eval(['vec_Y=Data.' YName ';'])
583            elseif numel(VarType.coord)==2 && ~isequal(VarType.coord,[0 0]);%coordinates defines by dimension variables
584                eval(['y=Data.' Data.ListVarName{VarType.coord(1)} ';'])
585                eval(['x=Data.' Data.ListVarName{VarType.coord(2)} ';'])
586                if numel(y)==2 % y defined by first and last values on aregular mesh
587                    y=linspace(y(1),y(2),size(vec_U,1));
588                end
589                if numel(x)==2 % y defined by first and last values on aregular mesh
590                    x=linspace(x(1),x(2),size(vec_U,2));
591                end
592                [vec_X,vec_Y]=meshgrid(x,y); 
593            else
594                errormsg='error in plot_field: invalid coordinate definition for vector field';
595                return
596            end
597            if ~isempty(ivar_C)
598                 eval(['vec_C=Data.' Data.ListVarName{ivar_C} ';']) ;
599                 vec_C=reshape(vec_C,1,numel(vec_C));
600                 test_C=1;
601            end
602            if ~isempty(ivar_F)%~(isfield(PlotParam.Vectors,'HideWarning')&& isequal(PlotParam.Vectors.HideWarning,1))
603                if test_vec
604                    eval(['vec_F=Data.' Data.ListVarName{ivar_F} ';']) % warning flags for  dubious vectors
605                    if  ~(isfield(PlotParam.Vectors,'HideWarning') && isequal(PlotParam.Vectors.HideWarning,1))
606                        test_black=1;
607                    end
608                end
609            end
610            if ~isempty(ivar_FF) %&& ~test_false
611                if test_vec% TODO: deal with FF for structured coordinates
612                    eval(['vec_FF=Data.' Data.ListVarName{ivar_FF} ';']) % flags for false vectors
613                end
614            end
615        end
616    elseif ~isempty(ivar_C) %scalar or image
617        if test_ima
618             errormsg='attempt to plot two scalar fields or images';
619            return
620        end
621        eval(['A=squeeze(Data.' Data.ListVarName{ivar_C} ');']) ;% scalar represented as color image
622        test_ima=1;
623        if ~isempty(ivar_X) && ~isempty(ivar_Y)% 2D field (with unstructured coordinates or structured ones (then ivar_X and ivar_Y empty)
624            XName=Data.ListVarName{ivar_X};
625            YName=Data.ListVarName{ivar_Y};
626            eval(['AX=Data.' XName ';'])
627            eval(['AY=Data.' YName ';'])
628            [A,AX,AY]=proj_grid(AX',AY',A',[],[],'np>256');  % interpolate on a grid 
629            if isfield(Data,'VarAttribute')
630                if numel(Data.VarAttribute)>=ivar_X && isfield(Data.VarAttribute{ivar_X},'units')
631                    x_units=[' (' Data.VarAttribute{ivar_X}.units ')'];
632                end
633                if numel(Data.VarAttribute)>=ivar_Y && isfield(Data.VarAttribute{ivar_Y},'units')
634                    y_units=[' (' Data.VarAttribute{ivar_Y}.units ')'];
635                end
636            end       
637        elseif numel(VarType.coord)==2 %structured coordinates
638            XName=Data.ListVarName{VarType.coord(2)};
639            YName=Data.ListVarName{VarType.coord(1)};
640            eval(['AY=Data.' Data.ListVarName{VarType.coord(1)} ';'])
641            eval(['AX=Data.' Data.ListVarName{VarType.coord(2)} ';'])
642            test_interp_X=0; %default, regularly meshed X coordinate
643            test_interp_Y=0; %default, regularly meshed Y coordinate
644            if isfield(Data,'VarAttribute')
645                if numel(Data.VarAttribute)>=VarType.coord(2) && isfield(Data.VarAttribute{VarType.coord(2)},'units')
646                    x_units=Data.VarAttribute{VarType.coord(2)}.units;
647                end
648                if numel(Data.VarAttribute)>=VarType.coord(1) && isfield(Data.VarAttribute{VarType.coord(1)},'units')
649                    y_units=Data.VarAttribute{VarType.coord(1)}.units;
650                end
651            end 
652            if numel(AY)>2
653                DAY=diff(AY);
654                DAY_min=min(DAY);
655                DAY_max=max(DAY);
656                if sign(DAY_min)~=sign(DAY_max);% =1 for increasing values, 0 otherwise
657                     errormsg=['errror in plot_field.m: non monotonic dimension variable ' Data.ListVarName{VarType.coord(1)} ];
658                      return
659                end
660                test_interp_Y=(DAY_max-DAY_min)> 0.0001*abs(DAY_max);
661            end
662            if numel(AX)>2
663                DAX=diff(AX);
664                DAX_min=min(DAX);
665                DAX_max=max(DAX);
666                if sign(DAX_min)~=sign(DAX_max);% =1 for increasing values, 0 otherwise
667                     errormsg=['errror in plot_field.m: non monotonic dimension variable ' Data.ListVarName{VarType.coord(2)} ];
668                      return
669                end
670                test_interp_X=(DAX_max-DAX_min)> 0.0001*abs(DAX_max);
671            end 
672            if test_interp_Y         
673                npxy(1)=max([256 floor((AY(end)-AY(1))/DAY_min) floor((AY(end)-AY(1))/DAY_max)]);
674                yI=linspace(AY(1),AY(end),npxy(1));
675                if ~test_interp_X
676                    xI=linspace(AX(1),AX(end),size(A,2));%default
677                    AX=xI;
678                end
679            end
680            if test_interp_X 
681                npxy(2)=max([256 floor((AX(end)-AX(1))/DAX_min) floor((AX(end)-AX(1))/DAX_max)]);
682                xI=linspace(AX(1),AX(end),npxy(2));   
683                if ~test_interp_Y
684                   yI=linspace(AY(1),AY(end),size(A,1));
685                   AY=yI;
686                end
687            end
688            if test_interp_X || test_interp_Y               
689                [AX2D,AY2D]=meshgrid(AX,AY);
690                A=interp2(AX2D,AY2D,double(A),xI,yI');
691            end
692            AX=[AX(1) AX(end)];% keep only the lower and upper bounds for image represnetation
693            AY=[AY(1) AY(end)];
694        else
695            errormsg='error in plot_field: invalid coordinate definition ';
696            return
697        end
698    end
699    %define coordinates as CoordUnits, if not defined as attribute for each variable
700    if isfield(Data,'CoordUnit')
701        if isempty(x_units)
702            x_units=Data.CoordUnit;
703        end
704        if isempty(y_units)
705            y_units=Data.CoordUnit;
706        end
707    end
708       
709end
710
711%%   image or scalar plot %%%%%%%%%%%%%%%%%%%%%%%%%%
712
713if ~isfield(PlotParam.Scalar,'Contours')
714    PlotParam.Scalar.Contours=0; %default
715end
716PlotParamOut=PlotParam; %default
717if test_ima
718    % distinguish B/W and color images
719    np=size(A);%size of image
720    siz=numel(np);
721    if siz>3
722       errormsg=['unrecognized scalar type: ' num2str(siz) ' dimensions'];
723            return
724    end
725    if siz==3
726        if np(3)==1
727            siz=2;%B W image
728        elseif np(3)==3
729            siz=3;%color image
730        else
731            errormsg=['unrecognized scalar type in plot_field: considered as 2D field with ' num2str(np(3)) ' color components'];
732            return
733        end
734    end
735   
736    %set the color map
737    if isfield(PlotParam.Scalar,'BW')
738        BW=PlotParam.Scalar.BW; %test for BW gray scale images
739    else
740        BW=(siz==2) && (isa(A,'uint8')|| isa(A,'uint16'));% non color images represented in gray scale by default
741    end
742   
743    %case of grey level images or contour plot
744    if siz==2
745        if ~isfield(PlotParam.Scalar,'FixScal')
746            PlotParam.Scalar.FixScal=0;%default
747        end
748        if ~isfield(PlotParam.Scalar,'MinA')
749            PlotParam.Scalar.MinA=[];%default
750        end
751        if ~isfield(PlotParam.Scalar,'MaxA')
752            PlotParam.Scalar.MaxA=[];%default
753        end
754        Aline=[];
755        if isequal(PlotParam.Scalar.FixScal,0)||isempty(PlotParam.Scalar.MinA)||~isa(PlotParam.Scalar.MinA,'double')  %correct if there is no numerical data in edit box
756            Aline=reshape(A,1,[]);
757            Aline=Aline(~isnan(A));
758            if isempty(Aline)
759                 errormsg='NaN input scalar or image in plot_field';
760                return
761            end
762            MinA=double(min(Aline));
763            %MinA=double(min(min(A)));
764        else
765            MinA=PlotParam.Scalar.MinA;
766        end;
767        if isequal(PlotParam.Scalar.FixScal,0)||isempty(PlotParam.Scalar.MaxA)||~isa(PlotParam.Scalar.MaxA,'double') %correct if there is no numerical data in edit box
768            if isempty(Aline)
769               Aline=reshape(A,1,[]);
770               Aline=Aline(~isnan(A));
771               if isempty(Aline)
772                 errormsg='NaN input scalar or image in plot_field';
773                return
774               end
775            end
776            MaxA=double(max(Aline));
777           % MaxA=double(max(max(A)));
778        else
779            MaxA=PlotParam.Scalar.MaxA; 
780        end;
781        PlotParamOut.Scalar.MinA=MinA;
782        PlotParamOut.Scalar.MaxA=MaxA;
783       
784        % case of contour plot
785        if isequal(PlotParam.Scalar.Contours,1)
786            if ~isempty(hima) && ishandle(hima)
787                delete(hima)
788            end
789            if ~isfield(PlotParam.Scalar,'IncrA')
790                PlotParam.Scalar.IncrA=NaN;
791            end
792            if isnan(PlotParam.Scalar.IncrA)% | PlotParam.Scalar.AutoScal==0
793                cont=colbartick(MinA,MaxA);
794                intercont=cont(2)-cont(1);%default
795                PlotParamOut.Scalar.IncrA=intercont;
796            else
797               intercont=PlotParam.Scalar.IncrA;
798            end
799            B=A;           
800            abscontmin=intercont*floor(MinA/intercont);
801            abscontmax=intercont*ceil(MaxA/intercont);
802            contmin=intercont*floor(min(min(B))/intercont);
803            contmax=intercont*ceil(max(max(B))/intercont);
804            cont_pos_plus=0:intercont:contmax;
805            cont_pos_min=double(contmin):intercont:-intercont;
806            cont_pos=[cont_pos_min cont_pos_plus];
807            sizpx=(AX(end)-AX(1))/(np(2)-1);
808            sizpy=(AY(1)-AY(end))/(np(1)-1);
809            x_cont=AX(1):sizpx:AX(end); % pixel x coordinates for image display
810            y_cont=AY(1):-sizpy:AY(end); % pixel x coordinates for image display
811           % axes(haxes)% set the input axes handle as current axis
812    txt=ver('MATLAB');
813    Release=txt.Release;
814            relnumb=str2double(Release(3:4));
815            if relnumb >= 14
816                    vec=linspace(0,1,(abscontmax-abscontmin)/intercont);%define a greyscale colormap with steps intercont
817                map=[vec' vec' vec'];
818                colormap(map);
819                [var,hcontour]=contour(x_cont,y_cont,B,cont_pos);       
820                set(hcontour,'Fill','on')
821                set(hcontour,'LineStyle','none')
822                hold on
823            end
824            [var_p,hcontour_p]=contour(x_cont,y_cont,B,cont_pos_plus,'k-');
825            hold on
826            [var_m,hcontour_m]=contour(x_cont,y_cont,B,cont_pos_min,':');
827            set(hcontour_m,'LineColor',[1 1 1])
828            hold off
829            caxis([abscontmin abscontmax])
830            colormap(map);
831        end
832       
833        % set  colormap for  image display
834        if ~isequal(PlotParam.Scalar.Contours,1) 
835            % rescale the grey levels with min and max, put a grey scale colorbar
836            B=A;
837            if BW
838                vec=linspace(0,1,255);%define a linear greyscale colormap
839                map=[vec' vec' vec'];
840                colormap(map);  %grey scale color map
841            else
842                colormap('default'); % standard faulse colors for div, vort , scalar fields
843            end
844        end
845       
846    % case of color images
847    else
848        if BW
849            B=uint16(sum(A,3));
850        else
851            B=uint8(A);
852        end
853        MinA=0;
854        MaxA=255;
855    end
856   
857    % display usual image
858    if ~isequal(PlotParam.Scalar.Contours,1)     
859        % interpolate field to increase resolution of image display
860        test_interp=1;
861        if max(np) <= 64
862            npxy=8*np;% increase the resolution 8 times
863        elseif max(np) <= 128
864            npxy=4*np;% increase the resolution 4 times
865        elseif max(np) <= 256
866            npxy=2*np;% increase the resolution 2 times
867        else
868            npxy=np;
869            test_interp=0; % no interpolation done
870        end
871        if test_interp==1%if we interpolate   
872            x=linspace(AX(1),AX(2),np(2));
873            y=linspace(AY(1),AY(2),np(1));
874            [X,Y]=meshgrid(x,y);
875            xi=linspace(AX(1),AX(2),npxy(2));
876            yi=linspace(AY(1),AY(2),npxy(1));
877            B = interp2(X,Y,double(B),xi,yi');
878        end           
879        % create new image if there  no image handle is found
880        if isempty(hima)
881            tag=get(haxes,'Tag');
882            if MinA<MaxA
883                hima=imagesc(AX,AY,B,[MinA MaxA]);
884            else % to deal with uniform field
885                hima=imagesc(AX,AY,B,[MaxA-1 MaxA]);
886            end
887            set(hima,'Tag','ima','HitTest','off')
888            set(haxes,'Tag',tag);%preserve the axes tag (removed by image fct !!!)     
889            uistack(hima, 'bottom')
890        % update an existing image
891        else
892            set(hima,'CData',B);
893            if MinA<MaxA
894                set(haxes,'CLim',[MinA MaxA])
895                %caxis([MinA MaxA])
896            else
897                set(haxes,'CLim',[MinA MaxA])
898                %caxis([MaxA-1 MaxA])
899            end
900            set(hima,'XData',AX);
901            set(hima,'YData',AY);
902        end
903    end
904%     if ~isstruct(AxeData)
905%         AxeData=[];
906%     end
907    test_ima=1;
908   
909    %display the colorbar code for B/W images if Poscolorbar not empty
910    if siz==2 && exist('PosColorbar','var')&& ~isempty(PosColorbar)
911        if isempty(hcol)||~ishandle(hcol)
912             hcol=colorbar;%create new colorbar
913        end
914        if length(PosColorbar)==4
915                 set(hcol,'Position',PosColorbar)           
916        end
917        %YTick=0;%default
918        if MaxA>MinA
919            if isequal(PlotParam.Scalar.Contours,1)
920                colbarlim=get(hcol,'YLim');
921                scale_bar=(colbarlim(2)-colbarlim(1))/(abscontmax-abscontmin);               
922                YTick=cont_pos(2:end-1);
923                YTick_scaled=colbarlim(1)+scale_bar*(YTick-abscontmin);
924                set(hcol,'YTick',YTick_scaled);
925            elseif (isfield(PlotParam.Scalar,'BW') && isequal(PlotParam.Scalar.BW,1))||isa(A,'uint8')|| isa(A,'uint16')%images
926                hi=get(hcol,'children');
927                if iscell(hi)%multiple images in colorbar
928                    hi=hi{1};
929                end
930                set(hi,'YData',[MinA MaxA])
931                set(hi,'CData',(1:256)')
932                set(hcol,'YLim',[MinA MaxA])
933                YTick=colbartick(MinA,MaxA);
934                set(hcol,'YTick',YTick)               
935            else
936                hi=get(hcol,'children');
937                if iscell(hi)%multiple images in colorbar
938                    hi=hi{1};
939                end
940                set(hi,'YData',[MinA MaxA])
941                set(hi,'CData',(1:64)')
942                YTick=colbartick(MinA,MaxA);
943                set(hcol,'YLim',[MinA MaxA])
944                set(hcol,'YTick',YTick)
945            end
946            set(hcol,'Yticklabel',num2str(YTick'));
947        end
948    elseif ishandle(hcol)
949        delete(hcol); %erase existing colorbar if not needed
950    end
951else%no scalar plot
952    if ~isempty(hima) && ishandle(hima)
953        delete(hima)
954    end
955    if ~isempty(hcol)&& ishandle(hcol)
956       delete(hcol)
957    end
958    PlotParamOut=rmfield(PlotParamOut,'Scalar');
959end
960
961%%   vector plot %%%%%%%%%%%%%%%%%%%%%%%%%%
962if test_vec
963   %vector scale representation
964    if size(vec_U,1)==numel(vec_Y) && size(vec_U,2)==numel(vec_X); % x, y  coordinate variables
965        [vec_X,vec_Y]=meshgrid(vec_X,vec_Y);
966    end   
967    vec_X=reshape(vec_X,1,numel(vec_X));%reshape in matlab vectors
968    vec_Y=reshape(vec_Y,1,numel(vec_Y));
969    vec_U=reshape(vec_U,1,numel(vec_U));
970    vec_V=reshape(vec_V,1,numel(vec_V));
971     MinMaxX=max(vec_X)-min(vec_X);
972    if  isfield(PlotParam.Vectors,'FixVec') && isequal(PlotParam.Vectors.FixVec,1)&& isfield(PlotParam.Vectors,'VecScale')...
973               &&~isempty(PlotParam.Vectors.VecScale) && isa(PlotParam.Vectors.VecScale,'double') %fixed vector scale
974        scale=PlotParam.Vectors.VecScale;  %impose the length of vector representation
975    else
976        if ~test_false %remove false vectors   
977            indsel=1:numel(vec_X);%
978        end
979        if isempty(vec_U)
980            scale=1;
981        else
982            if isempty(indsel)
983                MaxU=max(abs(vec_U));
984                MaxV=max(abs(vec_V));
985            else
986                MaxU=max(abs(vec_U(indsel)));
987                MaxV=max(abs(vec_V(indsel)));
988            end
989            scale=MinMaxX/(max(MaxU,MaxV)*50);
990            PlotParam.Vectors.VecScale=scale;%update the 'scale' display
991        end
992    end
993   
994    %record vectors on the plotting axes
995    if test_C==0
996        vec_C=ones(1,numel(vec_X));
997    end
998   
999    %decimate by a factor 2 in vector mesh(4 in nbre of vectors)
1000    if isfield(PlotParam.Vectors,'decimate4') && isequal(PlotParam.Vectors.decimate4,1)
1001        diffy=diff(vec_Y); %difference dy=vec_Y(i+1)-vec_Y(i)
1002        dy_thresh=max(abs(diffy))/2;
1003        ind_jump=find(abs(diffy) > dy_thresh); %indices with diff(vec_Y)> max/2, detect change of line
1004        ind_sel=1:ind_jump(1);%select the first line
1005        for i=2:2:length(ind_jump)-1
1006            ind_sel=[ind_sel (ind_jump(i)+1:ind_jump(i+1))];% select the odd lines
1007        end
1008        nb_sel=length(ind_sel);
1009        ind_sel=ind_sel(1:2:nb_sel);% take half the points on a line
1010        vec_X=vec_X(ind_sel);
1011        vec_Y=vec_Y(ind_sel);
1012        vec_U=vec_U(ind_sel);
1013        vec_V=vec_V(ind_sel);
1014        vec_C=vec_C(ind_sel);
1015        if ~isempty(ivar_F)
1016           vec_F=vec_F(ind_sel);
1017        end
1018        if ~isempty(ivar_FF)
1019           vec_FF=vec_FF(ind_sel);
1020        end
1021    end
1022   
1023    %get main level color code
1024    [colorlist,col_vec,PlotParamOut.Vectors]=set_col_vec(PlotParam.Vectors,vec_C);
1025   
1026    % take flags into account: add flag colors to the list of colors
1027    sizlist=size(colorlist);
1028    nbcolor=sizlist(1);
1029    if test_black
1030       nbcolor=nbcolor+1;
1031       colorlist(nbcolor,:)=[0 0 0]; %add black to the list of colors
1032       if ~isempty(ivar_FF)
1033          %  ind_flag=find(vec_F~=1 & vec_F~=0 & vec_FF==0);  %flag warning but not false
1034            col_vec(vec_F~=1 & vec_F~=0 & vec_FF==0)=nbcolor;
1035       else
1036           % ind_flag=find(vec_F~=1 & vec_F~=0);
1037            col_vec(vec_F~=1 & vec_F~=0)=nbcolor;
1038       end
1039      % col_vec(ind_flag)=nbcolor;   
1040    end
1041    nbcolor=nbcolor+1;
1042    if ~isempty(ivar_FF)
1043        %ind_flag=find(vec_FF~=0);
1044        if isfield(PlotParam.Vectors,'HideFalse') && PlotParam.Vectors.HideFalse==1
1045            colorlist(nbcolor,:)=[NaN NaN NaN];% no plot of false vectors
1046        else
1047            colorlist(nbcolor,:)=[1 0 1];% magenta color
1048        end
1049        col_vec(vec_FF~=0)=nbcolor;
1050    end
1051    %plot vectors:
1052    quiresetn(haxes,vec_X,vec_Y,vec_U,vec_V,scale,colorlist,col_vec);   
1053
1054else
1055    hvec=findobj(haxes,'Tag','vel');
1056    if ~isempty(hvec)
1057        delete(hvec);
1058    end
1059    PlotParamOut=rmfield(PlotParamOut,'Vectors');
1060end
1061
1062%listfields={'AY','AX','A','X','Y','U','V','C','W','F','FF'};
1063%listdim={'AY','AX',{'AY','AX'},'nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors'};
1064%Role={'coord_y','coord_x','scalar','coord_x','coord_y','vector_x','vector_y','scalar','vector_z','warnflag','errorflag'};
1065%ind_select=[];
1066nbvar=0;
1067
1068%store the coordinate extrema occupied by the field
1069if ~isempty(Data)
1070    fix_lim=isfield(PlotParam,'FixLimits') && PlotParam.FixLimits;
1071    if fix_lim
1072        if ~isfield(PlotParam,'MinX')||~isfield(PlotParam,'MaxX')||~isfield(PlotParam,'MinY')||~isfield(PlotParam,'MaxY')
1073            fix_lim=0; %free limits if lits are not set,
1074        end  %else PlotParamOut.XMin =PlotParam.XMin...
1075    end
1076    if ~fix_lim
1077        XMin=[];
1078        XMax=[];
1079        YMin=[];
1080        YMax=[];
1081        if test_ima%both background image and vectors coexist, take the wider bound
1082            XMin=min(AX);
1083            XMax=max(AX);
1084            YMin=min(AY);
1085            YMax=max(AY);
1086        end
1087        if test_vec
1088            XMin=[XMin min(vec_X)];
1089            XMax=[XMax max(vec_X)];
1090            YMin=[YMin min(vec_Y)];
1091            YMax=[YMax max(vec_Y)];
1092        end
1093        PlotParamOut.MinX=min(XMin);
1094        PlotParamOut.MaxX=max(XMax);
1095        PlotParamOut.MinY=min(YMin);
1096        PlotParamOut.MaxY=max(YMax);
1097        if XMax>XMin
1098            set(haxes,'XLim',[XMin XMax]);% set x limits of frame in axes coordinates
1099        end
1100        if YMax>YMin
1101            set(haxes,'YLim',[YMin YMax]);% set x limits of frame in axes coordinates
1102        end
1103    end
1104%     if Ylim(2)>Ylim(1)
1105%         set(haxes,'YLim',Ylim);% set y limits of frame in axes coordinate
1106%     end
1107%     end
1108    %    adjust the size of the plot to include the whole field, except if PlotParam.FixLimits=1
1109    %     if ~(isfield(PlotParam,'FixLimits') && PlotParam.FixLimits) && test_lim
1110    %         PlotParamOut.MinX=Xlim(1);
1111    %         PlotParamOut.MaxX=Xlim(2);
1112    %         PlotParamOut.MinY=Ylim(1);
1113    %         PlotParamOut.MaxY=Ylim(2);
1114%     if Xlim(2)>Xlim(1)
1115%         set(haxes,'XLim',Xlim);% set x limits of frame in axes coordinates
1116%     end
1117%     if Ylim(2)>Ylim(1)
1118%         set(haxes,'YLim',Ylim);% set y limits of frame in axes coordinate
1119%     end
1120    %     end
1121
1122    set(haxes,'YDir','normal')
1123    set(get(haxes,'XLabel'),'String',[XName ' (' x_units ')']);
1124    set(get(haxes,'YLabel'),'String',[YName ' (' y_units ')']);
1125    PlotParamOut.x_units=x_units;
1126    PlotParamOut.y_units=y_units;
1127end
1128%-------------------------------------------------------------------
1129% --- function for plotting vectors
1130%INPUT:
1131% haxes: handles of the plotting axes
1132% x,y,u,v: vectors coordinates and vector components to plot, arrays withb the same dimension
1133% scale: scaling factor for vector length representation
1134% colorlist(icolor,:): list of vector colors, dim (nbcolor,3), depending on color #i
1135% col_vec: matlab vector setting the color number #i for each velocity vector
1136function quiresetn(haxes,x,y,u,v,scale,colorlist,col_vec)
1137%-------------------------------------------------------------------
1138%define arrows
1139theta=0.5 ;%angle arrow
1140alpha=0.3 ;%length arrow
1141rot=alpha*[cos(theta) -sin(theta); sin(theta) cos(theta)]';
1142%find the existing lines
1143h=findobj(haxes,'Tag','vel');% search existing lines in the current axes
1144sizh=size(h);
1145set(h,'EraseMode','xor');
1146set(haxes,'NextPlot','replacechildren');
1147
1148%drawnow
1149%create lines (if no lines) or modify them
1150if ~isequal(size(col_vec),size(x))
1151    col_vec=ones(size(x));% case of error in col_vec input
1152end
1153sizlist=size(colorlist);
1154ncolor=sizlist(1);
1155
1156for icolor=1:ncolor
1157    %determine the line positions for each color icolor
1158    ind=find(col_vec==icolor);
1159    xc=x(ind);
1160    yc=y(ind);
1161    uc=u(ind)*scale;
1162    vc=v(ind)*scale;
1163    n=size(xc);
1164    xN=NaN*ones(size(xc));
1165    matx=[xc(:)-uc(:)/2 xc(:)+uc(:)/2 xN(:)]';
1166    %     matx=[xc(:) xc(:)+uc(:) xN(:)]';
1167    matx=reshape(matx,1,3*n(2));
1168    maty=[yc(:)-vc(:)/2 yc(:)+vc(:)/2 xN(:)]';
1169    %     maty=[yc(:) yc(:)+vc(:) xN(:)]';
1170    maty=reshape(maty,1,3*n(2));
1171   
1172    %determine arrow heads
1173    arrowplus=rot*[uc;vc];
1174    arrowmoins=rot'*[uc;vc];
1175    x1=xc+uc/2-arrowplus(1,:);
1176    x2=xc+uc/2;
1177    x3=xc+uc/2-arrowmoins(1,:);
1178    y1=yc+vc/2-arrowplus(2,:);
1179    y2=yc+vc/2;
1180    y3=yc+vc/2-arrowmoins(2,:);
1181    matxar=[x1(:) x2(:) x3(:) xN(:)]';
1182    matxar=reshape(matxar,1,4*n(2));
1183    matyar=[y1(:) y2(:) y3(:) xN(:)]';
1184    matyar=reshape(matyar,1,4*n(2));
1185    %draw the line or modify the existing ones
1186    tri=reshape(1:3*length(uc),3,[])';   
1187    isn=isnan(colorlist(icolor,:));%test if color NaN
1188    if 2*icolor > sizh(1) %if icolor exceeds the number of existing ones
1189        if ~isn(1) %if the vectors are visible color not nan
1190            if n(2)>0
1191                hold on
1192                line(matx,maty,'Color',colorlist(icolor,:),'Tag','vel');% plot new lines
1193                line(matxar,matyar,'Color',colorlist(icolor,:),'Tag','vel');% plot arrows
1194            end
1195        end
1196    else
1197        if isn(1)
1198            delete(h(2*icolor-1))
1199            delete(h(2*icolor))
1200        else
1201            set(h(2*icolor-1),'Xdata',matx,'Ydata',maty);
1202            set(h(2*icolor-1),'Color',colorlist(icolor,:));
1203            set(h(2*icolor-1),'EraseMode','xor');
1204            set(h(2*icolor),'Xdata',matxar,'Ydata',matyar);
1205            set(h(2*icolor),'Color',colorlist(icolor,:));
1206            set(h(2*icolor),'EraseMode','xor');
1207        end
1208    end
1209end
1210if sizh(1) > 2*ncolor
1211    for icolor=ncolor+1 : sizh(1)/2%delete additional objects
1212        delete(h(2*icolor-1))
1213        delete(h(2*icolor))
1214    end
1215end
1216
1217%-------------------------------------------------------------------
1218% ---- determine tick positions for colorbar
1219function YTick=colbartick(MinA,MaxA)
1220%-------------------------------------------------------------------
1221%determine tick positions with "simple" values between MinA and MaxA
1222YTick=0;%default
1223maxabs=max([abs(MinA) abs(MaxA)]);
1224if maxabs>0
1225ord=10^(floor(log10(maxabs)));%order of magnitude
1226div=1;
1227siz2=1;
1228while siz2<2
1229%     values=[-9:div:9];
1230    values=-10:div:10;
1231    ind=find((ord*values-MaxA)<0 & (ord*values-MinA)>0);%indices of 'values' such that MinA<ord*values<MaxA
1232    siz=size(ind);
1233    if siz(2)<4%if there are less than 4 selected values (4 levels)
1234        values=-9:0.5*div:9;
1235        ind=find((ord*values-MaxA)<0 & (ord*values-MinA)>0);
1236    end
1237    siz2=size(ind,2);
1238%     siz2=siz(2)
1239    div=div/10;
1240end
1241YTick=ord*values(ind);
1242end
1243
Note: See TracBrowser for help on using the repository browser.