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


 2  %


 3  %


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


 5  % appropriate to the input field structure:


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


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


[512]  8  % then split into blocks of related variables by find_field_cells.m.


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


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


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


 12  % dimension variables (case of matrices).


 13  %


[625]  14  % function [PlotType,PlotParamOut,haxes]= plot_field(Data,haxes,PlotParam,PosColorbar)


[201]  15  %


[822]  16  % OUPUT:


[866]  17  % PlotType: type of plot: 'text','line'(curve plot),'plane':2D view,'volume', or errormsg


[201]  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  %


[822]  21  %INPUT


[201]  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:


[748]  47  % .Axes: coordinate parameters:


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


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


 50  % .Axes.AspectRatio: imposed aspect ratio y/x of axis unit plots


[201]  51  % scalars


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


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


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


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


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


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


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


[201]  59  % .IncrA : contour interval


 60  %  vectors


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


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


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


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


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


[201]  66  %  vector color


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


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


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


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


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


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


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


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


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


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


 77  %


[546]  78  % PosColorbar: % if absent, no action on colorbar


 79  % % if empty, suppress any existing colorbar


 80  % % if not empty, display a colorbar for B&W images at position PosColorbar


 81  % expressed in figure relative unit (ex [0.821 0.471 0.019 0.445])


[201]  82 


[924]  83  %=======================================================================


[1061]  84  % Copyright 20082019, LEGI UMR 5519 / CNRS UGA GINP, Grenoble, France


[924]  85  % http://www.legi.grenobleinp.fr


 86  % Joel.Sommeria  Joel.Sommeria (A) legi.cnrs.fr


 87  %


[201]  88  % This file is part of the toolbox UVMAT.


[924]  89  %


[201]  90  % UVMAT is free software; you can redistribute it and/or modify


[924]  91  % it under the terms of the GNU General Public License as published


 92  % by the Free Software Foundation; either version 2 of the license,


 93  % or (at your option) any later version.


 94  %


[201]  95  % UVMAT is distributed in the hope that it will be useful,


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


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


[924]  98  % GNU General Public License (see LICENSE.txt) for more details.


 99  %=======================================================================


[201]  100 


[781]  101  function [PlotType,PlotParamOut,haxes]= plot_field(Data,haxes,PlotParam)


[247]  102 


[388]  103  %% default input and output


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


 105  PlotType='text'; %default


[749]  106  if ~isfield(PlotParam,'Axes')


[748]  107  PlotParam.Axes=[];


[582]  108  if isfield(Data,'CoordUnit')


[748]  109  PlotParam.Axes.CheckFixAspectRatio=1;


 110  PlotParam.Axes.AspectRatio=1; %set axes equal by default if CoordUnit is defined


[582]  111  end


[292]  112  end


[569]  113  PlotParamOut=PlotParam;%default


[201]  114 


[388]  115  %% check input structure


 116  % check the cells of fields :


[530]  117 


[1045]  118  % if ~isfield(PlotParam,'FieldName')


 119  % index_0D=[];


 120  % index_1D=1;


 121  % index_2D=[];%find 2D fields


 122  % index_3D=[];


 123  % else


 124  [CellInfo,NbDimArray,errormsg]=find_field_cells(Data);


 125  if ~isempty(errormsg)


 126  msgbox_uvmat('ERROR',['input of plot_field/find_field_cells: ' errormsg]);


 127  return


 128  end


 129  index_0D=find(NbDimArray==0);


 130  index_1D=find(NbDimArray==1);


 131  index_2D=find(NbDimArray==2);%find 2D fields


 132  index_3D=find(NbDimArray>2,1);


 133  if ~isempty(index_3D)


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


 135  return


 136  end


 137  % end


[294]  138 


[388]  139  %% test axes and figure


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


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


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


 143  if ishandle(haxes)


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


 145  testnewfig=0;


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


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


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


 149  testzoomaxes=1;


 150  zoomaxes=AxeData.ZoomAxes;


[388]  151  end


 152  end


 153  end


[201]  154  end


[429]  155  end


 156 


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


 158  if testnewfig


 159  hfig=figure;


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


 161  haxes=axes;


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


[569]  163  PlotParamOut.NextPlot='add'; %parameter for plot_profile


[429]  164  else


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


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


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


 168  end


 169 


 170  %% set axes properties


[748]  171  if isfield(PlotParamOut.Axes,'CheckFixLimits') && isequal(PlotParamOut.Axes.CheckFixLimits,1) %adjust the graph limits


[429]  172  set(haxes,'XLimMode', 'manual')


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


 174  else


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


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


 177  end


[760]  178 


 179  if isfield(PlotParamOut.Axes,'CheckFixAspectRatio') && isequal(PlotParamOut.Axes.CheckFixAspectRatio,1)&&isfield(PlotParamOut.Axes,'AspectRatio')


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


 181  set(haxes,'DataAspectRatio',[PlotParamOut.Axes.AspectRatio 1 1])


 182  else


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


 184  end


 185 


[429]  186  errormsg='';


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


 188 


 189  %% 2D plots


 190  if isempty(index_2D)


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


[429]  192  else %plot 2D field


[781]  193  % if ~exist('PosColorbar','var'),PosColorbar=[];end;


 194  [tild,PlotParamOut,PlotType,errormsg]=plot_plane(Data,CellInfo(index_2D),haxes,PlotParamOut);


[429]  195  AxeData.NbDim=2;


 196  if testzoomaxes && isempty(errormsg)


[781]  197  [zoomaxes,PlotParamOut,tild,errormsg]=plot_plane(Data,CellInfo(index_2D),zoomaxes,PlotParamOut);


[429]  198  AxeData.ZoomAxes=zoomaxes;


[201]  199  end


[429]  200  end


 201 


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


[1031]  203  if isempty(index_1D) ~isempty(index_2D)


[429]  204  if ~isempty(haxes)


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


[201]  206  end


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


[644]  208  CheckHold=0;


 209  if isfield(PlotParam,'CheckHold')


 210  CheckHold= PlotParam.CheckHold;


 211  end


[874]  212  PlotParamOut=plot_profile(Data,CellInfo(index_1D),haxes,PlotParamOut,CheckHold);%


[872]  213  if isempty(index_2D)


 214  if isfield(PlotParamOut,'Vectors')


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


 216  end


 217  if isfield(PlotParamOut,'Scalar')


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


 219  end


[871]  220  end


[429]  221  if testzoomaxes


[748]  222  [zoomaxes,PlotParamOut.Axes]=plot_profile(Data,CellInfo(index_1D),zoomaxes,PlotParamOut.Axes,CheckHold);


[429]  223  AxeData.ZoomAxes=zoomaxes;


[201]  224  end


