% Analyze Scaling Data
% Hannah Ross
% 4/20/2020, updated 9/26/2020

% Analyzes turbine performance and free surface deformation data collected
% in the University of Washington flume in September 2020.

clear

%% Load data
load('geometryData.mat')
load('turbineData.mat')
load('flowData.mat')
fieldNames = fields(tds);

%% Constants
D = 0.172;  % turbine diameter (m)
g = 9.81;  % acceleration of gravity (m/s^2)
c = 0.0405;  % blade chord length (m)

%% Calculate blockage ratio, Reynolds number, Froude numbers, and aspect ratio
for i = 1:length(fieldNames)
    blockage.(fieldNames{i}) = (turb.(fieldNames{i}).D*turb.(fieldNames{i}).H + channel.(fieldNames{i}).SS)/(channel.(fieldNames{i}).w*channel.(fieldNames{i}).d0);
    [~,nu,~] = getWaterProps(channel.(fieldNames{i}).T,channel.(fieldNames{i}).dt + turb.(fieldNames{i}).H/2,0);
    Reynolds.(fieldNames{i}) = mean(fds.(fieldNames{i}).Uinf)*c/nu;
    Froude_d0.(fieldNames{i}) = mean(fds.(fieldNames{i}).Uinf)/sqrt(g*channel.(fieldNames{i}).d0);
    Froude_ds.(fieldNames{i}) = mean(fds.(fieldNames{i}).Uinf)/sqrt(g*channel.(fieldNames{i}).dt);
    AR.(fieldNames{i}) = turb.(fieldNames{i}).H/turb.(fieldNames{i}).D;
    clear nu
end

%% Calculate TSR and CP
for i = 1:length(fieldNames)
    for j = 1:length(tds.(fieldNames{i}).cycleAvTSR)
        TSR.(fieldNames{i})(j) = median(tds.(fieldNames{i}).cycleAvTSR{j});
        TSRUCert.(fieldNames{i})(1,j) = quantile(tds.(fieldNames{i}).cycleAvTSR{j},0.25);
        TSRUCert.(fieldNames{i})(2,j) = quantile(tds.(fieldNames{i}).cycleAvTSR{j},0.75);
        CP.(fieldNames{i})(j) = median(tds.(fieldNames{i}).cycleAveff{j});
        CPUCert.(fieldNames{i})(1,j) = quantile(tds.(fieldNames{i}).cycleAveff{j},0.25);
        CPUCert.(fieldNames{i})(2,j) = quantile(tds.(fieldNames{i}).cycleAveff{j},0.75);
    end
end

%% Plot diagram of test rig
f1 = figure(1);
clf
im1 = imread('C:\Users\hanna\Dropbox (MREL)\Papers\ReynoldsNumber\Figures\ARTurbine.jpg');
im1R = imresize(im1,1);
imshow(im1R)
hold on
scatter(473,108,30,'k','Filled')
plot([473 590],[108 123],'k','LineWidth',1.5)
text(610,123,'Six-axis load cell (ATI Mini45)','FontName','Times New Roman','FontSize',15)
scatter(473,210,30,'k','Filled')
plot([473 590],[210 225],'k','LineWidth',1.5)
text(610,254,{'Servomotor and rotary encoder ','(Yaskawa SGMCS-05B3C41)'},'FontName','Times New Roman','FontSize',15)
scatter(473,858,30,'k','Filled')
plot([473 590],[858 843],'k','LineWidth',1.5)
text(610,872,{'Submersible six-axis load cell','(ATI Mini40)'},'FontName','Times New Roman','FontSize',15)
f1.Position = [320 118.3333 670 492.6667];

%% Plot top view diagram of experimental setup
f2 = figure(2);
clf
R = D/2;
ang = 0:0.1:2*pi+0.1;

% Plot measurement locations
scatter([-7 -2 -1 1 2],[0 0 0 0 0],100,'r','x','LineWidth',2);
hold on
scatter(-5,0,1000,'b','.');

% Plot end plates
xp = R*sin(ang)/D;
yp = R*cos(ang)/D;
plot(xp,yp,'k--','LineWidth',1.5)

