% Apply Blockage Corrections
% Hannah Ross

% Applies the Barnsley and Wellicome, Mikkelsen and Sorensen, Werle, and
% Houlsby et al. blockage corrections to performance data from a cross-flow
% turbine and an axial-flow turbine taken at high blockage. 

clear

%% Define constants
g = 9.81;

%% Load data
load Confined.mat
load Unconfined.mat

%% Calculate non-dimensional parameters
fieldNames = fieldnames(Confined);
for i = 1:length(fieldNames)
    Confined.(fieldNames{i}).beta = (Confined.(fieldNames{i}).AT+Confined.(fieldNames{i}).AS)./Confined.(fieldNames{i}).AC;  % equation (1.1)
    Unconfined.(fieldNames{i}).beta = (Unconfined.(fieldNames{i}).AT+Unconfined.(fieldNames{i}).AS)./Unconfined.(fieldNames{i}).AC;  % equation (1.1)
    Confined.(fieldNames{i}).Fr = mean(Confined.(fieldNames{i}).V0./(sqrt(g.*Confined.(fieldNames{i}).d0)));  % equation (2.2)
    Unconfined.(fieldNames{i}).Fr = mean(Unconfined.(fieldNames{i}).V0./(sqrt(g.*Unconfined.(fieldNames{i}).d0)));  % equation (2.2)
end

%% Apply Barnsley and Wellicome
u2u1Guess = [1.01 10];
for i = 1:length(fieldNames)
    for j = 1:length(Confined.(fieldNames{i}).CT)
        func = @(u2u1) u2u1-Confined.(fieldNames{i}).beta.*(-1+sqrt(1+Confined.(fieldNames{i}).beta.*(u2u1^2-1)))./ ...  % iterative solution to equation (3.9)
            (Confined.(fieldNames{i}).beta.*(u2u1-1)).*(u2u1-1)-(1./sqrt(Confined.(fieldNames{i}).CT(j)./(u2u1^2-1)));
        [BW.(fieldNames{i}).u2u1(j,1),BW.(fieldNames{i}).diff(j,1),BW.(fieldNames{i}).eFlag(j,1),~] = fzero(func,u2u1Guess);
        clear func
    end
    BW.(fieldNames{i}).uTu1(:,1) = (-1+sqrt(1+Confined.(fieldNames{i}).beta.*(BW.(fieldNames{i}).u2u1.^2-1)))./(Confined.(fieldNames{i}).beta.*(BW.(fieldNames{i}).u2u1-1));  % equation (3.9a)
    BW.(fieldNames{i}).u1V0(:,1) = sqrt(Confined.(fieldNames{i}).CT./(BW.(fieldNames{i}).u2u1.^2-1));  % equation (3.9c)
    BW.(fieldNames{i}).uTV0(:,1) = BW.(fieldNames{i}).uTu1.*BW.(fieldNames{i}).u1V0;  % equation (3.10)
    BW.(fieldNames{i}).V0(:,1) = Confined.(fieldNames{i}).V0.*((BW.(fieldNames{i}).uTV0.^2)+Confined.(fieldNames{i}).CT/4)./BW.(fieldNames{i}).uTV0;  % equation (3.8)
    BW.(fieldNames{i}).CP(:,1) = Confined.(fieldNames{i}).CP.*((Confined.(fieldNames{i}).V0./BW.(fieldNames{i}).V0).^3);  % equation (3.5a)
    BW.(fieldNames{i}).CT(:,1) = Confined.(fieldNames{i}).CT.*((Confined.(fieldNames{i}).V0./BW.(fieldNames{i}).V0).^2);  % equation (3.5b)
    BW.(fieldNames{i}).lambda(:,1) = Confined.(fieldNames{i}).lambda.*(Confined.(fieldNames{i}).V0./BW.(fieldNames{i}).V0);  % equation (3.5c)
end
clear u2u1Guess

