From 007058ef0a4ccf5b04049c41ba252df4c2562f07 Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Sun, 19 Nov 2023 16:55:48 -0500
Subject: [PATCH 1/6] Add some @view

---
 src/plot/utility.jl |  3 +--
 src/select.jl       | 14 +++++++-------
 2 files changed, 8 insertions(+), 9 deletions(-)

diff --git a/src/plot/utility.jl b/src/plot/utility.jl
index bf6dba5..37d03f6 100644
--- a/src/plot/utility.jl
+++ b/src/plot/utility.jl
@@ -3,8 +3,7 @@
 "Prepare 2D data arrays for passing to plotting functions."
 function getdata(data::BATLData, var::AbstractString, plotrange, plotinterval; griddim=1,
    innermask::Bool=false)
-   x, w = data.x, data.w
-   ndim = data.head.ndim
+   x, w, ndim = data.x, data.w, data.head.ndim
    @assert ndim == 2 "data must be in 2D!"
 
    varIndex_ = findindex(data, var)
diff --git a/src/select.jl b/src/select.jl
index d1b17ab..177d4d8 100644
--- a/src/select.jl
+++ b/src/select.jl
@@ -98,7 +98,7 @@ end
 	subsurface(x, y, data, limits)
 	subsurface(x, y, u, v, limits)
 
-Extract subset of 2D surface dataset. See also: [`subvolume`](@ref).
+Extract subset of 2D surface dataset in ndgrid format. See also: [`subvolume`](@ref).
 """
 function subsurface(x, y, data, limits)
    checkvalidlimits(limits)
@@ -188,24 +188,24 @@ function subdata(data, xind, yind, zind, sz)
    newsz = size(newdata)
 
    if length(sz) > 3
-      newdata = reshape(newdata, (newsz[1:3]..., sz[4:end]))
+      newdata = @views reshape(newdata, (newsz[1:3]..., sz[4:end]))
    end
 
    newdata
 end
 
 "Return variable data from string `var`."
-function getvar(data::BATLData, var)
+function getvar(data::BATLData, var::AbstractString)
    VarIndex_ = findfirst(x->x==lowercase(var), lowercase.(data.head.wnames))
-   isnothing(VarIndex_) && error("$(var) not found in file header variables!")
+   isnothing(VarIndex_) && error("$var not found in file header variables!")
 
    ndim = data.head.ndim
    if ndim == 1
-      w = data.w[:,VarIndex_]
+      w = @view data.w[:,VarIndex_]
    elseif ndim == 2
-      w = data.w[:,:,VarIndex_]
+      w = @view data.w[:,:,VarIndex_]
    elseif ndim == 3
-      w = data.w[:,:,:,VarIndex_]
+      w = @view data.w[:,:,:,VarIndex_]
    end
    w
 end

From e79ae73e8d54dd31ca75f9c7c8bfaa19576656ea Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Sun, 19 Nov 2023 17:56:23 -0500
Subject: [PATCH 2/6] Use capital cases

---
 README.md | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/README.md b/README.md
index 4624c6c..4343743 100644
--- a/README.md
+++ b/README.md
@@ -10,11 +10,11 @@ Fast [BATSRUS](https://github.com/MSTEM-QUDA/BATSRUS)/[SWMF](https://github.com/
 
 This package provides the following functionalities:
 
-* simulation data reader
+* Simulation data reader
 * 2D/3D region cut from the whole domain
-* interpolation from unstructured to structured data
-* data format conversion to VTK
-* simulation data visualization through multiple plotting libraries
+* Interpolation from unstructured to structured data
+* Data format conversion to VTK
+* Data visualization through multiple plotting libraries
 
 For more details, please check the [document][Batsrus-doc].
 

From b3274e18067b9a194446f6c29cb9972c09011cd4 Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Sun, 19 Nov 2023 18:46:05 -0500
Subject: [PATCH 3/6] Add Makie support as an extension

---
 Project.toml                           |  8 +++++-
 ext/BatsrusMakieExt/BatsrusMakieExt.jl | 10 +++++++
 ext/BatsrusMakieExt/typerecipe.jl      | 35 ++++++++++++++++++++++
 src/Batsrus.jl                         |  2 +-
 src/io.jl                              | 10 +++----
 src/plot/makie.jl                      | 40 --------------------------
 src/plot/plots.jl                      |  9 ++----
 src/plot/pyplot.jl                     | 28 +++++++++---------
 src/plot/utility.jl                    |  8 ++++++
 9 files changed, 84 insertions(+), 66 deletions(-)
 create mode 100644 ext/BatsrusMakieExt/BatsrusMakieExt.jl
 create mode 100644 ext/BatsrusMakieExt/typerecipe.jl
 delete mode 100644 src/plot/makie.jl

diff --git a/Project.toml b/Project.toml
index 70c8bec..12dcfd3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Batsrus"
 uuid = "e74ebddf-6ac1-4047-a0e5-c32c99e57753"
 authors = ["Hongyang Zhou <hyzhou@umich.edu>"]
-version = "0.3.11"
+version = "0.4.0"
 
 [deps]
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
@@ -15,6 +15,12 @@ Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
+[weakdeps]
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+
+[extensions]
+BatsrusMakieExt = "Makie"
+
 [compat]
 Dierckx = "0.5"
 FortranFiles = "0.6"
diff --git a/ext/BatsrusMakieExt/BatsrusMakieExt.jl b/ext/BatsrusMakieExt/BatsrusMakieExt.jl
new file mode 100644
index 0000000..17099f4
--- /dev/null
+++ b/ext/BatsrusMakieExt/BatsrusMakieExt.jl
@@ -0,0 +1,10 @@
+module BatsrusMakieExt
+
+using Batsrus, Printf
+import Batsrus: findindex, hasunit, getunit, getdata
+import Batsrus.UnitfulBatsrus
+import Makie
+
+include("typerecipe.jl")
+
+end
\ No newline at end of file
diff --git a/ext/BatsrusMakieExt/typerecipe.jl b/ext/BatsrusMakieExt/typerecipe.jl
new file mode 100644
index 0000000..1558e30
--- /dev/null
+++ b/ext/BatsrusMakieExt/typerecipe.jl
@@ -0,0 +1,35 @@
+# Type conversion from Batsrus to Makie
+
+"Conversion for 1D plots"
+function Makie.convert_arguments(P::Makie.PointBased, bd::BATLData, var::String)
+   var_ = findindex(bd, var)
+   if hasunit(bd)
+      unitx = getunit(bd, bd.head.variables[1])
+      unitw = getunit(bd, var)
+      x = bd.x .* unitx
+      y = bd.w[:,var_] .* unitw
+   else
+      x = bd.x
+      y = bd.w[:,var_]
+   end
+
+   ([Makie.Point2f(i, j) for (i, j) in zip(x, y)],)
+end
+
+"Conversion for 2D plots."
+function Makie.convert_arguments(P::Makie.SurfaceLike, bd::BATLData, var::String;
+   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1)
+   x, y, w = getdata(bd, var, plotrange, plotinterval)
+
+   unitx = getunit(bd, bd.head.variables[1])
+   unity = getunit(bd, bd.head.variables[2])
+   unitw = getunit(bd, var)
+
+   if unitx isa UnitfulBatsrus.Unitlike
+      x *= unitx
+      y *= unity
+      w *= unitw
+   end
+
+   (x, y, w)
+end
diff --git a/src/Batsrus.jl b/src/Batsrus.jl
index 77d2103..1758280 100644
--- a/src/Batsrus.jl
+++ b/src/Batsrus.jl
@@ -23,7 +23,7 @@ struct FileList
    lenhead::Int
 end
 
-"Primary BATLData storage type"
+"Primary Batsrus data storage type."
 struct BATLData{T<:AbstractFloat}
    "header information"
    head::NamedTuple
diff --git a/src/io.jl b/src/io.jl
index eedbb3e..d12213a 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -695,15 +695,15 @@ end
 function Base.show(io::IO, data::BATLData)
    showhead(io, data)
    if data.list.bytes ≥ 1e9
-      println(io, "filesize = $(data.list.bytes/1e9) GB")
+      println(io, "filesize: $(data.list.bytes/1e9) GB")
    elseif data.list.bytes ≥ 1e6
-      println(io, "filesize = $(data.list.bytes/1e6) MB")
+      println(io, "filesize: $(data.list.bytes/1e6) MB")
    elseif data.list.bytes ≥ 1e3
-      println(io, "filesize = $(data.list.bytes/1e3) KB")
+      println(io, "filesize: $(data.list.bytes/1e3) KB")
    else
-      println(io, "filesize = $(data.list.bytes) bytes")
+      println(io, "filesize: $(data.list.bytes) bytes")
    end
-   println(io, "snapshots = $(data.list.npictinfiles)")
+   println(io, "snapshots: $(data.list.npictinfiles)")
 end
 
 """
