; $Id: map_set.pro,v 1.2 1993/06/11 15:37:04 dave Exp $

; Copyright (c) 1993, Research Systems, Inc.  All rights reserved.
;	Unauthorized reproduction prohibited.

pro mapp_set, latmin, latmax, lonmin, lonmax, p0lon, p0lat, rot, $
	iproj, limit, t3d=t3d, zvalue=zvalue, $
	border = border, title = title, sat_p = sat_p,    $
	xmargin = xmargin, ymargin = ymargin, position=position, color=color

!x.s = [0,1]		;normal coverage for uv plane
!y.s = [0,1]

if (n_elements(zvalue) eq 0) then zvalue = 0
if (n_elements(color)  eq 0) then begin
	if (!d.flags and 512) ne 0 then color = 0 else color = !d.n_colors-1
	endif

if iproj eq 3 then begin		;Conic?
	!map.phioc = rot		;Parameter order is diff
	!map.p0lat = (p0lat + p0lon)/2.	;Center latitude
	if !map.p0lat lt 0.0 then !map.sino = -1 else !map.sino = 1
	chi1 = (90. - !map.sino * p0lat) * !dtor
	if p0lat eq p0lon then begin
		!map.coso = cos(chi1)
	   endif else begin
		chi2 = (90. - !map.sino * p0lon) * !dtor
		!map.coso = alog(sin(chi1)/sin(chi2))/ $
				alog(tan(0.5*chi1)/tan(0.5*chi2))
	   endelse
        !map.out(8)=p0lat * !dtor	;Ucen
        !map.out(9)=p0lon * !dtor	;Vcen
	goto, get_bounds
endif

;		Not conic, all other projections

x1 = rot * !dtor
x2 = p0lat * !dtor

if (iproj ge 8) then  BEGIN
 !map.out(8) = x1		;ucen
 !map.out(11) = x2		;vrng
ENDIF ELSE !map.out(8)=x2

!map.out(9) = p0lon * !dtor		;vcen
sinr = sin(x1)			;Stash cos / sin of map center lat and rotation
cosr = cos(x1)
sino = sin(x2)
coso = cos(x2)
del1 = .0001
!map.sat = 0.0		;Clear satellite proj params

if iproj ge 8 and iproj ne 14 then begin   ;Cylindrical projections
  if (abs(p0lat) ge del1) or $	;Simple case?
	((abs(rot) ge del1) and (abs(abs(rot)-180.) ge del1))  THEN BEGIN
	sint = coso * cosr	;NO...  complicated case
	cost = sqrt(1.-sint*sint)
	tmp1 =  sinr / cost
	tmp2 =  sino / cost
	!map.phioc = p0lon - atan(tmp1, - cosr * tmp2) * !radeg
	sinr = tmp1 * coso
	cosr = -tmp2
	sino = sint
	coso = cost
  endif else begin		;Simple cases
	if abs(rot) lt del1 then sino = 1.0 else begin 
	   sino = -1.
	   !map.phioc = p0lon + 180.
	   endelse
	coso = 0.0
	sinr = 0.0
	cosr = 1.0
	!map.projection = iproj+3	;Make a simple projection
  endelse
endif			;Cylindrical projections

if iproj eq 7 then begin	;Special params for satellite projection.
	!map.sat(0) = sat_p(0)	;Salt = sat_p(0)
		;Save em.  sat(1) = TRUE for simple case (Vertical perspective)
	!map.sat(1) = (sat_p(1) eq 0.0) and (sat_p(2) eq 0.0)
	!map.sat(2) = sin(!dtor * sat_p(1))	;Salpha = sat_p(1)
	!map.sat(3) = cos(!dtor * sat_p(1))	;calpha
	!map.sat(4) = sin(!dtor * sat_p(2))	;Sbeta = sat_p(2)
	!map.sat(5) = cos(!dtor * sat_p(2))	;cbeta
	endif	

!map.coso =  coso
!map.sinr = sinr
!map.sino = sino
!map.cosr = cosr
if (iproj eq 2) then !map.sat(0) =0;


GET_BOUNDS:

