clear all; close all; clc;

VOCALS_DataDir = '/home/disk/eos10/robwood/VOCALS/C130/HighRateData';
FlightFile = 'RF06.20081028.061800_151200.PNI.nc';
FlightPath = [VOCALS_DataDir filesep FlightFile];

temp = nc_varget(FlightPath,'Time');
tsec = length(temp);
temp = repmat(temp,[1,25]);
temp = double(temp) + repmat(0:.04:.96,[tsec,1]);

temp1 = nc_varget(FlightPath,'UIC');
temp2 = nc_varget(FlightPath,'VIC');
temp3 = nc_varget(FlightPath,'WIC');
temp4 = nc_varget(FlightPath,'HGM232');
temp5 = nc_varget(FlightPath,'THETA');

telapse = zeros(tsec*25,1);
Uvel = zeros(tsec*25,1);
Vvel = zeros(tsec*25,1);
Wvel = zeros(tsec*25,1);
HGM232 = zeros(tsec*25,1);

for j = 0:(tsec-1)
    telapse((1+25*j):(25*(j+1))) = temp((j+1),:);
    Uvel((1+25*j):(25*(j+1))) = temp1((j+1),:);
    Vvel((1+25*j):(25*(j+1))) = temp2((j+1),:);
    Wvel((1+25*j):(25*(j+1))) = temp3((j+1),:);
    HGM232((1+25*j):(25*(j+1))) = temp4((j+1),:);
    Theta((1+25*j):(25*(j+1))) = temp5((j+1),:);
end

clear temp*

sc_start = 10330*25;
sc_end = 12040*25;
cb_start = 12240*25;
cb_end = 14070*25;
cl_start = 14240*25;
cl_end = 16010*25;
saw_start = 17040*25;
saw_end = 18620*25;

Usc = Uvel(sc_start:sc_end);
Vsc = Vvel(sc_start:sc_end);
Wsc = Wvel(sc_start:sc_end);

Ucb = Uvel(cb_start:cb_end);
Vcb = Vvel(cb_start:cb_end);
Wcb = Wvel(cb_start:cb_end);

Ucl = Uvel(cl_start:cl_end);
Vcl = Vvel(cl_start:cl_end);
Wcl = Wvel(cl_start:cl_end);

Usaw = Uvel(saw_start:saw_end);
Vsaw = Vvel(saw_start:saw_end);
Wsaw = Wvel(saw_start:saw_end);

% Sample analysis block

X1 = Wcl(ceil(3/4*end):end); % Series goes here

x1_att.name = 'Wcl'; % Series attributes
x1_att.N = length(Wcl);
x1_att.delta_t = .04;
x1_att.delta_t_units = 'sec';
%x1_att.start = telapse(sc_start);
%x1_att.end = telapse(sc_end);
x1_att.units = 'm sec^{-1}';
delta_t1 = x1_att.delta_t;

[W1Jt, V1Jt] = modwt(X1, 'la8', 13, 'reflection');
[wvar1, CI_wvar1] = modwt_wvar(W1Jt);

X2 = Wcl(1:ceil(1/4*end)); % Series goes here

x2_att.name = 'Wcl'; % Series attributes
x2_att.N = length(Wcl);
x2_att.delta_t = .04;
x2_att.delta_t_units = 'sec';
%x2_att.start = telapse(sc_start);
%x2_att.end = telapse(sc_end);
x2_att.units = 'm sec^{-1}';

delta_t2 = x2_att.delta_t;

[W2Jt, V2Jt] = modwt(X2, 'la8', 13, 'reflection');
[wvar2, CI_wvar2] = modwt_wvar(W2Jt);

X3 = Wcb(1:ceil(1/4*end));

x3_att.name = 'Wcb'; % Series attributes
x3_att.N = length(Wcb);
x3_att.delta_t = .04;
x3_att.delta_t_units = 'sec';
%x3_att.start = telapse(cb_start);
%x3_att.end = telapse(cb_end);
x3_att.units = 'm sec^{-1}';

delta_t3 = x3_att.delta_t;

[W3Jt, V3Jt] = modwt(X3, 'la8', 13, 'reflection');
[wvar3, CI_wvar3] = modwt_wvar(W3Jt);

X4 = Wcb(ceil(3/4*end):end);

x4_att.name = 'Wcb'; % Series attributes
x4_att.N = length(Wcb);
x4_att.delta_t = .04;
x4_att.delta_t_units = 'sec';
%x4_att.start = telapse(cb_start);
%x4_att.end = telapse(cb_end);
x4_att.units = 'm sec^{-1}';

delta_t4 = x4_att.delta_t;

[W4Jt, V4Jt] = modwt(X4, 'la8', 13, 'reflection');
[wvar4, CI_wvar4] = modwt_wvar(W4Jt);

X5 = Wsc(ceil(3/4*end):end);

x5_att.name = 'Wsc'; % Series attributes
x5_att.N = length(Wsc);
x5_att.delta_t = .04;
x5_att.delta_t_units = 'sec';
%x5_att.start = telapse(sc_start);
%x5_att.end = telapse(sc_end);
x5_att.units = 'm sec^{-1}';

delta_t5 = x5_att.delta_t;

[W5Jt, V5Jt] = modwt(X5, 'la8', 13, 'reflection');
[wvar5, CI_wvar5] = modwt_wvar(W5Jt);

X6 = Wsc(1:ceil(1/4*end));

x6_att.name = 'Wcb'; % Series attributes
x6_att.N = length(Wcb);
x6_att.delta_t = .04;
x6_att.delta_t_units = 'sec';
%x6_att.start = telapse(cb_start);
%x6_att.end = telapse(cb_end);
x6_att.units = 'm sec^{-1}';

delta_t6 = x6_att.delta_t;

