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


[811]  16  % OUTPUT:


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


[811]  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 


[809]  83  %=======================================================================


 84  % Copyright 20082014, LEGI UMR 5519 / CNRS UJF GINP, Grenoble, France


 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.


[809]  89  %


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


[809]  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


[809]  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  [CellInfo,NbDimArray,errormsg]=find_field_cells(Data);


[388]  118  if ~isempty(errormsg)


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


[388]  120  return


 121  end


[530]  122 


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


[388]  124  if ~isempty(index_3D)


[546]  125  msgbox_uvmat('ERROR','volume plot not implemented yet');


 126  return


[388]  127  end


[530]  128  index_2D=find(NbDimArray==2);%find 2D fields


 129  index_1D=find(NbDimArray==1);


 130  index_0D=find(NbDimArray==0);


[294]  131 


[388]  132  %% test axes and figure


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


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


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


 136  if ishandle(haxes)


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


 138  testnewfig=0;


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


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


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


 142  testzoomaxes=1;


 143  zoomaxes=AxeData.ZoomAxes;


[388]  144  end


 145  end


 146  end


[201]  147  end


[429]  148  end


 149 


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


 151  if testnewfig


 152  hfig=figure;


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


 154  haxes=axes;


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


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


[429]  157  else


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


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


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


 161  end


 162 


 163  %% set axes properties


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


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


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


 167  else


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


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


 170  end


[760]  171 


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


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


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


 175  else


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


 177  end


 178 


[429]  179  errormsg='';


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


 181 


 182  %% 2D plots


 183  if isempty(index_2D)


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


[429]  185  else %plot 2D field


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


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


[429]  188  AxeData.NbDim=2;


 189  if testzoomaxes && isempty(errormsg)


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


[429]  191  AxeData.ZoomAxes=zoomaxes;


[201]  192  end


[429]  193  end


 194 


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


 196  if isempty(index_1D)


 197  if ~isempty(haxes)


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


[201]  199  end


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


[644]  201  CheckHold=0;


 202  if isfield(PlotParam,'CheckHold')


 203  CheckHold= PlotParam.CheckHold;


 204  end


[748]  205  PlotParamOut.Axes=plot_profile(Data,CellInfo(index_1D),haxes,PlotParamOut.Axes,CheckHold);%


[429]  206  if testzoomaxes


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


[429]  208  AxeData.ZoomAxes=zoomaxes;


[201]  209  end


[429]  210  PlotType='line';


 211  end


 212 


[760]  213  %% aspect ratio


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


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


 216 


[429]  217  %% text display


[690]  218  htext=findobj(hfig,'Tag','TableDisplay');


 219  hchecktable=findobj(hfig,'Tag','CheckTable');


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


 221  % htext=findobj(hfig,'Tag','TableDisplay');


 222  % else %text display added to plot


 223  % %htext=findobj(hfig,'Tag','text_display');


 224  % end


 225  if ~isempty(htext)&&~isempty(hchecktable)


[429]  226  if isempty(index_0D)


[690]  227  set(htext,'Data',{})


 228  set(htext,'visible','off')


 229  set(hchecktable,'visible','off')


 230  set(hchecktable,'Value',0)


 231  else


 232  errormsg=plot_text(Data,CellInfo(index_0D),htext);


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


 234  set(hchecktable,'visible','on')


 235  set(hchecktable,'Value',1)


 236  end


 237  set(hfig,'Unit','pixels');


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


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


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


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


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


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


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


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


[434]  246  else


[690]  247  set(haxes,'Visible','on')


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


[434]  249  end


[294]  250  end


[201]  251  end


 252 


[388]  253  %% display error message


[201]  254  if ~isempty(errormsg)


 255  msgbox_uvmat('ERROR', errormsg)


 256  end


 257 


 258  %% update the parameters stored in AxeData


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


 260  % AxeData=[];


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


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


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


[388]  264  end


 265  set(haxes,'UserData',AxeData)


 266  end


[201]  267 


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


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


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


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


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


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


 274  end


[388]  275  end


[292]  276 


[201]  277  %


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


 279  %


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


[690]  281 


[581]  282  errormsg='';


[201]  283  txt_cell={};


[429]  284  Data={};


