function [g, dg, ddg] = gaussian(x, a, x0, delta)
% gaussian Returns gaussian and its first two derivatives [PHYS210-fda]
%
% Utility function that returns gaussian(x; a,x0,delta) as well as its
% first and second derivatives.
%
% Input arguments
%
%      x:     (real vector) Spatial grid, length nx. 
%      a:     (real scalar) Overall amplitude of gaussian.
%      x0:    (real scalar) x location of centre of gaussian.
%      delta: (real scalar) Effective width of gaussian.
%
% Output arguments (all of length nx)
%
%      g:      (vector) Gaussian profile g(x(j))
%      dg:     (vector) First derivative of gaussian profile 
%              d[g(x(j))]/dx
%      dg:     (vector) Second derivative of gaussian profile,
%              d^2[g(x(j))]/dx^2
%
% Assumes that x is monotonically increasing.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Define default vectors to be 0 length so that *some* values are 
   % returned if an invalid input argument is detected.
   g   = zeros(0,0);
   dg  = zeros(0,0);
   ddg = zeros(0,0);
   if ~ isvector(x) 
      fprintf('gaussian: x is not a vector\n');
      return
   end

   % Define g, dg, dgg to be vectors of length nx.
   nx = length(x);
   g   = zeros(1,nx);
   dg  = zeros(1,nx);
   ddg = zeros(1,nx);
   if delta == 0
      fprintf('gaussian: delta = 0\n');
      return
   end

   if min(x) <= x0 & x0 <= max(x)
      delm2 = 1 / delta^2;
      m2delm2 = -2 / delta^2;
      g   = a .* exp( -delm2 .* (x - x0).^2 );
      dg  = m2delm2 .* (x - x0) .* g;
      ddg = m2delm2 .* (1.0 + m2delm2 .* (x - x0).^2) .* g;
   end

end
