;------------------------------------------------------------------
;+
; NAME:
;	HGRID
; PURPOSE:
;	     This program generates an IDL 512x384 pixel short integer
;	image with a solar coordinate grid on it.  It uses information
;       from a squared-up BBSO magnetogram image header to determine where
;       the grid should be placed.
;
;	This IDL version is designed to run by itself, without needing
;	a Fortran "core".
;
;	This IDL function produces an IDL short integer image (512x384),
;	which contains the grid (as very bright pixels) and is otherwise
;	zero.
;
;	Routines used:
;		sun  -- from Johns Hopkins IDL library to find P, B, and L0.
; CATEGORY:
;	POLAR
; CALLING SEQUENCE:
;	outimage = HGRID,fheader
; INPUTS:
;	fheader= BBSO VMG image FITS header.
;	         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:
;	outimage = IDL short integer image containing the grid.
; COMMON BLOCKS:
; NOTES:
;	Uses the solar ephemeris from the Johns Hopkins library
;	find P, and B.
; MODIFICATION HISTORY:
;	J. Varsik,	09 Nov. 1995 --- Port from IRAF task
;	J. Varsik,	24 Feb. 1999 --- Numerous changes.
;	J. Varsik,	25 Feb. 1999 --- Remove Fortran "core".
;       J. Varsik,      03 May  1999 --- Allow both old and new date formats.
;-
;-------------------------------------------------------------------

function hgrid,imh,help=help

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

    ; Create output image
    outimage = INTARR(512,384)

    ; Read input image info.
    ; Find center, radius, p angle, b angle, etc.
    ew2 = FXPAR(imh,'CRVAL1')
    ns2 = FXPAR(imh,'CRVAL2')
    datec = FXPAR(imh,'DATE-OBS')
    timec = FXPAR(imh,'TIME-OBS')
    h = FLOAT(FXPAR(imh,'CENTERX'))
    k = FLOAT(FXPAR(imh,'CENTERY'))
    rad = FLOAT(FXPAR(imh,'RADIUS'))

;   Note new Y2K FITS format
    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

    max = 31000
    br = b * 0.0174533
    pr = p * 0.0174533

    FOR j = 0, 383 DO BEGIN
	FOR i = 0, 511 DO BEGIN
            xfc = (i+1)-h
	    yfc = (j+1)-k
	    r1 = sqrt((xfc)^2 + (yfc)^2)
	    IF (r1 LE rad) THEN BEGIN
		IF (r1 GE (rad - 1 *rad/rsun)) THEN outimage[i,j] = 30000
		IF ((r1 GT (rad - 4*rad/rsun)) AND (r1 LT (rad-3*rad/rsun))) $
		    THEN outimage[i,j] = 30000
		r = ABS(ASIN(r1/rad) - ((r1*rsun/rad)/3600)*0.0174533)
		thetar = ATAN((xfc),(yfc))
		IF (thetar LT 0.0) THEN thetar = 6.283185 + thetar
		theta = thetar * 57.2957795
		latr = ASIN(sin(br)*cos(r) + $
			cos(br)*sin(r)*cos(pr - thetar))
		lat = latr * 57.2957795
		x = FIX(4*lat) + 360 + 1
		IF (COS(latr) NE 0) THEN BEGIN
		    tst = SIN(r)*SIN((p-theta)*0.0174533)/COS(latr)
		    IF (ABS(tst) LE 1.0) THEN BEGIN
			lon = (ASIN(tst) * 57.2957795) + 1.0
		    ENDIF ELSE BEGIN
			PRINT,tst,theta,lat
			lon = 20.0
		    ENDELSE
		ENDIF
		IF (lon LE 0.0) THEN lon = lon + 360.0
		IF (lon GT 360.0) THEN lon = lon - 360.0
		y = FIX(4.0*lon) + 1
		IF ((x MOD 40) EQ 0) THEN outimage[i,j] = max + x/4
		IF ((y MOD 40) EQ 0) THEN outimage[i,j] = max + y/4
	    ENDIF
	ENDFOR
    ENDFOR

    finishup:

    RETURN, outimage
END