diff --git a/src/plot/makie.jl b/src/plot/makie.jl
deleted file mode 100644
index 46d6827..0000000
--- a/src/plot/makie.jl
+++ /dev/null
@@ -1,40 +0,0 @@
-# Makie recipe
-
-using Batsrus
-using Makie
-
-include("utility.jl")
-
-# 1D example
-filename = "1d__raw_2_t25.60000_n00000258.out"
-data = load(filename, dir="test/data")
-
-var = "p"
-VarIndex_ = findindex(data, var)
-
-fig = Figure(resolution = (800, 600))
-ax = Axis(fig[1,1], title="1D shocktube", xlabel="x", ylabel="Pressure [nPa]")
-line1 = lines!(ax, data.x[:], data.w[:,VarIndex_], color = :red,
-   label="Pressure")
-axislegend()
-
-fig
-#=
-# 2D example
-#filename = "z=0_raw_1_t25.60000_n00000258.out"
-#data = load(filename, dir="test/data")
-
-var = "p"
-plotrange = [-Inf,Inf,-Inf,Inf]
-plotinterval = 0.1
-
-fig = Figure(resolution = (800, 600))
-
-xi, yi, wi = getdata(data, var, plotrange, plotinterval)
-
-fig[1,1] = Axis(fig, title="2D shocktube")
-c = heatmap!(fig[1,1], xi, yi, wi)
-cbar = Colorbar(fig[1,2], c, label="Pressure [nPa]", width=20)
-
-fig
-=#
\ No newline at end of file
diff --git a/src/plot/plots.jl b/src/plot/plots.jl
index e6978b2..398a3f0 100644
--- a/src/plot/plots.jl
+++ b/src/plot/plots.jl
@@ -8,17 +8,13 @@ using RecipesBase
 
    ndim = data.head.ndim
 
-   if startswith(data.head.headline, "normalized")
-      hasunits = false
-   else
-      hasunits = true
-      unitw = getunit(data, var)
-   end
+   hasunits = hasunit(bd)
 
    if ndim == 1
       VarIndex_ = findindex(data, var)
       if hasunits
          unitx = getunit(data, data.head.variables[1])
+         unitw = getunit(data, var)
          x = data.x .* unitx
          w = data.w
          y = w[:,VarIndex_] .* unitw
@@ -36,6 +32,7 @@ using RecipesBase
 
       unitx = getunit(data, data.head.variables[1])
       unity = getunit(data, data.head.variables[2])
+      unitw = getunit(data, var)
 
       if unitx isa UnitfulBatsrus.Unitlike
          x *= unitx
diff --git a/src/plot/pyplot.jl b/src/plot/pyplot.jl
index 447fef8..4f55f9a 100644
--- a/src/plot/pyplot.jl
+++ b/src/plot/pyplot.jl
@@ -53,22 +53,24 @@ Plot the variable from SWMF output.
 
 `plotdata(data, "p", plotmode="grid")`
 
-`plotdata(data, func, plotmode="trimesh",plotrange=plotrange, plotinterval=0.2)`
+`plotdata(data, func, plotmode="trimesh", plotrange=[-1.0, 1.0, -1.0, 1.0], plotinterval=0.2)`
 
-# Input arguments
+# Arguments
 - `data::BATLData`: BATSRUS data to be visualized.
 - `func::String`: variables for plotting.
