function [DfTest, F, U, xline] = CalcUFourier(dfline, mpp, k, A, f0)
%[DfTest, F, U, xline] = CalcUFourier(dfline, mpp, k, A, f0)
%Takes a vector dfline and scalars mpp: meters per pixel; k: sensor
%stiffness; A: Oscillation amplitude; and f0: Centre frequency. 
%Returns DfTest, which should reproduce dfline; F and U which are the
%forces and energies, as well as xline, a corresponding vector.
%F and U might require a low-pass filter to avoid high-frequency ringing
%A.J. Weymouth, 2019

maxit = round(length(dfline)/2)-1;

xline = mpp*(1:length(dfline));
L=xline(end);

%Estimate a, b, alpha and beta
for n = 1:maxit
    InnerArg = ((2*n*pi)/L).*xline;
    alpha(n) = trapz(xline, sin(InnerArg).*dfline) / (L/2);
    an(n) = -alpha(n)*(2*k/f0)*(A*L/4/pi/n)/besselj(1, 2*pi*n*A/L);
    beta(n) = trapz(xline, cos(InnerArg).*dfline) / (L/2);
    bn(n) = -beta(n)*(2*k/f0)*(A*L/4/pi/n)/besselj(1, 2*pi*n*A/L);
end

U=zeros(1,length(dfline));
F=zeros(1,length(dfline));
DfTest = zeros(1,length(dfline));
DfTest2 = zeros(1,length(dfline));

maxcomp = maxit;
%maxcomp = floor(L/2/A);

for n = 1:maxcomp
    InnerArg = ((2*n*pi)/L).*xline;
    U = U + an(n)*sin(InnerArg) + bn(n)*cos(InnerArg);
    F = F - an(n)*(2*pi*n/L)*cos(InnerArg) + bn(n)*(2*pi*n/L)*sin(InnerArg);
    DfTest = DfTest - (f0/2/k)*an(n)*(4*pi*n/A/L)*besselj(1, 2*pi*n*A/L)*sin(InnerArg) -(f0/2/k)*bn(n)*(4*pi*n/A/L)*besselj(1, 2*pi*n*A/L)*cos(InnerArg);
    DfTest2 = DfTest2 + alpha(n)*sin(InnerArg) + beta(n)*cos(InnerArg);
end