[690]  285  VarIndex=[];


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


[535]  287 


 288  % select types of variables to be projected


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


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


 291  for ilist=1:numel(ListProj)


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


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


 294  end


 295  end


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


 297  end


 298 


 299  % data need to be displayed in a table


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


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


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


 303  for ivar=1:numel(VarIndex)


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


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


[535]  306  end


[690]  307  set(htext,'ColumnName',VarNameCell)


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


 309 


 310  for ivar=1:numel(VarIndex)


 311  VarValue=FieldData.(VarNameCell{ivar});


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


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


[201]  314  end


[690]  315  set(htext,'Data',Data)


[201]  316  end


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


 318  % for ind=1:numel(VarValue)


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


 320  % end


 321  % else


 322  % VarNameCell={VarName};


 323  % end


 324  % if numel(VarValue)<10


 325  % if isempty(VarValue)


 326  % VarValueCell={'[]'};


 327  % else


 328  % VarValueCell=num2cell(VarValue);


 329  % end


 330  % if isempty(Data)


 331  % Data =[VarNameCell VarValueCell];


 332  % else


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


 334  % end


 335  % else


 336  % if isempty(Data)


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


 338  % else


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


 340  % end


 341  % end


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


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


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


 345  % end


 346  % end


 347  % end


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


 349  %


 350  %


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


 352  % else % display in a text edit box


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


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


 355  % end


[201]  356 


[429]  357 


[201]  358  %


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


 360  %


[644]  361  function CoordinatesOut=plot_profile(data,CellInfo,haxes,Coordinates,CheckHold)


[201]  362 


[644]  363  %% initialization


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


 365  Coordinates=[];


[201]  366  end


[292]  367  CoordinatesOut=Coordinates; %default


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


[644]  369  legend_str={};


 370 


 371  %% suppress existing plot if empty data


[201]  372  if isempty(data)


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


 374  if ~isempty(hplot)


 375  delete(hplot)


 376  end


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


 378  if ~isempty(hlegend)


 379  delete(hlegend)


 380  end


 381  return


 382  end


 383 


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


[651]  385  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];


 386  set(hfig,'DefaultAxesColorOrder',ColorOrder)


[644]  387  % if isfield(Coordinates,'NextPlot')


 388  % set(haxes,'NextPlot',Coordinates.NextPlot)


 389  % end


 390  if CheckHold


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


[644]  392  else


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


[201]  394  end


 395 


 396  %% prepare the string for plot command


 397  plotstr='hhh=plot(';


[644]  398  % coord_x_index=[];


[201]  399  xtitle='';


 400  ytitle='';


[644]  401  test_newplot=~CheckHold;


[569]  402  MinX=[];


 403  MaxX=[];


 404  MinY_cell=[];


 405  MaxY_cell=[];


[201]  406  %loop on input fields


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


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


 409  coord_x_index=CellInfo{icell}.CoordIndex;


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


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


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


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


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


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


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


 417  else


 418  xtitle=[xtitle ', '];


 419  end


[201]  420  end


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


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


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


 424  testplot(coord_x_index)=0;


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


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


 427  end


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


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


 430  end


 431  if isfield(data,'VarAttribute')


 432  VarAttribute=data.VarAttribute;


 433  for ivar=1:length(VarIndex)


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


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


 436  else


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


 438  end


[201]  439  end


 440  end


[729]  441  if isfield(CellInfo{icell},'VarIndex_discrete')


 442  charplot_0='''+''';


 443  else


 444  charplot_0='''''';


 445  end


 446  MinY=[];


 447  MaxY=[];%default


 448 


 449  nbplot=0;


 450  for ivar=1:length(VarIndex)


 451  if testplot(VarIndex(ivar))


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


 453  nbplot=nbplot+1;


 454  ytitle=[ytitle VarName];


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


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


 457  else


 458  ytitle=[ytitle ', '];


 459  end


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


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


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


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


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


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


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


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


 468  end


 469  if nbcomponent1==1 nbcomponent2==1


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


 471  else %variable with severals components


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


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


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


 475  end


[201]  476  end


 477  end


[729]  478  if ~isempty(MinY)


 479  MinY_cell(icell)=min(MinY);


 480  MaxY_cell(icell)=max(MaxY);


 481  end


