 function cont_fphc(datafile,Tmax,Tstep,PLOTS,MOVIE)

% function cont_fpc(datafile,Tmax,Tstep,PLOTS,MOVIE)
%
% Continues the previous experiment whose solution and setup
% are in contained in the data file.
% Time interval is linspace(0,Tmax,Tstep)
% PLOTS and MOVIE are flags (0/off | 1/on) to produce each
%
% No changes are made in parameter values

% Input check
switch nargin
   case 0
      error('You need to specify the data file');
   case 1
      Tmax =  1; Tstep = 3; PLOTS = 0; MOVIE = 0;
   case 2
      Tstep = 3; PLOTS = 0; MOVIE = 0;
   case 3
      PLOTS = 0; MOVIE = 0;
   otherwise
      MOVIE = 0;
end

load(datafile);

% set up some basic values
fem0 = fem;
oldtmax = fem.sol.tlist(end);
flowmax = fem.flowmax;
[Re,Pr,Gamma,T0,T1,D,d,B,Umax] = deal(fem.const{2:2:end});

% display current values of the constant
disp('The current values of the constants are:');
disp(fem.const);

% Map initial solution to current xmesh
fem0.sol.u(:,1:end-1)=[]; % clears old fem0.sol.u - not necessary here
umap={'T',1,1,1,1;'p',1,1,1,1;'u',1,1,1,1;'v',1,1,1,1};
init=asseminit(fem,'init', {fem0,umap});

% Restart and compute a time-dependent simulation

% Solve dynamic problem
tlist = linspace(0,Tmax,Tstep);
fem.sol = femtime(fem,'tlist',tlist,...
       'ode','fldaspk',...
       'atol',{'u' 1e-3*flowmax 'v' 1e-3*flowmax 'p' inf 'T' 1e-3},...
       'stop' , 'on',...
       'init',init,...
       'timeind','D',...
       'report','on');

% save for restart purposes
fem0=fem;

if PLOTS == 1
   figure;
 % Plot initial solution
   postplot(fem,'tridata','T',...
             'tribar','on',...
             'contdata','p',...
             'contbar','on',...
             'flowdata',{ 'u' 'v'},...
             'title', 't = 0 Surface: Temp (T) Contour: pressure (p) Flow: [x velocity (u),y velocity (v)]  ',...
             'solnum', 1,...
             'axisequal','on');
   figure;
   % Plot solution
   postplot(fem,'tridata','T',...
             'tribar','on',...
             'contdata','p',...
             'contbar','on',...
             'flowdata',{ 'u' 'v'},...
             'title', 'Surface: Temp (T) Contour: pressure (p)  Flow: [x velocity (u),y velocity (v)]  ',...
             'axisequal','on');

end

if MOVIE == 1
   figure;
   Mov = postmovie(fem,'tridata','T',...
                     'tribar','on',...
                     'contdata','p',...
                     'contbar','on',...
                     'flowdata',{'u' 'v'},...
                     'title', 'Surface: Temp (T) Contour: pressure (p) Flow: [x velocity (u),y velocity (v)]  ',...
                     'axisequal','on');
end

fname = sprintf('cfphc_R%d_G%d_B%d_t%f.mat',Re,Gamma,B*100,oldtmax+Tmax);
save(fname,'fem');


