;
; === Grid-based 1D hydro example code written in IDL ===
;
; written by
; Christoph Federrath - 2012-2020
;
forward_function eos, check_conservation, tickformat_exponential


pro advection_test, FLUXLIM=fluxlim, SCHEME=scheme, CFL=cfl

  if not keyword_set(fluxlim) then fluxlim = 'donor-cell'
  if not keyword_set(scheme) then scheme = 'upwind'

  ; basic user parameters
  nghost = 2
  n = 500
  nt = 200
  x0 = 0.d0
  x1 = 1.d0
  if not keyword_set(cfl) then cfl = 0.5
  gamma = 5./3.
  boundary_cond = 'periodic'

  ; total number of cell centers and cell interfaces
  nx  = n + 2 * nghost
  nxi = nx + 1

  ; cell width
  dx = (x1 - x0) / double(n)

  ; coordinates at cell interfaces
  xi = (dindgen(nxi)-nghost)*dx + x0

  ; coordinates at cell centers
  x = (dindgen(nx)-nghost)*dx + dx/2. + x0

  ; define containers to store density, momentum, energy, etc.
  rho     = dblarr(nx)
  rhou    = dblarr(nx)
  rhoetot = dblarr(nx)

  ; set the initial conditions
  rho[*] = 1. + exp(-(x-0.75)^2/0.1^2)
  index = where((x ge 0.25-0.05) and (x le 0.25+0.05))
  rho[index] = 2.
  rhou = rho ; set for constant velocity

  ; plot initial conditions
  y0 = 0.
  y1 = 3.
  window;, xpos=680
  plot, [x0, x1], [y0, y1], xstyle=1, ystyle=1, xtitle='x', ytitle='q', /nodata, charsize=1.5
  boundary, rho, rhou, rhoetot, BC=boundary_cond
  oplot, x, rho[*], psym=-6, symsize=0.7

  r = get_kbrd()
  wait, 2.

  time = 0.
  ui = dblarr(nx+1)
  for step = 0, nt do begin

     boundary, rho, rhou, rhoetot, BC=boundary_cond

     ; compute velocities at cell faces by averaging
     for i=1, nx-1 do begin
        ui[i] = 0.5 * ( rhou[i]/rho[i] + rhou[i-1]/rho[i-1] )
     endfor

     cs = 0.
     max_speed = max(abs(ui))+cs
     if max_speed gt 0. then dt = cfl*dx/max_speed else stop

     totmass = check_conservation(rho, xi, nghost)

     time = time + dt
     print, 'step=', step, ', time=', time, ', dt=', dt, ' totmass=', totmass

     advect, x, xi, rho, ui, dt, nghost, FLUXLIM=fluxlim, SCHEME=scheme
     rhou = rho ; set for constant velocity!
     boundary, rho, rhou, rhoetot, BC=boundary_cond

     erase
     plot, [x0, x1], [y0, y1], xstyle=1, ystyle=1, xtitle='x', ytitle='q', /nodata, charsize=2.0
     oplot, x, rho[*], psym=-6, symsize=0.7

     wait, 0.1

  endfor

end ; ================== advection test ===================