%% Apply Mikkelsen & Sorensen
for i = 1:length(fieldNames)
    if contains(fieldNames{i},'Wake')
        MS.(fieldNames{i}).uT(:,1) = (Confined.(fieldNames{i}).A1./Confined.(fieldNames{i}).AT.*(Confined.(fieldNames{i}).beta.*(Confined.(fieldNames{i}).A1./Confined.(fieldNames{i}).AT).^2-1))./ ...
            (Confined.(fieldNames{i}).A1./Confined.(fieldNames{i}).AT.*Confined.(fieldNames{i}).beta.*(3*Confined.(fieldNames{i}).A1./Confined.(fieldNames{i}).AT-2)-2*Confined.(fieldNames{i}).A1./Confined.(fieldNames{i}).AT+1);  % equation (3.11a)
        MS.(fieldNames{i}).u1(:,1) = MS.(fieldNames{i}).uT./(Confined.(fieldNames{i}).A1./Confined.(fieldNames{i}).AT);  % equation (3.11b)
        MS.(fieldNames{i}).u2(:,1) = (1-Confined.(fieldNames{i}).beta.*MS.(fieldNames{i}).uT)./(1-Confined.(fieldNames{i}).beta.*Confined.(fieldNames{i}).A1./Confined.(fieldNames{i}).AT);  % equation (3.11c)
        MS.(fieldNames{i}).CT(:,1) = MS.(fieldNames{i}).u2.^2-MS.(fieldNames{i}).u1.^2;  % equation (3.11d)
        MS.(fieldNames{i}).V0(:,1) = mean(Confined.(fieldNames{i}).V0).*(MS.(fieldNames{i}).uT+1/4*MS.(fieldNames{i}).CT./MS.(fieldNames{i}).uT);  % equation (3.8)
        MS.(fieldNames{i}).CP(:,1) = Confined.(fieldNames{i}).CP.*(Confined.(fieldNames{i}).V0./MS.(fieldNames{i}).V0).^3;  % equation (3.5a)
        MS.(fieldNames{i}).CT(:,1) = Confined.(fieldNames{i}).CT.*(Confined.(fieldNames{i}).V0./MS.(fieldNames{i}).V0).^2;  % equation (3.5b)
        MS.(fieldNames{i}).lambda(:,1) = Confined.(fieldNames{i}).lambda.*(Confined.(fieldNames{i}).V0./MS.(fieldNames{i}).V0);  % equation (3.5c)
    end
end

%% Apply Werle
for i = 1:length(fieldNames)
        Werle.(fieldNames{i}).CP(:,1) = Confined.(fieldNames{i}).CP.*(1-Confined.(fieldNames{i}).beta).^2;  % equation (3.15a)
        Werle.(fieldNames{i}).CT(:,1) = Confined.(fieldNames{i}).CT.*(1-Confined.(fieldNames{i}).beta).^2./(1+Confined.(fieldNames{i}).beta);  % equation (3.15b)
        Werle.(fieldNames{i}).lambda(:,1) = Confined.(fieldNames{i}).lambda.*(1-Confined.(fieldNames{i}).beta);  % equation (3.16)
end

