% Given a pointer to a function that represents the Fourier transform
% of a differential operator, return the free space Green's function.
function [G,t,fft_G,G2,fft_G2] = gf(FDfunc, N, dt)
% FDfunc = pointer to function that returns the FFT of some differential
%     operator, D, as a function of an input array of angular frequencies.
% G = Green's function with N elements, sampled at intervals of dt. This
%     is in the time domain.
%
% G will have zero offset. That is, the impulse is at the first element
% of the time array.
%
% It is essential that the units of dt match the time units used by FDfunc.
% That is, if FDfunc expects radians per fortnight, then dt should be in
% fortnights. 
%
% t = array of times, calculated here for convenient plotting by the end user.

omega = omega_array(N, dt); % angular frequency (radians per unit time)
fft_D = FDfunc(omega); % FFT of differential operator for all omega values
fft_G = 1 ./ fft_D; % FFT of Green's function

G = ifftn(fft_G); % Green's function

fft_G2 = anti_alias(fft_G);
G2 = ifftn(fft_G2); % Anti-aliased version of Green's function

t = ( 0:(N-1) ) * dt; % t, so that one can plot G(t).

