; These routines were extracted from the User Library routine
; MAP_SET.PRO and subsequently modified by D P Steele, 950302.

function map_grid_incr, span
ipow = 0
t = span
while t lt 5 do begin t = t * 10 & ipow = ipow +1 & endwhile
increments = [ 1., 2., 4., 5., 10., 15., 30., 45.]
i = 0
while t gt (increments(i) * 10) do i = i + 1
return, increments(i) / 10^ipow
end


; Put a grid on a previously established map projection.
pro mmap_grid, label=label,  latdel = latdel, no_grid=no_grid, zvalue=zvalue,$
       londel=londel, glinestyle = glinestyle, glinethick = glinethick,   $
       color = color, lonlab = lonlab, latlab = latlab, lonalign = lonalign,$
       latalign = latalign, charsize = charsize, t3d=t3d, whole_map=whole_map,$
       nice=nice

if (!x.type NE 2) THEN $  ; make sure we have mapping coordinates
   message, 'Current plotting device must have mapping coordinates'

; no grid? - in case someone wants just to put labels
if not keyword_set(no_grid) then no_grid = 0

if (n_elements(zvalue) eq 0) then zvalue = 0

lonmin = !map.out(2)		;Get lat/lon ranges from !MAP
lonmax = !map.out(3)

latmin = !map.out(4)
latmax = !map.out(5)

if n_elements(t3d) le 0 then t3d = 0

; if WHOLE_MAP specified, do not convert the limit into integers
if (not keyword_set(whole_map)) then begin
   if abs(latmax - latmin) gt 4. then begin	;Make range into integers
	latmin = floor(latmin)
	latmax = ceil(latmax)
	endif
   if abs(lonmax - lonmin) gt 4 then begin
	lonmin = floor(lonmin)
	lonmax = ceil(lonmax)
	endif
endif

			;Default grid spacings...
if n_elements(latdel) eq 0 then latdel = map_grid_incr(latmax-latmin)
if n_elements(londel) eq 0 then londel = map_grid_incr(lonmax-lonmin)
;; Old way:	londel = (lonmax-lonmin)/10. < 15 

if N_Elements(glinestyle) EQ 0 THEN glinestyle =1
if N_Elements(glinethick) EQ 0 THEN glinethick =1

if n_elements(color) le 0 then begin	;Default color?
	if (!d.flags and 512) ne 0 then color = 0 else color = !d.n_colors-1
	endif


if N_Elements(label) NE 0 OR (N_ELEMENTS(Latlab) ne 0) $
	OR (N_Elements(LonLab) NE 0) THEN BEGIN
	printno = 1  
	printno2 = 1
	if N_Elements(Latlab) eq 0 THEN Latlab = (lonmin + lonmax)/2
	if N_ELements(LonLab) eq 0 THEN LonLab = (latmin +latmax)/2
endif ELSE BEGIN
	printno = -1
	printno2 = -1
ENDELSE
                                    ; of grid numbers
if n_elements(latalign) eq 0 THEN latalign = .5	;Text alignment of lat labels
if n_elements(lonalign) eq 0 THEN lonalign = .5 ;Text alignment of lon labels
if n_elements(charsize) eq 0 THEN charsize =  1

step = 4 < (latmax - latmin)/10.
len = long((latmax-latmin) / step + 1)
lati = ((latmax-latmin) / (len-1)) * findgen(len) + latmin   ;Array of lats
;MESSAGE,'lati='+strcat(lati,/trim),/INFORMATIONAL

First = 1

; Following lines added by D P Steele, 950302
lonstart=lonmin
lonend=lonmax
IF N_ELEMENTS(nice) EQ 1 THEN BEGIN
    lonstart=CEIL(FLOAT(lonstart)/londel)*londel
    lonend=FLOOR(FLOAT(lonend)/londel)*londel
    ENDIF

