;+
;
; NAME       : odcenter (procedure)
; PURPOSE    :
;	obtain center of occulting disk & its radius
; CATEGORY :
;	idl/lib/nkr
; CALLING SEQUENCE :
;	odcenter,img,ox,oy,r,rxy=rxy
; INPUTS :
; 	img     -- 10-cm image data
; OPTIONAL INPUT PARAMETERS : 
;	none
; KEYWORD PARAMETERS :
;	rxy  --  ration between x & y-radius
; OUTPUTS :
;	ox,oy,r  -- center and radius
; COMMON BLOCKS : 	none
; SIDE EFFECTS :
;	if /ver is set, fitting curve is displayed
; RESTRICTIONS :	none
; PROCEDURE :
; MODIFICATION HISTORY :
;	K.I. 23 Jan.1993
;-

;*************************************************************************
;NAME       : get2edges (procedure)
;FUNCTION   : obtain 2 edges of occulting disk from a cross section
;=========================================================================
pro get2edges,cross,edge1,edge2,region=region

s=size(cross)
n=s(1)
if keyword_set(region) then n1=fix(n*region)	$
else n1=n/2

c1=fix(shift(cross(0:n1),1))-fix(cross(0:n1))
c2=fix(shift(cross(n-1-n1:n-1),-1))-fix(cross(n-1-n1:n-1))
mm=max(c1(1:n1),m1)
mm=max(c2(1:n1),m2)
edge1=fix(m1)+1
edge2=fix(m2)+n-1-n1+1+1	; <-- modified +1  '93/01/23

end

;*************************************************************************
;NAME       : odcenter (procedure)
;FUNCTION   : obtain center of occulting disk & its radius
;=========================================================================
pro odcenter,img,ox,oy,r,rxy=rxy
; rxy   --  ration between x & y-radius

nn=3	; no. of edge detection cross line = nn*2+1
x1=intarr(2*nn+1)
x2=intarr(2*nn+1)
yc=intarr(2*nn+1)
y1=intarr(2*nn+1)
y2=intarr(2*nn+1)
xc=intarr(2*nn+1)

s=size(img)	& nx=s(1)	& ny=s(2)
nx2=nx/2
ny2=ny/2

;----- O.D. edge on central cross section -----
get2edges,img(*,ny2:ny2),xmin,xmax,region=0.33
get2edges,transpose(img(nx2:nx2,*)),ymin,ymax,region=0.33
;----- O.D. edge on nn*2 cross sections -----
dx=(xmax-xmin)/(nn+1)/2
dy=(ymax-ymin)/(nn+1)/2
for n=0,nn*2 do begin
	yc(n)=ny2+dy*(n-nn)
	if (yc(n) gt 0.3*ny) and (yc(n) lt 0.7*ny) then	$
		get2edges,img(*,yc(n)),edge1,edge2,region=0.33	$
	else	get2edges,img(*,yc(n)),edge1,edge2
	x1(n)=edge1
	x2(n)=edge2
	xc(n)=nx2+dx*(n-nn)
	if (xc(n) gt 0.3*nx) and (xc(n) lt 0.7*nx) then	$
		get2edges,transpose(img(xc(n),*)),edge1,edge2,region=0.33 $
	else	get2edges,transpose(img(xc(n),*)),edge1,edge2
	y1(n)=edge1
	y2(n)=edge2
endfor
;----- obtain center pos. & average radius -----
ox = float(total(x1)+total(x2))/float(2*(nn*2+1))
oy = float(total(y1)+total(y2))/float(2*(nn*2+1))
r = ( total(sqrt((x1-ox)^2+(yc-oy)^2))+total(sqrt((x2-ox)^2+(yc-oy)^2))	  $ 
    + total(sqrt((xc-ox)^2+(y1-oy)^2))+total(sqrt((xc-ox)^2+(y2-oy)^2)) ) $ 
			/float(4*(nn*2+1))
;----- obtain ratio between x & y radius -----
;if keyword_set(rxy) then begin
    for n=0,nn*2 do begin
	get2edges,img(*,oy-nn+n:oy-nn+n),edge1,edge2,region=0.33
	x1(n)=edge1
	x2(n)=edge2
	get2edges,transpose(img(ox-nn+n:ox-nn+n,*)),edge1,edge2,region=0.33
	y1(n)=edge1
	y2(n)=edge2
    endfor
    rx=float(total(x2)-total(x1))/float(2*(nn*2+1))
    ry=float(total(y2)-total(y1))/float(2*(nn*2+1))
    rxy=rx/ry
;endif

end