%% Step 3: Assemble figures 

% Author: Eleanor Fadely, University of California Davis
% Last updated: 01/27/2026

% List of figures and associated processing:

% Figure 3: Color brightfield images of reactor inlet showing irradiation setup
    % (1) Load color brightfield images.
    % (2) Rotate and crop images.
    % (3) Save images.

% Figure 4: Schematic of image analysis steps
    % (1) Load and crop representative color brightfield images, binary
    % masks, and intensity-corrected matrix. 
    % (2) Mask cropped brightfield images for visualization. 
    % (3) Apply a colormap to the intensity-corrected matrix for
    % visualization. 
    % (4) Save images.

% Figure 5a-c: Reactor-scale mineral accumulation at the reference and final timepoints
    % (1) Load color brightfield images from the reference timepoint and
    % final timepoint, and intensity-corrected matrix from the final
    % timepoint.
    % (2) Apply a colormap to the intensity-corrected matrix for
    % visualization.
    % (3) Save images.

% Figure 5d: Mineral distribution along longitudinal dimension of pore space  
    % (1) Load the intensity-corrected matrix for the final timepoint
    % (triplicate experiments).
    % (2) Sum along the columns of each matrix to produce a row vector
    % (3) Calculate the moving average over 1 mm.
    % (4) Average the triplicate experiments and calculate the standard
    % deviation.
    % (5) Plot the total corrected intensity over distance.
    % (6) Save figure.

% Figure 6: Biofilm-scale mineral accumulation as a function of time
    % (1) Load representative color brightfield images and
    % intensity-corrected matrices.
    % (2) Crop color brightfield images.
    % (3) Apply colormaps to the intensity-corrected matrices for
    % visualization and crop.
    % (4) Save images.

% Figure 7: Rate of Mn oxide precipitation as a function of time
    % (1) Load spreadsheet of the endpoint Mn mass measured by ICP-MS
    % (triplicate experiments).
    % (2) Load vectors of the total corrected intensity over time 
    % (triplicate experiments).
    % (3) Rescale each vector by the final timepoint intensity.
    % (4) Multiply rescaled vector by the endpoint mass.
    % (5) Calculate the difference in mass between sequential timepoints
    % and divide by the time interval between images.
    % (6) Plot the rate over time.
    % (7) Save figure.

clear all; clc; close all;

%% Figure 3a, 3b, 3c: Color brightfield images of reactor inlet showing irradiation setup

% Load, rotate, and crop representative sample images

% Proper UV treatment (Figure 3a)
cd('C:\Users\ecfadely\OneDrive - University of California, Davis\Manuscripts\Mn Biomineralization Protocol\Figures\Figure_3\Proper_UV');
GUV = imread('WholeReactor__Color_XY1.tif'); figure; imshow(GUV);
GUV_Rotate = imrotate(GUV, 2); figure; imshow(GUV_Rotate);
GUV_RotCrop = imcrop(GUV_Rotate, [2096.5 1959.5 11104 6793]); figure; imshow(GUV_RotCrop);

% Incomplete UV treatment (Figure 3b)
cd('C:\Users\ecfadely\OneDrive - University of California, Davis\Manuscripts\Mn Biomineralization Protocol\Figures\Figure_3\Incomplete_UV');
InCompUV = imread('WholeReactor__Color_XY1.tif'); figure; imshow(InCompUV);
ICUV_Rotate = imrotate(InCompUV, 0.25); figure; imshow(ICUV_Rotate);
ICUV_RotCrop = imcrop(ICUV_Rotate, [1366.5 1054.5 11093 7085]); figure; imshow(ICUV_RotCrop);

% No UV treatment (Figure 3c)
cd('C:\Users\ecfadely\OneDrive - University of California, Davis\Manuscripts\Mn Biomineralization Protocol\Figures\Figure_3\No_UV');
NUV = imread('WholeReactor__Color_XY1.tif'); figure; imshow(NUV);
NUV_Rotate = imrotate(NUV, -0.8); figure; imshow(NUV_Rotate);
NUV_RotCrop = imcrop(NUV_Rotate, [1860.5 1305.5 11021 6743]); figure; imshow(NUV_RotCrop);