xlim([-7.5 2.5])
ylim([-4.42 4.42])
xlabel('$X/D$','Interpreter','latex')
ylabel('$Y/D$','Interpreter','latex')
text(-8.05,-3.95,'Flow','FontSize',16,'FontName','Times New Roman','BackgroundColor','w')
annotation('textarrow',[0.24 0.35],[0.17 0.17],'LineWidth',1,'HeadStyle','vback3');
text(-3.3,-2.1,'Turbine Swept Area','FontSize',16,'FontName','Times New Roman','BackgroundColor','w')
annotation('textarrow',[0.615 0.675],[0.37 0.47],'LineWidth',1,'HeadStyle','vback3');
set(gca,'FontName','Times New Roman','FontSize',16)
[legh,icons,plots,legend_text] = legend('Transducer Locations','ADV Location','location','northwest');
icons(3).Children.LineWidth = 2;
icons(3).Children.MarkerSize = 10;
icons(4).Children.MarkerSize = 30;
grid on
box on
axis equal

%% Plot CP
cmpD = get(groot,'DefaultAxesColorOrder');
k = [0 0 0];
cmpW = winter(7);
cmpP = parula(7);
cmp = [k; cmpD(2,:); cmpP(3,:); cmpP(5,:); cmpD(3,:)];

f3 = figure(3);
clf
s1 = subplot(2,2,1);
fill([TSR.Baseline fliplr(TSR.Baseline)],[CPUCert.Baseline(1,:) fliplr(CPUCert.Baseline(2,:))],cmp(1,:),'EdgeColor','none','FaceAlpha',0.2)
hold on
fill([TSR.Reynolds fliplr(TSR.Reynolds)],[CPUCert.Reynolds(1,:) fliplr(CPUCert.Reynolds(2,:))],cmp(3,:),'EdgeColor','none','FaceAlpha',0.3)
p1 = plot(TSR.Baseline,CP.Baseline,'-o','LineWidth',2,'Color',cmp(1,:),'MarkerEdgeColor','none','MarkerFaceColor',cmp(1,:),'MarkerSize',8);
p3 = plot(TSR.Reynolds,CP.Reynolds,'-d','LineWidth',2,'Color',cmp(3,:),'MarkerEdgeColor','none','MarkerFaceColor',cmp(3,:),'MarkerSize',8);
xlim([0.55 3.2])
ylim([0 0.32])
grid on
legend([p3 p1],'Reynolds','Baseline','location','northwest','FontName','Times New Roman')
xlabel('$\lambda$','Interpreter','latex')
ylabel('$C_P$','Interpreter','latex')
set(gca,'FontName','Times New Roman','FontSize',16)
text(0.6,0.3415,'(a)','FontName','Times New Roman','FontSize',16)

s2 = subplot(2,2,2);
fill([TSR.Baseline fliplr(TSR.Baseline)],[CPUCert.Baseline(1,:) fliplr(CPUCert.Baseline(2,:))],cmp(1,:),'EdgeColor','none','FaceAlpha',0.2)
hold on
fill([TSR.Blockage fliplr(TSR.Blockage)],[CPUCert.Blockage(1,:) fliplr(CPUCert.Blockage(2,:))],cmp(2,:),'EdgeColor','none','FaceAlpha',0.3)
p1 = plot(TSR.Baseline,CP.Baseline,'-o','LineWidth',2,'Color',cmp(1,:),'MarkerEdgeColor','none','MarkerFaceColor',cmp(1,:),'MarkerSize',8);
p2 = plot(TSR.Blockage,CP.Blockage,'-s','LineWidth',2,'Color',cmp(2,:),'MarkerEdgeColor','none','MarkerFaceColor',cmp(2,:),'MarkerSize',8);
xlim([0.55 3.2])
ylim([0 0.32])
grid on
legend([p2 p1],'Blockage','Baseline','location','northwest','FontName','Times New Roman')
xlabel('$\lambda$','Interpreter','latex')
ylabel('$C_P$','Interpreter','latex')
set(gca,'FontName','Times New Roman','FontSize',16)
text(0.6,0.3415,'(b)','FontName','Times New Roman','FontSize',16)

