Index: trunk/src/plot_field.m
===================================================================
--- trunk/src/plot_field.m	(revision 109)
+++ trunk/src/plot_field.m	(revision 128)
@@ -105,5 +105,5 @@
 % use htext: handles of the text edit box (uicontrol)
 % introduce PlotParam.Hold: 'on' or 'off' (for curves)
-  
+
 %default output
 if ~exist('PlotParam','var'),PlotParam=[];end;
@@ -127,5 +127,5 @@
     if ishandle(haxes)
         if isequal(get(haxes,'Type'),'axes')
-            axes(haxes)
+%             axes(haxes)
             testnewfig=0;
             AxeData=get(haxes,'UserData');
@@ -139,4 +139,5 @@
     end
 end
+
 if testnewfig% create a new figure and axes if the plotting axes does not exist
     hfig=figure;
@@ -157,6 +158,4 @@
     haxes=axes;
     set(haxes,'position',[0.13,0.2,0.775,0.73])
-else
-    hstack=findobj(allchild(0),'Type','figure');%current stack order of figures in matlab
 end
 if isfield(PlotParam,'text_display_1') && ishandle(PlotParam.text_display_1)
@@ -182,9 +181,11 @@
     [Data.NbDim]=max(NbDim);
 end
+errormsg=[];
 if isequal(Data.NbDim,0) 
         AxeData=plot_text(Data,htext);
         PlotType='text';
+        errormsg=[];
 elseif isequal(Data.NbDim,1)
-    [AxeData,haxes]=plot_profile(Data,CellVarIndex,VarType,haxes,PlotParam);% 
+    [AxeData]=plot_profile(Data,CellVarIndex,VarType,haxes,PlotParam);% 
     if testzoomaxes
         [AxeData,zoomaxes,PlotParamOut]=plot_profile(Data,CellVarIndex,VarType,zoomaxes,PlotParam);
@@ -192,22 +193,23 @@
     end
     PlotType='line';
+    errormsg=[];
 elseif isequal(Data.NbDim,2)
     ind_select=find(NbDim>=2);
     if numel(ind_select)>2
-        msgbox_uvmat('ERROR','more than two fields to map')
-        display('more than two fields to map')
-        return
-    end
-    [AxeData,haxes,PlotParamOut,PlotType]=plot_plane(Data,CellVarIndex(ind_select),VarType(ind_select),haxes,PlotParam,htext,PosColorbar);
-    if testzoomaxes
-        [AxeData,zoomaxes,PlotParamOut]=plot_plane(Data,CellVarIndex(ind_select),VarType(ind_select),zoomaxes,PlotParam,1,PosColorbar);
-        %AxeData.ZoomAxes=zoomaxes;
-        Data.ZoomAxes=zoomaxes;
+        errormsg='more than two fields to map';
+    else
+        [AxeData,xx,PlotParamOut,PlotType,errormsg]=plot_plane(Data,CellVarIndex(ind_select),VarType(ind_select),haxes,PlotParam,htext,PosColorbar);
+        if testzoomaxes && isempty(errormsg)
+            [AxeData,zoomaxes,PlotParamOut,xx,errormsg]=plot_plane(Data,CellVarIndex(ind_select),VarType(ind_select),zoomaxes,PlotParam,1,PosColorbar);
+            %AxeData.ZoomAxes=zoomaxes;
+            Data.ZoomAxes=zoomaxes;
+        end
     end
 elseif isequal(Data.NbDim,3)
+    errormsg='volume plot not implemented yet';
+end
+if ~isempty(errormsg)
     msgbox_uvmat('ERROR','volume plot not implemented yet')
     return
-% else
-%     testnbdim=0;
 end
 
@@ -234,8 +236,8 @@
 %set(haxes,'UserData',AxeData)
 set(haxes,'UserData',Data)
-if ~testnewfig
-set(0,'Children',hstack);%put back the initial figure stack after plot creation
-%set(haxes,'Tag','uvmat'); 
-end
+% if ~testnewfig
+% %set(0,'Children',hstack);%put back the initial figure stack after plot creation
+% 
+% end
 
 
@@ -303,5 +305,6 @@
 function [AxeData,haxes]=plot_profile(data,CellVarIndex,VarType,haxes,PlotParam)
 %-------------------------------------------------------------------
