% %% Code to load in glider data and calculate MLD averages
% 
% % load in data
% % Glider SN209
% load('\\atlas.shore.mbari.org\ProjectLibrary\901805_Coastal_Biogeochemical_Sensing\Spray_Data\25520301\25520301sat.mat');

function avgs = mldAverages(s)
    %% create vars
    % make a dive direction variable
    dnms = length(s.sdn) ;
    s.divedir = floor(diff(s.pres)) ;
    s.divedir(s.divedir>1) = 1;
    s.divedir(s.divedir<1) = -1;
    s.divedir = s.divedir .* -1; % change down to -1 and up to 1 because my brain works that way
    
    % choose MLD threshold:
    dens_thresh = 0.03; % kg/m3
    
    % preallocate variables
    s.mld = nan(size(s.sdn));
    
    % struct named 'sn209' is what I funnel into the ODSS code
    sn209.time = reshape(s.sdn_,1,[]);
    sn209.lat = reshape(s.lat_,1,[]);
    sn209.lon = reshape(s.lon_,1,[]);
    sn209.mld = [];
    sn209.ph = [];
    sn209.rhodamine = [];
    
    for i = 1:dnms
        % pull out dive and ascent profiles
        down = s.divedir(:,i) == -1 ;
        up = s.divedir(:,i) == 1 ;
    
        % for dive
        % calculate MLD
            pressure_prof = s.pres(down,i); % extract pressure profile
    
            % uncomment below for 10 dbar reference
            % [~,ref_idx] = min(abs(pressure_prof-10));
    
            % surface / 5 dbar reference
            [~,ref_idx] = min(abs(pressure_prof-5));
    
            density_prof = s.pdens(down,i); % extract density profile
    
            % define reference density
            density_ref = density_prof(ref_idx);
            MLDd_idx = find(s.pdens(down,i) > density_ref+dens_thresh, 1);
    
            if isempty(MLDd_idx) == 1
                MLDd(i) = NaN ;
            else
                MLDd(i) = pressure_prof(MLDd_idx) ; % dive MLD
            end
    
            % calculate pH average for dive
            phd = (s.pHin(down,i));
            phd_mean = nanmean(phd(1:MLDd_idx)) ;
            
            % Calculate rhodamine average for dive
      % +++ need to change the variable, it's a random input (i490) right now +++ %      
            rd = (s.i490(down,i));
            rd_mean = nanmean(rd(1:MLDd_idx)) ;
    
        % for ascent
        % calculate MLD
            pressure_prof = s.pres(up,i); % extract pressure profile
    
            % uncomment below for 10 dbar reference
            % [~,ref_idx] = min(abs(pressure_prof-10));
    
            % surface / 5 dbar reference
            [~,ref_idx] = min(abs(pressure_prof-5));
    
            density_prof = s.pdens(up,i); % extract density profile
            % define reference density
            density_ref = density_prof(ref_idx);
    
            MLDa_idx = find(flipud(s.pdens(up,i)) > density_ref+dens_thresh, 1);
            if isempty(MLDa_idx) == 1
                MLDa(i) = NaN ;
            else
                flp = flipud(pressure_prof);
                MLDa(i) = flp(MLDa_idx) ; % ascent MLD
            end
    
            % calculate pH average for ascent
            pha = flipud(s.pHin(up,i));
            pha_mean = nanmean(pha(1:MLDa_idx)) ;
            
            % calculate rhodamine average for ascent
      % +++ need to change the variable, it's a random input (i490) right now +++ %      
            ra = (s.i490(up,i));
            ra_mean = nanmean(ra(1:MLDa_idx)) ;
    
        % calculate ML means
            sn209.ph = [sn209.ph pha_mean] ;
            sn209.rhodamine = [sn209.rhodamine ra_mean] ;
    
        % save MLD 
            s.mld(i) = nanmean([MLDd(i) MLDa(i)]) ;
            sn209.mld = [sn209.mld MLDd(i) MLDa(i)] ;
            avgs = [sn209.mld MLDa(i)] ;
    end
end