warning off
figure(1);clf
figure(2);clf;hold on;
figure(3);clf;hold on;

% Define COVs from analysis results
sCOV{1}='6.3%';
sCOV{2}='6.9%';
sCOV{3}='8.6%';
sCOV{4}='8.7%';
sCOV{5}='9.4%';
sCOV{6}='14.9%';

% Iterate on each of 3 results files
N=6;
for i=1:N

    % Load the data
    if i==1
        sFile='(a) PP1, Capillary, Iso, U4';
        sStyle='kv';
    elseif i==2
        sFile='(b) PP2, Capillary, Iso, U4';
        sStyle='kx';
    elseif i==3
        sFile='(c) PP2, Parallel, Iso, U4';
        sStyle='ks';
    elseif i==4
        sFile='(d) PP2, Parallel, Non-iso, H500, U4';
        sStyle='ro';
    elseif i==5
        sFile='(e) PP2, Parallel, Non-iso, H50, U4';
        sStyle='rx';
    elseif i==6
        sFile='(f) PP2, Parallel, Non-iso, H50, U20';
        sStyle='r^';
    elseif i==7
        sFile='Outlet Velocities, PP2-Capillary (1)';
        sStyle='r-';
    %elseif i==4
    %    sFile='Newtonian Outlet Velocities';
    %    sStyle='r-';
    end
    data=load([sFile '.txt']);

    % Extract the columns, scale to mm, and manage orientation
    x = 1000*flipud(data(:,2));
    y = 1000*data(:,1);
    v = 1000*flipud(data(:,3));
    v_mean=mean(mean(v));
    v=v/v_mean;
    y = y-min(y);

    % Mirror the data about x=0 and y=0
    x_mirror = [x; x; ];
    y_mirror = [y; -y; ];
    v_mirror = [v; v; ];
    
    % Create a grid for contour plots
    grid_x = unique(x_mirror);
    grid_y = unique(y_mirror);
    [grid_X, grid_Y] = meshgrid(grid_x, grid_y);
    
    % Interpolate the velocities onto the grid
    grid_V = griddata(x_mirror, y_mirror, v_mirror, grid_X, grid_Y, 'nearest');
    
    % Calculate the areal flow rate, mm^2/s, across the y direction (gap) using trapz
    for j=1:size(grid_V,2)
        v_sum(j,i) = trapz(grid_Y(:,j),grid_V(:,j));
    end

    V_COV(i)=std(v_sum(:,i))/mean(v_sum(:,i));
    H_Mean(i)=mean(v_sum(:,i));
    norm_V=grid_V/H_Mean(i);

    % Create the contour plot
    %subplot(3,2,iPlot)
    %contourf(grid_X, grid_Y, grid_V);
    figure(1)
    subplot(3,6,3*i-1)
    %contourf(grid_Y, grid_X, norm_V);
    contourf(grid_Y, grid_X, grid_V, 15);
    colormap('parula');
    colormap('jet');
    %colorbar;
    %title(sFile);
    caxis([0.5 1.5]);  % Set the color scale
    axis off

    % Plot the total flow rate as a function of x
    figure(2);
    plot(grid_x, v_sum(:,i),sStyle);
    xlabel('Position across die width, mm');
    ylabel('Areal flow rate, mm^2/s');

    figure(3);hold on
    plot(grid_x, 100*v_sum(:,i)/mean(v_sum(:,i)),sStyle,'DisplayName',[sFile ', COV=' sCOV{i}]);
    xlabel('Position across die width, mm');
    ylabel('Relative thickness (% mean)');

end

grid on;
legend('location','south')
V_COV
H_Mean