;+
; NAME: hsi_ls_sin_cos.pro
;       
; PURPOSE: 
;       Returns the least-square fit to A*cos(omega*angle)+B*sin(omega*angle)
;       given scalar omegs and array angle.
;
; METHOD:
;  The method is to minimize the total squared deviation H of the fit to
;  y[i] = A * cos(omega*angle[i]) + B * sin(omega*angle[i]):
;
;  H = total( (y - A*cos(omega*angle[i]) - B * sin(omega*angle[i])))^2 )
;
; If the angle[i] are chosen to make the cos and sin terms uncorrelated
; then the matrix used herein is diagonal.  If, they are well correlated
; or nearly anti-correlated, the matrix is singular, and the method fails.
;
; INPUTS:
;       y=approximate sinusoid (vector)
;       angle=vector of angles, possibly irregularly spaced
; 
; OUTPUTS:
;        Best-fit cosine and sine coefficients to fit of A*cosine+B*sine
;        If sigma exists, the RMS deviation of the fit is returned in sigma.
;        If phi has only 2 elements, sigma is necessarily 0.
;        If yfit exits, the best fits to y will be returned in it.
;
; EXAMPLE:
;     Create an array of angles:  
;         angle=findgen(512)*!pi/512 -!pi/4
;     Define variable:
;         sigma=0.
;     Define function values (a simulated modulation profile):
;         omega=0.5
;         y=100*cos(omega*angle) + 10*randomn(seed,512)
;         z=hsi_ls_sin_cos(angle,y,omega,yfit,sigma,AMP_PHA=amp_pha)
;     Plot the function and the fit:
;         plot,angle,y & oplot,angle,yfit,psym=4
;
; RESTRICTIONS:
;       If the data set is enormous, this is slow.  For a way to speed it
;       up, see Press and Rybicki, Ap.J., 338, 277, 1989.
;
; REVISION HISTORY:
;  Version 1.0, Jun 21, 1999, EJS -- Derived from ls_sinusoid
;-

function hsi_ls_sin_cos,angle,y,omega,yfit,sigma,AMP_PHA=amp_pha

; Compute the sine and cosine arrays:

cosine=cos(omega*angle) 
sine=sin(omega*angle)   

            ; The regression matrix
m00=total(cosine^2)  & m01=total(sine*cosine)
m11= total(sine^2)

matr=[[m00, m01], $
      [m01, m11] ]

; Example: if omega*angle=[0,.25,.5,.75]*!pi, then the matrix is [[2,0],[0,2]]
; and its inverse is 0.5 * I, so AB = 0.5* vec

c0=total(y*cosine)
c1=total(y*sine)
del=m00*m11-m01^2

A=(c0*M11-c1*M01)/del
B=(c1*m00-c0*m01)/del

; Add the means back in to compute the goodness of fit
yfit=A*cosine+B*sine
if (n_elements(sigma) GT 0) then  $
   sigma=sqrt(variance(y-yfit))

AB=[A,B]

if keyword_set(amp_pha) then AB=[sqrt(A^2+B^2),atan(B,A)]

return,AB

end