[201]  482  end


 483  end


 484 


 485  %% activate the plot


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


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


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


 489  %%%


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


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


 492  %%%


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


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


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


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


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


 498  if length(legend_str)>=1


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


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


 501  end


 502  if ~isempty(legend_str)


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


 504  if isempty(hlegend)


 505  hlegend=legend(legend_str);


 506  txt=ver('MATLAB');


 507  Release=txt.Release;


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


 509  if relnumb >= 14


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


 511  end


 512  else


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


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


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


 516  end


 517  end


 518  end


 519  title_str='';


 520  if isfield(data,'filename')


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


 522  title_str=[title_str ext];


 523  end


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


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


 526  title_str=[title_str ', '];


 527  end


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


[201]  529  end


 530  htitle=title(title_str);


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


[201]  532  end


 533 


 534  %% determine axes bounds


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


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


[201]  537  if fix_lim


[569]  538  if ~check_lim


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


[221]  540  end


[201]  541  end


[294]  542  if fix_lim


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


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


[569]  545  else


 546  if ~isempty(MinX)


 547  if check_lim


 548  CoordinatesOut.MinX=min(min(MinX),CoordinatesOut.MinX);


 549  CoordinatesOut.MaxX=max(max(MaxX),CoordinatesOut.MaxX);


 550  else


 551  CoordinatesOut.MinX=min(MinX);


 552  CoordinatesOut.MaxX=max(MaxX);


 553  end


 554  end


 555  if ~isempty(MinY_cell)


 556  if check_lim


 557  CoordinatesOut.MinY=min(min(MinY_cell),CoordinatesOut.MinY);


 558  CoordinatesOut.MaxY=max(max(MaxY_cell),CoordinatesOut.MaxY);


 559  else


 560  CoordinatesOut.MinY=min(MinY_cell);


 561  CoordinatesOut.MaxY=max(MaxY_cell);


 562  end


 563  end


[201]  564  end


 565 


[428]  566  %% determine plot aspect ratio


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


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


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


 570  else


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


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


 573  CoordinatesOut.AspectRatio=AspectRatio(1)/AspectRatio(2);


 574  end


[428]  575 


[201]  576  %


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


[201]  578  %


[674]  579  PlotType='plane';


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


[411]  581 


[201]  582  %default plotting parameters


 583  if ~isfield(PlotParam,'Scalar')


 584  PlotParam.Scalar=[];


 585  end


 586  if ~isfield(PlotParam,'Vectors')


 587  PlotParam.Vectors=[];


 588  end


[674]  589  PlotParamOut=PlotParam;%default


 590  errormsg='';%default


[201]  591 


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


[781]  593  PosColorbar=[];


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


 595  if isfield(FigData,'PosColorbar')


 596  PosColorbar=FigData.PosColorbar;


 597  end


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


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


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


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


 602  test_black=0;


 603  test_false=0;


 604  test_C=0;


 605  XName='';


 606  x_units='';


 607  YName='';


 608  y_units='';


[674]  609 


 610  % loop on the input field cells


 611  for icell=1:numel(CellInfo)


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


[402]  613  continue


 614  end


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


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


 617  ivar_C=[];


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


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


 620  end


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


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


 623  end


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


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


 626  end


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


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


 629  end


[201]  630  if numel(ivar_C)>1


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


 632  return


 633  end


[530]  634  ivar_F=[];


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


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


[201]  637  end


[674]  638  ivar_FF_vec=[];


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


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


[515]  641  errormsg='error in plot_field: attempt to plot two vector fields: to get the difference project on a plane with mode interp';


[201]  642  return


 643  else


 644  test_vec=1;


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


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


 647  end


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


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


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


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


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


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


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


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


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


 657  x=Data.(Data.ListVarName{CellInfo{icell}.CoordIndex(end)});


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


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


 660  end


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


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


 663  end


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


[201]  665  end


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


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


[530]  668  if ~isempty(VarVal)


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


 670  test_C=1;


 671  end


[201]  672  end


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


 674  % if test_vec


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


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


[201]  677  test_black=1;


 678  end


[674]  679  % end


