% compare ID method b

%% load 10 Hz data for selected hour
clear
close all
Npts = 6000;
fsonic = 10;
dfl = '/Users/byronb/Raw_Data_Files.noindex/DYNAMO_2011/Revelle/flux/Processed_new/uvwStream_decorr_v3/UVWstream_v3_333_00_decorr.txt';
M = load(dfl,'-ascii');
U = M(:,2);
V = M(:,3);
W = M(:,4);
T = M(:,5);
Q = M(:,7);

Fxx = NaN(1,1+Npts/2);      % empty arrays for raw spectra
Puu = NaN(6,1+Npts/2);      % used for idiss calculations...
Pww = Puu; Ptt = Puu; Pqq = Puu;

Fx = NaN(1,48);               % empty arrays for 6x48-bin smoothed spectra
Pu = NaN(6,48);
Pw = Pu; Pt = Pu; Pq = Pu;

rwspd = NaN(6,1);

%% compute variance spectra and cospectra in 10-min windows
ii = 1;
for jj=1:6
    rwspd(jj) = mean(sqrt(U(ii:ii+Npts-1).^2 + V(ii:ii+Npts-1).^2 + W(ii:ii+Npts-1).^2));
    [Puu(jj,:),Fxx(1,:)] = psd2(detrend(U(ii:ii+Npts-1)),Npts,fsonic,hamming(Npts));
    Puu = real(Puu);
    [Pu(jj,:),Fx(1,:)] = specsmoo(Puu(jj,:),fsonic);
    [Pww(jj,:),~] = psd2(detrend(W(ii:ii+Npts-1))',Npts,fsonic,hamming(Npts));
    Pww = real(Pww);
    [Pw(jj,:),~] = specsmoo(Pww(jj,:),fsonic);
    [Ptt(jj,:),~] = psd2(detrend(T(ii:ii+Npts-1))',Npts,fsonic,hamming(Npts));
    Ptt = real(Ptt);
    [Pt(jj,:),~] = specsmoo(Ptt(jj,:),fsonic); % This is Tsonic...
    [Pqq(jj,:),~] = psd2(detrend(Q(ii:ii+Npts-1))',Npts,fsonic,hamming(Npts));
    Pqq = real(Pqq);
    [Pq(jj,:),~] = specsmoo(Pqq(jj,:),fsonic);
    ii = ii + Npts;
end;

%% METHOD B: compute Cx2 from spectral fit over inertial subrange in
% smoothed spectra using functionalized code

% Cu2 / Cw2
fminU = find(Fx>0.6*fsonic/20,1,'first');
if fminU<1
    fminU = 1;
end;
fmaxU = find(Fx<3.5*fsonic/20,1,'last');
if fmaxU>length(Fx)
    fmaxU = length(Fx);
end;
Cu2b = idiss_spec_fit(Pu, Fx, fminU, fmaxU, rwspd);
Cw2b = idiss_spec_fit(Pw, Fx, fminU, fmaxU, rwspd);

% Ct2
fminT = find(Fx>0.35*fsonic/20,1,'first');
if fminT<1
    fminT = 1;
end;
fmaxT = find(Fx<3*fsonic/20,1,'last');
if fmaxT>length(Fx)
    fmaxT = length(Fx);
end;
Ct2b = idiss_spec_fit(Pt, Fx, fminT, fmaxT, rwspd);

% Cq2
fminq = find(Fx>0.7*fsonic/20,1,'first');
if fminq<1
    fminq = 1;
end;
fmaxq = find(Fx<3.6*fsonic/20,1,'last');
if fmaxq>length(Fx)
    fmaxq = length(Fx);
end;
Cq2b = idiss_spec_fit(Pq, Fx, fminq, fmaxq, rwspd);

%% METHOD B: compute Cx2 from spectral fit over inertial subrange in
% smoothed spectra using original code

Cu2_b = NaN(6,1);
Cw2_b = Cu2_b;
Ct2_b = Cu2_b;
Cq2_b = Cu2_b;

fminU = min(find(Fx>.6*fsonic/20));
fmaxU = max(find(Fx<3.5*fsonic/20));

% U
cS=4;
Zpoly=polyfit(log(Fx(fminU:fmaxU)),log(cS*((2*pi/rwspd(1))^(2/3))*Pu(1,fminU:fmaxU).*(Fx(fminU:fmaxU).^(5/3))),4);
Zslope=Zpoly(4);
if Zslope<0.25,
   Cu2_b(1)=exp(Zpoly(5));
else
   plus=0;
   while 1
      zold=Zslope;
      Zslope=4*Zpoly(1).*(log(Fx(fminU+plus)).^3)+3*Zpoly(2).*(log(Fx(fminU+plus)).^2)+2*Zpoly(3).*log(Fx(fminU+plus))+Zpoly(4);
      plus=plus+1;
      if fminU+plus>=fmaxU,
         break
      end
      if Zslope>zold | abs(Zslope)<0.2,
         break,
      end
   end;
   Cu2_b(1)=exp(Zpoly(1).*(log(Fx(fminU+plus)).^4)+Zpoly(2).*(log(Fx(fminU+plus)).^3)+Zpoly(3).*(log(Fx(fminU+plus)).^2)+Zpoly(4).*log(Fx(fminU+plus))+Zpoly(5));
end;

Zpoly=polyfit(log(Fx(fminU:fmaxU)),log(cS*((2*pi/rwspd(2))^(2/3))*Pu(2,fminU:fmaxU).*(Fx(fminU:fmaxU).^(5/3))),4);
Zslope=Zpoly(4);
if Zslope<0.25,
   Cu2_b(2)=exp(Zpoly(5));
else
   plus=0;
   while 1
      zold=Zslope;
      Zslope=4*Zpoly(1).*(log(Fx(fminU+plus)).^3)+3*Zpoly(2).*(log(Fx(fminU+plus)).^2)+2*Zpoly(3).*log(Fx(fminU+plus))+Zpoly(4);
      plus=plus+1;
      if fminU+plus>=fmaxU,
         break
      end
      if Zslope>zold | abs(Zslope)<0.2,
         break,
      end
   end;
   Cu2_b(2)=exp(Zpoly(1).*(log(Fx(fminU+plus)).^4)+Zpoly(2).*(log(Fx(fminU+plus)).^3)+Zpoly(3).*(log(Fx(fminU+plus)).^2)+Zpoly(4).*log(Fx(fminU+plus))+Zpoly(5));
end;

Zpoly=polyfit(log(Fx(fminU:fmaxU)),log(cS*((2*pi/rwspd(3))^(2/3))*Pu(3,fminU:fmaxU).*(Fx(fminU:fmaxU).^(5/3))),4);
Zslope=Zpoly(4);
if Zslope<0.25,
   Cu2_b(3)=exp(Zpoly(5));
else
   plus=0;
   while 1
      zold=Zslope;
      Zslope=4*Zpoly(1).*(log(Fx(fminU+plus)).^3)+3*Zpoly(2).*(log(Fx(fminU+plus)).^2)+2*Zpoly(3).*log(Fx(fminU+plus))+Zpoly(4);
      plus=plus+1;
      if fminU+plus>=fmaxU,
         break
      end
      if Zslope>zold | abs(Zslope)<0.2,
         break,
      end
   end;
   Cu2_b(3)=exp(Zpoly(1).*(log(Fx(fminU+plus)).^4)+Zpoly(2).*(log(Fx(fminU+plus)).^3)+Zpoly(3).*(log(Fx(fminU+plus)).^2)+Zpoly(4).*log(Fx(fminU+plus))+Zpoly(5));
end;

Zpoly=polyfit(log(Fx(fminU:fmaxU)),log(cS*((2*pi/rwspd(4))^(2/3))*Pu(4,fminU:fmaxU).*(Fx(fminU:fmaxU).^(5/3))),4);
Zslope=Zpoly(4);
if Zslope<0.25,
   Cu2_b(4)=exp(Zpoly(5));
else
   plus=0;
   while 1
      zold=Zslope;
      Zslope=4*Zpoly(1).*(log(Fx(fminU+plus)).^3)+3*Zpoly(2).*(log(Fx(fminU+plus)).^2)+2*Zpoly(3).*log(Fx(fminU+plus))+Zpoly(4);
      plus=plus+1;
      if fminU+plus>=fmaxU,
         break
      end
      if Zslope>zold | abs(Zslope)<0.2,
         break,
      end
   end;
   Cu2_b(4)=exp(Zpoly(1).*(log(Fx(fminU+plus)).^4)+Zpoly(2).*(log(Fx(fminU+plus)).^3)+Zpoly(3).*(log(Fx(fminU+plus)).^2)+Zpoly(4).*log(Fx(fminU+plus))+Zpoly(5));
end;

Zpoly=polyfit(log(Fx(fminU:fmaxU)),log(cS*((2*pi/rwspd(5))^(2/3))*Pu(5,fminU:fmaxU).*(Fx(fminU:fmaxU).^(5/3))),4);
Zslope=Zpoly(4);
if Zslope<0.25,
   Cu2_b(5)=exp(Zpoly(5));
else
   plus=0;
   while 1
      zold=Zslope;
      Zslope=4*Zpoly(1).*(log(Fx(fminU+plus)).^3)+3*Zpoly(2).*(log(Fx(fminU+plus)).^2)+2*Zpoly(3).*log(Fx(fminU+plus))+Zpoly(4);
      plus=plus+1;
      if fminU+plus>=fmaxU,
         break
      end
      if Zslope>zold | abs(Zslope)<0.2,
         break,
      end
   end;
   Cu2_b(5)=exp(Zpoly(1).*(log(Fx(fminU+plus)).^4)+Zpoly(2).*(log(Fx(fminU+plus)).^3)+Zpoly(3).*(log(Fx(fminU+plus)).^2)+Zpoly(4).*log(Fx(fminU+plus))+Zpoly(5));
end;

Zpoly=polyfit(log(Fx(fminU:fmaxU)),log(cS*((2*pi/rwspd(6))^(2/3))*Pu(6,fminU:fmaxU).*(Fx(fminU:fmaxU).^(5/3))),4);
Zslope=Zpoly(4);
if Zslope<0.25,
   Cu2_b(6)=exp(Zpoly(5));
else
   plus=0;
   while 1
      zold=Zslope;
      Zslope=4*Zpoly(1).*(log(Fx(fminU+plus)).^3)+3*Zpoly(2).*(log(Fx(fminU+plus)).^2)+2*Zpoly(3).*log(Fx(fminU+plus))+Zpoly(4);
      plus=plus+1;
      if fminU+plus>=fmaxU,
         break
      end
      if Zslope>zold | abs(Zslope)<0.2,
         break,
      end
   end;
   Cu2_b(6)=exp(Zpoly(1).*(log(Fx(fminU+plus)).^4)+Zpoly(2).*(log(Fx(fminU+plus)).^3)+Zpoly(3).*(log(Fx(fminU+plus)).^2)+Zpoly(4).*log(Fx(fminU+plus))+Zpoly(5));
end;

% W
cS=4;
Zpoly=polyfit(log(Fx(fminU:fmaxU)),log(cS*((2*pi/rwspd(1))^(2/3))*Pw(1,fminU:fmaxU).*(Fx(fminU:fmaxU).^(5/3))),4);
Zslope=Zpoly(4);
if Zslope<0.25,
   Cw2_b(1)=exp(Zpoly(5));
else
   plus=0;
   while 1
      zold=Zslope;
      Zslope=4*Zpoly(1).*(log(Fx(fminU+plus)).^3)+3*Zpoly(2).*(log(Fx(fminU+plus)).^2)+2*Zpoly(3).*log(Fx(fminU+plus))+Zpoly(4);
      plus=plus+1;
      if fminU+plus>=fmaxU,
         break
      end
      if Zslope>zold | abs(Zslope)<0.2,
         break,
      end
   end;
   Cw2_b(1)=exp(Zpoly(1).*(log(Fx(fminU+plus)).^4)+Zpoly(2).*(log(Fx(fminU+plus)).^3)+Zpoly(3).*(log(Fx(fminU+plus)).^2)+Zpoly(4).*log(Fx(fminU+plus))+Zpoly(5));
end;

Zpoly=polyfit(log(Fx(fminU:fmaxU)),log(cS*((2*pi/rwspd(2))^(2/3))*Pw(2,fminU:fmaxU).*(Fx(fminU:fmaxU).^(5/3))),4);
Zslope=Zpoly(4);
if Zslope<0.25,
   Cw2_b(2)=exp(Zpoly(5));
else
   plus=0;
   while 1
      zold=Zslope;
      Zslope=4*Zpoly(1).*(log(Fx(fminU+plus)).^3)+3*Zpoly(2).*(log(Fx(fminU+plus)).^2)+2*Zpoly(3).*log(Fx(fminU+plus))+Zpoly(4);
      plus=plus+1;
      if fminU+plus>=fmaxU,
         break
      end
      if Zslope>zold | abs(Zslope)<0.2,
         break,
      end
   end;
   Cw2_b(2)=exp(Zpoly(1).*(log(Fx(fminU+plus)).^4)+Zpoly(2).*(log(Fx(fminU+plus)).^3)+Zpoly(3).*(log(Fx(fminU+plus)).^2)+Zpoly(4).*log(Fx(fminU+plus))+Zpoly(5));
end;

Zpoly=polyfit(log(Fx(fminU:fmaxU)),log(cS*((2*pi/rwspd(3))^(2/3))*Pw(3,fminU:fmaxU).*(Fx(fminU:fmaxU).^(5/3))),4);
Zslope=Zpoly(4);
if Zslope<0.25,
   Cw2_b(3)=exp(Zpoly(5));
else
   plus=0;
   while 1
      zold=Zslope;
      Zslope=4*Zpoly(1).*(log(Fx(fminU+plus)).^3)+3*Zpoly(2).*(log(Fx(fminU+plus)).^2)+2*Zpoly(3).*log(Fx(fminU+plus))+Zpoly(4);
      plus=plus+1;
      if fminU+plus>=fmaxU,
         break
      end
      if Zslope>zold | abs(Zslope)<0.2,
         break,
      end
   end;
   Cw2_b(3)=exp(Zpoly(1).*(log(Fx(fminU+plus)).^4)+Zpoly(2).*(log(Fx(fminU+plus)).^3)+Zpoly(3).*(log(Fx(fminU+plus)).^2)+Zpoly(4).*log(Fx(fminU+plus))+Zpoly(5));
end;

Zpoly=polyfit(log(Fx(fminU:fmaxU)),log(cS*((2*pi/rwspd(4))^(2/3))*Pw(4,fminU:fmaxU).*(Fx(fminU:fmaxU).^(5/3))),4);
Zslope=Zpoly(4);
if Zslope<0.25,
   Cw2_b(4)=exp(Zpoly(5));
else
   plus=0;
   while 1
      zold=Zslope;
      Zslope=4*Zpoly(1).*(log(Fx(fminU+plus)).^3)+3*Zpoly(2).*(log(Fx(fminU+plus)).^2)+2*Zpoly(3).*log(Fx(fminU+plus))+Zpoly(4);
      plus=plus+1;
      if fminU+plus>=fmaxU,
         break
      end
      if Zslope>zold | abs(Zslope)<0.2,
         break,
      end
   end;
   Cw2_b(4)=exp(Zpoly(1).*(log(Fx(fminU+plus)).^4)+Zpoly(2).*(log(Fx(fminU+plus)).^3)+Zpoly(3).*(log(Fx(fminU+plus)).^2)+Zpoly(4).*log(Fx(fminU+plus))+Zpoly(5));
end;

Zpoly=polyfit(log(Fx(fminU:fmaxU)),log(cS*((2*pi/rwspd(5))^(2/3))*Pw(5,fminU:fmaxU).*(Fx(fminU:fmaxU).^(5/3))),4);
Zslope=Zpoly(4);
if Zslope<0.25,
   Cw2_b(5)=exp(Zpoly(5));
else
   plus=0;
   while 1
      zold=Zslope;
      Zslope=4*Zpoly(1).*(log(Fx(fminU+plus)).^3)+3*Zpoly(2).*(log(Fx(fminU+plus)).^2)+2*Zpoly(3).*log(Fx(fminU+plus))+Zpoly(4);
      plus=plus+1;
      if fminU+plus>=fmaxU,
         break
      end
      if Zslope>zold | abs(Zslope)<0.2,
         break,
      end
   end;
   Cw2_b(5)=exp(Zpoly(1).*(log(Fx(fminU+plus)).^4)+Zpoly(2).*(log(Fx(fminU+plus)).^3)+Zpoly(3).*(log(Fx(fminU+plus)).^2)+Zpoly(4).*log(Fx(fminU+plus))+Zpoly(5));
end;

Zpoly=polyfit(log(Fx(fminU:fmaxU)),log(cS*((2*pi/rwspd(6))^(2/3))*Pw(6,fminU:fmaxU).*(Fx(fminU:fmaxU).^(5/3))),4);
Zslope=Zpoly(4);
if Zslope<0.25,
   Cw2_b(6)=exp(Zpoly(5));
else
   plus=0;
   while 1
      zold=Zslope;
      Zslope=4*Zpoly(1).*(log(Fx(fminU+plus)).^3)+3*Zpoly(2).*(log(Fx(fminU+plus)).^2)+2*Zpoly(3).*log(Fx(fminU+plus))+Zpoly(4);
      plus=plus+1;
      if fminU+plus>=fmaxU,
         break
      end
      if Zslope>zold | abs(Zslope)<0.2,
         break,
      end
   end;
   Cw2_b(6)=exp(Zpoly(1).*(log(Fx(fminU+plus)).^4)+Zpoly(2).*(log(Fx(fminU+plus)).^3)+Zpoly(3).*(log(Fx(fminU+plus)).^2)+Zpoly(4).*log(Fx(fminU+plus))+Zpoly(5));
end;

% H2O

fminq = min(find(Fx>0.7*fsonic/20));
fmaxq = max(find(Fx<3.6*fsonic/20));