-- `plotmode::String`: (optional) type of plotting ["cont","contbar"]...
-- `plotrange::Vector`: (optional) range of plotting.
-- `plotinterval`: (optional) interval for interpolation.
-- `level`: (optional) level of contour.
-- `innermask`: (optional) Bool for masking a circle at the inner boundary.
-- `dir`: (optional) 2D cut plane orientation from 3D outputs ["x","y","z"].
-- `sequence`: (optional) sequence of plane from - to + in that direction.
-- `multifigure`: (optional) 1 for multifigure display, 0 for subplots.
-- `verbose`: (optional) display additional information.
-- `density`: (optional) density for streamlines.
-- `stride`: (optional) quiver strides in number of cells.
+
+# Keywords
+- `plotmode::String`: type of plotting ["cont","contbar"]...
+- `plotrange::Vector`: range of plotting.
+- `plotinterval`: interval for interpolation.
+- `level`: level of contour.
+- `innermask`: Bool for masking a circle at the inner boundary.
+- `dir`: 2D cut plane orientation from 3D outputs ["x","y","z"].
+- `sequence`: sequence of plane from - to + in that direction.
+- `multifigure`: 1 for multifigure display, 0 for subplots.
+- `verbose`: display additional information.
+- `density`: density for streamlines.
+- `stride`: quiver strides in number of cells.
 Right now this can only deal with 2D plots or 3D cuts. Full 3D plots may be supported in the
 future.
 """
diff --git a/src/plot/utility.jl b/src/plot/utility.jl
index 37d03f6..8f643b0 100644
--- a/src/plot/utility.jl
+++ b/src/plot/utility.jl
@@ -84,4 +84,12 @@ function meshgrid(x, y)
    Y = [y for y in y, _ in x]
 
    X, Y
+end
+
+@inline function hasunit(bd::BATLData)
+   if startswith(bd.head.headline, "normalized")
+      return false
+   else
+      return true
+   end
 end
\ No newline at end of file

From 06cfe2be2c61336e9f18327208404d5414cc91fa Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Sun, 19 Nov 2023 21:22:49 -0500
Subject: [PATCH 4/6] Test for PyPlot, Plots, and Makie

---
 Project.toml                          |  10 +--
 src/plot/plots.jl                     |  29 ++++---
 src/plot/pyplot.jl                    |  28 +++----
 src/plot/utility.jl                   |   4 +-
 Artifacts.toml => test/Artifacts.toml |   0
 test/Project.toml                     |   8 ++
 test/runtests.jl                      | 110 ++++++++++++++++----------
 7 files changed, 108 insertions(+), 81 deletions(-)
 rename Artifacts.toml => test/Artifacts.toml (100%)
 create mode 100644 test/Project.toml

diff --git a/Project.toml b/Project.toml
index 12dcfd3..474813b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -31,12 +31,4 @@ RecipesBase = "1.1"
 Requires = "1.1"
 Unitful = "1.7"
 WriteVTK = "1.9"
-julia = "1.6"
-
-[extras]
-LazyArtifacts = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
-SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-
-[targets]
-test = ["Test", "LazyArtifacts", "SHA"]
+julia = "1.6"
\ No newline at end of file
diff --git a/src/plot/plots.jl b/src/plot/plots.jl
index 398a3f0..4933324 100644
--- a/src/plot/plots.jl
+++ b/src/plot/plots.jl
@@ -3,24 +3,23 @@
 using RecipesBase
 
 # Build a recipe which acts on a custom type.
-@recipe function f(data::BATLData, var::AbstractString;
+@recipe function f(bd::BATLData, var::AbstractString;
    plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1)
 
-   ndim = data.head.ndim
+   ndim = bd.head.ndim
 
    hasunits = hasunit(bd)
 
    if ndim == 1
-      VarIndex_ = findindex(data, var)
+      VarIndex_ = findindex(bd, var)
       if hasunits
-         unitx = getunit(data, data.head.variables[1])
-         unitw = getunit(data, var)
-         x = data.x .* unitx
-         w = data.w
-         y = w[:,VarIndex_] .* unitw
+         unitx = getunit(bd, bd.head.variables[1])
+         unitw = getunit(bd, var)
+         x = bd.x .* unitx
+         y = bd.w[:,VarIndex_] .* unitw
       else
-         x, w = data.x, data.w
-         y = w[:,VarIndex_]
+         x = bd.x
+         y = @view bd.w[:,VarIndex_]
       end
 
       @series begin
@@ -28,11 +27,11 @@ using RecipesBase
          x, y
       end
    elseif ndim == 2
-      x, y, w = getdata(data, var, plotrange, plotinterval)
+      x, y, w = getdata(bd, var, plotrange, plotinterval)
 
-      unitx = getunit(data, data.head.variables[1])
-      unity = getunit(data, data.head.variables[2])
-      unitw = getunit(data, var)
+      unitx = getunit(bd, bd.head.variables[1])
+      unity = getunit(bd, bd.head.variables[2])
+      unitw = getunit(bd, var)
 
       if unitx isa UnitfulBatsrus.Unitlike
          x *= unitx
@@ -42,7 +41,7 @@ using RecipesBase
 
       @series begin
          seriestype --> :contourf  # use := if you want to force it
-         x, y, w
+         x, y, w'
       end
    end
 end
\ No newline at end of file
diff --git a/src/plot/pyplot.jl b/src/plot/pyplot.jl
index 4f55f9a..5a27c3f 100644
--- a/src/plot/pyplot.jl
+++ b/src/plot/pyplot.jl
@@ -149,14 +149,14 @@ function plotdata(data::BATLData, func::AbstractString; dir="x", plotmode="contb
             # More robust method needed!
             if plotmode[ivar] ∈ ["contbar", "contbarlog"]
                if level == 0
-                  c = contourf(Xi, Yi, Wi; kwargs...)
+                  c = contourf(Xi, Yi, Wi'; kwargs...)
                else
-                  c = contourf(Xi, Yi, Wi, level; kwargs...)
+                  c = contourf(Xi, Yi, Wi', level; kwargs...)
                end
             elseif plotmode[ivar] ∈ ["cont", "contlog"]
-               c = contour(Xi, Yi, Wi; kwargs...)
+               c = contour(Xi, Yi, Wi'; kwargs...)
             elseif plotmode[ivar] ∈ ["surfbar", "surfbarlog"]
-               c = plot_surface(Xi, Yi, Wi; kwargs...)
+               c = plot_surface(Xi, Yi, Wi'; kwargs...)
             end
 
             occursin("bar", plotmode[ivar]) && colorbar()
@@ -182,9 +182,9 @@ function plotdata(data::BATLData, func::AbstractString; dir="x", plotmode="contb
                triang = matplotlib.tri.Triangulation(X, Y)
                c = ax.triplot(triang)
             elseif plotmode[ivar] == "trisurf"
-               c = ax.plot_trisurf(X, Y, W; kwargs...)
+               c = ax.plot_trisurf(X, Y, W'; kwargs...)
             elseif plotmode[ivar] == "tricont"
-               c = ax.tricontourf(X, Y, W; kwargs...)
+               c = ax.tricontourf(X, Y, W'; kwargs...)
                fig.colorbar(c,ax=ax)
             elseif plotmode[ivar] == "tristream"
                throw(ArgumentError("tristream not yet implemented!"))
@@ -352,9 +352,9 @@ function plotdata(data::BATLData, func::AbstractString; dir="x", plotmode="contb
          if plotmode[ivar] ∈ ("surf", "surfbar", "surfbarlog", "cont", "contbar", "contlog",
                "contbarlog")
             if level == 0
-               c = ax.contourf(cut1, cut2, W; kwargs...)
+               c = ax.contourf(cut1, cut2, W'; kwargs...)
             else
-               c = ax.contourf(cut1, cut2, W, level; kwargs...)
+               c = ax.contourf(cut1, cut2, W', level; kwargs...)
             end
             fig.colorbar(c, ax=ax)
             title(data.head.wnames[varIndex_])
@@ -409,7 +409,7 @@ function cutplot(data::BATLData, var::AbstractString, ax=nothing;
    Y = @view x[:,:,:,2]
    Z = @view x[:,:,:,3]
 
-   W = w[:,:,:,varIndex_]
+   W = @view w[:,:,:,varIndex_]
 
    if dir == "x"
       cut1 = @view X[sequence,:,:]
@@ -549,9 +549,9 @@ function PyPlot.contour(data::BATLData, var::AbstractString, levels=0; ax=nothin
    if isnothing(ax) ax = plt.gca() end
 
    if levels != 0
-      c = ax.contour(Xi, Yi, Wi, levels; kwargs...)
+      c = ax.contour(Xi, Yi, Wi', levels; kwargs...)
    else
-      c = ax.contour(Xi, Yi, Wi; kwargs...)
+      c = ax.contour(Xi, Yi, Wi'; kwargs...)
    end
 end
 
@@ -569,9 +569,9 @@ function PyPlot.contourf(data::BATLData, var::AbstractString, levels=0; ax=nothi
    if isnothing(ax) ax = plt.gca() end
 
    if levels != 0
-      c = ax.contourf(Xi, Yi, Wi, levels; kwargs...)
+      c = ax.contourf(Xi, Yi, Wi', levels; kwargs...)
    else
-      c = ax.contourf(Xi, Yi, Wi; kwargs...)
+      c = ax.contourf(Xi, Yi, Wi'; kwargs...)
    end
 end
 
@@ -645,7 +645,7 @@ function PyPlot.plot_surface(data::BATLData, var::AbstractString;
 
    Xi, Yi, Wi = getdata(data, var, plotrange, plotinterval; griddim=2, innermask)
 
-   plot_surface(Xi, Yi, Wi; kwargs...)
+   plot_surface(Xi, Yi, Wi'; kwargs...)
 end
 
 """
