function Cell_Coverage_Comparison

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Authors:
% Turner, Leo; Burbanks, Andrew; Cerasuolo, Marianna
%
% Affiliation:
% University of Portsmouth, UK
%
% Date:
% 13/05/2022
%
% Purpose:
% Create the figures comparing with TRAMP mice data.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Produce graphs to show comparison between TRAMP Mice data and
% Hybrid model simulations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

folder='./data/'
batch_names = [
    "20220906_0000",
    "20220906_0001",
    "20220906_0002", 
    "20220906_0003", 
    "20220906_0004",
    "20220906_0005",
    "20220906_0006",
    "20220906_0007",
    "20220906_0008",
    "20220906_0009",
    "20220906_0010",
    "20220906_0011",
    "20220906_0012", 
    "20220906_0013", 
    "20220906_0014",
    "20220906_0015",
    "20220906_0016",
    "20220906_0017",
    "20220906_0018",
    "20220906_0019",
    "20220906_0020",
    "20220906_0021",
    "20220906_0022", 
    "20220906_0023", 
    "20220906_0024",
    "20220906_0025",
    "20220906_0026",
    "20220906_0027",
    "20220906_0028",
    "20220906_0029",
    "20220906_0030",
    "20220906_0031",
    "20220906_0032", 
    "20220906_0033", 
    "20220906_0034",
    "20220906_0035",
    "20220906_0036",
    "20220906_0037",
    "20220906_0038",
    "20220906_0039",
    "20220906_0040",
    "20220906_0041",
    "20220906_0042", 
    "20220906_0043", 
    "20220906_0044",
    "20220906_0045",
    "20220906_0046",
    "20220906_0047",
    "20220906_0048",
    "20220906_0049"
    ]
dummy = [
];
batch_name = batch_names(1)

batchname = "20220906_0000"

% We make a comparison with the *18 week* therapies
schedule = "5_18_0"

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define Data Arrays
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load("./data/"+batchname+"_"+schedule+"_none.mat", 'Cell_Daily', 'T', 'T_start')

treatment_name = 'none'
TIME=1:T;
    %ADB start: handle multiple batches
    all_S_Daily = zeros(length(batch_names), length(TIME));
    all_R_Daily = zeros(length(batch_names), length(TIME));

    %ADB collect all the data on R and S
    original_batch_name = batch_name
    for i_batch=1:length(batch_names)
        single_batch_name = batch_names(i_batch)
        full_name = folder+single_batch_name+"_"+schedule+"_"+treatment_name;
        load(full_name+".mat");
            for k=1:T
        S_Daily(:,:,k)= (Cell_Daily(:,:,k)==1) + (Cell_Daily(:,:,k)==2);
        R_Daily(:,:,k)= (Cell_Daily(:,:,k)==3) + (Cell_Daily(:,:,k)==4);
            end
    length(batch_names)
    for k=1:T
        Total_S_Daily(k)=sum(sum(S_Daily(:,:,k)));
        Total_R_Daily(k)=sum(sum(R_Daily(:,:,k)));
    end
        all_S_Daily(i_batch,:)=Total_S_Daily;
        all_R_Daily(i_batch,:)=Total_R_Daily;
    end
    all_Daily_none = (all_S_Daily+all_R_Daily)/100;

load("./data/"+batchname+"_"+schedule+"_none.mat", 'Cell_Daily', 'T', 'T_start')

