% Plot the rate constant versus temperature for alanine dipeptide data.
% This script also generates an Arrhenius plot.

clear;

% Load data from NetCDF file.
filename = 'alanine-dipeptide-transition-rate-psi.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

%% Create figure.

%% Plot rate versus temperature
subplot(2,1,1);
hold on;

% Plot reweighting estimate with shaded 95% confidence interval.
x = vars.interpolated_temperature_k';
y = vars.rate_k';
dy = 2*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
x = vars.temperature_k;
y = vars.rate_single_k;
dy = 2*vars.drate_single_k;
% Remove those where rate estimate is zero.
indices = find(vars.rate_single_k ~= 0);
x = x(indices);
y = y(indices);
dy = dy(indices);
% plot
errorbar(x, y, dy, 'k.', 'MarkerSize', 10, 'LineWidth', 1.2);

% Adjust axes limits.
oldaxis = axis;
axis([vars.temperature_k(1) vars.temperature_k(end) min(y-dy) max(y+dy)]);

% 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('T (K)');
ylabel('k (1/ps)');

%% Plot Arrhenius plot (ln k vs 1/T):
subplot(2,1,2);
hold on;

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

% Plot single temperature estimates
x = (vars.temperature_k).^(-1);
y = log(vars.rate_single_k); 
dy = 2 * (vars.drate_single_k) ./ (vars.rate_single_k);
% Remove those where rate estimate is zero.
indices = find(vars.rate_single_k ~= 0);
x = x(indices);
y = y(indices);
dy = dy(indices);
% plot
errorbar(x, y, dy, 'k.', 'MarkerSize', 10, 'LineWidth', 1.2);

% Adjust axes limits.
oldaxis = axis;
axis([1/(vars.temperature_k(end)) 1/(vars.temperature_k(1)) min(y - dy), max(y + dy)]);

% 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('1/T (K)');
ylabel('ln (k * ps)');

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