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

  G       = 1.0
  Rstar   = 1.0
  Mstar   = 1.5815
  rhostar = 3.*Mstar/(4.*pi*Rstar**3)
  epsstar = G*Mstar/Rstar
  pstar   = rhostar*epsstar
  vstar   = sqrt(G*Mstar/Rstar)

  gamma = 5./3.
  p(U,rho) = (gamma-1.)*U*rho
  vrad(x,y,z,vx,vy,vz) = (vx*x+vy*y+vz*z)/sqrt(x**2+y**2+z**2)
  r(x,y,z) = sqrt(x**2+y**2+z**2)

  set xrange [0.005:3]
  set log x

  unset key
  set lmargin 10
  set tmargin 0
  set bmargin 0
  dy = 0.30
  y0 = 0.07

  do for [i=1:24] {
    ppmfile  = sprintf('PPM_%02d',i)
    snapfile = sprintf('output_%03d.txt',i)
    set multiplot
    set size 1,0.28
    unset xlabel
    set format x ''

    set origin 0,y0+dy*2
    set ylabel 'density' offset -1,0
    set log y
    set yrange [3e-3:3e3]
    set format y '10^{%T}'
    set title sprintf('snapshot %03d',i) offset 0,-0.5
    plot \
      snapfile u (r($6,$7,$8)):($4/rhostar) w p lt 1 pt 7 ps .5, \
      ppmfile every ::1 u 1:2 w l lt 3

    set origin 0,y0+dy*1
    unset title
    set ylabel 'pressure' offset -1,0
    set yrange [3e-5:3e3]
    plot \
      snapfile u (r($6,$7,$8)):(p($5,$4)/pstar) w p lt 1 pt 7 ps .5, \
      ppmfile every ::1 u 1:3 w l lt 3

    set origin 0,y0+dy*0
    set ylabel 'velocity'
    unset log y
    set yrange [-3:7]
    set format y '%g'
    set xlabel 'radius' offset 0,0.5
    set format x '%g'
    plot \
      snapfile u (r($6,$7,$8)):(vrad($6,$7,$8,$9,$10,$11)/vstar) w p lt 1 pt 7 ps .5, \
      ppmfile every ::1 u 1:4 w l lt 3

    unset multiplot
  }