[201]  680  end


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


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


 683  vec_FF=Data.(Data.ListVarName{ivar_FF_vec}); % flags for false vectors


 684  % end


[201]  685  end


 686  end


 687  elseif ~isempty(ivar_C) %scalar or image


 688  if test_ima


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


[201]  690  return


 691  end


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


[201]  693  test_ima=1;


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


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


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


 697  YName=Data.ListVarName{ivar_Y};


[782]  698  eval(['Coord_x=reshape(Data.' XName ',1,[]);'])


 699  eval(['Coord_y=reshape(Data.' YName ',1,[]);'])


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


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


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


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


 704  end


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


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


 707  end


[674]  708  end


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


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


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


 712  Coord_x=Data.(Data.ListVarName{CellInfo{icell}.CoordIndex(end)});


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


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


 715  if isfield(Data,'VarAttribute')


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


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


[201]  718  end


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


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


[201]  721  end


[674]  722  end


[782]  723  if numel(Coord_y)>2


 724  DCoord_y=diff(Coord_y);


 725  DCoord_y_min=min(DCoord_y);


 726  DCoord_y_max=max(DCoord_y);


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


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


 729  return


 730  end


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


[201]  732  end


[782]  733  if numel(Coord_x)>2


 734  DCoord_x=diff(Coord_x);


 735  DCoord_x_min=min(DCoord_x);


 736  DCoord_x_max=max(DCoord_x);


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


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


 739  return


 740  end


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


[674]  742  end


 743  if test_interp_Y


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


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


[201]  746  if ~test_interp_X


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


 748  Coord_x=xI;


[201]  749  end


 750  end


[674]  751  if test_interp_X


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


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


[201]  754  if ~test_interp_Y


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


 756  Coord_y=yI;


[201]  757  end


 758  end


[674]  759  if test_interp_X  test_interp_Y


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


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


[201]  762  end


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


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


[201]  765  end


 766  end


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


 768  if isfield(Data,'CoordUnit')


 769  if isempty(x_units)


 770  x_units=Data.CoordUnit;


 771  end


 772  if isempty(y_units)


 773  y_units=Data.CoordUnit;


 774  end


[674]  775  end


 776  end


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


[201]  778 


 779  %% image or scalar plot %%%%%%%%%%%%%%%%%%%%%%%%%%


 780  if test_ima


[688]  781 


[201]  782  % distinguish B/W and color images


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


 784  siz=numel(np);


 785  if siz>3


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


 787  return


[201]  788  end


 789  if siz==3


 790  if np(3)==1


 791  siz=2;%B W image


 792  elseif np(3)==3


 793  siz=3;%color image


 794  else


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


[201]  796  return


 797  end


 798  end


 799 


 800  %set the color map


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


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


 803  else % BW imposed automatically chosen


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


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


[428]  806  end


[688]  807 


 808  % determine the plot option 'image' or 'contours'


 809  if isfield(PlotParam.Scalar,'ListContour')


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


 811  else


 812  CheckContour=0; %default


 813  end


 814 


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


[428]  816  if siz==2


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


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


[201]  819  end


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


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


[201]  822  end


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


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


[201]  825  end


[688]  826 


 827  % determine the min scalar value


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


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


 830  else


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


[688]  832  end


 833 


 834  % error if the input scalar is NaN everywhere


 835  if isnan(MinA)


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


[206]  837  return


[688]  838  end


 839 


 840  % determine the max scalar value


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


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


[201]  843  else


[735]  844  MaxA=double(max(max(A))); % max value set as min of non NaN scalar values


[688]  845  end


 846 


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


 848  PlotParamOut.Scalar.MaxA=MaxA;


[569]  849  PlotParamOut.Scalar.Npx=size(A,2);


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


[688]  851 


[201]  852  % case of contour plot


[294]  853  if CheckContour


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


[688]  855  delete(hima) % delete existing image


[201]  856  end


[688]  857 


 858  % set the contour values


[201]  859  if ~isfield(PlotParam.Scalar,'IncrA')


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


[201]  861  end


[688]  862  if ~isempty(PlotParam.Scalar.IncrA) && isnumeric(PlotParam.Scalar.IncrA)


 863  interval=PlotParam.Scalar.IncrA;


 864  else % automatic contour interval


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