No_Therapy_6={34.7615,18.044,15.7451,32.6845,14.3223,10.0442,19.764,17.028,0,21.737,8.88654,19.6539,25.5023,19.135,23.3422,10.8765,16.1096,6.39154};
No_Therapy_6=cell2mat(No_Therapy_6);
No_Therapy_6_Time=zeros(1,length(No_Therapy_6))+(T_start);
No_Therapy_12={11.3004,13.285,87.1934,46.971,39.7019,69.9883,28.4091,23.206,6.76697,23.5508,27.172,0,15.6238,29.2909};
No_Therapy_12=cell2mat(No_Therapy_12);
No_Therapy_12_Time=zeros(1,length(No_Therapy_12))+(T_start+6*7);
No_Therapy_16={27.944,12.9866,18.246,35.51,36.0268,24.2356,14.6927,22.7556,35.991,84.9913,74.5916,46.001,23.9941,88.9205,61.5145,60.7755,30.0211,48.097,61.7402,21.6486,41.3649,35.8297,17.4095,57.4101,62.5677,15.5866,21.1429,22.139,22.1377,22.8326};
No_Therapy_16=cell2mat(No_Therapy_16);
No_Therapy_16_Time=zeros(1,length(No_Therapy_16))+(T_start+10*7);
No_Therapy_19={75.4369,70.5194,45.8653,60.407,83.1543,66.1963,76.3264,68.2305,72.6618,33.6864,35.6134,40.7182};
No_Therapy_19=cell2mat(No_Therapy_19);
No_Therapy_19_Time=zeros(1,length(No_Therapy_19))+(T_start+13*7);
No_Therapy_24={61.097,39.9855,59.5511,78.9933,69.4136,43.2717,49.7097,58.2118,44.1861,48.4874,79.5777,69.4906,95.0804,42.7433,56.7227,45.2716,41.5933};
No_Therapy_24=cell2mat(No_Therapy_24);
No_Therapy_24_Time=zeros(1,length(No_Therapy_24))+(T_start+18*7-1);

%No_Therapy_Means={mean(No_Therapy_6),mean(No_Therapy_12),mean(No_Therapy_16),mean(No_Therapy_19),mean(No_Therapy_24)};

%load('Hybrid_no_Enza_treat_24_0322.mat', 'Cell_Daily', 'T')

%CELL VOLUME OVER TIME
for k=1:T
    Total_Cell_Coverage_No_Therapy(k)= (sum(sum(((Cell_Daily(:,:,k)==1) + (Cell_Daily(:,:,k)==2)+(Cell_Daily(:,:,k)==3) + (Cell_Daily(:,:,k)==4))/10)))/10;
end

TIME=1:T;

treatment_name = 'enza';
    %ADB collect all the data on R and S
    original_batch_name = batch_name
    for i_batch=1:length(batch_names)
        single_batch_name = batch_names(i_batch)
        full_name = folder+single_batch_name+"_"+schedule+"_"+treatment_name;
        load(full_name+".mat");
            for k=1:T
        S_Daily(:,:,k)= (Cell_Daily(:,:,k)==1) + (Cell_Daily(:,:,k)==2);
        R_Daily(:,:,k)= (Cell_Daily(:,:,k)==3) + (Cell_Daily(:,:,k)==4);
            end
    length(batch_names)
    for k=1:T
        Total_S_Daily(k)=sum(sum(S_Daily(:,:,k)));
        Total_R_Daily(k)=sum(sum(R_Daily(:,:,k)));
    end
        all_S_Daily(i_batch,:)=Total_S_Daily;
        all_R_Daily(i_batch,:)=Total_R_Daily;
    end
    all_Daily_enza = (all_S_Daily+all_R_Daily)/100;

Enza_6={34.7615,18.0444,15.7451,32.6845,14.3223,10.0442,19.764,17.028,0,21.737,8.88654,19.6539,25.5023,19.135,23.3422,10.8765,16.1096,6.39154};
Enza_6=cell2mat(Enza_6);
Enza_6_Time=zeros(1,length(Enza_6))+(T_start);
Enza_12={31.3458,30.5736,30.1972,35.6338,9.84163,10.6647,23.9095,23.0273,20.5979,15.5657,1.49304,1.22729,42.3503,31.6328,19.8572,10.8212,0.830235,8.01625,7.69433,15.4685,19.0499,20.098};
Enza_12=cell2mat(Enza_12);
Enza_12_Time=zeros(1,length(Enza_12))+(T_start+6*7);
Enza_16={0,0.980107,61.4794,40.7001,4.53618,2.65753,8.69929,28.3359,32.8909,28.8386,32.8938,65.1652,6.07746,1.05155,6.43351,13.5993,22.821,0,0.260299,8.50278,10.9315};
Enza_16=cell2mat(Enza_16);
Enza_16_Time=zeros(1,length(Enza_16))+(T_start+10*7);
Enza_22={59.4857,59.0634,42.7101,38.4965,17.1829,38.0764,19.6511,21.1386,56.1201,35.5047,19.7523,40.6311};
Enza_22=cell2mat(Enza_22);
Enza_22_Time=zeros(1,length(Enza_22))+(T_start+13*7);
Enza_24={22.551,30.556,24.9437,46.9794,25.5014,39.7469,57.5567,23.0782,32.5606,40.5471,58.8412,52.8805,27.5649,29.6446};
Enza_24=cell2mat(Enza_24);
Enza_24_Time=zeros(1,length(Enza_24))+(T_start+18*7-1);

