function [x, e, n, g] = mesh_graded_line(c, l, g)
%MESH_GRADED_LINE Summary of this function goes here
%   Detailed explanation goes here

% Input:
%   c: 2 corner points
%   l: element size at corner points (smaller = stronger)
%   g: grading constant

% Output:
%   x: placed points
%   e: minimum and maximum element sizes
%   n: number of created elements

len = abs(diff(c));

if abs(diff(l)) < eps
    % Uniform spacing
    x = linspace(c(1), c(2), ceil(len/mean(l))+1).';
else
    % Non uniform spacing
    if l(1) > l(2)
        l = fliplr(l(:).');
        c = fliplr(c(:).');
        flipped = true;
    else
        flipped = false;
    end

    % Number of elements needed to reach grading
    n_grad = ceil(log(l(2) / l(1)) / log(g)) + 1;
    % Length achieved by exploiting full grading
    len_grad = l(1) * ((1 - g^n_grad) / (1 - g));

    % Check if full length can be exploited
    if len_grad < len
        % Divide the remainder length into elements
        x = cumsum([0; l(1)*g.^(0 : n_grad-2).']);
        % Add the remainder with equal spacing
        x = [x(1 : end-1); linspace(x(end), len, ceil((len-x(end))/(g*(x(end)-x(end-1))))+1).'];
    else
        % Divide the full length using the preferred grading
        n = log(1 - len*(1 - g)/l(1)) / log(g);
        % Create a logarithmic scale
        x = cumsum([0; l(1)*g.^(0 : ceil(n)-1).']);
        % Scale and add first element with predefined length
        if (isempty(x))
            x = [0; len];
        else
            x = x * len / x(end);
        end
    end

    % Map between end points
    x = (1-x/len)*c(1) + x/len*c(2);
    if flipped
        x = flipud(x);
    end
end
    
if nargout > 1
    e = [min(abs(diff(x))), max(abs(diff(x)))];
    if nargout > 2
        n = numel(x)-1;
    end
end

end % of function mesh_graded_line
