source: trunk/src/plot_field.m @ 364

Last change on this file since 364 was 364, checked in by sommeria, 12 years ago

civ2, fix2, patch2 introduced under Matlab (no interpolation and rotation done yet)
small corrections in plot_field and read_civdata

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