[688]  866  interval=cont(2)cont(1);%default


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


[201]  868  end


[688]  869  %B=A;


 870  abscontmin=interval*floor(MinA/interval);


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


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


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


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


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


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


 877 


[782]  878  sizpx=(Coord_x(end)Coord_x(1))/(np(2)1);


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


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


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


[688]  882 


 883  %axes(haxes)% set the input axes handle as current axis


 884 


 885  % colormap(map);


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


 887  Opacity=1;


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


 889  Opacity=PlotParam.Scalar.Opacity;


 890  end


 891  % fill the space between contours if opacity is undefined or =1


 892  if isequal(Opacity,1)


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


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


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


 896  hold on


 897  end


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


 899  hold on


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


 901  if isequal(Opacity,1)


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


 903  end


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


 905  hold off


 906 


 907  %determine the color scale and map


 908  caxis([abscontmin abscontmax])


 909  if BW


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


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


 912  colormap(map);


[688]  913  else


 914  colormap('default'); % default matlab colormap ('jet')


[201]  915  end


[688]  916 


[748]  917  if isfield(PlotParam.Axes,'CheckFixAspectRatio') && isequal(PlotParam.Axes.CheckFixAspectRatio,1)


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


[748]  919  if isfield(PlotParam.Axes,'AspectRatio')


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


[428]  921  else


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


 923  end


 924  end


[688]  925  else


[201]  926  % set colormap for image display


 927  % rescale the grey levels with min and max, put a grey scale colorbar


[688]  928  % B=A;