if n_elements(limit) ne 8 then begin	;Explicit limits?
  ;	  stereo	ortho	     conic
  bounds = [[-2,-2,2,2], [-1,-1,1,1], [-180,-1,180,1], $
  ;			lamb	gnomic		azimuthal
	    [-2,-2,2,2],[-2,-2,2,2],[-!pi,-!pi,!pi,!pi]]
  ;			satell	cylindrical	mercator
  bounds = [[bounds], [-1,-1,1,1],[-180,-90,180,90],[-!pi,-!pi,!pi,!pi], $
  ;			mollweide  Sinusoidal            
			[-2,-1,2,1], [-!pi,-!pi/2,!pi,!pi/2], $
  ;                       aitoff
			[-2*sqrt(2), -sqrt(2), 2*sqrt(2), sqrt(2)]]
  
  iproj1 = (iproj < 11) - 1
  plimit = bounds(*,iproj1)  ;Fudge factor to avoid clipping
  if n_elements(limit) eq 0 then begin
      umin = plimit(2)      ;MAXIMUM USEFUL AREA, xmax
      vmin = plimit(3)	;ymax
      umax = plimit(0)	;xmin
      vmax = plimit(1)	;ymin
      endif

  latdel = (latmax-latmin)/10.
  if latdel gt 10. THEN BEGIN  ;We dont want stepsizes too large
  	del = float(fix(latmax-latmin)/10)
  	latdel  = (latmax-latmin)/del
  	ENDIF
   
  londel = (lonmax - lonmin)/10.
  if londel gt 10. THEN BEGIN
  	del = float(fix(lonmax - lonmin)/10)
  	londel  = (lonmax - lonmin)/del
  	ENDIF
  
  eps = 1.e-5 
  latlim = latmax + (latmax-latmin)*eps
  lonlim = lonmax + (lonmax-lonmin)*eps
  
  ;Find max u and v coordinates
  
  for lat = float(latmin), latlim, latdel/2.  do $  
      for lon = float(lonmin), lonlim, londel/2. do BEGIN
          p00 = convert_coord(lon,lat, /to_norm)
          p0  = !map.out(0:1)	; (p,q in uv coordinates)
  	  if (p0(0) ge plimit(0) and p0(0) le plimit(2) and $
  		p0(1) ge plimit(1) and p0(1) le plimit(3)) THEN BEGIN
	    if n_elements(umin) le 0 then begin
		umin = p0(0) & umax = p0(0)
		vmin = p0(1) & vmax = p0(1)
		endif
            umin = umin < p0(0)	;Smallest area
            umax = umax > p0(0)
            vmin = vmin < p0(1)
            vmax = vmax > p0(1)
            ENDIF
      ENDFOR

ENDIF ELSE BEGIN  			;4 point (8 element) limit specified

; limit =  [ lat0, lon0, lat1, lon1, lat2, lon2, lat3, lon3] specify
; four points on the map which give, respectively,
; the leftmost, topmost, rightmost, and bottom
; most points to be shown.

; print,'Limit = ', limit
; print, latmin, latmax, lonmin, lonmax
	p0 = fltarr(2,4)		;for the 4 corners
	for i=0,3 do begin
		lat = limit(2*i)
		lon = limit(2*i+1)
		p1 = convert_coord(lon, lat, /to_norm)
		str = strtrim(lat,2) + ',' + strtrim(lon,2)
		if max(!map.out(0:1)) ge 1e3 then message, $
			'Map_set, limit point not mappable, lat,lon='+str
		p0(0,i) = !map.out(0:1)		;Get uv coordinates
;		print, i, p0(*,i)
		endfor
	umin = p0(0,0)
	umax = p0(0,2)
	vmin = p0(1,1)
	vmax = p0(1,3)
ENDELSE


if iproj eq 3 then  BEGIN		;Set tolerance for detecting wraps.
    ueps = 135.				;Conic?
    umin = umin / 180.			;Scale to units of pi
    umax = umax / 180.
endif  else ueps = 0.75 * ABS(umax-umin)
veps = 0.75 * ABS(vmax-vmin)
  
!map.out(6) = ueps
!map.out(7) = veps

if n_elements(xmargin) ne 2 THEN xmargin = [1,1]
if n_elements(ymargin) ne 2 THEN ymargin = [1,2]
if not N_Elements(title) THEN title =" "

eps = (umax - umin) * .02		;Extend the axes by a Fudge factor
umax = umax + eps
umin = umin - eps
eps = (vmax - vmin) * .02
vmax = vmax + eps
vmin = vmin - eps

;
; This establishes the plot scaling between uv coordinates and the 
;  plot range.

