source: trunk/src/plot_field.m @ 426

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

compiled matlab civ introduced in Batch mode,
minor bugs fixed for histogram of velocity

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