% Plot cumulative pseudowork.

clear;

timesteps = {0.25, 0.50, 0.75, 1.00, 1.25, 1.50, 1.75, 2.00}; % timesteps to plot
ntimesteps = length(timesteps);
colors = jet(ntimesteps);

clf;
hold on;
for i = 1:ntimesteps
	  timestep = timesteps{i};
	  filename = sprintf('cumulative-%04.2f.out', timestep)
	  x = textread(filename);
          color = colors(i);
          h = errorbar(x(:,1), x(:,2), 2*x(:,3), '.')
	  set(h, 'color', colors(i,:));	 
end

timesteps_labels = cell(1,ntimesteps);
for i = 1:ntimesteps
	  timestep = timesteps{i};
          timesteps_labels{i} = sprintf('%4.2f fs', timestep);
end

h = legend(timesteps_labels);
set(h, 'location', 'northwest');

xlabel('number of timesteps');
ylabel('average shadow work / kT');
axis([-10 4200 -1 +250]);

exportfig(gcf, 'shadow-work.eps', 'width', 8, 'height', 6, 'color', 'cmyk');
system('epstopdf shadow-work.eps');
