;------------------------------------------------------------------
;+
; NAME:
;	SBIGMAP
; PURPOSE:
;
;       This program takes a mosaic of BBSO VMG images (or a single
;       image) and projects it onto a square grid of latitude
;       vs. longitude.  Currently, the program is set up for polar images.
;
;	This program is based on the simpmap program.  It uses a grid
;	running from 45 to 90 deg (for north polar images) or 90 to
;	45 deg (for south polar images) in y.  In x, the grid
;       runs from c - 25 deg to c + 25 deg, where c is the central
;       Carrington longitude specified in the command line.
;
;	The scale is 16 pixels per degree in each direction.
;
;       This IDL version does not require a Fortran "core".
;
;	The program produces one FITS image (...... size),
;	which contains the map itself.  The map image name
;	is the original image name with xy appended.
;
;	Routines used:
;           sun (from Johns Hopkins library) to find P, B, and Carrington long.
; CATEGORY:
;	POLAR
; CALLING SEQUENCE:
;	map = SBIGMAP,im,imhdr,h,k,imgscl, ccentral
; INPUTS:
;       im     = BBSO VMG image mosaic (4000x2000 size).
;       imhdr  = A FITS header.  This is used to get date and time
;                information.
;       h,k    = X and Y position of the center of the Sun in pixels.
;                This is based on the lower left corner of the mosaic
;                being (0,0).  The nominal position for a north polar
;                image is (2000,0), for a south polar image
;                (2000,2000).
;       imgscl = image scale in pixels/arcsec.
;       ccentral = central Carrington longitude for map.
; KEYWORD PARAMETERS:
;       south    If set, image to map is a south pole image, otherwise
;                it is a north pole image.
; OUTPUTS:
;	A map.  It is 800 x 720 in size.
; COMMON BLOCKS:
; NOTES:
; MODIFICATION HISTORY:
;       J. Varsik,      12 Mar. 1999
;       J. Varsik,      17 Mar. 1999 Include secant correction for
;       mag. field.
;       J. Varsik,      19 Apr. 1999 Change mosaic size to 3000x1500.
;       J. Varsik,      03 May  1999 Change map size to 400 x 360.
;                                    Add central Carrington long.
;       J. Varsik,      07 May  1999 Change mosaic size to 4000x2000.
;                                    (sbigmaps handles 3000x1500)
;-
;-------------------------------------------------------------------
FUNCTION SBIGMAP,im,imh,h,k,imgscl,ccentral,SOUTH=south,SINGLE=single,HELP=help

    ; Display IDL header if help is required.
    IF (KEYWORD_SET(help)) THEN BEGIN
      GET_IDLHDR,'sbigmap.pro'
      GOTO,finishup
    ENDIF

    TRUE = 1
    FALSE = 0

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

    ; Check for single BBSO VMG image
    IF (KEYWORD_SET(SINGLE)) THEN bbsingle = TRUE ELSE bbsingle = FALSE

    ; Find center, radius, p angle, b angle, etc.
    datec = FXPAR(imh,'DATE-OBS')
    timec = FXPAR(imh,'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,long0=lnaught
    rad = imgscl * rsun

    ; Check parameters
    PRINT, 'p = ',p,' b = ',b,' lnaught = ',lnaught,' rsun = ',rsun
    PRINT, 'h = ',h,' k = ',k,' rad = ',rad
    ; First make a subimage that only has the points we want to use.

    ; Make output image
    cart = FLTARR(800,720)
    s = rsun
    rs2 = rsun^2
    TRUE = 1
    FALSE = 0

    IF (npole) THEN PRINT,'North pole image' ELSE PRINT,'South pole image'

    ; If image is single BBSO VMG, trim it
    IF (bbsingle) THEN BEGIN
        tmpim = im(11:503,*)
        h = h - 11.0
     ENDIF

    ; Omit flux conservation stuff for now

    PRINT,'Begin scanning output image'
    FOR i = 0, 799 DO BEGIN
        hi = (i / 16.0) + (ccentral - 25.0)
        IF (hi GT 360.0) THEN hi = hi - 360.0
        IF (hi LT - 0.0) THEN hi = hi + 360.0
        lambda = lnaught - hi
        IF (lambda LT -180.0) THEN lambda = lambda + 360.0
        IF (lambda GT 180.0) THEN lambda = lambda - 360.0
        FOR j = 0, 719 DO BEGIN
            IF (lambda LT -90.0 OR lambda GT 90.0) THEN BEGIN
                cart(i,j) = -32767
            ENDIF ELSE BEGIN
               IF (npole) THEN BEGIN
                   phi = (j/16.0) + 45.0
               ENDIF ELSE BEGIN
                   phi = -90.0 + (j/16.0)
               ENDELSE
                ; Convert from heliographic to geocentric
                suncor,phi,lambda,p,b,s,x,y

               ix = FIX(h - (x * imgscl)) - 1
                ; Note east is positive in pixel space, west positive
                ; in our geocentric coordinates
               iy = FIX((y * imgscl) + k) - 1
                ; Check input picture, apply to output
                ; picture and scale picture.
                ; This version uses nearest-neighbor interpolation
               IF (bbsingle) THEN BEGIN
                 IF (ix GE 0 AND ix LE 492 AND iy GE 0 AND iy LE 383) $
                  THEN BEGIN
                    pr2 = x^2 + y^2
                    IF (pr2 NE rs2) THEN BEGIN
                        cart(i,j) = tmpim(ix,iy) / SQRT(1.0 - (pr2/rs2))
                    ENDIF ELSE BEGIN
                        cart(i,j) = 0.0
                    ENDELSE
                 ENDIF ELSE BEGIN
                    cart(i,j) = -32767
                 ENDELSE
               ENDIF ELSE BEGIN
                 IF (ix GE 0 AND ix LE 3999 AND iy GE 0 AND iy LE 1999) $
                  THEN BEGIN
                    pr2 = x^2 + y^2
                    IF (pr2 NE rs2) THEN BEGIN
                        cart(i,j) = im(ix,iy) / SQRT(1.0 - (pr2/rs2))
                    ENDIF ELSE BEGIN
                        cart(i,j) = 0.0
                    ENDELSE
                 ENDIF ELSE BEGIN
                    cart(i,j) = -32767
                 ENDELSE
               ENDELSE
          ENDELSE
       ENDFOR
    ENDFOR

    
    finishup:

    RETURN,cart
end