# SM macros for contour plotting the density of discrete points.
  # Essentially a 2-d histogram.  Stellar astronomers refer to 
  # this as a "Hess diagram" when it is used to plot a color-magnitude
  # diagram of stars.    -bjw, 2/09/2006

  # How to use this:  save as file "docontour.sm" or whatever
  # to bin up and plot vectors mb and ub 
  #   (adjust the bin sizes and contour levels as needed)
  # macro read docontour.sm
  # set xbin = -25,-14,1
  # set ybin = -1.2, 0.7,0.1
  # set contlevels = 0,200,10
  # limits -14 -25 -1.2 0.7
  # box
  # plotbincontour xbin ybin mb ub contlevels
  
  # bin up 2 vectors into an SM image, essentially a 2-d histogram
  # Arguments: xbin ybin xdata ydata
  # can then be plotted with "contour"
  # xbin and ybin MUST be sorted increasing.
binupimage 4 #  xbins ybins xdata ydata
  set _xbin = $1
  set _ybin = $2
  set _xdata = $3
  set _ydata = $4
  define _nx $(dimen(_xbin))
  define _ny $(dimen(_ybin))
  define _nxm1 ($_nx-1)
  define _nym1 ($_ny-1)
  #define _xmin ((_xbin[1]+_xbin[0])/2)
  #define _xmax ((_xbin[$_nx-2]+_xbin[$_nx-1])/2)
  #define _ymin ((_ybin[1]+_ybin[0])/2)
  #define _ymax ((_ybin[$_ny-2]+_ybin[$_ny-1])/2)
  define _xmin (_xbin[0])
  define _xmax (_xbin[$_nx-1])
  define _ymin (_ybin[0])
  define _ymax (_ybin[$_ny-1])
  # define vectors that are the bin edges
  set dimen(_xedge) = $($_nx+1)
  set dimen(_yedge) = $($_ny+1)
  set _xedge[0] = _xbin[0] - (_xbin[1]-_xbin[0])/2
  set _yedge[0] = _ybin[0] - (_ybin[1]-_ybin[0])/2
  set _xedge[$_nx] = _xbin[$_nx-1] + (_xbin[$_nx-1]-_xbin[$_nx-2])/2
  set _yedge[$_ny] = _ybin[$_ny-1] + (_ybin[$_ny-1]-_ybin[$_ny-2])/2
  do _i=1,$_nx-1 {
    set _xedge[$_i] = (_xbin[$_i-1] + _xbin[$_i])/2
  }
  do _j=1,$_ny-1 {
    set _yedge[$_j] = (_ybin[$_j-1] + _ybin[$_j])/2
  }
  # clip the data to the edges, not the bin centers
  set _ifclip = (_xdata>_xedge[0] && _xdata<_xedge[$_nx] && \
    _ydata>_yedge[0] && _ydata<_yedge[$_ny]) ? 1 : 0
  set _xdata = _xdata if (_ifclip==1)
  set _ydata = _ydata if (_ifclip==1)
  image ( $_nx, $_ny ) $_xmin $_xmax $_ymin $_ymax
  do _i=0,$_nx-1 {
    set _dcol = _ydata if (_xdata>_xedge[$_i] && _xdata<_xedge[$_i+1])
    # the histogram is on the bin centers
    set _npcol = histogram(_dcol:_ybin)
    # This seems to actually work
    set image[$_i,*] = _npcol
    # do-loop way is unnecessary
    #do _j=0,$_nym1 {
    #  set image [$_i,$_j] = _npcol[$_j]
    #}
  }
  # use this if you want to save the image, it will only work
  # if you have file_type defined in your .sm file.
  #write image binimage.sm_img

  # to set the contour levels
  # set junk = 0,100,5
  # levels junk
  # limits _xbin _ybin
  # box
  # contour

  # plot density of points as contours
  # Arguments: xbin ybin xdata ydata contour-levels
plotbincontour 5
  binupimage $1 $2 $3 $4
  levels $5
  contour

  # plot density of points as greyscale
  # I haven't yet found sm greyscale settings that make it look really good.
  # Arguments: xbin ybin xdata ydata contour-levels
plotbingreyscale 5
  binupimage $1 $2 $3 $4
  glevels $5
  greyscale


 # Copyright 2006, Benjamin Weiner, licensed under GPL.
 # You may use or redistribute this code freely, please retain this
 # notice.  If you incorporate it into a package, you can use the
 # package freely for your own research.  (You can't sell it without
 # distributing the source freely.)