%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ploting U(x,y,t)
% call:
% plotU(x,y,u,ax,zl,title,mytext,plotop)
% plotop='phase','real','imag','height','mag2','mag','default'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plotU(x,y,u,ax,zl,title,mytext,plotop)

numcols=64;
n2=numcols/2;
xx=1:numcols/2;
blue_red=ones(numcols,3)*1;
% this gives NO white zone
blue_red(xx,1) =xx'/n2-1/n2;
blue_red(xx,2) =xx'/n2-1/n2;
blue_red(xx+n2,2) =1-xx'/n2;
blue_red(xx+n2,3) =1-xx'/n2;


 switch plotop
 case 'phase'
  surf(x,y,abs(u).^2,angle(u));
  shading interp
  colormap('jet')
%  colormap(blue_red)
 case 'real'
  surfl(x,y,real(u))
  shading interp
  colormap('gray')
 case 'imag'
  surfl(x,y,imag(u))
  shading interp
  colormap('gray')
 case 'height'
  surf(x,y,u)
  shading interp
  colormap('jet')
 case 'mag2'
  surfl(x,y,abs(u).^2)
  shading interp
  colormap('gray')
 case 'mag'
  surfl(x,y,u)
  shading interp
  colormap('gray')
 otherwise
  surfl(x,y,u)
  shading interp
  colormap('gray')
 end

 xlabel('x')
 ylabel('y')
 zlabel(zl)
 delta=1e-10;
% delta=10;
 if(isempty(ax))
  axis tight;
  ax=axis;
 else
  axis([ax(1)-delta ax(2)+delta ax(3)-delta ax(4)+delta]);
 end
 text(ax(2)-0.5*(ax(2)-ax(1)),ax(3)+0*(ax(4)-ax(3)),ax(6)+0*(ax(6)-ax(5)),mytext)
 axis(ax);


