clear all
close all

addpath(genpath('/Users/Mau/Documents/MATLABhelp_new')); 
% print_setup.m for saving plots
% cbdate, maximize
% spectra.m

%% Purpose: Compare w_ctd from SBE41cp on float (2-m binnned vs 1Hz data
% from .eng files) and dp/dt from 64Hz pressure sensor on MicroRider

% Directory for CTD data:
path_ctd='/Users/Mau/Documents/MBARIprojects/CPF/CPF003/Data/matfiles/ParagonTests/MICRORIDER/2021-10-8T16-19-7/'; % profile and w_ctd data in .mat files
% Directory for MR data:
path_mr='/Users/Mau/Documents/MBARIprojects/MR/TurbulenceData/ParagonTests/2021-10-08/';

path_ctd='/Users/Mau/Documents/MBARIprojects/CPF/CPF003/Data/matfiles/ParagonTests/MICRORIDER/2021-10-13/';
path_mr='/Users/Mau/Documents/MBARIprojects/MR/TurbulenceData/ParagonTests/2021-10-13/'; % only 1 prof returned data in .msg files

path_ctd='/Users/Mau/Documents/MBARIprojects/CPF/CPF003/Data/matfiles/ParagonTests/MICRORIDER/2021-10-15T16-10-6/';
path_mr='/Users/Mau/Documents/MBARIprojects/MR/TurbulenceData/ParagonTests/2021-10-15/';

mission_date='2021-10-19T16-3-9'; % 8profiles, half day; to be able to read CTD data keep time. mission_date='2021-10-19'; %
path_ctd=['/Users/Mau/Documents/MBARIprojects/CPF/CPF003/Data/matfiles/ParagonTests/MICRORIDER/',mission_date,'/'];
path_mr=['/Users/Mau/Documents/MBARIprojects/MR/TurbulenceData/ParagonTests/',mission_date(1:10),'/'];

% now specifying mission date instead, and code will go look into
% respective directories
mission_date='2021-11-16'; % mission_date='2021-10-16T?'; % Overnight mission, 18 profiles
% mission_date='2021-10-19'; % Half day mission, 8 profiles
% path_ctd=['/Users/Mau/Documents/MBARIprojects/CPF/CPF003/Data/matfiles/ParagonTests/MICRORIDER/',mission_date,'*/']; % added * to look into original file path but HR str doesn't have w_ctd
path_ctd=['/Users/Mau/Documents/MBARIprojects/CPF/CPF003/Data/matfiles_mr/ParagonTests/MICRORIDER/',mission_date,'/'];
path_mr=['/Users/Mau/Documents/MBARIprojects/MR/TurbulenceData/ParagonTests/',mission_date,'/matfiles/'];

%% Define Directories to save merged data from CPF+MR
path_out=['/Users/Mau/Documents/MBARIprojects/CPF-MR/ParagonTests/',mission_date,'/matfiles/'];
% Define Directories to save plots (with MR data)
path_plots=[path_out(1:end-9),'plots','/'];