[429]  225  PlotType='line';


 226  end


 227 


[760]  228  %% aspect ratio


 229  AspectRatio=get(haxes,'DataAspectRatio');


 230  PlotParamOut.Axes.AspectRatio=AspectRatio(1)/AspectRatio(2);


 231 


[429]  232  %% text display


[873]  233  if ~(isfield(PlotParamOut,'Axes')&&isfield(PlotParamOut.Axes,'TextDisplay')&&(PlotParamOut.Axes.TextDisplay)) % if text is not already given as statistics


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


[874]  235  if ~isempty(htext)%&&~isempty(hchecktable)


[873]  236  if isempty(index_0D)


[434]  237  else


[873]  238  errormsg=plot_text(Data,CellInfo(index_0D),htext);


 239  set(htext,'visible','on')


[434]  240  end


[873]  241  set(hfig,'Unit','pixels');


 242  set(htext,'Unit','pixels')


 243  PosFig=get(hfig,'Position');


 244  % case of no plot with view_field: only text display


 245  if strcmp(get(hfig,'Tag'),'view_field')


 246  if isempty(index_1D) && isempty(index_2D)% case of no plot: only text display


 247  set(haxes,'Visible','off')


 248  PosTable=get(htext,'Position');


 249  set(hfig,'Position',[PosFig(1) PosFig(2) PosTable(3) PosTable(4)])


 250  else


 251  set(haxes,'Visible','on')


 252  set(hfig,'Position',[PosFig(1) PosFig(2) 877 677])%default size for view_field


 253  end


 254  end


[294]  255  end


[201]  256  end


[388]  257  %% display error message


[201]  258  if ~isempty(errormsg)


[866]  259  PlotType=errormsg;


[201]  260  msgbox_uvmat('ERROR', errormsg)


 261  end


 262 


 263  %% update the parameters stored in AxeData


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


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


[569]  266  AxeData.RangeX=[PlotParamOut.MinX PlotParamOut.MaxX];


 267  AxeData.RangeY=[PlotParamOut.MinY PlotParamOut.MaxY];


[388]  268  end


 269  set(haxes,'UserData',AxeData)


 270  end


[201]  271 


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


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


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


[511]  275  if strcmp(get(hfig,'tag'),'view_field')strcmp(get(hfig,'tag'),'uvmat')


 276  FigData.(get(haxes,'tag'))=Data;


[429]  277  set(hfig,'UserData',FigData)


 278  end


[388]  279  end


[292]  280 


[201]  281  %


[690]  282  %  plot 0D fields: display data values without plot


 283  %


[530]  284  function errormsg=plot_text(FieldData,CellInfo,htext)


[690]  285 


[581]  286  errormsg='';


[201]  287  txt_cell={};


[429]  288  Data={};


[690]  289  VarIndex=[];


[530]  290  for icell=1:length(CellInfo)


[535]  291 


 292  % select types of variables to be projected


 293  ListProj={'VarIndex_scalar','VarIndex_image','VarIndex_color','VarIndex_vector_x','VarIndex_vector_y'};


 294  check_proj=false(size(FieldData.ListVarName));


 295  for ilist=1:numel(ListProj)


 296  if isfield(CellInfo{icell},ListProj{ilist})


 297  check_proj(CellInfo{icell}.(ListProj{ilist}))=1;


 298  end


 299  end


[690]  300  VarIndex=[VarIndex find(check_proj)];


 301  end


 302 


 303  % data need to be displayed in a table


[881]  304  % if strcmp(get(htext,'Type'),'uitable')% display data in a table


 305  % VarNameCell=cell(1,numel(VarIndex));% prepare list of variable names to display (titles of columns)


 306  % VarLength=zeros(1,numel(VarIndex)); % default number of values for each variable


 307  % for ivar=1:numel(VarIndex)


 308  % VarNameCell{ivar}=FieldData.ListVarName{VarIndex(ivar)};


 309  % VarLength(ivar)=numel(FieldData.(VarNameCell{ivar}));


 310  % end


 311  % set(htext,'ColumnName',VarNameCell)


 312  % Data=cell(max(VarLength),numel(VarIndex));% prepare the table of data display


 313  %


 314  % for ivar=1:numel(VarIndex)


 315  % VarValue=FieldData.(VarNameCell{ivar});


 316  % VarValue=reshape(VarValue,[],1);% reshape values array in a column


 317  % Data(1:numel(VarValue),ivar)=num2cell(VarValue);


 318  % end


 319  % set(htext,'Data',Data)


 320  % end


[690]  321  % if numel(VarValue)>1 && numel(VarValue)<10 % case of a variable with several values


 322  % for ind=1:numel(VarValue)


 323  % VarNameCell{1,ind}=[VarName '_' num2str(ind)];% indicate each value by an index


 324  % end


 325  % else


 326  % VarNameCell={VarName};


 327  % end


 328  % if numel(VarValue)<10


 329  % if isempty(VarValue)


 330  % VarValueCell={'[]'};


 331  % else


 332  % VarValueCell=num2cell(VarValue);


 333  % end


 334  % if isempty(Data)


 335  % Data =[VarNameCell VarValueCell];


 336  % else


 337  % Data =[Data [VarNameCell VarValueCell]];


 338  % end


 339  % else


 340  % if isempty(Data)


 341  % Data =[VarNameCell; num2cell(VarValue)];


 342  % else


 343  % Data =[Data [VarNameCell; {['size ' num2str(size(VarValue))]}]];


 344  % end


 345  % end


 346  % if size(VarValue,1)==1


 347  % txt=[VarName '=' num2str(VarValue)];


 348  % txt_cell=[txt_cell;{txt}];


 349  % end


 350  % end


 351  % end


 352  % if strcmp(get(htext,'Type'),'uitable')% display data in a table


 353  %


 354  %


 355  % set(htext,'Data',Data(2:end,:))


 356  % else % display in a text edit box


 357  % set(htext,'String',txt_cell)


 358  % set(htext,'UserData',txt_cell)% for temporary storage when the edit box is used for mouse display


 359  % end


[201]  360 


[429]  361 


[201]  362  %


[690]  363  %  plot 1D fields (usual x,y plots)


 364  %


[874]  365  function PlotParamOut=plot_profile(data,CellInfo,haxes,PlotParam,CheckHold)


[201]  366 


[644]  367  %% initialization


[874]  368  if ~(exist('PlotParam','var')&&~isempty(PlotParam.Axes))


[292]  369  Coordinates=[];


[874]  370  PlotParamOut.Axes=Coordinates;


 371  else


 372  Coordinates=PlotParam.Axes;


 373  PlotParamOut=PlotParam;


[201]  374  end


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


[644]  376  legend_str={};


 377 


 378  %% suppress existing plot if empty data