s3 = subplot(2,2,3);
fill([TSR.Baseline fliplr(TSR.Baseline)],[CPUCert.Baseline(1,:) fliplr(CPUCert.Baseline(2,:))],cmp(1,:),'EdgeColor','none','FaceAlpha',0.2)
hold on
fill([TSR.Froude_ds fliplr(TSR.Froude_ds)],[CPUCert.Froude_ds(1,:) fliplr(CPUCert.Froude_ds(2,:))],cmp(5,:),'EdgeColor','none','FaceAlpha',0.3)
p1 = plot(TSR.Baseline,CP.Baseline,'-o','LineWidth',2,'Color',cmp(1,:),'MarkerEdgeColor','none','MarkerFaceColor',cmp(1,:),'MarkerSize',8);
p4 = plot(TSR.Froude_ds,CP.Froude_ds,'-v','LineWidth',2,'Color',cmp(5,:),'MarkerEdgeColor','none','MarkerFaceColor',cmp(5,:),'MarkerSize',8);
xlim([0.55 3.2])
ylim([0 0.32])
grid on
legend([p4 p1],'Submergence','Baseline','location','northwest','FontName','Times New Roman')
xlabel('$\lambda$','Interpreter','latex')
ylabel('$C_P$','Interpreter','latex')
set(gca,'FontName','Times New Roman','FontSize',16)
text(0.6,0.3415,'(c)','FontName','Times New Roman','FontSize',16)

s4 = subplot(2,2,4);
fill([TSR.Baseline fliplr(TSR.Baseline)],[CPUCert.Baseline(1,:) fliplr(CPUCert.Baseline(2,:))],cmp(1,:),'EdgeColor','none','FaceAlpha',0.2)
hold on
fill([TSR.Froude_d0 fliplr(TSR.Froude_d0)],[CPUCert.Froude_d0(1,:) fliplr(CPUCert.Froude_d0(2,:))],cmp(4,:),'EdgeColor','none','FaceAlpha',0.3)
p1 = plot(TSR.Baseline,CP.Baseline,'-o','LineWidth',2,'Color',cmp(1,:),'MarkerEdgeColor','none','MarkerFaceColor',cmp(1,:),'MarkerSize',8);
p5 = plot(TSR.Froude_d0,CP.Froude_d0,'-^','LineWidth',2,'Color',cmp(4,:),'MarkerEdgeColor','none','MarkerFaceColor',cmp(4,:),'MarkerSize',8);
xlim([0.55 3.2])
ylim([0 0.32])
grid on
legend([p5 p1],'Froude','Baseline','location','northwest','FontName','Times New Roman')
xlabel('$\lambda$','Interpreter','latex')
ylabel('$C_P$','Interpreter','latex')
set(gca,'FontName','Times New Roman','FontSize',16)
text(0.6,0.3415,'(d)','FontName','Times New Roman','FontSize',16)

f3.Position = [317 778.3 1039.3 931.3];
s1.Position = [0.0800 0.5838 0.4047 0.3412];
s2.Position = [0.5803 0.5838 0.4047 0.3412]; 
s3.Position = [0.0800 0.1100 0.4047 0.3412]; 
s4.Position = [0.5803 0.1100 0.4047 0.3412]; 
  
%% Calculate cycle-averaged free surface height
for i = 1:length(fieldNames)-5
    for j = 1:length(fd.(fieldNames{i}))
        numSamples = floor(length(fd.(fieldNames{i})(j).FST.FrStr)/length(tds.(fieldNames{i}).cycleAvturbPos{j}));
        for k = 1:length(tds.(fieldNames{i}).cycleAvturbPos{j})
            for m = 1:size(fd.(fieldNames{i})(j).FST.TurbInterp,2)
                FST.(fieldNames{i}){j}(k,m) = mean(fd.(fieldNames{i})(j).FST.TurbInterp((k-1)*numSamples+1:k*numSamples,m));
            end
        end
    end
end

%% Calculate flow depth
for i = 1:length(fieldNames)-5
    for j = 1:length(FST.(fieldNames{i}))
        d.(fieldNames{i}){j} = median(FST.(fieldNames{i}){j});
        dUCert.(fieldNames{i}){j}(1,:) = quantile(FST.(fieldNames{i}){j},0.25);
        dUCert.(fieldNames{i}){j}(2,:) = quantile(FST.(fieldNames{i}){j},0.75);
    end
end

%% Plot free surface height at max CP point
for i = 1:length(fieldNames)
    [~,maxCPind.(fieldNames{i})] = max(CP.(fieldNames{i}));
