clear all;
close all;

MASTERPATH = '/Users/eliz/DATA/PISTON_2019/Sally_Ride/ship/';
mfiles = dir([MASTERPATH 'MET_SallyRide_PISTON_2019_thru_*.mat']);
nfiles = length(mfiles);

for i = 1:nfiles
    
    file = [MASTERPATH char(mfiles(i).name)];
    disp(['Loading file # ' strnd(i,0) '/' strnd(nfiles,0)  ': ' file]);
    load(file);
    
    if i ==1
%         disp('first file');
        f = fields(data);
            %%% load all variables
            for j = 1:length(f)
                %%% grab field from small array and put in big array
                MET.(f{j}) = data.(f{j});
            end
        
    else
%             disp('rest of variables');
            %%% grab and put rest of variables in big array
            for j = 1:length(f)
                MET.(f{j}) = [MET.(f{j}); data.(f{j})];
            end
    
    end
    
end % end file loop

% %%% problem with TSG
bad2 = find(MET.USW_flow_meter_2 > 3.5 | MET.USW_flow_meter_2 < 1.75); 
bad1 = find(MET.USW_flow_meter_1 > 3.5 | MET.USW_flow_meter_1 < 1.75); 

%%% 1 = good. 0 = bad
TSG_flow1_ID = ones(length(MET.dnum),1);
TSG_flow2_ID = ones(length(MET.dnum),1);
TSG_flow1_ID(bad1) = 0;
TSG_flow2_ID(bad2) = 0;
MET.TSG_goodflow1 = TSG_flow1_ID;
MET.TSG_goodflow2 = TSG_flow2_ID;

%%% adjust data
MET.TSG_conductivity_1(bad1) = nan;
MET.TSG_conductivity_2(bad1) = nan;
MET.TSG_temperature_1(bad1) = nan;
MET.TSG_temperature_2(bad1) = nan;
MET.salinity_1(bad1) = nan;
MET.salinity_2(bad1) = nan;
MET.sigma_t_1(bad1) = nan;
MET.sigma_t_2(bad1) = nan;

% calculate ship precipitation rate
% figure;hold on;
bad = find(diff(MET.precipitation) > 0.2 | diff(MET.precipitation) < -0.2);
% plot(MET.dnum, MET.precipitation,'.',MET.dnum(bad), MET.precipitation(bad),'.')
paccum = MET.precipitation;
paccum(bad) = nan;


%%% The original USPS data are every 5 sec but it needs to be standardized
%%% to a continuous level.

[yyy1, mmm1, ddd1, hhh1, MMM1, sss1] = datevec(MET.dnum(1));
[yyy2, mmm2, ddd2, hhh2, MMM2, sss2] = datevec(MET.dnum(end));

ndays = ddd2 - ddd1 +1;

newtime  = nan(ndays*1440,1);                        

% make a 1 min time series over same time period;
p = 1;
for d = 0:ndays % days
    for j = 0:23 % hrs
        for k = 0:59 % mins
%             for l = 0:5:55 % secs
            newtime(p) = datenum(2019,9,ddd1+d,j,k,0);
            p = p+1;
%             end
        end
    end
end
np = p-1;

[u_dnum, ID_u] = unique(MET.dnum);
paccum1 = interp1(MET.dnum(ID_u), paccum(ID_u), newtime);
x = [0;diff(paccum1)];
y = find(x < -0.2);
% plot(newtime(y), paccum1(y),'o');

paccum2 = paccum1;
paccum2(y) = nan;

MET.paccum1min = paccum2;
MET.dnum1min = newtime;

pr = cum2rate(MET.paccum1min', MET.dnum1min');

%%% this might be the best version, actually. 

%%% the effective resolution of paccum2 is really low... +/- 5 mm/hr. So,
%   use paccum2 instead to filter out bad data in MET.precipitation. 
% 
% ... that might not actually change. Because the effective resolution of the
%     actual measurement was only 0.2 mm. Look at specifications for
%     RMyoung model... what does documentation say? 

%%%%% find times when MET.precipitation is excessively higher or lower than
%%%%% the "low pass filtered" good data. Now that we have a better version
%%%%% of that... 

% x_id = nan(length(x),1);
% for i = 5:length(x)-5
%     if x(i) < -0.2 && x(i-2) > 0.2 && x(i-2) > 0.2
%         x_id(i) = 1;
%     end
%     
% end

% plot(newtime(abs(x) > 0.2), paccum1(abs(x) > 0.2),'x');

% %bad ones: 
% outlier = find(abs(pcd) > 0.6);
% pcc = pc;
% pcc(outlier) = nan;
% 
% 
% % MET.precipitation(MET.precipitation < 0.1) = nan;
% % good = find(diff(MET.precipitation) > 0.1);
% % y = diff(MET.precipitation);
% % bad = find(diff(MET.precipitation) > 0.2);
% pc = MET.precipitation;
% pdiff = [0; diff(pc)];
% bad2 = find(abs(pdiff) > 0.1);
% 
% % pc(bad) = nan;
% 
% %bad ones: 
% outlier = find(abs(pcd) > 0.6);
% pcc = pc;
% pcc(outlier) = nan;
% 
% % 
% % flag = nan(length(MET.dnum),1);
% % pcd = [0; diff(pc)];
% % for i = 3:length(MET.dnum)-3
% %     if pcd(i) < -1 & pcd(i-1) > -0.3 & pcd(i+1) > -0.3
% %         flag(i) = 1;
% %     end
% % end
% % 
% % wh_flag = find(flag ==1);
% % 
% % bad = find(diff(MET.precipitation) > 0.2);
% % 
% % % pacc_filt = f15(MET.precipitation);
% % 
% % mask = nan(length(MET.dnum),1);
% % 
% % % remove outliers
% % for i = 3:length(MET.dnum)-3
% %     if y(i) > 0 
% % end
% % 
% % pacc_filt_good = find(diff(pacc_filt) > 0.01);
% % pacc_filt_bad = find(diff(pacc_filt) <= 0.01);
% % pacc_filt2 = pacc_filt;
% % pacc_filt2(pacc_filt_bad) = nan;
% % 
% % figure;hold on;
% % plot(MET.dnum, pacc_filt,'.');
% % plot(MET.dnum, pacc_filt2,'o');
% % 
% % 
% % 
% % % MET.rain_rate = cum2rate(MET.precipitation',MET.dnum');

% MET.pcum = pcc;


%%% save the data
save([MASTERPATH 'MET_SallyRide_PISTON_2019'], 'MET');



% figure
% yyaxis left; plot(MET.dnum, MET.TSG_temperature_1);
% ylabel('T, deg C');
% 
% yyaxis right; plot(MET.dnum, MET.USW_flow_meter_1);
% ylabel('flow, L/min');
% 
% legend('flow','T')
% grid on;
% datetick('x','keeplimits');

% hold on;
% flow1 = MET.USW_flow_meter_1;
% flow1(bad1) = nan;
% plot(MET.dnum, flow1,'ok')