% Script to load data from transition data base and plot, including plots
% with correlations
% Davis and Atkins 2024

clear
close all
addpath('Correlations')

%% Load data

tab_data = readtable('Full_data_set.csv');
data = table2array(tab_data(:,1:3));

PG = data(:,1);
Re_theta = data(:,2);
Tu = data(:,3);

%% Sort

[~, Tu_ind] = sort(Tu);
[~, PG_ind] = sort(PG);
sorted_Tu = data(Tu_ind,:);
sorted_PG = data(PG_ind,:);

ZPG = sorted_Tu(sorted_Tu(:,1)==0,:);
PG = sorted_Tu(sorted_Tu(:,1)~=0,:);
tab_ZPG = tab_data(table2array(tab_data(:,1))==0,:);
tab_PG = tab_data(table2array(tab_data(:,1))~=0,:);

%% Evaluate correlations

Np = 1000;
Tu = linspace(0.01,10,Np);
lambda = linspace(-0.1,0.15,Np);

% Zero pressure gradient lines
AGS_ZPG = Abu_Ghannam_Shaw(zeros(1,Np),Tu);
Dunham_ZPG = Dunham(zeros(1,Np),Tu);
FJ_ZPG = Fasihfar_Johnson(zeros(1,Np),Tu);
HG_ZPG = Hall_Gibbings(zeros(1,Np),Tu);
Hourmouziadis_ZPG = Hourmouziadis(zeros(1,Np),Tu);
LM_ZPG = Langtry_Menter(zeros(1,Np),Tu);
Mayle_ZPG = Mayle(zeros(1,Np),Tu);
DA_ZPG = Davis_Atkins(zeros(1,Np),Tu);
Seyb_ZPG = Seyb(zeros(1,Np),Tu);
Suzen_ZPG = Suzen_et_al(zeros(1,Np),Tu);
VDB_ZPG = Van_Driest_Blumer(zeros(1,Np),Tu);

% Variable pressure gradient lines at set Tu
Tu_PG = [0.025 0.25 0.5 1 1.5 2.5 6];

AGS_PG = zeros(Np,length(Tu_PG));
Dunham_PG = zeros(Np,length(Tu_PG));
FJ_PG = zeros(Np,length(Tu_PG));
LM_PG = zeros(Np,length(Tu_PG));
Seyb_PG = zeros(Np,length(Tu_PG));
Suzen_PG = zeros(Np,length(Tu_PG));

for i = 1:length(Tu_PG)
    T = ones(1,Np)*Tu_PG(i);
    AGS_PG(:,i) = Abu_Ghannam_Shaw(lambda,T);
    Dunham_PG(:,i) = Dunham(lambda,T);
    FJ_PG(:,i) = Fasihfar_Johnson(lambda,T);
    LM_PG(:,i) = Langtry_Menter(lambda,T);
    Seyb_PG(:,i) = Seyb(lambda,T);
    Suzen_PG(:,i) = Suzen_et_al(lambda,T);
end

Dunham_PG(Dunham_PG<0) = nan;

%% Figures in paper

% ZPG data with sources
markers = {'o'; 'o'; 'square'; 'square'; '^'; '^'; 'diamond'; 'diamond'; '*';'x'; '>'; '>'; '_';'pentagram'; 'hexagram';'+'};
fill = [1 2 1 2 1 2 1 2 0 0 1 2 0 1 1 0];

