diff --git a/examples/euler_2d/euler_2d.py b/examples/euler_2d/euler_2d.py
index 6e897b394..c31a3c2e0 100755
--- a/examples/euler_2d/euler_2d.py
+++ b/examples/euler_2d/euler_2d.py
@@ -21,6 +21,7 @@
 """
 from clawpack import pyclaw
 from clawpack import riemann
+from clawpack.riemann.euler_4wave_2D_constants import *
 
 def load_frame(frame_number):
     from clawpack.pyclaw import Solution
@@ -31,7 +32,7 @@ def plot_frame(frame):
     import matplotlib.pyplot as plt
     q = frame.q
     x, y = frame.state.grid.c_centers
-    plt.pcolormesh(x, y, q[0,...])
+    plt.pcolormesh(x, y, q[density,...])
 
 def plot_results():
     from clawpack.visclaw import iplot
@@ -42,17 +43,17 @@ def plot_results():
 solver.all_bcs = pyclaw.BC.extrap
 
 domain = pyclaw.Domain([0.,0.],[1.,1.],[200,200])
-solution = pyclaw.Solution(solver.num_eqn,domain)
+solution = pyclaw.Solution(num_eqn,domain)
 gamma = 1.4
 solution.problem_data['gamma']  = gamma
 
 # Set initial data
 xx,yy = domain.grid.p_centers
 l = xx<0.5; r = xx>=0.5; b = yy<0.5; t = yy>=0.5
-solution.q[0,...] = 2.*l*t + 1.*l*b + 1.*r*t + 3.*r*b
-solution.q[1,...] = 0.75*t - 0.75*b
-solution.q[2,...] = 0.5*l  - 0.5*r
-solution.q[3,...] = 0.5*solution.q[0,...]*(solution.q[1,...]**2+solution.q[2,...]**2) + 1./(gamma-1.)
+solution.q[density,...] = 2.*l*t + 1.*l*b + 1.*r*t + 3.*r*b
+solution.q[x_momentum,...] = 0.75*t - 0.75*b
+solution.q[y_momentum,...] = 0.5*l  - 0.5*r
+solution.q[energy,...] = 0.5*solution.q[density,...]*(solution.q[x_momentum,...]**2+solution.q[y_momentum,...]**2) + 1./(gamma-1.)
 
 claw = pyclaw.Controller()
 claw.tfinal = 0.3