[201]  379  if isempty(data)


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


 381  if ~isempty(hplot)


 382  delete(hplot)


 383  end


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


 385  if ~isempty(hlegend)


 386  delete(hlegend)


 387  end


 388  return


 389  end


 390 


[644]  391  %% set the colors of the successive plots (designed to produce rgb for the three components of color images)


[651]  392  ColorOrder=[1 0 0;0 1 0;0 0 1;0 0.75 0.75;0.75 0 0.75;0.75 0.75 0;0.25 0.25 0.25];


 393  set(hfig,'DefaultAxesColorOrder',ColorOrder)


[644]  394  if CheckHold


[690]  395  set(haxes,'NextPlot','add')


[644]  396  else


 397  set(haxes,'NextPlot','replace')


[201]  398  end


 399 


 400  %% prepare the string for plot command


 401  plotstr='hhh=plot(';


 402  xtitle='';


 403  ytitle='';


[644]  404  test_newplot=~CheckHold;


[569]  405  MinX=[];


 406  MaxX=[];


 407  MinY_cell=[];


 408  MaxY_cell=[];


[871]  409  testplot=ones(size(data.ListVarName));%default test for plotted variables


[201]  410  %loop on input fields


[530]  411  for icell=1:numel(CellInfo)


[1048]  412  VarIndex=[CellInfo{icell}.YIndex CellInfo{icell}.YIndex_discrete];% indices of the selected variables in the list data.ListVarName


[1045]  413  coord_x_index=CellInfo{icell}.XIndex;


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


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


[569]  416  if isempty(find(strcmp(coord_x_name{icell},coord_x_name(1:end1)), 1)) %xtitle not already selected


[426]  417  xtitle=[xtitle coord_x_name{icell}];


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


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


 420  else


 421  xtitle=[xtitle ', '];


 422  end


[201]  423  end


[729]  424  if ~isempty(coord_x{icell})


 425  MinX(icell)=min(coord_x{icell});


 426  MaxX(icell)=max(coord_x{icell});


 427  testplot(coord_x_index)=0;


 428  if isfield(CellInfo{icell},'VarIndex_ancillary')


 429  testplot(CellInfo{icell}.VarIndex_ancillary)=0;


 430  end


 431  if isfield(CellInfo{icell},'VarIndex_warnflag')


 432  testplot(CellInfo{icell}.VarIndex_warnflag)=0;


 433  end


 434  if isfield(data,'VarAttribute')


 435  VarAttribute=data.VarAttribute;


 436  for ivar=1:length(VarIndex)


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


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


 439  else


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


 441  end


[201]  442  end


 443  end


[1048]  444  if isfield(CellInfo{icell},'YIndex_discrete')&& ~isempty(CellInfo{icell}.YIndex_discrete)


[729]  445  charplot_0='''+''';


 446  else


 447  charplot_0='''''';


 448  end


 449  MinY=[];


 450  MaxY=[];%default


 451 


 452  nbplot=0;


 453  for ivar=1:length(VarIndex)


 454  if testplot(VarIndex(ivar))


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


 456  nbplot=nbplot+1;


 457  ytitle=[ytitle VarName];


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


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


 460  else


 461  ytitle=[ytitle ', '];


 462  end


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


 464  MinY(ivar)=min(min(data.(VarName)));


 465  MaxY(ivar)=max(max(data.(VarName)));


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


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


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


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


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


 471  end


 472  if nbcomponent1==1 nbcomponent2==1


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


 474  else %variable with severals components


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


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


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


 478  end


[201]  479  end


 480  end


[729]  481  if ~isempty(MinY)


 482  MinY_cell(icell)=min(MinY);


 483  MaxY_cell(icell)=max(MaxY);


 484  end


[201]  485  end


 486  end


 487 


 488  %% activate the plot


[644]  489  if ~isequal(plotstr,'hhh=plot(')


[201]  490  set(hfig,'CurrentAxes',haxes)


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


 492  %%%


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


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


 495  %%%


[651]  496  set(haxes,'tag',tag)% restitute the axes tag (removed by the command plot)


 497  set(haxes,'ColorOrder',ColorOrder)% restitute the plot color order (to get red green blue for histograms or cuts of color images)


[201]  498  grid(haxes, 'on')


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


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


 501  if length(legend_str)>=1


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


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


 504  end


 505  if ~isempty(legend_str)


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


 507  if isempty(hlegend)


 508  hlegend=legend(legend_str);


 509  txt=ver('MATLAB');


 510  Release=txt.Release;


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


 512  if relnumb >= 14


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


 514  end


 515  else


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


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


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


 519  end


 520  end


 521  end


 522  title_str='';


 523  if isfield(data,'filename')


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


 525  title_str=[title_str ext];


 526  end


[477]  527  if isfield(data,'Action')&&isfield(data.Action,'ActionName')


[201]  528  if ~isequal(title_str,'')


 529  title_str=[title_str ', '];


 530  end


[477]  531  title_str=[title_str data.Action.ActionName];


[201]  532  end


 533  htitle=title(title_str);


[428]  534  set(htitle,'Interpreter','none')% desable tex interpreter


[201]  535  end


 536 


 537  %% determine axes bounds


[294]  538  fix_lim=isfield(Coordinates,'CheckFixLimits') && Coordinates.CheckFixLimits;


[569]  539  check_lim=isfield(Coordinates,'MinX')&&isfield(Coordinates,'MaxX')&&isfield(Coordinates,'MinY')&&isfield(Coordinates,'MaxY');


[201]  540  if fix_lim


[569]  541  if ~check_lim


 542  fix_lim=0; %free limits if limits are not set,


[221]  543  end


[201]  544  end


[294]  545  if fix_lim


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


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


[569]  548  else


 549  if ~isempty(MinX)


 550  if check_lim


[874]  551  Coordinates.MinX=min(min(MinX),Coordinates.MinX);


 552  Coordinates.MaxX=max(max(MaxX),Coordinates.MaxX);


[569]  553  else


[874]  554  Coordinates.MinX=min(MinX);


 555  Coordinates.MaxX=max(MaxX);


[569]  556  end


 557  end


 558  if ~isempty(MinY_cell)


 559  if check_lim


[874]  560  Coordinates.MinY=min(min(MinY_cell),Coordinates.MinY);


 561  Coordinates.MaxY=max(max(MaxY_cell),Coordinates.MaxY);


[569]  562  else


[874]  563  Coordinates.MinY=min(MinY_cell);


 564  Coordinates.MaxY=max(MaxY_cell);


[569]  565  end


 566  end


[201]  567  end


 568 


[428]  569  %% determine plot aspect ratio