% Save images
cd('C:\Users\ecfadely\OneDrive - University of California, Davis\Manuscripts\Mn Biomineralization Protocol\Figures\Figure_3');
imwrite(GUV_RotCrop, 'Rotated_Cropped_Inlet_GoodUV.tif');
imwrite(ICUV_RotCrop, 'Rotated_Cropped_Inlet_IncompleteUV.tif');
imwrite(NUV_RotCrop, 'Rotated_Cropped_Inlet_NoUV.tif');

%% Load variables for Figure 4, Figure 5, Figure 6, Figure 7

clear all; clc; close all;

% Load color brightfield image structure for one representative replicate
cd('F:\240919_Constit_mCh\Oxidation\Protocol_Variables\Processed_Images');
load('CroppedOxidationPhase.mat');

% Load grain mask for one representative replicate
cd('F:\240919_Constit_mCh\Growth\Protocol_Variables\Grain_Mask');
load('GrainMask.mat');

% Load pixel-intensity corrected (subtracted) images

% Replicate 1
cd('F:\240724_Constit_mCh\Oxidation\Protocol_Variables\Image_Subtraction');
load('Subtracted_Images_Ref.mat');
FinalTP_R1 = subtractedGray{1,58};
GrayMass1 = struct2array(load('Grayscale_Total_Amount.mat'));

% Replicate 3
clearvars subtractedGray 
cd('F:\240923_Constit_mCh\Oxidation\Protocol_Variables\Image_Subtraction');
load('Subtracted_Images_Ref.mat');
FinalTP_R3 = subtractedGray{1,58};
GrayMass3 = struct2array(load('Grayscale_Total_Amount.mat'));

% Replicate 2 (representative replicate in figures)
clearvars subtractedGray 
cd('F:\240919_Constit_mCh\Oxidation\Protocol_Variables\Image_Subtraction');
load('Subtracted_Images_Ref.mat');
FinalTP_R2 = subtractedGray{1,58};
GrayMass2 = struct2array(load('Grayscale_Total_Amount.mat'));

% Load image masks for one representative replicate
load('Brightfield_Image_Masks.mat');

%% Figure 4: Schematic of image analysis steps

addpath('C:\Users\ecfadely\OneDrive - University of California, Davis\Matlab_May2024\MatLab\Functions\Colormaps');

close all;

% Assign select timepoints
ReferenceImage = ProcOxidationBF{1,8};
FinalImage = ProcOxidationBF{1,58};
FinalMask = ImageMask{1,58};

% Dilate mask by one pixel for visualization
SE = strel('square', 1);
FMDilate = imdilate(FinalMask, SE); 

% Multiply brightfield images by mask 
MaskedRef = uint16(FMDilate) .* ReferenceImage;
MaskedFinal = uint16(FMDilate) .* FinalImage;

% Crop images
% Cropping coordinates
CropCoord = [1103.5 4447.5 554 555];

% Cropped panels
CropRef = imcrop(ReferenceImage, CropCoord); f1 = figure(1); imshow(CropRef);
CropFinal = imcrop(FinalImage, CropCoord); f2 = figure(2); imshow(CropFinal);
CropMask = imcrop(FinalMask, CropCoord); f3 = figure(3); imshow(CropMask);
CropMaskedRef = imcrop(MaskedRef, CropCoord); f4 = figure(4); imshow(CropMaskedRef);
CropMaskedFinal = imcrop(MaskedFinal, CropCoord); f5 = figure(5); imshow(CropMaskedFinal);

% Compile colormap for subtracted image 

GrainMask = imcomplement(GM_Final);
BiomassRef = imcomplement(ImageMask{1,8});
MaskedBiomassRef = BiomassRef .* GrainMask;
OxidesImageFinal = subtractedGray{1,58};

