;------------------------------------------------------------------
;+
; NAME:
;	PSLGRID
; PURPOSE:
;       This program generates an IDL 1540 x 360 pixel short integer
;       image with a solar coordinage grid on it.
;       The grid is based on the following parameters:
;       The solar radius is 1000 pixels.
;       Either the north pole or south pole is displayed.
;       The image is either very bright pixels (where the grid is)
;       or zero.
;
;       The grid is oriented with P = 0 (solar N/S meridian vertical).
;
;	Routines used:
;		sun  -- from Johns Hopkins IDL library to find P, B, and L0.
; CATEGORY:
;	POLAR
; CALLING SEQUENCE:
;	outimage = PSLGRID,fheader
; INPUTS:
;	fheader= BBSO VMG image FITS header (used to get date and time).
; 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,      28 May  2000
;-
;-------------------------------------------------------------------

function pslgrid,imh,help=help,south=south

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

    IF (KEYWORD_SET(south)) THEN southp = 1 ELSE southp = 0

    ; Create output image
    outimage = INTARR(1540,360)

    ; 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

    p = 0.0
    h = 769.0
    IF southp THEN BEGIN
        k = 1000.0
    ENDIF ELSE BEGIN
        k = -640.0
    ENDELSE
    rad = 1000.0
    imgscl = rad / rsun
    max = 31000
    br = b * 0.0174533
    pr = p * 0.0174533

    FOR j = 0, 359 DO BEGIN
	FOR i = 0, 1539 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