pro shocktube_test, FLUXlim=fluxlim

  if not keyword_set(fluxlim) then fluxlim = 'donor-cell'

  ; basic user parameters
  nghost = 2
  n = 500
  endtime = 5000.
  x0 = -50.
  x1 =  50.
  cfl = 0.8
  gamma = 5./3.
  boundary_cond = 'outflow'

  ; total number of cell centers and cell interfaces
  nx  = n + 2 * nghost
  nxi = nx + 1

  ; cell width
  dx = (x1 - x0) / double(n)

  ; coordinates at cell interfaces
  xi = (dindgen(nxi)-nghost)*dx + x0

  ; coordinates at cell centers
  x = (dindgen(nx)-nghost)*dx + dx/2. + x0

  ; define containers to store density, momentum, energy, etc.
  rho     = dblarr(nx)
  rhou    = dblarr(nx)
  rhoetot = dblarr(nx)

  ; set the initial conditions
  xmid = x0+0.5*(x1-x0)
  index = where(x le xmid) ; left state
  rho[index] = 1.d5
  rhoetot[index] = 0.5*rhou[index]^2/rho[index] + 1.0/(gamma-1.0) ; el=2.5d-5
  index = where(x gt xmid) ; right state
  rho[index] = 1.25d4
  rhoetot[index] = 0.5*rhou[index]^2/rho[index] + 0.1/(gamma-1.0) ; er=2.0d-5
  ;rhou = rho ; set for constant velocity

  ; plot initial conditions
  charsize = 1.2
  boundary, rho, rhou, rhoetot, BC=boundary_cond
  window, xsize=1400, ysize=1200;, xpos=690
  multiplot,/reset
  !p.multi=[0,1,3,0,1]
  multiplot
  y0 = 0.
  y1 = 1.2d5
  plot, [x0, x1], [y0, y1], xstyle=1, ystyle=1, ytitle='rho', /nodata, charsize=charsize
  oplot, x, rho, psym=-6, symsize=0.7
  multiplot
  y0 = 0.
  y1 = 1.1
  e = rhoetot/rho - 0.5*(rhou/rho)^2
  P = eos(rho, e, gamma)
  plot, [x0, x1], [y0, y1], xstyle=1, ystyle=1, ytitle='P', /nodata, charsize=charsize
  oplot, x, P, psym=-6, symsize=0.7
  multiplot
  y0 = -0.05
  y1 =  0.05
  u = rhou/rho
  plot, [x0, x1], [y0, y1], xstyle=1, ystyle=1, xtitle='x', ytitle='u', /nodata, charsize=charsize
  oplot, x, u, psym=-6, symsize=0.7

  e = rhoetot/rho - 0.5*(rhou/rho)^2
  P = eos(rho, e, gamma)
  cs = sqrt(gamma*P/rho)
  ;print, cs

  r = get_kbrd()
  wait, 2.

  time = 0.
  step = 0l

  while time lt endtime do begin

     ; impose CFL condition
     e = rhoetot/rho - 0.5*(rhou/rho)^2
     P = eos(rho, e, gamma)
     cs = sqrt(gamma*P/rho)
     max_speed = max(abs(rhou/rho))+max(cs)
     if max_speed gt 0. then dt = cfl*dx/max_speed else stop

     ; to get exactly to endtime in the last step
     if time+dt gt endtime then dt = endtime - time

     hydrostep, x, xi, rho, rhou, rhoetot, gamma, dt, nghost, BC=boundary_cond, FLUXLIM=fluxlim

     mass = check_conservation(rho,     xi, nghost)
     mom  = check_conservation(rhou,    xi, nghost)
     Etot = check_conservation(rhoetot, xi, nghost)

     time = time + dt
     print, 'step=', step, ', time=', time, ', dt=', dt, ' mass=', mass, ' mom=', mom, ' Etot=', Etot
     step++

     erase
     multiplot,/reset
     !p.multi=[0,1,3,0,1]
     multiplot
     y0 = 0.
     y1 = 1.2d5
     plot, [x0, x1], [y0, y1], xstyle=1, ystyle=1, ytitle='rho', /nodata, charsize=charsize
     oplot, x, rho, psym=-6, symsize=0.7
     multiplot
     y0 = 0.
     y1 = 1.1
     e = rhoetot/rho - 0.5*(rhou/rho)^2
     P = eos(rho, e, gamma)
     plot, [x0, x1], [y0, y1], xstyle=1, ystyle=1, ytitle='P', /nodata, charsize=charsize
     oplot, x, P, psym=-6, symsize=0.7
     multiplot
     y0 = -0.005
     y1 =  0.005
     u = rhou/rho
     plot, [x0, x1], [y0, y1], xstyle=1, ystyle=1, xtitle='x', ytitle='u', /nodata, charsize=charsize
     oplot, x, u, psym=-6, symsize=0.7

     wait, 0.1

  endwhile

end ; ================= shocktube test ===================


pro shockpiston_test, FLUXlim=fluxlim

  if not keyword_set(fluxlim) then fluxlim = 'donor-cell'

  ; basic user parameters
  nghost = 2
  n = 200
  endtime = 0.03
  x0 = -0.5
  x1 =  0.5
  cfl = 0.5
  gamma = 5./3.
  boundary_cond = 'outflow'

  ; total number of cell centers and cell interfaces
  nx  = n + 2 * nghost
  nxi = nx + 1

  ; cell width
  dx = (x1 - x0) / double(n)

  ; coordinates at cell interfaces
  xi = (dindgen(nxi)-nghost)*dx + x0

  ; coordinates at cell centers
  x = (dindgen(nx)-nghost)*dx + dx/2. + x0

  ; define containers to store density, momentum, energy, etc.
  rho      = dblarr(nx)
  rhou     = dblarr(nx)
  rhoetot  = dblarr(nx)
  pressure = dblarr(nx)
  velocity = dblarr(nx)

  ; set the initial conditions
  xmid = x0+0.5*(x1-x0)

  index = where(x gt xmid) ; right state
  velocity[index] = 0d0 ; gas at rest ahead of shock
  rho[index] = 1d0
  pressure[index] = 1d0/gamma ; set pressure to 1/gamma, such that cs = 1 for rho = 1

  index = where(x le xmid) ; left state
  velocity[index] = 10d0 ; piston speed; Mach 10 => results in a shock speed of about 4/3 * 10 ~ 13.33
  rho[index] = (gamma+1.0)*100d0 / ( (gamma-1.0)*100d0 + 2d0 )
  pressure[index] = ( 1.0 - gamma + 2.0*gamma*100d0 ) / (gamma+1.0)

  ; set remaining conserved variables
  rhou = rho * velocity
  rhoetot = 0.5*rhou^2/rho + pressure/(gamma-1d0) 

  ; plot initial conditions
  charsize = 1.4
  boundary, rho, rhou, rhoetot, BC=boundary_cond
  window, xsize=1400, ysize=1200;, xpos=690
  multiplot,/reset
  !p.multi=[0,1,3,0,1]
  multiplot
  y0 = 0.0
  y1 = 5.0
  plot, [x0, x1], [y0, y1], xstyle=1, ystyle=1, ytitle='rho', /nodata, charsize=charsize
  oplot, x, rho, psym=-6, symsize=0.7
  multiplot
  y0 = 0.1
  y1 = 500.0
  e = rhoetot/rho - 0.5*(rhou/rho)^2
  P = eos(rho, e, gamma)
  plot, [x0, x1], [y0, y1], xstyle=1, ystyle=1, ytitle='P', /nodata, charsize=charsize, ylog=1, ytickformat='tickformat_exponential'
  oplot, x, P, psym=-6, symsize=0.7
  multiplot
  y0 = -1.0
  y1 = 15.0
  u = rhou/rho
  plot, [x0, x1], [y0, y1], xstyle=1, ystyle=1, xtitle='x', ytitle='u', /nodata, charsize=charsize
  oplot, x, u, psym=-6, symsize=0.7

  e = rhoetot/rho - 0.5*(rhou/rho)^2
  P = eos(rho, e, gamma)
  cs = sqrt(gamma*P/rho)
  ;print, cs

  r = get_kbrd()
  wait, 2.

  time = 0.
  step = 0l

  while time lt endtime do begin

     ; impose CFL condition
     e = rhoetot/rho - 0.5*(rhou/rho)^2
     P = eos(rho, e, gamma)
     cs = sqrt(gamma*P/rho)
     max_speed = max(abs(rhou/rho))+max(cs)
     if max_speed gt 0. then dt = cfl*dx/max_speed else stop

     ; to get exactly to endtime in the last step
     if time+dt gt endtime then dt = endtime - time

     hydrostep, x, xi, rho, rhou, rhoetot, gamma, dt, nghost, BC=boundary_cond, FLUXLIM=fluxlim

     mass = check_conservation(rho,     xi, nghost)
     mom  = check_conservation(rhou,    xi, nghost)
     Etot = check_conservation(rhoetot, xi, nghost)

     time = time + dt
     print, 'step=', step, ', time=', time, ', dt=', dt, ' mass=', mass, ' mom=', mom, ' Etot=', Etot
     step++

     erase
     multiplot,/reset
     !p.multi=[0,1,3,0,1]
     multiplot
     y0 = 0.0
     y1 = 5.0
     plot, [x0, x1], [y0, y1], xstyle=1, ystyle=1, ytitle='rho', /nodata, charsize=charsize
     oplot, x, rho, psym=-6, symsize=0.7
     multiplot
     y0 = 0.1
     y1 = 500.0
     e = rhoetot/rho - 0.5*(rhou/rho)^2
     P = eos(rho, e, gamma)
     plot, [x0, x1], [y0, y1], xstyle=1, ystyle=1, ytitle='P', /nodata, charsize=charsize, ylog=1, ytickformat='tickformat_exponential'
     oplot, x, P, psym=-6, symsize=0.7
     multiplot
     y0 = -1.0
     y1 = 15.0
     u = rhou/rho
     plot, [x0, x1], [y0, y1], xstyle=1, ystyle=1, xtitle='x', ytitle='u', /nodata, charsize=charsize
     oplot, x, u, psym=-6, symsize=0.7

     wait, 0.1

  endwhile

