diff --git a/ext/UnderwaterAcousticsPlots.jl b/ext/UnderwaterAcousticsPlots.jl
index 20c0e24..1612185 100644
--- a/ext/UnderwaterAcousticsPlots.jl
+++ b/ext/UnderwaterAcousticsPlots.jl
@@ -4,6 +4,7 @@ using RecipesBase
 using Colors
 using UnderwaterAcoustics
 import UnderwaterAcoustics: AbstractAcousticSource, AbstractAcousticReceiver, RayArrival, value
+import UnderwaterAcoustics: SampledFieldZ, SampledFieldX, SampledFieldXZ, SampledFieldXY
 import DSP: amp2db
 
 @recipe function plot(env::UnderwaterEnvironment)
@@ -16,12 +17,12 @@ import DSP: amp2db
   xguide --> "x (m)"
   yguide --> "z (m)"
   @series begin
-    seriestype := :line
+    seriestype := :path
     linecolor := :royalblue
     xrange, [value(z1, (x, 0.0)) for x ∈ xrange]
   end
   @series begin
-    seriestype := :line
+    seriestype := :path
     linecolor := :brown
     xrange, [-value(z2, (x, 0.0)) for x ∈ xrange]
   end
@@ -107,11 +108,61 @@ end
   for (j, r) ∈ enumerate(rays)
     r = r.path
     @series begin
-      seriestype := :line
+      seriestype := :path
       linecolor := colors[cndx[j]]
       [r[i][1] for i ∈ 1:length(r)], [r[i][3] for i ∈ 1:length(r)]
     end
   end
 end
 
+@recipe function plot(fld::SampledFieldZ; npts=1000)
+  zr = range(minimum(fld.zrange), maximum(fld.zrange); length=npts)
+  ticks --> :native
+  legend --> false
+  yguide --> "z (m)"
+  @series begin
+    seriestype := :path
+    [fld(z) for z ∈ zr], zr
+  end
+end
+
+@recipe function plot(fld::SampledFieldX; npts=1000)
+  xr = range(minimum(fld.xrange), maximum(fld.xrange); length=npts)
+  ticks --> :native
+  legend --> false
+  xguide --> "x (m)"
+  @series begin
+    seriestype := :path
+    xr, [fld(x) for x ∈ xr]
+  end
+end
+
+@recipe function plot(fld::SampledFieldXY; npts=1000)
+  xr = range(minimum(fld.xrange), maximum(fld.xrange); length=npts)
+  yr = range(minimum(fld.yrange), maximum(fld.yrange); length=npts)
+  ticks --> :native
+  legend --> false
+  xguide --> "x (m)"
+  yguide --> "y (m)"
+  colorbar --> true
+  @series begin
+    seriestype := :heatmap
+    xr, yr, [fld(x, y) for x ∈ xr, y ∈ yr]
+  end
+end
+
+@recipe function plot(fld::SampledFieldXZ; npts=1000)
+  xr = range(minimum(fld.xrange), maximum(fld.xrange); length=npts)
+  zr = range(minimum(fld.zrange), maximum(fld.zrange); length=npts)
+  ticks --> :native
+  legend --> false
+  xguide --> "x (m)"
+  yguide --> "z (m)"
+  colorbar --> true
+  @series begin
+    seriestype := :heatmap
+    xr, zr, [fld(x, z) for x ∈ xr, z ∈ zr]
+  end
+end
+
 end # module
diff --git a/src/pm_stdlib.jl b/src/pm_stdlib.jl
index 5ce90d5..2070396 100644
--- a/src/pm_stdlib.jl
+++ b/src/pm_stdlib.jl
@@ -282,10 +282,12 @@ sampled on a regular grid.
 function SampledField(v; x=nothing, y=nothing, z=nothing, interp=:linear)
   interp === :linear || error("Only linear interpolation supported")
   if x === nothing && y === nothing && z !== nothing
-    f = linear_interpolation(z, v; extrapolation_bc=Flat())
+    ndx = sortperm(z)
+    f = linear_interpolation(z[ndx], v[ndx]; extrapolation_bc=Flat())
     SampledFieldZ(f, z)
   elseif x !== nothing && y === nothing && z === nothing
-    f = linear_interpolation(x, v; extrapolation_bc=Flat())
+    ndx = sortperm(x)
+    f = linear_interpolation(x[ndx], v[ndx]; extrapolation_bc=Flat())
     SampledFieldX(f, x)
   elseif x !== nothing && y === nothing && z !== nothing
     f = extrapolate(interpolate((x, z), v, Gridded(Linear())), Flat())