clear all
close all
fclose('all');
clc

font=16;
num=40;
x=0:9:360;

% lim1=0;
% lim2=360;
% ylim1=20;
% tickx=0:45:360;

lim1=50;
lim2=100;
ylim1=5;
tickx=0:5:360;

for txt=["Gas","Hexane","ACN"]
    
    init=zeros(1,num+1);
    ES0=init;
    ES1=init;
    ES2=init;
    ET1=init;
    ET2=init;
    
    %%
    for j=1:num
        for state=["S0","S1","S2","T1","T2"]
            
            filename = sprintf("C:/Users/Gebruiker/Desktop/Internship/plots/SOC-Pyr-DMA-dihedral-data/SOC-Pyr-DMA-dihedral-data/%s/%s/angle(%d)/logfile",txt,state,j);
            fid=fopen(filename);
            C = textscan(fid, '%s');
            C=C{1};
            Index = find(contains(C,'Bond'));
            
            if length(Index)==6
                shift=5;
            else
                shift=2;
            end
            
            En=str2double(C{Index(shift)+2});
            
            filename = sprintf("C:/Users/Gebruiker/Desktop/Internship/plots/SOC-Pyr-DMA-dihedral-data/SOC-Pyr-DMA-dihedral-data/%s/%s/angle(%d)/slurm.out",txt,"S1",j);
            fid2=fopen(filename);
            C = textscan(fid2, '%s');
            C=C{1};
            
            if state=="S0"
                ES0(j)=En;
            elseif state=="S1"
                Index = find(contains(C,'S1:'));
                ES1(j)=En+str2double(C{Index(1)+4});
            elseif state=="S2"
                Index = find(contains(C,'S2:'));
                ES2(j)=En+str2double(C{Index(1)+4});
            elseif state=="T1"
                Index = find(contains(C,'T1:'));
                ET1(j)=En+str2double(C{Index(1)+4});
            elseif state=="T2"
                Index = find(contains(C,'T2:'));
                ET2(j)=En+str2double(C{Index(1)+4});
            end
            
            fclose(fid);
            fclose('all');
        end
    end
    
    ES0(end)=ES0(1);
    ES1(end)=ES1(1);
    ES2(end)=ES2(1);
    ET1(end)=ET1(1);
    ET2(end)=ET2(1);
    %%
    figure
    tit=sprintf("$%s$",txt);
    [~,cal]=findmin(x,ES0);
    energyplot(x,ES0-cal,ES1-cal,ES2-cal,ET1-cal,ET2-cal,font,tit)
    xlim([lim1 lim2])
    ylim([0 ylim1])
    xticks(tickx)
    
    E0p=ES0-cal
    %     [xm, ~] = findmin(x,ES0-cal);
    %     p=polyfit(x(x>40 & x<140),ES0(x>40 & x<140)-cal,2);
    %     p(1)
    %     xline(xm);
    %     [xm, ~] = findmin(x,ES1-cal);
    %     p=polyfit(x(x>40 & x<140),ES1(x>40 & x<140)-cal,2);
    %     p(1)
    %     xline(xm);
    %     [xm, ~] = findmin(x,ET1-cal);
    %     p=polyfit(x(x>40 & x<140),ET1(x>40 & x<140)-cal,2);
    %     p(1)
    %     xline(xm);
    %     [xm, ym] = findmin(x,ET2-cal);
    %     xline(xm);
    %     p=polyfit(x(x>40 & x<140),ET2(x>40 & x<140)-cal,2);
    %     p(1)
    %
    
    legend('$E_{S_{0}}$','$E_{S_{2}}$','$E_{^1 CT}$','$E_{T_{1}}$','$E_{^3 CT}$','interpreter','latex')
end