% Index map of grains (= 2), pores (= 1), and biofilms (= 0)
GrainMaskIndex = GM_Final .* 2; 
indBio = find(MaskedBiomassRef == 1);
BiofilmGrains = GrainMaskIndex;
BiofilmGrains(indBio) = 1; 

% Index masks to distinguish grains, pores, biofilms
MaskedBiomassRefFinal = ImageMask{1,58} .* GrainMask;

% Final timepoint (t = 43.5 h)
indBioFinal = find(MaskedBiomassRefFinal == 1);
BiofilmGrainsFinal = GrainMaskIndex;
BiofilmGrainsFinal(indBioFinal) = 1; 

OxidesImage58 = double(OxidesImageFinal);
indNan58 = find(OxidesImage58 == 0);
OxidesImage58(indNan58) = NaN;

% Crop all select timepoints
CropGrainsFinal = imcrop(BiofilmGrainsFinal, CropCoord);
CropOxidesImageFinal = imcrop(OxidesImage58, CropCoord);

% t = 43.5 h (final timepoint)

% Create colormap for grains, pores, and biofilms
pp = [0.98 0.98 0.98];   % pore color
gg = [0.75 0.75 0.75];   % grain color
bb = [0.88 0.88 0.88];   % biofilm color

CM1 = [pp; bb; gg];
f6 = figure(6);
ax1 = axes;
imagesc(CropGrainsFinal);
view(2)
ax2 = axes;
imagesc(CropOxidesImageFinal, 'AlphaData',~isnan(CropOxidesImageFinal));
linkaxes([ax1,ax2]);
ax1.Visible = 'off';
ax2.Visible = 'off';
ax2.XTick = [];
ax2.YTick = [];
colormap(ax1, CM1);
colormap(ax2, flipud(viridis));
clim([min(min(OxidesImageFinal)) max(max(OxidesImageFinal))]);
set([ax1, ax2],'Position',[.17 .11 .685 .815]);

% cb2 = colorbar(ax2,'Position', [.9 .11 .03 .815]);
% cb2.FontSize = 18;
% cb2.Limits = [min(min(OxidesImage43)) max(max(OxidesImage43))];
% ylabel(cb2, 'Corrected pixel intensity (a.u.)', 'FontSize', 18);

% Save images/figures
cd('C:\Users\ecfadely\OneDrive - University of California, Davis\Manuscripts\Mn Biomineralization Protocol\Figures\Figure_4');
imwrite(CropRef, 'CroppedReferenceImage.tif');
imwrite(CropFinal, 'CroppedFinalImage.tif');
imwrite(CropMask, 'CroppedBinaryMask.tif');
imwrite(CropMaskedRef, 'CroppedMaskedReference.tif');
imwrite(CropMaskedFinal, 'CroppedMaskedFinal.tif');
exportgraphics(f6, 'Cropped_Subtracted_Matrix_Final.tif', 'Resolution', 1500);

%% Figure 5a, 5b, 5c: Reactor-scale mineral accumulation at the reference and final timepoints

close all;

% Plot colormap of final timepoint
f1 = figure(1);
set(gcf, 'Position', [100 100 1700 1000]);

CM11 = [bb; pp; gg]; 
ax1 = axes;
imagesc(BiofilmGrains);
view(2)
ax2 = axes;
imagesc(OxidesImage58, 'AlphaData', ~isnan(OxidesImage58));
linkaxes([ax1,ax2]);
ax1.Visible = 'off';
ax2.Visible = 'off';
ax2.XTick = [];
ax2.YTick = [];
colormap(ax1, CM11);
colormap(ax2, flipud(viridis));
clim([min(min(OxidesImageFinal)) max(max(OxidesImageFinal))]);

set([ax1, ax2],'Position',[.17 .11 .685 .815]);
cb2 = colorbar(ax2,'Position', [.12 .11 .015 .815]);
cb2.FontSize = 22;
cb2.Ticks = [0 1E4 2E4 3E4 4E4 5E4];
cb2.Limits = [min(min(OxidesImageFinal)) max(max(OxidesImageFinal))];
ylabel(cb2, 'Corrected pixel intensity (a.u.)', 'FontSize', 24);

