% building potential function
% same as mytest01 but now as a function with subfunctions
% Comparison of Lenard-Jones with Morse potential functions
% Milos Steinhart, Milky Way 
function pot
x=3.5:0.01:8;
% depth of the well
epsi=100;
% distance of the minimum of the potential
rmin=4;
% Lenard-Jones
%yLJ=epsi*((rmin./x).^12-2*(rmin./x).^6);
yLJ=local_lj(epsi,rmin,x);
plot(x,yLJ,'b')
hold on
% Morse
% adjustable parameter allegedly aa*rmin determines range of the forces;
% higher aa means shorter interactions !!!
x=3.4:0.01:8;
aa=1.3; aa=1.5;
yM=local_morse(aa,epsi,rmin,x);
plot(x,yM,'m.')
x0=aa*rmin;
yM0=local_morse(aa,epsi,rmin,x0);
%plot(x0,yM0,'go');
x=3.2:0.1:4.8;
y2=epsi*(2*(x-rmin).^2-1);
plot(x,y2,'g-')
legend('Lenard-Jones','Morse','Parabola');
xlabel('Radius [A]');
ylabel('Potential [au]');
title(['L-J vs Morse Potential (r_0 = 4, \epsilon = 100, a = 1.5)'])


function res=local_lj(epsi,rmin,x)
res=epsi*((rmin./x).^12-2*(rmin./x).^6);

function res=local_morse1(aa,epsi,rmin,x)
res=epsi*((1-exp(-aa*(x-rmin))).^2-1);

function res=local_morse(aa,epsi,rmin,x)
res=epsi*exp(aa*(rmin-x)).*(exp(aa*(rmin-x))-2);