diff --git a/src/plot/utility.jl b/src/plot/utility.jl
index 8f643b0..be6d07e 100644
--- a/src/plot/utility.jl
+++ b/src/plot/utility.jl
@@ -27,14 +27,14 @@ function getdata(data::BATLData, var::AbstractString, plotrange, plotinterval; g
       triang = @views matplotlib.tri.Triangulation(X[:], Y[:])
       interpolator = @views matplotlib.tri.LinearTriInterpolator(triang, W[:])
       Xi, Yi = meshgrid(xi, yi)
-      Wi = interpolator(Xi, Yi)
+      Wi = interpolator(Xi, Yi)'
    else # Cartesian coordinates
       xrange = range(x[1,1,1], x[end,1,1], length=size(x,1))
       yrange = range(x[1,1,2], x[1,end,2], length=size(x,2))
       if all(isinf.(plotrange))
          xi, yi = xrange, yrange
          Xi, Yi = meshgrid(xi, yi)
-         Wi = w[:,:,varIndex_]'
+         Wi = w[:,:,varIndex_] # Matplotlib does not accept view!
       else
          if plotrange[1] == -Inf plotrange[1] = minimum(xrange) end
          if plotrange[2] ==  Inf plotrange[2] = maximum(xrange) end
diff --git a/Artifacts.toml b/test/Artifacts.toml
similarity index 100%
rename from Artifacts.toml
rename to test/Artifacts.toml
diff --git a/test/Project.toml b/test/Project.toml
new file mode 100644
index 0000000..b10a922
--- /dev/null
+++ b/test/Project.toml
@@ -0,0 +1,8 @@
+[deps]
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+LazyArtifacts = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
+PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 5a26f4a..9f1c121 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,6 +2,10 @@
 
 using Batsrus, Test, SHA, LazyArtifacts
 using Batsrus.UnitfulBatsrus, Unitful
+using RecipesBase
+using CairoMakie
+using PyPlot
+ENV["MPLBACKEND"]="agg" # no GUI
 
 function filecmp(path1::AbstractString, path2::AbstractString)
    stat1, stat2 = stat(path1), stat(path2)
@@ -28,50 +32,50 @@ end
    datapath = artifact"testdata"
    @testset "Reading 1D ascii" begin
       file = "1d__raw_2_t25.60000_n00000258.out"
-      data = load(file, dir=datapath, verbose=true)
-      @test startswith(repr(data), "filename : 1d")
-      @test Batsrus.setunits(data.head, "NORMALIZED")
-      @test isa(data.head, NamedTuple)
-      @test extrema(data.x) == (-127.5, 127.5)
-      @test extrema(data.w) == (-0.79960780498, 1.9394335293)
+      bd = load(file, dir=datapath, verbose=true)
+      @test startswith(repr(bd), "filename : 1d")
+      @test Batsrus.setunits(bd.head, "NORMALIZED")
+      @test isa(bd.head, NamedTuple)
+      @test extrema(bd.x) == (-127.5, 127.5)
+      @test extrema(bd.w) == (-0.79960780498, 1.9394335293)
    end
 
    @testset "Reading 2D structured binary" begin
       file = "z=0_raw_1_t25.60000_n00000258.out"
-      data = load(file, dir=datapath)
-      @test data.head.time == 25.6f0
-      @test extrema(data.x) == (-127.5f0, 127.5f0)
-      @test extrema(data.w) == (-0.79985905f0, 1.9399388f0)
+      bd = load(file, dir=datapath)
+      @test bd.head.time == 25.6f0
+      @test extrema(bd.x) == (-127.5f0, 127.5f0)
+      @test extrema(bd.w) == (-0.79985905f0, 1.9399388f0)
    end
 
    @testset "Reading 2D unstructured binary" begin
       #file = "z=0_raw_1_t25.60000_n00000258.out"
-      #data = load(file)
+      #bd = load(file)
    end
 
    @testset "Reading 3D structured binary" begin
       file = "3d_raw.out"
-      data = load(file, dir=datapath)
+      bd = load(file, dir=datapath)
       plotrange = [-50.0, 50.0, -0.5, 0.5]
-      X, Z, p = cutdata(data, "p"; dir="y", sequence=1, plotrange)
+      X, Z, p = cutdata(bd, "p"; dir="y", sequence=1, plotrange)
       @test p[1] ≈ 0.560976f0
       @test p[2] ≈ 0.53704995f0
-      vars = getvars(data, ["p"]) 
+      vars = getvars(bd, ["p"]) 
       @test size(vars["p"]) == (8,8,8) 
    end
 
    @testset "Log" begin
       logfilename = joinpath(datapath, "log_n000001.log")
-      head, data = readlogdata(logfilename)
+      head, bd = readlogdata(logfilename)
       @test isa(head, NamedTuple)
-      @test extrema(data) == (-0.105, 258.0)
+      @test extrema(bd) == (-0.105, 258.0)
    end
 
    @testset "VTK" begin
       file = joinpath(datapath, "3d_bin.dat")
-      head, data, connectivity = readtecdata(file)
+      head, bd, connectivity = readtecdata(file)
       @test maximum(connectivity) ≤ head[:nNode] # check if it's read correctly
-      convertTECtoVTU(head, data, connectivity)
+      convertTECtoVTU(head, bd, connectivity)
       sha_str = bytes2hex(open(sha1, "out.vtu"))
       @test sha_str == "5b04747666542d802357dec183177f757754a254"
       rm("out.vtu")
@@ -87,43 +91,67 @@ end
    end
 
    @testset "Plotting" begin
-      using PyPlot
-      ENV["MPLBACKEND"]="agg" # no GUI
+      @testset "Plots" begin
+         RecipesBase.is_key_supported(k::Symbol) = true
+         # 1D
+         file = "1d__raw_2_t25.60000_n00000258.out"
+         bd = load(file, dir=datapath)
+         rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), bd, "Rho")
+         @test getfield(rec[1], 1)[:seriestype] == :path
+   
+         file = "z=0_raw_1_t25.60000_n00000258.out"
+         bd = load(file, dir=datapath)
+         rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), bd, "p")
+         @test getfield(rec[1], 1)[:seriestype] == :contourf
+      end
+
+      @testset "Makie" begin
+         file = "1d__raw_2_t25.60000_n00000258.out"
+         bd = load(file, dir=datapath)
+         fig, ax, plt = lines(bd, "Rho")
+         @test plt isa Lines
+   
+         file = "z=0_raw_1_t25.60000_n00000258.out"
+         bd = load(file, dir=datapath)
+         fig, ax, plt = heatmap(bd, "p")
+         @test plt isa Heatmap
+      end
+
       @testset "1D ascii" begin
          file = "1d__raw_2_t25.60000_n00000258.out"