% Save figures
cd('C:\Users\ecfadely\OneDrive - University of California, Davis\Manuscripts\Mn Biomineralization Protocol\Figures\Figure_5');
imwrite(ReferenceImage, 'Reference_PoreSpace.tif');
imwrite(FinalImage, 'Final_PoreSpace.tif');
exportgraphics(f1, 'WholePoreSpace_Final_Colormap.tif', 'Resolution', 1500);

%% Figure 5d: Mineral distribution along longitudinal dimension of pore space

close all;

% Sum final intensity-corrected matrix along the columns to produce 1D longitudinal signal
MassLong_R1 = sum(FinalTP_R1, 1);
MassLong_R2 = sum(FinalTP_R2, 1);
MassLong_R3 = sum(FinalTP_R3, 1);
MLMat = cat(1, MassLong_R1, MassLong_R2, MassLong_R3);

% Calculate a 1 mm moving average for each replicate signal (1 mm = 617 pixels at 1.62 um/pixel)
for j = 1:3
    MAX = movmean(MLMat(j,:), 617);
    MovAvgMassX(j, :) = MAX;
end

% Average replicate signals and calculate standard deviation
AvgMassLong = mean(MovAvgMassX);
STDMassLong = std(MovAvgMassX);

xvals = 0:1:length(AvgMassLong)-1;
xdim = (1.62 * xvals)/1000;

curve1Long = AvgMassLong - STDMassLong;
curve2Long = AvgMassLong + STDMassLong;

% Plot color
GRC = [0.0667    0.2118    0.5020];

f1 = figure(1);
set(gcf,'Position',[100 100 1200 220]);
hold on;
patch([xdim, flip(xdim)], [curve2Long, flip(curve1Long)], GRC, 'FaceAlpha', 0.2, 'EdgeColor','none');
pp1 = plot(xdim, AvgMassLong, '-', 'LineWidth', 2, 'Color', GRC, 'DisplayName', 'Composite');
xlabel('Distance along pore space (mm)');
ylabel(['Corrected pixel' newline 'intensity (a.u.)']);
ylim([0 4E7]);
xlim([0 15]); xticks([0 3 6 9 12 15]);
set(gca, 'FontSize', 17, 'LineWidth', 1.25);
box on;

% Save figure
cd('C:\Users\ecfadely\OneDrive - University of California, Davis\Manuscripts\Mn Biomineralization Protocol\Figures\Figure_5');
saveas(f1, 'Longitudinal_Profile.svg');

%% Figure 6a and 6b: Biofilm-scale mineral accumulation as a function of time

close all;

% Assign select timepoints of brightfield images
T19 = ProcOxidationBF{1,9};
T20 = ProcOxidationBF{1,11};
T21 = ProcOxidationBF{1,13};
T23 = ProcOxidationBF{1,17};
T33 = ProcOxidationBF{1,37};
T43 = ProcOxidationBF{1,57};

% Crop select timepoints of brightfield images
CropT19 = imcrop(T19, CropCoord); f1 = figure(1); imshow(CropT19);
CropT20 = imcrop(T20, CropCoord); f2 = figure(2); imshow(CropT20);
CropT21 = imcrop(T21, CropCoord); f3 = figure(3); imshow(CropT21);
CropT23 = imcrop(T23, CropCoord); f4 = figure(4); imshow(CropT23);
CropT33 = imcrop(T33, CropCoord); f5 = figure(5); imshow(CropT33);
CropT43 = imcrop(T43, CropCoord); f6 = figure(6); imshow(CropT43);

% Assign select timepoints of subtracted images
Sub19 = subtractedGray{1,9};
Sub20 = subtractedGray{1,11};
Sub21 = subtractedGray{1,13};
Sub23 = subtractedGray{1,17};
Sub33 = subtractedGray{1,37};
Sub43 = subtractedGray{1,57};

