clear;

nbins = 30; 

[E_t_complex] = textread('efield-complex.out');
[E_t_ligand] = textread('efield-ligand.out');

N_complex = length(E_t_complex);
N_ligand = length(E_t_ligand);

g_complex = statistical_inefficiency(E_t_complex)
g_ligand = statistical_inefficiency(E_t_ligand)

% Report average frequencies.
stark_tuning_rate = 0.87; % cm^{-1} / (MV/cm)
E_complex = mean(E_t_complex);
dE_complex = std(E_t_complex) / sqrt(N_complex/g_complex);
E_ligand = mean(E_t_ligand);
dE_ligand = std(E_t_ligand) / sqrt(N_ligand/g_ligand);
disp(sprintf('SOLVENT:  <E> =  %.1f +- %.1f MV/cm | Stark shift %.1f +- %.1f cm^{-1}', E_ligand, dE_ligand, E_ligand*stark_tuning_rate, dE_ligand*stark_tuning_rate));
disp(sprintf('COMPLEX:  <E> =  %.1f +- %.1f MV/cm | Stark shift %.1f +- %.1f cm^{-1}', E_complex, dE_complex, E_complex*stark_tuning_rate, dE_complex*stark_tuning_rate));
disp(sprintf('delta     <E> = %+.1f +- %.1f MV/cm | Stark shift %+.1f +- %.1f cm^{-1}', (E_complex-E_ligand), sqrt(dE_ligand^2 + dE_complex^2), (E_complex-E_ligand)*stark_tuning_rate, sqrt(dE_ligand^2+dE_complex^2)*stark_tuning_rate));

% PLOT
clf;

[Nx_complex,x_complex] = hist(E_t_complex, nbins); 
[Nx_ligand,x_ligand] = hist(E_t_ligand, nbins); 

Px_complex = Nx_complex / N_complex;
Px_ligand = Nx_ligand / N_ligand;

dPx_complex = sqrt(Px_complex .* (1-Px_complex)) / sqrt(N_complex/g_complex);
dPx_ligand = sqrt(Px_ligand .* (1-Px_ligand)) / sqrt(N_ligand/g_ligand);


hold on;

%errorbar(x_ligand, Px_ligand, dPx_ligand, 'k.-');
%errorbar(x_complex, Px_complex, dPx_complex, 'r.-');

fill([x_ligand fliplr(x_ligand)], [Px_ligand+dPx_ligand fliplr(Px_ligand-dPx_ligand)], [.5 .5 .5], 'linestyle', 'none');
fill([x_complex fliplr(x_complex)], [Px_complex+dPx_complex fliplr(Px_complex-dPx_complex)], [1 0.5 0.5], 'linestyle', 'none');

plot(x_ligand, Px_ligand, 'k-', 'linewidth', 2);
plot(x_complex, Px_complex, 'r-', 'linewidth', 2);

legend('bosutinib in solvent', 'bosutinib in complex')
xlabel('E (MV/cm)');
ylabel('P(E)');

Emin = min([E_t_complex; E_t_ligand]);
Emax = max([E_t_complex; E_t_ligand]);
Pmin = 0.0;
Pmax = max([(Px_complex+dPx_complex) (Px_ligand+dPx_ligand)]);
axis([Emin Emax Pmin Pmax]);

exportfig(gcf, 'src-bosutinib-efield.eps', 'width', 6, 'height', 4, 'color', 'cmyk');
system('epstopdf src-bosutinib-efield.eps');