[W6Jt, V6Jt] = modwt(X6, 'la8', 13, 'reflection');
[wvar6, CI_wvar6] = modwt_wvar(W6Jt);



% plot code
figure,

subplot(1,3,1)
semilogy(1:13,wvar1,'b'), hold on
semilogy(1:13,CI_wvar1(:,1),'b--')
semilogy(1:13,CI_wvar1(:,2),'b--')
semilogy(1:13,wvar2,'r'), hold on
semilogy(1:13,CI_wvar2(:,1),'r--')
semilogy(1:13,CI_wvar2(:,2),'r--')

subplot(1,3,2)
semilogy(1:13,wvar3,'b'), hold on
semilogy(1:13,CI_wvar3(:,1),'b--')
semilogy(1:13,CI_wvar3(:,2),'b--')
semilogy(1:13,wvar4,'r'), hold on
semilogy(1:13,CI_wvar4(:,1),'r--')
semilogy(1:13,CI_wvar4(:,2),'r--')

subplot(1,3,3)
semilogy(1:13,wvar5,'b'), hold on
semilogy(1:13,CI_wvar5(:,1),'b--')
semilogy(1:13,CI_wvar5(:,2),'b--')
semilogy(1:13,wvar6,'r'), hold on
semilogy(1:13,CI_wvar6(:,1),'r--')
semilogy(1:13,CI_wvar6(:,2),'r--')

% kolmogorov spectra
[CJ1, f1, CI_CJ1, f_band1] = modwt_wvar_psd(wvar1,delta_t1,CI_wvar1);
[CJ2, f2, CI_CJ2, f_band2] = modwt_wvar_psd(wvar2,delta_t2,CI_wvar2);
[CJ3, f3, CI_CJ3, f_band3] = modwt_wvar_psd(wvar3,delta_t3,CI_wvar3);
[CJ4, f4, CI_CJ4, f_band4] = modwt_wvar_psd(wvar4,delta_t4,CI_wvar4);
[CJ5, f5, CI_CJ5, f_band5] = modwt_wvar_psd(wvar5,delta_t5,CI_wvar5);
[CJ6, f6, CI_CJ6, f_band6] = modwt_wvar_psd(wvar6,delta_t6,CI_wvar6);

figure;
subplot(1,3,1)
loglog(f1,CJ1,'b*'), hold on
loglog(f2,CJ2,'r*')
for j = 1:length(f1)
    rectangle('Position',[f_band1(j,1),CI_CJ1(j,1),(f_band1(j,2)-f_band1(j,1)),CI_CJ1(j,2)-CI_CJ1(j,1)])
    rectangle('Position',[f_band2(j,1),CI_CJ2(j,1),(f_band2(j,2)-f_band2(j,1)),CI_CJ2(j,2)-CI_CJ2(j,1)])
    ylim([.00001 100])
    loglog((.1:.1:10),.1.*(.1:.1:10).^-(5/3),'k')
end

subplot(1,3,2)
loglog(f3,CJ3,'b*'), hold on
loglog(f4,CJ4,'r*')
for j = 1:length(f3)
    rectangle('Position',[f_band3(j,1),CI_CJ3(j,1),(f_band3(j,2)-f_band3(j,1)),CI_CJ3(j,2)-CI_CJ3(j,1)])
    rectangle('Position',[f_band4(j,1),CI_CJ4(j,1),(f_band4(j,2)-f_band4(j,1)),CI_CJ4(j,2)-CI_CJ4(j,1)])
    ylim([.00001 100])
    loglog((.1:.1:10),.1.*(.1:.1:10).^-(5/3),'k')
end

subplot(1,3,3)
loglog(f5,CJ5,'b*'), hold on
loglog(f6,CJ6,'r*')
for j = 1:length(f5)
    rectangle('Position',[f_band5(j,1),CI_CJ5(j,1),(f_band5(j,2)-f_band5(j,1)),CI_CJ5(j,2)-CI_CJ5(j,1)])
    rectangle('Position',[f_band6(j,1),CI_CJ6(j,1),(f_band6(j,2)-f_band6(j,1)),CI_CJ6(j,2)-CI_CJ6(j,1)])
    ylim([.00001 100])
    loglog((.1:.1:10),.1.*(.1:.1:10).^-(5/3),'k')
end

close all

% mra analysis

X = Wsaw;

wtfname = 'la8';
J0 = 10;
boundary = 'reflection';

[WJt, VJt, att] = modwt(X, wtfname, J0, boundary);
[DJt, SJt] = imodwt_mra(WJt, VJt, att);

% Setup x-axis for plotting.
delta_t = 1 / 25;   % sampling interval in seconds 
xaxis = (1:length(X))* delta_t;
xlabel_str = 't (sec)';

% Setup plotting parameters
title_str = {};
title_str(1) = {'MODWT multiresolution analysis of Vertical Velocity time series'};
title_str(2) = {['Wavelet: ' wtfname ',  NLevels: ', int2str(J0), ...
                 ',  Boundary: ', boundary]};

XplotAxesProp.XLim       = [0 1580];
XplotAxesProp.XTick      = [0:100:1580];
XplotAxesProp.XTickLabel = [0:100:1580];
XplotAxesProp.XMinorTick = 'off';

XplotAxesProp.YLim       = [-5 5];
XplotAxesProp.YTick      = [-5 -3 0 3 5];
XplotAxesProp.YTickLabel = [-5 -3 0 3 5];

% Plot MRA coefficients

[hMRAplotAxes, hXplotAxes] = plot_imodwt_mra(DJt, SJt, X, att);
hold(hXplotAxes), plot(hXplotAxes,HGM232(saw_start:saw_end)./2000,'k--');
plot(hXplotAxes,(Theta(saw_start:saw_end)-289)./5,'r--');