% Assign reference images for colormap
GrainMask = imcomplement(GM_Final);
BiomassRef = imcomplement(ImageMask{1,8});
MaskedBiomassRef = BiomassRef .* GrainMask;
OxidesImageFinal = subtractedGray{1,58};

% Index map of grains, pores, and biofilms
GrainMaskIndex = GM_Final .* 2; 
indBio = find(MaskedBiomassRef == 1);
BiofilmGrains = GrainMaskIndex;
BiofilmGrains(indBio) = 1; 

% Index masks to distinguish grains, pores, biofilms
MaskedBiomassRef19 = ImageMask{1,9} .* GrainMask;
MaskedBiomassRef20 = ImageMask{1,11} .* GrainMask;
MaskedBiomassRef21 = ImageMask{1,13} .* GrainMask;
MaskedBiomassRef23 = ImageMask{1,17} .* GrainMask;
MaskedBiomassRef33 = ImageMask{1,37} .* GrainMask;
MaskedBiomassRef43 = ImageMask{1,57} .* GrainMask;

% t = 19 h
indBio19 = find(MaskedBiomassRef19 == 1);
BiofilmGrains19 = GrainMaskIndex;
BiofilmGrains19(indBio19) = 1; 

OxidesImage19 = double(Sub19);
indNan19 = find(OxidesImage19 == 0);
OxidesImage19(indNan19) = NaN;

% t = 20 h
indBio20 = find(MaskedBiomassRef20 == 1);
BiofilmGrains20 = GrainMaskIndex;
BiofilmGrains20(indBio20) = 1; 

OxidesImage20 = double(Sub20);
indNan20 = find(OxidesImage20 == 0);
OxidesImage20(indNan20) = NaN;

% t = 21 h
indBio21 = find(MaskedBiomassRef21 == 1);
BiofilmGrains21 = GrainMaskIndex;
BiofilmGrains21(indBio21) = 1; 

OxidesImage21 = double(Sub21);
indNan21 = find(OxidesImage21 == 0);
OxidesImage21(indNan21) = NaN;

% t = 23 h
indBio23 = find(MaskedBiomassRef23 == 1);
BiofilmGrains23 = GrainMaskIndex;
BiofilmGrains23(indBio23) = 1; 

OxidesImage23 = double(Sub23);
indNan23 = find(OxidesImage23 == 0);
OxidesImage23(indNan23) = NaN;

% t = 33 h
indBio33 = find(MaskedBiomassRef33 == 1);
BiofilmGrains33 = GrainMaskIndex;
BiofilmGrains33(indBio33) = 1; 

OxidesImage33 = double(Sub33);
indNan33 = find(OxidesImage33 == 0);
OxidesImage33(indNan33) = NaN;

% t = 43 h
indBio43 = find(MaskedBiomassRef43 == 1);
BiofilmGrains43 = GrainMaskIndex;
BiofilmGrains43(indBio43) = 1; 

OxidesImage43 = double(Sub43);
indNan43 = find(OxidesImage43 == 0);
OxidesImage43(indNan43) = NaN;

% Crop all select timepoints
CropGrains = imcrop(BiofilmGrains, CropCoord);
CropGrains19 = imcrop(BiofilmGrains19, CropCoord);
CropGrains20 = imcrop(BiofilmGrains20, CropCoord);
CropGrains21 = imcrop(BiofilmGrains21, CropCoord);
CropGrains23 = imcrop(BiofilmGrains23, CropCoord);
CropGrains33 = imcrop(BiofilmGrains33, CropCoord);
CropGrains43 = imcrop(BiofilmGrains43, CropCoord);

CropOxidesImage19 = imcrop(OxidesImage19, CropCoord);
CropOxidesImage20 = imcrop(OxidesImage20, CropCoord);
CropOxidesImage21 = imcrop(OxidesImage21, CropCoord);
CropOxidesImage23 = imcrop(OxidesImage23, CropCoord);
CropOxidesImage33 = imcrop(OxidesImage33, CropCoord);
CropOxidesImage43 = imcrop(OxidesImage43, CropCoord);

