[201]  1  %'plot_field': plot any field with the structure defined in the uvmat package


 2  %


 3  %


 4  % This function is used by uvmat to plot fields. It automatically chooses the representation


 5  % appropriate to the input field structure:


[292]  6  % 2D vector fields are represented by arrows, 2D scalar fields by grey scale images or contour plots, 1D fields are represented by usual plot with (abscissa, ordinate).


[201]  7  % The input field structure is first tested by check_field_structure.m,


 8  % then split into blocks of related variables by find_field_indices.m.


[292]  9  % The dimensionality of each block is obtained by this function


[201]  10  % considering the presence of variables with the attribute .Role='coord_x'


 11  % and/or coord_y and/or coord_z (case of unstructured coordinates), or


 12  % dimension variables (case of matrices).


 13  %


 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  %


[294]  45  % PlotParam: structure containing the parameters for plotting, as read on the uvmat or view_field GUI (by function 'read_GUI.m').


 46  % Contains three substructures:


 47  % .Coordinates: coordinate parameters:


 48  % .CheckFixLimits:=0 (default) adjust axes limit to the X,Y data, =1: preserves the previous axes limits


[428]  49  % .Coordinates.CheckFixAspectRatio: =0 (default):automatic adjustment of the graph, keep 1 to 1 aspect ratio for x and y scales.


[429]  50  % .Coordinates.AspectRatio: imposed aspect ratio y/x of axis unit plots


[201]  51  % scalars


 52  % .Scalar.MaxA: upper bound (saturation color) for the scalar representation, max(field) by default


 53  % .Scalar.MinA: lower bound (saturation) for the scalar representation, min(field) by default


[292]  54  % .Scalar.CheckFixScal: =0 (default) lower and upper bounds of the scalar representation set to the min and max of the field


[201]  55  % =1 lower and upper bound imposed by .AMax and .MinA


[421]  56  % .Scalar.CheckBW= 1: black and white representation imposed, =0 color imposed (color scale or rgb),


 57  % =[]: automatic (B/W for integer positive scalars, color else)


[292]  58  % .Scalar.CheckContours= 1: represent scalars by contour plots (Matlab function 'contour'); =0 by default


[201]  59  % .IncrA : contour interval


 60  %  vectors


 61  % .Vectors.VecScale: scale for the vector representation


[292]  62  % .Vectors.CheckFixVec: =0 (default) automatic length for vector representation, =1: length set by .VecScale


 63  % .Vectors.CheckHideFalse= 0 (default) false vectors represented in magenta, =1: false vectors not represented;


 64  % .Vectors.CheckHideWarning= 0 (default) vectors marked by warnflag~=0 marked in black, 1: no warning representation;


 65  % .Vectors.CheckDecimate4 = 0 (default) all vectors reprtesented, =1: half of the vectors represented along each coordinate


[201]  66  %  vector color


 67  % .Vectors.ColorCode= 'black','white': imposed color (default ='blue')


 68  % 'rgb', : three colors red, blue, green depending


 69  % on thresholds .colcode1 and .colcode2 on the input scalar value (C)


 70  % 'brg': like rgb but reversed color order (blue, green, red)


 71  % '64 colors': continuous color from blue to red (multijet)


 72  % .Vectors.colcode1 : first threshold for rgb, first value for'continuous'


 73  % .Vectors.colcode2 : second threshold for rgb, last value (saturation) for 'continuous'


[292]  74  % .Vectors.CheckFixedCbounds; =0 (default): the bounds on C representation are min and max, =1: they are fixed by .Minc and .MaxC


[201]  75  % .Vectors.MinC = imposed minimum of the scalar field used for vector color;


 76  % .Vectors.MaxC = imposed maximum of the scalar field used for vector color;


 77  %


 78  % PosColorbar: if not empty, display a colorbar for B&W images


 79  % imposed position of the colorbar (ex [0.821 0.471 0.019 0.445])


 80 


 81  %AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA


 82  % Copyright Joel Sommeria, 2008, LEGI / CNRSUJFINPG, sommeria@coriolislegi.org.


 83  %AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA


 84  % This file is part of the toolbox UVMAT.


 85  %


 86  % UVMAT is free software; you can redistribute it and/or modify


 87  % it under the terms of the GNU General Public License as published by


 88  % the Free Software Foundation; either version 2 of the License, or


 89  % (at your option) any later version.


 90  %


 91  % UVMAT is distributed in the hope that it will be useful,


 92  % but WITHOUT ANY WARRANTY; without even the implied warranty of


 93  % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


 94  % GNU General Public License (file UVMAT/COPYING.txt) for more details.


 95  %AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA


 96 


 97  function [PlotType,PlotParamOut,haxes]= plot_field(Data,haxes,PlotParam,PosColorbar)


[247]  98 


[388]  99  %% default input and output


[201]  100  if ~exist('PlotParam','var'),PlotParam=[];end;


 101  if ~exist('PosColorbar','var'),PosColorbar=[];end;


 102  PlotType='text'; %default


 103  PlotParamOut=PlotParam;%default


[411]  104  if ~isfield(PlotParam,'Coordinates')


 105  PlotParam.Coordinates=[];


[292]  106  end


[201]  107 


[388]  108  %% check input structure


 109  index_2D=[];


 110  index_1D=[];


 111  index_0D=[];


[404]  112  errormsg=check_field_structure(Data);


[388]  113  if ~isempty(errormsg)


 114  msgbox_uvmat('ERROR',['input of plot_field/check_field_structure: ' errormsg])


 115  display(['input of plot_field/check_field_structure: ' errormsg])


 116  return


 117  end


 118  % check the cells of fields :


 119  [CellVarIndex,NbDim,VarType,errormsg]=find_field_indices(Data);


 120  if ~isempty(errormsg)


 121  msgbox_uvmat('ERROR',['input of plot_field/find_field_indices: ' errormsg]);


 122  return


 123  end


 124  index_2D=find(NbDim==2,2);%find 2D fields (at most 2)


 125  index_3D=find(NbDim>2,1);


 126  if ~isempty(index_3D)


 127  if isfield(Data,'NbDim')&& isequal(Data.NbDim,2)


 128  index_2D=[index_2D index_3D];


 129  else


 130  msgbox_uvmat('ERROR','volume plot not implemented yet');


 131  return


 132  end


 133  end


 134  index_1D=find(NbDim==1);


 135  index_0D=find(NbDim==0);


 136  %remove coordinates variables from 1D plot


 137  if ~isempty(index_2D)


 138  for ivar=1:length(index_1D)


 139  if isequal(CellVarIndex{index_1D(ivar)},VarType{index_1D(ivar)}.coord)


 140  index_1D(ivar)=0;


[201]  141  end


 142  end


[388]  143  index_1D=index_1D(index_1D>0);


[201]  144  end


[294]  145 


[388]  146  %% test axes and figure


[429]  147  testnewfig=1;%test to create a new figure (default)


 148  testzoomaxes=0;%test for the existence of a zoom secondary figure attached to the plotting axes


 149  if exist('haxes','var')


 150  if ishandle(haxes)


 151  if isequal(get(haxes,'Type'),'axes')


 152  testnewfig=0;


 153  AxeData=get(haxes,'UserData');


 154  if isfield(AxeData,'ZoomAxes')&& ishandle(AxeData.ZoomAxes)


 155  if isequal(get(AxeData.ZoomAxes,'Type'),'axes')


 156  testzoomaxes=1;


 157  zoomaxes=AxeData.ZoomAxes;


[388]  158  end


 159  end


 160  end


[201]  161  end