end
FSTLocs = [-2 -1 1 2];
figure(4)
clf
fill([FSTLocs fliplr(FSTLocs)],[dUCert.Baseline{maxCPind.Baseline}(1,:)./fds.Baseline.FST.FrStrInterp(maxCPind.Baseline) fliplr(dUCert.Baseline{maxCPind.Baseline}(2,:)./fds.Baseline.FST.FrStrInterp(maxCPind.Baseline))],cmp(1,:),'EdgeColor','none','FaceAlpha',0.15)
hold on
fill([FSTLocs fliplr(FSTLocs)],[dUCert.Froude_d0{maxCPind.Froude_d0}(1,:)./fds.Froude_d0.FST.FrStrInterp(maxCPind.Froude_d0) fliplr(dUCert.Froude_d0{maxCPind.Froude_d0}(2,:)./fds.Froude_d0.FST.FrStrInterp(maxCPind.Froude_d0))],cmp(4,:),'EdgeColor','none','FaceAlpha',0.25)
fill([FSTLocs fliplr(FSTLocs)],[dUCert.Froude_ds{maxCPind.Froude_ds}(1,:)./fds.Froude_ds.FST.FrStrInterp(maxCPind.Froude_ds) fliplr(dUCert.Froude_ds{maxCPind.Froude_ds}(2,:)./fds.Froude_ds.FST.FrStrInterp(maxCPind.Froude_ds))],cmp(5,:),'EdgeColor','none','FaceAlpha',0.25)
fill([FSTLocs fliplr(FSTLocs)],[dUCert.Blockage{maxCPind.Blockage}(1,:)./fds.Blockage.FST.FrStrInterp(maxCPind.Blockage) fliplr(dUCert.Blockage{maxCPind.Blockage}(2,:)./fds.Blockage.FST.FrStrInterp(maxCPind.Blockage))],cmp(2,:),'EdgeColor','none','FaceAlpha',0.1)
fill([FSTLocs fliplr(FSTLocs)],[dUCert.Reynolds{maxCPind.Reynolds}(1,:)./fds.Reynolds.FST.FrStrInterp(maxCPind.Reynolds) fliplr(dUCert.Reynolds{maxCPind.Reynolds}(2,:)./fds.Reynolds.FST.FrStrInterp(maxCPind.Reynolds))],cmp(3,:),'EdgeColor','none','FaceAlpha',0.1)
d1 = plot(FSTLocs,d.Baseline{maxCPind.Baseline}./fds.Baseline.FST.FrStrInterp(maxCPind.Baseline),'-o','LineWidth',2,'Color',cmp(1,:),'MarkerEdgeColor','none','MarkerFaceColor',cmp(1,:),'MarkerSize',8);
d2 = plot(FSTLocs,d.Froude_d0{maxCPind.Froude_d0}./fds.Froude_d0.FST.FrStrInterp(maxCPind.Froude_d0),'-^','LineWidth',2,'Color',cmp(4,:),'MarkerEdgeColor','none','MarkerFaceColor',cmp(4,:),'MarkerSize',8);
d3 = plot(FSTLocs,d.Froude_ds{maxCPind.Froude_ds}./fds.Froude_ds.FST.FrStrInterp(maxCPind.Froude_ds),'-v','LineWidth',2,'Color',cmp(5,:),'MarkerEdgeColor','none','MarkerFaceColor',cmp(5,:),'MarkerSize',8);
d4 = plot(FSTLocs,d.Blockage{maxCPind.Blockage}./fds.Blockage.FST.FrStrInterp(maxCPind.Blockage),'-s','LineWidth',2,'Color',cmp(2,:),'MarkerEdgeColor','none','MarkerFaceColor',cmp(2,:),'MarkerSize',8);
d5 = plot(FSTLocs,d.Reynolds{maxCPind.Reynolds}./fds.Reynolds.FST.FrStrInterp(maxCPind.Reynolds),'-d','LineWidth',2,'Color',cmp(3,:),'MarkerEdgeColor','none','MarkerFaceColor',cmp(3,:),'MarkerSize',8);
plot([0 0],[0.9 1.1],'k--','LineWidth',1)
grid on
ylim([0.955 1.015])
legend([d5 d4 d3 d2 d1],'Reynolds','Blockage','Submergence','Froude','Baseline','FontName','Times New Roman','location','southwest')
xlabel('$X/D$','Interpreter','latex')
ylabel('$d/d_0$','Interpreter','latex')
set(gca,'FontName','Times New Roman','FontSize',16)