if keyword_set(position) then  $
    plot, [umin,umax], [vmin,vmax], xsty=5, ysty=5, /nodata,    $
        xmargin= xmargin, ymargin= ymargin, title=title, font=0,/noerase,$
        position=position $
else plot, [umin,umax], [vmin,vmax], xsty=5, ysty=5, /nodata,    $
    xmargin= xmargin, ymargin= ymargin, title=title, font=0, /noerase

!x.type = 2			;Reset mapping type
if keyword_set(border) then  $
	plots ,!x.window([0,1,1,0,0]), color=color,   $
	       !y.window([0,0,1,1,0]), zvalue, /norm, /noclip, t3d=t3d

return
end



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 map_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

if (!x.type NE 2) THEN $  ; make sure we have mapping coordinates
   message, 'map_grid---Current ploting 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

First = 1

for lon = lonmin, lonmax, londel do begin
    if (lon lt -180) then lon2 =lon +360  $
    	else if (lon gt 180) then lon2 = -360 +lon $
    	else lon2 = lon
    pres = convert_coord(lon,latmin,/to_norm)
    pres = !map.out(0:1)
    pres1 = convert_coord(lon,latmax,/to_norm)
    pres1 = !map.out(0:1)
    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
    
    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

for lat = float(latmin), latmax, 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
   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



pro map_continents,  $
	mlinestyle = mlinestyle, color = color, t3d=t3d, $
	mlinethick = mlinethick, USA = usa, CONTINENTS = cont, zvalue=zvalue

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

if (!x.type NE 2) THEN message,'Map transform not established.'
latmin = !map.out(4)
latmax = !map.out(5)
lonmin = !map.out(2)
lonmax = !map.out(3)

if N_Elements(mlinestyle) EQ 0 THEN mlinestyle=0
if N_Elements(mlinethick) EQ 0 THEN mlinethick=1
if n_elements(t3d) le 0 then t3d = 0
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

maxlat=0. & minlat=0. & maxlon=0. & minlon=0.

openr, lun, /get, FILEPATH('supmap.dat',subdir = "maps"),/xdr,/stream
npts = 0L
; 	cont us_only  both
fbyte = [ 0, 71612L, 165096L]
nsegs = [ 283, 325, 594 ]


if n_elements(usa) eq 0 then usa = 0
if n_elements(cont) eq 0 then cont = 0

ij = 0			;Default = continents only
if usa ne 0 then begin
	if cont then ij = 2 else ij = 1
	endif
point_lun, lun, fbyte(ij)

for i=1,nsegs(ij) do begin
	READU, lun, npts,maxlat,minlat,maxlon,minlon
	npts = npts / 2		;# of points
	xy = fltarr(2,npts)
	READU, lun, xy
	if (maxlat lt latmin) or (minlat gt latmax) then goto,skipit
	
	if (maxlon lt lonmin) or (minlon gt lonmax) then BEGIN
          if (lonmax gt 180 and maxlon + 360 ge lonmin) then goto,goon
          if ( lonmin lt -180 and minlon -360 le lonmax) then goto,goon
          goto, skipit
        END
     

goon:	plots, xy(1,*), xy(0,*), zvalue, NOCLIP=0, t3d=t3d, $
		THICK=mlinethick,linestyle=mlinestyle, color = color

skipit:

       	endfor
FREE_LUN, lun
end



pro checkparam, latmin, latmax, lonmin, lonmax, p0lon,p0lat, limit ,  $
                iproj, Rot

NOLIMIT = 0
if (N_Elements(limit) NE 4) and (N_Elements(limit) ne 8) then BEGIN
  lonmin = -180.0 + p0lon
  lonmax  = 180.0 + p0lon
  if (iproj NE 3) THEN BEGIN
    latmax  = 90 
    latmin = -90 
  ENDIF
  NOLIMIT = 1
ENDIF

; If entire globe, center about P0lon
if (lonmin EQ -180 AND lonmax EQ 180) THEN BEGIN  
  if iproj ne 3 THEN BEGIN
    lonmin = -180 + P0lon
    lonmax = 180 + P0lon
  ENDIF ELSE BEGIN
    lonmin = -180 + Rot
    lonmax = 180 + rot
  ENDELSE
ENDIF