%% Create directory to save data & plots (if it doesn't exist already)
if isfolder(path_out)==0 & numel(dir(path_out))<=2 % will create directory only if it is empty!
    mkdir(path_out)
elseif isfolder(path_plots)==0 & numel(dir(path_plots))<=2 % will create directory only if it is empty!
    mkdir(path_plots)
end

%%
% NOTE: should consolidate folder naming conventions so data from the same
% mission matches. "done" for MR data.

%% Define colors for plots
cbar=cbrewer('seq','Blues',101);
% cbar=cbar(5:5:end,:);
lblue=cbar(20,:);
blue=cbar(60,:);

Reds=cbrewer('seq','Reds',101);
lred=Reds(20,:);
red=Reds(50,:);

Greens=cbrewer('seq','Greens',101);
lgreen=Greens(20,:);
green1=Greens(40,:);
green2=Greens(60,:);

lgrey=[0.7 0.7 0.7];

% Get default color order
for p=1:8
h = plot(1:10,1:10); hold on
Colors = get(h,'Color');
end
close all


%% List files to process and preallocate memory
files=dir([path_ctd,'*.mat']); % list matfiles with CPF data for a given mission
nfiles=length(files);
nfiles=nfiles-1; % Note, last CPF profile usually is garbage (if mission was recovered before writing data to file)

Lat=nan(nfiles,1);
Lon=nan(nfiles,1);
Time=nan(nfiles,1);

W_1hz=nan(800,nfiles);
P_1hz=nan(800,nfiles);

W_2mb=nan(200,nfiles);
D_2mb=nan(200,nfiles);
Nbin=nan(200,nfiles);

RHO=nan(200,nfiles);
CT=nan(200,nfiles);
SA=nan(200,nfiles);

N2gsw=nan(200,nfiles);
Ngsw=nan(200,nfiles);
p_mid=nan(200,nfiles);
N2=nan(200,nfiles);
N=nan(200,nfiles);

epsilon2=nan(300,nfiles);
mr_paxis=nan(300,nfiles);

W_mr=nan(300,nfiles);
P_mr=nan(300,nfiles);

W_1hzlp=nan(800,nfiles); %added from Gene's ENG file. units fixed and cutoff filter known
P_1hzlp=nan(800,nfiles);

ppl=1; % 0 do not plot profiles, 1 for plotting

%% Process each file/profile

for f=1:nfiles

    % NOte: should fix how file names when saved in .mat format to match
    % CPF matfile data ('.i.e. '0' extra)
    
    % Load CTD and MR data
    if f<10
        load([files(f).folder,'/','CPF003.00',num2str(f),'.mat']); % CPF data
        load([path_mr,'CPF003.00',num2str(f),'.MR.mat']); % MR data
    elseif (f>=10 && f<=99)
        load([files(f).folder,'/','CPF003.0',num2str(f),'.mat']); % CPF data
        load([path_mr,'CPF003.00',num2str(f),'.MR.mat']); % MR data
    end
    % Note: should fix how file names when saved in .mat format to match
    % CPF matfile data ('.i.e. '0' extra)

    maxp=100;
    
    % Compute W from SBE41cp 1Hz pressure signal (wihout smoothing, the RSI way)
    P1hz=ENG.P1hz; % LowPass press, ask Gene what type of filter
    P1hz=P1hz(~isnan(P1hz));
    
    dP1Hz=diff(P1hz); dt=1;
    w1hz=-1*dP1Hz./dt;
    P1hz(w1hz==0)=[]; 
    w1hz(w1hz==0)=[]; %check w/Gene, is this because of the sensitivity of sensor or writing values twice in the .eng files?
    
    % NOTE: it makes sense to remove zeros (consecutive values with equal
    % pressure. Otherwise, a discontinuity will show up in the w profile at~80m)
    W_1hz(1:length(w1hz),f)=flipud(w1hz); %flip, such that surface is on top of matrix
    P_1hz(1:length(w1hz),f)=flipud(P1hz(1:end-1));


    % Alternatively, use W from SBE41cp computed by Gene
    W1hzlp_g=ENG.Wlp; % LowPass press, ask Gene what type of filter. looks like fcutoof/pass =0.1000
    W1hzlp_g=W1hzlp_g(~isnan(W1hzlp_g));
    P1hzlp_g=ENG.Plp;
    P1hzlp_g=P1hzlp_g(~isnan(P1hzlp_g));
    W1hzlp_g=3.*W1hzlp_g; % fixing bug in code
    
%     W1hz_g=ENG.W; % Slope of LP Press, right units but too smoothed
%     
    % How does RSI calculate platform speed?
%     fc = 0.5;
%     fs_slow=63.99;
%     [b,a]=butter(4,fc/(fs_slow/2));
%     dP_dt = gradient(P_slow,1/fs_slow);
%     dP_dt_LP = filtfilt(b,a,dP_dt);
    
    % Save all profiles from a given mission in a matrix
    W_2mb(1:length(HR.w_ctd1),f)=HR.w_ctd1;
    D_2mb(1:length(HR.w_ctd1),f)=HR.d_axis;
    Nbin(1:length(HR.N),f)=HR.N;
    % NOTE: wctd1 compares much better with w from MR P

    % Save all profiles from 1Hz data  given mission in a matrix
    W_1hzlp(1:length(W1hzlp_g),f)=W1hzlp_g;
    P_1hzlp(1:length(W1hzlp_g),f)=P1hzlp_g;

    % NOTE: wctd1 compares much better with w from MR P
    
    
    if ppl==1
    plot(w1hz,P1hz(1:end-1),'-o','color',green1); % computed doing the difference between samples at depth
    hold on
%    plot(ENG.W,ENG.Plp,'-o','color',green2); % computed by Gene, SlopeLPFpressure, units OK but too smooth
%    plot(3*ENG.Wlp,ENG.Plp,'-o','color','r'); % wrong units? these platVel don't make much sense
%     plot(3*ENG.Wdiff,ENG.Plp,'-o','color','g'); % wrong units? off by of 0.333?
%   Mistery solved (see notes in Evernote/notebook). SlopeLPFpressure has
%   the correct units (m/s).
%   All other PV estimates should be multiplied by 3 to get the units
%   right. ENG.Wdiff looks pretty much like what I computed above
    plot(W1hzlp_g,P1hzlp_g,'-o','color',green2);
    
    plot(MR.speed,MR.P,'-','color',red); %cbar(10,:)); % from MicroRider (LowPassed)
    plot(HR.w_ctd1,HR.d_axis,'color',blue); %cbar(10,:)); %SBE41cp 2-m binned
    
    axis ij

    legend('SBE41cp 1Hz','SBE41cp 1Hz LP0.1','MR P 64Hz LP0.5','SBE41cp 2m binned',...
        'location','W')

    xlim([range(W_2mb(:))])
    ylim([0 maxp])
    xlabel('w_{dpdt} [m/s]');

    print_setup(20,40)
    print([path_plots,'w_2sensors_p_',num2str(f),'_wctd1.jpg'],'-djpeg');

    close all
    
    end
    
    % Compute N2 (using gsw function)
    %     if (LAT{f}(p)<90 && LAT{f}(p)>-90) % quick fix to compute MLD where lat/lon is within range only
    g=9.81; % [m/s]
    [N2gsw(1:length(HR.P)-1,f), p_mid(1:length(HR.P)-1,f)] = gsw_Nsquared(HR.SA,HR.CT,HR.P,HR.LAT); % Note this is in [rad2 s-2 ]. 
    % To convert frequency to hertz, cycles sec-1, divide the frequency by 2π, ie N/(2π).
    
    nnan=find(~isnan(N2gsw(:,f)));
    Ngsw(nnan,f)=sqrt(N2gsw(nnan,f)); % in rad/s
    
    % Compute N2 from density (i.e. to sort profile and avoid
    % negative N^2)
    drho=diff(sort(HR.RHO));
    dz=diff(HR.DEPTH);
    drhodz=drho./dz;
    rho_ref=nanmean(HR.RHO); % mean density for each profile
    N2(1:length(HR.P)-1,f)=drhodz.*[g./rho_ref];
    
    N(1:length(HR.P)-1,f)=sqrt(N2(1:length(HR.P)-1,f));
    % NOTE: already in HR! local version
    
    % Save epsilon and platform velocity from the MR
    ne=length(MR.e(2,:));
    epsilon2(1:ne,f)=flipud(MR.e(2,:)');
    mr_paxis(1:ne,f)=flipud(MR.P);
    
    W_mr(1:ne,f)=flipud(MR.speed);
    
    % Save Coordinates for plotting (from HR CPF structure)
    Lat(f)=HR.LAT;
    Lon(f)=HR.LON;
    Time(f)=HR.TIME;
    
    % SAVE other sensors variables:
    RHO(1:length(HR.P),f)=HR.RHO;
    CT(1:length(HR.P),f)=HR.CT;
    SA(1:length(HR.P),f)=HR.SA;
    
end

%% Plot Buoyancy Frequency, N
close all

LocalTime=Time-8/24;

[X,Y]= meshgrid(LocalTime,D_2mb(:,1));

h=pcolor(X,Y,log10(N2)); % in Hz,cycles per sec
h.EdgeColor = 'none'; axis ij
% hold on

% ylim([1 50])
colormap(flipud(cbrewer('div','Spectral',50)));
% colormap(flipud(cbrewer('div','RdYlBu',21)));
hc=colorbar; 
set(get(hc,'Title'),'String','[s^{-2}]','fontsize',20);
datetick('x','HH:MM','keeplimits');

ylabel('Depth [m]')
xlabel('Local Time')
% xlabel('Profile Number')
% xlabel('Longitude [^{\circ} W]');

title('log(N^2)')
caxis([-6 -3])

print_setup(50,30)
print([path_plots,'N2_logscale_all.jpg'],'-djpeg');

% regresion(N2gsw(:),N2(:),2); great agreement between the gsw code and calculating from scratch
% corr of 1!
% will use N2 because it doesn't have any negative values (computed from
% sorted RHO profiles)

%% Plot Buoyancy Frequency, N2 in (rad/s)^2
close all

LocalTime=Time-8/24;

[X,Y]= meshgrid(LocalTime,D_2mb(:,1));

h=pcolor(X,Y,log10(N2.*(2*pi)^2)); % in radians per second (to compare with GSW)
h.EdgeColor = 'none'; axis ij
% hold on

% ylim([1 50])
colormap(flipud(cbrewer('div','Spectral',50)));
% colormap(flipud(cbrewer('div','RdYlBu',21)));
hc=colorbar; 
set(get(hc,'Title'),'String','[(rad/s)^{2}]','fontsize',20);
datetick('x','HH:MM','keeplimits');

ylabel('Depth [m]')
xlabel('Local Time')
% xlabel('Profile Number')
% xlabel('Longitude [^{\circ} W]');

title('log(N^2)')
% caxis([-6 -3])
caxis([-3 -1])

print_setup(50,30)
print([path_plots,'N2_logscale_all_radpsec2.jpg'],'-djpeg');

%% Plot Buoyancy Frequency (in Hz, 1/seconds)
close all

h=pcolor(X,Y,N);
h.EdgeColor = 'none'; axis ij
% hold on

% ylim([1 50])
colormap(flipud(cbrewer('div','Spectral',50)));
% colormap(flipud(cbrewer('div','RdYlBu',21)));
hc=colorbar; 
set(get(hc,'Title'),'String','[s^{-1}]','fontsize',20);

datetick('x','HH:MM','keeplimits');

ylabel('Depth [m]')
xlabel('Local Time')

title('N')
print_setup(50,30)
% print([path_plots,'N_all.jpg'],'-djpeg');

Nmin=nanmin(N(:))
Nmax=nanmax(N(:))

%% Plot Buoyancy period (in minutes!)
close all

h=pcolor(X,Y,(1./N)./60);
h.EdgeColor = 'none'; axis ij
% hold on

% ylim([1 50])
colormap(flipud(cbrewer('div','Spectral',50)));
% colormap(flipud(cbrewer('div','RdYlBu',21)));
hc=colorbar; 
set(get(hc,'Title'),'String','[s]','fontsize',20);
set(get(hc,'Title'),'String','[min]','fontsize',20);

datetick('x','HH:MM','keeplimits');

ylabel('Depth [m]')
xlabel('Local Time')

% caxis([1 1000])

title('Buoyancy Period, T_{N} (minutes)')
print_setup(50,30)
print([path_plots,'Tn_minutes_all.jpg'],'-djpeg');

% Plot w profile data

%% Plot W (raw ascent rates)
close all

LocalTime=Time-8/24;
[X,Y]= meshgrid(LocalTime,D_2mb(:,1));
Z=W_2mb-nanmean(W_2mb);


% [X,Y]= meshgrid(LocalTime,P_1hzlp(:,1));
% Z=W_1hzlp-nanmean(W_1hzlp);

% NOTE this plots are an approx.. Paxis is not the same for every profile.

h=pcolor(X,Y,Z);
h.EdgeColor = 'none'; axis ij
% hold on

caxis([-0.14 0.14])
% ylim([1 50])
% colormap(flipud(cbrewer('div','Spectral',50)));
% colormap(flipud(cbrewer('div','RdYlBu',21)));
colormap(flipud(cbrewer('div','RdBu',21)));
hc=colorbar; 
set(get(hc,'Title'),'String','[m/s]','fontsize',20);

datetick('x','HH:MM','keeplimits');

ylabel('Depth [m]')
xlabel('Local Time')

% caxis([1 1000])

print_setup(50,30)
print([path_plots,'w_CTD_2mb_all.jpg'],'-djpeg');
% print([path_plots,'w_CTD_2mb_1Hz_MR_all.jpg'],'-djpeg');

%% Compare Spectra from all 3 signals
%% Spectral estimates for w_ctd

% Take all data to the same high res grid for meaningful comparisons of
% spectral estimates!

dz=0.1; %m @ PV ~0.18 m/s, dt~0.5542
dz=0.002 %m @ PV ~0.18 m/s, dt~1/64;
dz=0.18*(1/64); dst=1./64;
% dz=0.35 % meandiff in mr_paxis. NOTE: data resolution had been reduced to 2s, or 0.35dbar
% dz=0.175 % meandiff in CTD 1Hz data.

% use the sam zbounds and dz to get all data with the same size
pmr_interp=[round(nanmin(mr_paxis)):dz:round(nanmax(mr_paxis(:)))]';
p1hz_interp=pmr_interp;
p2mb_interp=pmr_interp;
p1hzlp_interp=pmr_interp;

% dz=2;
% p2mb_interp=[round(nanmin(D_2mb(:))):dz:round(nanmax(D_2mb(:)))]';
W_2mb_interp=nan(length(p2mb_interp),f);
N_interp=nan(length(p2mb_interp),f);

% dz=0.15;
% p1hz_interp=[round(nanmin(P_1hz(:))):dz:round(nanmax(P_1hz(:)))]';
W_1hz_interp=nan(length(p1hz_interp),f);

% dz=0.0156; %1/64
% pmr_interp=[round(nanmin(mr_paxis)):dz:round(nanmax(mr_paxis(:)))]';
W_mr_interp=nan(length(pmr_interp),f);

% p1hzlp_interp=[round(nanmin(P_1hzlp(:))):dz:round(nanmax(P_1hzlp(:)))]';
W_1hzlp_interp=nan(length(pmr_interp),f);

% Interpolating in P space
for p=1:f
    
    indnani=find(~isnan(W_2mb(:,p)),1,'first');
    indnanf=find(~isnan(W_2mb(:,p)),1,'last');
    y=interp1(D_2mb(indnani:indnanf,p),W_2mb(indnani:indnanf,p),p2mb_interp,'linear');
    W_2mb_interp(1:length(y),p)=y;
    
    [xu, uind]=unique(P_1hz(:,p));
    wu=W_1hz(uind,p);
    indnani=find(~isnan(xu),1,'first');
    indnanf=find(~isnan(xu),1,'last');
    y=interp1(xu(indnani:indnanf),wu(indnani:indnanf),p1hz_interp,'linear');
    W_1hz_interp(1:length(y),p)=y;

    [xu, uind]=unique(P_1hzlp(:,p)); %idem as above for Gene's computed W and LP filtered
    wu=W_1hzlp(uind,p);
    indnani=find(~isnan(xu),1,'first');
    indnanf=find(~isnan(xu),1,'last');
    y=interp1(xu(indnani:indnanf),wu(indnani:indnanf),p1hzlp_interp,'linear');
    W_1hzlp_interp(1:length(y),p)=y;
    
    indnani=find(~isnan(W_mr(:,p)),1,'first');
    indnanf=find(~isnan(W_mr(:,p)),1,'last');
    y=interp1(mr_paxis(indnani:indnanf,p),W_mr(indnani:indnanf,p),pmr_interp,'linear');
    W_mr_interp(1:length(y),p)=y;

    indnani=find(~isnan(N(:,p)),1,'first');
    indnanf=find(~isnan(N(:,p)),1,'last');
    y=interp1(D_2mb(indnani:indnanf,p),N(indnani:indnanf,p),p2mb_interp,'linear');
    N_interp(1:length(y),p)=y;
    
end

%% Check that interpolation is working 

close all
plot(p2mb_interp,W_2mb_interp(:,p),'*r'); hold on; plot(D_2mb(:,p),W_2mb(:,p),'--k');

close all
subplot(211)
plot(p1hz_interp,W_1hz_interp(:,p),'*r'); hold on; plot(P_1hz(:,p),W_1hz(:,p),'--k'); % spikes in the 1Hz Pres data from CPF when going through the thermocline
xlim([0 100])
subplot(212)
plot(HR.P,HR.T);

close all
plot(pmr_interp,W_mr_interp(:,p),'*r'); hold on; plot(mr_paxis(:,p),W_mr(:,p),'--k'); % spikes in the 1Hz Pres data from CPF when going through the thermocline
xlim([0 100])

close all
plot(pmr_interp,W_1hzlp_interp(:,p),'*r'); hold on; plot(mr_paxis(:,p),W_1hzlp(:,p),'--k'); 
xlim([0 100])


%% SAVE all w data and corresponding P in one .mat file
save([path_out,'w_CTD_MR_allprofiles_v2.mat'],...
    'W_1hz','W_2mb','W_mr','W_1hzlp','P_1hz','D_2mb','mr_paxis','P_1hzlp',...
    'W_1hz_interp','W_2mb_interp','W_mr_interp','p1hz_interp','p2mb_interp','pmr_interp','N_interp', ...
    'W_1hzlp_interp','p1hzlp_interp',...
    'N2','N','RHO','CT','SA','Nbin','Time', ... %);
    'epsilon2');
