% function [ECens, ECiter] = ECens( times, celP, IRGS)

%	ENENS is a wrapper function for ECFlux.m to calculate eddy covariance from ensemble arrays
%   Jason Kelley
%	Last modified 21Apr2018
%	Questions?  <a href="mailto:jasonrk@uidaho.edu">Contact</a>

%   INPUTS : IRGS(n) IRGASON data from n sensors, with four fields: .wind, .Ts, .rhov, .Thmp

cd('C:\Users\kelleyja\Data\eclipse\matfiles')
% 20AUG starts at 8
% 21AUG starts at 32
% 22AUG starts at 56
clear DATA
start = 2; fin = 66;
fn = ls('TS*'); fn = fn(:,3:end);

f = start;
fns = ['*',fn(f,:)]; fnn = ls(fns); vn = fnn(:,1:2);
for g = 1:size(fnn,1)
    %DATA.(vn(g,:)) = load(fnn(g,:),vn(g,:));
    load(fnn(g,:));
end
vnn = {'TS';'A1';'A2';'B1';'B2';'C1';'C2';'D1';'D2';'D3';'P';'fieldNames'};
DATA = v2struct(vnn);
clear(vnn{1:end-1})
for f = start+1:fin
    fns = ['*',fn(f,:)]; fnn = ls(fns);
    for g = 1:size(fnn,1)
        load(fnn(g,:))
    end
    DATAt = v2struct(vnn);
    DATA.TS = [DATA.TS; DATAt.TS];
    DATA.P  = [DATA.P;  DATAt.P ];
    for h = 2:size(vnn)-2
        DATA.(vnn{h,:}) = catStruct(DATA.(vnn{h,:}),DATAt.(vnn{h,:}));
    end
    
    clear DATAt
    clear(vnn{1:end-1})
clear g h
end

for h = 2:length(vnn)-2
    DATA.(vnn{h,:}).wind = [DATA.(vnn{h,:}).Ux DATA.(vnn{h,:}).Uy DATA.(vnn{h,:}).Uz];
    DATA.(vnn{h,:}) = rmfield(DATA.(vnn{h,:}),'Ux'); 
    DATA.(vnn{h,:}) = rmfield(DATA.(vnn{h,:}),'Uy');
    DATA.(vnn{h,:}) = rmfield(DATA.(vnn{h,:}),'Uz');
    
end
DATA.P = (DATA.P+800)./10; % add 800 hPa sensor offset, convert to kPa

clear f fn start fin vnn vn fns fnn h

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% iterative flux calculation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
cd('C:\Users\kelleyja\Data\eclipse');

clearvars -except ECiter* DATA H*
clear ECiter
fns = fieldnames(DATA);

% set time parameters
TIME.tavEC = 20;                                % EC averaging time, minutes
TIME.t_smp = 0.05;                              % sample time spacing (s)
TIME.freq  = 1/TIME.t_smp;                      % sampling frequency (Hz)
TIME.n     = length(DATA.TS);
TIME.ttot  = TIME.n*TIME.t_smp/60;              % length of record (minutes)
TIME.start = datevec(DATA.TS(1));
TIME.fin   = datevec(DATA.TS(end));
TIME.npEC   = floor(TIME.ttot/TIME.tavEC);      % number of averaging periods
TIME.nEC   = floor(60/TIME.t_smp*TIME.tavEC);   % samples per period

for j = 1:TIME.npEC
    jj = (j-1)*TIME.nEC+1:j*TIME.nEC;
    for i = 1:length(fns)-2
        [EC(i,j), WIND(i,j),T_p] = EC_Hflux( DATA.TS(jj), DATA.(fns{i+1}).wind(jj,:),...
             DATA.(fns{i+1}).Ts(jj), DATA.(fns{i+1}).e(jj), DATA.(fns{i+1}).Th(jj), DATA.P(jj));
    end
end

clear i j jj T_p

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ensemble flux calculation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

reps = size(EC,1);
TIME.Nens  = length(DATA.TS)*reps;            % number of samples / ensemble

% collate wind vectors and decompose fluctuations

