PRO dust_stamp, lon, lat, xsize, ysize, image, scale=scale_in, $
        equinox=equinox, grid=grid
;------------------------------------------------------------------------------
;+
; NAME:
;   dust_stamp
;
; PURPOSE:
;   Make postage-stamp images from BH files or our dust maps.
;
; CALLING SEQUENCE:
;   dust_stamp, lon, lat, xsize, ysize, image, scale=scale_in, $
;        equinox=equinox, grid=grid
;
; INPUTS:
;   lon, lat      longitude and latitude of image center 
;                   Galactic (l,b) assumed UNLESS equinox keyword is used
;   xsize, ysize  size of image to produce
;
;
; OPTIONAL INPUTS:
;   scale         pixel scale of desired image (in arcminutes per
;                   pixel)
;   equinox       equinox of equatorial coordinates (RA,dec)
;   grid          set keyword to overplot coordinate grid
;
;
; NOTE:  
;   All images are produced in a gnomic projection with north up
;     and east to the left.  (As viewed in IDL, or with SAOimage, etc...)
;   


  ; determine coordinate system
  IF keyword_set(equinox) THEN galcoords = 0 ELSE galcoords = 1
  IF galcoords THEN print, 'Galactic (l,b)' ELSE $
    print, 'RA, dec ('+string(equinox, format='(I4)')+')'


  ; default scale is 1.5 arcmin
  IF keyword_set(scale_in) THEN pixscale = scale_in ELSE pixscale = 1.5
  print, 'Pixel size: ', pixscale, ' arcmin'
  arcminscale = 3437.7468
  scale = arcminscale/pixscale


  xbox = lindgen(xsize, ysize) MOD long(xsize) - xsize/2.
  ybox = lindgen(xsize, ysize) / long(xsize) - ysize/2.
  xs = lindgen(xsize) - xsize/2.
  ys = lindgen(ysize) - ysize/2.


  IF galcoords THEN BEGIN 
      l0 = lon/!radeg
      b0 = lat/!radeg
      ; this procedure requires l,b in radians
      issa_proj_ignom, xbox, ybox, l0, b0, scale, l, b
      l = l*!radeg
      b = b*!radeg
  ENDIF ELSE BEGIN 
      ra0  = lon/!radeg
      dec0 = lat/!radeg
      ; this procedure requires l,b in radians
      issa_proj_ignom, xbox, ybox, ra0, dec0, scale, ra, dec
      ra  = ra*!radeg
      dec = dec*!radeg
      ra1  = reform(ra, long(xsize)*long(ysize))
      dec1 = reform(dec, long(xsize)*long(ysize))
      precess, ra1, dec1, equinox, 1950.
      euler, ra1, dec1, l1, b1, 1
      l = reform(l1, xsize, ysize)
      b = reform(b1, xsize, ysize)
  ENDELSE 

  image = dust_getval(l, b, map='I100', /noloop, /verbose)
  
  display,image,xs,ys, min=-3, max=80

; this is another way to draw coordinates
;  nlong = 24
;  nlat = 9
;  stamp_coord_grid, nlong, nlat, l0, b0, scale

  lspc = [1, 2, 4, 5, 10, 15, 30, 45, 60, 90]
  bspc = [1, 2, 5, 10, 15, 30]


  minb = mean(abs(b))
  imgsizedeg = pixscale*(xsize >ysize)/60.
  lind = (where((imgsizedeg/(4*cos(minb/!radeg)))>lspc EQ lspc))[0]
  IF lind EQ -1 THEN lind = n_elements(lspc)
  intl = lspc[lind]

  bind = (where((imgsizedeg/4)>bspc EQ bspc))[0]
  IF bind EQ -1 THEN bind = n_elements(bspc)
  intb = bspc[bind]

  IF galcoords THEN BEGIN   
      contour,l,xs,ys,/over,level=findgen(360/intl)*intl, c_charsize=1.5
      contour,b,xs,ys,/over,level=(findgen(180/intb-1)-90/intb+1)*intb, $
        c_charsize=1.5
  ENDIF ELSE BEGIN 
      contour,ra,xs,ys,/over,level=findgen(360/intl)*intl, c_charsize=1.5
      contour,dec,xs,ys,/over,level=(findgen(180/intb-1)-90/intb+1)*intb, $
        c_charsize=1.5
  ENDELSE 

  return
END


PRO stamp_coord_grid, nlong, nlat, l0, b0, scale

  x = findgen(361)/360.     ; use 360 to make 2x2 plots still look good.
  x1 = findgen(721)/360 -1.

  z = fltarr(361)
  z1 = fltarr(721)

  FOR i=1, nlong DO BEGIN
      ang = i*(360/nlong)
      IF (ang MOD 90) EQ 0 THEN BEGIN
          stamp_gridover, ang+z1, x1*90, l0, b0, scale
      ENDIF ELSE BEGIN 
          IF (ang MOD 30) EQ 0 THEN BEGIN
              stamp_gridover, ang+z1, x1*80, l0, b0, scale
          ENDIF ELSE BEGIN 
              stamp_gridover, ang+z1, x1*60, l0, b0, scale
          ENDELSE 
      ENDELSE

  ENDFOR

  FOR i=1-nlat, nlat-1 DO BEGIN 
      stamp_gridover, x*360, z+(90/nlat)*i, l0, b0, scale
  ENDFOR

  return
END


PRO stamp_gridover, ldeg, bdeg, l0, b0, scale
  lrad = ldeg/!radeg
  brad = bdeg/!radeg
  r0   = ll2uv([[l0], [b0]]*!radeg)
  rdeg = ll2uv([[ldeg], [bdeg]])
  dot = rdeg#transpose(r0)

  good = where(dot GE 0)

  IF good[0] NE -1 THEN BEGIN 
      issa_proj_gnom, lrad(good), brad(good), l0, b0, scale, xgrid, ygrid
      oplot, xgrid, ygrid
  ENDIF 
  return
END