case 1 of
 iproj eq 1 or iproj eq 2 or iproj eq 4 :BEGIN ;stereo, ortho, lambert

        if (NOLIMIT eq 1) and p0lat eq 0.0 then BEGIN
           lonmax = p0lon + 90 	; center on equator
           lonmin = p0lon - 90
        ENDIF

        if (NOLIMIT eq 1) and p0lat eq 90.0   then  begin   ;center at a pole
            latmin = 0.0 
            latmax = 90.0
        ENDIF else if (NOLIMIT eq 1) and p0lat eq -90.0 then BEGIN
            latmin = -90.0
            latmax = 0.0
        ENDIF
    ENDIF			;Iproj eq 1 or 2 or 4
 iproj eq 5: BEGIN		;gnomic
  if (NOLIMIT eq 1) and p0lat eq 0.0 then BEGIN ; center on equator
           lonmax = p0lon + 60
           lonmin = p0lon - 60
           latmin = p0lat - 60
           latmax = p0lat + 60
  ENDIF
  if (NOLIMIT eq 1) and p0lat eq 90.0  then  begin
     latmin = 30
     latmax = 90.0
  ENDIF else if (NOLIMIT eq 1) and p0lat eq -90.0 then BEGIN
     latmin = -90.0
     latmax = -30.0
     ENDIF
  END

 iproj eq 3: BEGIN		;Conic
       if NOLIMIT eq 1 THEN BEGIN
	sgn = (p0lat ge 0) * 2 - 1	; + 1 or -1
	t0 = sgn * 15.
	t1 = sgn * 75.
        latmin = t0 < t1
        latmax = t0 > t1
        lonmin = -180 + Rot
        lonmax = 180  + Rot
        ENDIF		;Nolimit
      ENDCASE

 iproj eq 9 or iproj eq 12: if NOLIMIT eq 1 THEN BEGIN  ;mercator, merc/simple
		latmin = -80
		latmax = 80
		ENDIF
          
      ELSE:
ENDCASE
RETURN
END



pro map_set, p0lat, p0lon, rot, proj = proj, $
        Query = interg,                $
	CYLINDRICAL = cyl,      MERCATOR = merc,  $
	MOLLWEIDE = moll,       STEREOGRAPHIC = stereo,  $
        ORTHOGRAPHIC = ortho,   Conic = cone,            $
	LAMBERT = lamb,         GNOMIC = gnom,           $
        AZIMUTHAL = azim,       SATELLITE = satel,       $
        SINUSOIDAL = sinu,      AITOFF = aitoff,          $
  latdel = latdel, londel = londel, limit = limit,       $
  Sat_P = Sat_p, title = title, noborder = noborder,        $
  noerase = noerase, label = label, $
  glinestyle = glinestyle, glinethick=glinethick,       $
  mlinestyle=mlinestyle, mlinethick=mlinethick,		$
  color = color, con_color = con_color,                 $
  continents = continent, grid = grid, xmargin = xmargin, $
  ymargin = ymargin, lonlab = lonlab, latlab = latlab, $
  lonalign = lonalign, latalign = latalign, charsize = charsize, $
  advance = advance, usa=usa, t3d=t3d, position=position, zvalue=zvalue, $
  whole_map = whole_map