[429]  162  end


 163 


 164  %% create a new figure and axes if the plotting axes does not exist


 165  if testnewfig


 166  hfig=figure;


 167  set(hfig,'Units','normalized')


 168  haxes=axes;


 169  set(haxes,'position',[0.13,0.2,0.775,0.73])


 170  PlotParam.NextPlot='add'; %parameter for plot_profile and plot_his


 171  else


 172  hfig=get(haxes,'parent');


 173  set(0,'CurrentFigure',hfig)% the parent of haxes becomes the current figure


 174  set(hfig,'CurrentAxes',haxes)% haxes becomes the current axes of the parent figure


 175  end


 176 


 177  %% set axes properties


 178  if isfield(PlotParam.Coordinates,'CheckFixLimits') && isequal(PlotParam.Coordinates.CheckFixLimits,1) %adjust the graph limits


 179  set(haxes,'XLimMode', 'manual')


 180  set(haxes,'YLimMode', 'manual')


 181  else


 182  set(haxes,'XLimMode', 'auto')


 183  set(haxes,'YLimMode', 'auto')


 184  end


 185  % if ~isfield(PlotParam.Coordinates,'CheckFixAspectRatio')&& isfield(Data,'CoordUnit')


 186  % PlotParam.Coordinates.CheckFixAspectRatio=1;% if CoordUnit is defined, the two coordiantes should be plotted with equal scale by default


 187  % end


 188  errormsg='';


 189  PlotParamOut.Coordinates=[]; %default output


 190  AxeData=get(haxes,'UserData');


 191 


 192  %% 2D plots


 193  if isempty(index_2D)


 194  plot_plane([],[],[],haxes);%removes images or vector plots in the absence of 2D field plot


 195  else %plot 2D field


 196  [tild,PlotParamOut,PlotType,errormsg]=plot_plane(Data,CellVarIndex(index_2D),VarType(index_2D),haxes,PlotParam,PosColorbar);


 197  AxeData.NbDim=2;


 198  if testzoomaxes && isempty(errormsg)


 199  [zoomaxes,PlotParamOut,tild,errormsg]=plot_plane(Data,CellVarIndex(index_2D),VarType(index_2D),zoomaxes,PlotParam,PosColorbar);


 200  AxeData.ZoomAxes=zoomaxes;


[201]  201  end


[429]  202  end


 203 


 204  %% 1D plot (usual graph y vs x)


 205  if isempty(index_1D)


 206  if ~isempty(haxes)


 207  plot_profile([],[],[],haxes);%removes usual praphs y vs x in the absence of 1D field plot


[201]  208  end


[429]  209  else %plot 1D field (usual graph y vs x)


 210  Coordinates=plot_profile(Data,CellVarIndex(index_1D),VarType(index_1D),haxes,PlotParam.Coordinates);%


 211  if testzoomaxes


 212  [zoomaxes,Coordinates]=plot_profile(Data,CellVarIndex(index_1D),VarType(index_1D),zoomaxes,PlotParam.Coordinates);


 213  AxeData.ZoomAxes=zoomaxes;


[201]  214  end


[429]  215  if ~isempty(Coordinates)


 216  PlotParamOut.Coordinates=Coordinates;


[201]  217  end


[429]  218  PlotType='line';


 219  end


 220 


 221  %% text display


 222  if isempty(index_2D) && isempty(index_1D)%text display alone


 223  htext=findobj(hfig,'Tag','TableDisplay');


 224  else %text display added to plot


[388]  225  htext=findobj(hfig,'Tag','text_display');


[429]  226  end


 227  if ~isempty(htext)


 228  if isempty(index_0D)


[434]  229  if strcmp(get(htext,'Type'),'uitable')


 230  set(htext,'Data',{})


 231  else


 232  set(htext,'String',{''})


 233  end


[429]  234  else


 235  [errormsg]=plot_text(Data,CellVarIndex(index_0D),VarType(index_0D),htext);


[294]  236  end


[201]  237  end


 238 


[388]  239  %% display error message


[201]  240  if ~isempty(errormsg)


 241  msgbox_uvmat('ERROR', errormsg)


 242  end


 243 


 244  %% update the parameters stored in AxeData


[429]  245  if ishandle(haxes)&&( ~isempty(index_2D) ~isempty(index_1D))


 246  % AxeData=[];


[388]  247  if isfield(PlotParamOut,'MinX')


 248  AxeData.RangeX=[PlotParamOut.MinX PlotParamOut.MaxX];%'[PlotParamOut.MinX PlotParamOut.MaxX];


 249  AxeData.RangeY=[PlotParamOut.MinY PlotParamOut.MaxY];%[PlotParamOut.MinY PlotParamOut.MaxY]


 250  end


 251  set(haxes,'UserData',AxeData)


 252  end


[201]  253 


[206]  254  %% update the plotted field stored in parent figure


[429]  255  if ~isempty(index_2D) ~isempty(index_1D)


 256  FigData=get(hfig,'UserData');


 257  % if strcmp(get(hfig,'tag'),'view_field')


 258  % set(hfig,'UserData',[]); % refresh user data in view_field (set by civ/TestCiv )


 259  % end


 260  tagaxes=get(haxes,'tag');% tag of the current plot axis


 261  if isfield(FigData,tagaxes)


 262  FigData.(tagaxes)=Data;


 263  set(hfig,'UserData',FigData)


 264  end


[388]  265  end


[292]  266 


[201]  267  %


[389]  268  function errormsg=plot_text(FieldData,CellVarIndex,VarTypeCell,htext)


[201]  269  %


 270  errormsg=[];


 271  txt_cell={};


[429]  272  Data={};


[201]  273  for icell=1:length(CellVarIndex)


 274  VarIndex=CellVarIndex{icell};% indices of the selected variables in the list data.ListVarName


 275  for ivar=1:length(VarIndex)


[389]  276  checkancillary=0;


 277  if length(FieldData.VarAttribute)>=VarIndex(ivar)


 278  VarAttribute=FieldData.VarAttribute{VarIndex(ivar)};


[397]  279  if isfield(VarAttribute,'Role')&&(strcmp(VarAttribute.Role,'ancillary')strcmp(VarAttribute.Role,'coord_tps')...


 280  strcmp(VarAttribute.Role,'vector_x_tps')strcmp(VarAttribute.Role,'vector_y_tps'))


[389]  281  checkancillary=1;


 282  end


 283  end


 284  if ~checkancillary% does not display variables with attribute '.Role=ancillary'


 285  VarName=FieldData.ListVarName{VarIndex(ivar)};


 286  VarValue=FieldData.(VarName);


[429]  287  Data =[Data [{VarName}; num2cell(VarValue)]];


[389]  288  if size(VarValue,1)~=1


[429]  289  VarValue=VarValue';% put the different values on a line


[389]  290  end


[397]  291  if size(VarValue,1)==1


[389]  292  txt=[VarName '=' num2str(VarValue)];


 293  txt_cell=[txt_cell;{txt}];


[397]  294  end


[389]  295  end


[201]  296  end


 297  end


[429]  298  if strcmp(get(htext,'Type'),'uitable')


 299  get(htext,'ColumnName')


 300  set(htext,'ColumnName',Data(1,:))


 301  set(htext,'Data',Data(2:end,:))


 302  else


 303  set(htext,'String',txt_cell)


 304  set(htext,'UserData',txt_cell)% for storage during mouse display


 305  end


[201]  306 


[429]  307 


[201]  308  %


[292]  309  function CoordinatesOut=plot_profile(data,CellVarIndex,VarType,haxes,Coordinates)


[201]  310  %


 311 


[292]  312  if ~exist('Coordinates','var')


 313  Coordinates=[];


[201]  314  end


[292]  315  CoordinatesOut=Coordinates; %default


[201]  316  hfig=get(haxes,'parent');


 317  %suppress existing plot isf empty data


 318  if isempty(data)


 319  hplot=findobj(haxes,'tag','plot_line');


 320  if ~isempty(hplot)


 321  delete(hplot)


 322  end


 323  hlegend=findobj(hfig,'tag','legend');


 324  if ~isempty(hlegend)


 325  delete(hlegend)


 326  end


 327  return


 328  end


 329 


 330  ColorOrder=[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];


 331  set(haxes,'ColorOrder',ColorOrder)


[292]  332  if isfield(Coordinates,'NextPlot')


 333  set(haxes,'NextPlot',Coordinates.NextPlot)


[201]  334  end


 335  % adjust the size of the plot to include the whole field,


 336 


 337  legend_str={};


 338 


 339  %% prepare the string for plot command


 340  plotstr='hhh=plot(';


 341  coord_x_index=[];


 342  xtitle='';


 343  ytitle='';


 344  test_newplot=1;


 345 


 346  %loop on input fields


 347  for icell=1:length(CellVarIndex)


 348  VarIndex=CellVarIndex{icell};% indices of the selected variables in the list data.ListVarName


 349  if ~isempty(VarType{icell}.coord_x)


 350  coord_x_index=VarType{icell}.coord_x;


 351  else


 352  coord_x_index_cell=VarType{icell}.coord(1);


 353  if isequal(coord_x_index_cell,0)


 354  continue % the cell has no abscissa, skip it


 355  end


 356  coord_x_index=coord_x_index_cell;


 357  end


 358  testplot=ones(size(data.ListVarName));%default test for plotted variables


