;------------------------------------------------------------------
;+
; NAME:
;	MHGRID
; PURPOSE:
;       This program generates an IDL 4000x2000 pixel short integer
;       image with a solar coordinate grid on it.
;
;	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,h,k,imgscl,fheader
; INPUTS:
;       h      = x-coordinate of center of Sun
;       k      = y-coordinage of center of Sun
;       imgscl = image scale (in pixels/arcsec).
;	fheader= BBSO VMG image FITS header (used to get date and time).
;       h and k are measured with (0,0) at the lower left corner of
;       the image.  Correct values are around (2000,0) for north pole
;       mosaics and (2000,2000) for south pole mosaics.
; 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,      11 Mar. 1999
;       J. Varsik,      06 May  1999 -- modify FITS date processing
;-
;-------------------------------------------------------------------

function mhgrid,h,k,imgscl,imh,help=help

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

    ; Create output image
    outimage = INTARR(4000,2000)

    ; Read input image info.
    ; 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
;    pandb,mon,day,yrs,hr,p,b,rsun,mjd

    rad = imgscl * rsun
    max = 31000
    br = b * 0.0174533
    pr = p * 0.0174533

    FOR j = 0, 1999 DO BEGIN
	FOR i = 0, 3999 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 *imgscl)) THEN outimage[i,j] = 30000
		IF ((r1 GT (rad - 4*imgscl)) AND (r1 LT (rad-3*imgscl))) $
		    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