end ; ================= shockpiston test ===================


pro sedov_test, FLUXlim=fluxlim

  if not keyword_set(fluxlim) then fluxlim = 'donor-cell'

  ; basic user parameters
  nghost = 2
  n = 500
  endtime = 1d-3
  x0 = -0.5
  x1 =  0.5
  cfl = 0.5
  gamma = 5./3.
  boundary_cond = 'outflow'

  ; total number of cell centers and cell interfaces
  nx  = n + 2 * nghost
  nxi = nx + 1

  ; cell width
  dx = (x1 - x0) / double(n)

  ; coordinates at cell interfaces
  xi = (dindgen(nxi)-nghost)*dx + x0

  ; coordinates at cell centers
  x = (dindgen(nx)-nghost)*dx + dx/2. + x0

  ; define containers to store density, momentum, energy, etc.
  rho      = dblarr(nx)
  rhou     = dblarr(nx)
  rhoetot  = dblarr(nx)
  pressure = dblarr(nx)
  velocity = dblarr(nx)

  ; set the initial conditions
  xmid = x0+0.5*(x1-x0)

  rho[*] = 1d0
  pressure[*] = 1d0/gamma ; set pressure to 1/gamma, such that cs = 1 for rho = 1
  velocity[*] = 0d0

  index = where(x ge xmid-0.05 and x le xmid+0.05)
  pressure[index] = 1d5/gamma

  ; set remaining conserved variables
  rhou = rho * velocity
  rhoetot = 0.5*rhou^2/rho + pressure/(gamma-1d0)

  ; plot initial conditions
  boundary, rho, rhou, rhoetot, BC=boundary_cond
  window, xsize=1400, ysize=1200;, xpos=690
  multiplot,/reset
  !p.multi=[0,1,3,0,1]
  multiplot
  y0 = 0.0
  y1 = 5.0
  plot, [x0, x1], [y0, y1], xstyle=1, ystyle=1, ytitle='rho', /nodata, charsize=1.4
  oplot, x, rho, psym=-6, symsize=0.7
  multiplot
  y0 = 0.1
  y1 = 2d5
  e = rhoetot/rho - 0.5*(rhou/rho)^2
  P = eos(rho, e, gamma)
  plot, [x0, x1], [y0, y1], xstyle=1, ystyle=1, ytitle='P', /nodata, charsize=1.4, ylog=1, ytickformat='tickformat_exponential'
  oplot, x, P, psym=-6, symsize=0.7
  multiplot
  y0 = -200.0
  y1 =  200.0
  u = rhou/rho
  plot, [x0, x1], [y0, y1], xstyle=1, ystyle=1, xtitle='x', ytitle='u', /nodata, charsize=1.4
  oplot, x, u, psym=-6, symsize=0.7

  r = get_kbrd()
  wait, 2.

  time = 0.
  step = 0l

  while time lt endtime do begin

     ; impose CFL condition
     e = rhoetot/rho - 0.5*(rhou/rho)^2
     P = eos(rho, e, gamma)
     cs = sqrt(gamma*P/rho)
     max_speed = max(abs(rhou/rho))+max(cs)
     if max_speed gt 0. then dt = cfl*dx/max_speed else stop

     ; to get exactly to endtime in the last step
     if time+dt gt endtime then dt = endtime - time

     hydrostep, x, xi, rho, rhou, rhoetot, gamma, dt, nghost, BC=boundary_cond, FLUXLIM=fluxlim

     mass = check_conservation(rho,     xi, nghost)
     mom  = check_conservation(rhou,    xi, nghost)
     Etot = check_conservation(rhoetot, xi, nghost)

     time = time + dt
     print, 'step=', step, ', time=', time, ', dt=', dt, ' mass=', mass, ' mom=', mom, ' Etot=', Etot
     step++

     erase
     multiplot,/reset
     !p.multi=[0,1,3,0,1]
     multiplot
     y0 = 0.0
     y1 = 5.0
     plot, [x0, x1], [y0, y1], xstyle=1, ystyle=1, ytitle='rho', /nodata, charsize=1.4
     oplot, x, rho, psym=-6, symsize=0.7
     multiplot
     y0 = 0.1
     y1 = 2d5
     e = rhoetot/rho - 0.5*(rhou/rho)^2
     P = eos(rho, e, gamma)
     plot, [x0, x1], [y0, y1], xstyle=1, ystyle=1, ytitle='P', /nodata, charsize=1.4, ylog=1, ytickformat='tickformat_exponential'
     oplot, x, P, psym=-6, symsize=0.7
     multiplot
     y0 = -200.0
     y1 =  200.0
     u = rhou/rho
     plot, [x0, x1], [y0, y1], xstyle=1, ystyle=1, xtitle='x', ytitle='u', /nodata, charsize=1.4
     oplot, x, u, psym=-6, symsize=0.7

     wait, 0.1

  endwhile

