function [xFine, fFine] = m20120528_03_FourierInterpolation(f, xMin, xMax, NFine, boundary)

% Interpolate a periodic function from one uniform grid to a finer grid,
% using Fourier transforms to achieve spectral accuracy.

% Possible values for boundary:
% 1 = first element is xMin, last element is xMax-dx.
% 2 = first element is xMin+dx, last element is xMax.
% 3 = first element is xMin, last element is xMax, so f(1)=f(end)

if NFine==numel(f)
    fFine=f;
    switch boundary
        case 1
            xFine = linspace(xMin,xMax,NFine+1);
            xFine(end)=[];
        case 2
            xFine = linspace(xMin,xMax,NFine+1);
            xFine(1)=[];
        case 3
            xFine = linspace(xMin,xMax,NFine);
        otherwise
            error('Invalid value of boundary')
    end
    return
end
if NFine < numel(f)
    error('Not set up yet for reducing resolution')
end

column=1;
row=2;
if iscolumn(f)
    orientation=column;
    f=f';
else
    orientation=row;
end
% Internally, everything will be done with rows.  We'll transpose at the
% end if needed;

switch boundary
    case 1
        xFine = linspace(xMin,xMax,NFine+1);
        xFine(end)=[];
    case 2
        xFine = linspace(xMin,xMax,NFine+1);
        xFine(1)=[];
        f = [f(end), f(1:(end-1))];
    case 3
        if abs(f(1)- f(end))> 1e-13
            error('For boundary=3, f(1) must equal f(end)')
        end
        xFine = linspace(xMin,xMax,NFine);
        f(end)=[];
        NFine=NFine-1;
    otherwise
        error('Invalid value of boundary')
end

NCoarse = numel(f);

ff=fft(f);

if mod(NCoarse,2)==0
    halfN = ceil(NCoarse/2);
    ffFine = [ff(1:(halfN)), 0.5*ff(halfN+1), zeros(1,NFine-NCoarse-1), 0.5*ff(halfN+1), ff((halfN+2):end)] * NFine/NCoarse;
else
    halfN = ceil(NCoarse/2);
    ffFine = [ff(1:halfN), zeros(1,NFine-NCoarse), ff((halfN+1):end)] * NFine/NCoarse;
end
fFine = ifft(ffFine);


switch boundary
    case 2
        fFine = [fFine(2:end), fFine(1)];
    case 3
        fFine = [fFine, fFine(1)];
end

if orientation==column
    xFine=xFine';
    fFine=fFine';
end

end