;+
; NAME:
;      map_set
; PURPOSE:
;        The procedure map_set establishes
;        the axis type and coordinate conversion mechanism
;        for mapping  points on the earth's surface, expressed
;        in latitude and longitude, to points on a plane, according
;        to one of ten possible map projections. The user may 
;        may select the map projection, the map center, polar rotation
;        and geographical limits.
;
;        The user may request map_set to plot the
;        graticule and continental boundaries by setting
;        the keywords Grid and Continent.
;
; CATEGORY:
;           Mapping
; CALLING SEQUENCE:
;                 map_set, p0lat, p0lon, rot, proj = proj,            $
;                    Query = interg,                                  $
;              	     CYLINDRICAL = cyl,      MERCATOR = merc,         $
;   	             MOLLWEIDE = moll,       STEREOGRAPHIC = stereo,  $
;                    ORTHOGRAPHIC = ortho,   Conic = cone,            $
;   	             LAMBERT = lamb,         GNOMIC = gnom,           $
;                    AZIMUTHAL = azim,       SATELLITE = satel,       $
;                    SINUSOIDAL = sinu,      AITOFF  = aitoff,         $
;                    latdel = latdel, londel = londel, limit = limit, $
;                    Sat_p = Sat_p, $
;                     title = title, noborder = noborder, $
;                    noerase = noerase, label = label,                $
;                    glinestyle = glinestyle, glinethick=glinethick,  $
;                    mlinestyle=mlinestyle, mlinethick=mlinethick,    $
;                    color = color, con_color = con_color, 	      $
;		     continents = continent, grid = grid, lonlab = lonlab,$
;                    latlab = latlab, lonalign = lonalign,           $
;                    latalign = latalign, charsize = charsize
;
; OPTIONAL INPUT:
;               P0lat--- For all but Lambert's conformal conic projection
;                        with two standard parallels, P0lat should be set
;                        to the latitude of the point on the earth's surface
;                        to be mapped to the center of the projection plane.
;                        If the projection type is sinusoidal, P0lat must be
;                        0. -90 <= P0lat <= 90. If P0lat is not set, the
;                        default value is 0. If the user has selected
;                        Lambert's conformal conic projection with two
;                        standard parallels, P0lat should be set to the
;                        latitude in degrees of one of the parallels.
;                        If not both P0lat and P0lon are defined, P0lat is
;                        20 by default.
;
;               P0lon--- For all but Lambert's conformal conic projection with
;                        two standard parallels, P0lon should be set
;                        by the user to the longitude of the point on the
;                        earth's surface to be mapped to the center of
;                        the map projection. -180 <= P0lon <= 180. If P0lon
;                        is not set, the default value is zero.
;
;              Rot-----  The user should set the argument Rot to the angle
;                        through which the North direction should be
;                        rotated around the line L between the earth's center
;                        and the point (P0lat, P0lon). Rot is measured in 
;                        degrees with the positive direction being clockwise
;                         rotation around L. Default value is 0.
;
; KEYWORDS:
;
;         Advance       = if set, advance to the next frame when there
;                         are multiple plots on the screen.
;         Azimuthal     = azimuthal equidistant projection
;         Conic         = Conic projection
;         Cylindrical   = cylindrcal equidistant projection
;         Gnomic        = gnomonic projection
;         Lambert       = Lambert's equal area projection
;         Mercator      = Mercator's projection
;         Mollweide     = Mollweide type projection
;         Orthographic  = Orthographic projection
;	  Satellite	= Satellite (General perspective) projection
;         Sinusoidal    = sinsoidal projection
;         Stereographic = stereographic projection
;         Aitoff        = Aitoff's projection
;         charsize      = size of characters in labels.
;         Color         = color of the map boundary.
;	  con_color	= color of continents.
;         Continents    = flag to signal that continental boundaries
;                         should be drawn.         
;         Glinestyle    = linestyle of gridlines. (default = 1 = dotted)
;         Glinethick    = thickness of gridlines. (default = 1 = normal)
;         Grid          = flag to signal that gridlines should be drawn.
;         Label         = flag to signal labeling of parallels and 
;                         meridians in the grid.
;
;         latalign     = the aligment of the text baseline for
;                        latitude labels. A value of 0.0 left justifies
;                        the label, 1.0 right justifies it and
;                         0.5 centers. The default value is .5.
;                         
;         Latdel        = if set, spacing of parallels drawn on grid. 
;			  Default is 10 <  (latmin - latmax) / 10).
;         LatLab       =  if set, longitude at which to place latitude
;                         labels. Default is longitude of center.
;         Limit         =  A four or eight element vector.
;			  If a four element vector, [latmin, lonmin, latmax,
;                          lonmax] specifying the boundaries of the region
;                          to be mapped. (latmin, lonmin) and (latmax,
;                          lonmax) are the latitudes and longitudes of
;                           two diagonal points on the boundary with
;                           latmin < latmax and lonmin <lonmax.
;			  If an eight element vector: [ lat0, lon0, 
;			    lat1, lon1, lat2, lon2, lat3, lon3] specify
;			    four points on the map which give, respectively,
;			    the location of a point on the left edge, 
;			    top edge, right edge, and bottom edge of
;			    the map extent.
;         lonalign     = the aligment of the text baseline for
;                        longitude labels. A value of 0.0 left justifies
;                        the label, 1.0 right justifies it and
;                         0.5 centers it. Default value is .5.
;        Londel         = if set, spacing of meridians drawn on grid.
;                         Y axis, and the screen border
;        LonLab       =  if set, latitude at which to place longitude
;                         labels. Default is longitude of center.
;        Mlinethick     = thickness of continental boundaries. Default is 1.
;        Mlinestyle     = linestyle of continental boundaries. Default is
;                         0 = normal.
;        Noborder       = if set, no border is drawn around the map.
;        Noerase        = flag to signal that map_set should not erase the
;                         current plot window. Default is to erase.
; 
;	Position	= Optional POSITION keyword for PLOT to specify the
;			  location of the map in the drawing area. See manual.
;	Sat_p		= an three element array of additional parameters
;			  required by the satellite projection.  The parameters
;			  are:  [ P, omega, gamma], where:  P = distance of
;			  viewpoint from center of globe, in units of the 
;			  globe's radius (P may be negative).
;			  Gamma = rotation of 
;			  projection plane clockwise from north.
;			  Omega = downward tilt of the 
;			  viewing plane in degrees.  The
;			  projection plane is first rotated gamma degrees 
;			  clockwise from north, and then tilted down
;			  omega degrees.  If both gamma and omega are zero,
;			  a Vertical Perspective projection results.
;	T3D 		= 1 to use the existing 3D transformation, 0 
;			  or omitted for normal.
;       Title           = Title of plot to be centered over plot window.
;       XMargin         = if set, a two element vector that specifies
;                         in character units the vertical margin between
;                         the map and the screen border.
;       YMargin         = if set, a two element vector that specifies
;                         in character units the horizontal margin between
;                         the map and the screen border.
;       Usa             = If set draw the state boundries
;	zvalue		= This keyword goes with T3D and specifies
;			  the height of the mapping plane.
;	whole_map	= in some cases where the limit is small, the gridlines
;			  will not cover the whole map.  This keyword can be
;			  specified to eliminate that problem.
;
; Modification history:
;			Written, Ann Bateson, 1991
;	SMR, 10-'91	Added a check for additional satelite parameters
;			   when /SATELITE is set.
;       UCSB, ESRG,5-92 Added USA flag to check and map state boundries
;	DMS,  6-92	Fixed bug that incorrectly printed grid values with
;				small lat/lon ranges.
;	DMS, 8-92	Added 8 element lat/lon keyword, fixed problems
;				with Satellite, conics with 1 std parallel,
;				smaller and/or wierd maps
;	DMS, 3-93	Extended axes by a fudge factor, improved labelling.
;	JIY, 4-93	Added ZVALUE keyword to work with T3D.
;			Added NO_GRID keyword to MAP_GRID for displaying
;			labels only.
;			Added WHOLE_MAP keyword to make the gridlines cover
;			the entire map area when the limit is small
;			COLOR now works with map border and labels.
;-
; !Map.out coorespondence:
;     0 p, 1 q, 2 umin, 3 umax, 4 vmin, 5 vmax, 6 ueps, 7 veps,
;      8 ucen, 9 vcen, 10 urng, 11 vrng;