end ; ================= sedov test ===================


pro hydro_test, FLUXLIM=fluxlim

  if not keyword_set(fluxlim) then fluxlim = 'vanleer'

  ; basic user parameters
  nghost = 2
  n = 100
  nt = 1000
  x0 = 0.d0
  x1 = 1.d0
  cfl = 0.5
  gamma = 7./5.
  boundary_cond = 'periodic'

  ; total number of cell centers and cell interfaces
  nx  = n + 2 * nghost
  nxi = nx + 1

  ; cell width
  dx = (x1 - x0) / double(n)

  ; coordinates at cell interfaces
  xi = (dindgen(nxi)-nghost)*dx + x0

  ; coordinates at cell centers
  x = (dindgen(nx)-nghost)*dx + dx/2. + x0

  ; define containers to store density, momentum, energy, etc.
  rho     = dblarr(nx)
  rhou    = dblarr(nx)
  rhoetot = dblarr(nx)

  ; set the initial conditions
  rho[*] = 1. + exp(-(x-0.75)^2/0.1^2)
  index = where((x ge 0.25-0.05) and (x le 0.25+0.05))
  rho[index] = 2.
  rhou = 0.01*rho ; set for constant velocity
  rhoetot = 0.5*rhou^2/rho + 1.*rho
 
  ; plot initial conditions
  y0 = 0.
  y1 = 3.
  window, xsize=2000, ysize=1000;, xpos=690
  ;window;, xpos=680
  plot, [x0, x1], [y0, y1], xstyle=1, ystyle=1, xtitle='x', ytitle='q', /nodata, charsize=2.0
  boundary, rho, rhou, rhoetot, BC=boundary_cond
  oplot, x, rho[*], psym=-6, symsize=0.7

  r = get_kbrd()
  wait, 2.

  time = 0.

  for step = 0, nt do begin

     ; impose CFL condition
     e = rhoetot/rho - 0.5*(rhou/rho)^2
     P = eos(rho, e, gamma)
     cs = sqrt(gamma*P/rho)
     max_speed = max(abs(rhou/rho))+max(cs)
     if max_speed gt 0. then dt = cfl*dx/max_speed else stop

     hydrostep, x, xi, rho, rhou, rhoetot, gamma, dt, nghost, BC=boundary_cond, FLUXLIM=fluxlim

     mass = check_conservation(rho,     xi, nghost)
     mom  = check_conservation(rhou,    xi, nghost)
     Etot = check_conservation(rhoetot, xi, nghost)

     time = time + dt
     print, 'step=', step, ', time=', time, ', dt=', dt, ' mass=', mass, ' mom=', mom, ' Etot=', Etot

     erase
     plot, [x0, x1], [y0, y1], xstyle=1, ystyle=1, xtitle='x', ytitle='q', /nodata, charsize=2.0
     oplot, x, rho[*], psym=-6, symsize=0.7

     wait, 0.1

  endfor

end ; ================= hydro test ===================


