
REstarA = [-1800  2400 ];
nR = length(REstarA);
tstop = [1020.13 1020.97 1021.22];


Bo=1000;
scalarz='z';
scalarp='p';

[Xb,Zb] = return_patch_bathy();

pflag = 1;

for k=1:1,

REstar = REstarA(k);
index_stop=return_basilisk_index(tstop(k));


figure(1),clf


 cmap = cbrewer('div','RdBu',80);
 colormap(flipud(cmap));
xSize = 16; ySize = 6;  %size of the figure to be generated
set(gcf,'PaperUnits','centimeters')
set(gcf,'PaperPosition',[0 0  xSize ySize])

index2=143;
for index=581:582,
    %for t=1011.0:1011.0, %0.5:1020.5,
    %  [Z,x,y,t]  = read_binary_file(9,REstar,Bo,index,scalarz);
  [P,x,y,t]  = read_binary_file(9,REstar,Bo,index,scalarp);  

  minx = min(x)+4.5; maxx=max(x)-2;
  
  [Xs,Ys,R] = read_processed_surface_extract(9,REstar,Bo,t);

  ii=find(Xs>=minx-1 & Xs<=maxx+1);
  Xsn = [Xs(ii(1)); Xs(ii); Xs(ii(end))];
 Ysn = [-1; Ys(ii); -1];

 if Bo==4001,
     Bo_n = 4000;
 else
     Bo_n = Bo;
 end;
 
 if (pflag==1),
     S = P;
     clim = [-0.08 0.08];
     titlestr =sprintf('pressure: $\\mathrm{Re}^*=%d$, $\\mathrm{Bo}=%d$: $t=%.2f$',REstar,Bo,t-1001);
     imgfile=sprintf('tmp/focusP_REstar%d_Bo%04d_%03d.png',REstar,Bo,index)     
 else
     S=Z;
     clim = [-10 10];
     titlestr =sprintf('vorticity: $\\mathrm{Re}^*=%d$, $\\mathrm{Bo}=%d$: $t=%.2f$',REstar,Bo,t-1001);
     imgfile=sprintf('tmp/focusZ_REstar%d_Bo%04d_%03d.png',REstar,Bo,index)     
 end;
 
    clf
   sand_color = [194 178 128]/255;
   water_color = [100 208 216]/255;   

   hh=subplot('position',[0.06 0.17 0.925 0.8])


            pcolorcen(x,y,S),shading interp
            axis xy
            axis([minx maxx -0.3 0.85])
            caxis(clim)
            hold on
            plot(Xsn,Ysn,'k','linewidth',3)
            h1=patch(Xsn,Ysn,water_color);
            %                       set(h1,'FaceAlpha',0.1)
            h2=patch(Xb,Zb,sand_color)
            colorbar('vert','position',[0.10 0.65 0.03 0.25])
            ylabel('$z/h_0$','interpreter','latex')
            xlabel('$x/h_0$','interpreter','latex')
            set(gca,'TickDir','out')
            text(minx+2.5,1.6,titlestr,'interpreter','latex')
            drawnow

           set(gcf,'PaperUnits','centimeters')
           set(gcf,'papersize',[xSize ySize],'paperposition',[ 0 0 xSize ySize]);

                                 print('-dpng','-r400',imgfile)
           index2=index2+1;            
end;
    end;

            return
               
%nn = length(Xn);
%Cn = [1:nn]'/nn;
    
% figure(1)
% hold on
%      scatter(Xs,Ys+(t-1002)/5,4,C,'d')
%     drawnow
%     pause(0.02)
%     %    plot(Spc.Location(:,1),Spc.Location(:,2),'kd')
% 
%     figure(2)
%     hold on
%     %for t=1011.0:1011.0, %0.5:1020.5,
% 
% nraw = length(Xraw);
% Cr = [1:nraw]'/nraw;
% %nn = length(Xn);
% %Cn = [1:nn]'/nn;
%     
% figure(2)
% hold on
%      scatter(Xraw,Yraw+(t-1002)/5,4,Cr,'d')
%     drawnow
%     pause(0.02)
%      
% end;
    
    
  %[F,x,y,t]  = read_binary_file(9,REstar,Bo,index,'f');
%[V,x,y,t]  = read_binary_file(9,REstar,Bo,index,'z')
%[P,x,y,t]  = read_binary_file(9,REstar,Bo,index,'p');
