;********************************************************************
;NAME       : cgspot (procedure)
;FUNCTION   : get the center of gravity of sunspot
;					K.I. '92/09/02 
;====================================================================
pro cgspot,img0,xc,yc,ww=ww

s=size(img0)
nx=s(1)	&	ny=s(2)

; ------- find reference spot in img's ----------
dd=30	; distance for taking 2-nd derivative

wx=fix(nx*0.8)	&	wy=fix(ny*0.8)	; area for serching spot 
xs=nx/2-wx/2	&	xe=xs+wx-1
ys=ny/2-wy/2	&	ye=ys+wy-1
img=img0(xs:xe,ys:ye)/2
deriv = shift(img,dd,0)+shift(img,-dd,0)-2*img	$
	+shift(img,0,dd)+shift(img,0,-dd)-2*img
max1=max(deriv(dd:wx-dd-1,dd:wy-dd-1),ii)
yp1=ii/(wx-2*dd)+dd
xp1=ii-(yp1-dd)*(wx-2*dd)+dd+xs
yp1=yp1+ys

; ------- get gravity center of the spot ---------
if keyword_set(ww) then begin
	wxf=ww
	wyf=ww
endif else begin
	wxf=fix(64*nx/512)
	wyf=fix(64*ny/480)	; area for calculating G-C around the spot
endelse

xsf1=xp1-wxf/2
ysf1=yp1-wyf/2
part1=float(img0(xsf1:xsf1+wxf-1,ysf1:ysf1+wyf-1))
cx1=rebin(part1,wxf,1)
xc=total((max(cx1)-cx1)*lindgen(wxf))/total((max(cx1)-cx1))+xsf1
cy1=transpose(rebin(part1,1,wyf))
yc=total((max(cy1)-cy1)*lindgen(wyf))/total((max(cy1)-cy1))+ysf1

end