% NOTE: % _v2.mat shoud have all vars using HR.w_ctd1 and saving epsilon2
% data to bin and compare in test_pwech.m, and saving now LP filtered W1hz
% computed by Gene.

%% Compute Spectra for w_ctd (2-m binned data)
% NOTE: similar analysis done in test_pwelch.m

close all

% dst=0.55;
% see note above: for dz=0.1; %m @ mean PV ~0.18 m/s, dt~0.5542

% dst=10; % round(nanmean(nanmean(Nbin(1:40,1:end)))); % % in sec for 1Hz sampling rate.

% Spectra code uses a mean dt
time=dst; %10; %/3600;
% data=W_2mb(:,~isnan(W_2mb(1,:))); % remove profiles w/NaN
% indpi=1; indpf=40;
% data=data(indpi:indpf,p); %  truncate to omit bottom layer (BC ON)

% Nyquist and fundamental freq
fNy=1./(2*dst);

z1=10; z2=80; % bounds to cutoff profile data
fTn=1./((z2-z1)./0.18);

% spectra using interpolated data instead
data=W_2mb_interp; %(:,~isnan(W_2mb_interp(1,:))); % remove profiles w/NaN
indpi=searchclosest_m(p2mb_interp,z1);
indpf=searchclosest_m(p2mb_interp,z2);
data=data(indpi:indpf,:); % % truncate to omit bottom layer (BC ON)
method='interp';

