% Plot cumulative distribution function from analyze-distances.py.

clear;

%r = textread('restrained-distances-000');
[r, U, x, y, z] = textread('distances.txt');

ny = 3;
nx = 2;


clf;
subplot(ny,nx,1);
plot(r,'.');

% Compute CDF.
rmax = max(r);
r_sorted = sort(r);
% Normalize so that integral up to 'rmax' is unity.
rnorm = rmax;
Pr = cumsum((1:length(r))');
Pr = Pr / Pr(end);

% Plot.
subplot(ny,nx,2);
plot(r_sorted,Pr,'.');
hold on;

npoints = 10000;
rpoints = linspace(0,rnorm,npoints);
dr = rpoints(2) - rpoints(1);
%Pr_anal = (4/3)*pi*rpoints.^3;
temperature = 298.; % kelvin
kB = 1.381e-23 * 6.022e+23 / 4184.0; % kJ/mol/K
kT = kB * temperature;
beta = 1.0 / kT; % 1 / (kJ/mol)
K = 404.408; % kJ/mol/nm^2
K = K / 10^2; % kJ/mol/A^2
sigma = (beta * K).^(-1/2);
sigma = 5.0; % A
sigma = 3.816371; % A

sigma = 10.0;
r0 = 0.0;

%sigma = 0.2;
%r0 = 3.5;


K = 1.0 / (beta * sigma^2);

dPr_anal = exp(-0.5*((rpoints-r0)/sigma).^2) .* (rpoints.^2) * dr;
Pr_anal = cumsum(dPr_anal);
Pr_anal = Pr_anal / Pr_anal(end);
%Pr_anal = sqrt(pi/2) .* sigma^3 .* erf(rpoints / (sqrt(2) * sigma)) - sigma^2 .* rpoints .* exp(-0.5*(rpoints/sigma).^2);
%Pr_anal = Pr_anal / Pr_anal(end);

plot(rpoints,Pr_anal,'r-');
axis([0 rnorm 0 1]);

% DEBUG
cla;
nbins = 40;
[Nr, rv] = hist(r,nbins);
plot(rv, Nr, 'k.');
hold on;
dPr_anal = exp(-0.5*((rpoints-r0)/sigma).^2) .* (rpoints.^2) * dr;
dPr_anal = dPr_anal / max(dPr_anal) * max(Nr);
plot(rpoints, dPr_anal, 'r-');
axis([0 rmax 0 max(Nr)]);


% Compare U(r).
subplot(ny,nx,3);
plot(r, U, 'k.');
hold on;
Upoints = (K/2) * (rpoints - r0).^2;
plot(rpoints, Upoints, 'r-');

subplot(ny,nx,4);
nbins = 40;
[Nx,x] = hist(x,nbins);
plot(x,Nx,'.');
hold on;
Px = exp(-0.5 * (x/sigma).^2);
Px = Px / max(Px) * max(Nx);
plot(x,Px,'r-');
xlabel('x');
ylabel('p(x)');
set(gca,'YTick',[]);

subplot(ny,nx,5);
[Ny,y] = hist(y,nbins);
plot(y,Ny,'.');
hold on;
Py = exp(-0.5 * (y/sigma).^2);
Py = Py / max(Py) * max(Ny);
plot(y,Py,'r-');
xlabel('y');
ylabel('p(y)');
set(gca,'YTick',[]);

subplot(ny,nx,6);
[Nz,z] = hist(z,nbins);
plot(z,Nz,'.');
hold on;
Pz = exp(-0.5 * (z/sigma).^2);
Pz = Pz / max(Pz) * max(Nz);
plot(z,Pz,'r-');
xlabel('z');
ylabel('p(z)');
set(gca,'YTick',[]);
