% Script to read in the Station M trap data and then compare it to the SES
% mean attenuance flux

%% Set up the data directory and load the files
datadir = '/Users/mestapa/Documents/Research_Maine/NSF-STTR/MBARI_SES_data/';
load([datadir 'SESatnsumm.mat']);
load([datadir 'SESdata.mat'],'bgr');
trapfnam = 'P64_Estapa_sourcefileMAIN_MergedDaily50_600.csv';
rawtrapdata = readmatrix(trapfnam,'NumHeaderLines',1,'OutputType','string');

%% Parse the trap data file
startcollect = datenum(rawtrapdata(:,1),'mm/dd/yyyy');
endcollect = datenum(rawtrapdata(:,2),'mm/dd/yyyy');
massflux_mgm2d = str2double(rawtrapdata(:,8));
[~,uniquerowinds] = unique(massflux_mgm2d);
uniquerowinds = sort(uniquerowinds);
startcollect = startcollect(uniquerowinds);
endcollect = endcollect(uniquerowinds);
massflux_mgm2d = massflux_mgm2d(uniquerowinds);
OCflux_mgm2d = str2double(rawtrapdata(uniquerowinds,9));
TCflux_mgm2d = str2double(rawtrapdata(uniquerowinds,10));
ICflux_mgm2d = str2double(rawtrapdata(uniquerowinds,11));

%% Convert SES attenuance to attenuance flux
SEStraparea = 0.5;
% SESimgarea = length(find(isfinite(bgr(:))))./(0.07194.^2)./(1000000.^2);  % pix / (pix/um^2) / (um^2/m^2)
SESimgarea = pi.*(25./1000./2).^2;  % diameter of lower end of funnel is 25 mm, compute area in m^2
SEScolltime = 2/24;
meanatnflux = meanatn./SEScolltime.*SESimgarea./SEStraparea;

%% Compare the two datasets
startwindow = datenum([2014 10 26]); % start of first trap collection period
endwindow = datenum([2015 05 14]);   % end of last trap collection period
goodtrapinds = intersect(find(startcollect>=startwindow),find(endcollect<=endwindow));
goodSESinds = intersect(find(imgdates>=startwindow),find(imgdates<=endwindow));

atnflux_trap = NaN(size(OCflux_mgm2d(goodtrapinds),1),2);
for i = 1:length(goodtrapinds)
    these = intersect(find(imgdates>=startcollect(goodtrapinds(i))),...
        find(imgdates<endcollect(goodtrapinds(i))));
    atnflux_trap(i,1) = mean(meanatnflux(these),'omitnan');
    atnflux_trap(i,2) = std(meanatnflux(these),'omitnan');
    meantime_trap(i,1) = (startcollect(goodtrapinds(i))+endcollect(goodtrapinds(i)))./2;
    OCflux_comp_mgm2d(i,1) = OCflux_mgm2d(goodtrapinds(i));
    massflux_comp_mgm2d(i,1) = massflux_mgm2d(goodtrapinds(i));
end

%% Make a plot of OC flux
figure;scatter(atnflux_trap(:,1),OCflux_mgm2d(goodtrapinds),40,meantime_trap,'filled');
ch = colorbar;
set(ch,'TickLabels',datestr(get(ch,'Ticks')));
ylabel('POC flux, mg-C/m2/d');
xlabel('Atn flux, m2/m2/d');
% set(gca,'xscale','log','yscale','log')
% ylim([10^0,10^2]);xlim([10^-3 10^-2]);

%% Make a plot of mass flux
figure;scatter(atnflux_trap(:,1),massflux_mgm2d(goodtrapinds),40,meantime_trap,'filled');
ch = colorbar;
set(ch,'TickLabels',datestr(get(ch,'Ticks')));
ylabel('Mass flux, mg-C/m2/d');
xlabel('Atn flux, m2/m2/d');
% set(gca,'xscale','log','yscale','log')
% ylim([10^0,10^2]);xlim([10^-3 10^-2]);
%% Save it
save('SESatnvsPOC.mat','atnflux_trap','OCflux_comp_mgm2d','massflux_comp_mgm2d','meantime_trap');