pro rll_grid, dlong, dlat, label = label, charsize = charsize, shift = $
              shift, _extra = _extra

;+
; NAME:
;	RLL_GRID
;
;
; PURPOSE:
;	Add a latitude/longitude grid to a rectangular image
;
;
; CATEGORY:
;	Utils
;
;
; CALLING SEQUENCE:
;	rll_grid[, dlong, dlat]
;
;
; OPTIONAL INPUTS:
;	dlong	int	The step in longitude
;	dlat	int	The step in latitude
;
;
; KEYWORD PARAMETERS:
;	/label			If set, then label the grid.
;	charsize	Float	Set a size for the labels (ignored if
;				label is not set.
;	max_elongation	float	Set the largest elongation to show
;				(default=135)
;	shift		float	The shift in degrees for the display.
;	Any key accepted by PLOTS will be passed through.
;
;
; RESTRICTIONS:
;	Steps will be rounded to the nearest whole degree
;
;
; MODIFICATION HISTORY:
;	Original (after fll_grid): 9/7/08; SJT
;	Support shifting (and make it work): 23/11/09; SJT
;-

if n_params() ne 2 then dlat = 30 else dlat = round(dlat)
if n_params() eq 0 then dlong = 30 else dlong = round(dlong)
if n_elements(shift) eq 0 then shift = 0.

lng = findgen(361)-180
for lat = -90, 90, dlat do begin
    elon = acos(cos(lng*!dtor)*cos(lat*!dtor))*!radeg
    pat = 90-asin(sin(lat*!dtor)/sin(elon*!dtor))*!radeg
    locs = where(lng lt 0)
    pat[locs] = -pat[locs]
    pat = mod_pm(pat-shift, 180.)
    lbr = [-1, where(abs(pat-pat[1:*])+1 gt 90., nseg), $
           n_elements(pat)-1]
    if nseg eq 0 then lbr = [-1, n_elements(pat)-1]
    for j = 0, nseg do $
      plots, pat[lbr[j]+1:lbr[j+1]], elon[lbr[j]+1:lbr[j+1]], _extra = _extra
endfor

lat = findgen(181)-90
for lng = -180, 180, dlong do begin
    if abs(lng) mod 180 eq 0 then continue
    elon = acos(cos(lng*!dtor)*cos(lat*!dtor))*!radeg
    pat = 90-asin(sin(lat*!dtor)/sin(elon*!dtor))*!radeg
    if lng lt 0 then pat = -pat
    pat = mod_pm(pat-shift, 180.)
    lbr = [-1, where(abs(pat-pat[1:*])+1 gt 90., nseg), $
           n_elements(pat)-1]
    if nseg eq 0 then lbr = [-1, n_elements(pat)-1]
    for j = 0, nseg do $
      plots, pat[lbr[j]+1:lbr[j+1]], elon[lbr[j]+1:lbr[j+1]], _extra = _extra
endfor

elon = [0., 180.]
pat = mod_pm(replicate(-shift, 2), 180.)
plots, pat, elon, _extra = _extra
pat = mod_pm(replicate(180.-shift, 2), 180.)
plots, pat, elon, _extra = _extra

if keyword_set(label) then begin
     dy = 0.5
     dx = 0.5
     lat = 0.0
     for lng = -180+dlong, 180, dlong do begin
         y = acos(cos(lng*!dtor)*cos(lat*!dtor))*!radeg
         x = 90-asin(sin(lat*!dtor)/sin(lng*!dtor))*!radeg
        
         if lng lt 0 then x = -x
         x = mod_pm(x-shift, 180.)
         if x eq -180. then begin
             x = x+dx
             algn = 1.
         endif else begin
             x = x-dx
             algn = 0.
         endelse
         xyouts, x, y+dy, string(lng, format = "(I0)"), $
           align = algn, charsize = charsize
     endfor

     lng = -90.0
     dx = 0.5

     for lat = -90, 90, dlat do begin
         y = acos(cos(lng*!dtor)*cos(lat*!dtor))*!radeg
         x = 90-asin(-1. > sin(-lat*!dtor)/sin(lng*!dtor) < 1.)*!radeg
         if lng lt 0 then x = -x
         x = mod_pm(x-shift, 180)
         if x eq -180. then begin
             x = x+dx
             algn = 1.
         endif else begin
             x = x-dx
             algn = 0.
         endelse
         xyouts, x, y, string(lat, format = "(I0)"), $
           align = algn, charsize = charsize
    endfor
endif

end