% Plot steady-state power.

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('timestep-%04.2f.out', timestep)
	  x = textread(filename);
          color = colors(i);
          h = errorbar(x(:,1), x(:,5) ./ (timestep*x(:,1)), 2*x(:,6) ./ (timestep*x(:,1)), '.')
	  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('steady-state shadow work / (kT/fs)');
axis([-10 2100 -0.001 +0.03]);
plot([-10 2100], [0 0], 'k-');
title('steady-state power');

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