[names, ~, IC] = unique(table2array(tab_ZPG(:,4)));
occ = cumsum(histcounts(IC)');

figure(1);
h = zeros(size(occ));
scatter(table2array(tab_ZPG(1:occ(1),3)),table2array(tab_ZPG(1:occ(1),2)),30,'k',cell2mat(markers(1)),'filled')
hold on
h(1) = scatter(nan,30,'k',cell2mat(markers(1)),'filled','DisplayName',cell2mat(names(1)));
for i = 2:length(occ)
    if fill(i)
        if fill(i) == 1
            scatter(table2array(tab_ZPG(occ(i-1)+1:occ(i),3)),table2array(tab_ZPG(occ(i-1)+1:occ(i),2)),30,'k',cell2mat(markers(i)),'filled')
            h(i) = scatter(nan,30,'k',cell2mat(markers(i)),'filled','DisplayName',cell2mat(names(i)));
        else
            scatter(table2array(tab_ZPG(occ(i-1)+1:occ(i),3)),table2array(tab_ZPG(occ(i-1)+1:occ(i),2)),30,'w',cell2mat(markers(i)),'filled','MarkerEdgeColor','k','MarkerFaceColor','w')
            h(i) = scatter(nan,30,'k',cell2mat(markers(i)),'filled','DisplayName',cell2mat(names(i)),'MarkerEdgeColor','k','MarkerFaceColor','w');
        end
    else
        scatter(table2array(tab_ZPG(occ(i-1)+1:occ(i),3)),table2array(tab_ZPG(occ(i-1)+1:occ(i),2)),30,'k',cell2mat(markers(i)))
        h(i) = scatter(nan,30,'k',cell2mat(markers(i)),'DisplayName',cell2mat(names(i)));
    end
end
hold off
grid on
ylim([0 1500])
xticks([0 2 4 6 8 10])
title('Zero pressure gradient transition data labelled by source')
xlabel('{\it Tu} (%)')
ylabel('{\it Re}_{\theta,t}')
legend(h,'Location','northoutside','NumColumns',3)

% PG data with sources
markers = {'o'; 'square'; 'o'; '^'; 'diamond'; 'diamond'; '*'; 'square'; 'x'; '>'; '+'};
fill = [1 1 2 2 1 2 0 2 0 1 0];

[names, ~, IC] = unique(table2array(tab_PG(:,4)));
occ = cumsum(histcounts(IC)');

figure(2)
clear h
h = zeros(size(occ));
scatter(table2array(tab_PG(1:occ(1),1)),table2array(tab_PG(1:occ(1),2)),30,'k',cell2mat(markers(1)),'filled')
hold on
h(1) = scatter(nan,30,'k',cell2mat(markers(1)),'filled','DisplayName',cell2mat(names(1)));
for i = 2:length(occ)
    if fill(i)
        if fill(i)==1
            scatter(table2array(tab_PG(occ(i-1)+1:occ(i),1)),table2array(tab_PG(occ(i-1)+1:occ(i),2)),30,'k',cell2mat(markers(i)),'filled')
            h(i) = scatter(nan,30,'k',cell2mat(markers(i)),'filled','DisplayName',cell2mat(names(i)));
        else
            scatter(table2array(tab_PG(occ(i-1)+1:occ(i),1)),table2array(tab_PG(occ(i-1)+1:occ(i),2)),30,'k',cell2mat(markers(i)),'filled','MarkerEdgeColor','k','MarkerFaceColor','w')
            h(i) = scatter(nan,30,'k',cell2mat(markers(i)),'filled','DisplayName',cell2mat(names(i)),'MarkerEdgeColor','k','MarkerFaceColor','w');
        end
    else
        scatter(table2array(tab_PG(occ(i-1)+1:occ(i),1)),table2array(tab_PG(occ(i-1)+1:occ(i),2)),30,'k',cell2mat(markers(i)))
        h(i) = scatter(nan,30,'k',cell2mat(markers(i)),'DisplayName',cell2mat(names(i)));
    end
end
hold off
grid on
ylim([0 1500])
xlim([-0.1 0.15])
title('Pressure gradient transition data labelled by source. Symbols are matched to Fig. 2 where appropriate.')
xlabel('{\it\lambda}_\theta')
ylabel('{\it Re}_{\theta,t}')
legend(h,'Location','northoutside','NumColumns',3)

% ZPG data with correlations
figure(3)
scatter(ZPG(:,3),ZPG(:,2),25,'k','.')
hold on
plot(Tu,Mayle_ZPG,'k-.')
plot(Tu,DA_ZPG,'k-')
hold off
grid on
ylim([0 1500])
xticks([0 2 4 6 8 10])
title('Zero pressure gradient transition data with correlations')
h = {'Data' 'Mayle 1991' 'New best fit'};
xlabel('{\it Tu} (%)')
ylabel('{\it Re}_{\theta,t}')
legend(h,'Location','northoutside','NumColumns',3)

% PG data coloured by Tu
Nc = 10000;
cmap = turbo(Nc);
cm_Tu = round(interp1([0 6],[1 Nc],sorted_PG(:,3)));
cm_Tu(isnan(cm_Tu)) = Nc;
cm_Tu = cmap(cm_Tu',:);

% PD data with correlations
Tu_PG = [0.025 0.25 0.5 1 1.5 2.5 6];

AGS_PG = zeros(Np,length(Tu_PG));
LM_PG = zeros(Np,length(Tu_PG));

for i = 1:length(Tu_PG)
    T = ones(1,Np)*Tu_PG(i);
    AGS_PG(:,i) = Abu_Ghannam_Shaw(lambda,T);
    LM_PG(:,i) = Langtry_Menter(lambda,T);
end

line_cs = round(interp1([0 6],[1 Nc],Tu_PG));
line_cs = num2cell(cmap(line_cs,:),2);

figure(4)
scatter(sorted_PG(:,1),sorted_PG(:,2),10,cm_Tu,'filled')
hold on
h = plot(lambda,AGS_PG,'k-.','LineWidth',1);set(h, {'color'}, line_cs);
h = plot(lambda,LM_PG,'k-','LineWidth',1);set(h, {'color'}, line_cs);
clear h
h(2) = plot(nan,'k-.','DisplayName','Abu-Ghannam and Shaw 1980');
h(1) = plot(nan,'k-','DisplayName','Langtry and Menter 2009');
hold off
grid on
ylim([0 1500])
xlim([-0.1 0.15])
title('Pressure gradient transition data coloured by Tu. Correlations from Abu-Ghannam and Shaw [6] (dashed) and Langtry and Menter [43] (solid) are plotted with lines coloured to the same scale.')
xlabel('{\it\lambda}_\theta')
ylabel('{\it Re}_{\theta,t}')

% Histograms
Tu_bins = linspace(0.25,9.75,20);
Re_bins = linspace(50,1450,15);
pg_bins = linspace(-0.1,0.15,26);

x_tu = data(:,3);
y_pg = data(:,1);
z_re = data(:,2);

X = [y_pg,x_tu];

figure(5)
hist3(X,'ctrs',{pg_bins Tu_bins},'CdataMode','auto')
view(2)
c = colorbar;
clim([0 10])
a = gca;
cm = a.Colormap;
a.Colormap = cm(1:24:end,:);
xlim([-0.1 0.15])
ylim([0 10])
yticks([0 2 4 6 8 10])
xticks([-0.1 -0.05 0 0.05 0.1 0.15])
c.Label.String = 'Number of data points';
c.Ticks = [0.5 2.5 4.5 6.5 8.5 10.5]*10/11;
c.TickLabels = {'0','2','4','6','8','\geq 10'};
c.Location = 'northoutside';c.Label.Color = 'k'; c.Color = 'k';
a = gca;
a.SortMethod = 'childorder';
title('Histogram showing of data points for Tu and \lambda_\theta')
xlabel('{\it\lambda}_\theta')
ylabel('{\it Tu} (%)')

X = [y_pg,z_re];

figure(6)
hist3(X,'ctrs',{pg_bins Re_bins},'CdataMode','auto')
view(2)
c = colorbar;
clim([0 10])
a = gca;
cm = a.Colormap;
a.Colormap = cm(1:24:end,:);
xlim([-0.1 0.15])
ylim([0 1500])
yticks([0 500 1000 1500])
xticks([-0.1 -0.05 0 0.05 0.1 0.15])
c.Label.String = 'Number of data points';
c.Ticks = [0.5 2.5 4.5 6.5 8.5 10.5]*10/11;
c.TickLabels = {'0','2','4','6','8','\geq 10'};
c.Location = 'northoutside';c.Label.Color = 'k'; c.Color = 'k';
a.SortMethod = 'childorder';
title('Histogram showing of data points for Re_{\theta,t} and \lambda_\theta')
xlabel('{\it\lambda}_\theta')
ylabel('{\it Re}_{\theta,t}')

% Error at each point

tu = data(:,3);
pg = data(:,1);
re = data(:,2);
re_lm = Langtry_Menter(pg,tu);

e_re = 100*abs(re-re_lm)./re_lm;

Nc = 10000;
cmap = turbo(Nc);
cm_e = round(interp1([0 50],[1 Nc],e_re));
cm_e(isnan(cm_e)) = Nc;
cm_e = cmap(cm_e',:);

figure(7)
scatter(pg,re,10,cm_e,'filled')
grid on
ylim([0 1500])
xlim([-0.1 0.15])
c = colorbar; clim([0 50]); colormap turbo; c.Label.String = 'Relative error in Reθ,t calculated by the correlation of Langtry and Menter [43] for each data point.';
c.Ticks = [0 10 20 30 40 50];
c.TickLabels = {'0','10','20','30','40','\geq 50'};
c.Location = 'northoutside';c.Label.Color = 'k'; c.Color = 'k';
xlabel('{\it\lambda}_\theta')
ylabel('{\it Re}_{\theta,t}')

%% Additional plots for interest

% ZPG data with AGS and LM correlations

figure(9)
scatter(ZPG(:,3),ZPG(:,2),25,'k','filled')
hold on
plot(Tu,AGS_ZPG,'k-')
plot(Tu,LM_ZPG,'k--')
hold off
grid on
ylim([0 1500])
xlabel('{\it Tu} (%)')
ylabel('{\it Re}_{\theta,t}')
h = {'Data' 'Abu-Ghannam and Shaw 1980' 'Langtry and Menter 2009'};
legend(h,'Location','northoutside','NumColumns',2)

% ZPG data with other ZPG correlations

figure(10)
scatter(ZPG(:,3),ZPG(:,2),25,'k','filled')
hold on
plot(Tu,Hourmouziadis_ZPG,'k-')
plot(Tu,HG_ZPG,'k--')
plot(Tu,VDB_ZPG,'k:')
hold off
grid on
ylim([0 1500])
xlabel('{\it Tu} (%)')
ylabel('{\it Re}_{\theta,t}')
h = {'Data' 'Hourmouziadis 1989' 'Hall and Gibbings 1972' 'Van Driest and Blumer 1963'};
legend(h,'Location','northoutside','NumColumns',2)

% ZPG data with other PG correlations

figure(11)
scatter(ZPG(:,3),ZPG(:,2),25,'k','filled')
hold on
plot(Tu,FJ_ZPG,'k-')
plot(Tu,Dunham_ZPG,'k--')
plot(Tu(1:10:end),Seyb_ZPG(1:10:end),'k:','MarkerSize',2)
plot(Tu(1:10:end),Suzen_ZPG(1:10:end),'k-x','MarkerSize',4)
hold off
grid on
ylim([0 1500])
xlabel('{\it Tu} (%)')
ylabel('{\it Re}_{\theta,t}')
h = {'Data' 'Fasihfar and Johnson 1992' 'Dunham 1972' 'Seyb 1972' 'Suzen et al 2002'};
legend(h,'Location','northoutside','NumColumns',2)

% PG data with other PG correlations
figure(12)
scatter(sorted_PG(:,1),sorted_PG(:,2),30,cm_Tu,'filled')
hold on
h = plot(lambda,FJ_PG,'k-');set(h, {'color'}, line_cs);
h = plot(lambda,Dunham_PG,'k--');set(h, {'color'}, line_cs);
h = plot(lambda,Seyb_PG,'k:');set(h, {'color'}, line_cs);
h = plot(lambda,Suzen_PG,'k-.');set(h, {'color'}, line_cs);

clear h
h(1) = plot(nan,'k','DisplayName','Fasihfar and Johnson 1992');
h(2) = plot(nan,'k--','DisplayName','Dunham 1972');
h(3) = plot(nan,'k:','DisplayName','Seyb 1972');
h(4) = plot(nan,'k-.','DisplayName','Suzen et al 2002');

hold off
grid on
ylim([0 1500])
xlim([-0.1 0.15])
c = colorbar; clim([0 5]); colormap turbo; c.Ticks = [0 1 2 3 4 5]; c.Label.String = '{\it Tu} (%)';
legend(h,'Location','northoutside','NumColumns',2)
xlabel('{\it\lambda}_\theta')
ylabel('{\it Re}_{\theta,t}')