[201]  929  if BW


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


 931  map=[vec' vec' vec'];


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


[201]  933  else


[688]  934  colormap('default'); % standard false colors for div, vort , scalar fields


[201]  935  end


 936  end


 937 


[428]  938  % case of color images


 939  else


[201]  940  if BW


[688]  941  A=uint16(sum(A,3));


[201]  942  else


[688]  943  A=uint8(A);


[201]  944  end


 945  MinA=0;


 946  MaxA=255;


 947  end


 948 


 949  % display usual image


[428]  950  if ~CheckContour


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


[652]  952  test_interp=0;


[688]  953  if size(A,3)==1 % scalar of B/W image


[652]  954  test_interp=1;


 955  if max(np) <= 64


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


 957  elseif max(np) <= 128


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


 959  elseif max(np) <= 256


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


 961  else


 962  npxy=np;


 963  test_interp=0; % no interpolation done


 964  end


[201]  965  end


[652]  966  if test_interp%if we interpolate


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


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


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


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


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


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


[428]  973  end


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


[428]  975  if isempty(hima)


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


 977  if MinA<MaxA


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


[201]  979  else % to deal with uniform field


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


[201]  981  end


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


[411]  983  % requested:


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


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


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


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


[428]  988  % update an existing image


[201]  989  else


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


[201]  991  if MinA<MaxA


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


 993  else


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


[201]  995  end


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


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


[201]  998  end


[546]  999 


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


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


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


[405]  1003  else


[428]  1004  if test_vec


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


 1006  PlotParamOut.Scalar.Opacity=0.5;


 1007  else


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


 1009  end


[405]  1010  end


[201]  1011  end


 1012  test_ima=1;


 1013 


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


[546]  1015  if ~isempty(PosColorbar)


 1016  if siz==2 && exist('PosColorbar','var')


 1017  if isempty(hcol)~ishandle(hcol)


 1018  hcol=colorbar;%create new colorbar


 1019  end


 1020  if length(PosColorbar)==4


 1021  set(hcol,'Position',PosColorbar)


 1022  end


 1023  %YTick=0;%default


 1024  if MaxA>MinA


 1025  if CheckContour


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


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


 1028  YTick=cont_pos(2:end1);


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


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


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


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


 1033  if iscell(hi)%multiple images in colorbar


 1034  hi=hi{1};


 1035  end


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


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


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


 1039  YTick=colbartick(MinA,MaxA);


 1040  set(hcol,'YTick',YTick)


 1041  else


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


 1043  if iscell(hi)%multiple images in colorbar


 1044  hi=hi{1};


 1045  end


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


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


 1048  YTick=colbartick(MinA,MaxA);


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


 1050  set(hcol,'YTick',YTick)


[201]  1051  end


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


[201]  1053  end


[546]  1054  elseif ishandle(hcol)


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


[201]  1056  end


 1057  end


 1058  else%no scalar plot


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


[201]  1060  delete(hima)


 1061  end


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


[428]  1063  delete(hcol)


[201]  1064  end


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


 1066  end


 1067 


 1068  %% vector plot %%%%%%%%%%%%%%%%%%%%%%%%%%


 1069  if test_vec


 1070  %vector scale representation


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


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


 1073  end


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


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


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


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


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


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


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


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


 1082  else


 1083  if ~test_false %remove false vectors


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


 1085  end


 1086  if isempty(vec_U)


 1087  scale=1;


 1088  else


 1089  if isempty(indsel)


 1090  MaxU=max(abs(vec_U));


 1091  MaxV=max(abs(vec_V));


 1092  else


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


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


 1095  end


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


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


 1098  end


[210]  1099  end


[201]  1100 


 1101  %record vectors on the plotting axes


 1102  if test_C==0


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


 1104  end


 1105 


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


[581]  1107  check_decimate=0;


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


[581]  1109  check_decimate=1;


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


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


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


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


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


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


 1116  end


 1117  nb_sel=length(ind_sel);


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


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


 1120  check_decimate=1;


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


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


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


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


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


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


 1127  end


 1128  nb_sel=length(ind_sel);


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


 1130  end


 1131  if check_decimate


[201]  1132  vec_X=vec_X(ind_sel);


 1133  vec_Y=vec_Y(ind_sel);


 1134  vec_U=vec_U(ind_sel);


 1135  vec_V=vec_V(ind_sel);


 1136  vec_C=vec_C(ind_sel);


 1137  if ~isempty(ivar_F)


 1138  vec_F=vec_F(ind_sel);


 1139  end


[674]  1140  if ~isempty(ivar_FF_vec)


[201]  1141  vec_FF=vec_FF(ind_sel);


 1142  end


 1143  end


 1144 


 1145  %get main level color code


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


 1147 


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


 1149  sizlist=size(colorlist);


 1150  nbcolor=sizlist(1);


 1151  if test_black


 1152  nbcolor=nbcolor+1;


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


[674]  1154  if ~isempty(ivar_FF_vec)


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


[201]  1156  else


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


[201]  1158  end


 1159  end


 1160  nbcolor=nbcolor+1;


[674]  1161  if ~isempty(ivar_FF_vec)


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


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


 1164  else


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


 1166  end


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


[201]  1168  end


 1169  %plot vectors:


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


 1171 


 1172  else


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


 1174  if ~isempty(hvec)


 1175  delete(hvec);


 1176  end


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


 1178  end


[546]  1179  % nbvar=0;


[201]  1180 


 1181  %store the coordinate extrema occupied by the field


 1182  if ~isempty(Data)


[569]  1183  MinX=[];


 1184  MaxX=[];


 1185  MinY=[];


 1186  MaxY=[];


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


[201]  1188  if fix_lim


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


 1190  MinX=PlotParam.Axes.MinX;


 1191  MaxX=PlotParam.Axes.MaxX;


 1192  MinY=PlotParam.Axes.MinY;


 1193  MaxY=PlotParam.Axes.MaxY;


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


[252]  1195  else


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


[782]  1197  MinX=min(Coord_x);


 1198  MaxX=max(Coord_x);


 1199  MinY=min(Coord_y);


 1200  MaxY=max(Coord_y);


[252]  1201  if test_vec


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


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


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


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


[252]  1206  end


 1207  elseif test_vec


[569]  1208  MinX=min(vec_X);


 1209  MaxX=max(vec_X);


 1210  MinY=min(vec_Y);


 1211  MaxY=max(vec_Y);


[201]  1212  end


[252]  1213  end


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


 1215  PlotParamOut.Axes.MaxX=MaxX;


 1216  PlotParamOut.Axes.MinY=MinY;


 1217  PlotParamOut.Axes.MaxY=MaxY;


[569]  1218  if MaxX>MinX


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


[546]  1220  end


[569]  1221  if MaxY>MinY


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


[546]  1223  end


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


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


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


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


 1228  PlotParamOut.Axes.y_units=y_units;


[201]  1229  end


[760]  1230  % if isfield(PlotParam,'Axes') && isfield(PlotParam.Axes,'CheckFixAspectRatio') && isequal(PlotParam.Axes.CheckFixAspectRatio,1)


 1231  % set(haxes,'DataAspectRatioMode','manual')


 1232  % if isfield(PlotParam.Axes,'AspectRatio')


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


 1234  % end


 1235  % else


 1236  % set(haxes,'DataAspectRatioMode','auto')


 1237  % end


 1238 


[201]  1239  %


 1240  %  function for plotting vectors


 1241  %INPUT:


 1242  % haxes: handles of the plotting axes


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


 1244  % scale: scaling factor for vector length representation


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


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


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


 1248  %


 1249  %define arrows


 1250  theta=0.5 ;%angle arrow


 1251  alpha=0.3 ;%length arrow


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


 1253  %find the existing lines


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


 1255  sizh=size(h);


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


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


 1258 


 1259  %drawnow


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


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


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


 1263  end


 1264  sizlist=size(colorlist);


 1265  ncolor=sizlist(1);


 1266 


 1267  for icolor=1:ncolor


 1268  %determine the line positions for each color icolor


 1269  ind=find(col_vec==icolor);


 1270  xc=x(ind);


 1271  yc=y(ind);


 1272  uc=u(ind)*scale;


 1273  vc=v(ind)*scale;


 1274  n=size(xc);


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


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


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


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


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


 1280 


 1281  %determine arrow heads


 1282  arrowplus=rot*[uc;vc];


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


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


 1285  x2=xc+uc/2;


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


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


 1288  y2=yc+vc/2;


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


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


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


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


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


 1294  %draw the line or modify the existing ones


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


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


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


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


 1299  if n(2)>0


 1300  hold on


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


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


 1303  end


 1304  end


 1305  else


 1306  if isn(1)


 1307  delete(h(2*icolor1))


 1308  delete(h(2*icolor))


 1309  else


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


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


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


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


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


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


 1316  end


 1317  end


 1318  end


 1319  if sizh(1) > 2*ncolor


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


 1321  delete(h(2*icolor1))


 1322  delete(h(2*icolor))


 1323  end


 1324  end


 1325 


 1326  %


 1327  %  determine tick positions for colorbar


 1328  function YTick=colbartick(MinA,MaxA)


 1329  %


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


 1331  YTick=0;%default


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


 1333  if maxabs>0


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


 1335  div=1;


 1336  siz2=1;


 1337  while siz2<2


 1338  values=10:div:10;


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


 1340  siz=size(ind);


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


 1342  values=9:0.5*div:9;


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


 1344  end


 1345  siz2=size(ind,2);


 1346  div=div/10;


 1347  end


 1348  YTick=ord*values(ind);


 1349  end


 1350 


[356]  1351  % 


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


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


[389]  1354  % 


[356]  1355  if length(vec_Y)<2


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


 1357  return;


 1358  end


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


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


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


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


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


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


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


 1366  rangx=rgx_in; % first and last positions


 1367  rangy=rgy_in;


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


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


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


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


 1372  else % interpolation grid automatically determined


 1373  rangx(1)=min(vec_X);


 1374  rangx(2)=max(vec_X);


 1375  rangy(2)=min(vec_Y);


 1376  rangy(1)=max(vec_Y);


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


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


 1379  end


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


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


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


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


 1384  else


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


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


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


 1388  ymin=vec_Y(index(1));


 1389  y=vec_Y(index);


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


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


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


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


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


 1395 


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


 1397  %to prescribed positions


 1398  test_interp=1;


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


 1400  rangx=rgx_in; % first and last positions


 1401  rangy=rgy_in;


 1402  npxy=npxy_in;


 1403  else


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


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


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


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


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


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


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


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


 1412  else


 1413  npxy=[ny nx];


 1414  test_interp=0; % no interpolation done


 1415  end


 1416  end


 1417  if test_interp==1%if we interpolate


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


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


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


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


 1422  else %no interpolation for a resolution higher than 256


 1423  A=B;


 1424  end


[801]  1425  end

