function [V,Frab,Frbc,Frac] = lepspoint(drab,drbc,theta,Drab,Drbc,Drac,Brab,Brbc,Brac,lrab,lrbc,lrac,H,grad)

% Gives the Energy and first derivative at any point on LEPS surface.

drac = ((drab.^2)+(drbc.^2)-2*drab*drbc*cos(theta)).^0.5;

%ENERGY
% morse (y) and anti morse (z) functions
ydrab=Drab*(exp(-2*Brab*(drab-lrab))-(2*exp(-Brab*(drab-lrab))));
zdrab=(Drab/2)*(exp(-2*Brab*(drab-lrab))+(2*exp(-Brab*(drab-lrab))));
ydrbc=Drbc*(exp(-2*Brbc*(drbc-lrbc))-(2*exp(-Brbc*(drbc-lrbc))));
zdrbc=(Drbc/2)*(exp(-2*Brbc*(drbc-lrbc))+(2*exp(-Brbc*(drbc-lrbc))));
ydrac=Drac*(exp(-2*Brac*(drac-lrac))-(2*exp(-Brac*(drac-lrac))));
zdrac=(Drac/2)*(exp(-2*Brac*(drac-lrac))+(2*exp(-Brac*(drac-lrac))));

% Coulomb (Q) and Exchange (J) integrals
k=0.18; % In TRIATOMICS this is 0.18 but should be sato parameter?
Qdrab=(ydrab+zdrab+k*(ydrab-zdrab))/2;
Jdrab=(ydrab-zdrab+k*(ydrab+zdrab))/2;
Qdrbc=(ydrbc+zdrbc+k*(ydrbc-zdrbc))/2;
Jdrbc=(ydrbc-zdrbc+k*(ydrbc+zdrbc))/2;
Qdrac=(ydrac+zdrac+k*(ydrac-zdrac))/2;
Jdrac=(ydrac-zdrac+k*(ydrac+zdrac))/2;

% Potential Energy
sq=H.^2;
ist=1; % 1 for ground state, >1 for excited states
FKK=0.5*((Jdrab-Jdrbc).^2+(Jdrbc-Jdrac).^2+(Jdrac-Jdrab).^2); 
V=1/(1+sq)*(Qdrab+Qdrbc+Qdrac-ist*sqrt(FKK));

% FIRST DERIVATIVE
if (grad ~= 0)
    r(1) = drab;
    r(2) = drbc;
    r(3) = drac;
    dd(1) = Drab;
    dd(2) = Drbc;
    dd(3) = Drac;
    fa(1) = Brab;
    fa(2) = Brbc;
    fa(3) = Brac;
    re(1) = lrab;
    re(2) = lrbc;
    re(3) = lrac;
    term = 1/drac;

    for m = 1:3;
       dr(m) = r(m) - re(m); 
       ex1 = exp(-fa(m)*dr(m));
       ex2 = ex1*ex1;

       dmdr(m) = dd(m)*((-2.*fa(m)*ex2)+(2.*fa(m)*ex1));% morse 1st deriv
       dmsdr(m)=dd(m)*((-fa(m)*ex2)-(fa(m)*ex1));% anti-morse 1st deriv

       dQdr(m) = 0.5*((dmdr(m)+dmsdr(m))+0.18*(dmdr(m)-dmsdr(m))); % coulomb 1st deriv
       dJdr(m) = 0.5*((dmdr(m)-dmsdr(m))+0.18*(dmdr(m)+dmsdr(m))); % exchange 1st deriv

    end

    xdQdr1 = dQdr(3)*(r(1)-r(2)*cos(theta))*term; 
    xdQdr2 = dQdr(3)*(r(2)-r(1)*cos(theta))*term;
    xdQdro = dQdr(3)*(r(1)*r(2)*sin(theta))*term; 

    xdJdr1 = dJdr(3)*(r(1)-r(2)*cos(theta))*term; 
    xdJdr2 = dJdr(3)*(r(2)-r(1)*cos(theta))*term;
    xdJdro = dJdr(3)*(r(1)*r(2)*sin(theta))*term;

    xmu = xdJdr1-dJdr(1);
    xnew = (Jdrbc-Jdrac)*(-xdJdr1);
    ymu = dJdr(2)-xdJdr2;
    ynew = (Jdrac-Jdrab)*xdJdr2;

    tm1 = ((Jdrab-Jdrbc)*dJdr(1)+(Jdrac-Jdrab)*xmu+xnew)*0.5*(1/sqrt(FKK));
    tm2 = (-(Jdrab-Jdrbc)*dJdr(2)+(Jdrbc-Jdrac)*ymu+ynew)*0.5*(1/sqrt(FKK));
    tm3 = (-(Jdrbc-Jdrac)*xdJdro+(Jdrac-Jdrab)*xdJdro)*0.5*(1/sqrt(FKK));

    Frab = -1/(1+sq)*(dQdr(1)+xdQdr1-ist*tm1);
    Frbc = -1/(1+sq)*(dQdr(2)+xdQdr2-ist*tm2);
    Frac = -1/(1+sq)*(xdQdro-ist*tm3);
end