-         data = load(file, dir=datapath, verbose=false)
-         plotdata(data, "p", plotmode="line")
+         bd = load(file, dir=datapath, verbose=false)
+         plotdata(bd, "p", plotmode="line")
          line = get(gca().lines, 0)
-         @test line.get_xdata() ≈ data.x
-         @test line.get_ydata() ≈ data.w[:,10]
+         @test line.get_xdata() ≈ bd.x
+         @test line.get_ydata() ≈ bd.w[:,10]
       end
 
       @testset "2D structured binary" begin
          file = "z=0_raw_1_t25.60000_n00000258.out"
-         data = load(file, dir=datapath)
-         plotdata(data, "p bx;by", plotmode="contbar streamover")
+         bd = load(file, dir=datapath)
+         plotdata(bd, "p bx;by", plotmode="contbar streamover")
          @test isa(gca(), PyPlot.PyObject)
-         contourf(data, "p")
+         PyPlot.contourf(bd, "p")
          @test isa(gca(), PyPlot.PyObject)
-         @test_throws ErrorException contourf(data, "rho", innermask=true)
-         contour(data, "rho")
+         @test_throws ErrorException PyPlot.contourf(bd, "rho", innermask=true)
+         PyPlot.contour(bd, "rho")
          @test isa(gca(), PyPlot.PyObject)
-         tricontourf(data, "rho")
+         PyPlot.tricontourf(bd, "rho")
          @test isa(gca(), PyPlot.PyObject)
-         streamplot(data, "bx;by")
+         PyPlot.streamplot(bd, "bx;by")
          @test isa(gca(), PyPlot.PyObject)
          plt.close()
          fig = plt.figure()
          ax = fig.add_subplot(111, projection="3d")
-         plot_surface(data, "rho")
+         plot_surface(bd, "rho")
          @test isa(gca(), PyPlot.PyObject)
          plt.close()
       end
 
       @testset "2D AMR Cartesian" begin
          file = "bx0_mhd_6_t00000100_n00000352.out"
-         data = load(file, dir=datapath)
-         plotdata(data, "P", plotmode="contbar")
+         bd = load(file, dir=datapath)
+         plotdata(bd, "P", plotmode="contbar")
          ax = gca()
          @test isa(ax, PyPlot.PyObject)
       end
@@ -132,18 +160,18 @@ end
    @testset "Units" begin
       @test 1.0bu"R" > 2.0bu"Rg"
       file = "y=0_var_1_t00000000_n00000000.out"
-      data = load(file, dir=datapath)
-      varunit = getunit(data, "Rho")
+      bd = load(file, dir=datapath)
+      varunit = getunit(bd, "Rho")
       @test varunit == bu"amucc"
-      varunit = getunit(data, "Ux")
+      varunit = getunit(bd, "Ux")
       @test varunit == u"km/s"
-      varunit = getunit(data, "Bx")
+      varunit = getunit(bd, "Bx")
       @test varunit == u"nT"