%load('Hybrid_Enza_treat_24_1221_0.0017.mat', 'Cell_Daily', 'T')
%load('Hybrid_Enza_treat_24_0322.mat', 'Cell_Daily', 'T')
load("./data/"+batchname+"_"+schedule+"_enza.mat")

%CELL VOLUME OVER TIME
for k=1:T
    Total_Cell_Coverage_Enza(k)= (sum(sum(((Cell_Daily(:,:,k)==1) + (Cell_Daily(:,:,k)==2)+(Cell_Daily(:,:,k)==3) + (Cell_Daily(:,:,k)==4))/10)))/10;
end

TIMEt=1:T;

treatment_name = 'caba';
    %ADB collect all the data on R and S
    original_batch_name = batch_name
    for i_batch=1:length(batch_names)
        single_batch_name = batch_names(i_batch)
        full_name = folder+single_batch_name+"_"+schedule+"_"+treatment_name;
        load(full_name+".mat");
            for k=1:T
        S_Daily(:,:,k)= (Cell_Daily(:,:,k)==1) + (Cell_Daily(:,:,k)==2);
        R_Daily(:,:,k)= (Cell_Daily(:,:,k)==3) + (Cell_Daily(:,:,k)==4);
            end
    length(batch_names)
    for k=1:T
        Total_S_Daily(k)=sum(sum(S_Daily(:,:,k)));
        Total_R_Daily(k)=sum(sum(R_Daily(:,:,k)));
    end
        all_S_Daily(i_batch,:)=Total_S_Daily;
        all_R_Daily(i_batch,:)=Total_R_Daily;
    end
    all_Daily_caba = (all_S_Daily+all_R_Daily)/100;
%Enza_Means={mean(Enza_6),mean(Enza_12),mean(Enza_16),mean(Enza_22),mean(Enza_24),};

Caba_6={34.7615,18.0444,15.7451,32.6845,14.3223,10.0442,19.764,17.028,0,21.737,8.88654,19.6539,25.5023,19.135,23.3422,10.8765,16.1096,6.39154};
Caba_6=cell2mat(Caba_6);
Caba_6_Time=zeros(1,length(Caba_6))+(T_start);
Caba_24={52,23,21.5,8.2,3.7,2.5};
Caba_24=cell2mat(Caba_24);
Caba_24_Time=zeros(1,length(Caba_24))+(T_start+18*7-1);

%load('Hybrid_caba_treat_24_0422.mat', 'Cell_Daily', 'T')
load("./data/"+batchname+"_"+schedule+"_caba.mat")

%CELL VOLUME OVER TIME
for k=1:T
    Total_Cell_Coverage_Caba(k)= sum(sum(((Cell_Daily(:,:,k)==1) + (Cell_Daily(:,:,k)==2)+(Cell_Daily(:,:,k)==3) + (Cell_Daily(:,:,k)==4))/10))/10;
end

%Cab_Means={mean(Caba_6),mean(Caba_24)};

treatment_name = 'enza_with_caba';
    %ADB collect all the data on R and S
    original_batch_name = batch_name
    for i_batch=1:length(batch_names)
        single_batch_name = batch_names(i_batch)
        full_name = folder+single_batch_name+"_"+schedule+"_"+treatment_name;
        load(full_name+".mat");
            for k=1:T
        S_Daily(:,:,k)= (Cell_Daily(:,:,k)==1) + (Cell_Daily(:,:,k)==2);
        R_Daily(:,:,k)= (Cell_Daily(:,:,k)==3) + (Cell_Daily(:,:,k)==4);
            end
    length(batch_names)
    for k=1:T
        Total_S_Daily(k)=sum(sum(S_Daily(:,:,k)));
        Total_R_Daily(k)=sum(sum(R_Daily(:,:,k)));
    end
        all_S_Daily(i_batch,:)=Total_S_Daily;
        all_R_Daily(i_batch,:)=Total_R_Daily;
    end
    all_Daily_both = (all_S_Daily+all_R_Daily)/100;

