;+
; NAME:
;       quadrupole
;
; PURPOSE:
;
;       The purpose of this program is to compute the ellipticity of an object
;       from its quadrupole moments
;       
; AUTHOR:
; 
;       M. James Jee, Ph.D. (Johns Hopkins University)
;
; CALLING SEQUENCE:
;
;       test=quadrupole(img,e1,e2,width,xc=xc,yc=yc,sigma=sigma,q11,q22,q12)
;
;
; INPUTS:
;       IMG - two-dimensional astronomical image
;       
; OUTPUTS:
;       The function returns 0 if it fails to compute the ellipticity.
;       Otherwise it returns 1.
;
;       e1,e2 - two ellipticity components. The value e=sqrt(e1^2+e2^2)
;               correspond to the ellipticity (a-b)/(a+b), where
;               a and b are the major and minor axis, respectively.
;
;       width - contains the measure of the width of the object
;
;       q11,q22,q12 - contains the quadrupole moments of the object. See the paper.
;
; OPTIONAL KEYWORDS:
;
;       xc,yc - the center of the object. Remember that the subscript of the first pixel
;               is 0 not 1.
;       sigma - the sigma for the weighting Gaussian function. default value is 2.
;
function quadrupole,img,e1,e2,width,xc=xc,yc=yc,sigma=sigma,q11,q22,q12

s=size(img)
xsize=s[1]
ysize=s[2]

makenxy,1,xsize,xsize,1,ysize,ysize,xa,ya
xa=xa-1.d
ya=ya-1.d

tot=total(img)

if not keyword_set(xc) then xc=total(xa*img)/tot
if not keyword_set(yc) then yc=total(ya*img)/tot
if not keyword_set(sigma) then sigma=2.d


xa=xa-xc
ya=ya-yc

rmap=sqrt(xa^2+ya^2)
gmap=exp(-rmap^2/(2.*sigma^2))
gmap=gmap/total(gmap)

tot=total(img*gmap)

q11=total(xa*xa*img*gmap)/tot 
q22=total(ya*ya*img*gmap)/tot 
q12=total(xa*ya*img*gmap)/tot

xc=total(xa*img*gmap)/tot
yc=total(ya*img*gmap)/tot

test=q11*q22 - q12^2


e1=(q11-q22) /  (  q11+q22 + 2.d * sqrt(test) )
e2=2.0d * q12 / (  q11+q22 + 2.d * sqrt(test) )

width=sqrt(q11+q22)

if test lt 0. then return, 0 else return,1

end