function StepLim_values = StepLimFunction(x_val, number_of_cusp_points,...
    init_final_StepLim_values, xLimit, StepLim_functional_form)
    % Output - 
    %    StepLim_values - StepLim given the x-coordinate, basal StepLim specified,
    %                     if there is a viscosity gradient
    % 
    % Inputs - 
    %    x_val - x-coordinate of particle
    %    number_of_cusp_points - Number of cusps in case of a step-wise
    %    viscosity gradient, only matters if StepLim_functional_form is set
    %    as 'stepwise'
    %    init_final_StepLim_values - 1x2 array with values for StepLim at left and right
    %                                ends of the box.
    %    xLimit - box length 
    % StepLim_functional_form - 'default' for no StepLim gradient, 'stepwise'
    %                           for a stepwise StepLim gradient, 'linear' for a 
    %                           linear StepLim gradient 

    if strcmpi(StepLim_functional_form,'default') && ...
            init_final_StepLim_values(1) ~= init_final_StepLim_values(2)
        error('Error:incorrectInput',...
            'Error. \nStepLim functional form should be stepwise or linear and not default, since the StepLim varies along the box length')
    end
    StepLim_values = zeros(size(x_val));
    if strcmpi(StepLim_functional_form,'stepwise')
        cusps = 0:xLimit/(number_of_cusp_points+1):xLimit; % 5
        StepLim_possibilities = linspace(init_final_StepLim_values(1),...
            init_final_StepLim_values(2),number_of_cusp_points+1);
    end
    
    for i = 1:length(x_val)
        if strcmpi(StepLim_functional_form,'stepwise')
            for j = 1:length(cusps)
                if x_val(i) < cusps(end)
                    if x_val(i) >= cusps(j) && x_val(i) < cusps(j+1)
                        StepLim_values(i) = StepLim_possibilities(j);
                    end
                else
                    StepLim_values(i) = StepLim_possibilities(end);
                end
            end
        elseif strcmpi(StepLim_functional_form,'linear')
            StepLim_values(i) = (((init_final_StepLim_values(2) - init_final_StepLim_values(1))...
                /(xLimit - 0))*(x_val(i)-0)) + init_final_StepLim_values(1);
        end
    end
end