if exist('WINDens','var')
    clear WINDens;
    WINDens(TIME.npEC).w = [];
else
    WINDens(TIME.npEC).w = [];
end

% collate rotated wind speeds for each averaging period
for i = 1:TIME.npEC
    WINDens(i).w = WIND(1,i).rot(:,3);
    WINDens(i).v = WIND(1,i).rot(:,2);
    WINDens(i).u = WIND(1,i).rot(:,1);
    %collate all sensors
    for r = 2:reps
        WINDens(i).w = [WINDens(i).w; WIND(r,i).rot(:,3)];
        WINDens(i).v = [WINDens(i).v; WIND(r,i).rot(:,2)];
        WINDens(i).u = [WINDens(i).u; WIND(r,i).rot(:,1)];
    end
    % calculate means and decompose flucuations
    WINDens(i).w_m = nanmean(WINDens(i).w);
    WINDens(i).w_p = WINDens(i).w-WINDens(i).w_m;
    WINDens(i).v_m = nanmean(WINDens(i).v);
    WINDens(i).v_p = WINDens(i).v-WINDens(i).v_m;
    WINDens(i).u_m = nanmean(WINDens(i).u);
    WINDens(i).u_p = WINDens(i).u-WINDens(i).u_m;

end
clear i r

% collate SUSAN data

if exist('ECens','var')
    clear ECens
    ECens(TIME.npEC).Ts = [];
else
    ECens(TIME.npEC).Ts = [];
end

R = 8.31451;                           % [J /mol /K] Gas constant
for i = 1:TIME.npEC
    ii = (i-1)*TIME.nEC+1:i*TIME.nEC;
    ECens(i).Ts = DATA.(fns{2}).Ts(ii);
    ECens(i).e  = DATA.(fns{2}).e(ii);
    ECens(i).Th = DATA.(fns{2}).Th(ii);
    for r = 2:reps
        ECens(i).Ts = [ECens(i).Ts; DATA.(fns{r+1}).Ts(ii)];
        ECens(i).e  = [ECens(i).e;  DATA.(fns{r+1}).e(ii) ];
        ECens(i).Th = [ECens(i).Th; DATA.(fns{r+1}).Th(ii)];        
    end
% Thermodynamic calcs
    ECens(i).Tm   = nanmean(ECens(i).Th);
    ECens(i).TK   = ECens(i).Tm + 273.15;
    ECens(i).P    = repmat(DATA.P(ii),reps,1);
    rhov = ECens(i).e.*18.01./(R.*(ECens(i).Th+273.15));       % [kg/m^3] vapor density
    rhod = (ECens(i).P-ECens(i).e)./0.287./ECens(i).TK; % [kg/m^3] dry air dens.
    chiv = rhov./rhod;
    rcp = rho_cp(ECens(i).Tm, nanmean([ECens(i).P]), nanmean(rhov));
    
% Temperature flucuations
    Ts_m = nanmean(ECens(i).Ts);         % sonic temperature
    Ts_p = ECens(i).Ts-Ts_m;
    T   = ECens(i).Ts./(1+0.32.*chiv);   % true temperature corrected for density of water vapor
    T_m = nanmean(T);                           % Leuning (2004) provides this eq.as approx., referring to
    T_p = T-T_m;                                % Kaimal & Gaynor (1991) which states Ts = T*[1+0.32(e/p)]

% Covariances
    ECens(i).uw  = nanmean(WINDens(i).u_p.*WINDens(i).w_p);     % momentum flux
    ECens(i).vw  = nanmean(WINDens(i).v_p.*WINDens(i).w_p);     % 
    ECens(i).wT  = nanmean(WINDens(i).w_p.*T_p);                % ~Sensible flux
    ECens(i).wTs = nanmean(WINDens(i).w_p.*Ts_p);               % Buoyancy flux
    
    ECens(i).ustar = (ECens(i).uw^2+ECens(i).vw^2)^0.25;        % Friction velocity
    ECens(i).ObuL  = -1*(ECens(i).ustar^3)/(0.4*9.81/(Ts_m+273.15)*ECens(i).wTs);
                                                                % Obukhov Length