%% Apply Houlsby et al. and Maskell
for i = 1:length(fieldNames)
    for j = 1:length(Confined.(fieldNames{i}).CT)
        u2Guess = BW.(fieldNames{i}).u2u1(j,1).*BW.(fieldNames{i}).u1V0(j,1).*Confined.(fieldNames{i}).V0(j,1);
        func = @(u2) (Confined.(fieldNames{i}).Fr.^2.*u2.^4-(4+2*Confined.(fieldNames{i}).Fr.^2).*Confined.(fieldNames{i}).V0(j,1).^2.*u2.^2+ ... % iterative solution to equation (3.18)
            8*Confined.(fieldNames{i}).V0(j,1).^3.*u2-4*Confined.(fieldNames{i}).V0(j,1).^4+4*Confined.(fieldNames{i}).beta.* ...
            Confined.(fieldNames{i}).CT(j,1).*Confined.(fieldNames{i}).V0(j,1).^4+Confined.(fieldNames{i}).Fr.^2.* ...
            Confined.(fieldNames{i}).V0(j,1).^4)./(-4*Confined.(fieldNames{i}).Fr.^2.*u2.^3+8*Confined.(fieldNames{i}).V0(j,1).^2.*u2+ ...
            4*Confined.(fieldNames{i}).Fr.^2.*Confined.(fieldNames{i}).V0(j,1).^2.*u2-8*Confined.(fieldNames{i}).V0(j,1).^3)- ...
            ((u2.^2-Confined.(fieldNames{i}).CT(j,1).*Confined.(fieldNames{i}).V0(j,1).^2)^(1/2));
        [Houlsby.(fieldNames{i}).u2(j,1),Houlsby.(fieldNames{i}).diff(j,1),Houlsby.(fieldNames{i}).eFlag(j,1),~] = fzero(func,u2Guess);
        clear u2Guess func
        Houlsby.(fieldNames{i}).u1(j,1) = (Houlsby.(fieldNames{i}).u2(j,1).^2-Confined.(fieldNames{i}).CT(j,1).*Confined.(fieldNames{i}).V0(j,1).^2)^(1/2);  % equation (3.18b)
    end
    for j = 1:length(Houlsby.(fieldNames{i}).u1)
        Houlsby.(fieldNames{i}).uT(j,1) = (Houlsby.(fieldNames{i}).u1(j,1).*(Houlsby.(fieldNames{i}).u2(j,1)-Confined.(fieldNames{i}).V0(j,1)).*(2*g.*Confined.(fieldNames{i}).d0- ... % equation (3.19)
            Houlsby.(fieldNames{i}).u2(j,1).^2-Houlsby.(fieldNames{i}).u2(j,1).*Confined.(fieldNames{i}).V0(j,1)))./(2*Confined.(fieldNames{i}).beta.*g.* ...
            Confined.(fieldNames{i}).d0.*(Houlsby.(fieldNames{i}).u2(j,1)-Houlsby.(fieldNames{i}).u1(j,1)));
    end
    Houlsby.(fieldNames{i}).V0(:,1) = Confined.(fieldNames{i}).V0.*(Confined.(fieldNames{i}).CT/4+ ...
        (Houlsby.(fieldNames{i}).uT./Confined.(fieldNames{i}).V0).^2)./(Houlsby.(fieldNames{i}).uT./Confined.(fieldNames{i}).V0);  % equation (3.8)
    Houlsby.(fieldNames{i}).CP(:,1) = Confined.(fieldNames{i}).CP.*((Confined.(fieldNames{i}).V0./Houlsby.(fieldNames{i}).V0).^3);  % equation (3.5a)
    Houlsby.(fieldNames{i}).CT(:,1) = Confined.(fieldNames{i}).CT.*((Confined.(fieldNames{i}).V0./Houlsby.(fieldNames{i}).V0).^2);  % equation (3.5b)
    Houlsby.(fieldNames{i}).lambda(:,1) = Confined.(fieldNames{i}).lambda.*(Confined.(fieldNames{i}).V0./Houlsby.(fieldNames{i}).V0);  % equation (3.5c)
    Maskell.(fieldNames{i}).CP(:,1) = Confined.(fieldNames{i}).CP.*((Confined.(fieldNames{i}).V0./Houlsby.(fieldNames{i}).u2).^3);  % equation (5.5)        
    Maskell.(fieldNames{i}).CT(:,1) = Confined.(fieldNames{i}).CT.*((Confined.(fieldNames{i}).V0./Houlsby.(fieldNames{i}).u2).^2);  % equation (5.4a)
    Maskell.(fieldNames{i}).lambda(:,1) = Confined.(fieldNames{i}).lambda.*(Confined.(fieldNames{i}).V0./Houlsby.(fieldNames{i}).u2);  % equation (5.4b)