[789]  570  if isfield(Coordinates,'CheckFixAspectRatio') && isequal(Coordinates.CheckFixAspectRatio,1)&&isfield(Coordinates,'AspectRatio')


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


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


 573  else


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


 575  AspectRatio=get(haxes,'DataAspectRatio');


[874]  576  Coordinates.AspectRatio=AspectRatio(1)/AspectRatio(2);


[789]  577  end


[874]  578  PlotParamOut.Axes= Coordinates;


[428]  579 


[871]  580  %% give statistics for pdf


[880]  581  %ind_var=find(testplot);


[881]  582  TableData={'Variable';'SampleNbr';'bin size';'Mean';'RMS';'Skewness';'Kurtosis';...


[873]  583  'Min';'FirstCentile';'FirstDecile';'Median';'LastDecile';'LastCentile';'Max'};


[880]  584 


[874]  585  TextDisplay=0;


[873]  586  for icell=1:numel(CellInfo)


[1058]  587  if isfield(CellInfo{icell},'VarIndex_histo')% case of histogram plot


[874]  588  TextDisplay=1;


[873]  589  VarName=data.ListVarName{CellInfo{icell}.CoordIndex};


 590  pdf_val=data.(data.ListVarName{CellInfo{icell}.VarIndex_histo});


 591  x=coord_x{icell};


[881]  592  if isrow(x)


 593  x=x';


 594  end


 595  if ~isequal(size(x,1),size(pdf_val,1))


 596  pdf_val=pdf_val';


 597  end


 598  Val=pdf2stat(x,pdf_val);


 599  Column=mat2cell(Val,ones(13,1),ones(1,size(Val,2)));


 600  if size(Val,2)==1%single component


 601  TitleBar={VarName};


 602  else


 603  TitleBar=cell(1,size(Val,2));


 604  for icomp=1:size(Val,2)


 605  TitleBar{icomp}=[VarName '_' num2str(icomp)];


 606  end


 607  end


 608  Column=[TitleBar;Column];


[873]  609  TableData=[TableData Column];


 610  end


 611  end


[874]  612  if TextDisplay;


[873]  613  disp(TableData);


[874]  614  PlotParamOut.TableDisplay=TableData;


 615  else


 616  if isfield(PlotParamOut,'TableDisplay')


 617  PlotParamOut=rmfield(PlotParamOut,'TableDisplay');


 618  end


[871]  619  end


 620 


[201]  621  %


[781]  622  function [haxes,PlotParamOut,PlotType,errormsg]=plot_plane(Data,CellInfo,haxes,PlotParam)


[201]  623  %


[674]  624  PlotType='plane';


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


[411]  626 


[201]  627  %default plotting parameters


 628  if ~isfield(PlotParam,'Scalar')


 629  PlotParam.Scalar=[];


 630  end


 631  if ~isfield(PlotParam,'Vectors')


 632  PlotParam.Vectors=[];


 633  end


[674]  634  PlotParamOut=PlotParam;%default


 635  errormsg='';%default


[201]  636 


[674]  637  hfig=get(haxes,'parent');%handle of the figure containing the plot axes


[781]  638  PosColorbar=[];


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


 640  if isfield(FigData,'PosColorbar')


 641  PosColorbar=FigData.PosColorbar;


 642  end


[201]  643  hcol=findobj(hfig,'Tag','Colorbar'); %look for colorbar axes


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


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


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


 647  test_black=0;


 648  test_false=0;


 649  test_C=0;


 650  XName='';


 651  x_units='';


 652  YName='';


 653  y_units='';


[674]  654 


 655  % loop on the input field cells


[1031]  656  for icell=1:numel(CellInfo)


[530]  657  if strcmp(CellInfo{icell}.CoordType,'tps') %do not plot directly tps data (used for projection only)


[402]  658  continue


 659  end


[530]  660  ivar_X=CellInfo{icell}.CoordIndex(end); % defines (unique) index for the variable representing unstructured x coordinate (default =[])


 661  ivar_Y=CellInfo{icell}.CoordIndex(end1); % defines (unique)index for the variable representing unstructured y coordinate (default =[])


 662  ivar_C=[];


 663  if isfield(CellInfo{icell},'VarIndex_scalar')


 664  ivar_C=[ivar_C CellInfo{icell}.VarIndex_scalar];


 665  end


 666  if isfield(CellInfo{icell},'VarIndex_image')


 667  ivar_C=[ivar_C CellInfo{icell}.VarIndex_image];


 668  end


 669  if isfield(CellInfo{icell},'VarIndex_color')


 670  ivar_C=[ivar_C CellInfo{icell}.VarIndex_color];


 671  end


 672  if isfield(CellInfo{icell},'VarIndex_ancillary')


 673  ivar_C=[ivar_C CellInfo{icell}.VarIndex_ancillary];


 674  end


[201]  675  if numel(ivar_C)>1


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


 677  return


 678  end


[530]  679  ivar_F=[];


 680  if isfield(CellInfo{icell},'VarIndex_warnflag')


[674]  681  ivar_F=CellInfo{icell}.VarIndex_warnflag; %defines index (unique) for warning flag variable


[201]  682  end


[674]  683  ivar_FF_vec=[];


[530]  684  if isfield(CellInfo{icell},'VarIndex_vector_x')&&isfield(CellInfo{icell},'VarIndex_vector_y') % vector components detected


[674]  685  if test_vec% a vector field has been already detected


[866]  686  errormsg='error in plot_field: attempt to plot two vector fields: to get the difference project on a plane with ProjMode= interp_lin or interp_tps';


[201]  687  return


 688  else


[864]  689  if numel(CellInfo{icell}.VarIndex_vector_x)>1


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


[1031]  691  return


[864]  692  end


[201]  693  test_vec=1;


[674]  694  if isfield(CellInfo{icell},'VarIndex_errorflag')


 695  ivar_FF_vec=CellInfo{icell}.VarIndex_errorflag; %defines index (unique) for error flag variable


 696  end


 697  vec_U=Data.(Data.ListVarName{CellInfo{icell}.VarIndex_vector_x});


[530]  698  vec_V=Data.(Data.ListVarName{CellInfo{icell}.VarIndex_vector_y});


[1001]  699  XName=Data.ListVarName{CellInfo{icell}.CoordIndex(end)};


 700  YName=Data.ListVarName{CellInfo{icell}.CoordIndex(end1)};


[674]  701  if strcmp(CellInfo{icell}.CoordType,'scattered')%2D field with unstructured coordinates


[517]  702  vec_X=reshape(Data.(XName),[],1); %transform vectors in column matlab vectors


 703  vec_Y=reshape(Data.(YName),[],1);


[674]  704  elseif strcmp(CellInfo{icell}.CoordType,'grid')%2D field with structured coordinates