[nz,np]=size(data);

nu=np*2; % !! to compute spectra for each profile, and average across profiles
         % it needs to be 2*Nprofiles
varunits='m/s';
timedt='s';
color='k';
[kn_2mbin,meanSE_2mbin,dk,varcheckSat,h,SE_2mbin]=spectra(time,data,nu,varunits,timedt,color);
hold on

xlim([fTn fNy])

% Compute Spectra for 64Hz dp/dt (smoothed by RSI !!)
% dst=1/64; % % in sec for 64Hz sampling rate.
% % dst=0.0982; % % Effective resolution? for a generic low-pass filter: pi* DeltaR_raw/w_cuttoff
% dst=2; %s 0.35/0.15;

% Spectra code uses a mean dt
time=dst; %10; %/3600;
% data=MR.speed(:,~isnan(MR.speed(1,:))); % remove profiles w/NaN
% data=data(:,:); % truncate to omit bottom layer (BC ON)

% if using interpolated data
data=W_mr_interp; %(:,~isnan(W_mr_interp(100,:))); % remove profiles w/NaN
indpi=searchclosest_m(pmr_interp,z1);
indpf=searchclosest_m(pmr_interp,z2);
data=data(indpi:indpf,:); % truncate to omit bottom layer (BC ON)

[nz,np]=size(data);

nu=np*2; % to compute spectra for each profile, and average across profiles
varunits='m/s';
timedt='s';
color='r';
[kn_MRw,meanSE_MRw,dk,varcheckSat,h,SE_MRw]=spectra(time,data,nu,varunits,timedt,color);

% Compute Spectra for 1Hz dp/dt
% dst=1; % % in sec for 1Hz sampling rate.

% Spectra code uses a mean dt
time=dst; %10; %/3600;
% data=W_1hz(:,~isnan(W_1hz(1,:))); % remove profiles w/NaN
% data=data(1:500,:); % truncate to omit bottom layer (BC ON)

data=W_1hz_interp; %(:,~isnan(W_1hz_interp(50,:))); % remove profiles w/NaN
indpi=searchclosest_m(p1hz_interp,z1);
indpf=searchclosest_m(p1hz_interp,z2);
data=data(indpi:indpf,:); % truncate to omit bottom layer (BC ON)

[nz,np]=size(data);

nu=np*2; % to compute spectra for each profile, and average across profiles
varunits='m/s';
timedt='s';
color='g';
[kn_W1hz,meanSE_W1hz,dk,varcheckSat,h,SE_W1hz]=spectra(time,data,nu,varunits,timedt,color);


% Axis
% xlim([1*10^(-3) 4*10^(1)])

% Change color of lines overlaying plot

kn_W1hz_eff=1;
kn_2mbin_eff=1/(2*10);
kn_MRw_eff=1/(2*2.3); % effective res for butterworthfilter?
%using ~2s (estimated from dp=0.35dbar and mean PV from MR of 0.15 m/s)

h3=loglog(kn_W1hz,meanSE_W1hz,'color',lgreen);
h2=loglog(kn_2mbin,meanSE_2mbin,'color',lblue);
h1=loglog(kn_MRw,meanSE_MRw,'color',lred);

h3=loglog(kn_W1hz(kn_W1hz<kn_W1hz_eff),meanSE_W1hz(kn_W1hz<kn_W1hz_eff),'color',green1);
h2=loglog(kn_2mbin(kn_2mbin<=kn_2mbin_eff),meanSE_2mbin(kn_2mbin<=kn_2mbin_eff),'color',blue);
h1=loglog(kn_MRw(kn_MRw<=kn_MRw_eff),meanSE_MRw(kn_MRw<=kn_MRw_eff),'color',red);

legend('','','','','','','CTD 1Hz','2-m binned CTD 1Hz','P from MR 64Hz','');

% maximize;
print_setup(40,30)
% print([path_plots,'Spectra_wctd_2mbinnedCTD_MRp_z1-80m_',method,'.eps'],'-depsc2');
print([path_plots,'Spectra_wctd_2mbinnedCTD_MRp_z',num2str(z1),'-',num2str(z2),'m_',method,'.eps'],'-depsc2');

% % To plot all spectra (not just teh average)
% hold on
% loglog(knSat,SESat,'linewidth',1); %plot all spectra in random colors
% 
% % plot all spectra in gray shades
% % cbar=colormap(gray); %gray shades
% % cbar=cbar(5:10:end,:); cbar=cbar(1:4,:);
% % for p=1:np
% % loglog(knSat,SESat(:,p),'color',cbar(p,:),'linestyle','--','linewidth',0.5);
% % end
% 
% [knSat,meanSESat,dk,varcheckSat,h,SESat]=spectra(time,data,nu,varunits,timedt,color); % jut to overlay

% print([path_plots,'Spectra_wctd_raw_z1-80m.eps'],'-depsc2');