end

%% Calculate Error

spacing = 0.001;
    % Interpolate CP (CFT)
for i = 1:length(fieldNames)
    if strcmp(fieldNames{i},'CFT')
        ConfinedInterp.(fieldNames{i}).lambda = Confined.(fieldNames{i}).lambda(1):spacing:Confined.(fieldNames{i}).lambda(end);
        UnconfinedInterp.(fieldNames{i}).lambda = Unconfined.(fieldNames{i}).lambda(1):spacing:Unconfined.(fieldNames{i}).lambda(end);
        BWInterp.(fieldNames{i}).lambda = BW.(fieldNames{i}).lambda(1):spacing:BW.(fieldNames{i}).lambda(end);
        HoulsbyInterp.(fieldNames{i}).lambda = Houlsby.(fieldNames{i}).lambda(1):spacing:Houlsby.(fieldNames{i}).lambda(end);
        WerleInterp.(fieldNames{i}).lambda = Werle.(fieldNames{i}).lambda(1):spacing:Werle.(fieldNames{i}).lambda(end);
        MaskellInterp.(fieldNames{i}).lambda = Maskell.(fieldNames{i}).lambda(1):spacing:Maskell.(fieldNames{i}).lambda(end);
        ConfinedInterp.(fieldNames{i}).CP = interp1(Confined.(fieldNames{i}).lambda,Confined.(fieldNames{i}).CP,ConfinedInterp.(fieldNames{i}).lambda,'pchip');
        UnconfinedInterp.(fieldNames{i}).CP = interp1(Unconfined.(fieldNames{i}).lambda,Unconfined.(fieldNames{i}).CP,UnconfinedInterp.(fieldNames{i}).lambda,'pchip');
        BWInterp.(fieldNames{i}).CP = interp1(BW.(fieldNames{i}).lambda,BW.(fieldNames{i}).CP,BWInterp.(fieldNames{i}).lambda,'pchip');
        HoulsbyInterp.(fieldNames{i}).CP = interp1(Houlsby.(fieldNames{i}).lambda,Houlsby.(fieldNames{i}).CP,HoulsbyInterp.(fieldNames{i}).lambda,'pchip');
        WerleInterp.(fieldNames{i}).CP = interp1(Werle.(fieldNames{i}).lambda,Werle.(fieldNames{i}).CP,WerleInterp.(fieldNames{i}).lambda,'pchip');
        MaskellInterp.(fieldNames{i}).CP = interp1(Maskell.(fieldNames{i}).lambda,Maskell.(fieldNames{i}).CP,MaskellInterp.(fieldNames{i}).lambda,'pchip');
    end
end
    % Interpolate CP (AFT)
