function Vmat = plotlepsgui(a,b,c,theta,xrabi,xrbci,prabi,prbci,H,listno,cutoff,stepno,mep)

% Produces output from GUI input.

% get parameters 
[ma,mb,mc,Drab,Drbc,Drac,lrab,lrbc,lrac,Brab,Brbc,Brac,mina,maxa,minb,maxb,cont,atma,atmb,atmc] = params(a,b,c);

% Reduced masses.
mab = (ma*mb)/(ma+mb);
mbc = (mb*mc)/(mb+mc);
mac = (ma*mc)/(ma+mc);

% Set gap between points on surface. If surf slow change top val. TRIATOMICS was 0.15 anyway. 
if listno == 2;
    resl = 0.02;
else
    resl = 0.02;
end

% Always keep as zero here: Stops gradients being determined in surface loop
grad = 0;

% Calculate Surface

drbccount = 0;
drabcount = 0;

h = waitbar(0,'Calculating Surface...');

for drab = mina:resl:maxa
    drabcount = drabcount + 1;
    for drbc = minb:resl:maxb
        drbccount = drbccount + 1;

        V = lepspoint(drab,drbc,theta,Drab,Drbc,Drac,Brab,Brbc,Brac,lrab,lrbc,lrac,H,grad);
        
        if V < cutoff
          Vmat(drbccount, drabcount) = V;
        else
          Vmat(drbccount, drabcount) = cutoff;
        end
        
    end
    drbccount = 0;
    waitbar(drab/maxa)
end
close(h)

% Calculate Trajectory 
itlimit = stepno;
dt = 0.005;
s = 5E-4;
% Define Step method:
%   0 = Euler step with numerical gradients
%   1 = Euler step with analytical gradients
%   2 = Gradients with Hessians (Helgaker et al, Chem. Phys. Lett., 2, 3,
%   1990)
grad = 2;

% Initial values derivied from GUI input
vrabi = prabi/mab;
vrbci = prbci/mbc;
xraci = ((xrabi.^2)+(xrbci.^2)-2*xrabi*xrbci*cos(theta)).^0.5;
vraci=0;
% if (vrabi+vrbci < 0)
%    vraci = -((vrabi.^2)+(vrbci.^2)-2*vrabi*vrbci*cos(theta)).^0.5;
% else
%    vraci = ((vrabi.^2)+(vrbci.^2)-2*vrabi*vrbci*cos(theta)).^0.5;
% end
thetai = theta;
ti = 0;

xrab(1,:) = xrabi;
xrbc(1,:) = xrbci;
xrac(1,:) = xraci;
vrab(1,:) = vrabi;
vrbc(1,:) = vrbci;
vrac(1,:) = vraci;
theta(1,:) = theta;
t(1,:) = ti;

h = waitbar(0,'Calculating Trajectory...');

