pro drawcoordgrid,hdr, DEBUG=debug, SYSTEM=system, RSUN=rsun, DX=dx, DY=dy, COORDS=coords, units=units, position= position, $
 CARRINGTON=carrington, _EXTRA=_extra
;+
; $Id: drawcoordgrid.pro,v 1.15 2019/11/07 20:13:00 nathan Exp $
;
; Project   : STEREO SECCHI
;                   
; Name      : DrawCoordGrid
;               
; Purpose   : Draw contours on image corresponding to  a given projection
;               
; Explanation: 
;               
; Use       : IDL> im=scc_mk_image(filename,outsize=512,outhdr=h)
;   	      IDL> drawcoordgrid,h, system='hcr'
;    
; Inputs    : 
;   	hdr 	FITS header structure OR WCS structure
;
; Outputs   : None
;
; Keywords  : 
;   	Any keywords accepted by CONTOUR procedure
;   	SYSTEM= STRING designation of desired system, default is primary system of header
;   	    'HCR' - elongation, position angle (Heliocentric Radial)
;   	    'HAE' - RA-Dec (Helioprojective Ares Ecliptic)
;   	/RSUN	If system='hcr', show height units Rsun not elongation
;   	DX= 	set space between x-coord contours
;   	DY= 	set space between y-coord contours
;   	COORDS= Use this [2,nx,ny] array to draw contours instead of from header
;       units = units of coords
;   	/CARRINGTON 	Plot carrington coords
;               
; Calls from LASCO : NONE
;
; Common    : 
;               
; Restrictions: Assumes input header matches the image displayed in current window, 
;   	    	and the image starts at [0,0] in this window.
;               
; Side effects: Draws on current display window
;               
; Category    : Astrometry coordinates display wcs grid
;               
; Prev. Hist. : None.
;
; Written     : N.Rich, I2, 6/2008
;               
; $Log: drawcoordgrid.pro,v $
; Revision 1.15  2019/11/07 20:13:00  nathan
; add /carrington
;
; Revision 1.14  2011/09/26 18:33:22  mcnutt
; if coords defined no longer needs wcs and add keyword units for coords
;
; Revision 1.13  2011/07/25 22:16:46  nathan
; fix nlines
;
; Revision 1.12  2010/10/15 19:51:57  nathan
; auto-change to deg based on nlines not xran
;
; Revision 1.11  2010/07/15 15:12:51  nathan
; allow hdr input to be undefined
;
; Revision 1.10  2010/06/23 12:22:50  mcnutt
; returns dx and dy
;
; Revision 1.9  2010/01/20 13:59:38  mcnutt
; added keyword position for ab movies
;
; Revision 1.8  2009/03/12 18:26:46  nathan
; fix setting of units if input is WCS
;
; Revision 1.7  2008/12/19 22:48:47  nathan
; if value range exceeds 3600 change to degr; print units if /loud
;
; Revision 1.6  2008/11/14 14:42:26  mcnutt
; change x and y style to 5
;
; Revision 1.5  2008/07/14 18:44:15  nathan
; do ceil on params
;
; Revision 1.4  2008/06/25 18:24:07  nathan
; add COORDS= option
;
; Revision 1.3  2008/06/24 17:03:44  nathan
; fix some problems
;
; Revision 1.2  2008/06/24 15:25:32  nathan
; details...
;
; Revision 1.1  2008/06/24 15:18:19  nathan
; so far works on hi2-b
;
;
;-            
IF n_params() EQ 0 THEN BEGIN
    ;units='undefined' 
    IF ~keyword_set(COORDS) THEN BEGIN
    	message,'If no hdr input, COORDS= must be set. Returning.',/info
	return
    ENDIF
ENDIF 
IF keyword_set(COORDS) THEN BEGIN
    c=coords 
    if keyword_set(units)then units=units else units='degrees'
ENDIF ELSE BEGIN
    IF tag_exist(hdr,'CUNIT1') THEN units=hdr.cunit1 ELSE units=hdr.cunit[0]
    IF keyword_set(SYSTEM) THEN sys=strupcase(system) ELSE sys=''
    IF sys  EQ 'HAE' THEN wcs=fitshead2wcs(hdr,system='A') ELSE $
    IF tag_exist(hdr,'COORD_TYPE') THEN wcs=hdr ELSE $
    wcs=fitshead2wcs(hdr)

    IF sys EQ 'HCR' THEN BEGIN
	c=pix2coord(wcs) 
	IF keyword_set(RSUN) THEN BEGIN
    	    IF hdr.cunit1 EQ 'deg' THEN degfac=3600. ELSE degfac=1.
	    c[0,*,*]=c[0,*,*]*degfac/hdr.rsun
	    units='Rsun'
	ENDIF 
    ENDIF ELSE $
    c = wcs_get_coord(wcs)	
