;------------------------------------------------------------------
;+
; NAME:
;	BIGMAP
; PURPOSE:
;
;	     This program takes a squared, clipped 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, 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:
;	BIGMAP,imname
; INPUTS:
;	imname = BBSO VMG image file to map.  This image must be "squared";
;	         that is, it must have square pixels.
;	         The actual coordinates of Sun center (h and k in pixels)
;	         must be in the FITS header keywords CENTERX and CENTERY.
;	         The radius of the Sun in pixels must be in the keyword
;	         RADIUS.  The pixel coordinates are based on the lower
;	         left corner of the picture being (0,0).
; KEYWORD PARAMETERS:
; OUTPUTS:
;	A map with the name imnamexy.  It is 2880 x 360 in size.
; COMMON BLOCKS:
; NOTES:
;	The first 11 and last 8 columns in the VMGs are bad.
;	The procedure does not use them and subtracts 11 from the
;	h value before mapping.  Uses the procedure sunlong to
;	find P, B, and the Carrington longitude of the central meridian.
; MODIFICATION HISTORY:
;	J. Varsik,	19 Mar. 1996 --- Port from CMAP procedure.
;       J. Varsik,      26 Feb. 1999 --- Remove Fortran core.
;-
;-------------------------------------------------------------------
pro bigmap,imname,help=help

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

    ; Make output file names.
    dotpos = STRPOS(imname,".")
    iflen = dotpos
    ifin1 = STRMID(imname,0,iflen)
    exlen = STRLEN(imname) - iflen
    ext = STRMID(imname,dotpos,exlen)
    hname = ifin1 + 'xy' + ext

    ; Read input image
    im = FLOAT(READFITS(imname,imh,/NOSCALE))
    ; Find center, radius, p angle, b angle, etc.
    ew2 = LONG(FXPAR(imh,'CRVAL1'))
    ns2 = LONG(FXPAR(imh,'CRVAL2'))
    h = FLOAT(FXPAR(imh,'CENTERX'))
    k = FLOAT(FXPAR(imh,'CENTERY'))
    rad = FLOAT(FXPAR(imh,'RADIUS'))
    datec = FXPAR(imh,'DATE-OBS')
    timec = FXPAR(imh,'TIME-OBS')
;   Note new Y2K FITS format
    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
;    pandb,mon,day,yrs,hr,p,b,rsun,mjd

    ; Check parameters
    print, 'ew2 = ',ew2,' ns2 = ',ns2
    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.
    temp = im(11:503,*)
    ; temp = im(11:501,*)

    ; Make output image
    cart = fltarr(2880,360)
    s = rsun
    TRUE = 1
    FALSE = 0

    imgscl = rad/s
    ; which pole are we doing?
    thetar = ATAN((256-h),(192-k))
    IF (thetar LT 0.0) THEN thetar = (2.0 * !pi) + thetar
    theta = thetar * 57.2958
    br = b * 0.0174533
    pr = p * 0.0174533
    r1 = SQRT((256-h)^2 + (192-k)^2)
    r = ABS(ASIN(r1/rad) - (((r1*s/rad)/3600) * 0.0174533))
    latr = ASIN((SIN(br) * COS(r)) + (COS(br)*SIN(r)*COS(pr - thetar)))
    lat = latr * 57.2958
    IF (lat GT 30.0) THEN npole = TRUE
    IF (lat LT -30.0) THEN npole = FALSE
    IF (lat LE 30.0 AND lat GE -30.0) THEN BEGIN
        PRINT,'Not a polar image',lat
        RETURN
    ENDIF
    IF (npole) THEN PRINT,'North pole image' ELSE PRINT,'South pole image'
    PRINT,'image scale = ',imgscl

    ; 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 492 AND iy GE 0 AND iy LE 383) THEN BEGIN
                    cart(i,j) = temp(ix,iy)
                ENDIF ELSE BEGIN
                    cart(i,j) = -32767
                ENDELSE
          ENDELSE
       ENDFOR
    ENDFOR

    ; Place arrays into images, place information in FITS headers,
    ; write FITS files.
    FXHMAKE,imouth,cart
    FXWRITE,hname,imouth,cart
    
    finishup:

    return
end