-%axes(haxes)
+%TODO: modify existing plot if it exists
+
 hfig=get(haxes,'parent');
 AxeData=data;
@@ -322,8 +325,8 @@
 legend_str={};
 
-%initiates string of the plot command
+%% prepare the string for plot command
+%initiate  the plot command
 plotstr='hhh=plot(';
 textmean={};
-% abscissa_name='';
 coord_x_index=[];
 test_newplot=1;
@@ -331,12 +334,8 @@
 % num_curve=numel(hh);
 % icurve=0;
+
+%loop on input  fields
 for icell=1:length(CellVarIndex)
-    testfalse=0;
     VarIndex=CellVarIndex{icell};%  indices of the selected variables in the list data.ListVarName
-    DimCell=data.VarDimName{VarIndex(1)};
-    if ischar(DimCell)
-        DimCell={DimCell};
-    end
-    XName=DimCell{1}; %first dimension considered as abscissa
     if ~isempty(VarType{icell}.coord_x)
         coord_x_index=VarType{icell}.coord_x;
@@ -353,94 +352,82 @@
     end
     testplot=ones(size(data.ListVarName));%default test for plotted variables
-    testcoordvar=0;
-    charplot_0='''-''';%default
-    if isfield(data,'ObjectProjMode')&& isequal(data.ObjectProjMode,'projection')
-        charplot_0='+';
-    end
-    xtitle='';  
-    
     xtitle=data.ListVarName{coord_x_index};
     eval(['coord_x{icell}=data.' data.ListVarName{coord_x_index} ';']);%coordinate variable set as coord_x
     if isfield(data,'VarAttribute')&& numel(data.VarAttribute)>=coord_x_index && isfield(data.VarAttribute{coord_x_index},'units')
-         xtitle=[xtitle '(' data.VarAttribute{coord_x_index}.units ')'];
+        xtitle=[xtitle '(' data.VarAttribute{coord_x_index}.units ')'];
     end
     eval(['coord_x{icell}=data.' data.ListVarName{coord_x_index} ';']);%coordinate variable set as coord_x
