;------------------------------------------------------------------
;+
; NAME:
;	MOSTOPSLNS
; PURPOSE:
;       Rescale BBSO polar mosaic to constant size Sun, correct for
;       P angle, and display from pole to 40 deg.  This version
;       includes an option to turn off the conversion to 8 bits.
; CATEGORY:
;	POLAR
; CALLING SEQUENCE:
;       newimg = mostopslns(stk,stkhdr,imgscl,xcen,ycen,sclout,xcout,ycout)
; INPUTS:
;       stk    = BBSO polar mosaic
;       stkhdr = FITS header (as string array) for date/time info
;       imgscl = BBSO mosaic image scale, pixels/arcsec
;       xcen   = position of Sun center for mosaic (lower left corner
;                of image is (0,0)
;       ycen   = position of Sun center for mosaic
; KEYWORD PARAMETERS:
;       noscale = if set, do not scale to 8 bits
; OUTPUTS:
;       newimg    = rescaled, rotated image
;       sclout    = image scale of output image
;       xcout     = Sun center position in output image
;       ycout     = Sun center position in output image
;
; COMMON BLOCKS:
;       SUNRAD contains:
;          R      = radius of Sun in pixels.
;          Q      = quadrant of circle to fit
; NOTES:
; MODIFICATION HISTORY:
;       J. Varsik,      11 Nov. 1999
;-
;-------------------------------------------------------------------

FUNCTION mostopslns,stk,stkhdr,imgscl,xcen,ycen,sclout,xcout,ycout,noscale=nscl
COMMON SUNRAD,sprad,quad

    ; Display IDL header if help is required.
    IF (KEYWORD_SET(help)) THEN BEGIN
      GET_IDLHDR,'mostopslns.pro'
      RETURN, 0.0
    ENDIF

    TRUE = 1
    FALSE = 0

    IF (KEYWORD_SET(nscl)) THEN ns = 1 ELSE ns = 0

    ; Obtain FITS parameters from V image
    datec = FXPAR(stkhdr,'DATE-OBS')
    timec = FXPAR(stkhdr,'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
    PRINT,'P = ', p, 'B = ', b, 'SD = ', rsun
;   sprad is radius of Sun in pixels.

;   Use fixed radius of 1000 pixels.  Final image is 1500 x 300, with
;   the projected pole (the B angle is not removed) centered at
;   (750,300) for north pole images and (750, 0) for south pole images.
;   so we make a temporary picture of 3000x3000 pixels.

;   Clean up mosaic image to remove noise beyond limb
    isz = SIZE(stk)
    ixz = isz(1)
    iyz = isz(2)
    filled = WHERE(stk NE -32000, fc)
    IF fc EQ 0 THEN PRINT,'empty mosaic'
    i = filled MOD ixz
    j = filled / ixz
    xfc = (i+1)-xcen
    yfc = (j+1)-ycen
    r1 = SQRT((xfc)^2 + (yfc)^2)
    tmpimg = FLTARR(ixz,iyz)
    tmpimg[filled] = r1
    big = WHERE(tmpimg GE sprad)
    stkc = stk
    stkc[big] = -32000

;   Rescale BBSO mosaic to fixed solar radius, make 8-bit version
    stksiz = SIZE(stk)
    sfactor = 1000.0 / sprad
    PRINT, sfactor
    stkmidsize = CONGRID(stkc,FIX(sfactor*stksiz(1)),FIX(sfactor*stksiz(2)))
    stkmidx = xcen * sfactor
    stkmidy = ycen * sfactor
    PRINT, stkmidx, stkmidy
    IF (NOT ns) THEN stkbyt = BYTSCL(stkmidsize,MIN=-30,MAX=30) $
        ELSE stkbyt = stkmidsize

;   Make blank picture, place rescaled mosaic on picture
    IF (NOT ns) THEN stkmask = BYTARR(3000, 3000) ELSE $
         stkmask = MAKE_ARRAY(3000, 3000,/INT,VALUE=-32000)
    stkmask[1500.0-stkmidx,1500.0-stkmidy] = stkbyt
    stkmf = ROTATE(stkmask, 5)
;   Rotate image to correct for P angle about correct center point
;    stkrot = ROT(stkmask,p,1.0,1500.0,1500.0)
    stkrot = ROT(stkmf,p,1.0,1500.0,1500.0)
;   Mask non-image pixels
    picmask = REPLICATE(-32000,3000,3000)
    picmask[1500.0-stkmidx,1500.0-stkmidy] = stkmidsize
;    picmf = ROTATE(picmask, 5)
;    picrot = ROT(picmf,p,1.0,1500,1500)
;    picrot = ROT(picmf,p,1.0)
    picrot = ROT(picmask,p,1.0)
    good = WHERE (picrot NE -32000)

    picmask(good) = stkrot(good)

;   Output image scale (pixels/arcsec)
    sclout = 1000.0 / rsun

;   Extract output image
    IF ycen LT 0 THEN BEGIN
;       north pole image
;        outimg = picmask[542:1957,1958:2249]
        outimg = stkrot[730:2269, 2140:2499]
        xcout = 769
        ycout = -641
    ENDIF ELSE BEGIN
;        outimg = picmask[542:1957,250:542]
        outimg = stkrot[730:2269, 500:859]
        xcout = 769
        ycout = 999
    ENDELSE

    RETURN,outimg
    END