% Overlay Buoyancy frequency in Spectra plot
legend('AutoUpdate','off')

Nmin=nanmin(N(:)); % minimum buoyancy freq (s^(-1))
Nmax=nanmax(N(:)); % maximum buoyancy freq (s^(-1))
Nmean=nanmean(N(:));% maximum buoyancy freq (s^(-1))

% Plot Mean N vertical Bar
xlocation=Nmean;
ylocation=10^(0);

err_high=1;
err_low=10^(-12);

semilogy([xlocation xlocation],[err_low err_high].*ylocation,'color',lgrey,'linewidth',4)
xlocation=Nmax;
semilogy([xlocation xlocation],[err_low err_high].*ylocation,'color',lgrey,'linewidth',4)
xlocation=Nmin;
semilogy([xlocation xlocation],[err_low err_high].*ylocation,'color',lgrey,'linewidth',4)

ylim([err_low 1]);

ylim([10^(-10) 1]);
xlim([fTn 1]);

% print([path_plots,'Spectra_wctd_2mbinnedCTD_MRp_z1-80m_',method,'_wfNrange.eps'],'-depsc2');
print([path_plots,'Spectra_wctd_2mbinnedCTD_MRp_',num2str(z1),'-',num2str(z2),'m_',method,'_wfNrange.eps'],'-depsc2');


%%  Pre-processing of W profile data for epsilon estimates
%% ------------------------------------------------------------------------

% indzf=35;
% 
% data_truncated=W_2mb(1:indzf,:); % Define W_raw data to use (3 different signals)
% p_truncated=D_2mb(1:indzf,:); % Define W_raw data to use (3 different signals)

% data_truncated=W_2mb(1:zf,:); % Define W_raw data to use (3 different signals)
% p_truncated=D_2mb(1:zf,:); % Define W_raw data to use (3 different signals)

% data=W_2mb(:,:); % Define W_raw data to use (3 different signals)
% data_p=D_2mb(:,:); % will truncate after filtering, etc
% sensor='CTD2mbin';

% data=W_2mb_interp(:,:); % Define W_raw data to use (3 different signals)
% data_p=p2mb_interp(:,:); % will truncate after filtering, etc
% sensor='CTD2mbin_interp';

% data=W_1hz_interp(:,:); % Define W_raw data to use (3 different signals)
% data_p=p1hz_interp(:,:); % will truncate after filtering, etc
% sensor='CTD1hz_interp';

data=W_1hzlp_interp(:,:); % Define W_raw data to use (3 different signals)
data_p=p1hzlp_interp(:,:); % will truncate after filtering, etc
sensor='CTD1hzlp_interp';

% data=W_1hz(:,:); % Define W_raw data to use (3 different signals)
% data_p=P_1hz(:,:); % will truncate after filtering, etc
% sensor='CTD1Hz';

% data=W_mr_interp(:,:); % Define W_raw data to use (3 different signals)
% data_p=pmr_interp(:,:); % will truncate after filtering, etc
% sensor='MRp';

% Preallocate memory for processed profiles
[nz,np]=size(data);
data_hp=nan(nz,np);
rms_w_hp=nan(nz,np);
e_1=nan(nz,np);
ce=1;

Loz=nan(nz,np);
Kz_1=nan(nz,np);

%%  Following processing steps for Spray glider data

dt=10; % sampling period for 2m bin data
fs=1/dt; % sampling frequency
T=640; % floor(nanmean(nansum(Nbin))); % approx length of record in sec

fc_lpass=0.04; % 0.02when noise level starts to increase in spectra for w
fc_hpass=1/T; % no high pass if using legth of record
fc_hpass=1/150; % ~30m @ a PV~0.2

% Defining filters
Fd_hp=designfilt('highpassiir', 'FilterOrder', 4, 'PassbandFrequency', fc_hpass, 'PassbandRipple', 1, 'SampleRate', fs);
Fd_lp=designfilt('lowpassiir', 'FilterOrder', 4, 'PassbandFrequency', fc_lpass, 'PassbandRipple', 1, 'SampleRate', fs);

% data=nan(nz,np);
[nz,np]=size(data);
data_filt=nan(nz,np);
data_filt_lp=nan(nz,np);
data_filt_lp_hp=nan(nz,np);
rms_w_hp=nan(nz,np);

W_detrend=cell(3,1);
w_meas_dtrend=nan(nz,np);

d=1;
for p=1:np
    
%         % Remove mean/trend, to get w around a ~zero mean
%         y=detrend(data_truncated(:,p));
%         data(:,p)=detrend(y); % detrend
    
        indnani=find(~isnan(data(:,p)),1,'first');
        indnanf=find(~isnan(data(:,p)),1,'last');
        
        % Quadratic detrending
        x=data_p(indnani:indnanf,1); %,p);
        y=data(indnani:indnanf,p);
        
%        if length(y)>15
            
            % % Quadratic polynomial fit
            [p_coef, S, mu] = polyfit(x,y,4);
            pfit(indnani:indnanf,p)=polyval(p_coef,x,[],mu);
            
            w_meas_dtrend(indnani:indnanf,p)=y-pfit(indnani:indnanf,p); % detrend low-frequencies introduced by buoyancy control
            
            W_detrend{d}=w_meas_dtrend;
%         else
%             continue
%         end
    
    
    % Apply additional filters
    y=w_meas_dtrend(indnani:indnanf,p);
%    data_filt(:,p)=filtfilt(Fd_hp,y); % to remove low freq signals (ommiting)
    
    data_filt_lp(indnani:indnanf,p)=filtfilt(Fd_lp,y); % to remove high freq sensor noise
    data_filt_lp_hp(indnani:indnanf,p)=filtfilt(Fd_hp,data_filt_lp(indnani:indnanf,p)); % to remove low freq internal wave (non dissipative)
 
    
    % Compute rms_velocity;
%     y=data(indnani:indnanf,p); % Signal to use for epsilon estimate,
%     method='raw'; 

%     y=data(indnani:indnanf,p)-nanmean(data(indnani:indnanf,p)); % Signal to use for epsilon estimate,
%     method='detrend'; 

%     y=w_meas_dtrend(indnani:indnanf,p); % Signal to use for epsilon estimate,
%     method='detrend_polyfit'; 
    
%     y=data_filt_lp(indnani:indnanf,p); % Signal to use for epsilon estimate,
%     method=['lp_fc',num2str(fc_lpass)];
    
    y=data_filt_lp_hp(indnani:indnanf,p); % Signal to use for epsilon estimate,
    method=['lp_hp_fc',fc_hpass];
    
    k=11; % 11; % moving std over k points; k=2 (for 2mb data)
    k=6757; % ; % moving std over k points; note that for interpolated data 20m~6757
   %  k=3201; % ; % moving std over k points; note that for interpolated data 10m~3201
    k=find(data_p==10,1,'first'); % the size of the window should be representative of sevral eddy length scales

   % rms_w_hp(indnani:indnanf,p)=movvar(y,k); % in units of velocity squared (moving variance, makes sense when using raw data!)
    rms_w_hp(indnani:indnanf,p)=movmean(y.^2,k); % RMS!! of high-pass w