% Flux calcs                                                                                
    ECens(i).H  = rcp .* ECens(i).wT;     % Sensible heat flux [W/m^2]
    ECens(i).Hb  = rcp .* ECens(i).wTs;   % Buoyant heat flux [W/m^2]
    
end
clear i ii r R rhov rhod Ts_m Ts_p T T_m T_p chiv rcp reps

%% save output structures
ECens = rmfield(ECens,{'Ts';'e';'Th';'P'});
%savestr = ['EC_', num2str(TIME.tavEC),'min'];
savestr = 'EC_30s';
save(savestr,'EC','ECens','TIME')

%% load EC structures
cd('C:\Users\kelleyja\Data\eclipse\matfiles')
EC_20m = load('EC_20min.mat');
% EC_10m = load('EC_10min.mat');
% EC_5m  = load('EC_05min.mat');
% EC_2m  = load('EC_02min.mat');
EC_1m  = load('EC_01min.mat');
% EC_30s = load('EC_30s.mat');

Hb_20m = reshape([EC_20m.EC(:,:).Hb],9,EC_20m.TIME.npEC);
Hb_ave = nanmean(Hb_20m,1);
%%  EC plots
figure, plot([EC_20m.EC(1,:).ctr],Hb_ave,'k','Linewidth',2)
hold('on')
xlr = [[EC_20m.EC(1,:).ctr] fliplr([EC_20m.EC(1,:).ctr])];
yhl = [nanmax(Hb_20m,[],1) fliplr(nanmin(Hb_20m,[],1))];
patch(xlr,yhl,'r','Facealpha',0.25,'EdgeColor','none')
hhmmTick(240,gca); grid('on');

clear xlr yhl

%% Rnet
curdir = cd;
cd('C:\Users\kelleyja\Data\ANN-ECL')
EB = load('ANN-ECL_EB_DATA.mat'); EB = EB.EB_DATA;
ANNEC = load('20mEC_19Aug-31Aug.mat');
cd(curdir); clear curdir

EB.Rnet = EB.SRUp_Avg+EB.IRUp_Avg - (EB.SRDn_Avg+EB.IRDn_Avg);
Rnet = interp1(EB.TIMESTAMP,EB.Rnet,[SR_1m.SR(1,:).ctr]);
G = interp1(EB.TIMESTAMP,EB.GrFlux_Avg,[SR_1m.SR(1,:).ctr]);
plot([SR_1m.SR(1,:).ctr],Rnet-G,'Color',[0.9 0.5 0],'Linewidth',1)

