function [yall, ypos] = f_localbndry(x,h,flag)

% function solves for Mahalanobis D^2 or Euclidean distance squared (E^2)
% (flag dependent) used to delineating local boundaries
%
% INPUT ARGUMENTS:
% z is the dependent variable of a traverse or time series
% h is a number denoting half the window size
% flag: 0 for D^2 or 1 for E^2
%
% OUTPUT:
% yall is either D^2 or E^2 for entire seriesdepending on vlaue of flag
% ypos is D^2 or E^2 calculated when the mean of the right side of the 
% window (X2mean as calculated below) is greater than the right side mean
% (X1mean)
%
% NOTES:
% The window size (2*h) is the entire number of points used in calculating
% means and standard deviations.  Calculations from each half of the window
% at every possible point are used in determining E^2 and D^2
%
% any flag other than 0 or 1 will result in an error

fill = NaN*ones(1,h);
switch flag
    case 0
        for i = 1:length(x)-2*h
            X1mean(i) = mean(x(i:i+h));
            X2mean(i) = mean(x(i+h:i+2*h));
            s1(i) = std(x(i:i+h));
            s2(i) = std(x(i+h:i+2*h));
            D2(i) = (X1mean(i) - X2mean(i))^2/(s1(i)+s2(i));
            if X2mean(i) >= X1mean(i)
                D2pos(i) = (X1mean(i) - X2mean(i))^2/(s1(i)+s2(i));
            else
                D2pos(i) = 0;
            end
        end
        yall = [fill,D2,fill];
        ypos = [fill,D2pos,fill];
        
    case 1
        for i = 1:length(x)-2*h
            X1mean(i) = mean(x(i:i+h));
            X2mean(i) = mean(x(i+h:i+2*h));
            E2(i) = (X1mean(i) - X2mean(i))^2;
            if X2mean(i) >= X1mean(i)
                E2pos(i) = (X1mean(i) - X2mean(i))^2;
            else
                E2pos(i) = 0;
            end
        end
        yall = [fill,E2,fill];
        ypos = [fill,E2pos,fill];
        
    otherwise
        display('improper flag used to call f_localbndry function')
end