set term pngcairo size 1536,768 background rgb '#ffffee'
# set output 'blast.png'
  set output '| ffmpeg -y -r 15 -f image2pipe -i - -vcodec libx264 -profile:v baseline -pix_fmt yuv420p blast.mp4'
  set colors classic
  set encoding utf8
  set minussign

# Box edges are 0...2, the blast center is at (1,1,1).
  r(x,y,z) = sqrt((x-1.)**2+(y-1.)**2+(z-1.)**2)
  TimeBetSnapshot = 0.1
  gamma = 5./3.
  p(eps,rho) = (gamma-1.)*eps*rho
  rho0 = 1e6/8.

# These are the scaling relations from the lecture.
# Proportionality coefficients are roughly fitted "by eye".
  rhomax = (gamma+1.)/(gamma-1.)*rho0
  Rscal(t) = 0.216*t**0.4
  pscal(t) = t==0?NaN:481.*t**-1.2
  vscal(t) = t==0?NaN:0.4*0.216*t**-0.6*2./(gamma+1.)
  escal(t) = pscal(t)/(gamma-1.)/rhomax

  unset key
  set tics front
  set object rectangle from graph 0,0 to graph 1,1 fs solid noborder fc rgb 'white'
  set style line 1 lt 2 lc rgb '#00cc00'
  set style line 2 lt 1 lc rgb '#ff0000'

  set palette model HSV functions 5./6.*gray, 1, 1
  set colorbox horizontal user size graph 1,0.01 origin graph 0,1.015
  set cbtics scale 0.5 offset 0,2.4
  set cblabel 'density' offset 0,5.8
  set cbrange [0:5e5]

# We set all margins to zero and position the plots manually instead.
  set tmargin 0
  set bmargin 0
  set lmargin 0
  set rmargin 0
  ar = 768./1536.
  wd = 0.39*ar
  ht = 0.39
  sx = 0.02*ar
  sy = 0.02
  x0 = 0.14*ar
  y0 = 0.1
  x1 = x0+wd+sx
  y1 = y0+ht+sy

  n = 101
  do for [i=0:n] {
    filename=sprintf('output_%03d.txt',i)
    print filename
    t = i*TimeBetSnapshot
    lbl = sprintf('t = %4.1f',t)
    set multiplot
    set size   wd,ht
    set origin x0,y1
    unset y2tics
    unset y2label
    set ytics mirror
    unset log y
    set yrange [0:8e5]
    set ylabel 'density'
    set xrange [0:1]
    unset xlabel
    set format x ''
    set label lbl at screen 0,1 left offset 3,-1.5 font 'Helvetica-Bold'
    set arrow 1 from first Rscal(t),graph 0 to first Rscal(t),graph 1 nohead ls 2
#   set arrow 2 from graph 0,first rhomax   to graph 1,first rhomax   nohead ls 2
    plot \
      filename u (r($6,$7,$8)):4 w d ls 1
    unset label

    set origin x0,y0
    set log y
    set yrange [1e-21:10]
    set ylabel 'internal energy'
    set format x '%g'
    set xlabel 'radius'
    set arrow 2 from graph 0,first escal(t) to graph 1,first escal(t) nohead ls 2
    plot \
      filename u (r($6,$7,$8)):5 w d ls 1

    set origin x1,y1
    unset ytics
    unset ylabel
    set y2tics mirror
    set y2range [0:0.5]
    set y2label 'velocity'
    unset xlabel
    set format x ''
    set arrow 2 from graph 0,second vscal(t) to graph 1,second vscal(t) nohead ls 2
    plot \
      filename u (r($6,$7,$8)):(sqrt($9**2+$10**2+$11**2)) w d ls 1 axes x1y2

    set origin x1,y0
    set y2tics mirror
    set log y2
    set y2range [1e-16:1e6]
    set y2label 'pressure'
    set format x '%g'
    set xlabel 'radius'
    set arrow 2 from graph 0,second pscal(t) to graph 1,second pscal(t) nohead ls 2
    plot \
      filename u (r($6,$7,$8)):(p($5,$4)) w d ls 1 axes x1y2

    unset arrow
    set origin 0.55,y0
    set size wd*2.,ht*2.
    set y2tics mirror
    unset log y2
    set y2range [0:2]
    set y2label 'y'
    set xrange [0:2]
    set format x '%g'
    set xlabel 'x'
    plot \
      filename u (abs($8-1.)<0.1?$6:NaN):7:4 w d axes x1y2 lc palette

    unset multiplot
  }