% This script uses the Sum of Squared Differences (SSD) method to
% compute the value of the curve-fitting parameter y0 in
% a fit of the form u'= u/vp = 1-exp(-y/y0).
%
% INPUT ARGUMENTS
% vp: flow velocity outside the boundary layer (scalar)
% u: flow velocities, calculated from pressures (row vector)
% y: y-values at which velocity was measured (row vector)

% VARIABLE CONVENTIONS used in this code
% u* (ustar): normalized velocity data (input)
% u' (uprime): velocity data calculated using y0

% The function finds the y0 that minimizes the sum of the squared
% difference between u* and u', sum((u*(i)-u'(i))^2)

function SSD_fit(vp,u,y)
ustar = u/vp;       % calculate the normalized velocities to be fitted 

% *** INPUT ***
                                                % values that work well
% initial guess for y0 (upper bound) 
y01U = input('Initial y0 (upper bound): ');     % 0.001 (two zeros)

% initial guess for y0 (lower bound)
y01L = input('Initial y0 (lower bound): ');     % 0.0001 (three zeros)

% amount to increment y0 by each step
step = input('y0 step size: ');                 % 0.000001 (five zeros)

% *** MAIN CALCULATION ***
L = length(y);          % length of input array, used in for loop and
                        % to create variable arrays used in calculation

% Call SSD fitting functions
y0L = SSDfitUp(ustar,y01L,L,step,y);        % returns lower bound for y0
y0U = SSDfitDown(ustar,y01U,L,step,y);      % returns upper bound for y0

% Calculate magnitude of error between upper and lower estimates
diff = y0L - y0U;
diff = abs(diff);

% *** OUTPUT ***
% Display results
if diff < 0.000001             % error allowed between upper and lower fits
    y0 = sprintf('Fitted y0: %d', y0L);
    disp(y0)
else 
    disp('The upper and lower fits do not match to within the given margin of error')
end

% *** FUNCTIONS USED IN CALCULATION ***

% SSDfitUp
% SSD fit iterating up from lower bound for y0
function y0L = SSDfitUp(ustar,y01L,L,step,y)
uprime = zeros(2,L);            % matrix to store computed u'
squareDiff = zeros(2,L);        % vector to store [(u*-u')^2]i

curDiff = 1;                    % current and next sum of squared
nextDiff = 0;                   % differences for iteration

while curDiff > nextDiff
    for i = 1:L
        % set "current" y0 to last iteration's "next" y0
        y0L = y01L;     
        % calculate "current" u' values
        uprime(1,i) = 1 - exp(-y(i)/y0L); 
    end
    for i = 1:L
        % increment y0 to the next value
        y01L = y0L + step;
        % calculate "next" u' values
        uprime(2,i) = 1 - exp(-y(i)/y01L);
    end
    for i = 1:L
        % vector of differences for u* - "current" u'
      	squareDiff(1,i) = (ustar(i) - uprime(1,i))^2;
    end
    for i = 1:L
        % vector of differences for u* - "next" u'
        squareDiff(2,i) = (ustar(i) - uprime(2,i))^2;
    end
        % sum of "current" squared differences
        curDiff = sum(squareDiff(1,:));
        % sum of "next" squared differences
        nextDiff = sum(squareDiff(2,:));
end

% SSDfitDown
% SSD fit iterating down from upper bound for y0
function y0U = SSDfitDown(ustar,y01U,L,step,y)
uprime = zeros(2,L);            % matrix to store computed u'
squareDiff = zeros(2,L);        % vector to store [(u*-u')^2]i

curDiff = 1;                    % current and next sum of squared
nextDiff = 0;                   % differences for iteration

while curDiff > nextDiff
    for i = 1:L
        % set "current" y0 to last iteration's "next" y0
        y0U = y01U;
        % calculate "current" u' values
        uprime(1,i) = 1 - exp(-y(i)/y0U); 
    end
    for i = 1:L
        % increment y0 to the next value
        y01U = y0U - step;
        % calculate "next" u' values
        uprime(2,i) = 1 - exp(-y(i)/y01U);
    end
    for i = 1:L
        % vector of differences for u* - "current" u'
      	squareDiff(1,i) = (ustar(i) - uprime(1,i))^2;
    end
    for i = 1:L
        % vector of differences for u* - "next" u'
        squareDiff(2,i) = (ustar(i) - uprime(2,i))^2;
    end
        % sum of "current" squared differences
        curDiff = sum(squareDiff(1,:));
        % sum of "next" squared differences
        nextDiff = sum(squareDiff(2,:));
end