; === START subroutine advect ===
;
pro advect, x, xi, q, ui, dt, nghost, FLUXLIM=fluxlim, SCHEME=scheme

  nx = n_elements(xi)-1
  if n_elements(x) ne nx then stop
  if n_elements(xi) ne nx+1 then stop
  if n_elements(q) ne nx then stop
  if n_elements(ui) ne nx+1 then stop
  if nghost lt 1 then stop
  if not keyword_set(fluxlim) then fluxlim = 'mc'
  if not keyword_set(scheme) then scheme = 'upwind'

  ; determine the r_{i-1/2} for the flux limiter
  r = dblarr(nx+1)
  for i=2, nx-2 do begin
     dq = (q[i]-q[i-1])
     if abs(dq) gt 0.d0 then begin
        if(ui[i] ge 0.d0) then begin
           r[i] = (q[i-1]-q[i-2])/dq
        endif else begin
           r[i] = (q[i+1]-q[i])/dq
           endelse
     endif
  endfor

  ; Determine the flux limiter
  ; (many other flux limiters can be implemented here!)
  case fluxlim of
     'donor-cell': begin
        phi = dblarr(nx+1)
     end
     'lax-wendroff': begin
        phi = dblarr(nx+1)
        phi[*] = 1.
     end
     'beam-warming': begin
        phi = r
     end
     'fromm': begin
        phi = 0.5*(r+1.)
     end
     'minmod': begin
        phi = dblarr(nx+1)
        for i=1, nx-1 do begin
           a = 1.
           b = r[i]
           if (abs(a) lt abs(b)) and (a*b gt 0.) then phi[i] = a
           if (abs(a) gt abs(b)) and (a*b gt 0.) then phi[i] = b
           if a*b le 0. then phi[i] = 0.
        endfor
     end
     'superbee': begin
        phi = dblarr(nx+1)
        for i=1, nx-1 do begin
           a = min([1.d0,2.d0*r[i]])
           b = min([2.d0,r[i]])
           phi[i] = max([0.d0,a,b])
        endfor
     end
     'mc': begin
        phi = dblarr(nx+1)
        for i=1, nx-1 do begin
           a = min( [ (1.+r[i])/2., 2., 2.*r[i] ] )
           phi[i] = max([0.,a])
        endfor
     end
     'vanleer': begin
        phi = dblarr(nx+1)
        for i=1, nx-1 do begin
           phi[i] = ( r[i] + abs(r[i]) ) / ( 1. + abs(r[i]) )
        endfor
     end
     else: stop
  endcase

  ; now construct the flux
  flux = dblarr(nx+1)
  for i=1, nx-1 do begin
     case scheme of
      'upwind': begin
       if ui[i] ge 0. then begin
          flux[i] = ui[i] * q[i-1]
       endif else begin
          flux[i] = ui[i] * q[i]
       endelse
      end
      'centred': begin
           flux[i] = ui[i] * (q[i]+q[i-1])/2.
      end
      'downwind': begin
       if ui[i] ge 0. then begin
          flux[i] = ui[i] * q[i]
       endif else begin
          flux[i] = ui[i] * q[i-1]
       endelse
      end
     endcase
     fluxlimiter_correction = 0.5 * abs(ui[i]) * $
                              (1-abs(ui[i]*dt/(x[i]-x[i-1]))) * $
                              phi[i] * (q[i]-q[i-1])
     flux[i] = flux[i] + fluxlimiter_correction
  endfor

  ; update the cells, except the ghost cells
  for i=nghost, nx-1-nghost do begin
     q[i] = q[i] - dt * ( flux[i+1] - flux[i] ) / ( xi[i+1] - xi[i] )
  endfor

end ; === END subroutine advect ===

;
; === START soubroutine hydrostep ===
;
pro hydrostep, x, xi, rho, rhou, rhoetot, gamma, dt, nghost, BC=bc, FLUXLIM=fluxlim

  if not keyword_set(fluxlim) then fluxlim = 'minmod'
  nx = n_elements(x)

  ; impose boundary conditions
  boundary, rho, rhou, rhoetot, BC=bc

  ; compute the velocity at the cell interfaces (by averaging)
  ui = dblarr(nx+1)
  for i=1, nx-1 do begin
     ui[i] = 0.5 * ( rhou[i]/rho[i] + rhou[i-1]/rho[i-1] )
  endfor

  ; advect rho
  advect, x, xi, rho, ui, dt, nghost, FLUXLIM=fluxlim

  ; advect rho*u
  advect, x, xi, rhou, ui, dt, nghost, FLUXLIM=fluxlim

  ; advect rho*etot
  advect, x, xi, rhoetot, ui, dt, nghost, FLUXLIM=fluxlim

  ; re-impose BCs
  boundary, rho, rhou, rhoetot, BC=bc

  ; compute pressure (call equation of state)
  e = rhoetot/rho - 0.5*(rhou/rho)^2
  P = eos(rho, e, gamma)

  ; Now add the pressure force and energy, for all cells except the ghost cells
  ; operator-splitting
  for i=nghost, nx-1-nghost do begin
     rhou   [i] = rhou   [i] - dt * (P[i+1]-P[i-1])                 / (x[i+1]-x[i-1])
     rhoetot[i] = rhoetot[i] - dt * (P[i+1]*ui[i+1]-P[i-1]*ui[i-1]) / (x[i+1]-x[i-1])
  endfor

  ; re-impose BCs (not strictly necessary)
  boundary, rho, rhou, rhoetot, BC=bc

end ; === END soubroutine hydrostep ===