-%     testcoordvar=1;
     testplot(coord_x_index)=0;
     if ~isempty(VarType{icell}.ancillary')
-            testplot(VarType{icell}.ancillary)=0;
+        testplot(VarType{icell}.ancillary)=0;
     end
     if ~isempty(VarType{icell}.warnflag')
-            testplot(VarType{icell}.warnflag)=0;
-    end
- 
+        testplot(VarType{icell}.warnflag)=0;
+    end
     if isfield(data,'VarAttribute')
         VarAttribute=data.VarAttribute;
-        for ivar=1:length(VarIndex) 
-             if length(VarAttribute)>=VarIndex(ivar) && isfield(VarAttribute{VarIndex(ivar)},'long_name')
-                 plotname{VarIndex(ivar)}=VarAttribute{VarIndex(ivar)}.long_name;
-             else
-                 plotname{VarIndex(ivar)}=data.ListVarName{VarIndex(ivar)};%name for display in plot A METTRE
-             end
-        end
-    end  
-%     test_newplot=0;%default
-%     if num_curve>=icurve+numel(find(testplot(VarIndex)))%update existing curves
-%         if ~isempty(VarType{icell}.discrete')
-%             charplot_0='+';
-%             LineStyle='none';
-%         else
-%             charplot_0='none';
-%             LineStyle='-';
-%         end
-%         for ivar=1:length(VarIndex)
-%             if testplot(VarIndex(ivar))
-%                 icurve=icurve+1;
-%                 VarName=data.ListVarName{VarIndex(ivar)};
-%                 eval(['data.' VarName '=squeeze(data.' VarName ');'])
-%                 set(hh(icurve),'LineStyle',LineStyle)
-%                 set(hh(icurve),'Marker',charplot_0)
-%                 set(hh(icurve),'XData',coord_x{icell})
-%                 eval(['yy=data.' VarName ';'])
-%                 set(hh(icurve),'YData',yy);
-%             end
-%         end
-%     else% new plot
-        if ~isempty(VarType{icell}.discrete')
-            charplot_0='''+''';
-        else
-            charplot_0='''-''';
-        end
         for ivar=1:length(VarIndex)
-            if testplot(VarIndex(ivar))
-                VarName=data.ListVarName{VarIndex(ivar)};
-                eval(['data.' VarName '=squeeze(data.' VarName ');'])
-%                 if isequal(VarName,'A')
-%                     charplot='''-''';
-%                 else
-%                     charplot=charplot_0;
-%                 end
-                plotstr=[plotstr 'coord_x{' num2str(icell) '},data.' VarName ',' charplot_0 ','];
-                eval(['nbcomponent2=size(data.' VarName ',2);']);
-                eval(['nbcomponent1=size(data.' VarName ',1);']);
-                if numel(coord_x{icell})==2
-                    coord_x{icell}=linspace(coord_x{icell}(1),coord_x{icell}(2),nbcomponent1);
-                end
-                eval(['varmean=mean(double(data.' VarName '));']);%mean value
-                textmean=[textmean; {[VarName 'mean= ' num2str(varmean,4)]}];
-                if nbcomponent1==1|| nbcomponent2==1
-                    legend_str=[legend_str {VarName}]; %variable with one component
-                else  %variable with severals  components
-                    for ic=1:min(nbcomponent1,nbcomponent2)
-                        legend_str=[legend_str [VarName '_' num2str(ic)]]; %variable with severals  components 
-                    end                                                   % labeled by their index (e.g. color component)
-                end
-            end
-        end
-%     end
-end
+            if length(VarAttribute)>=VarIndex(ivar) && isfield(VarAttribute{VarIndex(ivar)},'long_name')
+                plotname{VarIndex(ivar)}=VarAttribute{VarIndex(ivar)}.long_name;
+            else
+                plotname{VarIndex(ivar)}=data.ListVarName{VarIndex(ivar)};%name for display in plot A METTRE
+            end
+        end
+    end
+    %     test_newplot=0;%default
+    %     if num_curve>=icurve+numel(find(testplot(VarIndex)))%update existing curves
+    %         if ~isempty(VarType{icell}.discrete')
+    %             charplot_0='+';
+    %             LineStyle='none';
+    %         else
+    %             charplot_0='none';
+    %             LineStyle='-';
+    %         end
+    %         for ivar=1:length(VarIndex)
+    %             if testplot(VarIndex(ivar))
+    %                 icurve=icurve+1;
+    %                 VarName=data.ListVarName{VarIndex(ivar)};
+    %                 eval(['data.' VarName '=squeeze(data.' VarName ');'])
+    %                 set(hh(icurve),'LineStyle',LineStyle)
+    %                 set(hh(icurve),'Marker',charplot_0)
+    %                 set(hh(icurve),'XData',coord_x{icell})
+    %                 eval(['yy=data.' VarName ';'])
+    %                 set(hh(icurve),'YData',yy);
+    %             end
+    %         end
+    %     else% new plot
+    if ~isempty(VarType{icell}.discrete')
+        charplot_0='''+''';
+    else
+        charplot_0='''-''';
+    end
+    for ivar=1:length(VarIndex)
+        if testplot(VarIndex(ivar))
+            VarName=data.ListVarName{VarIndex(ivar)};
+            eval(['data.' VarName '=squeeze(data.' VarName ');'])
+            plotstr=[plotstr 'coord_x{' num2str(icell) '},data.' VarName ',' charplot_0 ','];
+            eval(['nbcomponent2=size(data.' VarName ',2);']);
+            eval(['nbcomponent1=size(data.' VarName ',1);']);
+            if numel(coord_x{icell})==2
+                coord_x{icell}=linspace(coord_x{icell}(1),coord_x{icell}(2),nbcomponent1);
+            end
+            eval(['varmean=mean(double(data.' VarName '));']);%mean value
+            textmean=[textmean; {[VarName 'mean= ' num2str(varmean,4)]}];
+            if nbcomponent1==1|| nbcomponent2==1
+                legend_str=[legend_str {VarName}]; %variable with one component
+            else  %variable with severals  components
+                for ic=1:min(nbcomponent1,nbcomponent2)
+                    legend_str=[legend_str [VarName '_' num2str(ic)]]; %variable with severals  components
+                end                                                   % labeled by their index (e.g. color component)
+            end
+        end
+    end
+end
+
+%% activate the plot
 if test_newplot && ~isequal(plotstr,'hhh=plot(')
     plotstr=[plotstr '''tag'',''plot_line'');'];
                 %execute plot (instruction  plotstr)  
+    axes(haxes)% select the plotting axes for plot operation
     eval(plotstr)
    
@@ -502,6 +489,7 @@
 
 %-------------------------------------------------------------------
-function [AxeData,haxes,PlotParamOut,PlotType]=plot_plane(Data,CellVarIndex,VarTypeCell,haxes,PlotParam,htext,PosColorbar)
+function [AxeData,haxes,PlotParamOut,PlotType,errormsg]=plot_plane(Data,CellVarIndex,VarTypeCell,haxes,PlotParam,htext,PosColorbar)
 %-------------------------------------------------------------------
+
 grid(haxes, 'off')
 %default plotting parameters
@@ -517,6 +505,4 @@
 end
 PlotParamOut=PlotParam;%default
-
-%plotting axes
 hfig=get(haxes,'parent');
 hcol=findobj(hfig,'Tag','Colorbar'); %look for colorbar axes
@@ -530,5 +516,5 @@
     AxeData.ObjectCoord=Data.ObjectCoord;
 end
-
+errormsg=[];%default
 test_ima=0; %default: test for image or map plot
 test_vec=0; %default: test for vector plots
@@ -548,5 +534,5 @@
     ivar_C=[VarType.scalar VarType.image VarType.color VarType.ancillary]; %defines index (indices) for the scalar or ancillary fields
     if numel(ivar_C)>1
-        msgbox_uvmat('ERROR','error in plot_field: too many scalar inputs')
+        errormsg= 'error in plot_field: too many scalar inputs';
         return
     end
@@ -558,8 +544,8 @@
     end
 %     idim_Y=[];  
-    test_grid=0;
+%     test_grid=0;
     if ~isempty(ivar_U) && ~isempty(ivar_V)% vector components detected
         if test_vec
-            msgbox_uvmat('ERROR','error in plot_field: attempt to plot two vector fields')
+            errormsg='error in plot_field: attempt to plot two vector fields';
             return
         else
@@ -570,5 +556,5 @@
                 eval(['vec_X=Data.' Data.ListVarName{ivar_X} ';']) 
                 eval(['vec_Y=Data.' Data.ListVarName{ivar_Y} ';'])
-            elseif numel(VarType.coord)==2 & VarType.coord~=[0 0];%coordinates defines by dimension variables
+            elseif numel(VarType.coord)==2 && ~isequal(VarType.coord,[0 0]);%coordinates defines by dimension variables
                 eval(['y=Data.' Data.ListVarName{VarType.coord(1)} ';']) 
                 eval(['x=Data.' Data.ListVarName{VarType.coord(2)} ';'])
@@ -581,5 +567,5 @@
                 [vec_X,vec_Y]=meshgrid(x,y);  
             else
-                msgbox_uvmat('ERROR','error in plot_field: invalid coordinate definition for vector field')
+                errormsg='error in plot_field: invalid coordinate definition for vector field';
                 return
             end
@@ -605,5 +591,5 @@
     elseif ~isempty(ivar_C) %scalar or image
         if test_ima
-             msgbox_uvmat('ERROR','attempt to plot two scalar fields or images')
+             errormsg='attempt to plot two scalar fields or images';
             return
         end
@@ -635,5 +621,5 @@
                     x_units=['(' Data.VarAttribute{VarType.coord(2)}.units ')'];
                 end
-                if numel(Data.VarAttribute)>=VarType.coord(1) & isfield(Data.VarAttribute{VarType.coord(1)},'units')
+                if numel(Data.VarAttribute)>=VarType.coord(1) && isfield(Data.VarAttribute{VarType.coord(1)},'units')
                     y_units=['(' Data.VarAttribute{VarType.coord(1)}.units ')'];
                 end
@@ -644,5 +630,5 @@
                 DAY_max=max(DAY);
                 if sign(DAY_min)~=sign(DAY_max);% =1 for increasing values, 0 otherwise
-                     errormsg=['errror in plot_field.m: non monotonic dimension variable # ' ListVarName{VarType.coord(1)} ];
+                     errormsg=['errror in plot_field.m: non monotonic dimension variable ' Data.ListVarName{VarType.coord(1)} ];
                       return
                 end 
@@ -654,5 +640,5 @@
                 DAX_max=max(DAX);
                 if sign(DAX_min)~=sign(DAX_max);% =1 for increasing values, 0 otherwise
-                     errormsg=['errror in plot_field.m: non monotonic dimension variable # ' ListVarName{VarType.coord(2)} ];
+                     errormsg=['errror in plot_field.m: non monotonic dimension variable ' Data.ListVarName{VarType.coord(2)} ];
                       return
                 end 
@@ -682,12 +668,12 @@
             AY=[AY(1) AY(end)];
         else
-            msgbox_uvmat('ERROR','error in plot_field: invalid coordinate definition ')
+            errormsg='error in plot_field: invalid coordinate definition ';
             return
         end
-          x_label=[Data.ListVarName{ivar_X} '(' x_units ')'];
+          %x_label=[Data.ListVarName{ivar_X} '(' x_units ')'];
     end          
 end 
 
-%%%%%%%%%%%%%%%%%%%%%   image or scalar plot %%%%%%%%%%%%%%%%%%%%%%%%%%
+%%   image or scalar plot %%%%%%%%%%%%%%%%%%%%%%%%%%
 
 if ~isfield(PlotParam.Scalar,'Contours')
@@ -700,5 +686,5 @@
     siz=numel(np);
     if siz>3
-        msgbox_uvmat('ERROR',['unrecognized scalar type: ' num2str(siz) ' dimensions'])
+       errormsg=['unrecognized scalar type: ' num2str(siz) ' dimensions'];
             return
     end
@@ -709,8 +695,9 @@
             siz=3;%color image
         else
-             msgbox_uvmat('ERROR',['unrecognized scalar type: ' num2str(np(3)) ' color components'])
+            errormsg=['unrecognized scalar type: ' num2str(np(3)) ' color components'];
             return
         end
     end
+    
     %set the color map
     if isfield(PlotParam.Scalar,'BW')
@@ -719,5 +706,7 @@
         BW=(siz==2) && (isa(A,'uint8')|| isa(A,'uint16'));% non color images represented in gray scale by default
     end
-    if siz==2 %for black and white images
+    
+    %case of grey level images or contour plot
+    if siz==2 
         if ~isfield(PlotParam.Scalar,'AutoScal')
             PlotParam.Scalar.AutoScal=0;%default
@@ -729,5 +718,5 @@
             PlotParam.Scalar.MaxA=[];%default
         end
-        if isequal(PlotParam.Scalar.AutoScal,0)|isempty(PlotParam.Scalar.MinA)|~isa(PlotParam.Scalar.MinA,'double')  %correct if there is no numerical data in edit box
+        if isequal(PlotParam.Scalar.AutoScal,0)||isempty(PlotParam.Scalar.MinA)||~isa(PlotParam.Scalar.MinA,'double')  %correct if there is no numerical data in edit box
             MinA=double(min(min(A)));
         else
@@ -741,7 +730,8 @@
         PlotParamOut.Scalar.MinA=MinA;
         PlotParamOut.Scalar.MaxA=MaxA;
-        axes(haxes)
+        
+        % case of contour plot
         if isequal(PlotParam.Scalar.Contours,1)
-            if ~isempty(hima) & ishandle(hima)
+            if ~isempty(hima) && ishandle(hima)
                 delete(hima)
             end
@@ -770,5 +760,5 @@
             txt=ver;%version of Matlab
             Release=txt(1).Release;
-            relnumb=str2num(Release(3:4));
+            relnumb=str2double(Release(3:4));
             if relnumb >= 14
                     vec=linspace(0,1,(abscontmax-abscontmin)/intercont);%define a greyscale colormap with steps intercont
@@ -788,4 +778,6 @@
             colormap(map);
         end
+        
+        % set  colormap for  image display
         if ~isequal(PlotParam.Scalar.Contours,1)  
             % rescale the grey levels with min and max, put a grey scale colorbar
@@ -799,6 +791,7 @@
             end
         end
-    else %color images
-        axes(haxes)
+        
+    % case of color images 
+    else 
         if BW
             B=uint16(sum(A,3));
@@ -809,6 +802,8 @@
         MaxA=255;
     end
-    if ~isequal(PlotParam.Scalar.Contours,1)
-        %interpolate to increase resolution
+    
+    % display usual image
+    if ~isequal(PlotParam.Scalar.Contours,1)      
+        % interpolate field to increase resolution of image display
         test_interp=1;
         if max(np) <= 64 
@@ -829,14 +824,18 @@
             yi=linspace(AY(1),AY(2),npxy(1));
             B = interp2(X,Y,double(B),xi,yi');
-        end        
+        end           
+        % create new image if there  no image handle is found
         if isempty(hima)
+           % axes(haxes)% set haxes the current axes for image creation
+            set(hfig,'CurrentAxes',haxes) % set haxes the current axes for image creation 
             tag=get(haxes,'Tag');
             if MinA<MaxA
                 hima=imagesc(AX,AY,B,[MinA MaxA]);
-            else
+            else % to deal with uniform field
                 hima=imagesc(AX,AY,B,[MaxA-1 MaxA]);
             end
             set(hima,'Tag','ima','HitTest','off')
-            set(haxes,'Tag',tag);%preserve the axes tag (removed by image fct !!!)       
+            set(haxes,'Tag',tag);%preserve the axes tag (removed by image fct !!!)            
+        % update an existing image
         else
             set(hima,'CData',B);
@@ -857,4 +856,5 @@
     AxeData.AY=[AY(1) AY(end)];
     test_ima=1;
+    
     %display the colorbar code for B/W images if Poscolorbar not empty
     if siz==2 && exist('PosColorbar','var')&& ~isempty(PosColorbar)
@@ -865,5 +865,5 @@
                  set(hcol,'Position',PosColorbar)           
         end 
-        YTick=0;%default
+        %YTick=0;%default
         if MaxA>MinA
             if isequal(PlotParam.Scalar.Contours,1)
@@ -873,5 +873,5 @@
                 YTick_scaled=colbarlim(1)+scale_bar*(YTick-abscontmin);
                 set(hcol,'YTick',YTick_scaled);
-            elseif (isfield(PlotParam.Scalar,'BW') & isequal(PlotParam.Scalar.BW,1))|isa(A,'uint8')| isa(A,'uint16')%images
+            elseif (isfield(PlotParam.Scalar,'BW') && isequal(PlotParam.Scalar.BW,1))||isa(A,'uint8')|| isa(A,'uint16')%images
                 hi=get(hcol,'children');
                 if iscell(hi)%multiple images in colorbar
@@ -879,5 +879,5 @@
                 end
                 set(hi,'YData',[MinA MaxA])
-                set(hi,'CData',[1:256]')
+                set(hi,'CData',(1:256)')
                 set(hcol,'YLim',[MinA MaxA])
                 YTick=colbartick(MinA,MaxA);
@@ -889,5 +889,5 @@
                 end
                 set(hi,'YData',[MinA MaxA])
-                set(hi,'CData',[1:64]')
+                set(hi,'CData',(1:64)')
                 YTick=colbartick(MinA,MaxA); 
                 set(hcol,'YLim',[MinA MaxA])
@@ -912,5 +912,5 @@
 end
 
-%%%%%%%%%%%%%%%%%%%%%   vector plot %%%%%%%%%%%%%%%%%%%%%%%%%%
+%%   vector plot %%%%%%%%%%%%%%%%%%%%%%%%%%
 if test_vec
    %vector scale representation
@@ -927,5 +927,5 @@
     if  ~isfield(PlotParam.Vectors,'AutoVec') || isequal(PlotParam.Vectors.AutoVec,0)|| ~isfield(PlotParam.Vectors,'VecScale')...
                ||isempty(PlotParam.Vectors.VecScale)||~isa(PlotParam.Vectors.VecScale,'double') %automatic vector scale
-        scale=[];
+%         scale=[];
         if test_false %remove false vectors
             indsel=find(AxeData.FF==0);%indsel =indices of good vectors
@@ -1045,5 +1045,5 @@
 listdim={'AY','AX',{'AY','AX'},'nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors'};
 Role={'coord_y','coord_x','scalar','coord_x','coord_y','vector_x','vector_y','scalar','vector_z','warnflag','errorflag'};
-ind_select=[];
+%ind_select=[];
 nbvar=0;
 AxeData.ListVarName={};
@@ -1168,5 +1168,9 @@
     isn=isnan(colorlist(icolor,:));%test if color NaN
     if 2*icolor > sizh(1) %if icolor exceeds the number of existing ones
-        axes(haxes)
+        %axes(haxes)
+        hfig=get(haxes,'parent');
+%         axes(haxes)
+        set(0,'CurrentFigure',hfig)
+        set(hfig,'CurrentAxes',haxes)
         if ~isn(1) %if the vectors are visible color not nan
             if n(2)>0
