function profile = tm_dissipation_profile(dudz,Pr,w,T,Ax,Ay,Az,fs,binsize,...
    inverse_ovlap,numensembles,spectra_inverse_ovlap,name,do_plot)

% profile = tm_dissipation_profile(dudz,Pr,w,sbeT,Ax,Ay,Az,fs,binsize,...
%            inverse_ovlap,numensembles,spectra_inverse_ovlap,name,do_plot)
%
% modified from vmp_dissipation_profile.m for Turbomap data (KJR 26 July 2011)
% included default spectra_inverse_ovlap

%% defaults
if ~exist('binsize','var') || isempty(binsize)
    binsize = 4;
end
if ~exist('inverse_ovlap','var') || isempty(inverse_ovlap)
    inverse_ovlap = 2;
end
if ~exist('numensembles','var') || isempty(numensembles)
    numensembles = 8;
end
if ~exist('spectra_inverse_ovlap','var') || isempty(spectra_inverse_ovlap)
    spectra_inverse_ovlap= 2;
end
if ~exist('name','var') || isempty(name)
    str = 'Calculating epsilon profile';
    name = 'N/A';
else
    str = ['Calculating epsilon profile (' name ')'];
end
if ~exist('do_plot','var') || isempty(do_plot)
    do_plot = 0;
end

%% Figure out overlap and array sizes
if inverse_ovlap == 0
    bin_offset = binsize;
    bo = binsize;
else
    bin_offset = binsize/inverse_ovlap;
    bo = bin_offset;
end

pmin = min(Pr);
if mod(pmin,bo) ~= 0
    pmin = pmin + bin_offset-mod(pmin,bo);
end
clear bo
    
pmax = floor(max(Pr));
total_depth = pmax - pmin;

if inverse_ovlap == 0
    N = fix(total_depth/binsize);
else
    N = fix(total_depth/bin_offset) - (inverse_ovlap - 1);
end

%% Calculated the NFFT for the spectra, and for storing the spectra
%% try and account for spikes any spikes in w
N1 = binsize/max([min(w) mean(w)-3*std(w)])*fs;
if spectra_inverse_ovlap == 0
    ovlap = 0;
else
    ovlap = 1/spectra_inverse_ovlap;
end
M = floor(N1/(numensembles - ovlap*numensembles + ovlap));
NFFT = 2^nextpow2(M);
clear N1 M ovlap

%% Gravity
g = 9.81;

%% Initialize arrays 
profile.epsilon = nan*ones(N,1);
profile.z_db = nan*ones(N,1);
profile.nu = nan*ones(N,1);
profile.kmin = nan*ones(N,1);
profile.kc = nan*ones(N,1);
profile.kolmogorov = nan*ones(N,1);
profile.rms = nan*ones(N,1);
profile.used_unclean = nan*ones(N,1);
profile.k = nan*ones(N,NFFT/2+1);
profile.spectra = nan*ones(N,NFFT/2+1);

%% Main loop
counter = 0;
if ~do_plot
    h_wb=waitbar(0,str);
