;------------------------------------------------------------------
;+
; NAME:
;	MBIGMAP
; PURPOSE:
;
;       This program takes a mosaic of BBSO VMG images 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, and running from 0 to
;	360 in Carrington longitude in x.
;
;	The scale is 8 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 = MBIGMAP,im,imhdr,h,k,imgscl
; 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.
; KEYWORD PARAMETERS:
; OUTPUTS:
;	A map.  It is 2880 x 360 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.
;-
;-------------------------------------------------------------------
FUNCTION MBIGMAP,im,imh,h,k,imgscl,help=help

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

    ; Find center, radius, p angle, b angle, etc.
    datec = FXPAR(imh,'DATE-OBS')
    timec = FXPAR(imh,'TIME-OBS')
    READS,datec,yrs,mon,day,FORMAT='(i4,1x,i2,1x,i2)'
    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(2880,360)
    s = rsun
    rs2 = rsun^2
    TRUE = 1
    FALSE = 0

    IF (k LT 1000) THEN npole = TRUE ELSE npole = FALSE
    IF (npole) THEN PRINT,'North pole image' ELSE PRINT,'South pole image'

    ; Omit flux conservation stuff for now

    PRINT,'Begin scanning output image'
    FOR i = 0, 2879 DO BEGIN
        hi = i / 8.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, 359 DO BEGIN
            IF (lambda LT -90.0 OR lambda GT 90.0) THEN BEGIN
                cart(i,j) = -32767
            ENDIF ELSE BEGIN
               IF (npole) THEN phi = (j/8.0) + 45.0  ELSE phi = -90.0 + (j/8.0)
                ; 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 (ix GE 0 AND ix LE 2999 AND iy GE 0 AND iy LE 1499) $
                  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
       ENDFOR
    ENDFOR

    
    finishup:

    RETURN,cart
end