set term post enh color solid
  set output 'histogram.ps'
  set encoding iso_8859_1
  set colors classic

  rstar = 1.
  mstar = 1.5815
  rhomean = mstar/(4.*pi/3.*rstar**3)
# a truncated 1/r density profile should have rho(R)=2/3*rhomean:
  rho1 = 2./3.*rhomean
# for comparison with our histogram, which is in units of mstar:
  f = rho1/mstar

  set grid lt 6 lc rgb '#eeeeee' xtics mxtics ytics mytics
  set samples 1001
  unset key
  set xlabel 'radius'
  set ylabel 'density'

  set xrange [0:1.1]
  set yrange [0:18]
  plot \
    'histogram.txt' u 3:($6==0?1/0:$6):7 w errorbars lt 2 pt 7, \
    'histogram.txt' u ($4==0?1e-3:$4):6 w steps lt 1, \
    f/x w l lt 3

  set yrange [0:2]
  replot

  set yrange [0.1:30]
  set xrange [5e-3:2]
  set log x
  set log y
  plot \
    'histogram.txt' u 3:6:7 w errorbars lt 2 pt 7, \
    'histogram.txt' u ($4==0?1e-3:$4):6 w steps lt 1, \
    f/x w l lt 3