Both_6={34.7615,18.0444,15.7451,32.6845,14.3223,10.0442,19.764,17.028,0,21.737,8.88654,19.6539,25.5023,19.135,23.3422,10.8765,16.1096,6.39154};
Both_6=cell2mat(Both_6);
Both_6_Time=zeros(1,length(Both_6))+(T_start);
Both_24={5.2,2.3,1.5,2.2,1.7,2.5};
Both_24=cell2mat(Both_24);
Both_24_Time=zeros(1,length(Both_24))+(T_start+18*7-1);

%load('Hybrid_Enza_Caba_treat_30_63_0322.mat', 'Cell_Daily', 'T')
%Hybrid_Enza_Caba_treat_24_1221_final
load("./data/"+batchname+"_"+schedule+"_enza_with_caba.mat")

%CELL VOLUME OVER TIME
for k=1:T
    Total_Cell_Coverage_Both(k)= sum(sum(((Cell_Daily(:,:,k)==1) + (Cell_Daily(:,:,k)==2)+(Cell_Daily(:,:,k)==3) + (Cell_Daily(:,:,k)==4))/10))/10;
end

% load('Hybrid_Patel_Enza_Caba_24_alternate.mat', 'Cell_Daily', 'T')
%
%     %CELL VOLUME OVER TIME
%     for k=1:T
%     Total_Cell_Coverage_Alternate(k)= sum(sum(((Cell_Daily(:,:,k)==1) + (Cell_Daily(:,:,k)==2)+(Cell_Daily(:,:,k)==3) + (Cell_Daily(:,:,k)==4))/10));
%     end

%Both_Means={mean(Both_6),mean(Both_24)};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Common parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

my_font_size = 18;
my_linewidth = 2;
my_scatterwidth = 1;
my_facecolor = [0.75,0.75,1];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Cell coverage with no treatment
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

h1 = figure(1);
ax = gca; ax.FontSize = my_font_size; ax.Box = 'on'; ax.TickDir = 'out'; ax.XMinorTick = 'off'; ax.YMinorTick = 'off'; ax.XLabel.Interpreter = 'latex'; ax.YLabel.Interpreter = 'latex';
%title('Cell Coverage with no Drug Therapy')
%set(gca,'FontSize',14)
xlabel('Time (Days)') %,'FontSize',20)
ylabel('Cell Coverage (\%)') %,'FontSize',20)
hold on

%ADB plotting min and max over runs - can amend to plot +/- 1sd from mean
    fill([TIME, fliplr(TIME)], [min(all_Daily_none), fliplr(max(all_Daily_none))],'k','EdgeColor','k','FaceAlpha', 0.2, 'EdgeAlpha', 0.3, 'linewidth', my_linewidth/2)
    %ADB plotting 95% confidence intervals
    N = size(all_Daily_none, 1)
    ci_95 = tinv([0.025, 0.975], N-1)'; %tinv
    s_mean = mean(all_Daily_none);
    s_sem = std(all_Daily_none)/sqrt(N);
    sci_95 = bsxfun(@times, s_sem, ci_95);
    fill([TIME, fliplr(TIME)], [s_mean+sci_95(1,:), fliplr(s_mean+sci_95(2,:))],'k','EdgeColor','k', 'FaceAlpha', 0.5, "EdgeAlpha", 1.0, 'linewidth', my_linewidth/2)