for itcounter = 1:itlimit;
    if mep == true;
        vrabi = 0;
        vrbci = 0;
        vraci = 0;
    end
    % GET FORCES
    if grad == 0
        % NUMERICAL GRADIENTS (grad=0)
        Vrint(itcounter,:) = lepspoint(xrabi,xrbci,thetai,Drab,Drbc,Drac,Brab,Brbc,Brac,lrab,lrbc,lrac,H,grad);
        Vrab(itcounter,:) = lepspoint(xrabi + s,xrbci,thetai,Drab,Drbc,Drac,Brab,Brbc,Brac,lrab,lrbc,lrac,H,grad);
        Vrbc(itcounter,:) = lepspoint(xrabi,xrbci - s,thetai,Drab,Drbc,Drac,Brab,Brbc,Brac,lrab,lrbc,lrac,H,grad);
        Frab(itcounter,:) = (Vrab(itcounter,:) - Vrint(itcounter,:))/((xrabi - s) - xrabi);
        Frbc(itcounter,:) = (Vrbc(itcounter,:) - Vrint(itcounter,:))/((xrbci + s) - xrbci);
        Frac(itcounter,:) = ((Frab(itcounter,:).^2)+(Frbc(itcounter,:).^2)-2*Frab(itcounter,:)*Frbc(itcounter,:)*cos(thetai)).^0.5;
        % EULER STEP
        [xrabf,xrbcf,xracf,thetaf,vrabf,vrbcf,vracf,arab(itcounter,:),arbc(itcounter,:),arac(itcounter,:),Ktot(itcounter,:),tf]=eulerstep(xrabi,xrbci,xraci,thetai,vrabi,vrbci,vraci,ti,dt,Frab(itcounter,:),Frbc(itcounter,:),Frac(itcounter,:),mab,mbc,mac);
    elseif grad == 1
        % ANALYTICAL GRADIENTS (grad=1)
        [Vrint(itcounter,:),Frab(itcounter,:),Frbc(itcounter,:),Frac(itcounter,:)] = lepspoint(xrabi,xrbci,thetai,Drab,Drbc,Drac,Brab,Brbc,Brac,lrab,lrbc,lrac,H,grad);
        % EULER STEPvracf = -((vrabf.^2)+(vrbcf.^2)-2*vrabf*vrbcf*cos(theta + XX(3))).^0.5;
        [xrabf,xrbcf,xracf,thetaf,vrabf,vrbcf,vracf,arab(itcounter,:),arbc(itcounter,:),arac(itcounter,:),Ktot(itcounter,:),tf]=eulerstep(xrabi,xrbci,xraci,thetai,vrabi,vrbci,vraci,ti,dt,Frab(itcounter,:),Frbc(itcounter,:),Frac(itcounter,:),mab,mbc,mac);
    else
        % GRADIENTS WITH HESSIANS (grad=2)
        [Vrint(itcounter,:),Frab(itcounter,:),Frbc(itcounter,:),Frac(itcounter,:),hr1r1(itcounter,:),hr1r2(itcounter,:),hr1r3(itcounter,:),hr2r2(itcounter,:),hr2r3(itcounter,:),hr3r3(itcounter,:)] = lepspoint(xrabi,xrbci,thetai,Drab,Drbc,Drac,Brab,Brbc,Brac,lrab,lrbc,lrac,H,grad);
        % GRADIENTS WITH HESSIANS STEP
        [xrabf,xrbcf,xracf,thetaf,vrabf,vrbcf,vracf,tf,arab(itcounter,:),arbc(itcounter,:),arac(itcounter,:),Ktot(itcounter,:)] = lepnorm(xrabi,xrbci,thetai,Frab(itcounter,:),Frbc(itcounter,:),Frac(itcounter,:),vrabi,vrbci,vraci,hr1r1(itcounter,:),hr1r2(itcounter,:),hr1r3(itcounter,:),hr2r2(itcounter,:),hr2r3(itcounter,:),hr3r3(itcounter,:),ma,mb,mc,ti,dt,mep);
    end
    
    etot(itcounter,:) = Vrint(itcounter,:)+Ktot(itcounter,:);
    
    if mep == true;
        vrabf = 0;
        vrbcf = 0;
        vracf = 0;
    end
        
    if itcounter ~= itlimit
        xrab(itcounter+1,:) = xrabf;
        xrbc(itcounter+1,:) = xrbcf;
        xrac(itcounter+1,:) = xracf;
        vrab(itcounter+1,:) = vrabf;
        vrbc(itcounter+1,:) = vrbcf;
        vrac(itcounter+1,:) = vracf;
        theta(itcounter+1,:) = thetaf;
        t(itcounter+1,:) = tf;
    
    end

    xrabi = xrabf;
    xrbci = xrbcf;
    xraci = xracf;
    vrabi = vrabf;
    vrbci = vrbcf;
    vraci = vracf;
    thetai = thetaf;
    ti = tf;
    
    print = thetai/pi;
    
    waitbar(itcounter/itlimit)
end
close(h)

% OUTPUT
if listno == 1;
    % Contour plot
    contour(mina:resl:maxa, minb:resl:maxb, Vmat, cont); axis tight
    xlabel('A-B bond distance (\AA)','interpreter','latex');
    ylabel('B-C bond distance (\AA)','interpreter','latex');
    hold on;
    plot(xrab,xrbc,'r');
    hold off;
    axis([mina maxa minb maxb])