%     % Estimate Dissipation Rates
%     e_1(:,p)=N(:,p).*rms_w_hp(:,p); % where ce is a calibration coefficient
    
end

%% Plot Filtering
%%  Example plot: binned data HR

Z=rms_w_hp; Zname='w_rms'; Zunits='[m/s]';
% Z=SAint; Zname='Salinity'; Zunits='[g kg^{-1}]';
% Z=RHOint-1000; Zname='Sigma'; Zunits='[mg m^{-3}]';
% Z=sqrt(N2int); Zname='N'; Zunits='[s^{-1}]';
% Z=N2int; Zname='N2'; Zunits='[s^{-2}]';
LocalTime=Time-8/24;
[X,Y]=meshgrid(LocalTime,data_p(:,1));

close all

figure(1)
h=pcolor(X,Y,Z);
h.EdgeColor = 'none';
hold on
axis ij

hc=colorbar; 
set(get(hc,'Title'),'String',Zunits,'fontsize',20);
cmin=floor(nanmin(Z(:))); cmax=floor(nanmax(Z(:)));
cmin=nanmin(Z(:)); cmax=nanmax(Z(:));
caxis([cmin cmax])
if strcmp(Zname,'w_rms')
    cmap=flipud(cbrewer('div','RdBu',25));
end
colormap(cmap)

% % % Overlay Density Contours
%     hold on
%     v=24:0.02:27;
%     v1=[floor(nanmin(RHOint(:)-1000)):0.05:26.7];
%     v2=[26.7:0.05:ceil(nanmax(RHOint(:)-1000))];
%     v=[v1 v2];
%     [C,h]=contour(X,Y,RHOint-1000,v); %contourf
%     % v=floor(nanmin(nanmin(Temp))):0.4:ceil(nanmax(nanmax(Temp)));
%     % [C,h]=contour(Time,Depth,Temp,v); %contourf
%     set(h,'color','k');
%     clabel(C,h,v,'FontSize',18,'LabelSpacing',1200); %,'FontSize',16,'LabelSpacing',1200); %,'manual);
%     
% %   % Overlay MLD (largest threshold)
%     hold on
%     plot(X(1,:),MLD(:,3),'--o','color',[0.9 0.9 0.9],'linewidth',2,'markersize',10,...
%         'markerfacecolor',[0.9 0.9 0.9],'markeredgecolor',[0.2 0.2 0.2]);
%     plot(X(1,:),MLD_n2,'--o','color',[0.9 0.9 0.9],'linewidth',2,'markersize',10,...
%         'markerfacecolor',[0.9 0.9 0.9],'markeredgecolor',[0.2 0.2 0.2]);

ylabel('Depth [m]');
title(Zname)

% round date limits to the hour
tmin=nanmin(X(:)); tmax=nanmax(X(:));
dmin=datevec(tmin); dmax=datevec(tmax);
dmin=[dmin(1:4) 0 0]; dmax=[dmax(1:4) 0 0];
xtick([datenum(dmin):1/24:datenum(dmax)])
xlim([datenum(dmin) datenum(dmax)])
datetick('keeplimits','keepticks')
day=datetime(datevec(LocalTime),'Format','MM/dd');
xlabel([char(day(1)),' - ',char(day(end-1))]);
% xlabel('Longitude [^{\circ} W]');
set(gca,'xMinortick','on')
xlim([nanmin(X(:)) nanmax(X(:))])
% xticklabel_rotate([],45,[]);

print_setup(60,30)

caxis([-0.001 0.001])

%% 
% % %% Estimate epsilon from the Large Eddy Method (LEM)
% % % e=ce*N*rms_w_hp; where ce is a calibration coefficient
% % 
% % % Filter out non-dissipative linear internal waves
% % % e.g. with abutterworth filter and cutoff wavenumber of ~ lambda=30m or
% % % freq of ~300s (~40 data points, ~8s and ~1m per bin)
% % 
% % Fd=designfilt('highpassiir', 'FilterOrder', 4, 'PassbandFrequency', 0.003, 'PassbandRipple', 0.1, 'SampleRate', 0.125);
% % 
% % w_hp=nan(nZ-1,nP);
% % for p=1:nP
% %     
% %     indnani=find(~isnan(w_meas_filt(:,p)),1,'first');
% %     indnanf=find(~isnan(w_meas_filt(:,p)),1,'last');
% %     
% %     % Signal to filter
% %     y=w_meas_filt(indnani:indnanf,p);
% %     
% %     w_hp(indnani:indnanf,p)=filtfilt(Fd,w_meas_filt(indnani:indnanf,p)); %low-pass to remove high-freq noise
% % 
% % end


%%  Cal coefficient

e_03=nan(size(W_2mb));
e_1=nan(size(W_2mb));

e_03=nan(size(rms_w_hp));
e1=nan(size(rms_w_hp));

% ce=0.37; % From Beaird et al. (2012), from straight comparison against VMP (stratified ocean)
% ce=1; % Garrett (1988), from straight comparison against VMP.. first application of Eq to the ocean (tidal channel)
%     % % ce=3-5; % in both convective and stratified atm boundary layer
%     % e_03=ce.*N.*rms_w_hp; % where ce is a calibration coefficient
%     % ce=1; % Garrett (1988), from straight comparison against VMP.. first application of Eq to the ocean (tidal channel)
%     % e_1=ce.*N.*rms_w_hp; % where ce is a calibration coefficient
%     % 
%     % e=e_1; cen='ce1';

% e_1=ce.*N.*rms_w_hp; 

e_1=ce.*N_interp.*rms_w_hp; 

%% Bin MR epsilon estimates in 2-mbins

zres=2;
zi=1; zf=100;
edges=[zi:zres:zf];
p2mb=edges-zres/2;

[val,indb]=histc(mr_paxis,edges);

Nb=unique(indb(:));
epsilon2_2m=nan(length(Nb),np);

indzi=searchclosest_m(mr_paxis,zi); % to use only from 10-70m
indzf=searchclosest_m(mr_paxis,zf);

for p=1:f

    [val,indb]=histc(mr_paxis(indzi(p):indzf(p),p),edges);
     Nb=unique(indb(:));
    
    for b=1:length(Nb)
        ind=find(indb==Nb(b));
        epsilon2_2m(b,p)=nanmedian(epsilon2(ind,p));
    end
    
end


%% Bin e from LEM estimates in 2-mbins
% 
% zres=2;
% zi=1; zf=100;
% edges=[zi:zres:zf];
% p2mb=edges-zres/2;
% 
% [val,indb]=histc(D_2mb,edges);
% 
% Nb=unique(indb(:));
% e1_2m=nan(length(Nb),np);
% 
% indzi=searchclosest_m(D_2mb,zi); % to use only from 10-70m
% indzf=searchclosest_m(D_2mb,zf);
% 
% for p=1:f
% 
%     [val,indb]=histc(D_2mb(indzi(p):indzf(p),p),edges);
%      Nb=unique(indb(:));
%     
%     for b=1:length(Nb)
%         ind=find(indb==Nb(b));
%         e1_2m(b,p)=nanmedian(e_1(ind,p));
%     end
%     
% end

