;------------------------------------------------------------------
;+
; NAME:
;	MCENLS
; PURPOSE:
;	Find center Sun position in a mosaic of BBSO polar magnetograms.
;	This procedure estimates the Sun center position in pixel
;	coordinates.  It also provides a revised estimate of the true
;       image scale.
; CATEGORY:
;	POLAR
; CALLING SEQUENCE:
;	imgsclnew = mcenls(jimg,jih,imgscl,gcut, xout, yout)
; INPUTS:
;       ji     = IDL image for input direct image mosiac.
;       jih    = FITS header for the V image taken closest to the
;                middle of the frames used for the mosaic.  This is
;                used to get the time for the ephemeris.
;	imgscl = image scale (pixels / arcsec).
;       gcut   = gradient cutoff value
; KEYWORD PARAMETERS:
;       south   = if set, the image is a south polar mosaic.  If not,
;                 it's a north pole mosaic.
; OUTPUTS:
;       imgsclnew = new image scale value.
;       xout      = position of center Sun in pixels.
;       yout      = position of center Sun in pixels. 
;
; COMMON BLOCKS:
;       SUNRAD contains:
;          R      = radius of Sun in pixels.
;          Q      = quadrant of circle to fit
; NOTES:
;       The position of center Sun is given in pixels.  The image
;       scales are in pixels/arcsec.  The origin of the coordinates
;       is with (0,0) at the lower left corner of the image(s).
;       The initial coordinates for Sun center are (1500,0) for the
;       north polar mosaics and (1500,1500) for the south polar
;       mosaics.
; MODIFICATION HISTORY:
;	J. Varsik,	11 Mar. 1999
;       J. Varsik,      19 Apr. 1999 changed mosaic size to 3000x1500.
;       J. Varsik,      06 May  1999 changed back to 4000x2000, allow
;                                    old and new FITS dates to be used
;-
;-------------------------------------------------------------------

FUNCTION mcenls,ji,jih,imgscl,gcut,xout,yout,SOUTH=south
COMMON SUNRAD,sprad,quad

    ; Display IDL header if help is required.
    if (keyword_set(help)) then begin
      GET_IDLHDR,'mcenls.pro'
      RETURN, 0.0
    endif

    TRUE = 1
    FALSE = 0
    ; Check for south pole keyword
    IF (KEYWORD_SET(SOUTH)) THEN npole = FALSE ELSE npole = TRUE

    ; Obtain FITS parameters from V image
    datec = FXPAR(jih,'DATE-OBS')
    timec = FXPAR(jih,'TIME-OBS')
    IF (STRMID(datec, 2, 1) EQ '/') THEN BEGIN
        READS,datec,day,mon,yrs,FORMAT='(i2,1x,i2,1x,i2)'
        yrs = yrs + 1900
    ENDIF ELSE BEGIN
        READS,datec,yrs,mon,day,FORMAT='(i4,1x,i2,1x,i2)'
    ENDELSE
    READS,timec,hour,min,sec,FORMAT='(i2,1x,i2,1x,i2)'
    hr = hour + ((min + (sec / 60.0)) / 60.0)
;   Use solar ephemeris from Johns Hopkins library
    sun,yrs,mon,day,hr,pa=p,lat0=b,sd=rsun
;    pandb,mon,day,yrs,hr,p,b,rsun,mjd
    sprad = rsun * imgscl

;    h = 1500.0
    h = 2000.0
    IF (npole) THEN BEGIN
        k = 0.0
        quad = 1
    ENDIF ELSE BEGIN
;        k = 1500.0
        k = 2000.0
        quad = 3
    ENDELSE
; 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)

; 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
    p = r(r(maxreg):r(maxreg+1)-1)
    xl = INTARR(N_ELEMENTS(p))
    yl = INTARR(N_ELEMENTS(p))
    FOR i = 0, N_ELEMENTS(p)-1 DO BEGIN
      kpix = p(i)
      xl(i) = kpix MOD xax
      yl(i) = kpix / xax
    ENDFOR

