# 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.)