function modprof2visibility,phi,delta,sigma
;+
; NAME: modprof2visibility.pro
;       
; PURPOSE: Returns the best-fit visibility phi_c + i * phi_s from N 
;           independent sinuoidal transforms.
;
; METHOD:
;  The method is to minimize the total squared deviation H of the fit to
;  phi[i] = phi_c * cos(delta[i]) - phi_s * sin(delta[i]),
;
;  H = total( (phi - phi_c*cos(delta) + phi_s * sin(delta))^2 )
;
; If the delta are chosen to make sin(delta) uncorrelated with cos(delta)
; then the matrix used below is diagonal.  If, they are well correlated
; or nearly anti-correlated, the matrix is singular, and the method fails.
;
; INPUTS:
;       Sinusoidal transforms phi=vector of N floating point numbers (N GE 2)
;       phases delta=vector of N floating point numbers
;
;    phi[i] is presumed to be the sinusoidal transform of some function f(x,y):
;    phi = \int dx \int dy f(x,y) cos(k(x+cos(th[i])+y*sin(th[i]))+delta[i])
;    where the x and y are relative to the phase center of the map (0 at ctr)
;
; 
; OUTPUTS:
;        Best-fit cosine and sine transform in complex form.
;        If sigma exists, the RMS deviation of the fit is returned in sigma.
;        If phi has only 2 elements, sigma is necessarily 0.
;     The output is the exponential transform of the same f(x,y):
;     phi_LS= \int dx \int dy f(x,y)*exp(ik(x*cos(th_bar)+y*sin(th_bar)))
;     where th_bar is the mean of th[*],
;
;
; EXAMPLE:
;       
;  delta=!pi*[0,.25,.5,.75]
;  phic=0.5 & phis=.75
;  phi=phic*cos(delta)-phis*sin(delta)+0.02*randomn(seed,4)
;  sigma=0.
;  print,modprof2visibility(phi,delta,sigma),sigma
;  
; RESTRICTIONS:
;  The values of phi must be independent, and cos(delta) and sin(delta)
;  must be uncorrelated arrays.  Phi must have at least 2 elements.
;
; REVISION HISTORY:
;  Version 1.0, March 25, 1999, EJS
;
;-
matr=fltarr(2,2)
matr(0,0)=total(cos(delta)^2)
matr(1,1)=total(sin(delta)^2)
matr(0,1)=-total(sin(delta)*cos(delta))
matr(1,0)=matr(0,1)

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

vec=[total(phi*cos(delta)),-total(phi*sin(delta))]



phi_LS=invert(matr)#vec

sigma=sqrt(total((phi-phi_LS[0]*cos(delta)+phi_LS[1]*sin(delta))^2))


return,phi_LS[0]+complex(0,1)*phi_LS[1]

end