Rnet_p =  Rnet - convsame(Rnet,20);
G_p    =  G - convsame(G,20);
%% SR, EC mw_std
SR_std = mw_std(Hsr_1a',20)';
EC_std = mw_std([EC_1m.ECens.Hb]',20)';
Hsr_p = Hsr_1a - convsame(Hsr_1a,20);
Hec_p = [EC_1m.ECens.Hb] - convsame([EC_1m.ECens.Hb],20);
figure, crosscorr(naninterp(Hec_p),naninterp(Hsr_p),1440)
title('Cross Correlation Function [H_{EC}'' x H_{SR}''] - 1 minute flucuations)')
%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Cross Correlation plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%      H x Rnet
dayi = Rnet > 50;
trnsi = Rnet < 50 & Rnet > 0;
nighti = Rnet <= 5;

figure, crosscorr(naninterp(Hec_p(dayi)),naninterp(Rnet_p(dayi)),360)
title('Cross Correlation [H_{EC}'' x R_{net}''] - 1 minute fluctuations (daytime)')
figure, crosscorr(naninterp(Hsr_p(dayi)),naninterp(Rnet_p(dayi)),360)
title('Cross Correlation [H_{SR}'' x R_{net}''] - 1 minute fluctuations(daytime)')

figure, crosscorr(naninterp(Hec_p(trnsi)),naninterp(Rnet_p(trnsi)),60)
title('Cross Correlation [H_{EC}'' x R_{net}''] - 1 minute fluctuations (dawn/dusk)')
figure, crosscorr(naninterp(Hsr_p(trnsi)),naninterp(Rnet_p(trnsi)),60)
title('Cross Correlation [H_{SR}'' x R_{net}''] - 1 minute fluctuations (dawn/dusk)')

figure, crosscorr(naninterp(Hec_p(nighti)),naninterp(Rnet_p(nighti)),360)
title('Cross Correlation [H_{EC}'' x R_{net}''] - 1 minute fluctuations (night)')
figure, crosscorr(naninterp(Hsr_p(nighti)),naninterp(Rnet_p(nighti)),360)
title('Cross Correlation [H_{SR}'' x R_{net}''] - 1 minute fluctuations (night)')
%%      H x G
figure, crosscorr(naninterp(Hec_p(dayi)),naninterp(G_p(dayi)),360)
title('Cross Correlation [H_{EC}'' x G''] - 1 minute fluctuations (daytime)')
figure, crosscorr(naninterp(Hsr_p(dayi)),naninterp(G_p(dayi)),360)
title('Cross Correlation [H_{SR}'' x G''] - 1 minute fluctuations(daytime)')

figure, crosscorr(naninterp(Hec_p(trnsi)),naninterp(G_p(trnsi)),60)
title('Cross Correlation [H_{EC}'' x G''] - 1 minute fluctuations (dawn/dusk)')
figure, crosscorr(naninterp(Hsr_p(trnsi)),naninterp(G_p(trnsi)),60)
title('Cross Correlation [H_{SR}'' x G''] - 1 minute fluctuations (dawn/dusk)')

figure, crosscorr(naninterp(Hec_p(nighti)),naninterp(G_p(nighti)),360)
title('Cross Correlation [H_{EC}'' x G''] - 1 minute fluctuations (night)')
figure, crosscorr(naninterp(Hsr_p(nighti)),naninterp(G_p(nighti)),360)
title('Cross Correlation [H_{SR}'' x G''] - 1 minute fluctuations (night)')
%%      Rnet x G
figure, crosscorr(naninterp(Rnet_p(dayi)),naninterp(G_p(dayi)),360)
title('Cross Correlation [R_{net}'' x G''] - 1 minute fluctuations (daytime)')

figure, crosscorr(naninterp(Rnet_p(trnsi)),naninterp(G_p(trnsi)),60)
title('Cross Correlation [R_{net}'' x G''] - 1 minute fluctuations (dawn/dusk)')

figure, crosscorr(naninterp(Rnet_p(nighti)),naninterp(G_p(nighti)),360)
title('Cross Correlation [R_{net}'' x G''] - 1 minute fluctuations (night)')
%%      LW x G and SW x G
LW = EB.IRUp_Avg - EB.IRDn_Avg;
LW = interp1(EB.TIMESTAMP,LW,[SR_1m.SR(1,:).ctr]);
LW_p = LW - convsame(LW,20);
SW = EB.SRUp_Avg - EB.SRDn_Avg;
SW = interp1(EB.TIMESTAMP,SW,[SR_1m.SR(1,:).ctr]);
SW_p = SW - convsame(SW,20);
%%      SW' x G' & LW' x G'
figure, subplot(2,2,1)
crosscorr(naninterp(LW_p(dayi)),naninterp(G_p(dayi)),360)
title('Cross Correlation [LW'' x G''] - 1 minute fluctuations (day)')
xlim([-60 60]), ylim([-0.45 0.45]), zeroline(gcf,'y')
xlabel('')
ax = gca; ax.Position = [0.075 0.55 0.4 0.375];

subplot(2,2,2), crosscorr(naninterp(LW_p(nighti)),naninterp(G_p(nighti)),360)
title('Cross Correlation [LW'' x G''] - 1 minute fluctuations (night)')
xlim([-60 60]), ylim([-0.45 0.45]), zeroline(gcf,'y')
xlabel(''), ylabel('')
ax = gca; ax.Position = [0.525 0.55 0.4 0.375];

subplot(2,2,4), crosscorr(naninterp(SW_p(nighti)),naninterp(G_p(nighti)),360)
title('Cross Correlation [SW'' x G''] - 1 minute fluctuations (night)')
xlim([-60 60]), ylim([-0.45 0.45]), zeroline(gcf,'y')
ylabel(''), xlabel('Lag (minutes)')
ax = gca; ax.Position = [0.525 0.075 0.4 0.375];

subplot(2,2,3), crosscorr(naninterp(SW_p(dayi)),naninterp(G_p(dayi)),360)
title('Cross Correlation [SW'' x G''] - 1 minute fluctuations (day)')
xlim([-60 60]), ylim([-0.45 0.45]), zeroline(gcf,'y')
xlabel('Lag (minutes)')
ax = gca; ax.Position = [0.075 0.075 0.4 0.375];

clear ax
%%      LW x Hsr & SW x Hsr
figure, subplot(2,2,1)
crosscorr(naninterp(LW_p(dayi)),naninterp(Hsr_p(dayi)),360)
title('Cross Correlation [LW'' x H_{SR}''] - 1 minute fluctuations (day)')
xlim([-60 60]), ylim([-0.2 0.3]), zeroline(gcf,'y'); yticks(-0.1:0.1:0.3);
xlabel('')
ax = gca; ax.Position = [0.075 0.55 0.4 0.375];

subplot(2,2,2), crosscorr(naninterp(LW_p(nighti)),naninterp(Hsr_p(nighti)),360)
title('Cross Correlation [LW'' x H_{SR}''] - 1 minute fluctuations (night)')
xlim([-60 60]), ylim([-0.2 0.3]), zeroline(gcf,'y'); yticks(-0.1:0.1:0.3);
xlabel(''), ylabel('')
ax = gca; ax.Position = [0.525 0.55 0.4 0.375];

subplot(2,2,4), crosscorr(naninterp(SW_p(nighti)),naninterp(Hsr_p(nighti)),360)
title('Cross Correlation [SW'' x H_{SR}''] - 1 minute fluctuations (night)')
xlim([-60 60]), ylim([-0.2 0.3]), zeroline(gcf,'y'); yticks(-0.1:0.1:0.3);
ylabel(''), xlabel('Lag (minutes)')
ax = gca; ax.Position = [0.525 0.075 0.4 0.375];

subplot(2,2,3), crosscorr(naninterp(SW_p(dayi)),naninterp(Hsr_p(dayi)),360)
title('Cross Correlation [SW'' x H_{SR}''] - 1 minute fluctuations (day)')
xlim([-60 60]), ylim([-0.2 0.3]), zeroline(gcf,'y'); yticks(-0.1:0.1:0.3);
xlabel('Lag (minutes)')
ax = gca; ax.Position = [0.075 0.075 0.4 0.375];

clear ax
%%      LW x Hec & SW x Hec
figure, subplot(2,2,1)
crosscorr(naninterp(LW_p(dayi)),naninterp(Hec_p(dayi)),360)
title('Cross Correlation [LW'' x H_{EC}''] - 1 minute fluctuations (day)')
xlim([-60 60]), ylim([-0.2 0.3]), zeroline(gcf,'y'); yticks(-0.1:0.1:0.3);
xlabel('')
ax = gca; ax.Position = [0.075 0.55 0.4 0.375]

subplot(2,2,2), crosscorr(naninterp(LW_p(nighti)),naninterp(Hec_p(nighti)),360)
title('Cross Correlation [LW'' x H_{EC}''] - 1 minute fluctuations (night)')
xlim([-60 60]), ylim([-0.2 0.3]), zeroline(gcf,'y'); yticks(-0.1:0.1:0.3);
xlabel(''), ylabel('')
ax = gca; ax.Position = [0.525 0.55 0.4 0.375];

subplot(2,2,4), crosscorr(naninterp(SW_p(nighti)),naninterp(Hec_p(nighti)),360)
title('Cross Correlation [SW'' x H_{EC}''] - 1 minute fluctuations (night)')
xlim([-60 60]), ylim([-0.2 0.3]), zeroline(gcf,'y'); yticks(-0.1:0.1:0.3);
ylabel(''), xlabel('Lag (minutes)')
ax = gca; ax.Position = [0.525 0.075 0.4 0.375];

subplot(2,2,3), crosscorr(naninterp(SW_p(dayi)),naninterp(Hec_p(dayi)),360)
title('Cross Correlation [SW'' x H_{EC}''] - 1 minute fluctuations (day)')
xlim([-60 60]), ylim([-0.2 0.3]), zeroline(gcf,'y'); yticks(-0.1:0.1:0.3);
xlabel('Lag (minutes)')
ax = gca; ax.Position = [0.075 0.075 0.4 0.375];

clear ax

%%      (Rnet'-G') x Hec & (Rnet'-G') x Hsr
figure, subplot(2,2,1)
crosscorr(naninterp(Rnet_p(dayi)-G_p(dayi)),naninterp(Hec_p(dayi)),360)
title('Cross Correlation [(R_n''-G'') x H_{EC}''] - 1 minute fluctuations (day)')
xlim([-60 60]), ylim([-0.2 0.3]), zeroline(gcf,'y'); yticks(-0.1:0.1:0.3);
xlabel('')
ax = gca; ax.Position = [0.075 0.55 0.4 0.375]

subplot(2,2,2), crosscorr(naninterp(Rnet_p(nighti)-G_p(nighti)),naninterp(Hec_p(nighti)),360)
title('Cross Correlation [(R_n''-G'') x H_{EC}''] - 1 minute fluctuations (night)')
xlim([-60 60]), ylim([-0.2 0.3]), zeroline(gcf,'y'); yticks(-0.1:0.1:0.3);
xlabel(''), ylabel('')
ax = gca; ax.Position = [0.525 0.55 0.4 0.375];

subplot(2,2,4), crosscorr(naninterp(Rnet_p(nighti)-G_p(nighti)),naninterp(Hsr_p(nighti)),360)
title('Cross Correlation [(R_n''-G'') x H_{SR}''] - 1 minute fluctuations (night)')
xlim([-60 60]), ylim([-0.2 0.3]), zeroline(gcf,'y'); yticks(-0.1:0.1:0.3);
ylabel(''), xlabel('Lag (minutes)')
ax = gca; ax.Position = [0.525 0.075 0.4 0.375];

subplot(2,2,3), crosscorr(naninterp(Rnet_p(dayi)-G_p(dayi)),naninterp(Hsr_p(dayi)),360)
title('Cross Correlation [(R_n''-G'') x H_{SR}''] - 1 minute fluctuations (day)')
xlim([-60 60]), ylim([-0.2 0.3]), zeroline(gcf,'y'); yticks(-0.1:0.1:0.3);
xlabel('Lag (minutes)')
ax = gca; ax.Position = [0.075 0.075 0.4 0.375];

clear ax
%% Time series of cross correlation
clear RnG_Hsr % Hsr_SW

nnR = naninterp(Rnet_p-G_p);
nnH = naninterp(Hsr_p);

RnG_Hsr(SR_1m.TIME.npSR) = struct();
RnG_Hsr = deal2SA(nan(SR_1m.TIME.npSR,1),RnG_Hsr,'XCmax'); 
RnG_Hsr = deal2SA(nan(SR_1m.TIME.npSR,1),RnG_Hsr,'XCmin');
RnG_Hsr = deal2SA(nan(SR_1m.TIME.npSR,1),RnG_Hsr,'maxlag');
RnG_Hsr = deal2SA(nan(SR_1m.TIME.npSR,1),RnG_Hsr,'minlag');

for i = 6:SR_1m.TIME.npSR-6
    ii = i-5:i+5;
    [XCF,lags,bounds] = crosscorr(nnR(ii),nnH(ii),10);
    [RnG_Hsr(i).XCmax,j] = max(XCF);
    RnG_Hsr(i).maxlag = lags(j);
    [RnG_Hsr(i).XCmin,j] = min(XCF);
    RnG_Hsr(i).minlag = lags(j);
end

clear i ii j nnR nnH XCF lags bounds