;
; === START soubroutine boundary (FOR 2 ghost cells) ===
;
pro boundary, rho, rhou, rhoetot, BC=bc
;
; Get the number of grid points including the ghost cells
;
  nx = n_elements(rho)
  case bc of
     'periodic': begin
        rho    [0]    = rho    [nx-4]
        rho    [1]    = rho    [nx-3]
        rho    [nx-2] = rho    [2]
        rho    [nx-1] = rho    [3]
        rhou   [0]    = rhou   [nx-4]
        rhou   [1]    = rhou   [nx-3]
        rhou   [nx-2] = rhou   [2]
        rhou   [nx-1] = rhou   [3]
        rhoetot[0]    = rhoetot[nx-4]
        rhoetot[1]    = rhoetot[nx-3]
        rhoetot[nx-2] = rhoetot[2]
        rhoetot[nx-1] = rhoetot[3]
     end
     'reflect': begin
        rho    [0]    =  rho    [3]
        rho    [1]    =  rho    [2]
        rho    [nx-2] =  rho    [nx-3]
        rho    [nx-1] =  rho    [nx-4]
        rhou   [0]    = -rhou   [3]
        rhou   [1]    = -rhou   [2]
        rhou   [nx-2] = -rhou   [nx-3]
        rhou   [nx-1] = -rhou   [nx-4]
        rhoetot[0]    =  rhoetot[3]
        rhoetot[1]    =  rhoetot[2]
        rhoetot[nx-2] =  rhoetot[nx-3]
        rhoetot[nx-1] =  rhoetot[nx-4]
     end
     'outflow': begin
        rho    [0]    =  rho    [2]
        rho    [1]    =  rho    [2]
        rho    [nx-2] =  rho    [nx-3]
        rho    [nx-1] =  rho    [nx-3]
        rhou   [0]    =  rhou   [2]
        rhou   [1]    =  rhou   [2]
        rhou   [nx-2] =  rhou   [nx-3]
        rhou   [nx-1] =  rhou   [nx-3]
        rhoetot[0]    =  rhoetot[2]
        rhoetot[1]    =  rhoetot[2]
        rhoetot[nx-2] =  rhoetot[nx-3]
        rhoetot[nx-1] =  rhoetot[nx-3]
     end
     else: stop
  endcase

end ; === END soubroutine boundary ===

;
; === START function eos ===
;
function eos, rho, e, gamma

  P = (gamma-1.)*rho*e

  return, P

end ; === END function eos ===

;
; === START function check_conservation ===
;
function check_conservation, q, xi, nghost

  nx = n_elements(xi)-1
  dx = xi[1:nx]-xi[0:nx-1]

  ; compute integral over conserved q (excluding ghost zones)
  istart = nghost
  iend   = nx-1-nghost
  integral = total(q[istart:iend]*dx[istart:iend],/double)

  return, integral

end ; === END function check_conservation ===

;+
; NAME:
;	MULTIPLOT
; PURPOSE:
;	Create multiple plots with shared axes.
; EXPLANATION:
;	This procedure makes a matrix of plots with *SHARED AXES*, either using
;	parameters passed to multiplot or !p.multi in a non-standard way.
;	It is good for data with one or two shared axes and retains all the
;	versatility of the plot commands (e.g. all keywords and log scaling).
;	The plots are connected with the shared axes, which saves space by
;	omitting redundant ticklabels and titles.  Multiplot does this by
;	setting !p.position, !x.tickname and !y.tickname automatically.
;	A call (multiplot,/reset) restores original values.
;
;	Note: This method may be superseded by future improvements in !p.multi
;	by RSI.  For now, it's a good way to gang plots together.
; CALLING SEQUENCE:
;	multiplot[pmulti][,/help][,/initialize][,/reset][,/rowmajor]
; EXAMPLES:
;	multiplot,/help			; print this header.
;	; Then copy & paste, from your xterm, the following lines to test:
;
;	x = findgen(100)		;	   MULTIPLOT
;	t=exp(-(x-50)^2/300)		;	 -------------------------
;	erase				;	 |           |           |
;	u=exp(-x/30)			;	 |           |           |
;	y = sin(x)			;	 |  UL plot  |  UR plot  |
;	r = reverse(y*u)		;	 |           |           |
;	!p.multi=[0,2,2,0,0]		;	 |           |           |
;	multiplot 	 		;	y-------------------------
;	plot,x,y*u,title='MULTIPLOT'	;	l|           |           |
;	multiplot & plot,x,r 		;	a|           |           |
;	multiplot 			;	b|  LL plot  |  LR plot  |
;	plot,x,y*t,ytit='ylabels'	;	e|           |           |
;	multiplot 			;	l|           |           |
;	plot,x,y*t,xtit='xlabels'	;	s-------------------------
;	multiplot,/reset		;		        xlabels
;					 
;	wait,2 & erase			;		 TEST
;	multiplot,[1,3]			;	H------------------------
;	plot,x,y*u,title='TEST'		;	E|	plot #1		|
;	multiplot			;	I------------------------
;	plot,x,y*t,ytit='HEIGHT'	;	G|	plot #2 	|
;	multiplot			;	H------------------------
;	plot,x,r,xtit='PHASE'		;	T|	plot #3		|
;	multiplot,/reset		;	 ------------------------
;					;		 PHASE
;
;	multiplot,[1,1],/init,/verbose	; one way to return to single plot
;	% MULTIPLOT: Initialized for 1x1, plotted across then down (column major).
; OPTIONAL INPUTS:
;	pmulti = 2-element or 5-element vector giving number of plots, e.g.,
;	  multiplot,[1,6]		; 6 plots vertically
;	  multiplot,[0,4,2,0,0]		; 4 plots along x and 2 along y
;	  multiplot,[0,4,2,0,1]		; ditto, except rowmajor (down 1st)
;	  multiplot,[4,2],/rowmajor 	; identical to previous line
; OPTIONAL KEYWORDS:
;	help = flag to print header
;	initialize = flag to begin only---no plotting, just setup,
;	  e.g., multiplot,[4,2],/init,/verbose & multiplot & plot,x,y
;	reset = flag to reset system variables to values prior to /init
;	default = flag to restore IDL's default value for system variables
;	rowmajor = flag to number plots down column first (D=columnmajor)
;	verbose = flag to output informational messages
; Outputs:
;	!p.position = 4-element vector to place a plot
;	!x.tickname = either '' or else 30 ' ' to suppress ticknames
;	!y.tickname = either '' or else 30 ' ' to suppress ticknames
;	!p.noerase = 1
; Common blocks:
;	multiplot---to hold saved variables and plot counter.  See code.
; Side Effects:
;	Multiplot sets a number of system variables: !p.position, !p.multi,
;	!x.tickname, !y.tickname, !P.noerase---but all can be reset with
;	the call: multiplot,/reset
; RESTRICTIONS:
;	1. If you use !p.multi as the method of telling how many plots
;	are present, you have to set !p.multi at the beginning each time you
;	use multiplot or call multiplot with the /reset keyword.
;	2. There's no way to make an xtitle or ytitle span more than one plot,
;	except by adding spaces to shift it or to add it manually with xyouts.
;	3. There is no way to make plots of different sizes; each plot
;	covers the same area on the screen or paper.
; PROCEDURE:
;	This routine makes a matrix of plots with common axes, as opposed to
;	the method of !p.multi where axes are separated to allow labels.
;	Here the plots are joined and labels are suppressed, except at the
;	left edge and the bottom.  You tell multiplot how many plots to make
;	using either !p.multi (which is then reset) or the parameter pmulti.
;	However, multiplot keeps track of the position by itself because
;	!p.multi interacts poorly with !p.position.
; MODIFICATION HISTORY:
;	write, 21-23 Mar 94, Fred Knight (knight@ll.mit.edu)
;	alter plot command that sets !x.window, etc. per suggestion of
;	  Mark Hadfield (hadfield@storm.greta.cri.nz), 7 Apr 94, FKK
;	add a /default keyword restore IDL's default values of system vars,
;	  7 Apr 94, FKK
;	modify two more sys vars !x(y).tickformat to suppress user-formatted
;	  ticknames, per suggestion of Mark Hadfield (qv), 8 Apr 94, FKK
;	Converted to IDL V5.0   W. Landsman   September 1997
;-
pro multiplot,help=help,pmulti $
  ,initialize=initialize,reset=reset,default=default $
  ,rowmajor=rowmajor,verbose=verbose