%% Bin e from LEM estimates in 2-mbins, if uding interpolated data

zres=2;
zi=1; zf=100;
edges=[zi:zres:zf];
p2mb=edges-zres/2;

[val,indb]=histc(data_p,edges);

Nb=unique(indb(:));
e1_2m=nan(length(Nb),np);

indzi=searchclosest_m(data_p,zi); % to use only from 10-70m
indzf=searchclosest_m(data_p,zf);

for p=1:f

   %[val,indb]=histc(data_p(indzi(p):indzf(p),p),edges);
    [val,indb]=histc(data_p(indzi:indzf,1),edges); %for interpolated data
     Nb=unique(indb(:));
    
    for b=1:length(Nb)
        ind=find(indb==Nb(b));
        e1_2m(b,p)=nanmedian(e_1(ind,p));
    end
    
end

%% Regression: estimated epsilon (e_1) vs epsilon2_2mb from the MR
close all

lgrey=[0.9 0.9 0.9];
dgrey=[0.4 0.4 0.4];

figure
for p=1:np
% regresion(e_1(10:end,p),epsilon2_2m(10:35,p),2);
% hold on
% scatter(e_1(:,p),epsilon2_2m(1:35,p),50,edges(:)+1,'filled');

subplot(121)

semilogx(epsilon2_2m(:,p),p2mb,'-o','color',lblue); axis ij
hold on
h1=semilogx(nanmedian(epsilon2_2m,2),p2mb,'-o','color',blue); axis ij

semilogx(e1_2m(:,p),p2mb,'--o','color',lred); axis ij
hold on
h2=semilogx(nanmedian(e1_2m,2),p2mb,'--o','color',red); axis ij
ylim([1 100])
grid on

legend([h1,h2],'\epsilon_{MR}','e_{LEM}','location','NW');

subplot(122)
semilogx(epsilon2_2m(:,p)./e1_2m(:,p),p2mb,'-','color',dgrey); axis ij
hold on
ylim([1 100])
grid on

print_setup(50,40)

% print([path_plots,'Epsilon_mr_vs_wctd_2mbinCTD_z1-70m_',method,'_p',num2str(p),'.eps'],'-depsc2');

end

% % print([path_plots,'Epsilon_mr_vs_wctd_2mbinCTD_z1-70m_',method,'_all.eps'],'-depsc2');
% print([path_plots,'Epsilon_mr_vs_wctd_2mbinCTD_z',num2str(zi),'-',num2str(zf),'m_',method,'_all.eps'],'-depsc2');

print([path_plots,'Epsilon_mr_vs_wctd_',sensor,'_z',num2str(zi),'-',num2str(zf),'m_',method,'_all.eps'],'-depsc2');

%% Overlay Average profiles

figure
subplot(121)
% h1=semilogx(nanmedian(epsilon2_2m(:,:),2),p2mb,'-o'); axis ij
% hold on
% h2=semilogx(nanmedian(e1_2m(:,:),2),p2mb,'-o'); axis ij
h1=semilogx(nanmean(epsilon2_2m(:,:),2),p2mb,'-o'); axis ij
hold on
h2=semilogx(nanmean(e1_2m(:,:),2),p2mb,'-o'); axis ij
hold on
grid on
ylim([1 100])

% legend(h1,h2,'\epsilon_{MR}','e_{LEM}','location','SE');
legend('\epsilon_{MR}','e_{LEM}','location','NW');

subplot(122)
% mCe=nanmedian(epsilon2_2m(:,:),2)./nanmedian(e1_2m(:,:),2);
% mCe=nanmedian((epsilon2_2m./e1_2m),2);% new estimate, first calculate Ce, then average
mCe=nanmean(epsilon2_2m(:,:),2)./nanmean(e1_2m(:,:),2);
semilogx(mCe,p2mb,'ko'); axis ij
hold on
% semilogx(mCe-1,D_2mb(indzi:indzf,p),'--g'); hold on
% semilogx(mCe+1,D_2mb(indzi:indzf,p),'--m'); axis ij
ylim([1 100])
grid on

xlabel('C_e = \epsilon_{MR}/e_{LEM}')

print_setup(50,40)

% % print([path_plots,'MeanEpsilon_mr_vs_wctd_2mbinCTD_z1-70m_',method,'.eps'],'-depsc2');
% print([path_plots,'MeanEpsilon_mr_vs_wctd_2mbinCTD_z1-100m_',method,'.jpg'],'-djpeg');

%print([path_plots,'MeanEpsilon_mr_vs_wctd_2mbinCTD_z',num2str(zi),'-',num2str(zf),'m_',method,'.jpg'],'-djpeg');
print([path_plots,'MeanEpsilon_mr_vs_wctd_',sensor,'_z',num2str(zi),'-',num2str(zf),'m_',method,'.jpg'],'-djpeg');

%% Overlay Average profiles (linear scale)

figure
subplot(121)
% h1=plot(nanmedian(epsilon2_2m(:,:),2),p2mb,'-o'); axis ij
% hold on
% h2=plot(nanmedian(e1_2m(:,:),2),p2mb,'-o'); axis ij

h1=plot(nanmean(epsilon2_2m(:,:),2),p2mb,'-o'); axis ij
hold on
h2=plot(nanmean(e1_2m(:,:),2),p2mb,'-o'); axis ij

hold on
grid on
% ylim([1 100])

% legend(h1,h2,'\epsilon_{MR}','e_{LEM}','location','SE');
legend('\epsilon_{MR}','e_{LEM}','location','NW');

title(sensor)

subplot(122)
% mCe=nanmedian(epsilon2_2m(:,:),2)./nanmedian(e1_2m(:,:),2);
% mCe=nanmedian((epsilon2_2m./e1_2m),2);% new estimate, first calculate Ce, then average
mCe=nanmean(epsilon2_2m(:,:),2)./nanmean(e1_2m(:,:),2);


plot(mCe,p2mb,'ko'); axis ij
hold on
% semilogx(mCe-2,D_2mb(indzi:indzf,p),'--g'); %axis ij
% semilogx(mCe+2,D_2mb(indzi:indzf,p),'--m'); axis ij
% ylim([1 100])
grid on

xlim([0 1])

print_setup(50,40)

% % % print([path_plots,'MeanEpsilon_mr_vs_wctd_2mbinCTD_z1-70m_',method,'.eps'],'-depsc2');
% % print([path_plots,'MeanEpsilon_mr_vs_wctd_2mbinCTD_z1-100m_',method,'.jpg'],'-djpeg');
% print([path_plots,'MeanEpsilon_mr_vs_wctd_2mbinCTD_z',num2str(zi),'-',num2str(zf),'m_',method,'_linear.jpg'],'-djpeg');

print([path_plots,'MeanEpsilon_mr_vs_wctd_',sensor,'_z',num2str(zi),'-',num2str(zf),'m_',method,'_linear.jpg'],'-djpeg');