ENDELSE

; account for only lower left corner being displayed.
; If another subfield, then must be indicated in header.

IF keyword_set(CARRINGTON) THEN BEGIN
    wcs_convert_from_coord,wcs,c,'HG',lon,lat,/carrington
    c[0,*,*]=lon
    c[1,*,*]=lat
ENDIF

csz=size(c)
wx=!d.x_size<(csz[2])
wy=!d.y_size<(csz[3])
c=c[*,0:wx-1,0:wy-1]

; units are in arcsec, degrees, or rsun
sz=size(c)
; assuming 3-d array
nx=sz[2]
ny=sz[3]
cx= reform(c[0,*,*])
okx=where(finite(cx))
xext=minmax(cx[okx])
xran=xext[1]-xext[0]

nlines=nx/50 	; lines in image default
IF keyword_set(DX) THEN BEGIN
    intrv=dx
    nlines=ceil(xran/intrv)
ENDIF ELSE intrv=xran/nlines
IF nlines GT 1000 THEN BEGIN
; should convert from arcsec to degrees
;; If range exceeds 3000 arcsec, change to degrees
;IF xran GT 3000 THEN BEGIN
    cx=cx/3600
    c=c/3600.
    xran=xran/3600.
    xext=xext/3600.
    units='Degrees'
    nlines=ceil(xran/intrv)
ENDIF
IF keyword_set(DEBUG) THEN help,units,nlines

IF ~keyword_set(DX) THEN BEGIN
    intrv0=xran/nlines	
    IF intrv0 GT 20 THEN intrv=10*(fix(intrv0/10)) ELSE $
    IF intrv0 GT 15 THEN intrv=20 ELSE $
    IF intrv0 GT 7  THEN intrv=10 ELSE $
    IF intrv0 GT 3  THEN intrv=5  ELSE $
    IF intrv0 GT 1.5 THEN intrv=2 ELSE $
    IF intrv0 GT 0.5 THEN intrv=1 ELSE $
    IF intrv0 GT 0.3 THEN intrv=0.5 ELSE $
    IF intrv0 GT 0.15 THEN intrv=0.2 ELSE $
	intrv=0.1
ENDIF 

;'Draw Carrington Grid'
   if 'x' eq 'Carrington' then begin
       plevels=findgen(360/gridstr.grspace(0)+1)*gridstr.grspace(0)+(-180.)
       ; create labels for longitute 0-360 not (-180) - 180
       lonlevels=findgen(360/gridstr.grspace(1)+1)*gridstr.grspace(1)+(0.)
   endif

if keyword_set(position) ne 0 then xypos=position else xypos=[0,0,nx-1,ny-1]

x0=ceil(xext[0]/intrv)*intrv
nlines=fix(xran/intrv)

IF nlines LT 2 THEN BEGIN
    help,dx,dy
    message,'Grid spacing must be LT '+trim(xran/2.)+' deg.',/info
    return
ENDIF
contour,cx,/follow,levels=x0+indgen(nlines)*intrv,/noerase,position=xypos, /device, $
    xstyle=5, ystyle=5, ticklen=0, _EXTRA=_extra
IF keyword_set(DEBUG) THEN BEGIN
    print,'Done with x contour.'
    wait,2
ENDIF
DX=intrv

cy= reform(c[1,*,*])
oky=where(finite(cy))
yext=minmax(cy[oky])
yran=yext[1]-yext[0]
nlines=ny/50    ; lines in image default
IF keyword_set(DY) THEN BEGIN
    intrv=dy
    nlines=ceil(yran/intrv)
ENDIF ELSE BEGIN
    intrv0=yran/nlines	
    IF intrv0 GT 20 THEN intrv=10*(fix(intrv0/10)) ELSE $
    IF intrv0 GT 15 THEN intrv=20 ELSE $
    IF intrv0 GT 7  THEN intrv=10 ELSE $
    IF intrv0 GT 3  THEN intrv=5  ELSE $
    IF intrv0 GT 1.5 THEN intrv=2 ELSE $
    IF intrv0 GT 0.5 THEN intrv=1 ELSE $
    IF intrv0 GT 0.3 THEN intrv=0.5 ELSE $
    IF intrv0 GT 0.15 THEN intrv=0.2 ELSE $
	intrv=0.1
ENDELSE
y0=ceil(yext[0]/intrv)*intrv
nlines=fix(yran/intrv)
contour,cy,/follow,levels=y0+indgen(nlines)*intrv,/noerase,position=xypos,/device, $
    xstyle=5, ystyle=5, ticklen=0, _EXTRA=_extra
IF keyword_set(DEBUG) THEN wait,2
DY=intrv

END