function quadratic_basis,x,i,basis_extra=basis_extra
  xx=x[*,0]
  yy=x[*,1]
  case i of
    0:result=0*xx+1
    1:result=xx
    2:result=yy
    3:result=xx*yy
    4:result=xx^2
    5:result=yy^2
  end
  return,result
end
;Linear-combination fitter. Fits data to a linear combination of arbitrary functions
;
;y=A0*f0(x)+A1*f1(x)+A2*f2(x)...
;
;Each function fn can be wildly nonlinear, just the combination has to be linear.
;
;Input:
;  x: Independent variable. May be a 1D or 2D array. Use a 2D array for multiple regression,
;     or in other words if the independent variable is multi-dimensional. Each row is
;     a single (possibly multidimensional) value of the independent variable
;  y: Dependent variable. Must be a 1D array, with the same number of elements as x has rows
;  n: Number of bases. For instance, 3 for a quadratic fit
;  basis: Name of basis function.
;Keyword input:
;  All keyword parameters will be passed to the basis function
;Returns:
;  An array of linear coefficients, such that
;  y_model=A[0]*f0(x)+A[1]*f1(x)+A[2]*f2(x)...
;  and such that the total variance is minimized (least-squares fit)
;
;Basis function:
;  Write a function which can evaluate all the bases. Its header should be as follows
;  function basis,x,i,...
;  where basis can be any function name, x is the independent variable, and i is the
;  basis index to calculate.
function lin_comb_fit,x,y,n,basis,basis_extra=basis_extra
  A=fltarr(n_elements(y),n)
  for i=0,n-1 do begin
    A[*,i]=call_function(basis,x,i,basis_extra=basis_extra)
  end
  AT=transpose(A)
  B=y
  alpha=(AT#A)
  beta=(AT#B)
  return,invert(alpha)#beta
end