[426]  359  coord_x_name{icell}=data.ListVarName{coord_x_index};


 360  coord_x{icell}=data.(data.ListVarName{coord_x_index});%coordinate variable set as coord_x


 361  if isempty(find(strcmp(coord_x_name{icell},coord_x_name(1:end1)))) %xtitle not already selected


 362  xtitle=[xtitle coord_x_name{icell}];


 363  if isfield(data,'VarAttribute')&& numel(data.VarAttribute)>=coord_x_index && isfield(data.VarAttribute{coord_x_index},'units')


 364  xtitle=[xtitle '(' data.VarAttribute{coord_x_index}.units '), '];


 365  else


 366  xtitle=[xtitle ', '];


 367  end


[201]  368  end


 369  XMin(icell)=min(coord_x{icell});


 370  XMax(icell)=max(coord_x{icell});


 371  testplot(coord_x_index)=0;


 372  if ~isempty(VarType{icell}.ancillary')


 373  testplot(VarType{icell}.ancillary)=0;


 374  end


 375  if ~isempty(VarType{icell}.warnflag')


 376  testplot(VarType{icell}.warnflag)=0;


 377  end


 378  if isfield(data,'VarAttribute')


 379  VarAttribute=data.VarAttribute;


 380  for ivar=1:length(VarIndex)


 381  if length(VarAttribute)>=VarIndex(ivar) && isfield(VarAttribute{VarIndex(ivar)},'long_name')


 382  plotname{VarIndex(ivar)}=VarAttribute{VarIndex(ivar)}.long_name;


 383  else


 384  plotname{VarIndex(ivar)}=data.ListVarName{VarIndex(ivar)};%name for display in plot A METTRE


 385  end


 386  end


 387  end


 388  if ~isempty(VarType{icell}.discrete')


 389  charplot_0='''+''';


 390  else


 391  charplot_0='''''';


 392  end


 393  YMin=0;


 394  YMax=1;%default


 395  for ivar=1:length(VarIndex)


 396  if testplot(VarIndex(ivar))


 397  VarName=data.ListVarName{VarIndex(ivar)};


 398  ytitle=[ytitle VarName];


 399  if isfield(data,'VarAttribute')&& numel(data.VarAttribute)>=VarIndex(ivar) && isfield(data.VarAttribute{VarIndex(ivar)},'units')


 400  ytitle=[ytitle '(' data.VarAttribute{VarIndex(ivar)}.units '), '];


 401  else


 402  ytitle=[ytitle ', '];


 403  end


 404  eval(['data.' VarName '=squeeze(data.' VarName ');'])


 405  %eval(['min(data.' VarName ')'])


[294]  406  YMin(ivar)=min(min(data.(VarName)));


 407  YMax(ivar)=max(max(data.(VarName)));


[201]  408  plotstr=[plotstr 'coord_x{' num2str(icell) '},data.' VarName ',' charplot_0 ','];


 409  eval(['nbcomponent2=size(data.' VarName ',2);']);


 410  eval(['nbcomponent1=size(data.' VarName ',1);']);


 411  if numel(coord_x{icell})==2


 412  coord_x{icell}=linspace(coord_x{icell}(1),coord_x{icell}(2),nbcomponent1);


 413  end


 414  if nbcomponent1==1 nbcomponent2==1


 415  legend_str=[legend_str {VarName}]; %variable with one component


 416  else %variable with severals components


 417  for ic=1:min(nbcomponent1,nbcomponent2)


 418  legend_str=[legend_str [VarName '_' num2str(ic)]]; %variable with severals components


 419  end % labeled by their index (e.g. color component)


 420  end


 421  end


 422  end


 423  YMin_cell(icell)=min(YMin);


 424  YMax_cell(icell)=max(YMax);


 425  end


 426 


 427  %% activate the plot


 428  if test_newplot && ~isequal(plotstr,'hhh=plot(')


 429  set(hfig,'CurrentAxes',haxes)


 430  tag=get(haxes,'tag');


 431  %%%


 432  plotstr=[plotstr '''tag'',''plot_line'');'];


 433  eval(plotstr) %execute plot (instruction plotstr)


 434  %%%


 435  set(haxes,'tag',tag)


 436  grid(haxes, 'on')


 437  hxlabel=xlabel(xtitle(1:end2));% xlabel (removes ', ' at the end)


 438  set(hxlabel,'Interpreter','none')% desable tex interpreter


 439  if length(legend_str)>=1


 440  hylabel=ylabel(ytitle(1:end2));% ylabel (removes ', ' at the end)


 441  set(hylabel,'Interpreter','none')% desable tex interpreter


 442  end


 443  if ~isempty(legend_str)


 444  hlegend=findobj(hfig,'Tag','legend');


 445  if isempty(hlegend)


 446  hlegend=legend(legend_str);


 447  txt=ver('MATLAB');


 448  Release=txt.Release;


 449  relnumb=str2double(Release(3:4));% should be changed to Version for better compatibility


 450  if relnumb >= 14


 451  set(hlegend,'Interpreter','none')% desable tex interpreter


 452  end


 453  else


 454  legend_old=get(hlegend,'String');


 455  if isequal(size(legend_old,1),size(legend_str,1))&&~isequal(legend_old,legend_str)


 456  set(hlegend,'String',[legend_old legend_str]);


 457  end


 458  end


 459  end


 460  title_str='';


 461  if isfield(data,'filename')


 462  [Path, title_str, ext]=fileparts(data.filename);


 463  title_str=[title_str ext];


 464  end


 465  if isfield(data,'Action')


 466  if ~isequal(title_str,'')


 467  title_str=[title_str ', '];


 468  end


 469  title_str=[title_str data.Action];


 470  end


 471  htitle=title(title_str);


[428]  472  % txt=ver('MATLAB');


 473  % Release=txt.Release;


 474  % relnumb=str2double(Release(3:4));


 475  % if relnumb >= 14


 476  set(htitle,'Interpreter','none')% desable tex interpreter


 477  % end


[201]  478  end


 479 


 480  %% determine axes bounds


[294]  481  %CoordinatesOut.RangeX=[min(XMin) max(XMax)];


 482  %CoordinatesOut.RangeY=[min(YMin_cell) max(YMax_cell)];


 483  fix_lim=isfield(Coordinates,'CheckFixLimits') && Coordinates.CheckFixLimits;


[201]  484  if fix_lim


[292]  485  if ~isfield(Coordinates,'MinX')~isfield(Coordinates,'MaxX')~isfield(Coordinates,'MinY')~isfield(Coordinates,'MaxY')


[201]  486  fix_lim=0; %free limits if lits are not set,


[221]  487  end


[201]  488  end


[294]  489  if fix_lim


 490  set(haxes,'XLim',[Coordinates.MinX Coordinates.MaxX])


 491  set(haxes,'YLim',[Coordinates.MinY Coordinates.MaxY])


 492  else


[292]  493  CoordinatesOut.MinX=min(XMin);


 494  CoordinatesOut.MaxX=max(XMax);


 495  CoordinatesOut.MinY=min(YMin_cell);


 496  CoordinatesOut.MaxY=max(YMax_cell);


[201]  497  end


 498 


[428]  499  %% determine plot aspect ratio


[429]  500  if isfield(Coordinates,'CheckFixAspectRatio') && isequal(Coordinates.CheckFixAspectRatio,1)&&isfield(Coordinates,'AspectRatio')


[428]  501  set(haxes,'DataAspectRatioMode','manual')


 502  set(haxes,'DataAspectRatio',[Coordinates.AspectRatio 1 1])


 503  else


 504  set(haxes,'DataAspectRatioMode','auto')%automatic aspect ratio


[429]  505  AspectRatio=get(haxes,'DataAspectRatio');


[428]  506  CoordinatesOut.AspectRatio=AspectRatio(1)/AspectRatio(2);


 507  end


 508 


[201]  509  %


 510  function [haxes,PlotParamOut,PlotType,errormsg]=plot_plane(Data,CellVarIndex,VarTypeCell,haxes,PlotParam,PosColorbar)


 511  %


[411]  512 


 513  grid(haxes, 'off')% remove grid (possibly remaining from other graphs)


[201]  514  %default plotting parameters


 515  PlotType='plane';%default


 516  if ~exist('PlotParam','var')


 517  PlotParam=[];


 518  end


[411]  519 


[201]  520  if ~isfield(PlotParam,'Scalar')


 521  PlotParam.Scalar=[];


 522  end


 523  if ~isfield(PlotParam,'Vectors')


 524  PlotParam.Vectors=[];


 525  end


 526 


 527  PlotParamOut=PlotParam;%default


 528  hfig=get(haxes,'parent');


 529  hcol=findobj(hfig,'Tag','Colorbar'); %look for colorbar axes


 530  hima=findobj(haxes,'Tag','ima');% search existing image in the current axes


 531  errormsg='';%default


 532  test_ima=0; %default: test for image or map plot


 533  test_vec=0; %default: test for vector plots


 534  test_black=0;


 535  test_false=0;


 536  test_C=0;


 537  XName='';


 538  x_units='';


 539  YName='';


 540  y_units='';


 541  for icell=1:length(CellVarIndex) % length(CellVarIndex) =1 or 2 (from the calling function)


 542  VarType=VarTypeCell{icell};


[402]  543  if ~isempty(VarType.coord_tps)


 544  continue


 545  end


[201]  546  ivar_X=VarType.coord_x; % defines (unique) index for the variable representing unstructured x coordinate (default =[])


 547  ivar_Y=VarType.coord_y; % defines (unique)index for the variable representing unstructured y coordinate (default =[])


 548  ivar_U=VarType.vector_x; % defines (unique) index for the variable representing x vector component (default =[])


 549  ivar_V=VarType.vector_y; % defines (unique) index for the variable representing y vector component (default =[])


 550  ivar_C=[VarType.scalar VarType.image VarType.color VarType.ancillary]; %defines index (indices) for the scalar or ancillary fields


 551  if numel(ivar_C)>1


 552  errormsg= 'error in plot_field: too many scalar inputs';


 553  return


 554  end


 555  ivar_F=VarType.warnflag; %defines index (unique) for warning flag variable


 556  ivar_FF=VarType.errorflag; %defines index (unique) for error flag variable


 557  ind_coord=find(VarType.coord);


 558  if numel(ind_coord)==2


 559  VarType.coord=VarType.coord(ind_coord);


 560  end


 561  if ~isempty(ivar_U) && ~isempty(ivar_V)% vector components detected


 562  if test_vec


 563  errormsg='error in plot_field: attempt to plot two vector fields';


 564  return


 565  else


 566  test_vec=1;


[313]  567  vec_U=Data.(Data.ListVarName{ivar_U});


 568  vec_V=Data.(Data.ListVarName{ivar_V});


[201]  569  if ~isempty(ivar_X) && ~isempty(ivar_Y)% 2D field (with unstructured coordinates or structured ones (then ivar_X and ivar_Y empty)


 570  XName=Data.ListVarName{ivar_X};


 571  YName=Data.ListVarName{ivar_Y};


[227]  572  eval(['vec_X=reshape(Data.' XName ',[],1);'])


 573  eval(['vec_Y=reshape(Data.' YName ',[],1);'])


[201]  574  elseif numel(VarType.coord)==2 && ~isequal(VarType.coord,[0 0]);%coordinates defines by dimension variables


 575  eval(['y=Data.' Data.ListVarName{VarType.coord(1)} ';'])


 576  eval(['x=Data.' Data.ListVarName{VarType.coord(2)} ';'])


 577  if numel(y)==2 % y defined by first and last values on aregular mesh


 578  y=linspace(y(1),y(2),size(vec_U,1));


 579  end


 580  if numel(x)==2 % y defined by first and last values on aregular mesh


 581  x=linspace(x(1),x(2),size(vec_U,2));


 582  end


 583  [vec_X,vec_Y]=meshgrid(x,y);


 584  else


 585  errormsg='error in plot_field: invalid coordinate definition for vector field';


 586  return


 587  end


 588  if ~isempty(ivar_C)


 589  eval(['vec_C=Data.' Data.ListVarName{ivar_C} ';']) ;


 590  vec_C=reshape(vec_C,1,numel(vec_C));


 591  test_C=1;


 592  end


 593  if ~isempty(ivar_F)%~(isfield(PlotParam.Vectors,'HideWarning')&& isequal(PlotParam.Vectors.HideWarning,1))


 594  if test_vec


[313]  595  vec_F=Data.(Data.ListVarName{ivar_F}); % warning flags for dubious vectors


[294]  596  if ~(isfield(PlotParam.Vectors,'CheckHideWarning') && isequal(PlotParam.Vectors.CheckHideWarning,1))


[201]  597  test_black=1;


 598  end


 599  end


 600  end


 601  if ~isempty(ivar_FF) %&& ~test_false


 602  if test_vec% TODO: deal with FF for structured coordinates


[313]  603  vec_FF=Data.(Data.ListVarName{ivar_FF}); % flags for false vectors


[201]  604  end


 605  end


 606  end


 607  elseif ~isempty(ivar_C) %scalar or image


 608  if test_ima


 609  errormsg='attempt to plot two scalar fields or images';


 610  return


 611  end


 612  eval(['A=squeeze(Data.' Data.ListVarName{ivar_C} ');']) ;% scalar represented as color image


 613  test_ima=1;


[227]  614  if ~isempty(ivar_X) && ~isempty(ivar_Y)% 2D field (with unstructured coordinates (then ivar_X and ivar_Y not empty)


 615  A=reshape(A,1,[]);


[201]  616  XName=Data.ListVarName{ivar_X};


 617  YName=Data.ListVarName{ivar_Y};


[227]  618  eval(['AX=reshape(Data.' XName ',1,[]);'])


 619  eval(['AY=reshape(Data.' YName ',1,[]);'])


[201]  620  [A,AX,AY]=proj_grid(AX',AY',A',[],[],'np>256'); % interpolate on a grid


 621  if isfield(Data,'VarAttribute')


 622  if numel(Data.VarAttribute)>=ivar_X && isfield(Data.VarAttribute{ivar_X},'units')


 623  x_units=[' (' Data.VarAttribute{ivar_X}.units ')'];


 624  end


 625  if numel(Data.VarAttribute)>=ivar_Y && isfield(Data.VarAttribute{ivar_Y},'units')


 626  y_units=[' (' Data.VarAttribute{ivar_Y}.units ')'];


 627  end


 628  end


 629  elseif numel(VarType.coord)==2 %structured coordinates


 630  XName=Data.ListVarName{VarType.coord(2)};


 631  YName=Data.ListVarName{VarType.coord(1)};


 632  eval(['AY=Data.' Data.ListVarName{VarType.coord(1)} ';'])


 633  eval(['AX=Data.' Data.ListVarName{VarType.coord(2)} ';'])


 634  test_interp_X=0; %default, regularly meshed X coordinate


 635  test_interp_Y=0; %default, regularly meshed Y coordinate


 636  if isfield(Data,'VarAttribute')


 637  if numel(Data.VarAttribute)>=VarType.coord(2) && isfield(Data.VarAttribute{VarType.coord(2)},'units')


 638  x_units=Data.VarAttribute{VarType.coord(2)}.units;


 639  end


 640  if numel(Data.VarAttribute)>=VarType.coord(1) && isfield(Data.VarAttribute{VarType.coord(1)},'units')


 641  y_units=Data.VarAttribute{VarType.coord(1)}.units;


 642  end


 643  end


 644  if numel(AY)>2


 645  DAY=diff(AY);


 646  DAY_min=min(DAY);


 647  DAY_max=max(DAY);


 648  if sign(DAY_min)~=sign(DAY_max);% =1 for increasing values, 0 otherwise


 649  errormsg=['errror in plot_field.m: non monotonic dimension variable ' Data.ListVarName{VarType.coord(1)} ];


 650  return


 651  end


 652  test_interp_Y=(DAY_maxDAY_min)> 0.0001*abs(DAY_max);


 653  end


 654  if numel(AX)>2


 655  DAX=diff(AX);


 656  DAX_min=min(DAX);


 657  DAX_max=max(DAX);


 658  if sign(DAX_min)~=sign(DAX_max);% =1 for increasing values, 0 otherwise


 659  errormsg=['errror in plot_field.m: non monotonic dimension variable ' Data.ListVarName{VarType.coord(2)} ];


 660  return


 661  end


 662  test_interp_X=(DAX_maxDAX_min)> 0.0001*abs(DAX_max);


 663  end


 664  if test_interp_Y


 665  npxy(1)=max([256 floor((AY(end)AY(1))/DAY_min) floor((AY(end)AY(1))/DAY_max)]);


 666  yI=linspace(AY(1),AY(end),npxy(1));


 667  if ~test_interp_X


 668  xI=linspace(AX(1),AX(end),size(A,2));%default


 669  AX=xI;


 670  end


 671  end


 672  if test_interp_X


 673  npxy(2)=max([256 floor((AX(end)AX(1))/DAX_min) floor((AX(end)AX(1))/DAX_max)]);


 674  xI=linspace(AX(1),AX(end),npxy(2));


 675  if ~test_interp_Y


 676  yI=linspace(AY(1),AY(end),size(A,1));


 677  AY=yI;


 678  end


 679  end


 680  if test_interp_X  test_interp_Y


 681  [AX2D,AY2D]=meshgrid(AX,AY);


 682  A=interp2(AX2D,AY2D,double(A),xI,yI');


 683  end


 684  AX=[AX(1) AX(end)];% keep only the lower and upper bounds for image represnetation


 685  AY=[AY(1) AY(end)];


 686  else


 687  errormsg='error in plot_field: invalid coordinate definition ';


 688  return


 689  end


 690  end


 691  %define coordinates as CoordUnits, if not defined as attribute for each variable


 692  if isfield(Data,'CoordUnit')


 693  if isempty(x_units)


 694  x_units=Data.CoordUnit;


 695  end


 696  if isempty(y_units)


 697  y_units=Data.CoordUnit;


 698  end


 699  end


 700 


 701  end


 702 


 703  %% image or scalar plot %%%%%%%%%%%%%%%%%%%%%%%%%%


 704 


[294]  705  if isfield(PlotParam.Scalar,'ListContour')


[364]  706  CheckContour=strcmp(PlotParam.Scalar.ListContour,'contours');


[294]  707  else


 708  CheckContour=0; %default


[201]  709  end


 710  PlotParamOut=PlotParam; %default


 711  if test_ima


 712  % distinguish B/W and color images


 713  np=size(A);%size of image


 714  siz=numel(np);


 715  if siz>3


[428]  716  errormsg=['unrecognized scalar type: ' num2str(siz) ' dimensions'];


 717  return


[201]  718  end


 719  if siz==3


 720  if np(3)==1


 721  siz=2;%B W image


 722  elseif np(3)==3


 723  siz=3;%color image


 724  else


[206]  725  errormsg=['unrecognized scalar type in plot_field: considered as 2D field with ' num2str(np(3)) ' color components'];


[201]  726  return


 727  end


 728  end


 729 


 730  %set the color map


[421]  731  if isfield(PlotParam.Scalar,'CheckBW') && ~isempty(PlotParam.Scalar.CheckBW)


 732  BW=PlotParam.Scalar.CheckBW; %BW=0 color imposed, else gray scale imposed.


 733  else % BW imposed automatically chosen


[201]  734  BW=(siz==2) && (isa(A,'uint8') isa(A,'uint16'));% non color images represented in gray scale by default


[421]  735  PlotParamOut.Scalar.CheckBW=BW;


[428]  736  end


[201]  737  %case of grey level images or contour plot


[428]  738  if siz==2


[294]  739  if ~isfield(PlotParam.Scalar,'CheckFixScalar')


 740  PlotParam.Scalar.CheckFixScalar=0;%default


[201]  741  end


 742  if ~isfield(PlotParam.Scalar,'MinA')


 743  PlotParam.Scalar.MinA=[];%default


 744  end


 745  if ~isfield(PlotParam.Scalar,'MaxA')


 746  PlotParam.Scalar.MaxA=[];%default


 747  end


[206]  748  Aline=[];


[294]  749  if ~PlotParam.Scalar.CheckFixScalar isempty(PlotParam.Scalar.MinA)~isa(PlotParam.Scalar.MinA,'double') %correct if there is no numerical data in edit box


[206]  750  Aline=reshape(A,1,[]);


 751  Aline=Aline(~isnan(A));


 752  if isempty(Aline)


[428]  753  errormsg='NaN input scalar or image in plot_field';


[206]  754  return


 755  end


[210]  756  MinA=double(min(Aline));


[201]  757  else


[206]  758  MinA=PlotParam.Scalar.MinA;


[428]  759  end;


[294]  760  if ~PlotParam.Scalar.CheckFixScalarisempty(PlotParam.Scalar.MaxA)~isa(PlotParam.Scalar.MaxA,'double') %correct if there is no numerical data in edit box


[206]  761  if isempty(Aline)


[428]  762  Aline=reshape(A,1,[]);


 763  Aline=Aline(~isnan(A));


 764  if isempty(Aline)


 765  errormsg='NaN input scalar or image in plot_field';


 766  return


 767  end


[206]  768  end


[210]  769  MaxA=double(max(Aline));


[201]  770  else


[428]  771  MaxA=PlotParam.Scalar.MaxA;


 772  end;


[201]  773  PlotParamOut.Scalar.MinA=MinA;


 774  PlotParamOut.Scalar.MaxA=MaxA;


 775  % case of contour plot


[294]  776  if CheckContour


[201]  777  if ~isempty(hima) && ishandle(hima)


 778  delete(hima)


 779  end


 780  if ~isfield(PlotParam.Scalar,'IncrA')


 781  PlotParam.Scalar.IncrA=NaN;


 782  end


[294]  783  if isempty(PlotParam.Scalar.IncrA) isnan(PlotParam.Scalar.IncrA)%  PlotParam.Scalar.AutoScal==0


[201]  784  cont=colbartick(MinA,MaxA);


 785  intercont=cont(2)cont(1);%default


 786  PlotParamOut.Scalar.IncrA=intercont;


 787  else


[428]  788  intercont=PlotParam.Scalar.IncrA;


[201]  789  end


[428]  790  B=A;


[201]  791  abscontmin=intercont*floor(MinA/intercont);


 792  abscontmax=intercont*ceil(MaxA/intercont);


 793  contmin=intercont*floor(min(min(B))/intercont);


 794  contmax=intercont*ceil(max(max(B))/intercont);


 795  cont_pos_plus=0:intercont:contmax;


 796  cont_pos_min=double(contmin):intercont:intercont;


 797  cont_pos=[cont_pos_min cont_pos_plus];


 798  sizpx=(AX(end)AX(1))/(np(2)1);


 799  sizpy=(AY(1)AY(end))/(np(1)1);


[428]  800  x_cont=AX(1):sizpx:AX(end); % pixel x coordinates for image display


[201]  801  y_cont=AY(1):sizpy:AY(end); % pixel x coordinates for image display


[428]  802  % axes(haxes)% set the input axes handle as current axis


 803  txt=ver('MATLAB');


 804  Release=txt.Release;


[201]  805  relnumb=str2double(Release(3:4));


 806  if relnumb >= 14


[428]  807  vec=linspace(0,1,(abscontmaxabscontmin)/intercont);%define a greyscale colormap with steps intercont


[201]  808  map=[vec' vec' vec'];


 809  colormap(map);


[428]  810  [var,hcontour]=contour(x_cont,y_cont,B,cont_pos);


[201]  811  set(hcontour,'Fill','on')


 812  set(hcontour,'LineStyle','none')


 813  hold on


 814  end


 815  [var_p,hcontour_p]=contour(x_cont,y_cont,B,cont_pos_plus,'k');


 816  hold on


 817  [var_m,hcontour_m]=contour(x_cont,y_cont,B,cont_pos_min,':');


 818  set(hcontour_m,'LineColor',[1 1 1])


 819  hold off


[428]  820  caxis([abscontmin abscontmax])


[201]  821  colormap(map);


[428]  822  if isfield(PlotParam.Coordinates,'CheckFixAspectRatio') && isequal(PlotParam.Coordinates.CheckFixAspectRatio,1)


[415]  823  set(haxes,'DataAspectRatioMode','manual')


[428]  824  if isfield(PlotParam.Coordinates,'AspectRatio')


[429]  825  set(haxes,'DataAspectRatio',[PlotParam.Coordinates.AspectRatio 1 1])


[428]  826  else


 827  set(haxes,'DataAspectRatio',[1 1 1])


 828  end


 829  end


[201]  830  end


 831 


 832  % set colormap for image display


[294]  833  if ~CheckContour


[201]  834  % rescale the grey levels with min and max, put a grey scale colorbar


 835  B=A;


 836  if BW


 837  vec=linspace(0,1,255);%define a linear greyscale colormap


 838  map=[vec' vec' vec'];


[428]  839  colormap(map); %grey scale color map


[201]  840  else


[428]  841  colormap('default'); % standard faulse colors for div, vort , scalar fields


[201]  842  end


 843  end


 844 


[428]  845  % case of color images


 846  else


[201]  847  if BW


 848  B=uint16(sum(A,3));


 849  else


 850  B=uint8(A);


 851  end


 852  MinA=0;


 853  MaxA=255;


 854  end


 855 


 856  % display usual image


[428]  857  if ~CheckContour


[201]  858  % interpolate field to increase resolution of image display


 859  test_interp=1;


[428]  860  if max(np) <= 64


[201]  861  npxy=8*np;% increase the resolution 8 times


[428]  862  elseif max(np) <= 128


[201]  863  npxy=4*np;% increase the resolution 4 times


[428]  864  elseif max(np) <= 256


[201]  865  npxy=2*np;% increase the resolution 2 times


 866  else


 867  npxy=np;


 868  test_interp=0; % no interpolation done


 869  end


[428]  870  if test_interp==1%if we interpolate


[201]  871  x=linspace(AX(1),AX(2),np(2));


 872  y=linspace(AY(1),AY(2),np(1));


 873  [X,Y]=meshgrid(x,y);


 874  xi=linspace(AX(1),AX(2),npxy(2));


 875  yi=linspace(AY(1),AY(2),npxy(1));


 876  B = interp2(X,Y,double(B),xi,yi');


[428]  877  end


[201]  878  % create new image if there no image handle is found


[428]  879  if isempty(hima)


[201]  880  tag=get(haxes,'Tag');


 881  if MinA<MaxA


 882  hima=imagesc(AX,AY,B,[MinA MaxA]);


 883  else % to deal with uniform field


 884  hima=imagesc(AX,AY,B,[MaxA1 MaxA]);


 885  end


[428]  886  % the function imagesc reset the axes 'DataAspectRatioMode'='auto', change if .CheckFixAspectRatio is


[411]  887  % requested:


 888  set(hima,'Tag','ima')


 889  set(hima,'HitTest','off')


[428]  890  set(haxes,'Tag',tag);%preserve the axes tag (removed by image fct !!!)


[201]  891  uistack(hima, 'bottom')


[428]  892  % update an existing image


[201]  893  else


 894  set(hima,'CData',B);


 895  if MinA<MaxA


 896  set(haxes,'CLim',[MinA MaxA])


 897  else


[238]  898  set(haxes,'CLim',[MinA MaxA+1])


[201]  899  end


 900  set(hima,'XData',AX);


 901  set(hima,'YData',AY);


 902  end


[429]  903 


[405]  904  % set the transparency to 0.5 if vectors are also plotted


[428]  905  if isfield(PlotParam.Scalar,'Opacity')&& ~isempty(PlotParam.Scalar.Opacity)


 906  set(hima,'AlphaData',PlotParam.Scalar.Opacity)


[405]  907  else


[428]  908  if test_vec


 909  set(hima,'AlphaData',0.5)%set opacity to 0.5 by default in the presence of vectors


 910  PlotParamOut.Scalar.Opacity=0.5;


 911  else


 912  set(hima,'AlphaData',1)% full opacity (no transparency) by default


 913  end


[405]  914  end


[201]  915  end


 916  test_ima=1;


 917 


 918  %display the colorbar code for B/W images if Poscolorbar not empty


 919  if siz==2 && exist('PosColorbar','var')&& ~isempty(PosColorbar)


 920  if isempty(hcol)~ishandle(hcol)


[428]  921  hcol=colorbar;%create new colorbar


[201]  922  end


 923  if length(PosColorbar)==4


[428]  924  set(hcol,'Position',PosColorbar)


 925  end


[201]  926  %YTick=0;%default


 927  if MaxA>MinA


[294]  928  if CheckContour


[201]  929  colbarlim=get(hcol,'YLim');


[428]  930  scale_bar=(colbarlim(2)colbarlim(1))/(abscontmaxabscontmin);


[201]  931  YTick=cont_pos(2:end1);


 932  YTick_scaled=colbarlim(1)+scale_bar*(YTickabscontmin);


 933  set(hcol,'YTick',YTick_scaled);


[294]  934  elseif (isfield(PlotParam.Scalar,'CheckBW') && isequal(PlotParam.Scalar.CheckBW,1))isa(A,'uint8') isa(A,'uint16')%images


[201]  935  hi=get(hcol,'children');


 936  if iscell(hi)%multiple images in colorbar


 937  hi=hi{1};


 938  end


 939  set(hi,'YData',[MinA MaxA])


 940  set(hi,'CData',(1:256)')


 941  set(hcol,'YLim',[MinA MaxA])


 942  YTick=colbartick(MinA,MaxA);


[428]  943  set(hcol,'YTick',YTick)


[201]  944  else


 945  hi=get(hcol,'children');


 946  if iscell(hi)%multiple images in colorbar


 947  hi=hi{1};


 948  end


 949  set(hi,'YData',[MinA MaxA])


 950  set(hi,'CData',(1:64)')


[428]  951  YTick=colbartick(MinA,MaxA);


[201]  952  set(hcol,'YLim',[MinA MaxA])


 953  set(hcol,'YTick',YTick)


 954  end


 955  set(hcol,'Yticklabel',num2str(YTick'));


 956  end


 957  elseif ishandle(hcol)


[428]  958  delete(hcol); %erase existing colorbar if not needed


[201]  959  end


 960  else%no scalar plot


[428]  961  if ~isempty(hima) && ishandle(hima)


[201]  962  delete(hima)


 963  end


 964  if ~isempty(hcol)&& ishandle(hcol)


[428]  965  delete(hcol)


[201]  966  end


 967  PlotParamOut=rmfield(PlotParamOut,'Scalar');


 968  end


 969 


 970  %% vector plot %%%%%%%%%%%%%%%%%%%%%%%%%%


 971  if test_vec


 972  %vector scale representation


 973  if size(vec_U,1)==numel(vec_Y) && size(vec_U,2)==numel(vec_X); % x, y coordinate variables


 974  [vec_X,vec_Y]=meshgrid(vec_X,vec_Y);


 975  end


 976  vec_X=reshape(vec_X,1,numel(vec_X));%reshape in matlab vectors


 977  vec_Y=reshape(vec_Y,1,numel(vec_Y));


 978  vec_U=reshape(vec_U,1,numel(vec_U));


 979  vec_V=reshape(vec_V,1,numel(vec_V));


 980  MinMaxX=max(vec_X)min(vec_X);


[294]  981  if isfield(PlotParam.Vectors,'CheckFixVectors') && isequal(PlotParam.Vectors.CheckFixVectors,1)&& isfield(PlotParam.Vectors,'VecScale')...


[210]  982  &&~isempty(PlotParam.Vectors.VecScale) && isa(PlotParam.Vectors.VecScale,'double') %fixed vector scale


 983  scale=PlotParam.Vectors.VecScale; %impose the length of vector representation


 984  else


 985  if ~test_false %remove false vectors


[201]  986  indsel=1:numel(vec_X);%


 987  end


 988  if isempty(vec_U)


 989  scale=1;


 990  else


 991  if isempty(indsel)


 992  MaxU=max(abs(vec_U));


 993  MaxV=max(abs(vec_V));


 994  else


 995  MaxU=max(abs(vec_U(indsel)));


 996  MaxV=max(abs(vec_V(indsel)));


 997  end


 998  scale=MinMaxX/(max(MaxU,MaxV)*50);


 999  PlotParam.Vectors.VecScale=scale;%update the 'scale' display


 1000  end


[210]  1001  end


[201]  1002 


 1003  %record vectors on the plotting axes


 1004  if test_C==0


 1005  vec_C=ones(1,numel(vec_X));


 1006  end


 1007 


 1008  %decimate by a factor 2 in vector mesh(4 in nbre of vectors)


[294]  1009  if isfield(PlotParam.Vectors,'CheckDecimate4') && PlotParam.Vectors.CheckDecimate4


[201]  1010  diffy=diff(vec_Y); %difference dy=vec_Y(i+1)vec_Y(i)


 1011  dy_thresh=max(abs(diffy))/2;


 1012  ind_jump=find(abs(diffy) > dy_thresh); %indices with diff(vec_Y)> max/2, detect change of line


 1013  ind_sel=1:ind_jump(1);%select the first line


 1014  for i=2:2:length(ind_jump)1


 1015  ind_sel=[ind_sel (ind_jump(i)+1:ind_jump(i+1))];% select the odd lines


 1016  end


 1017  nb_sel=length(ind_sel);


 1018  ind_sel=ind_sel(1:2:nb_sel);% take half the points on a line


 1019  vec_X=vec_X(ind_sel);


 1020  vec_Y=vec_Y(ind_sel);


 1021  vec_U=vec_U(ind_sel);


 1022  vec_V=vec_V(ind_sel);


 1023  vec_C=vec_C(ind_sel);


 1024  if ~isempty(ivar_F)


 1025  vec_F=vec_F(ind_sel);


 1026  end


 1027  if ~isempty(ivar_FF)


 1028  vec_FF=vec_FF(ind_sel);


 1029  end


 1030  end


 1031 


 1032  %get main level color code


 1033  [colorlist,col_vec,PlotParamOut.Vectors]=set_col_vec(PlotParam.Vectors,vec_C);


 1034 


 1035  % take flags into account: add flag colors to the list of colors


 1036  sizlist=size(colorlist);


 1037  nbcolor=sizlist(1);


 1038  if test_black


 1039  nbcolor=nbcolor+1;


 1040  colorlist(nbcolor,:)=[0 0 0]; %add black to the list of colors


 1041  if ~isempty(ivar_FF)


[210]  1042  % ind_flag=find(vec_F~=1 & vec_F~=0 & vec_FF==0); %flag warning but not false


 1043  col_vec(vec_F~=1 & vec_F~=0 & vec_FF==0)=nbcolor;


[201]  1044  else


[210]  1045  col_vec(vec_F~=1 & vec_F~=0)=nbcolor;


[201]  1046  end


 1047  end


 1048  nbcolor=nbcolor+1;


 1049  if ~isempty(ivar_FF)


[294]  1050  if isfield(PlotParam.Vectors,'CheckHideFalse') && PlotParam.Vectors.CheckHideFalse==1


[201]  1051  colorlist(nbcolor,:)=[NaN NaN NaN];% no plot of false vectors


 1052  else


 1053  colorlist(nbcolor,:)=[1 0 1];% magenta color


 1054  end


[210]  1055  col_vec(vec_FF~=0)=nbcolor;


[201]  1056  end


 1057  %plot vectors:


 1058  quiresetn(haxes,vec_X,vec_Y,vec_U,vec_V,scale,colorlist,col_vec);


 1059 


 1060  else


 1061  hvec=findobj(haxes,'Tag','vel');


 1062  if ~isempty(hvec)


 1063  delete(hvec);


 1064  end


 1065  PlotParamOut=rmfield(PlotParamOut,'Vectors');


 1066  end


 1067 


 1068  %listfields={'AY','AX','A','X','Y','U','V','C','W','F','FF'};


 1069  %listdim={'AY','AX',{'AY','AX'},'nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors'};


 1070  %Role={'coord_y','coord_x','scalar','coord_x','coord_y','vector_x','vector_y','scalar','vector_z','warnflag','errorflag'};


 1071  %ind_select=[];


 1072  nbvar=0;


 1073 


 1074  %store the coordinate extrema occupied by the field


 1075  if ~isempty(Data)


[252]  1076  XMin=[];


 1077  XMax=[];


 1078  YMin=[];


 1079  YMax=[];


[334]  1080  fix_lim=isfield(PlotParam.Coordinates,'CheckFixLimits') && PlotParam.Coordinates.CheckFixLimits;


[201]  1081  if fix_lim


[334]  1082  if isfield(PlotParam.Coordinates,'MinX')&&isfield(PlotParam.Coordinates,'MaxX')&&isfield(PlotParam.Coordinates,'MinY')&&isfield(PlotParam.Coordinates,'MaxY')


 1083  XMin=PlotParam.Coordinates.MinX;


 1084  XMax=PlotParam.Coordinates.MaxX;


 1085  YMin=PlotParam.Coordinates.MinY;


 1086  YMax=PlotParam.Coordinates.MaxY;


[201]  1087  end %else PlotParamOut.XMin =PlotParam.XMin...


[252]  1088  else


 1089  if test_ima %both background image and vectors coexist, take the wider bound


[201]  1090  XMin=min(AX);


 1091  XMax=max(AX);


 1092  YMin=min(AY);


 1093  YMax=max(AY);


[252]  1094  if test_vec


 1095  XMin=min(XMin,min(vec_X));


 1096  XMax=max(XMax,max(vec_X));


 1097  YMin=min(YMin,min(vec_Y));


 1098  YMax=max(YMax,max(vec_Y));


 1099  end


 1100  elseif test_vec


 1101  XMin=min(vec_X);


 1102  XMax=max(vec_X);


 1103  YMin=min(vec_Y);


 1104  YMax=max(vec_Y);


[201]  1105  end


[252]  1106  end


 1107  % PlotParamOut.RangeX=[XMin XMax]; %range of x, to be stored in the user data of the plot axes


 1108  % PlotParamOut.RangeY=[YMin YMax]; %range of x, to be stored in the user data of the plot axes


 1109  % if ~fix_lim


[334]  1110  PlotParamOut.Coordinates.MinX=XMin;


 1111  PlotParamOut.Coordinates.MaxX=XMax;


 1112  PlotParamOut.Coordinates.MinY=YMin;


 1113  PlotParamOut.Coordinates.MaxY=YMax;


[201]  1114  if XMax>XMin


 1115  set(haxes,'XLim',[XMin XMax]);% set x limits of frame in axes coordinates


 1116  end


 1117  if YMax>YMin


 1118  set(haxes,'YLim',[YMin YMax]);% set x limits of frame in axes coordinates


 1119  end


[252]  1120  % end


[201]  1121  set(haxes,'YDir','normal')


 1122  set(get(haxes,'XLabel'),'String',[XName ' (' x_units ')']);


 1123  set(get(haxes,'YLabel'),'String',[YName ' (' y_units ')']);


[334]  1124  PlotParamOut.Coordinates.x_units=x_units;


 1125  PlotParamOut.Coordinates.y_units=y_units;


[201]  1126  end


[429]  1127  if isfield(PlotParam,'Coordinates') && isfield(PlotParam.Coordinates,'CheckFixAspectRatio') && isequal(PlotParam.Coordinates.CheckFixAspectRatio,1)


 1128  set(haxes,'DataAspectRatioMode','manual')


 1129  if isfield(PlotParam.Coordinates,'AspectRatio')


 1130  set(haxes,'DataAspectRatio',[PlotParam.Coordinates.AspectRatio 1 1])


 1131  end


 1132  else


 1133  set(haxes,'DataAspectRatioMode','auto')


 1134  end


[201]  1135  %


 1136  %  function for plotting vectors


 1137  %INPUT:


 1138  % haxes: handles of the plotting axes


 1139  % x,y,u,v: vectors coordinates and vector components to plot, arrays withb the same dimension


 1140  % scale: scaling factor for vector length representation


 1141  % colorlist(icolor,:): list of vector colors, dim (nbcolor,3), depending on color #i


 1142  % col_vec: matlab vector setting the color number #i for each velocity vector


 1143  function quiresetn(haxes,x,y,u,v,scale,colorlist,col_vec)


 1144  %


 1145  %define arrows


 1146  theta=0.5 ;%angle arrow


 1147  alpha=0.3 ;%length arrow


 1148  rot=alpha*[cos(theta) sin(theta); sin(theta) cos(theta)]';


 1149  %find the existing lines


 1150  h=findobj(haxes,'Tag','vel');% search existing lines in the current axes


 1151  sizh=size(h);


 1152  set(h,'EraseMode','xor');


 1153  set(haxes,'NextPlot','replacechildren');


 1154 


 1155  %drawnow


 1156  %create lines (if no lines) or modify them


 1157  if ~isequal(size(col_vec),size(x))


 1158  col_vec=ones(size(x));% case of error in col_vec input


 1159  end


 1160  sizlist=size(colorlist);


 1161  ncolor=sizlist(1);


 1162 


 1163  for icolor=1:ncolor


 1164  %determine the line positions for each color icolor


 1165  ind=find(col_vec==icolor);


 1166  xc=x(ind);


 1167  yc=y(ind);


 1168  uc=u(ind)*scale;


 1169  vc=v(ind)*scale;


 1170  n=size(xc);


 1171  xN=NaN*ones(size(xc));


 1172  matx=[xc(:)uc(:)/2 xc(:)+uc(:)/2 xN(:)]';


 1173  matx=reshape(matx,1,3*n(2));


 1174  maty=[yc(:)vc(:)/2 yc(:)+vc(:)/2 xN(:)]';


 1175  maty=reshape(maty,1,3*n(2));


 1176 


 1177  %determine arrow heads


 1178  arrowplus=rot*[uc;vc];


 1179  arrowmoins=rot'*[uc;vc];


 1180  x1=xc+uc/2arrowplus(1,:);


 1181  x2=xc+uc/2;


 1182  x3=xc+uc/2arrowmoins(1,:);


 1183  y1=yc+vc/2arrowplus(2,:);


 1184  y2=yc+vc/2;


 1185  y3=yc+vc/2arrowmoins(2,:);


 1186  matxar=[x1(:) x2(:) x3(:) xN(:)]';


 1187  matxar=reshape(matxar,1,4*n(2));


 1188  matyar=[y1(:) y2(:) y3(:) xN(:)]';


 1189  matyar=reshape(matyar,1,4*n(2));


 1190  %draw the line or modify the existing ones


[421]  1191  % tri=reshape(1:3*length(uc),3,[])';


[201]  1192  isn=isnan(colorlist(icolor,:));%test if color NaN


 1193  if 2*icolor > sizh(1) %if icolor exceeds the number of existing ones


 1194  if ~isn(1) %if the vectors are visible color not nan


 1195  if n(2)>0


 1196  hold on


 1197  line(matx,maty,'Color',colorlist(icolor,:),'Tag','vel');% plot new lines


 1198  line(matxar,matyar,'Color',colorlist(icolor,:),'Tag','vel');% plot arrows


 1199  end


 1200  end


 1201  else


 1202  if isn(1)


 1203  delete(h(2*icolor1))


 1204  delete(h(2*icolor))


 1205  else


 1206  set(h(2*icolor1),'Xdata',matx,'Ydata',maty);


 1207  set(h(2*icolor1),'Color',colorlist(icolor,:));


 1208  set(h(2*icolor1),'EraseMode','xor');


 1209  set(h(2*icolor),'Xdata',matxar,'Ydata',matyar);


 1210  set(h(2*icolor),'Color',colorlist(icolor,:));


 1211  set(h(2*icolor),'EraseMode','xor');


 1212  end


 1213  end


 1214  end


 1215  if sizh(1) > 2*ncolor


 1216  for icolor=ncolor+1 : sizh(1)/2%delete additional objects


 1217  delete(h(2*icolor1))


 1218  delete(h(2*icolor))


 1219  end


 1220  end


 1221 


 1222  %


 1223  %  determine tick positions for colorbar


 1224  function YTick=colbartick(MinA,MaxA)


 1225  %


 1226  %determine tick positions with "simple" values between MinA and MaxA


 1227  YTick=0;%default


 1228  maxabs=max([abs(MinA) abs(MaxA)]);


 1229  if maxabs>0


 1230  ord=10^(floor(log10(maxabs)));%order of magnitude


 1231  div=1;


 1232  siz2=1;


 1233  while siz2<2


 1234  values=10:div:10;


 1235  ind=find((ord*valuesMaxA)<0 & (ord*valuesMinA)>0);%indices of 'values' such that MinA<ord*values<MaxA


 1236  siz=size(ind);


 1237  if siz(2)<4%if there are less than 4 selected values (4 levels)


 1238  values=9:0.5*div:9;


 1239  ind=find((ord*valuesMaxA)<0 & (ord*valuesMinA)>0);


 1240  end


 1241  siz2=size(ind,2);


 1242  div=div/10;


 1243  end


 1244  YTick=ord*values(ind);


 1245  end


 1246 


[356]  1247  % 


[389]  1248  %  'proj_grid': project fields with unstructured coordinantes on a regular grid


[356]  1249  function [A,rangx,rangy]=proj_grid(vec_X,vec_Y,vec_A,rgx_in,rgy_in,npxy_in)


[389]  1250  % 


[356]  1251  if length(vec_Y)<2


 1252  msgbox_uvmat('ERROR','less than 2 points in proj_grid.m');


 1253  return;


 1254  end


 1255  diffy=diff(vec_Y); %difference dy=vec_Y(i+1)vec_Y(i)


 1256  index=find(diffy);% find the indices of vec_Y after wich a change of horizontal line occurs(diffy non zero)


 1257  if isempty(index); msgbox_uvmat('ERROR','points aligned along abscissa in proj_grid.m'); return; end;%points aligned% A FAIRE: switch to line plot.


 1258  diff2=diff(diffy(index));% diff2 = fluctuations of the detected vertical grid mesh dy


 1259  if max(abs(diff2))>0.001*abs(diffy(index(1))) % if max(diff2) is larger than 1/1000 of the first mesh dy


 1260  % the data are not regularly spaced and must be interpolated on a regular grid


 1261  if exist('rgx_in','var') & ~isempty (rgx_in) & isnumeric(rgx_in) & length(rgx_in)==2% positions imposed from input


 1262  rangx=rgx_in; % first and last positions


 1263  rangy=rgy_in;


 1264  dxy(1)=1/(npxy_in(1)1);%grid mesh in y


 1265  dxy(2)=1/(npxy_in(2)1);%grid mesh in x


 1266  dxy(1)=(rangy(2)rangy(1))/(npxy_in(1)1);%grid mesh in y


 1267  dxy(2)=(rangx(2)rangx(1))/(npxy_in(2)1);%grid mesh in x


 1268  else % interpolation grid automatically determined


 1269  rangx(1)=min(vec_X);


 1270  rangx(2)=max(vec_X);


 1271  rangy(2)=min(vec_Y);


 1272  rangy(1)=max(vec_Y);


 1273  dxymod=sqrt((rangx(2)rangx(1))*(rangy(1)rangy(2))/length(vec_X));


 1274  dxy=[dxymod/4 dxymod/4];% increase the resolution 4 times


 1275  end


 1276  xi=[rangx(1):dxy(2):rangx(2)];


 1277  yi=[rangy(1):dxy(1):rangy(2)];


 1278  A=griddata_uvmat(vec_X,vec_Y,vec_A,xi,yi');


 1279  A=reshape(A,length(yi),length(xi));


 1280  else


 1281  x=vec_X(1:index(1));% the set of abscissa (obtained on the first line)


 1282  indexend=index(end);% last vector index of line change


 1283  ymax=vec_Y(indexend+1);% y coordinate AFTER line change


 1284  ymin=vec_Y(index(1));


 1285  y=vec_Y(index);


 1286  y(length(y)+1)=ymax;


 1287  nx=length(x); %number of grid points in x


 1288  ny=length(y); % number of grid points in y


 1289  B=(reshape(vec_A,nx,ny))'; %vec_A reshaped as a rectangular matrix


 1290  [X,Y]=meshgrid(x,y);% positions X and Y also reshaped as matrix


 1291 


 1292  %linear interpolation to improve the image resolution and/or adjust


 1293  %to prescribed positions


 1294  test_interp=1;


 1295  if exist('rgx_in','var') & ~isempty (rgx_in) & isnumeric(rgx_in) & length(rgx_in)==2% positions imposed from input


 1296  rangx=rgx_in; % first and last positions


 1297  rangy=rgy_in;


 1298  npxy=npxy_in;


 1299  else


 1300  rangx=[vec_X(1) vec_X(nx)];% first and last position found for x


 1301  rangy=[max(ymax,ymin) min(ymax,ymin)];


 1302  if max(nx,ny) <= 64 & isequal(npxy_in,'np>256')


 1303  npxy=[8*ny 8*nx];% increase the resolution 8 times


 1304  elseif max(nx,ny) <= 128 & isequal(npxy_in,'np>256')


 1305  npxy=[4*ny 4*nx];% increase the resolution 4 times


 1306  elseif max(nx,ny) <= 256 & isequal(npxy_in,'np>256')


 1307  npxy=[2*ny 2*nx];% increase the resolution 2 times


 1308  else


 1309  npxy=[ny nx];


 1310  test_interp=0; % no interpolation done


 1311  end


 1312  end


 1313  if test_interp==1%if we interpolate


 1314  xi=[rangx(1):(rangx(2)rangx(1))/(npxy(2)1):rangx(2)];


 1315  yi=[rangy(1):(rangy(2)rangy(1))/(npxy(1)1):rangy(2)];


 1316  [XI,YI]=meshgrid(xi,yi);


 1317  A = interp2(X,Y,B,XI,YI);


 1318  else %no interpolation for a resolution higher than 256


 1319  A=B;


 1320  end


 1321  end 