for i = 1:length(fieldNames)
    if strcmp(fieldNames{i},'AFT')
        ConfinedInterp.(fieldNames{i}).lambda = Confined.(fieldNames{i}).lambda(1)*0.3:spacing:Confined.(fieldNames{i}).lambda(end);
        UnconfinedInterp.(fieldNames{i}).lambda = Unconfined.(fieldNames{i}).lambda(1)*0.3:spacing:Unconfined.(fieldNames{i}).lambda(end)*1.1;
        BWInterp.(fieldNames{i}).lambda = BW.(fieldNames{i}).lambda(1)*0.3:spacing:BW.(fieldNames{i}).lambda(end)*1.1;
        HoulsbyInterp.(fieldNames{i}).lambda = Houlsby.(fieldNames{i}).lambda(1)*0.3:spacing:Houlsby.(fieldNames{i}).lambda(end)*1.1;
        WerleInterp.(fieldNames{i}).lambda = Werle.(fieldNames{i}).lambda(1)*0.3:spacing:Werle.(fieldNames{i}).lambda(end)*1.1;
        MaskellInterp.(fieldNames{i}).lambda = Maskell.(fieldNames{i}).lambda(1)*0.3:spacing:Maskell.(fieldNames{i}).lambda(end)*1.1;
        ConfinedInterp.(fieldNames{i}).CP = interp1(Confined.(fieldNames{i}).lambda,Confined.(fieldNames{i}).CP,ConfinedInterp.(fieldNames{i}).lambda,'linear','extrap');
        UnconfinedInterp.(fieldNames{i}).CP = interp1(Unconfined.(fieldNames{i}).lambda,Unconfined.(fieldNames{i}).CP,UnconfinedInterp.(fieldNames{i}).lambda,'linear','extrap');
        BWInterp.(fieldNames{i}).CP = interp1(BW.(fieldNames{i}).lambda,BW.(fieldNames{i}).CP,BWInterp.(fieldNames{i}).lambda,'linear','extrap');
        HoulsbyInterp.(fieldNames{i}).CP = interp1(Houlsby.(fieldNames{i}).lambda,Houlsby.(fieldNames{i}).CP,HoulsbyInterp.(fieldNames{i}).lambda,'linear','extrap');
        WerleInterp.(fieldNames{i}).CP = interp1(Werle.(fieldNames{i}).lambda,Werle.(fieldNames{i}).CP,WerleInterp.(fieldNames{i}).lambda,'linear','extrap');
        MaskellInterp.(fieldNames{i}).CP = interp1(Maskell.(fieldNames{i}).lambda,Maskell.(fieldNames{i}).CP,MaskellInterp.(fieldNames{i}).lambda,'linear','extrap');
    end
end
    % Select TSR range of positive CP
for i = 1:length(fieldNames)
    if strcmp(fieldNames{i},'CFT') || strcmp(fieldNames{i},'AFT')
        ConfinedInterp.(fieldNames{i}).lambda = ConfinedInterp.(fieldNames{i}).lambda(ConfinedInterp.(fieldNames{i}).CP > 0);
        UnconfinedInterp.(fieldNames{i}).lambda = UnconfinedInterp.(fieldNames{i}).lambda(UnconfinedInterp.(fieldNames{i}).CP > 0);
        BWInterp.(fieldNames{i}).lambda = BWInterp.(fieldNames{i}).lambda(BWInterp.(fieldNames{i}).CP > 0);
        HoulsbyInterp.(fieldNames{i}).lambda = HoulsbyInterp.(fieldNames{i}).lambda(HoulsbyInterp.(fieldNames{i}).CP > 0);
        WerleInterp.(fieldNames{i}).lambda = WerleInterp.(fieldNames{i}).lambda(WerleInterp.(fieldNames{i}).CP > 0);
        MaskellInterp.(fieldNames{i}).lambda = MaskellInterp.(fieldNames{i}).lambda(MaskellInterp.(fieldNames{i}).CP > 0);
    end
end
    % Create new TSR values with equal numbers of points
n = 1500;
for i = 1:length(fieldNames)
    if strcmp(fieldNames{i},'CFT') || strcmp(fieldNames{i},'AFT')
        ConfinedInterp.(fieldNames{i}).lambda = linspace(ConfinedInterp.(fieldNames{i}).lambda(1),ConfinedInterp.(fieldNames{i}).lambda(end),n);
        UnconfinedInterp.(fieldNames{i}).lambda = linspace(UnconfinedInterp.(fieldNames{i}).lambda(1),UnconfinedInterp.(fieldNames{i}).lambda(end),n);
        BWInterp.(fieldNames{i}).lambda = linspace(BWInterp.(fieldNames{i}).lambda(1),BWInterp.(fieldNames{i}).lambda(end),n);
        HoulsbyInterp.(fieldNames{i}).lambda = linspace(HoulsbyInterp.(fieldNames{i}).lambda(1),HoulsbyInterp.(fieldNames{i}).lambda(end),n);
        WerleInterp.(fieldNames{i}).lambda = linspace(WerleInterp.(fieldNames{i}).lambda(1),WerleInterp.(fieldNames{i}).lambda(end),n);
        MaskellInterp.(fieldNames{i}).lambda = linspace(MaskellInterp.(fieldNames{i}).lambda(1),MaskellInterp.(fieldNames{i}).lambda(end),n);
    end