; Now we use the non-linear least squares curve-fitting function
; to locate the new Sun center points.
    a = FLTARR(3)
    a(0) = h
    a(1) = k
    a(2) = sprad
    w = MAKE_ARRAY(N_ELEMENTS(xl),/FLOAT,VALUE=1.0)
    yfit = CURVEFIT(xl,yl,w,a,siga,FUNCTION_NAME='NSLIMB2')
    h = a(0)
    k = a(1)
    sprad = a(2)
    sigh = siga(0)
    sigk = siga(1)
    sigr = siga(2)
    PRINT,h,sigh,k,sigk,sprad,sigr

; Now return the variables we want to return.
    xout = h
    yout = k
    newimgscl = sprad/rsun
    RETURN,newimgscl
END

PRO NSLIMB,x,a,f,pder
  COMMON SUNRAD,r,q
  ; function for use with CURVEFIT. 
  ; r is set to solar radius in pixels.
  CASE q OF

  1: BEGIN
    f = a(1) + SQRT(r^2 - (x - a(0))^2)
    IF N_PARAMS() GE 4 THEN BEGIN
      pder = FLTARR(N_ELEMENTS(x),2)
      pder(*,0) = (x - a(0)) / SQRT(r^2 - (x - a(0))^2)
      pder(*,1) = 1.0
    ENDIF
  END
  3: BEGIN
    f = a(1) - SQRT(r^2 - (x - a(0))^2)
    IF N_PARAMS() GE 4 THEN BEGIN
      pder = FLTARR(N_ELEMENTS(x),2)
      pder(*,0) = -1.0*(x - a(0)) / SQRT(r^2 - (x - a(0))^2)
      pder(*,1) = 1.0
    ENDIF
  END
  ENDCASE
END
PRO EWLIMB,y,a,f,pder
  COMMON SUNRAD,r,q
  ; function for use with CURVEFIT. 
  ; r is set to solar radius in pixels.
  CASE q OF

  2: BEGIN
    f = a(0) + SQRT(r^2 - (y - a(1))^2)
    IF N_PARAMS() GE 4 THEN BEGIN
      pder = FLTARR(N_ELEMENTS(y),2)
      pder(*,0) = 1.0
      pder(*,1) = (y - a(1)) / SQRT(r^2 - (y - a(1))^2)
    ENDIF
  END
  4: BEGIN
    f = a(0) - SQRT(r^2 - (y - a(1))^2)
    IF N_PARAMS() GE 4 THEN BEGIN
      pder = FLTARR(N_ELEMENTS(x),2)
      pder(*,0) = 1.0
      pder(*,1) = -1.0*(y - a(1)) / SQRT(r^2 - (y - a(1))^2)
    ENDIF
  END
  ENDCASE
END
PRO NSLIMB2,x,a,f,pder
  COMMON SUNRAD,r,q
  ; function for use with CURVEFIT. 
  ; r is set to solar radius in pixels.
  CASE q OF

  1: BEGIN
    f = a(1) + SQRT(a(2)^2 - (x - a(0))^2)
    IF N_PARAMS() GE 4 THEN BEGIN
      pder = FLTARR(N_ELEMENTS(x),3)
      pder(*,0) = (x - a(0)) / SQRT(a(2)^2 - (x - a(0))^2)
      pder(*,1) = 1.0
      pder(*,2) = a(2) / SQRT(a(2)^2 - (x - a(0))^2)
    ENDIF
  END
  3: BEGIN
    f = a(1) - SQRT(a(2)^2 - (x - a(0))^2)
    IF N_PARAMS() GE 4 THEN BEGIN
      pder = FLTARR(N_ELEMENTS(x),3)
      pder(*,0) = -1.0*(x - a(0)) / SQRT(a(2)^2 - (x - a(0))^2)
      pder(*,1) = 1.0
      pder(*,2) = -1.0*a(2) / SQRT(a(2)^2 - (x - a(0))^2)
    ENDIF
  END
  ENDCASE
END