;
;	=====>> COMMON
;
common multiplot $
  ,nplots $	; [# of plots along x, # of plots along y]
  ,nleft $	; # of plots remaining---like the first element of !p.multi
  ,pdotmulti $	; saved value of !p.multi
  ,margins $	; calculated margins based on !p.multi or pmulti
  ,pposition $	; saved value of !p.position
  ,colmajor $	; flag for column major order
  ,noerase $	; saved value of !p.noerase
  ,xtickname $	; Original value
  ,ytickname $	; Original value
  ,xtickformat $; Original value
  ,ytickformat 	; Original value
;
;	=====>> HELP
;
;on_error,2
if keyword_set(help) then begin & doc_library,'multiplot' & return & endif
;
;	=====>> RESTORE IDL's DEFAULT VALUES (kill multiplot's influence)
;
if keyword_set(default) then begin
  !p.position = 0
  !x.tickname = ''
  !y.tickname = ''
  !x.tickformat = ''
  !y.tickformat = ''
  !p.multi = 0
  !p.noerase = 0
  nleft = 0
  nplots = [1,1]
  pdotmulti = !p.multi
  margins = 0
  pposition = !p.position
  noerase = !p.noerase
  xtickname = !x.tickname
  ytickname = !y.tickname
  xtickformat = !x.tickformat
  ytickformat = !y.tickformat
  if keyword_set(verbose) then begin
    message,/inform,'Restore IDL''s defaults for affected system variables.'
    message,/inform,'Reset multiplot''s common to IDL''s defaults.'
    endif
  return
  endif
;
;	=====>> RESTORE SAVED SYSTEM VARIABLES
;
if keyword_set(reset) then begin
  if n_elements(pposition) gt 0 then begin
    !p.position = pposition
    !x.tickname = xtickname
    !y.tickname = ytickname
    !x.tickformat = xtickformat
    !y.tickformat = ytickformat
    !p.multi = pdotmulti
    !p.noerase = noerase
    endif
  nleft = 0
  if keyword_set(verbose) then begin
    coords = '['+string(!p.position,form='(3(f4.2,","),f4.2)')+']'
    multi = '['+string(!p.multi,form='(4(i2,","),i2)')+']'
    message,/inform,'Reset.  !p.position='+coords+', !p.multi='+multi
    endif
  return
  endif
;
;	=====>> SETUP: nplots, MARGINS, & SAVED SYSTEM VARIABLES
;
if n_elements(nleft) eq 1 then init = (nleft eq 0) else init = 1
if (n_elements(pmulti) eq 2) or (n_elements(pmulti) eq 5) then init = 1
if (n_elements(!p.multi) eq 5) then begin
  if (!p.multi[1] gt 0) and (!p.multi[2] gt 0) then init = (!p.multi[0] eq 0) 
  endif
