clear all;
close all;

nplot = readv('fort.nplot');
x = readv('fort.x');
y = readv('fort.y');
[mx,my,dx,dy,ax,ay,bx,by,iprob,fname] = read_dat;
[ha, hf] = get_figs(iprob);  %  Figure and axis handles

set(ha,'XLim',[ax bx]);
set(ha,'YLim',[ay by]);
set(ha,'XTick',[ax:bx]);
set(ha,'YTick',[ay:by]);


% Get mask for plotting contour plots.
[t,mxs,mys,data] = read_data([fname,'fort.kappa']);
kappa = reshape(data,mx,my);
kpa = sparse(kappa);
[si,sj] = find(kpa);
nz = nnz(kpa);
mask = nan*kappa;  % Mask out entries
for i = 1:nz
  mask(si(i),sj(i)) = 1.0;
end;

xe = ax:dx:bx;
ye = ay:dy:by;

[xs, ys] = init_shapes(iprob);

% 1 = pseudocolor plot; 2 = contour plot; 3 = mesh plot

if (iprob <= 2)
  p_choice = 2;
else
  p_choice = 1;
end;

for n = 0:nplot,
  [t,mxs,mys,data] = read_data([fname,'fort'],100+n);
  q = reshape(data(:,1),mxs,mys);
  if (p_choice == 1)
    qaug = [[q zeros(mx,1)]; zeros(1,my+1)];
    phandle = pcolor(xe',ye',qaug');
    set(phandle,'EdgeColor','none')
    shading('flat');
    yrbcolormap;
    hold on;
    % plot_shapes(iprob);
    plot(xs,ys,'r-');
  elseif (p_choice == 2)
    contour(x,y,(mask.*q)', [0.01 0.1:0.1:0.9],'b');
    hold on;
    plot(xs,ys,'r-');
    % plot_shapes(iprob);
  elseif (p_choice == 3)
    phandle = mesh(x',y',(q.*mask)');
    colormap([0 0 1]);
    view(20,30);
  end;
  set(ha,'XLim',[ax bx]);
  set(ha,'YLim',[ay by]);
  set(ha,'XTick',[ax:bx]);
  set(ha,'YTick',[ay:by]);
  if (p_choice == 3)
    set(ha,'ZLim',[0 1]);
    set(ha,'ZTick',[0:0.5:1]);
  end;
  caxis([0 1]);
  title(['t = ',num2str(t)]);
  if (n < nplot)
    query1;
  end;
  drawnow;
  hold off;
end;
