IDL

PRO plotc,  x0,  y0,  ctable = ctable, every = every, buffer = buffer, $
            overplot = overplot, _extra = eee

compile_opt idl2

IF ~keyword_set(ctable) THEN ctable = 33

IF keyword_set(every) THEN BEGIN
   message,  "Every is not working right!", /inf
   n = n_elements(x0)
   w = indgen(floor(1.*n/every))*every < n-1
   nw = n_elements(w) 
   x = x0[w]
   y = y0[w]
   message, /inf, "Plotting "+n2s(100.*nw/n )+"% of the points. That's "+n2s(nw)+" points" 
ENDIF ELSE BEGIN
   x = x0
   y = y0
ENDELSE  

h2=sshist_2d(x,y, re=ri1, cost=co,  outbin = bin)
;rad = 20
;d= 255b-SHIFT(DIST(2*rad+1), rad,rad)
;h2 = convol(h2,d, /edge_w)
;h2 = smooth(h2,rad, /edg)
xmin = min(X) &  ymin = min(y)
h2size = size(h2, /dimen)
col = h2[ floor((x-xmin)/bin[0]) + floor((y-ymin)/bin[1])*h2size[0] ]

cgloadct, ctable

IF keyword_set(buffer) THEN BEGIN
   set_plot, 'z'
   device, z_buff=0, set_res=[!D.X_SIZE,!D.Y_SIZE]
ENDIF 

col = bytscl(col)

IF ~keyword_set(overplot) THEN cgplot, x, y, /noda, back=cgcolor('black'), $
color=cgcolor('white'), chars=1.5,  _extra = eee



cmin = min(fix(col),  max = cmax)

for c=cmax, cmin, -1 do begin
   w=where(col EQ c)
   if w[0] ne -1 THEN BEGIN 
      cgplot, x[w], y[w], syms=.1,  /ov, col=c, _extra=eee
   ENDIF 
ENDFOR 

IF keyword_set(buffer) THEN BEGIN
   a=tvrd(/tr)
   set_plot, 'x'
   tv, 255b-a, /tr
ENDIF 
;plot_hist,  col,  /l

RETURN
END