%% Statistics for Ce

nanmean(mCe)


%% lagged correlation?

close all
x=nanmedian(e1_2m(:,:),2);
y=nanmedian(epsilon2_2m(:,:),2);

maxlag=10;
[c,lags] = xcorr(x,y,maxlag,'normalized');
stem(lags,c)

xlim([-maxlag maxlag])
xlabel('Lags (1 lag = 2m bin)')

ylabel('Correlation')

print_setup(40,30)
% print([path_plots,'LaggedCorr_e1_vs_eosilon_2mbin_z1-100m_',method,'.jpg'],'-djpeg');


%%

close all
for p=1:np
% regresion(e_1(10:end,p),epsilon2_2m(10:35,p),0);
% hold on

loglog(epsilon2_2m(:,p),e1_2m(:,p),'o');
hold on
grid on
scatter(log10(epsilon2_2m(:,p)),log10(e1_2m(:,p)),50,'filled');
hold on

[m(p)]=regresion(log10(epsilon2_2m(:,p)),log10(e1_2m(:,p)),0);

end

% regresion(log10(e_1(10:end,p)),log10(epsilon2_2m(10:35,p)),0);

close all
pdffit6(m,0.025,'k',0.6);
xlabel('C_e, calibration coefficient')
% print([path_plots,'Ce_wctd_2mbinnedCTD_z1-70m_',method,'.eps'],'-depsc2');

%% Regression in loglog scale
ce=nanmean(m); % 0.6405

close all
for p=1:np
% regresion(e_1(10:end,p),epsilon2_2m(10:35,p),0);
% hold on

loglog(epsilon2_2m(:,p),e1_2m(:,p),'o');
hold on
grid on
scatter(log10(epsilon2_2m(:,p)),log10(e1_2m(:,p)),50,'filled');
hold on

[m(p)]=regresion(log10(epsilon2_2m(:,p)),log10(e1_2m(:,p)),0);

end

%% Regrssion prfile by profile
close all

for p=1:np
s(p)=regresion(e1_2m(10:end-5,p),epsilon2_2m(10:end-5,p),0);
scatter(e1_2m(:,p),epsilon2_2m(:,p),50,'filled');
hold on
end

%%
ms=regresion(nanmedian(e1_2m,2),nanmedian(epsilon2_2m,2),0);

%%
close all 
pdffit2(s,0.1,'k')

%%
close all
pdffit2(log10(epsilon2_2m(:)),0.5,'k')

hold on
pdffit2(10^(0.645).*log10(e_1(:)),0.5,'r')
% %% semilogyfit
% semilogy(x,y,'o','color',[0.2 0.2 0.2],'markerfacecolor',[0.2 0.2 0.2])
% hold on
% grid on
% 
% P = polyfit(x, log(y), 1);
% fit = exp(polyval(P, x));
% plot(x,fit,'k','linewidth',2)
% 
% a=P(1); b=P(2);
% [r,pval]=corr(x,log(y));
% 
% label = ['y = ',num2str(exp(b), '%6.2f'),' exp^{',num2str(a, '%6.3f'),'x} , r =',...
%     num2str(r, '%6.3f'),' , p =',num2str(pval, '%6.3f')];
% lgnd=legend(label,'location','NW');
% set(lgnd,'box', 'off','fontsize',20);



%%

coeff1=e_1./epsilon2_2m(1:35,:);

coeff=nanmedian(e_1,2)./nanmedian(epsilon2_2m(1:35,:),2);


figure(3)
pdffit4(coeff(:),'k',0.6)

% % %%% From Glider processing code:
% % % Fd=designfilt('highpassiir', 'PassbandFrequency', 0.011, 'StopbandFrequency', 0.009, 'PassbandRipple', 1, 'StopbandAttenuation', 60, 'SampleRate', 0.125);
% % Fd_bp=designfilt('bandpassiir', 'FilterOrder', 4, 'PassbandFrequency1', 1/5400, 'PassbandFrequency2', 1/100, 'PassbandRipple', 0.1, 'SampleRate', 0.125);
% % 
% % Fd_lp=designfilt('lowpassiir', 'FilterOrder', 4, 'PassbandFrequency', 0.01,'PassbandRipple', 0.1, 'SampleRate', 0.125); %low-pass filter to remove high-freq noise <100s
% % Fd_hp=designfilt('highpassiir', 'FilterOrder', 4, 'PassbandFrequency', 1/5400,'PassbandRipple', 0.1, 'SampleRate', 0.125); %low-pass filter to remove low freq from inacuracies of the model 5400s cutoff
% % 
% % pfit=nan(nZ-1,nP);
% % w_meas_dtrend=nan(nZ-1,nP);
% % w_meas_lp=nan(nZ-1,nP);
% % w_meas_filt=nan(nZ-1,nP);
% % w_water_hpfilt=nan(nZ-1,nP);
% % 
% % for p=1:nP
% %     
% %     indnani=find(~isnan(w_meas_interp(:,p)),1,'first');
% %     indnanf=find(~isnan(w_meas_interp(:,p)),1,'last');
% %     
% %     % Quadratic detrending
% %     x=du.depth(indnani:indnanf,p);
% %     y=w_meas_interp(indnani:indnanf,p);
% %     
% %     if length(y)>15
% %      
% % %     % Quadratic polynomial fit
% %     [p_coef, S, mu] = polyfit(x,y,4);
% %     pfit(indnani:indnanf,p)=polyval(p_coef,x,[],mu);
% %     
% %     w_meas_dtrend(indnani:indnanf,p)=y-pfit(indnani:indnanf,p); % detrend low-frequencies introdced by buoyancy control
% %     
% %     % High Pass (to remove high freq noise >1/100s)
% %      w_meas_lp(indnani:indnanf,p)=filtfilt(Fd_lp,w_meas_dtrend(indnani:indnanf,p)); %low-pass to remove high-freq noise (T>100s)
% %      w_meas_filt(indnani:indnanf,p)=filtfilt(Fd_hp,w_meas_lp(indnani:indnanf,p)); % high-pass with cutoff 1/5400s (~1.5h), after low-pass was applied
% %      
% % %     % Band-Pass filter (1/100s < f < 1/5400);
% % %      w_water_hpfilt(indnani:indnanf,p)=filtfilt(Fd_bp,w_meas_dtrend(indnani:indnanf,p));
% % %      % agrees with above , except on edges but above method agrees better
% % %      with the model based w
% % 
% % 
% % %     % 4) Finally, data above 125 m and within 75m of the maximum dive depth are excluded to remove accelerations/decelarations not due to water motion 
% % %     % and to avoid filter edge effects
% % %     w_meas_filt(indnani:indnani+125,p)=NaN;
% % %     w_meas_filt(indnanf-25:indnanf,p)=NaN;
% %     
% %     else
% %         continue
% %     end
% % end
% % 
% % % % Estimated error from R2013, w_water=0.005 m/s
% % % % should compute spectra to estimate the noise level
% % % w_error=zeros(nZ-1,nP);
% % % w_error(w_meas_dtrend<=0.002)=1;