scatter(No_Therapy_6_Time, No_Therapy_6,'b', 'LineWidth', my_scatterwidth,'MarkerFaceColor',my_facecolor)
scatter(No_Therapy_12_Time, No_Therapy_12,'b', 'LineWidth', my_scatterwidth,'MarkerFaceColor',my_facecolor)
scatter(No_Therapy_16_Time, No_Therapy_16,'b', 'LineWidth', my_scatterwidth,'MarkerFaceColor',my_facecolor)
scatter(No_Therapy_19_Time, No_Therapy_19,'b', 'LineWidth', my_scatterwidth,'MarkerFaceColor',my_facecolor)
scatter(No_Therapy_24_Time, No_Therapy_24,'b', 'LineWidth', my_scatterwidth,'MarkerFaceColor',my_facecolor)
%plot(TIME, Total_Cell_Coverage_No_Therapy, 'LineWidth', my_linewidth, 'Color', 'k')

xlim([0,168-7])
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Cell coverage with enza
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

h2 = figure(2);
ax = gca; ax.FontSize = my_font_size; ax.Box = 'on'; ax.TickDir = 'out'; ax.XMinorTick = 'off'; ax.YMinorTick = 'off'; ax.XLabel.Interpreter = 'latex'; ax.YLabel.Interpreter = 'latex';
%title('Cell Coverage with Enzalutamide')
%set(gca,'FontSize',14)
hold on
xlabel('Time (Days)'); %,'FontSize',20)
ylabel('Cell Coverage (\%)'); %,'FontSize',20)
hold on

%ADB plotting min and max over runs - can amend to plot +/- 1sd from mean
    fill([TIME, fliplr(TIME)], [min(all_Daily_enza), fliplr(max(all_Daily_enza))],'k','EdgeColor','k','FaceAlpha', 0.2, 'EdgeAlpha', 0.3, 'linewidth', my_linewidth/2)
    %ADB plotting 95% confidence intervals
    N = size(all_Daily_enza, 1)
    ci_95 = tinv([0.025, 0.975], N-1)'; %tinv
    s_mean = mean(all_Daily_enza);
    s_sem = std(all_Daily_enza)/sqrt(N);
    sci_95 = bsxfun(@times, s_sem, ci_95);
    fill([TIME, fliplr(TIME)], [s_mean+sci_95(1,:), fliplr(s_mean+sci_95(2,:))],'k','EdgeColor','k', 'FaceAlpha', 0.5, "EdgeAlpha", 1.0, 'linewidth', my_linewidth/2)

scatter(Enza_6_Time, Enza_6,'b', 'LineWidth', my_scatterwidth,'MarkerFaceColor',my_facecolor)
scatter(Enza_12_Time, Enza_12,'b', 'LineWidth', my_scatterwidth,'MarkerFaceColor',my_facecolor)
scatter(Enza_16_Time, Enza_16,'b', 'LineWidth', my_scatterwidth,'MarkerFaceColor',my_facecolor)
scatter(Enza_22_Time, Enza_22,'b', 'LineWidth', my_scatterwidth,'MarkerFaceColor',my_facecolor)
scatter(Enza_24_Time, Enza_24,'b', 'LineWidth', my_scatterwidth,'MarkerFaceColor',my_facecolor)
%plot(TIMEt, Total_Cell_Coverage_Enza, 'LineWidth', my_linewidth, 'Color', 'k')

xlim([0,168-7])
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Cell coverage with caba
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

h3 = figure(3);
ax = gca; ax.FontSize = my_font_size; ax.Box = 'on'; ax.TickDir = 'out'; ax.XMinorTick = 'off'; ax.YMinorTick = 'off'; ax.XLabel.Interpreter = 'latex'; ax.YLabel.Interpreter = 'latex';
%set(gca,'FontSize',14)
hold on
%title('Cell Coverage with Cabazitaxel')
xlabel('Time (Days)'); %,'FontSize',20)
ylabel('Cell Coverage (\%)'); %,'FontSize',20)
hold on