if (N_ELements(interg) eq 0) THEN BEGIN                ; DONT INTERROGATE 

if n_elements(proj) ne 0 then iproj = proj $        
else if keyword_set(stereo) then iproj =1		$
else if keyword_set(ortho ) then iproj =2		$
else if keyword_set(cone  ) then iproj =3		$
else if keyword_set(lamb  ) then iproj =4		$
else if keyword_set(gnom  ) then iproj =5		$
else if keyword_set(azim  ) then iproj =6		$
else if keyword_set(satel ) then iproj =7		$
else if keyword_set(merc  ) then iproj =9		$
else if keyword_set(moll  ) then iproj =10		$
else if keyword_set(sinu  ) then iproj =14              $
else if keyword_set(aitoff) then iproj =15              $
else iproj=8		                         	;Assume cylindrical

IF (iproj eq 7) and (NOT(KEYWORD_SET(sat_p))) THEN $
    MESSAGE, "Satellite parameters must be set for satellite projection."

if keyword_set(noborder) eq 0 then border =1
if n_elements(limit) eq 4 then BEGIN
	latmin = limit(0)
	latmax = limit(2)
	lonmin = limit(1)
	lonmax = limit(3)                        
ENDIF ELSE if n_elements(limit) eq 8 then BEGIN
	latmin = min(limit([0,2,4,6]))
	latmax = max(limit([0,2,4,6]))
	lonmin = min(limit([1,3,5,7]))
	lonmax = max(limit([1,3,5,7]))
