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

  gamma = 5./3.
  pressure(U,rho) = (gamma-1.)*U*rho
  entropy(U,rho)  = pressure(U,rho)/rho**gamma

# Make sure this corresponds to your initial conditions.
  U1 = 1.; rho1 = 8000.
  U5 = 1.; rho5 = 1000.

  p1 = pressure(U1,rho1)
  p5 = pressure(U5,rho5)
  lambda = p1/p5
# This is the formula from the lecture.
  f1(P) = rho1/rho5/lambda*(1.-P)**2/(gamma*(1.+P)-1.+P)
  f2(P) = 2.*gamma/(gamma-1.)**2*(1.-(P/lambda)**((gamma-1.)/(2.*gamma)))**2

# Plot the two functions to find the intersection.
  set samples 1001
  set grid lt 7 lc rgb '#eeeeee'
  set xrange [0:4]
  set yrange [0:1]
  set xlabel 'P'
  set ylabel 'lhs, rhs'
  plot \
    f1(x) w l lt 1 t 'lhs', \
    f2(x) w l lt 3 t 'rhs'

# Zoom in to the intersection.
  set xrange [2.5:2.6]
  set xtics 0.01
  set auto y
  replot

# Zoom in again.
  set xrange [2.53:2.54]
  set xtics 0.001
  set auto y
  replot

# Read this number from the previous plot.
  P = 2.5329

# Analytic values.
  pc = P*p5
  rho4 = rho5*((gamma-1.)*p5+(gamma+1.)*pc)/((gamma-1.)*pc+(gamma+1.)*p5)
  rho3 = rho1*(pc/p1)**(1./gamma)
  vc = sqrt((pc-p5)*(1./rho5-1./rho4))
  vs = vc*rho4/(rho4-rho5)
  U4 = U5+0.5*vc**2*(pc+p5)/(pc-p5)
  U3 = pc/((gamma-1.)*rho3)
  S4 = entropy(U4,rho4)

  set samples 2
  set auto x
  set xtics auto
  unset key

  set multiplot
  set size 1,.20
  dy = 0.23
  y0 = 0.09
  set lmargin 10
  set tmargin 0
  set bmargin 0
  set arrow from 10.       ,graph 0 rto 0,graph 1 lt 7 dt 4 nohead front
  set arrow from 10.-1.0*vc,graph 0 rto 0,graph 1 lt 7 dt 2 nohead front
  set arrow from 10.-2.0*vc,graph 0 rto 0,graph 1 lt 7 dt 3 nohead front
  set arrow from 10.-1.0*vs,graph 0 rto 0,graph 1 lt 7 dt 2 nohead front
  set arrow from 10.-2.0*vs,graph 0 rto 0,graph 1 lt 7 dt 3 nohead front

  set ylabel 'density'
  set yrange [0:10000]
  set format x ''
  unset xlabel
  set origin 0,y0+3*dy
  plot \
    'output_000.txt' u 6:4 w d lt 1, \
    'output_010.txt' u 6:4 w d lt 2, \
    'output_020.txt' u 6:4 w d lt 3, \
    rho4 w l lt 7 dt 2, \
    rho3 w l lt 7 dt 2

  set ylabel 'pressure'
  set yrange [0:6000]
  set origin 0,y0+2*dy
  plot \
    'output_000.txt' u 6:(pressure($5,$4)) w d lt 1, \
    'output_010.txt' u 6:(pressure($5,$4)) w d lt 2, \
    'output_020.txt' u 6:(pressure($5,$4)) w d lt 3, \
    pc w l lt 7 dt 2

  set ylabel 'internal energy'
  set yrange [0:2]
  set origin 0,y0+1*dy
  plot \
    'output_000.txt' u 6:5 w d lt 1, \
    'output_010.txt' u 6:5 w d lt 2, \
    'output_020.txt' u 6:5 w d lt 3, \
    U4 w l lt 7 dt 2, \
    U3 w l lt 7 dt 2

  set ylabel 'entropy'
  set yrange [0:1e-2]
  set xlabel 'x'
  set format x '%g'
  set origin 0,y0+0*dy
  plot \
    'output_000.txt' u 6:(entropy($5,$4)) w d lt 1, \
    'output_010.txt' u 6:(entropy($5,$4)) w d lt 2, \
    'output_020.txt' u 6:(entropy($5,$4)) w d lt 3, \
    S4 w l lt 7 dt 2

  unset arrow
  unset multiplot