elseif listno == 2;
    % Surface plot
    set(gcf,'renderer','opengl');
    surfc(mina:resl:maxa, minb:resl:maxb, Vmat, 'FaceAlpha',0.7);
    hold on
    plot3(xrab,xrbc,Vrint+0.5,'-k','Linewidth',4)%,'Markersize',8)
    hold off
    axis tight;
    shading interp;
    view([135,30]);
    xlabel('A-B bond distance (\AA)','interpreter','latex');
    ylabel('B-C bond distance (\AA)','interpreter','latex');
    zlabel('Potential Energy (Kcal/mol)','interpreter','latex');
    axis([mina maxa minb maxb])
elseif listno == 3;
    % Internuclear distance vs time
    plot(t,xrab,'-b');
    xlabel('time');
    ylabel('internuclear distance (Angstrom)');
    hold on;
    plot(t,xrbc,'-r');
    plot(t,xrac,'-g');
    h = legend('A-B bond distance (Ang)','B-C bond distance (Ang)','A-C bond distance (Ang)');
    hold off;
elseif listno == 4;
    % Potential energy vs time
    plot(t,Vrint);
    xlabel('time');
    ylabel('Potential Energy (kcal/mol)');
elseif listno == 5;
    % Kinetic energy vs time
    plot(t,Ktot,'-r');
    xlabel('time');
    ylabel('Kinetic Energy (kcal/mol)');
elseif listno == 6;
    % Internuclear velocities vs internuclear distance
    h1 = plot(vrab(1),xrab(1),'ob','Erasemode','none');
    hold on;
    h2 = plot(vrbc(1),xrbc(1),'dr','Erasemode','none');
    h3 = plot(vrac(1),xrac(1),'sg','Erasemode','none');
    hold off;
    axis([-8,8,0,4.5]); % axis required for unfixed plots but how best to do this?
    xlabel('internuclear velocity (m/s)');
    ylabel('internuclear distance (m)');
    h = legend('A-B','B-C','A-C');
    nframes = 1000;
    for frame = 1:nframes;
        vrabt = vrab(frame,:);
        xrabt = xrab(frame,:);
        vrbct = vrbc(frame,:);
        xrbct = xrbc(frame,:);
        vract = vrac(frame,:);
        xract = xrac(frame,:);
        set(h1,'XData',vrabt,'YData',xrabt);
        set(h2,'XData',vrbct,'YData',xrbct);
        set(h3,'XData',vract,'YData',xract);
        drawnow
        pause(0.05)
    end
elseif listno == 7;
    % Internuclear acceleration vs internuclear distance
    h1 = plot(arab(1),xrab(1),'ob','Erasemode','none');
    hold on;
    h2 = plot(arbc(1),xrbc(1),'dr','Erasemode','none');
    h3 = plot(arac(1),xrac(1),'sg','Erasemode','none');
    hold off;
    axis([-60,60,0,4.5]); % axis required for unfixed plots but how best to do this?
    xlabel('internuclear acceleration (m/s^2)');
    ylabel('internuclear distance (m)');
    h = legend('A-B','B-C','A-C');
    nframes = 1000;
    for frame = 1:nframes;
        arabt = arab(frame,:);
        xrabt = xrab(frame,:);
        arbct = arbc(frame,:);
        xrbct = xrbc(frame,:);
        aract = arac(frame,:);
        xract = xrac(frame,:);
        set(h1,'XData',arabt,'YData',xrabt);
        set(h2,'XData',arbct,'YData',xrbct);
        set(h3,'XData',aract,'YData',xract);
        drawnow
        pause(0.05)
    end
    elseif listno == 8;
    % Internuclear momenta vs internuclear distance
    h1 = plot(vrab(1)*mab,xrab(1),'ob','Erasemode','none');
    hold on;
    h2 = plot(vrbc(1)*mbc,xrbc(1),'dr','Erasemode','none');
    h3 = plot(vrac(1)*mac,xrac(1),'sg','Erasemode','none');
    hold off;
    axis([-8,8,0,4.5]); % axis required for unfixed plots but how best to do this?
    xlabel('internuclear momentum (N s)');
    ylabel('internuclear distance (m)');
    h = legend('A-B','B-C','A-C');
    nframes = 1000;
    for frame = 1:nframes;
        prabt = vrab(frame,:)*mab;
        xrabt = xrab(frame,:);
        prbct = vrbc(frame,:)*mbc;
        xrbct = xrbc(frame,:);
        pract = vrac(frame,:)*mac;
        xract = xrac(frame,:);
        set(h1,'XData',prabt,'YData',xrabt);
        set(h2,'XData',prbct,'YData',xrbct);
        set(h3,'XData',pract,'YData',xract);
        drawnow
        pause(0.05)
    end
