% Plot the rate constant versus temperature for alanine dipeptide data.

clear;

% Load data from NetCDF file.
filename = 'alanine-dipeptide-transition-rate.nc';
nc = netcdf(filename);

% Extract variables into a structure.
nvars = length(nc.VarArray);
vars = struct();
for i = 1:nvars  
  vars = setfield(vars, nc.VarArray(i).Str, nc.VarArray(i).Data);
end

% Plot rate versus temperature
clf;
hold on;

% Plot reweighting estimate with shaded 95% confidence interval.
x = vars.interpolated_temperature_k';
y = vars.rate_k';
dy = vars.drate_k';
fill([x fliplr(x)], [(y+dy) fliplr(y-dy)], 0.7 * [1 1 1]);
plot(vars.interpolated_temperature_k, vars.rate_k, 'k-', 'LineWidth', 1.5);

% Plot single temperature estimates
errorbar(vars.temperature_k, vars.rate_single_k, vars.drate_single_k, 'k.', 'MarkerSize', 10, 'LineWidth', 1.2);

% Adjust axes limits.
oldaxis = axis;
axis([vars.temperature_k(1) vars.temperature_k(end) 0 oldaxis(4)]);

% Turn on box surrounding plot.
set(gca, 'box', 'on');

% Set axes ticks.
set(gca, 'YTick', [0.1, 0.2, 0.3, 0.4]);

% Label plot axes.
xlabel('temperature / K');
ylabel('rate / ps^{-1}');

% Generate plot
outfilename = 'alanine-dipeptide-transition-rate';
exportfig(gcf, sprintf('%s.eps', outfilename), 'width', 3.25);
unix(sprintf('epstopdf %s.eps', outfilename));