if init or keyword_set(initialize) then begin
  case n_elements(pmulti) of
  0:begin
    if n_elements(!p.multi) eq 1 then return	; NOTHING TO SET
    if n_elements(!p.multi) ne 5 then message,'Bogus !p.multi; aborting.'
    nplots = !p.multi[1:2] > 1
    if keyword_set(rowmajor) then colmajor = 0 else colmajor = !p.multi[4] eq 0
    end
  2:begin
    nplots = pmulti
    colmajor = not keyword_set(rowmajor)	; D=colmajor: left to rt 1st
    end
  5:begin
    nplots = pmulti[1:2]
    if keyword_set(rowmajor) then colmajor = 0 else colmajor = pmulti[4] eq 0
    end
  else: message,'pmulti can only have 0, 2, or 5 elements.'
  endcase
  pposition = !p.position			; save sysvar to be altered
  xtickname = !x.tickname
  ytickname = !y.tickname
  xtickformat = !x.tickformat
  ytickformat = !y.tickformat
  pdotmulti = !p.multi
  nleft = nplots[0]*nplots[1]			; total # of plots
  !p.position = 0				; reset
  !p.multi = 0
  plot,/nodata,xstyle=4,ystyle=4,!x.range,!y.range,/noerase	; set window & region
  margins = [min(!x.window)-min(!x.region) $	; in normlized coordinates
    ,min(!y.window)-min(!y.region) $
    ,max(!x.region)-max(!x.window) $
    ,max(!y.region)-max(!y.window)]
  noerase = !p.noerase
  !p.noerase = 1				; !p.multi does the same
  if keyword_set(verbose) then begin
    major = ['across then down (column major).','down then across (row major).']
    if colmajor then index = 0 else index = 1
    message,/inform,'Initialized for '+strtrim(nplots[0],2) $
    +'x'+strtrim(nplots[1],2)+', plotted '+major[index]
    endif
;  print,margins,'=margins'
  if keyword_set(initialize) then return
endif
;
;	=====>> Define the plot region without using !p.multi.
;
cols = nplots[0]			; for convenience
rows = nplots[1]
nleft = nleft - 1			; decrement plots remaining
cur = cols*rows - nleft			; current plot #: 1 to cols*rows
idx = [(1.-margins[0]-margins[2])/cols $
  ,(1.-margins[1]-margins[3])/rows]	; normalized coords per plot
if colmajor then begin			; location in matrix of plots
  col = cur mod cols
  if col eq 0 then col = cols
  row = (cur-1)/cols + 1
endif else begin			; here (1,2) is 1st col, 2nd row
  row = cur mod rows
  if row eq 0 then row = rows
  col = (cur-1)/rows + 1
endelse
pos = [(col-1)*idx[0],(rows-row)*idx[1],col*idx[0],(rows-row+1)*idx[1]] $
  + [margins[0],margins[1],margins[0],margins[1]]
;print,row,col,rows,cols,pos
;
;	=====>> Finally set the system variables; user shouldn't change them.
;
!p.position = pos
onbottom = (row eq rows) or (rows eq 1)
onleft = (col eq 1) or (cols eq 1)
if onbottom then !x.tickname = xtickname else !x.tickname = replicate(' ',30)
if onleft then !y.tickname = ytickname else !y.tickname = replicate(' ',30)
if onbottom then !x.tickformat = xtickformat else !x.tickformat = ''
if onleft then !y.tickformat = ytickformat else !y.tickformat = ''
if keyword_set(verbose) then begin
  coords = '['+string(pos,form='(3(f4.2,","),f4.2)')+']'
  plotno = 'Setup for plot ['+strtrim(col,2)+','+strtrim(row,2)+'] of ' $
    +strtrim(cols,2)+'x'+strtrim(rows,2)
  message,/inform,plotno+' at '+coords
  endif
;stop
return
end

;
; format string to produce exponential display
; in tickmarks
;
; assume tickmark on every decade
;
function tickformat_exponential, axis, index, value 

  ex    = long(alog10(value))
  if (ex LT 0) then begin
      if (ex GT -10) then begin
          case ex of
              -1: fstr = string(value,FORMAT='(f3.1)')
              -2: fstr = string(value,FORMAT='(f4.2)')
              -3: fstr = string(value,FORMAT='(f5.3)')
              else: fstr = string(ex,FORMAT='("10!U",i2.1,"!N")')
;               fstr = string(ex,FORMAT='("10!U",i2.1,"!N")')
          endcase
      endif else begin
          if (ex GT -100) then begin
              fstr = string(ex,FORMAT='("10!U",i3.2,"!N")')
          endif else fstr = string(ex,FORMAT='("10!U",i4.3,"!N")')
      endelse
  endif else begin
      if (ex LT 10) then begin
          case ex of
              0: fstr = '1'
              1: fstr = '10'
              2: fstr = '100'
              else: fstr = string(ex,FORMAT='("10!U",i1.1,"!N")')
;               fstr = string(ex,FORMAT='("10!U",i1.1,"!N")')
          endcase
      endif else begin
          if (ex LT 100) then begin
              fstr = string(ex,FORMAT='("10!U",i2.2,"!N")')
          endif else fstr = string(ex,FORMAT='("10!U",i3.3,"!N")')
      endelse
  endelse

  return, fstr
END