%% Calculate percent changes in CP at max CP point
for i = 1:length(fieldNames)
    if ~contains(fieldNames{i},char('Baseline'))
        CPPerc.(fieldNames{i}) = (CP.(fieldNames{i})(maxCPind.(fieldNames{i})) - CP.Baseline(maxCPind.Baseline))./CP.Baseline(maxCPind.Baseline)*100;
        blockagePerc.(fieldNames{i}) = (blockage.(fieldNames{i}) - blockage.Baseline)./blockage.Baseline*100;
        ReynoldsPerc.(fieldNames{i}) = (Reynolds.(fieldNames{i}) - Reynolds.Baseline)./Reynolds.Baseline*100;
        Frouded0Perc.(fieldNames{i}) = (Froude_d0.(fieldNames{i}) - Froude_d0.Baseline)./Froude_d0.Baseline*100;
        FroudedsPerc.(fieldNames{i}) = (Froude_ds.(fieldNames{i}) - Froude_ds.Baseline)./Froude_ds.Baseline*100;
    end
end

%% Compare Fr_s and Fr_0 CP
f5 = figure(5);
clf
s1 = subplot(2,2,1);
p5 = plot(TSR.Froude_ds_5,CP.Froude_ds_5,'-v','LineWidth',2,'Color',cmp(5,:),'MarkerEdgeColor','none','MarkerFaceColor',cmp(5,:),'MarkerSize',8);
hold on
p4 = plot(TSR.Froude_ds_4,CP.Froude_ds_4,'-p','LineWidth',2,'Color',cmpW(2,:),'MarkerEdgeColor','none','MarkerFaceColor',cmpW(2,:),'MarkerSize',10);
p3 = plot(TSR.Froude_ds_3,CP.Froude_ds_3,'-o','LineWidth',2,'Color',cmp(1,:),'MarkerEdgeColor','none','MarkerFaceColor',cmp(1,:),'MarkerSize',8);
p2 = plot(TSR.Froude_ds_2,CP.Froude_ds_2,'->','LineWidth',2,'Color',cmpW(4,:),'MarkerEdgeColor','none','MarkerFaceColor',cmpW(4,:),'MarkerSize',8);
p1 = plot(TSR.Froude_ds_1,CP.Froude_ds_1,'-<','LineWidth',2,'Color',cmpW(5,:),'MarkerEdgeColor','none','MarkerFaceColor',cmpW(5,:),'MarkerSize',8);
xlim([0.55 3.2])
ylim([0 0.32])
grid on
legend([p5 p4 p3 p2 p1],'$Fr_\mathrm{s} = 0.80$','$Fr_\mathrm{s} = 0.71$','$Fr_\mathrm{s} = 0.60$','$Fr_\mathrm{s} = 0.44$','$Fr_\mathrm{s} = 0.35$','location','south','FontName','Times New Roman','Interpreter','Latex')
xlabel('$\lambda$','Interpreter','latex')
ylabel('$C_P$','Interpreter','latex')
set(gca,'FontName','Times New Roman','FontSize',16)
text(0.6,0.3415,'(a)','FontName','Times New Roman','FontSize',16)

s2 = subplot(2,2,2);
p1 = plot(TSR.Baseline,CP.Baseline,'-o','LineWidth',2,'Color',cmp(1,:),'MarkerEdgeColor','none','MarkerFaceColor',cmp(1,:),'MarkerSize',8);
hold on
p2 = plot(TSR.Froude_d0,CP.Froude_d0,'-^','LineWidth',2,'Color',cmp(4,:),'MarkerEdgeColor','none','MarkerFaceColor',cmp(4,:),'MarkerSize',8);
p3 = plot(TSR.Froude_d0_High,CP.Froude_d0_High,'-x','LineWidth',2,'Color',cmpW(3,:),'MarkerEdgeColor',cmpW(3,:),'MarkerFaceColor',cmpW(3,:),'MarkerSize',10);
xlim([0.55 3.2])
ylim([0 0.32])
grid on
legend([p3 p2 p1],'$Fr_0 = 0.39$','$Fr_0 = 0.32$','$Fr_0 = 0.25$','location','south','FontName','Times New Roman','Interpreter','Latex')
xlabel('$\lambda$','Interpreter','latex')
ylabel('$C_P$','Interpreter','latex')
set(gca,'FontName','Times New Roman','FontSize',16)
text(0.6,0.3415,'(b)','FontName','Times New Roman','FontSize',16)

f5.Position = [317 778.3 1039.3 931.3];
s1.Position = [0.0800 0.5838 0.4047 0.3412];
s2.Position = [0.5803 0.5838 0.4047 0.3412];