;for lon = lonmin, lonmax, londel do begin
for lon = lonstart, lonend, londel do begin
    if (lon lt -180) then lon2 =lon +360  $
    	else if (lon gt 180) then lon2 = -360 +lon $
    	else lon2 = lon
    ; MESSAGE,'lon='+STRTRIM(lon,2)+' lon2='+STRTRIM(lon2,2),/INFORMATIONAL
    pres = convert_coord(lon,latmin,/to_norm)
    ; MESSAGE,'pres='+strcat(pres,/trim),/INFORMATIONAL
    pres = !map.out(0:1)
    ; MESSAGE,'pres='+strcat(pres,/trim),/INFORMATIONAL
    pres1 = convert_coord(lon,latmax,/to_norm)
    ; MESSAGE,'pres1='+strcat(pres1,/trim),/INFORMATIONAL
    pres1 = !map.out(0:1)
    ; MESSAGE,'pres1='+strcat(pres1,/trim),/INFORMATIONAL
    lon1 = lon
    if First eq 1  THEN First = 0 else $
      if abs(pres(0) - past(0)) GE !map.out(6)   OR  $
         abs(pres(1) - past(1)) GE !map.out(7)   OR  $
         abs(pres1(0) - past1(0)) GE !map.out(6) OR  $
         abs(pres1(1) - past1(1)) GE !map.out(7)     $
      THEN BEGIN
       if(lon ge 0) then dd = .0001 else dd = -.0001
       lon1 = lon - dd
      ENDIF  
    past = pres
    past1 = pres1
    
    ; MESSAGE,'lon1='+STRTRIM(lon1,2),/INFORMATIONAL
    if (not no_grid) then plots, Replicate(lon1,len), lati, zvalue, $
	color = color, t3d=t3d, NOCLIP=0,linestyle=glinestyle,thick=glinethick

    if lon2 ne long(lon2) then fmt = '(f7.2)' else fmt = '(i4)'
    if (printno eq 1) and $  ;Dont repeat -180....
	((lonmin ne -180) or (lonmax ne 180) or (lon ne -180)) then $
	xyouts,lon, LonLab, z=zvalue, ali=lonalign, t3d=t3d, color=color,$
          strtrim(string(lon2,format=fmt),2), charsize = charsize
    printno  = 1 - printno
 endfor

  !map.out(3) = lonmax 


step = 4 < (lonmax - lonmin)/10.
len = (lonmax-lonmin)/step + 1
loni = findgen(len)*step + lonmin


if (loni(len-1) NE lonmax) THEN  BEGIN
    loni = [loni, lonmax]
    len = len + 1
ENDIF
;MESSAGE,'loni='+strcat(loni,/trim),/INFORMATIONAL

; Lines below added by D P Steele, 950302
latstart=latmin
latend=latmax
IF N_ELEMENTS(nice) GT 0 THEN BEGIN
    latstart=CEIL(FLOAT(latstart)/latdel)*latdel
    latend=FLOOR(FLOAT(latend)/latdel)*latdel
    ENDIF
;MESSAGE,'latstart='+STRTRIM(latstart,2)+' latend='+STRTRIM(latend,2),/INFORMATIONAL

;for lat = float(latmin), latmax, latdel do begin
for lat = FLOAT(latstart), latend, latdel do begin
   if lat ne long(lat) then fmt = '(f7.2)' else fmt = '(i4)'
   if printno2 eq 1 then xyouts,latlab,lat, z=zvalue, ali=latalign, t3d=t3d, $
	strtrim(string(lat, format=fmt),2), charsize = charsize, color=color
   printno2 = 1 - printno2
   ;MESSAGE,'lat='+STRTRIM(lat,2),/INFORMATIONAL
   if (not no_grid) then plots,loni, Replicate(lat,len), zvalue, $
	NOCLIP=0,linestyle=glinestyle,color = color, thick=glinethick, t3d=t3d
 endfor
    !map.out(5) = latmax

return
end