% Compile colormaps for each timepoint
% t = 19 h
f7 = figure(7);
ax1 = axes;
imagesc(CropGrains19);
view(2)
ax2 = axes;
imagesc(CropOxidesImage19, 'AlphaData',~isnan(CropOxidesImage19));
linkaxes([ax1,ax2]);
ax1.Visible = 'off';
ax2.Visible = 'off';
ax2.XTick = [];
ax2.YTick = [];
colormap(ax1, CM1);
colormap(ax2, flipud(viridis));
clim([min(min(OxidesImageFinal)) max(max(OxidesImageFinal))]);
set([ax1, ax2],'Position',[.17 .11 .685 .815]);

% t = 20 h
f8 = figure(8);
ax1 = axes;
imagesc(CropGrains20);
view(2)
ax2 = axes;
imagesc(CropOxidesImage20, 'AlphaData',~isnan(CropOxidesImage20));
linkaxes([ax1,ax2]);
ax1.Visible = 'off';
ax2.Visible = 'off';
ax2.XTick = [];
ax2.YTick = [];
colormap(ax1, CM1);
colormap(ax2, flipud(viridis));
clim([min(min(OxidesImageFinal)) max(max(OxidesImageFinal))]);
set([ax1, ax2],'Position',[.17 .11 .685 .815]);

% t = 21 h
f9 = figure(9);
ax1 = axes;
imagesc(CropGrains21);
view(2)
ax2 = axes;
imagesc(CropOxidesImage21, 'AlphaData',~isnan(CropOxidesImage21));
linkaxes([ax1,ax2]);
ax1.Visible = 'off';
ax2.Visible = 'off';
ax2.XTick = [];
ax2.YTick = [];
colormap(ax1, CM1);
colormap(ax2, flipud(viridis));
clim([min(min(OxidesImageFinal)) max(max(OxidesImageFinal))]);
set([ax1, ax2],'Position',[.17 .11 .685 .815]);

% t = 23 h
f10 = figure(10);
ax1 = axes;
imagesc(CropGrains23);
view(2)
ax2 = axes;
imagesc(CropOxidesImage23, 'AlphaData',~isnan(CropOxidesImage23));
linkaxes([ax1,ax2]);
ax1.Visible = 'off';
ax2.Visible = 'off';
ax2.XTick = [];
ax2.YTick = [];
colormap(ax1, CM1);
colormap(ax2, flipud(viridis));
clim([min(min(OxidesImageFinal)) max(max(OxidesImageFinal))]);
set([ax1, ax2],'Position',[.17 .11 .685 .815]);

% t = 33 h
f11 = figure(11);
ax1 = axes;
imagesc(CropGrains33);
view(2)
ax2 = axes;
imagesc(CropOxidesImage33, 'AlphaData',~isnan(CropOxidesImage33));
linkaxes([ax1,ax2]);
ax1.Visible = 'off';
ax2.Visible = 'off';
ax2.XTick = [];
ax2.YTick = [];
colormap(ax1, CM1);
colormap(ax2, flipud(viridis));
clim([min(min(OxidesImageFinal)) max(max(OxidesImageFinal))]);
set([ax1, ax2],'Position',[.17 .11 .685 .815]);

% t = 43 h
f12 = figure(12);
ax1 = axes;
imagesc(CropGrains43);
view(2)
ax2 = axes;
imagesc(CropOxidesImage43, 'AlphaData',~isnan(CropOxidesImage43));
linkaxes([ax1,ax2]);
ax1.Visible = 'off';
ax2.Visible = 'off';
ax2.XTick = [];
ax2.YTick = [];
colormap(ax1, CM1);
colormap(ax2, flipud(viridis));
clim([min(min(OxidesImageFinal)) max(max(OxidesImageFinal))]);
set([ax1, ax2],'Position',[.17 .11 .685 .815]);