end
for pres = pmin+bin_offset:bin_offset:pmin+(N-1)*bin_offset
    
    counter = counter + 1;
    
    start_pr = pres - bin_offset;
    end_pr = pres + bin_offset;
    idx_f = find(Pr >= start_pr & Pr < end_pr);



    
    W = nanmean(w(idx_f)); 
    nu = visc35(nanmean(T(idx_f)));
    dudz_int=dudz(idx_f);  
    Ax_int=Ax(idx_f);  
    Ay_int=Ay(idx_f);
    Az_int=Az(idx_f);
        
    [epsilon,details] = vmp_epsilon(dudz_int,g*Ax_int,g*Ay_int,g*Az_int, ...
        W,nu,fs,numensembles,spectra_inverse_ovlap,NFFT);
    
    profile.epsilon(counter) = epsilon;
    profile.z_db(counter) = pres;
    profile.nu(counter) = nu;
    profile.kmin(counter) = details.kmin;
    profile.kc(counter) = details.kc;
    profile.kolmogorov(counter) = details.kolmogorov;
    profile.rms(counter) = details.rms;
    profile.used_unclean(counter) = details.use_unclean;
    profile.k(counter,:) = details.k(:)';
    profile.spectra(counter,:) = details.spectra(:)';
    
    if do_plot
        clf
        Nasmyth_values = [1e-11 1e-10 1e-9 1e-8 1e-7 1e-6 1e-5 1e-4];
        Nasmyth_spectra = nasmyth(Nasmyth_values, nu, details.k);
        xx = [details.kmin details.kc details.kc details.kmin details.kmin];
        yy = [10^-15 10^-15 10^-1 10^-1 10^-15]; 
        Pn = nasmyth(epsilon,nu,details.k);
        
        subplot(2,1,1)
        plot(idx_f,dudz(idx_f),'b')
        title(sprintf('depth: %d to %d, \\epsilon: %e, rms: %e',start_pr,end_pr,epsilon,details.rms))
        subplot(2,1,2)
        loglog(details.k,details.spectra,'r')
        hold on
        patch(xx,yy,'b','facecolor',0.8*[1 1 1],'edgecolor','none')
        loglog(details.k,Nasmyth_spectra,'g')
        loglog(details.k,details.spectra,'r',details.k,Pn,'b')
        axis tight
        pause
    else
        waitbar(counter/(N-1),h_wb);
    end
end
close(h_wb);   


%% Attributes
profile.epsilon_att(1).Name = 'Description';
profile.epsilon_att(1).Value = sprintf('Turbulent kinetic energy dissipation from shear probe %s',name);
profile.epsilon_att(2).Name = 'units';
profile.epsilon_att(2).Value = 'W/kg';

profile.z_db_att(1).Name = 'Description';
profile.z_db_att(1).Value = sprintf('Z for epsilon calculation (binsize=%4.2f, offset=%4.2f)',binsize,bin_offset);
profile.z_db_att(2).Name = 'units';
profile.z_db_att(2).Value = 'dB';

profile.nu_att(1).Name = 'Description';
profile.nu_att(1).Value = 'Viscosity';
profile.nu_att(2).Name = 'units';
profile.nu_att(2).Value = 'm^2/s';

profile.kmin_att(1).Name = 'Description';
profile.kmin_att(1).Value = sprintf('Minimum wavenumber used in the integration for shear probe %s',name);
profile.kmin_att(2).Name = 'units';
profile.kmin_att(2).Value = 'cpm';

profile.kc_att(1).Name = 'Description';
profile.kc_att(1).Value = sprintf('Cutoff wavenumber used in the integration for shear probe %s',name);
profile.kc_att(2).Name = 'units';
profile.kc_att(2).Value = 'cpm';

profile.kol_att(1).Name = 'Description';
profile.kol_att(1).Value = sprintf('Kolmogorov wavenumber for epsilon calculated with shear probe %s',name);
profile.kol_att(2).Name = 'units';
profile.kol_att(2).Value = 'cpm';

profile.rms_att(1).Name = 'Description';
profile.rms_att(1).Value = 'RMS difference between log10 of data spectra and log10 of corresponding Nasmyth';
profile.rms_att(2).Name = 'units';
profile.rms_att(2).Value = 'None';

profile.used_unclean_att(1).Name = 'Description';
profile.used_unclean_att(1).Value = 'Flag that the acceleration corrected spectra could not be used';
profile.used_unclean_att(2).Name = 'units';
profile.used_unclean_att(2).Value = 'None';

profile.k_att(1).Name = 'Description';
profile.k_att(1).Value = 'Wavenumber vector for dudz spectra';
profile.k_att(2).Name = 'units';
profile.k_att(2).Value = 'cpm';

profile.spectra_att(1).Name = 'Description';
profile.spectra_att(1).Value = sprintf('dudz (%s) spectra',name);
profile.spectra_att(2).Name = 'units';
profile.spectra_att(2).Value = 's^2/cpm';