%ADB plotting min and max over runs - can amend to plot +/- 1sd from mean
    fill([TIME, fliplr(TIME)], [min(all_Daily_caba), fliplr(max(all_Daily_caba))],'k','EdgeColor','k','FaceAlpha', 0.2, 'EdgeAlpha', 0.3, 'linewidth', my_linewidth/2)
    %ADB plotting 95% confidence intervals
    N = size(all_Daily_caba, 1)
    ci_95 = tinv([0.025, 0.975], N-1)'; %tinv
    s_mean = mean(all_Daily_caba);
    s_sem = std(all_Daily_caba)/sqrt(N);
    sci_95 = bsxfun(@times, s_sem, ci_95);
    fill([TIME, fliplr(TIME)], [s_mean+sci_95(1,:), fliplr(s_mean+sci_95(2,:))],'k','EdgeColor','k', 'FaceAlpha', 0.5, "EdgeAlpha", 1.0, 'linewidth', my_linewidth/2)

scatter(Caba_6_Time, Caba_6,'b', 'LineWidth', my_scatterwidth,'MarkerFaceColor',my_facecolor)
scatter(Caba_24_Time, Caba_24,'b', 'LineWidth', my_scatterwidth,'MarkerFaceColor',my_facecolor)
%plot(TIMEt, Total_Cell_Coverage_Caba, 'LineWidth', my_linewidth, 'Color', 'k')


xlim([0,168-7])
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Cell coverage with enza and caba
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

TIMEt1=1:T; %1:210; %?????????????????????????

h4 = figure(4);
ax = gca; ax.FontSize = my_font_size; ax.Box = 'on'; ax.TickDir = 'out'; ax.XMinorTick = 'off'; ax.YMinorTick = 'off'; ax.XLabel.Interpreter = 'latex'; ax.YLabel.Interpreter = 'latex';
%set(gca,'FontSize',14)
hold on
%title('Cell Coverage with both drugs')
xlabel('Time (Days)'); %,'FontSize',20)
ylabel('Cell Coverage (\%)'); %,'FontSize',20)
hold on

%ADB plotting min and max over runs - can amend to plot +/- 1sd from mean
    fill([TIME, fliplr(TIME)], [min(all_Daily_both), fliplr(max(all_Daily_both))],'k','EdgeColor','k','FaceAlpha', 0.2, 'EdgeAlpha', 0.3, 'linewidth', my_linewidth/2)
    %ADB plotting 95% confidence intervals
    N = size(all_Daily_both, 1)
    ci_95 = tinv([0.025, 0.975], N-1)'; %tinv
    s_mean = mean(all_Daily_both);
    s_sem = std(all_Daily_both)/sqrt(N);
    sci_95 = bsxfun(@times, s_sem, ci_95);
    fill([TIME, fliplr(TIME)], [s_mean+sci_95(1,:), fliplr(s_mean+sci_95(2,:))],'k','EdgeColor','k', 'FaceAlpha', 0.5, "EdgeAlpha", 1.0, 'linewidth', my_linewidth/2)

scatter(Both_6_Time, Both_6,'b', 'LineWidth', my_scatterwidth,'MarkerFaceColor',my_facecolor)
scatter(Both_24_Time, Both_24,'b', 'LineWidth', my_scatterwidth,'MarkerFaceColor',my_facecolor)
%plot(TIMEt1, Total_Cell_Coverage_Both, 'LineWidth', my_linewidth, 'Color', 'k')
axis([0 168-7 0 35])
%plot(TIMEt, Total_Cell_Coverage_Alternate, 'LineWidth', 2, 'Color', 'g')



hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Adjust figure positions and save figures
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

set(h1,'Units','Inches');
set(h2,'Units','Inches');
set(h3,'Units','Inches');
set(h4,'Units','Inches');

pos1 = get(h1,'Position');
pos2 = get(h2,'Position');
pos3 = get(h3,'Position');
pos4 = get(h4,'Position');

set(h1,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos1(3), pos1(4)])
set(h2,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos2(3), pos2(4)])
set(h3,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos3(3), pos3(4)])
set(h4,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos4(3), pos4(4)])

print(h1,"./data/"+batchname+"_"+schedule+"_none_comp.pdf",'-dpdf','-r0')
print(h2,"./data/"+batchname+"_"+schedule+"_enza_comp.pdf",'-dpdf','-r0')
print(h3,"./data/"+batchname+"_"+schedule+"_caba_comp.pdf",'-dpdf','-r0')
print(h4,"./data/"+batchname+"_"+schedule+"_enza_with_caba_comp.pdf",'-dpdf','-r0')

end