% Save images/figures
cd('C:\Users\ecfadely\OneDrive - University of California, Davis\Manuscripts\Mn Biomineralization Protocol\Figures\Figure_6');

imwrite(CropT19, 'CroppedT19Image.tif');
imwrite(CropT20, 'CroppedT20Image.tif');
imwrite(CropT21, 'CroppedT21Image.tif');
imwrite(CropT23, 'CroppedT23Image.tif');
imwrite(CropT33, 'CroppedT33Image.tif');
imwrite(CropT43, 'CroppedT43Image.tif');

exportgraphics(f7, 'Cropped_Subtracted_Matrix_19.tif', 'Resolution', 1500);
exportgraphics(f8, 'Cropped_Subtracted_Matrix_20.tif', 'Resolution', 1500);
exportgraphics(f9, 'Cropped_Subtracted_Matrix_21.tif', 'Resolution', 1500);
exportgraphics(f10, 'Cropped_Subtracted_Matrix_23.tif', 'Resolution', 1500);
exportgraphics(f11, 'Cropped_Subtracted_Matrix_33.tif', 'Resolution', 1500);
exportgraphics(f12, 'Cropped_Subtracted_Matrix_43.tif', 'Resolution', 1500);

%% Figure 7: Rate of Mn oxide precipitation as a function of time

close all;

cd('C:\Users\ecfadely\OneDrive - University of California, Davis\Manuscripts\Mn Biomineralization Protocol\Data');

% Load ICP-MS reactor digestate measurements
ICPMS_Digest = readtable('ICP_MS_Reactor_Digestate_Triplicates.csv');
MassVals = table2array(ICPMS_Digest(:,2));

% Rescale pixel intensity data by measured ICP-MS reactor digestate
GM1_Rescale = (GrayMass1 ./ GrayMass1(1,51)) * MassVals(1,1);
GM2_Rescale = (GrayMass2 ./ GrayMass2(1,51)) * MassVals(2,1);
GM3_Rescale = (GrayMass3 ./ GrayMass3(1,51)) * MassVals(3,1);

% Concatenate rescaled data
Gray_Mat = cat(1, GM1_Rescale, GM2_Rescale, GM3_Rescale);

% Take derivative of mass 
Gray_Deriv = cat(1, diff(Gray_Mat(1,:))/0.5, diff(Gray_Mat(2,:))/0.5, diff(Gray_Mat(3,:))/0.5);

% Calculate average
AvgGray_DT = mean(Gray_Deriv);

% Calculate standard deviation
STDGray_DT = std(Gray_Deriv);

% Data for patch plot
curve1Gray_DT = AvgGray_DT - STDGray_DT;
curve2Gray_DT = AvgGray_DT + STDGray_DT;

% Plotting parameters
tdiff = 4:0.5:28.5;
GRC = [0.0667    0.2118    0.5020];

% Plot data
f1 = figure(1);
set(gcf,'Position', [100 100 525 300]); 
box on;
hold on;
patch([tdiff, flip(tdiff)], [curve2Gray_DT, flip(curve1Gray_DT)], GRC, 'FaceAlpha', 0.2, 'EdgeColor','none');
plot(tdiff, AvgGray_DT, '-o', 'MarkerSize', 4, 'MarkerFaceColor', GRC, 'Color', GRC, 'LineWidth', 1.5, 'DisplayName', 'Composite');
hold off;
xticks([0 5 10 15 20 25 30]); xticklabels({'0' '20' '25' '30' '35' '40' '45'});
xlabel('Time (h)');
ylabel('Precipitation rate (μg h^{-1})');
set(gca, 'FontSize', 14, 'LineWidth', 1.25);

% Save figure
cd('C:\Users\ecfadely\OneDrive - University of California, Davis\Manuscripts\Mn Biomineralization Protocol\Figures\Figure_7');
saveas(f1, 'PrecipitationRate.svg');
