;+
; NAME:
;   SFITCUSP
;
; PURPOSE:
;   This function determines a "cusp" polynomial weighted fit to a surface.
;
; CATEGORY:
;   Curve and surface fitting.
;
; CALLING SEQUENCE:
;   Result = SFITCUSP( Degree, X, Y, Z, W)
;
; INPUTS:
;
;   Degree:   The maximum degree of fit (in one dimension), [integer].
;
;      X,Y:   The x- and y-coordinates for each data point, [double(npt)].
;
;        Z:   The surface value for each data point, [double(npt)].
;
;        W:   The weight associated with each data point, [double(npt)].
;         For poisson statistics use 1/sigma2 = 1/Z.
;
; INPUT KEYWORDS:
;
; DEADTIME:   Set this keyword to include a dead time correction in the fit.
;
; OUTPUT:
;   This function returns a fitted array, [double(npt)].
;
; OUTPUT KEYWORDS:
;    COEFF:   The array of coefficients for a cusp-polynomial function
;         of x and y to fit data.
;         This parameter is returned as a (Degree+1)^2 element [double] array.
;
;      ERR:   The error matrix of dimension (Degree+1)^2 x (Degree+1)^2
;         containing the variances and covariances of the fitted coefficients.
;
; PROCEDURE:
;   Fit a 1D array Z as a polynomial function of x and y.
;   The function fitted is:
;       F(x,y) = Sum over i and j of coeff(i*(degree+1)+j) * abs( x^i * y^j )
;   where coeff is returned as a keyword.
;
;   The deadtime correction is of the form
;       D(z)   = - coeff * cts^2
;   where coeff = (deadtime)/(timing mode=(320 or 5 ms)
;   WARNING: This functional form for the deadtime correction is only an
;
;
; MODIFICATION HISTORY:
;   July, 1994               H.C. Wen
;   07-SEP-1994              Added the DEADTIME keyword.
;
;-
function SFITCUSP, degree, x, y, z, w, DEADTIME=Cts, COEFF=coeff, ERR=err

         on_error, 2
         npts = N_ELEMENTS( x )
         ndim = (degree+1)^2
         if keyword_set( Cts ) then ndim = ndim + 1

         beta = dblarr( ndim, /nozero )
         alpha= dblarr( ndim, ndim, /nozero )

;   Enter matrix and vector elements
;   for the cusp-polynomial

         for jj=0,degree do for kk=0,degree do begin
              beta( jj*(degree+1)+kk) = TOTAL( z*w*abs( x^jj*y^kk ) )
              for j=0,degree do for k=0,degree do $
                   alpha( jj*(degree+1)+kk, j*(degree+1)+k ) =$
                   TOTAL( w*abs( x^(j+jj)*y^(k+kk) ) )
         endfor

;   for the dead time correction

         if keyword_set( Cts ) then begin
              beta( ndim-1 )      = TOTAL( z*w* (-cts^2) )
              for j=0,degree do for k=0,degree do $
                   alpha( ndim-1, j*(degree+1)+k ) = $
                        TOTAL( w*abs( x^j * y^k )*( -cts^2 ) )
              alpha( *, ndim-1 )  = alpha( ndim-1, * )
              alpha( ndim-1, ndim-1 ) = TOTAL( w* (-cts^2)^2 )
         endif

;   Invert matrix to determine regression coefficients

         detA = DETERM( alpha )
         if detA eq 0 then message, 'Determinant = 0!'

         err   = INVERT( alpha )       ;error matrix
         coeff = err # beta

         fit = dblarr( npts)
         for j=0,degree do for k=0,degree do $
              fit = fit+coeff( j*(degree+1)+k )*abs( x^j*y^k )

         if keyword_set( Cts ) then $
              fit = fit + coeff( ndim-1 )*( -cts^2 )

         return, fit
end