ENDIF

if !P.multi(0) eq 0 and keyword_set(advance) then erase

if n_elements(rot) eq 0 then rot = 0.	;Default rotation
if n_elements(t3d) le 0 then t3d = 0

if n_elements(p0lat) eq 0 then  $	;Defaults values for center of proj
  if iproj NE 3 THEN p0lat = 0.0 $
    else  BEGIN 
          p0lat = 20
          p0lon = 60
    ENDELSE

if n_elements(p0lon) eq 0 then   $
 if iproj NE 3 THEN p0lon = 0.0 $
 ELSE BEGIN
     p0lon = 60
     p0lat = 20
 ENDELSE
ENDIF

if (iproj eq 14) THEN BEGIN		;Sinusoidal
  if N_Elements(P0lat) NE 0 THEN $
    if P0lat NE 0 THEN  $
      message,'map_set--- Sinusoidal projection must have P0lat = 0'
  if N_Elements(Rot) NE 0 THEN  $
    if Rot NE 0 THEN   $
     message,'map_set--- Sinusoidal projection must have Rot = 0'
 p0lat = 0
 rot   = 0
endif

if N_Elements(interg) ne 0 THEN BEGIN
	iproj = 0
	read,'Projection number ',iproj
	read,'p0lat, p0lon, rot:  ',p0lat,p0lon,rot
	read,'latmin, latmax: ', latmin, latmax
	read,'lonmin, lonmax: ', lonmin, lonmax
ENDIF

if KEYWORD_SET(xmargin) THEN BEGIN
   csize = !d.X_SIZE/!d.X_CH_SIZE	;1/Char width in normalized units
   if total(xmargin) gt csize THEN BEGIN
      message,'map_set--- xmargin too large'
      return
   ENDIF
ENDIF

if KEYWORD_SET(ymargin) THEN BEGIN
   csize = !d.Y_SIZE/!d.Y_CH_SIZE
   if total(ymargin) gt csize THEN BEGIN
      message,'map_set--- ymargin too large'
      return
   ENDIF
ENDIF



!map.projection = iproj		;Set mapping parameters
!x.type = 2
!map.phioc = p0lon		;Save Center longitude /lat for internals
!map.p0lat = p0lat

if not N_Elements(noerase) and not KEYWORD_SET(advance) THEN erase

checkparam, latmin, latmax, lonmin, lonmax, p0lon, p0lat, limit, iproj, Rot
; print, latmin, latmax, lonmin, lonmax

!map.out(2) = lonmin		;Limits
!map.out(3) = lonmax
!map.out(4) = latmin
!map.out(5) = latmax


mapp_set, latmin, latmax, lonmin, lonmax, p0lon, p0lat,rot,iproj, limit,  $
      border=border, title =title, sat_p = sat_p, xmargin = xmargin, $
      ymargin = ymargin, t3d=t3d, position = position, zvalue=zvalue, $
      color=color

if keyword_set(grid) or keyword_set(label) then $		;Grid?
	map_grid, label=label, latdel=latdel, zvalue=zvalue, $
          londel = londel , glinestyle = glinestyle, $
          glinethick = glinethick, color = color, latlab = latlab,$
          lonlab = lonlab, lonalign = lonalign, latalign = latalign, $
          charsize = charsize, t3d=t3d, whole_map=whole_map, $
	  no_grid = (not keyword_set(grid))

if keyword_set(continent) then begin		;Continents?
	if n_elements(con_color) le 0 then begin
	    if n_elements(color) gt 0 then con_color = color $
		else begin
		  if (!d.flags and 512) ne 0 then color = 0 $
		    else color = !d.n_colors-1
		endelse
	    endif		;con_color
	    map_continents, color = con_color, t3d=t3d, $
		mlinestyle = mlinestyle, mlinethick = mlinethick, zvalue=zvalue
	endif

if keyword_set(usa) ne 0 then $			;US State outlines?
        map_continents, color = con_color, t3d=t3d, $
        mlinestyle = mlinestyle, mlinethick = mlinethick, us=usa, zvalue=zvalue

if KEYWORD_SET(ADVANCE) and !P.Multi(0) gt 0 THEN $
    !P.Multi(0) = !P.Multi(0) - 1 $
else !p.multi(0) = !p.multi(1) * !p.multi(2) - 1
       
RETURN
END