;------------------------------------------------------------------
;+
; NAME:
;	GETL
; PURPOSE:
;	Get limb points
;	Find center Sun position in BBSO polar magnetograms.
;	This procedure estimates the Sun center position in pixel
;	coordinates.  This procedure uses the center X, Y position
;	(provided by the guider) along with image scale information
;	provided by the user.
; CATEGORY:
;	POLAR
; CALLING SEQUENCE:
;	limbimg = centerls(vfile,jfile,imgscl,gcut)
; INPUTS:
;	vfile = file name for input magnetogram (FITS format).
;	jfile = file name for input direct image (FITS format).
;	imgscl = image scale (pixels / arcsec).
;       gcut   = gradient cutoff value
;	gcorr  = correction to picture center coordinates.  The
;	         real coordinates are the file coordinates divided
;	         by gcorr.  If gcorr is zero or not present, no correction
;	         is applied.
; KEYWORD PARAMETERS:
; OUTPUTS:
;	an image includling the limb points.
;
; COMMON BLOCKS:
; NOTES:
;	The time, date, and position information are taken from the
;	V image.  The FITS keywords are added in BOTH images.
;	This assumes that the time and position are the same for both
;	images.
; MODIFICATION HISTORY:
;	J. Varsik,	18 Mar. 1996
;-
;-------------------------------------------------------------------

function getl,vfinam,jfinam,imgscl,gcut,c_corr=gcorrk

    ; Display IDL header if help is required.
    if (keyword_set(help)) then begin
      get_idlhdr,'centere.pro'
      goto,finishup
    endif

    ; Check for coordinate correction keyword
    if (keyword_set(gcorrk)) then gcorr=gcorrk else gcorr = 0.0

    ; Read V file
    vi = READFITS(vfinam,vih,/NOSCALE)
    ji = READFITS(jfinam,jih,/NOSCALE)

    ; Obtain FITS parameters from V image
    ew = FXPAR(vih,'CRVAL1')
    ns = FXPAR(vih,'CRVAL2')
    datec = FXPAR(vih,'DATE-OBS')
    timec = FXPAR(vih,'TIME-OBS')
    READS,datec,day,mon,yrs,FORMAT='(i2,1x,i2,1x,i2)'
    READS,timec,hour,min,sec,FORMAT='(i2,1x,i2,1x,i2)'
    hr = hour + ((min + (sec / 60.0)) / 60.0)
    pandb,mon,day,yrs,hr,p,b,rsun,mjd
    sprad = rsun * imgscl
    sdrad = 0.0
    if (gcorr EQ 0.0) then begin
        h = imgscl * ew + 256
        k = (-1.0 * imgscl * ns) + 192
    endif else begin
        h = imgscl * (ew / gcorr) + 256
        k = (-1.0 * imgscl * (ns / gcorr)) + 192
    endelse

    print,h,k
; now we have our first estimate of the Sun center position.  Next
; we isolate the limb points.  We use the Sobel function to
; give the gradient of the direct image.
    jgrad = SOBEL(ji)

; Next we use the region-coloring function to identify all the
; regions of the direct image with large gradients.
    jgcut = jgrad GT gcut
    jglabel = LABEL_REGION(jgcut)
    jgmag = HISTOGRAM(jglabel,REVERSE_INDICES=r)
    print,jgmag(0:10)

; Now we select the largest region.  This should be the set of limb
; points.
    imsiz = SIZE(ji)
    xax = imsiz(1)
    jgmagmax = MAX(jgmag(1:*),maxreg)
    maxreg = maxreg + 1
    print,jgmagmax,maxreg
    p = r(r(maxreg):r(maxreg+1)-1)
    xl1 = INTARR(N_ELEMENTS(p))
    yl1 = INTARR(N_ELEMENTS(p))
    FOR i = 0, N_ELEMENTS(p)-1 DO BEGIN
      kpix = p(i)
      xl1(i) = kpix MOD xax
      yl1(i) = kpix / xax
    ENDFOR

; Be sure to avoid bad points.
    ig = WHERE(xl1 GT 11 AND xl1 LT 504)
    xl = xl1(ig)
    yl = yl1(ig)

    limbi = intarr(xax, imsiz(2))
    for i = 0, n_elements(xl)-1 do begin
       limbi(xl(i),yl(i)) = 1
    endfor
    finishup:

    RETURN,limbi
END