-      varunit = getunit(data, "P")
+      varunit = getunit(bd, "P")
       @test varunit == u"nPa"
-      varunit = getunit(data, "jx")
+      varunit = getunit(bd, "jx")
       @test dimension(varunit) == dimension(Unitful.A/Unitful.m^2)
-      varunit = getunit(data, "ex")
+      varunit = getunit(bd, "ex")
       @test varunit == u"mV/m"
    end
 

From ef77cb3347eace99821b0ffb4222c2d05d520887 Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Sun, 19 Nov 2023 21:31:08 -0500
Subject: [PATCH 5/6] Add whitespace

---
 src/plot/pyplot.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/plot/pyplot.jl b/src/plot/pyplot.jl
index 5a27c3f..d39ed92 100644
--- a/src/plot/pyplot.jl
+++ b/src/plot/pyplot.jl
@@ -32,7 +32,7 @@ function plotlogdata(data, head::NamedTuple, func::AbstractString; plotmode="lin
       if plotmode[ivar] == "line"
          plot(data[1,:],data[varIndex_,:])
       elseif plotmode[ivar] == "scatter"
-         scatter(data[1,:],data[varIndex_,:])
+         scatter(data[1,:], data[varIndex_,:])
       else
          throw(ArgumentError("unknown plot mode $(plotmode[ivar])!"))
       end

From a1dcc042d265b551fd795dd766319341d7f54fa4 Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Sun, 19 Nov 2023 21:35:45 -0500
Subject: [PATCH 6/6] Rename variables

---
 src/plot/pyplot.jl | 116 ++++++++++++++++++++++-----------------------
 1 file changed, 58 insertions(+), 58 deletions(-)

diff --git a/src/plot/pyplot.jl b/src/plot/pyplot.jl
index d39ed92..4fdff64 100644
--- a/src/plot/pyplot.jl
+++ b/src/plot/pyplot.jl
@@ -56,7 +56,7 @@ Plot the variable from SWMF output.
 `plotdata(data, func, plotmode="trimesh", plotrange=[-1.0, 1.0, -1.0, 1.0], plotinterval=0.2)`
 
 # Arguments
-- `data::BATLData`: BATSRUS data to be visualized.
+- `bd::BATLData`: BATSRUS data to be visualized.
 - `func::String`: variables for plotting.
 
 # Keywords
@@ -74,26 +74,26 @@ Plot the variable from SWMF output.
 Right now this can only deal with 2D plots or 3D cuts. Full 3D plots may be supported in the
 future.
 """
-function plotdata(data::BATLData, func::AbstractString; dir="x", plotmode="contbar",
+function plotdata(bd::BATLData, func::AbstractString; dir="x", plotmode="contbar",
    plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, sequence=1, multifigure=true,
    getrangeonly=false, level=0, innermask=false, verbose=false, density=1.0,
    stride=10, kwargs...)
 
-   x, w = data.x, data.w
+   x, w = bd.x, bd.w
    plotmode = split(plotmode)
    vars     = split(func)
-   ndim     = data.head.ndim
+   ndim     = bd.head.ndim
    nvar     = length(vars)
 
    if verbose || getrangeonly
       @info "============ PLOTTING PARAMETERS ==============="
-      @info "wnames = $(data.head.wnames)"
+      @info "wnames = $(bd.head.wnames)"
       wmin = Vector{Float64}(undef,nvar)
       wmax = Vector{Float64}(undef,nvar)
       # Display min & max for each variable
       for (ivar,var) in enumerate(vars)
          if occursin(";",var) continue end # skip the vars for streamline
-         varIndex_ = findindex(data, var)
+         varIndex_ = findindex(bd, var)
          if ndim == 1
             wmin[ivar] = minimum(w[:,varIndex_])
             wmax[ivar] = maximum(w[:,varIndex_])
@@ -114,7 +114,7 @@ function plotdata(data::BATLData, func::AbstractString; dir="x", plotmode="contb
    ## Plot
    if ndim == 1
       for (ivar,var) in enumerate(vars)
-         varIndex_ = findindex(data, var)
+         varIndex_ = findindex(bd, var)
          if ivar == 1 || multifigure fig, ax = subplots() else ax = gca() end
          if !occursin("scatter", plotmode[ivar])
             plot(x, w[:,varIndex_]; kwargs...)
@@ -126,7 +126,7 @@ function plotdata(data::BATLData, func::AbstractString; dir="x", plotmode="contb
          end
          xlabel("x"); ylabel("$(var)")
          dim = [0.125, 0.013, 0.2, 0.045]
-         str = @sprintf "it=%d, time=%4.2f" data.head.it data.head.time
+         str = @sprintf "it=%d, time=%4.2f" bd.head.it bd.head.time
          at = matplotlib.offsetbox.AnchoredText(str,
             loc="lower left", prop=Dict("size"=>8), frameon=true,
             bbox_to_anchor=(0., 1.), bbox_transform=ax.transAxes)
@@ -138,13 +138,13 @@ function plotdata(data::BATLData, func::AbstractString; dir="x", plotmode="contb
          occursin("over", plotmode[ivar]) && (multifigure = false)
          if ivar == 1 || multifigure fig, ax = subplots() else ax = gca() end
          if !occursin(";",var)
-            varIndex_ = findindex(data, var)
+            varIndex_ = findindex(bd, var)
          end
 
          if plotmode[ivar] ∈ ("surf","surfbar","surfbarlog","cont","contbar",
             "contlog","contbarlog")
 
-            Xi, Yi, Wi = getdata(data, var, plotrange, plotinterval; griddim=2, innermask)
+            Xi, Yi, Wi = getdata(bd, var, plotrange, plotinterval; griddim=2, innermask)
 
             # More robust method needed!
             if plotmode[ivar] ∈ ["contbar", "contbarlog"]
@@ -162,7 +162,7 @@ function plotdata(data::BATLData, func::AbstractString; dir="x", plotmode="contb
             occursin("bar", plotmode[ivar]) && colorbar()
             occursin("log", plotmode[ivar]) &&
             ( c.locator = matplotlib.ticker.LogLocator() )
-            title(data.head.wnames[varIndex_])
+            title(bd.head.wnames[varIndex_])
 
          elseif plotmode[ivar] ∈ ("trimesh","trisurf","tricont","tristream")
             X = vec(x[:,:,1])
@@ -190,14 +190,14 @@ function plotdata(data::BATLData, func::AbstractString; dir="x", plotmode="contb
                throw(ArgumentError("tristream not yet implemented!"))
             end
 
-            title(data.head.wnames[varIndex_])
+            title(bd.head.wnames[varIndex_])
 
          elseif plotmode[ivar] ∈ ("stream","streamover")
             VarStream  = split(var,";")
-            VarIndex1_ = findindex(data, VarStream[1])
-            VarIndex2_ = findindex(data, VarStream[2])
+            VarIndex1_ = findindex(bd, VarStream[1])
+            VarIndex2_ = findindex(bd, VarStream[2])
 
-            if data.head.gencoord # Generalized coordinates
+            if bd.head.gencoord # Generalized coordinates
                xrange = @view x[:,:,1]
                yrange = @view x[:,:,2]
                if any(isinf.(plotrange))
@@ -258,8 +258,8 @@ function plotdata(data::BATLData, func::AbstractString; dir="x", plotmode="contb
 
          elseif occursin("quiver", plotmode[ivar])
             VarQuiver  = split(var, ";")
-            VarIndex1_ = findindex(data, VarQuiver[1])
-            VarIndex2_ = findindex(data, VarQuiver[2])
+            VarIndex1_ = findindex(bd, VarQuiver[1])
+            VarIndex2_ = findindex(bd, VarQuiver[2])
 
             X, Y = x[:,1,1], x[1,:,2]
             v1, v2 = w[:,:,VarIndex1_]', w[:,:,VarIndex2_]'
@@ -279,9 +279,9 @@ function plotdata(data::BATLData, func::AbstractString; dir="x", plotmode="contb
             throw(ArgumentError("unknown plot mode: $(plotmode[ivar])"))
          end
 
-         xlabel(data.head.variables[1]); ylabel(data.head.variables[2])
+         xlabel(bd.head.variables[1]); ylabel(bd.head.variables[2])
          dim = [0.125, 0.013, 0.2, 0.045]
-         str = @sprintf "it=%d, time=%4.2f" data.head.it data.head.time
+         str = @sprintf "it=%d, time=%4.2f" bd.head.it bd.head.time
          at = matplotlib.offsetbox.AnchoredText(str,
             loc="lower left", prop=Dict("size"=>8), frameon=true,
             bbox_to_anchor=(0., 1.), bbox_transform=ax.transAxes)
@@ -299,7 +299,7 @@ function plotdata(data::BATLData, func::AbstractString; dir="x", plotmode="contb
          if plotmode[ivar] ∈ ("surf","surfbar","surfbarlog","cont","contbar", "contlog",
                "contbarlog")
 
-            varIndex_ = findindex(data, var)
+            varIndex_ = findindex(bd, var)
 
             if ivar == 1 || multifigure fig, ax = subplots() else ax = gca() end
 
@@ -320,8 +320,8 @@ function plotdata(data::BATLData, func::AbstractString; dir="x", plotmode="contb
             end
          elseif plotmode[ivar] ∈ ("stream","streamover")
             varStream  = split(var,";")
-            VarIndex1_ = findindex(data, varStream[1])
-            VarIndex2_ = findindex(data, varStream[2])
+            VarIndex1_ = findindex(bd, varStream[1])
+            VarIndex2_ = findindex(bd, varStream[2])
 
             v1 = @view w[:,:,:,VarIndex1_]
             v2 = @view w[:,:,:,VarIndex2_]
@@ -357,7 +357,7 @@ function plotdata(data::BATLData, func::AbstractString; dir="x", plotmode="contb
                c = ax.contourf(cut1, cut2, W', level; kwargs...)
             end
             fig.colorbar(c, ax=ax)
-            title(data.head.wnames[varIndex_])
+            title(bd.head.wnames[varIndex_])
 
          elseif plotmode[ivar] ∈ ("stream", "streamover")
             xi = range(cut1[1,1], stop=cut1[1,end],
@@ -380,7 +380,7 @@ function plotdata(data::BATLData, func::AbstractString; dir="x", plotmode="contb
 
          ax = gca()
          dim = [0.125, 0.013, 0.2, 0.045]
-         str = @sprintf "it=%d, time=%4.2f" data.head.it data.head.time
+         str = @sprintf "it=%d, time=%4.2f" bd.head.it bd.head.time
          at = matplotlib.offsetbox.AnchoredText(str,
             loc="lower left", prop=Dict("size"=>8), frameon=true,
             bbox_to_anchor=(0., 1.), bbox_transform=ax.transAxes)
@@ -399,11 +399,11 @@ end
 
 2D plane cut contourf of 3D box data.
 """
-function cutplot(data::BATLData, var::AbstractString, ax=nothing;
+function cutplot(bd::BATLData, var::AbstractString, ax=nothing;
    plotrange=[-Inf,Inf,-Inf,Inf], dir="x", sequence=1, level=20)
 
-   x, w = data.x, data.w
-   varIndex_ = findindex(data, var)
+   x, w = bd.x, bd.w
+   varIndex_ = findindex(bd, var)
 
    X = @view x[:,:,:,1]
    Y = @view x[:,:,:,2]
@@ -431,7 +431,7 @@ function cutplot(data::BATLData, var::AbstractString, ax=nothing;
    if isnothing(ax) ax = plt.gca() end
    c = ax.contourf(cut1, cut2, W, level)
 
-   title(data.head.wnames[varIndex_])
+   title(bd.head.wnames[varIndex_])
 
    if dir == "x"
       xlabel("y"); ylabel("z")
@@ -452,13 +452,13 @@ end
 Plot streamlines on 2D slices of 3D box data. Variable names in `var` string must be
 separated with `;`.
 """
-function streamslice(data::BATLData, var::AbstractString, ax=nothing;
+function streamslice(bd::BATLData, var::AbstractString, ax=nothing;
    plotrange=[-Inf,Inf,-Inf,Inf], dir="x", sequence=1, kwargs...)
 
-   x,w = data.x, data.w
+   x,w = bd.x, bd.w
    varStream  = split(var, ";")
-   VarIndex1_ = findindex(data, varStream[1])
-   VarIndex2_ = findindex(data, varStream[2])
+   VarIndex1_ = findindex(bd, varStream[1])
+   VarIndex2_ = findindex(bd, varStream[2])
 
    X = @view x[:,:,:,1]
    Y = @view x[:,:,:,2]
@@ -514,9 +514,9 @@ end
 
 Wrapper over `plot` in matplotlib.
 """
-function PyPlot.plot(data::BATLData, var::AbstractString, ax=nothing; kwargs...)
-   x, w = data.x, data.w
-   varIndex_ = findindex(data, var)
+function PyPlot.plot(bd::BATLData, var::AbstractString, ax=nothing; kwargs...)
+   x, w = bd.x, bd.w
+   varIndex_ = findindex(bd, var)
    if isnothing(ax) ax = plt.gca() end
 
    c = ax.plot(x, w[:,varIndex_]; kwargs...)
@@ -527,9 +527,9 @@ end
 
 Wrapper over `scatter` in matplotlib.
 """
-function PyPlot.scatter(data::BATLData, var::AbstractString, ax=nothing; kwargs...)
-   x, w = data.x, data.w
-   varIndex_ = findindex(data, var)
+function PyPlot.scatter(bd::BATLData, var::AbstractString, ax=nothing; kwargs...)
+   x, w = bd.x, bd.w
+   varIndex_ = findindex(bd, var)
    if isnothing(ax) ax = plt.gca() end
 
    c = ax.scatter(x, w[:,varIndex_]; kwargs...)
@@ -541,10 +541,10 @@ end
 
 Wrapper over `contour` in matplotlib.
 """
-function PyPlot.contour(data::BATLData, var::AbstractString, levels=0; ax=nothing,
+function PyPlot.contour(bd::BATLData, var::AbstractString, levels=0; ax=nothing,
    plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, innermask=false, kwargs...)
 
-   Xi, Yi, Wi = getdata(data, var, plotrange, plotinterval; innermask)
+   Xi, Yi, Wi = getdata(bd, var, plotrange, plotinterval; innermask)
 
    if isnothing(ax) ax = plt.gca() end
 
@@ -561,10 +561,10 @@ end
 
 Wrapper over `contourf` in matplotlib.
 """
-function PyPlot.contourf(data::BATLData, var::AbstractString, levels=0; ax=nothing,
+function PyPlot.contourf(bd::BATLData, var::AbstractString, levels=0; ax=nothing,
    plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, innermask=false, kwargs...)
 
-   Xi, Yi, Wi = getdata(data, var, plotrange, plotinterval; innermask)
+   Xi, Yi, Wi = getdata(bd, var, plotrange, plotinterval; innermask)
 
    if isnothing(ax) ax = plt.gca() end
 
@@ -581,11 +581,11 @@ end
 
 Wrapper over `tricontourf` in matplotlib.
 """
-function PyPlot.tricontourf(data::BATLData, var::AbstractString, ax=nothing;
+function PyPlot.tricontourf(bd::BATLData, var::AbstractString, ax=nothing;
    plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, kwargs...)
 
-   x, w = data.x, data.w
-   varIndex_ = findindex(data, var)
+   x, w = bd.x, bd.w
+   varIndex_ = findindex(bd, var)
 
    X = vec(x[:,:,1])
    Y = vec(x[:,:,2])
@@ -610,11 +610,11 @@ end
 
 Wrapper over `plot_trisurf` in matplotlib.
 """
-function PyPlot.plot_trisurf(data::BATLData, var::AbstractString, ax=nothing;
+function PyPlot.plot_trisurf(bd::BATLData, var::AbstractString, ax=nothing;
    plotrange=[-Inf,Inf,-Inf,Inf], kwargs...)
 
-   x, w = data.x, data.w
-   varIndex_ = findindex(data, var)
+   x, w = bd.x, bd.w
+   varIndex_ = findindex(bd, var)
 
    X = vec(x[:,:,1])
    Y = vec(x[:,:,2])
@@ -640,10 +640,10 @@ end
 
 Wrapper over `plot_surface` in matplotlib.
 """
-function PyPlot.plot_surface(data::BATLData, var::AbstractString;
+function PyPlot.plot_surface(bd::BATLData, var::AbstractString;
    plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, innermask=false, kwargs...)
 
-   Xi, Yi, Wi = getdata(data, var, plotrange, plotinterval; griddim=2, innermask)
+   Xi, Yi, Wi = getdata(bd, var, plotrange, plotinterval; griddim=2, innermask)
 
    plot_surface(Xi, Yi, Wi'; kwargs...)
 end
@@ -655,16 +655,16 @@ end
 Wrapper over `streamplot` in matplotlib. `streamplot` does not have **kwargs in the API, but
 it supports `density`, `color`, and some other keywords.
 """
-function PyPlot.streamplot(data::BATLData, var::AbstractString, ax=nothing;
+function PyPlot.streamplot(bd::BATLData, var::AbstractString, ax=nothing;
    plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, kwargs...)
 
-   x, w = data.x, data.w
+   x, w = bd.x, bd.w
    VarStream  = split(var,";")
-   wnames = lowercase.(data.head.wnames)
+   wnames = lowercase.(bd.head.wnames)
    VarIndex1_ = findfirst(x->x==lowercase(VarStream[1]), wnames)
    VarIndex2_ = findfirst(x->x==lowercase(VarStream[2]), wnames)
 
-   if data.head.gencoord # generalized coordinates
+   if bd.head.gencoord # generalized coordinates
       X, Y = vec(x[:,:,1]), vec(x[:,:,2])
       if any(isinf.(plotrange))
          if plotrange[1] == -Inf plotrange[1] = minimum(X) end
@@ -725,12 +725,12 @@ end
 
 Wrapper over `quiver` in matplotlib. Only supports Cartesian grid for now.
 """
-function PyPlot.quiver(data::BATLData, var::AbstractString, ax=nothing;
+function PyPlot.quiver(bd::BATLData, var::AbstractString, ax=nothing;
    stride::Integer=10, kwargs...)
-   x, w = data.x, data.w
+   x, w = bd.x, bd.w
    VarQuiver  = split(var, ";")
-   VarIndex1_ = findindex(data, VarQuiver[1])
-   VarIndex2_ = findindex(data, VarQuiver[2])
+   VarIndex1_ = findindex(bd, VarQuiver[1])
+   VarIndex2_ = findindex(bd, VarQuiver[2])
 
    @views X, Y = x[:,1,1], x[1,:,2]
    @views v1, v2 = w[:,:,VarIndex1_]', w[:,:,VarIndex2_]'