function Vmat = plotlepsgui(a,b,c,theta,xrabi,xrbci,vrabi,vrbci,H,listno,cutoff,stepno)

% 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/3.5)
end
close(h)

% Calculate Trajectory 
itlimit = stepno;
dt = 0.005;
s = 5E-4;
% Define Gradient method: 0 is numerical, 1 is analytical
grad = 1;

% Initial values derivied from GUI input
xraci = ((xrabi.^2)+(xrbci.^2)-2*xrabi*xrbci*cos(theta)).^0.5;
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
ti = 0;

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

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

for itcounter = 1:itlimit;
    
    % GET FORCES
    if grad == 0
        % NUMERICAL GRADIENTS (grad=0)
        Vrint(itcounter,:) = lepspoint(xrabi,xrbci,theta,Drab,Drbc,Drac,Brab,Brbc,Brac,lrab,lrbc,lrac,H,grad);
        Vrab(itcounter,:) = lepspoint(xrabi + s,xrbci,theta,Drab,Drbc,Drac,Brab,Brbc,Brac,lrab,lrbc,lrac,H,grad);
        Vrbc(itcounter,:) = lepspoint(xrabi,xrbci - s,theta,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(theta)).^0.5;
    else
        % ANALYTICAL GRADIENTS (grad=1)
        [Vrint(itcounter,:),Frab(itcounter,:),Frbc(itcounter,:),Frac(itcounter,:)] = lepspoint(xrabi,xrbci,theta,Drab,Drbc,Drac,Brab,Brbc,Brac,lrab,lrbc,lrac,H,grad);
    end

    % EULER STEP: maybe extend to leapfrog or Runge-kutta.
    arab(itcounter,:) = (Frab(itcounter,:))/mab;
    arbc(itcounter,:) = (Frbc(itcounter,:))/mbc;
    if (arab(itcounter,:)+arbc(itcounter,:) < 0)
        arac(itcounter,:) = -((arab(itcounter,:).^2)+(arbc(itcounter,:).^2)-2*arab(itcounter,:)*arbc(itcounter,:)*cos(theta)).^0.5;%(Frac(itcounter,:))/mac;
    else
        arac(itcounter,:) = ((arab(itcounter,:).^2)+(arbc(itcounter,:).^2)-2*arab(itcounter,:)*arbc(itcounter,:)*cos(theta)).^0.5;
    end
    

    Krab(itcounter,:) = 0.5*mab*vrabi.^2;
    Krbc(itcounter,:) = 0.5*mbc*vrbci.^2;
    Krac(itcounter,:) = 0.5*mac*vraci.^2;
    Ktot(itcounter,:) = Krab(itcounter,:)+Krbc(itcounter,:);
    % energy not quite consistent with large oscilations.
%    Etot(itcounter,:) = Ktot(itcounter,:) + Vrint(itcounter,:);
        
    vrabf = vrabi + (arab(itcounter,:)*dt);
    vrbcf = vrbci + (arbc(itcounter,:)*dt);
    vracf = vraci + (arac(itcounter,:)*dt);

    xrabf = xrabi + (vrabf*dt);
    xrbcf = xrbci + (vrbcf*dt);
    xracf = xraci + (vracf*dt);

    tf = ti + dt;

    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;
        t(itcounter+1,:) = tf;
    
    end

    xrabi = xrabf;
    xrbci = xrbcf;
    xraci = xracf;
    vrabi = vrabf;
    vrbci = vrbcf;
    vraci = vracf;
    ti = tf;
    
    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;
    % Total energy vs time
    plot(t,Vrint,'-b');
    xlabel('time');
    ylabel('Energy (kcal/mol)');
    hold on;
    plot(t,Ktot,'-r');
%    plot(t,Etot,'-g');
    h = legend('potential energy (kcal/mol)','kinetic energy (kcal/mol)');%,'total energy (kcal/mol)');
    hold off;
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 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 == 9;
    % 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 == 10; 
    % 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