%%%%%%%%%%%%%%%%%%%%%
% radial_function.m %
%%%%%%%%%%%%%%%%%%%%%

clf;
A_ZZ = 2.9296875;
B_ZZ = 2.58787395638939;
C_ZZ = 1.32308467741935;
D_ZZ = 4.3672464262563;
RCUT = 7.0;
K = 0.1;

A_ZC = 2.01006334248919;
B_ZC = 3.13774768336519;
C_ZC = 1.62561120468171;
D_ZC = 3.33014286910615;
W = 3.5;
U = 0.1;

a0 = 4.711804;
r1 = a0 / 2;
r2 = a0 / sqrt(2);
r3 = a0 / 2 * sqrt(3);
orig = 3.1744;

r = 2 : .001 : RCUT-.001;
repulsion_ZZ = exp( A_ZZ * ( B_ZZ - r ) + K ./ ( r - RCUT ) );
attraction_ZZ = - exp( C_ZZ / 2 * ( D_ZZ - r ) + K ./ ( r - RCUT ) );

s = 2 : .001 : W-.001;
repulsion_ZC = exp( A_ZC * ( B_ZC - s ) + U ./ ( s - W ) );
attraction_ZC = - exp( C_ZC * ( D_ZC - s ) + U ./ ( s - W ) );

plot (s, repulsion_ZC, 'r');
axis ([2 4 -6 6]);
hold on;
plot (s, attraction_ZC, 'r--');
plot (r, repulsion_ZZ, '-.');
plot (r, attraction_ZZ, ':');

v = axis;
arrowdraw ([r1,v(3)*0.9], [r1,v(3)], 0.3, 30, 'r');
set(text(r1,v(3)*0.85,'1st'),'Color','red');
arrowdraw ([r2,v(3)*0.9], [r2,v(3)], 0.3, 30, 'b');
set(text(r2,v(3)*0.85,'2nd'),'Color','blue');

arrowdraw ([orig,v(3)*0.9], [orig,v(3)], 0.3, 30, 'k');
set(text(orig*0.975,v(3)*0.85,'hcp Zr'),'Color','black');

legend ('\Phi_{ZrC}(r)', '-h_{ZrC}(r)', ...
        '\Phi_{ZrZr}(r)', '-h_{ZrZr}(r)', 1);

xlabel('r [A]');
ylabel('Energy [eV]');
title(['Radial Functions']);
grid on;

print -depsc radial_function.eps;