end
    % Interpolate CP and CT using new TSR ranges (CFT)
errorNames = {'CT','CP'};
for i = 1:length(fieldNames)
    if strcmp(fieldNames{i},'CFT')
        for j = 1:length(errorNames)
            ConfinedInterp.(fieldNames{i}).(errorNames{j}) = interp1(Confined.(fieldNames{i}).lambda,Confined.(fieldNames{i}).(errorNames{j}),ConfinedInterp.(fieldNames{i}).lambda,'pchip');
            UnconfinedInterp.(fieldNames{i}).(errorNames{j}) = interp1(Unconfined.(fieldNames{i}).lambda,Unconfined.(fieldNames{i}).(errorNames{j}),UnconfinedInterp.(fieldNames{i}).lambda,'pchip');
            BWInterp.(fieldNames{i}).(errorNames{j}) = interp1(BW.(fieldNames{i}).lambda,BW.(fieldNames{i}).(errorNames{j}),BWInterp.(fieldNames{i}).lambda,'pchip');
            HoulsbyInterp.(fieldNames{i}).(errorNames{j}) = interp1(Houlsby.(fieldNames{i}).lambda,Houlsby.(fieldNames{i}).(errorNames{j}),HoulsbyInterp.(fieldNames{i}).lambda,'pchip');
            WerleInterp.(fieldNames{i}).(errorNames{j}) = interp1(Werle.(fieldNames{i}).lambda,Werle.(fieldNames{i}).(errorNames{j}),WerleInterp.(fieldNames{i}).lambda,'pchip');
            MaskellInterp.(fieldNames{i}).(errorNames{j}) = interp1(Maskell.(fieldNames{i}).lambda,Maskell.(fieldNames{i}).(errorNames{j}),MaskellInterp.(fieldNames{i}).lambda,'pchip');
        end
    end
end
    % Interpolate CP and CT using new TSR ranges (AFT)
errorNames = {'CT','CP'};
for i = 1:length(fieldNames)
    if strcmp(fieldNames{i},'AFT')
        for j = 1:length(errorNames)
            ConfinedInterp.(fieldNames{i}).(errorNames{j}) = interp1(Confined.(fieldNames{i}).lambda,Confined.(fieldNames{i}).(errorNames{j}),ConfinedInterp.(fieldNames{i}).lambda,'linear','extrap');
            UnconfinedInterp.(fieldNames{i}).(errorNames{j}) = interp1(Unconfined.(fieldNames{i}).lambda,Unconfined.(fieldNames{i}).(errorNames{j}),UnconfinedInterp.(fieldNames{i}).lambda,'linear','extrap');
            BWInterp.(fieldNames{i}).(errorNames{j}) = interp1(BW.(fieldNames{i}).lambda,BW.(fieldNames{i}).(errorNames{j}),BWInterp.(fieldNames{i}).lambda,'linear','extrap');
            HoulsbyInterp.(fieldNames{i}).(errorNames{j}) = interp1(Houlsby.(fieldNames{i}).lambda,Houlsby.(fieldNames{i}).(errorNames{j}),HoulsbyInterp.(fieldNames{i}).lambda,'linear','extrap');
            WerleInterp.(fieldNames{i}).(errorNames{j}) = interp1(Werle.(fieldNames{i}).lambda,Werle.(fieldNames{i}).(errorNames{j}),WerleInterp.(fieldNames{i}).lambda,'linear','extrap');
            MaskellInterp.(fieldNames{i}).(errorNames{j}) = interp1(Maskell.(fieldNames{i}).lambda,Maskell.(fieldNames{i}).(errorNames{j}),MaskellInterp.(fieldNames{i}).lambda,'linear','extrap');
        end
    end