elseif listno == 9;
    % Internuclear velocities vs time 
    plot(t,vrab,'-b');
    xlabel('time');
    ylabel('internuclear velocity');
    hold on;
    plot(t,vrbc,'-r');
    plot(t,vrac,'-g');
    h = legend('A-B internuclear velocity','B-C internuclear velocity','A-C internuclear velocity');
    hold off;
elseif listno == 10;
    % Internuclear acceleration vs time
    plot(t,arab,'-b');
    xlabel('time (s)');
    ylabel('internuclear acceleration (m/s^2)');
    hold on;
    plot(t,arbc,'-r');
    plot(t,arac,'-g');
    h = legend('A-B internuclear acceleration (m/s^2)','B-C internuclear acceleration (m/s^2)','A-C internuclear acceleration (m/s^2)');
    hold off;
elseif listno == 11;
    % Internuclear momenta vs time 
    plot(t,vrab*mab,'-b');
    xlabel('time');
    ylabel('internuclear momentum');
    hold on;
    plot(t,vrbc*mbc,'-r');
    plot(t,vrac*mac,'-g');
    h = legend('A-B internuclear momentum','B-C internuclear momentum','A-C internuclear momentum');
    hold off;
elseif listno == 12; 
    % animation
    % problem with the distance display when it is >= 10: formatting?
    % can prevent problem by adjusting number of steps.
    nframes = itlimit;
    for frame = 1:nframes;
        axpos(frame,:) = xrab(frame,:);
        aypos(frame,:) = 0;
        bxpos(frame,:) = 0;
        bypos(frame,:) = 0;
        cypos(frame,:) = sin(theta) * xrbc(frame,:);
        cxpos(frame,:) = cos(theta) * xrbc(frame,:);
        strab(frame,:) = ['A-B = ',num2str(xrab(frame,:),'%-7.4f')];
        strbc(frame,:) = ['B-C = ',num2str(xrbc(frame,:),'%-7.4f')];
        strac(frame,:) = ['A-C = ',num2str(xrac(frame,:),'%-7.4f')];
    end
    h1 = image([axpos(1)-0.25,axpos(1)+0.25],[aypos(1)-0.25,aypos(1)+0.25],atma,'Erasemode','xor');
    hold on
    h2 = image([bxpos(1)-0.25,bxpos(1)+0.25],[bypos(1)-0.25,bypos(1)+0.25],atmb,'Erasemode','xor');
    h3 = image([cxpos(1)-0.25,cxpos(1)+0.25],[cypos(1)-0.25,cypos(1)+0.25],atmc,'Erasemode','xor');
    h4 = text(2.0,-2.8,strab(1,:),'Edgecolor','blue','Erasemode','xor');
    h5 = text(2.0,-2.5,strbc(1,:),'Edgecolor','red','Erasemode','xor');
    h6 = text(2.0,-2.2,strac(1,:),'Edgecolor','green','Erasemode','xor');
    hold off
    h = legend('A','B','C');
    axis([-3.5,3.5,-3.5,3.5])
    xlabel('distance (\AA)','interpreter','Latex');
    ylabel('distance (\AA)','interpreter','Latex');
    while 1;
    for frame = 1:nframes; 
        axpost = axpos(frame,:);
        aypost = aypos(frame,:);
        bxpost = bxpos(frame,:);
        bypost = bypos(frame,:);
        cxpost = cxpos(frame,:);
        cypost = cypos(frame,:);
        strabt = strab(frame,:);
        strbct = strbc(frame,:);
        stract = strac(frame,:);
        set(h1,'XData',[axpost-0.25,axpost+0.25],'YData',[aypost-0.25,aypost+0.25]);
        set(h2,'XData',[bxpost-0.25,bxpost+0.25],'YData',[bypost-0.25,bypost+0.25]);
        set(h3,'XData',[cxpost-0.25,cxpost+0.25],'YData',[cypost-0.25,cypost+0.25]);
        set(h4,'string',strabt);
        set(h5,'string',strbct);
        set(h6,'string',stract);
        drawnow
        pause(0.05)
    end
    end
end