[1001]  705  y=Data.(YName);


 706  x=Data.(XName);


[201]  707  if numel(y)==2 % y defined by first and last values on aregular mesh


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


 709  end


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


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


 712  end


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


[201]  714  end


[517]  715  if isfield(PlotParam.Vectors,'ColorScalar') && ~isempty(PlotParam.Vectors.ColorScalar)


 716  [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp([],Data,PlotParam.Vectors.ColorScalar);


[530]  717  if ~isempty(VarVal)


 718  vec_C=reshape(VarVal{1},1,numel(VarVal{1}));


 719  test_C=1;


 720  end


[201]  721  end


[674]  722  if ~isempty(ivar_F)%~(isfield(PlotParam.Vectors,'HideWarning')&& isequal(PlotParam.Vectors.HideWarning,1))


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


 724  if ~(isfield(PlotParam.Vectors,'CheckHideWarning') && isequal(PlotParam.Vectors.CheckHideWarning,1))


 725  test_black=1;


 726  end


[201]  727  end


[674]  728  if ~isempty(ivar_FF_vec) %&& ~test_false


[1031]  729  vec_FF=Data.(Data.ListVarName{ivar_FF_vec}); % flags for false vectors


[201]  730  end


 731  end


 732  elseif ~isempty(ivar_C) %scalar or image


 733  if test_ima


[674]  734  errormsg='attempt to plot two scalar fields or images';


[201]  735  return


 736  end


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


[201]  738  test_ima=1;


[674]  739  if strcmp(CellInfo{icell}.CoordType,'scattered')%2D field with unstructured coordinates


[227]  740  A=reshape(A,1,[]);


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


 742  YName=Data.ListVarName{ivar_Y};


[1001]  743  Coord_x=reshape(Data.(XName),1,[]);


 744  Coord_y=reshape(Data.(YName),1,[]);


[782]  745  [A,Coord_x,Coord_y]=proj_grid(Coord_x',Coord_y',A',[],[],'np>256'); % interpolate on a grid


[201]  746  if isfield(Data,'VarAttribute')


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


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


 749  end


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


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


 752  end


[674]  753  end


 754  elseif strcmp(CellInfo{icell}.CoordType,'grid')%2D field with structured coordinates


[530]  755  YName=Data.ListVarName{CellInfo{icell}.CoordIndex(end1)};


[1001]  756  XName=Data.ListVarName{CellInfo{icell}.CoordIndex(end)};


[782]  757  Coord_y=Data.(YName);


[1001]  758  Coord_x=Data.(XName);


[201]  759  test_interp_X=0; %default, regularly meshed X coordinate


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


 761  if isfield(Data,'VarAttribute')


[530]  762  if numel(Data.VarAttribute)>=CellInfo{icell}.CoordIndex(end) && isfield(Data.VarAttribute{CellInfo{icell}.CoordIndex(end)},'units')


 763  x_units=Data.VarAttribute{CellInfo{icell}.CoordIndex(end)}.units;


[201]  764  end


[530]  765  if numel(Data.VarAttribute)>=CellInfo{icell}.CoordIndex(end1) && isfield(Data.VarAttribute{CellInfo{icell}.CoordIndex(end1)},'units')


 766  y_units=Data.VarAttribute{CellInfo{icell}.CoordIndex(end1)}.units;


[201]  767  end


[674]  768  end


[782]  769  if numel(Coord_y)>2


 770  DCoord_y=diff(Coord_y);


 771  DCoord_y_min=min(DCoord_y);


 772  DCoord_y_max=max(DCoord_y);


 773  if sign(DCoord_y_min)~=sign(DCoord_y_max);% =1 for increasing values, 0 otherwise


[1001]  774  errormsg=['errror in plot_field.m: non monotonic dimension variable ' YName ];


[674]  775  return


 776  end


[782]  777  test_interp_Y=(DCoord_y_maxDCoord_y_min)> 0.0001*abs(DCoord_y_max);


[201]  778  end


[782]  779  if numel(Coord_x)>2


 780  DCoord_x=diff(Coord_x);


 781  DCoord_x_min=min(DCoord_x);


 782  DCoord_x_max=max(DCoord_x);


 783  if sign(DCoord_x_min)~=sign(DCoord_x_max);% =1 for increasing values, 0 otherwise


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


 785  return


 786  end


[782]  787  test_interp_X=(DCoord_x_maxDCoord_x_min)> 0.0001*abs(DCoord_x_max);


[674]  788  end


 789  if test_interp_Y


[782]  790  npxy(1)=max([256 floor((Coord_y(end)Coord_y(1))/DCoord_y_min) floor((Coord_y(end)Coord_y(1))/DCoord_y_max)]);


 791  yI=linspace(Coord_y(1),Coord_y(end),npxy(1));


[201]  792  if ~test_interp_X


[782]  793  xI=linspace(Coord_x(1),Coord_x(end),size(A,2));%default


 794  Coord_x=xI;


[201]  795  end


 796  end


[674]  797  if test_interp_X


[782]  798  npxy(2)=max([256 floor((Coord_x(end)Coord_x(1))/DCoord_x_min) floor((Coord_x(end)Coord_x(1))/DCoord_x_max)]);


 799  xI=linspace(Coord_x(1),Coord_x(end),npxy(2));


[201]  800  if ~test_interp_Y


[782]  801  yI=linspace(Coord_y(1),Coord_y(end),size(A,1));


 802  Coord_y=yI;


[201]  803  end


 804  end


[674]  805  if test_interp_X  test_interp_Y


[782]  806  [Coord_x2D,Coord_y2D]=meshgrid(Coord_x,Coord_y);


 807  A=interp2(Coord_x2D,Coord_y2D,double(A),xI,yI');


[201]  808  end


[782]  809  Coord_x=[Coord_x(1) Coord_x(end)];% keep only the lower and upper bounds for image represnetation


 810  Coord_y=[Coord_y(1) Coord_y(end)];


[201]  811  end


 812  end


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


 814  if isfield(Data,'CoordUnit')


 815  if isempty(x_units)


 816  x_units=Data.CoordUnit;


 817  end


 818  if isempty(y_units)


 819  y_units=Data.CoordUnit;


 820  end


[1031]  821  end


[674]  822  end


[1031]  823 


[688]  824  PlotParamOut=PlotParam; % output plot parameters equal to input by default


[201]  825 


 826  %% image or scalar plot %%%%%%%%%%%%%%%%%%%%%%%%%%


 827  if test_ima


 828  % distinguish B/W and color images


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


 830  siz=numel(np);


 831  if siz>3


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


 833  return


[201]  834  end


 835  if siz==3


 836  if np(3)==1


 837  siz=2;%B W image


 838  elseif np(3)==3


 839  siz=3;%color image


 840  else


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


[201]  842  return


 843  end


 844  end


 845 


[822]  846  %set for grey scale setting


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


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


 849  else % BW imposed automatically chosen


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


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


[428]  852  end


[688]  853 


[822]  854  % determine the plot option 'image' or 'contours'


 855  CheckContour=0; %default


[688]  856  if isfield(PlotParam.Scalar,'ListContour')


 857  CheckContour=strcmp(PlotParam.Scalar.ListContour,'contours');% =1 for contour plot option


[822]  858  end


 859 


 860  %case of grey level images or contour plot


 861  if ~isfield(PlotParam.Scalar,'CheckFixScalar')


 862  PlotParam.Scalar.CheckFixScalar=0;% free scalar threshold value scale (from min to max) by default


 863  end


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


 865  PlotParam.Scalar.MinA=[];%no min scalar threshold value set


 866  end


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


 868  PlotParam.Scalar.MaxA=[];%no max scalar threshold value set


 869  end


 870 


 871  % determine the min scalar value


 872  if PlotParam.Scalar.CheckFixScalar && ~isempty(PlotParam.Scalar.MinA) && isnumeric(PlotParam.Scalar.MinA)


 873  MinA=double(PlotParam.Scalar.MinA); % min value set as input


[688]  874  else


[822]  875  MinA=double(min(min(min(A)))); % min value set as min of non NaN scalar values


[688]  876  end


 877 


[822]  878  % error if the input scalar is NaN everywhere


 879  if isnan(MinA)


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


 881  return


 882  end


 883 


 884  % determine the max scalar value


 885  CheckFixScalar=0;


 886  if PlotParam.Scalar.CheckFixScalar && ~isempty(PlotParam.Scalar.MaxA) && isnumeric(PlotParam.Scalar.MaxA)


 887  MaxA=double(PlotParam.Scalar.MaxA); % max value set as input


 888  CheckFixScalar=1;


 889  else


 890  MaxA=double(max(max(max(A)))); % max value set as min of non NaN scalar values


 891  end


 892 


 893  PlotParamOut.Scalar.MinA=MinA;


 894  PlotParamOut.Scalar.MaxA=MaxA;


 895  PlotParamOut.Scalar.Npx=size(A,2);


 896  PlotParamOut.Scalar.Npy=size(A,1);


 897  % if siz==2


 898  % case of contour plot


 899  if CheckContour


 900  if ~isempty(hima) && ishandle(hima)


 901  delete(hima) % delete existing image


[201]  902  end


[822]  903 


 904  % set the contour values


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


 906  PlotParam.Scalar.IncrA=[];% automatic contour interval


[201]  907  end


[822]  908  if ~isempty(PlotParam.Scalar.IncrA) && isnumeric(PlotParam.Scalar.IncrA)


 909  interval=PlotParam.Scalar.IncrA;


 910  else % automatic contour interval


 911  cont=colbartick(MinA,MaxA);


 912  interval=cont(2)cont(1);%default


 913  PlotParamOut.Scalar.IncrA=interval;% set the interval as output for display on the GUI


[201]  914  end


[822]  915  abscontmin=interval*floor(MinA/interval);


 916  abscontmax=interval*ceil(MaxA/interval);


 917  contmin=interval*floor(min(min(A))/interval);


 918  contmax=interval*ceil(max(max(A))/interval);


 919  cont_pos_plus=0:interval:contmax;% zero and positive contour values (plotted as solid lines)


 920  cont_pos_min=double(contmin):interval:interval;% negative contour values (plotted as dashed lines)


 921  cont_pos=[cont_pos_min cont_pos_plus];% set of all contour values


[688]  922 


[822]  923  sizpx=(Coord_x(end)Coord_x(1))/(np(2)1);


 924  sizpy=(Coord_y(1)Coord_y(end))/(np(1)1);


 925  x_cont=Coord_x(1):sizpx:Coord_x(end); % pixel x coordinates for image display


 926  y_cont=Coord_y(1):sizpy:Coord_y(end); % pixel x coordinates for image display


[688]  927 


[822]  928  %axes(haxes)% set the input axes handle as current axis


 929 


 930  % colormap(map);


 931  tag_axes=get(haxes,'Tag');% axes tag


 932  Opacity=1;


 933  if isfield(PlotParam.Scalar,'Opacity')&&~isempty(PlotParam.Scalar.Opacity)


 934  Opacity=PlotParam.Scalar.Opacity;


[688]  935  end


[822]  936  % fill the space between contours if opacity is undefined or =1


 937  if isequal(Opacity,1)


 938  [var,hcontour]=contour(haxes,x_cont,y_cont,A,cont_pos);% determine all contours


 939  set(hcontour,'Fill','on')% fill the space between contours


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


 941  hold on


 942  end


 943  [var_p,hcontour_p]=contour(haxes,x_cont,y_cont,A,cont_pos_plus,'k');% draw the contours for positive values as solid lines


 944  hold on


 945  [var_m,hcontour_m]=contour(haxes,x_cont,y_cont,A,cont_pos_min,'');% draw the contours for negative values as dashed lines


 946  if isequal(Opacity,1)


 947  set(hcontour_m,'LineColor',[1 1 1])% draw negative contours in white (better visibility in dark background)


 948  end


 949  set(haxes,'Tag',tag_axes);% restore axes tag (removed by the matlab fct contour !)


 950  hold off


[688]  951 


[822]  952  %determine the color scale and map


 953  caxis([abscontmin abscontmax])


 954  if BW


 955  vec=linspace(0,1,(abscontmaxabscontmin)/interval);%define a greyscale colormap with steps interval


 956  map=[vec' vec' vec'];


 957  colormap(map);


[201]  958  else


[943]  959  colormap('jet'); % default matlab colormap ('jet')


[822]  960  end


[688]  961 


[822]  962  if isfield(PlotParam.Axes,'CheckFixAspectRatio') && isequal(PlotParam.Axes.CheckFixAspectRatio,1)


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


 964  if isfield(PlotParam.Axes,'AspectRatio')


 965  set(haxes,'DataAspectRatio',[PlotParam.Axes.AspectRatio 1 1])


[688]  966  else


[822]  967  set(haxes,'DataAspectRatio',[1 1 1])


[201]  968  end


[822]  969  end


 970  else %usual images (no contour)


[201]  971  % set colormap for image display


[822]  972  if BW


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


 974  map=[vec' vec' vec'];


 975  colormap(map); %grey scale color map


 976  if siz==3% true color images visualized in BW


 977  A=uint16(sum(A,3));%sum the three color components for color images displayed with BW option


[201]  978  end


[822]  979  else


 980  if siz==3 && CheckFixScalar % true color images rescaled by MaxA


 981  A=uint8(255*double(A)/double(MaxA));


 982  end


 983  colormap('default'); % standard false colors for div, vort , scalar fields


[201]  984  end


 985 


 986  % interpolate field to increase resolution of image display


[652]  987  test_interp=0;


[822]  988  if size(A,3)==1 % scalar or B/W image


[652]  989  test_interp=1;


 990  if max(np) <= 64


 991  npxy=8*np;% increase the resolution 8 times


 992  elseif max(np) <= 128


 993  npxy=4*np;% increase the resolution 4 times


 994  elseif max(np) <= 256


 995  npxy=2*np;% increase the resolution 2 times


 996  else


 997  npxy=np;


 998  test_interp=0; % no interpolation done


 999  end


[201]  1000  end


[652]  1001  if test_interp%if we interpolate


[782]  1002  x=linspace(Coord_x(1),Coord_x(2),np(2));


 1003  y=linspace(Coord_y(1),Coord_y(2),np(1));


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


[782]  1005  xi=linspace(Coord_x(1),Coord_x(2),npxy(2));


 1006  yi=linspace(Coord_y(1),Coord_y(2),npxy(1));


[688]  1007  A = interp2(X,Y,double(A),xi,yi');


[428]  1008  end


[848]  1009  % create new image if no image handle is found


[428]  1010  if isempty(hima)


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


 1012  if MinA<MaxA


[782]  1013  hima=imagesc(Coord_x,Coord_y,A,[MinA MaxA]);


[201]  1014  else % to deal with uniform field


[782]  1015  hima=imagesc(Coord_x,Coord_y,A,[MaxA1 MaxA]);


[201]  1016  end


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


[411]  1018  % requested:


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


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


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


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


[428]  1023  % update an existing image


[201]  1024  else


[688]  1025  set(hima,'CData',A);


[201]  1026  if MinA<MaxA


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


 1028  else


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


[201]  1030  end


[782]  1031  set(hima,'XData',Coord_x);


 1032  set(hima,'YData',Coord_y);


[201]  1033  end


[546]  1034 


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


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


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


[405]  1038  else


[428]  1039  if test_vec


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


 1041  PlotParamOut.Scalar.Opacity=0.5;


 1042  else


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


 1044  end


[405]  1045  end


[201]  1046  end


 1047  test_ima=1;


 1048 


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


[546]  1050  if ~isempty(PosColorbar)


[822]  1051  if size(A,3)==1 && exist('PosColorbar','var')


[546]  1052  if isempty(hcol)~ishandle(hcol)


 1053  hcol=colorbar;%create new colorbar


 1054  end


 1055  if length(PosColorbar)==4


 1056  set(hcol,'Position',PosColorbar)


 1057  end


 1058  %YTick=0;%default


 1059  if MaxA>MinA


 1060  if CheckContour


 1061  colbarlim=get(hcol,'YLim');


 1062  scale_bar=(colbarlim(2)colbarlim(1))/(abscontmaxabscontmin);


 1063  YTick=cont_pos(2:end1);


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


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


 1066  elseif (isfield(PlotParam.Scalar,'CheckBW') && isequal(PlotParam.Scalar.CheckBW,1))isa(A,'uint8') isa(A,'uint16')%images


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


 1068  if iscell(hi)%multiple images in colorbar


 1069  hi=hi{1};


 1070  end


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


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


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


 1074  YTick=colbartick(MinA,MaxA);


 1075  set(hcol,'YTick',YTick)


 1076  else


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


 1078  if iscell(hi)%multiple images in colorbar


 1079  hi=hi{1};


 1080  end


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


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


 1083  YTick=colbartick(MinA,MaxA);


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


 1085  set(hcol,'YTick',YTick)


[201]  1086  end


[546]  1087  set(hcol,'Yticklabel',num2str(YTick'));


[201]  1088  end


[546]  1089  elseif ishandle(hcol)


 1090  delete(hcol); %erase existing colorbar if not needed


[201]  1091  end


 1092  end


 1093  else%no scalar plot


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


[201]  1095  delete(hima)


 1096  end


[546]  1097  if ~isempty(PosColorbar) && ~isempty(hcol)&& ishandle(hcol)


[428]  1098  delete(hcol)


[201]  1099  end


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


 1101  end


 1102 


 1103  %% vector plot %%%%%%%%%%%%%%%%%%%%%%%%%%


 1104  if test_vec


 1105  %vector scale representation


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


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


 1108  end


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


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


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


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


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


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


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


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


 1117  else


 1118  if ~test_false %remove false vectors


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


 1120  end


 1121  if isempty(vec_U)


 1122  scale=1;


 1123  else


 1124  if isempty(indsel)


 1125  MaxU=max(abs(vec_U));


 1126  MaxV=max(abs(vec_V));


 1127  else


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


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


 1130  end


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


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


 1133  end


[210]  1134  end


[201]  1135 


 1136  %record vectors on the plotting axes


 1137  if test_C==0


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


 1139  end


 1140 


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


[581]  1142  check_decimate=0;


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


[581]  1144  check_decimate=1;


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


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


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


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


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


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


 1151  end


 1152  nb_sel=length(ind_sel);


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


[581]  1154  elseif isfield(PlotParam.Vectors,'CheckDecimate16') && PlotParam.Vectors.CheckDecimate16


 1155  check_decimate=1;


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


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


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


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


 1160  for i=2:4:length(ind_jump)1


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


 1162  end


 1163  nb_sel=length(ind_sel);


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


 1165  end


 1166  if check_decimate


[201]  1167  vec_X=vec_X(ind_sel);


 1168  vec_Y=vec_Y(ind_sel);


 1169  vec_U=vec_U(ind_sel);


 1170  vec_V=vec_V(ind_sel);


 1171  vec_C=vec_C(ind_sel);


 1172  if ~isempty(ivar_F)


 1173  vec_F=vec_F(ind_sel);


 1174  end


[674]  1175  if ~isempty(ivar_FF_vec)


[201]  1176  vec_FF=vec_FF(ind_sel);


 1177  end


 1178  end


 1179 


 1180  %get main level color code


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


 1182 


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


 1184  sizlist=size(colorlist);


 1185  nbcolor=sizlist(1);


 1186  if test_black


 1187  nbcolor=nbcolor+1;


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


[674]  1189  if ~isempty(ivar_FF_vec)


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


[201]  1191  else


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


[201]  1193  end


 1194  end


 1195  nbcolor=nbcolor+1;


[674]  1196  if ~isempty(ivar_FF_vec)


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


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


 1199  else


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


 1201  end


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


[201]  1203  end


 1204  %plot vectors:


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


 1206 


 1207  else


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


 1209  if ~isempty(hvec)


 1210  delete(hvec);


 1211  end


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


 1213  end


[546]  1214  % nbvar=0;


[201]  1215 


 1216  %store the coordinate extrema occupied by the field


 1217  if ~isempty(Data)


[569]  1218  MinX=[];


 1219  MaxX=[];


 1220  MinY=[];


 1221  MaxY=[];


[748]  1222  fix_lim=isfield(PlotParam.Axes,'CheckFixLimits') && PlotParam.Axes.CheckFixLimits;


[201]  1223  if fix_lim


[748]  1224  if isfield(PlotParam.Axes,'MinX')&&isfield(PlotParam.Axes,'MaxX')&&isfield(PlotParam.Axes,'MinY')&&isfield(PlotParam.Axes,'MaxY')


 1225  MinX=PlotParam.Axes.MinX;


 1226  MaxX=PlotParam.Axes.MaxX;


 1227  MinY=PlotParam.Axes.MinY;


 1228  MaxY=PlotParam.Axes.MaxY;


[569]  1229  end %else PlotParamOut.MinX =PlotParam.MinX...


[252]  1230  else


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


[782]  1232  MinX=min(Coord_x);


 1233  MaxX=max(Coord_x);


 1234  MinY=min(Coord_y);


 1235  MaxY=max(Coord_y);


[252]  1236  if test_vec


[569]  1237  MinX=min(MinX,min(vec_X));


 1238  MaxX=max(MaxX,max(vec_X));


 1239  MinY=min(MinY,min(vec_Y));


 1240  MaxY=max(MaxY,max(vec_Y));


[252]  1241  end


 1242  elseif test_vec


[569]  1243  MinX=min(vec_X);


 1244  MaxX=max(vec_X);


 1245  MinY=min(vec_Y);


 1246  MaxY=max(vec_Y);


[201]  1247  end


[252]  1248  end


[748]  1249  PlotParamOut.Axes.MinX=MinX;


 1250  PlotParamOut.Axes.MaxX=MaxX;


 1251  PlotParamOut.Axes.MinY=MinY;


 1252  PlotParamOut.Axes.MaxY=MaxY;


[569]  1253  if MaxX>MinX


 1254  set(haxes,'XLim',[MinX MaxX]);% set x limits of frame in axes coordinates


[546]  1255  end


[569]  1256  if MaxY>MinY


 1257  set(haxes,'YLim',[MinY MaxY]);% set x limits of frame in axes coordinates


[546]  1258  end


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


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


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


[748]  1262  PlotParamOut.Axes.x_units=x_units;


 1263  PlotParamOut.Axes.y_units=y_units;


[201]  1264  end


[760]  1265 


[201]  1266  %


 1267  %  function for plotting vectors


 1268  %INPUT:


 1269  % haxes: handles of the plotting axes


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


 1271  % scale: scaling factor for vector length representation


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


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


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


 1275  %


 1276  %define arrows


 1277  theta=0.5 ;%angle arrow


 1278  alpha=0.3 ;%length arrow


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


 1280  %find the existing lines


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


 1282  sizh=size(h);


[932]  1283  %set(h,'EraseMode','xor');


[201]  1284  set(haxes,'NextPlot','replacechildren');


 1285 


 1286  %drawnow


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


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


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


 1290  end


 1291  sizlist=size(colorlist);


 1292  ncolor=sizlist(1);


 1293 


 1294  for icolor=1:ncolor


 1295  %determine the line positions for each color icolor


 1296  ind=find(col_vec==icolor);


 1297  xc=x(ind);


 1298  yc=y(ind);


 1299  uc=u(ind)*scale;


 1300  vc=v(ind)*scale;


 1301  n=size(xc);


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


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


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


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


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


 1307 


 1308  %determine arrow heads


 1309  arrowplus=rot*[uc;vc];


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


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


 1312  x2=xc+uc/2;


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


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


 1315  y2=yc+vc/2;


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


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


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


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


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


 1321  %draw the line or modify the existing ones


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


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


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


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


 1326  if n(2)>0


 1327  hold on


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


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


 1330  end


 1331  end


 1332  else


 1333  if isn(1)


 1334  delete(h(2*icolor1))


 1335  delete(h(2*icolor))


 1336  else


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


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


[932]  1339  % set(h(2*icolor1),'EraseMode','xor');


[201]  1340  set(h(2*icolor),'Xdata',matxar,'Ydata',matyar);


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


[932]  1342  %set(h(2*icolor),'EraseMode','xor');


[201]  1343  end


 1344  end


 1345  end


 1346  if sizh(1) > 2*ncolor


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


 1348  delete(h(2*icolor1))


 1349  delete(h(2*icolor))


 1350  end


 1351  end


 1352 


 1353  %


 1354  %  determine tick positions for colorbar


 1355  function YTick=colbartick(MinA,MaxA)


 1356  %


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


 1358  YTick=0;%default


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


 1360  if maxabs>0


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


 1362  div=1;


 1363  siz2=1;


 1364  while siz2<2


 1365  values=10:div:10;


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


 1367  siz=size(ind);


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


 1369  values=9:0.5*div:9;


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


 1371  end


 1372  siz2=size(ind,2);


 1373  div=div/10;


 1374  end


 1375  YTick=ord*values(ind);


 1376  end


 1377 


[356]  1378  % 


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


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


[389]  1381  % 


[356]  1382  if length(vec_Y)<2


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


 1384  return;


 1385  end


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


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


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


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


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


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


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


 1393  rangx=rgx_in; % first and last positions


 1394  rangy=rgy_in;


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


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


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


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


 1399  else % interpolation grid automatically determined


 1400  rangx(1)=min(vec_X);


 1401  rangx(2)=max(vec_X);


 1402  rangy(2)=min(vec_Y);


 1403  rangy(1)=max(vec_Y);


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


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


 1406  end


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


 1408  yi=rangy(1):dxy(1):rangy(2);


 1409  A=griddata(vec_X,vec_Y,vec_A,xi,yi');


[356]  1410  A=reshape(A,length(yi),length(xi));


 1411  else


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


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


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


 1415  ymin=vec_Y(index(1));


 1416  y=vec_Y(index);


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


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


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


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


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


 1422 


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


 1424  %to prescribed positions


 1425  test_interp=1;


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


 1427  rangx=rgx_in; % first and last positions


 1428  rangy=rgy_in;


 1429  npxy=npxy_in;


 1430  else


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


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


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


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


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


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


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


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


 1439  else


 1440  npxy=[ny nx];


 1441  test_interp=0; % no interpolation done


 1442  end


 1443  end


 1444  if test_interp==1%if we interpolate


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


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


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


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


 1449  else %no interpolation for a resolution higher than 256


 1450  A=B;


 1451  end


[801]  1452  end