end

% Calculate mean percentage error (MPE)
errorNames = {'lambda','CT','CP'};
for i = 1:length(fieldNames)
    if strcmp(fieldNames{i},'CFT') || strcmp(fieldNames{i},'AFT')
        for j = 1:length(errorNames)        
            ConfinedError.(fieldNames{i}).(errorNames{j}) = mean(abs(UnconfinedInterp.(fieldNames{i}).(errorNames{j})-ConfinedInterp.(fieldNames{i}).(errorNames{j}))./abs(UnconfinedInterp.(fieldNames{i}).(errorNames{j})))*100;
            BWError.(fieldNames{i}).(errorNames{j}) = mean(abs(UnconfinedInterp.(fieldNames{i}).(errorNames{j})-BWInterp.(fieldNames{i}).(errorNames{j}))./abs(UnconfinedInterp.(fieldNames{i}).(errorNames{j})))*100;
            HoulsbyError.(fieldNames{i}).(errorNames{j}) = mean(abs(UnconfinedInterp.(fieldNames{i}).(errorNames{j})-HoulsbyInterp.(fieldNames{i}).(errorNames{j}))./abs(UnconfinedInterp.(fieldNames{i}).(errorNames{j})))*100;
            WerleError.(fieldNames{i}).(errorNames{j}) = mean(abs(UnconfinedInterp.(fieldNames{i}).(errorNames{j})-WerleInterp.(fieldNames{i}).(errorNames{j}))./abs(UnconfinedInterp.(fieldNames{i}).(errorNames{j})))*100;
            MaskellError.(fieldNames{i}).(errorNames{j}) = abs(UnconfinedInterp.(fieldNames{i}).(errorNames{j})-MaskellInterp.(fieldNames{i}).(errorNames{j}))./abs(UnconfinedInterp.(fieldNames{i}).(errorNames{j}));
        end
    end
end

% Calculate wake mean percentage error (MPE)
for i = 1:length(fieldNames)
    if strcmp(fieldNames{i},'CFTWake') || strcmp(fieldNames{i},'AFTWake')
        for j = 1:length(errorNames)
            ConfinedError.(fieldNames{i}).(errorNames{j}) = abs(Unconfined.(fieldNames{i}).(errorNames{j})-Confined.(fieldNames{i}).(errorNames{j}))./abs(Unconfined.(fieldNames{i}).(errorNames{j}))*100;
            BWError.(fieldNames{i}).(errorNames{j}) = abs(Unconfined.(fieldNames{i}).(errorNames{j})-BW.(fieldNames{i}).(errorNames{j}))./abs(Unconfined.(fieldNames{i}).(errorNames{j}))*100;
            HoulsbyError.(fieldNames{i}).(errorNames{j}) = abs(Unconfined.(fieldNames{i}).(errorNames{j})-Houlsby.(fieldNames{i}).(errorNames{j}))./abs(Unconfined.(fieldNames{i}).(errorNames{j}))*100;
            WerleError.(fieldNames{i}).(errorNames{j}) = abs(Unconfined.(fieldNames{i}).(errorNames{j})-Werle.(fieldNames{i}).(errorNames{j}))./abs(Unconfined.(fieldNames{i}).(errorNames{j}))*100;        
            MSError.(fieldNames{i}).(errorNames{j}) = abs(Unconfined.(fieldNames{i}).(errorNames{j})-MS.(fieldNames{i}).(errorNames{j}))./abs(Unconfined.(fieldNames{i}).(errorNames{j}))*100;
        end
    end
end

% A1 at 2.25D gives lowest error for CFT devices and A1 at 1.75D gives
% lowest error for AFT devices