close all
fclose('all');
clc

font=16;
tickx=0:45:360;
num=40;
x=0:9:360;
conv=8065.6;
init=zeros(1,num+1);
SOCS1T1=init;

xfilt=(x>40) & (x<330);
xinterp=linspace(40,330,2000);
xSOCS1filt=(xinterp>50) & (xinterp<120);
xCTS2filt=(xinterp>140) & (xinterp<220);
xSOC=xinterp(xSOCS1filt);
xCT=xinterp(xCTS2filt);

%% Hexane

CT=JH;

for j=1:num
    filename = sprintf("C:/Users/Gebruiker/Desktop/Internship/plots/SOC-Pyr-DMA-dihedral-data/SOC-Pyr-DMA-dihedral-data/Hexane/S1/angle(%d)/slurm.out",j);
    fid=fopen(filename);
    C = textscan(fid, '%s');
    C=C{1};
    Index = find(contains(C,'S1:'));
    SOCS1T1(j)=str2double(C{Index(2)+1});
end

thS1=65.7;
thS2=67.1;

yCTeq=interp1(x(xfilt),conv*abs(CT(xfilt)),thS2);
ySOCeq=interp1(x(xfilt),SOCS1T1(xfilt),thS1);

yinterpCT=interp1(x(xfilt),conv*abs(CT(xfilt)),xinterp);
yinterpSOC=interp1(x(xfilt),SOCS1T1(xfilt),xinterp);

[SOCmax, ind]=max(yinterpSOC(xSOCS1filt));
SOCmax;
thSOCmax=xSOC(ind);
[CTmax, ind]=max(yinterpCT(xCTS2filt));
CTmax;
thCTmax=xCT(ind);

% ytable=[yCTeq ySOCeq SOCmax CTmax thSOCmax thCTmax]
% ytable=[yCTeq  CTmax  thCTmax]
yCTeq
CTmax
thCTmax
figure
yyaxis left
CT_plot(x(xfilt),conv^2*abs(CT(xfilt)).^2,font','','-d')
xlim([30 300])
ylim([-0.1 1.075*CTmax^2])
ylabel('$|V_{CT}|^2\:(cm^{-2})$ ','interpreter','latex')
hold on

yyaxis right
SOCplot(x(xfilt),SOCS1T1(xfilt).^2,font,'','-o')
xlim([30 300])
ylim([-0.1 1.075*SOCmax^2])
ylabel('$V_{SOC}^2\:(cm^{-2})$','interpreter','latex')
title('$NHX$','interpreter','latex')
xticks(tickx)
% Hexane
xline(thS1) % S1
xline(thS2) % S2

%% gas

CT=JG;

for j=1:num
    filename = sprintf("C:/Users/Gebruiker/Desktop/Internship/plots/SOC-Pyr-DMA-dihedral-data/SOC-Pyr-DMA-dihedral-data/Gas/S1/angle(%d)/slurm.out",j);
    fid=fopen(filename);
    C = textscan(fid, '%s');
    C=C{1};
    Index = find(contains(C,'S1:'));
    SOCS1T1(j)=str2double(C{Index(2)+1});
end

%Gas
thS1=69.9; % S1
thS2=65.7; % S2

yCTeq=interp1(x(xfilt),conv*abs(CT(xfilt)),thS2);
ySOCeq=interp1(x(xfilt),SOCS1T1(xfilt),thS1);

yinterpCT=interp1(x(xfilt),conv*abs(CT(xfilt)),xinterp);
yinterpSOC=interp1(x(xfilt),SOCS1T1(xfilt),xinterp);

[SOCmax, ind]=max(yinterpSOC(xSOCS1filt));
SOCmax;
thSOCmax=xSOC(ind);
[CTmax, ind]=max(yinterpCT(xCTS2filt));
CTmax;
thCTmax=xCT(ind);


%%
% ytable=[yCTeq ySOCeq SOCmax CTmax thSOCmax thCTmax]
% ytable=[yCTeq  CTmax  thCTmax]
yCTeq
CTmax
thCTmax

figure
yyaxis left
CT_plot(x(xfilt),conv^2*abs(CT(xfilt)).^2,font','','-d')
xlim([30 300])
ylim([-0.1 1.075*CTmax^2])
ylabel('$|V_{CT}|^2\:(cm^{-2})$ ','interpreter','latex')
hold on

yyaxis right
SOCplot(x(xfilt),SOCS1T1(xfilt).^2,font,'','-o')
xlim([30 300])
ylim([-0.1 1.075*SOCmax^2])
ylabel('$V_{SOC}^2\:(cm^{-2})$','interpreter','latex')
title('$Gas$','interpreter','latex')
xticks(tickx)
%Gas
xline(thS1) % S1
xline(thS2) % S2

%% ACN

CT=JA;

for j=1:num
    filename = sprintf("C:/Users/Gebruiker/Desktop/Internship/plots/SOC-Pyr-DMA-dihedral-data/SOC-Pyr-DMA-dihedral-data/ACN/S1/angle(%d)/slurm.out",j);
    fid=fopen(filename);
    C = textscan(fid, '%s');
    C=C{1};
    Index = find(contains(C,'S1:'));
    SOCS1T1(j)=str2double(C{Index(2)+1});
end

thS1=71.6;  % S1
thS2=70.4; % S2

yCTeq=interp1(x(xfilt),conv*abs(CT(xfilt)),thS2);
ySOCeq=interp1(x(xfilt),SOCS1T1(xfilt),thS1);

yinterpCT=interp1(x(xfilt),conv*abs(CT(xfilt)),xinterp);
yinterpSOC=interp1(x(xfilt),SOCS1T1(xfilt),xinterp);

[SOCmax, ind]=max(yinterpSOC(xSOCS1filt));
SOCmax;
thSOCmax=xSOC(ind);
[CTmax, ind]=max(yinterpCT(xCTS2filt));
CTmax;
thCTmax=xCT(ind);

%%
% ytable=[yCTeq ySOCeq SOCmax CTmax thSOCmax thCTmax]
% ytable=[yCTeq  CTmax  thCTmax]
yCTeq
CTmax
thCTmax

figure
yyaxis left
CT_plot(x(xfilt),conv^2*abs(CT(xfilt)).^2,font','','-d')
xlim([30 300])
ylim([-0.1 1.05]*CTmax^2)
ylabel('$|V_{CT}|^2\:(cm^{-2})$ ','interpreter','latex')
hold on

yyaxis right
SOCplot(x(xfilt),SOCS1T1(xfilt).^2,font,'','-o')
xlim([30 300])
ylim([-0.1 1.075*SOCmax^2])
ylabel('$V_{SOC}^2\:(cm^{-2})$','interpreter','latex')
title('$ACN$','interpreter','latex')
xticks(tickx)
%ACN
xline(thS1) % S1
xline(thS2) % S2

