function [y_smooth] = quad_smooth(x, y, d)

% [y_smooth] = quad_smooth(x, y, n)

temp = length(x);
ii = find(~isnan(y));
y = y(ii);
x = x(ii);
temp = temp - length(x);

% Make sure d is odd
d2 = floor(d./2);
x_start = x(1);
x = x-x_start;

% Initialize variables...

n = length(x);
x0 = 0 * x;
x1 = x0;
x2 = x0;
x3 = x0;
x4 = x0;
y1 = x0;
x1y1 = x0;
x2y1 = x0;
x00 = 1 + x0;

% Accumulate powers of x and y

for dd = -d2:d2
    % Set limits of index...
    index1a = 1;
    index1b = n;
    index2a = dd+1;
    if index2a < 1
        index2a = 1;
        index1b = n + dd;
    end
    index2b = n+dd;
    if index2b > n
        index2b = n;
        index1a = 1 + dd;
    end;
    
    % Accumulate powers of x and y...
    index2 = index2a:index2b;
    index1 = index1a:index1b;
   
    x0(index1) = x0(index1) + x00(index2);
    x1(index1) = x1(index1) + x(index2);
    x2(index1) = x2(index1) + x(index2).^2;
    x3(index1) = x3(index1) + x(index2).^3;
    x4(index1) = x4(index1) + x(index2).^4;
    
    y1(index1) = y1(index1) + y(index2);
    x1y1(index1) = x1y1(index1) + y(index2).*x(index2);
    x2y1(index1) = x2y1(index1) + y(index2).*x(index2).^2;
end

% Solve for quadratic fit
% -- First load matrix and vector
a1 = x0; a2 = x1; a3 = x2;
a4 = x1; a5 = x2; a6 = x3;
a7 = x2; a8 = x3; a9 = x4;
b1 = y1; b2 = x1y1; b3 = x2y1;

% clear x0 x1 x2 x3 x4 y1 x1y1 x2y1;

% -- Now solve for quadratic terms...
c1 = -((a5.*a7.*b1 - a4.*a8.*b1 - a3.*a7.*b2 + a2.*a8.*b2 + a3.*a4.*b3 - a2.*a5.*b3)./(-(a3.*a4.*a6) + a2.*a5.*a6 + a3.*a4.*a7 - a1.*a5.*a7 + a1.*a4.*a8 - a2.*a4.*a8));
c2 = -((-(a5.*a6.*b1) + a4.*a8.*b1 + a3.*a6.*b2 - a1.*a8.*b2 - a3.*a4.*b3 + a1.*a5.*b3)./(-(a3.*a4.*a6) + a2.*a5.*a6 + a3.*a4.*a7 - a1.*a5.*a7 + a1.*a4.*a8 - a2.*a4.*a8));
c3 = -((a4.*a6.*b1 - a4.*a7.*b1 - a2.*a6.*b2 + a1.*a7.*b2 - a1.*a4.*b3 + a2.*a4.*b3)./(-(a3.*a4.*a6) + a2.*a5.*a6 + a3.*a4.*a7 - a1.*a5.*a7 + a1.*a4.*a8 - a2.*a4.*a8));

% clear a1 a2 a3 a4 a5 a6 a7 a8 a9 b1 b2 b3;

% -- Finally, solve for smoothed variable from quadratic fit...
y_smooth = c1 + c2 .* x + c3 .* x.^2;
