diff --git a/Project.toml b/Project.toml
index e9576de3..10517668 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.3"
+version = "0.3.4"
 
 [deps]
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
@@ -30,9 +30,8 @@ WriteVTK = "1.9"
 julia = "1.6"
 
 [extras]
-PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "SHA", "PyPlot"]
+test = ["Test", "SHA"]
diff --git a/docs/src/man/examples.md b/docs/src/man/examples.md
index 3730b57c..b9d0982b 100644
--- a/docs/src/man/examples.md
+++ b/docs/src/man/examples.md
@@ -102,6 +102,8 @@ filenames = vcat(filenames, filesfound)
 end
 ```
 
+More examples can be found in [examples](https://github.com/henry2004y/Batsrus.jl/tree/master/examples).
+
 ## Data visualization
 
 We provide plot recipes for Plots.jl and wrappers for PyPlot.jl.
@@ -161,7 +163,7 @@ using PyPlot
 plt.axis("scaled")
 
 subplot(2,2,(1,3))
-cutplot(data, "Ex", cut='y', cutPlaneIndex=128, plotrange=plotrange)
+cutplot(data, "Ex", cut="y", cutPlaneIndex=128, plotrange=plotrange)
 ```
 
 #### Finding indexes
@@ -234,7 +236,7 @@ set_cmap("RdBu_r")
 
 contourf(data, "uxS0", 50, plotrange=[-3,3,-3,3], plotinterval=0.05, norm=DN(0))
 colorbar()
-streamplot(data, "uxS0;uzS0", density=2.0, color="g",plotrange=[-3,3,-3,3])
+streamplot(data, "uxS0;uzS0", density=2.0, color="g", plotrange=[-3,3,-3,3])
 xlabel("x"); ylabel("y"); title("Ux [km/s]")
 
 contourf(data,"uxS0", 50, plotinterval=0.05, norm=DN(0))
diff --git a/src/convert2vtk_distributed.jl b/examples/convert2vtk_distributed.jl
similarity index 88%
rename from src/convert2vtk_distributed.jl
rename to examples/convert2vtk_distributed.jl
index 2f869d82..2c8554bb 100644
--- a/src/convert2vtk_distributed.jl
+++ b/examples/convert2vtk_distributed.jl
@@ -10,11 +10,9 @@
 using Distributed
 using Pkg
 if "VisAna" ∈ keys(Pkg.installed())
-   @everywhere using VisAna
+   @everywhere using Batsrus
 else
-   @warn "VisAna not installed. Consider installing the package by
-   `Pkg.add(PackageSpec(url=\"https://github.com/henry2004y/VisAnaJulia\",
-    rev=\"master\"))`"
+   @warn "Batsrus.jl not installed. Install the package by `Pkg.add("Batsrus")`"
 end
 @everywhere using Glob
 
diff --git a/src/convert2vtk_threads.jl b/examples/convert2vtk_threads.jl
similarity index 89%
rename from src/convert2vtk_threads.jl
rename to examples/convert2vtk_threads.jl
index 9cbd235f..4ec6a034 100644
--- a/src/convert2vtk_threads.jl
+++ b/examples/convert2vtk_threads.jl
@@ -10,11 +10,9 @@
 
 using Pkg, Glob
 if "VisAna" ∈ keys(Pkg.installed())
-   using VisAna
+   using Batsrus
 else
-   @warn "VisAna not installed. Consider installing the package by
-   `Pkg.add(PackageSpec(url=\"https://github.com/henry2004y/VisAnaJulia\",
-    rev=\"master\"))`"
+   @warn "Batsrus.jl not installed. Install the package by `Pkg.add("Batsrus")`"
 end
 
 dir = "IO2"
diff --git a/src/convert2vtm.jl b/examples/convert2vtm.jl
similarity index 100%
rename from src/convert2vtm.jl
rename to examples/convert2vtm.jl
diff --git a/src/Batsrus.jl b/src/Batsrus.jl
index 0ac63775..e3f2feff 100644
--- a/src/Batsrus.jl
+++ b/src/Batsrus.jl
@@ -7,29 +7,27 @@ using Printf, Requires
 
 export Data
 
-"""
-	FileList
-
-Type for the file information. Contains filename `name`, file type `type` of
-ascii or binary, file size `bytes`, and number of snapshots `npictinfiles`.
-"""
+"Type for the file information."
 struct FileList
+   "filename"
    name::String
+   "file type"
    type::String
+   "file size"
    bytes::Int64
+   "number of snapshots"
    npictinfiles::Int64
 end
 
-"""
-    Data
-
-Primary data storage type, with fields `head` of header info, grid `x`, value
-`w`, and file info `list`.
-"""
+"Primary data storage type"
 struct Data{T<:AbstractFloat}
+   "header information"
    head::NamedTuple
+   "grid"
    x::Array{T}
+   "variables"
    w::Array{T}
+   "file information"
    list::FileList
 end
 
diff --git a/src/constants.jl b/src/constants.jl
deleted file mode 100644
index 57675020..00000000
--- a/src/constants.jl
+++ /dev/null
@@ -1,16 +0,0 @@
-# Commonly used constants
-
-const q = 1.6021765e-19 #[C]
-const dx = 1/100 #[Rg]
-const Rg = 2634000. #[m], Ganymede's radius
-const Rm = 2444000. #[m], Mercury's radius
-
-const ionElectronMassRatio = 100 # ion-electron mass ratio in PIC simulation
-const ionMass = 14 # single fluid ion mass in Ganymede's simulation
-
-const me = 9.10938356e-31 #[kg]
-const mp = 1.6726219e-27  #[kg]
-
-const c = 3e8 #[m/s]
-const Rₑ = 6.38e6 #[m]
-const μ = 4*π*1e-7  # Vacuum permeability, [V*s/(A*m)]
diff --git a/src/plot/pyplot.jl b/src/plot/pyplot.jl
index 66a33d8c..f3e944ea 100644
--- a/src/plot/pyplot.jl
+++ b/src/plot/pyplot.jl
@@ -3,11 +3,11 @@
 using PyPlot
 using Dierckx: Spline2D
 
-export plotdata, plotlogdata, animatedata, plot, scatter, contour, contourf,
-	   plot_surface, tricontourf, plot_trisurf, streamplot, streamslice, cutplot
+export plotdata, plotlogdata, animatedata, plot, scatter, contour, contourf, plot_surface,
+   tricontourf, plot_trisurf, streamplot, streamslice, cutplot
 
-import PyPlot: plot, scatter, contour, contourf, plot_surface, tricontourf,
-	   plot_trisurf, streamplot
+import PyPlot: plot, scatter, contour, contourf, plot_surface, tricontourf, plot_trisurf,
+   streamplot
 
 include("utility.jl")
 
@@ -70,12 +70,12 @@ Plot the variable from SWMF output.
 - `multifigure`: (optional) 1 for multifigure display, 0 for subplots.
 - `verbose`: (optional) display additional information.
 - `density`: (optional) density for streamlines.
-Right now this can only deal with 2D plots or 3D cuts. Full 3D plots may be
-supported in the future.
+Right now this can only deal with 2D plots or 3D cuts. Full 3D plots may be supported in the
+future.
 """
 function plotdata(data::Data, func::AbstractString; cut="", plotmode="contbar",
-   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, cutPlaneIndex=1,
-   multifigure=true, getrangeOnly=false, level=0, verbose=false, density=1.0)
+   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, cutPlaneIndex=1, multifigure=true,
+   getrangeOnly=false, level=0, verbose=false, density=1.0)
 
    x, w = data.x, data.w
    plotmode = split(plotmode)
@@ -126,9 +126,8 @@ function plotdata(data::Data, func::AbstractString; cut="", plotmode="contbar",
          dim = [0.125, 0.013, 0.2, 0.045]
          str = @sprintf "it=%d, time=%4.2f" data.head.it data.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)
+            loc="lower left", prop=Dict("size"=>8), frameon=true,
+            bbox_to_anchor=(0., 1.), bbox_transform=ax.transAxes)
          at.patch.set_boxstyle("round,pad=0.,rounding_size=0.2")
          ax.add_artist(at)
       end
@@ -171,7 +170,7 @@ function plotdata(data::Data, func::AbstractString; cut="", plotmode="contbar",
             # This needs to be modified!!!
             if !all(isinf.(plotrange))
                xyIndex = X .> plotrange[1] .& X .< plotrange[2] .&
-               Y .> plotrange[3] .& Y .< plotrange[4]
+                  Y .> plotrange[3] .& Y .< plotrange[4]
                X = X[xyIndex]
                Y = Y[xyIndex]
                W = W[xyIndex]
@@ -263,12 +262,12 @@ function plotdata(data::Data, func::AbstractString; cut="", plotmode="contbar",
             X, Y = x[:,1,1], x[1,:,2]
             v1, v2 = w[:,:,VarIndex1_]', w[:,:,VarIndex2_]'
 
-            q = quiver(X,Y,v1,v2,color="w")
+            q = quiver(X, Y, v1, v2, color="w")
 
          elseif occursin("grid", plotmode[ivar])
             # This does not take subdomain plot into account!
             X, Y = x[:,:,1], x[:,:,2]
-            scatter(X,Y,marker=".",alpha=0.6)
+            scatter(X, Y, marker=".", alpha=0.6)
             title("Grid illustration")
          else
             error("unknown plot mode: $(plotmode[ivar])")
@@ -278,9 +277,8 @@ function plotdata(data::Data, func::AbstractString; cut="", plotmode="contbar",
          dim = [0.125, 0.013, 0.2, 0.045]
          str = @sprintf "it=%d, time=%4.2f" data.head.it data.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)
+            loc="lower left", prop=Dict("size"=>8), frameon=true,
+            bbox_to_anchor=(0., 1.), bbox_transform=ax.transAxes)
          at.patch.set_boxstyle("round,pad=0.,rounding_size=0.2")
          ax.add_artist(at)
          # recover status
@@ -292,8 +290,8 @@ function plotdata(data::Data, func::AbstractString; cut="", plotmode="contbar",
       Y = @view x[:,:,:,2]
       Z = @view x[:,:,:,3]
       for (ivar,var) in enumerate(vars)
-         if plotmode[ivar] ∈ ("surf","surfbar","surfbarlog","cont","contbar",
-            "contlog","contbarlog")
+         if plotmode[ivar] ∈ ("surf","surfbar","surfbarlog","cont","contbar", "contlog",
+               "contbarlog")
 
             VarIndex_ = findindex(data, var)
 
@@ -345,8 +343,8 @@ function plotdata(data::Data, func::AbstractString; cut="", plotmode="contbar",
             cut1, cut2, v1, v2 = subsurface(cut1, cut2, v1, v2, plotrange)
          end
 
-         if plotmode[ivar] ∈ ("surf","surfbar","surfbarlog","cont","contbar",
-            "contlog","contbarlog")
+         if plotmode[ivar] ∈ ("surf", "surfbar", "surfbarlog", "cont", "contbar", "contlog",
+               "contbarlog")
             if level == 0
                c = ax.contourf(cut1, cut2, W)
             else
@@ -355,19 +353,15 @@ function plotdata(data::Data, func::AbstractString; cut="", plotmode="contbar",
             fig.colorbar(c, ax=ax)
             title(data.head.wnames[VarIndex_])
 
-         elseif plotmode[ivar] ∈ ("stream","streamover")
-            # Surprisingly, some box outputs do not have equal spaces???
-            #xi = range(cut1[1,1], stop=cut1[1,end], length=size(cut1)[2])
-            #yi = range(cut2[1,1], stop=cut2[end,1], length=size(cut2)[1])
-
+         elseif plotmode[ivar] ∈ ("stream", "streamover")
             xi = range(cut1[1,1], stop=cut1[1,end],
-            step=(cut1[1,end] - cut1[1,1]) / (size(cut1,2) - 1))
+               step=(cut1[1,end] - cut1[1,1]) / (size(cut1,2) - 1))
             yi = range(cut2[1,1], stop=cut2[end,1],
-            step=(cut2[end,1] - cut2[1,1]) / (size(cut2,1) - 1))
+               step=(cut2[end,1] - cut2[1,1]) / (size(cut2,1) - 1))
 
             Xi, Yi = meshgrid(xi, yi)
 
-            s = streamplot(Xi,Yi,v1,v2; color="w", linewidth=1.0, density)
+            s = streamplot(Xi, Yi, v1, v2; color="w", linewidth=1.0, density)
          end
 
          if cut == "x"
@@ -382,9 +376,8 @@ function plotdata(data::Data, func::AbstractString; cut="", plotmode="contbar",
          dim = [0.125, 0.013, 0.2, 0.045]
          str = @sprintf "it=%d, time=%4.2f" data.head.it data.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)
+            loc="lower left", prop=Dict("size"=>8), frameon=true,
+            bbox_to_anchor=(0., 1.), bbox_transform=ax.transAxes)
          at.patch.set_boxstyle("round,pad=0.,rounding_size=0.2")
          ax.add_artist(at)
       end
@@ -394,13 +387,12 @@ end
 
 
 """
-    cutplot(data, var; plotrange=[-Inf,Inf,-Inf,Inf], cut=' ',
-       plotinterval=0.1, density=1.0, cutPlaneIndex=1,level=20)
+    cutplot(data, var; plotrange=[-Inf,Inf,-Inf,Inf], cut="x", plotinterval=0.1,
+       density=1.0, cutPlaneIndex=1, level=20)
 
 2D plane cut contourf of 3D box data.
 """
-function cutplot(data::Data, var::AbstractString;
-   plotrange=[-Inf,Inf,-Inf,Inf], cut=' ',
+function cutplot(data::Data, var::AbstractString; plotrange=[-Inf,Inf,-Inf,Inf], cut="x",
    plotinterval=0.1, density=1.0, cutPlaneIndex=1,level=20)
 
    x, w = data.x, data.w
@@ -412,15 +404,15 @@ function cutplot(data::Data, var::AbstractString;
 
    W = w[:,:,:,VarIndex_]
 
-   if cut ∈ ('x',' ')
+   if cut == "x"
       cut1 = @view X[cutPlaneIndex,:,:]
       cut2 = @view Y[cutPlaneIndex,:,:]
       W    = @view W[cutPlaneIndex,:,:]
-   elseif cut ==  'y'
+   elseif cut == "y"
       cut1 = @view X[:,cutPlaneIndex,:]
       cut2 = @view Z[:,cutPlaneIndex,:]
       W    = @view W[:,cutPlaneIndex,:]
-   elseif cut == 'z'
+   elseif cut == "z"
       cut1 = @view X[:,:,cutPlaneIndex]
       cut2 = @view Y[:,:,cutPlaneIndex]
       W    = @view W[:,:,cutPlaneIndex]
@@ -442,24 +434,22 @@ function cutplot(data::Data, var::AbstractString;
       xlabel("x"); ylabel("y")
    end
 
-   return c::PyCall.PyObject
+   c
 end
 
 
 """
-    streamslice(data::Data, var::String;
-       plotrange=[-Inf,Inf,-Inf,Inf], cut=' ',
+    streamslice(data::Data, var::String; plotrange=[-Inf,Inf,-Inf,Inf], cut="x",
        plotinterval=0.1, density=1.0, cutPlaneIndex=1, color="w", linewidth=1.0)
 
-Plot streamlines on 2D slices of 3D box data. Variable string must be separated
-with `;`. Tranposes aree needed because of `meshgrid` and `ndgrid` conversion.
+Plot streamlines on 2D slices of 3D box data. Variable names in `var` string must be
+separated with `;`.
 """
 function streamslice(data::Data, var::AbstractString;
-   plotrange=[-Inf,Inf,-Inf,Inf], cut=' ', cutPlaneIndex=1,
-   plotinterval=0.1, kwargs...)
+   plotrange=[-Inf,Inf,-Inf,Inf], cut="x", cutPlaneIndex=1, plotinterval=0.1, kwargs...)
 
    x,w = data.x, data.w
-   varStream  = split(var,";")
+   varStream  = split(var, ";")
    VarIndex1_ = findindex(data, varStream[1])
    VarIndex2_ = findindex(data, varStream[2])
 
@@ -470,17 +460,17 @@ function streamslice(data::Data, var::AbstractString;
    v1 = @view w[:,:,:,VarIndex1_]
    v2 = @view w[:,:,:,VarIndex2_]
 
-   if cut ∈ ('x',' ')
+   if cut == "x"
       cut1 = @view X[cutPlaneIndex,:,:]
       cut2 = @view Y[cutPlaneIndex,:,:]
       v1   = v1[cutPlaneIndex,:,:]
       v2   = v2[cutPlaneIndex,:,:]
-   elseif cut ==  'y'
+   elseif cut ==  "y"
       cut1 = @view X[:,cutPlaneIndex,:]
       cut2 = @view Z[:,cutPlaneIndex,:]
       v1   = v1[:,cutPlaneIndex,:]
       v2   = v2[:,cutPlaneIndex,:]
-   elseif cut == 'z'
+   elseif cut == "z"
       cut1 = @view X[:,:,cutPlaneIndex]
       cut2 = @view Y[:,:,cutPlaneIndex]
       v1   = v1[:,:,cutPlaneIndex]
@@ -492,13 +482,13 @@ function streamslice(data::Data, var::AbstractString;
    end
 
    xi = range(cut1[1,1], stop=cut1[end,1],
-   step = (cut1[end,1] - cut1[1,1]) / (size(cut1,1) - 1))
+      step = (cut1[end,1] - cut1[1,1]) / (size(cut1,1) - 1))
    yi = range(cut2[1,1], stop=cut2[1,end],
-   step = (cut2[1,end] - cut2[1,1]) / (size(cut2,2) - 1))
+      step = (cut2[1,end] - cut2[1,1]) / (size(cut2,2) - 1))
 
    Xi, Yi = meshgrid(xi, yi)
 
-   s = streamplot(Xi,Yi,v1',v2'; kwargs...)
+   s = streamplot(Xi, Yi, v1', v2'; kwargs...)
 
    if cut == 'x'
       xlabel("y"); ylabel("z")
@@ -521,8 +511,6 @@ function plot(data::Data, var::AbstractString; kwargs...)
    VarIndex_ = findindex(data, var)
 
    c = plot(x, w[:,VarIndex_]; kwargs...)
-
-   return c::Vector{PyCall.PyObject}
 end
 
 """
@@ -535,8 +523,6 @@ function scatter(data::Data, var::AbstractString; kwargs...)
    VarIndex_ = findindex(data, var)
 
    c = scatter(x, w[:,VarIndex_]; kwargs...)
-
-   return c::Vector{PyCall.PyObject}
 end
 
 """
@@ -602,7 +588,7 @@ end
 
 """
     plot_trisurf(data::Data, var::String;
-       plotrange::Vector{Float64}=[-Inf,Inf,-Inf,Inf], kwargs::Dict=Dict())
+       plotrange::Vector{Float64}=[-Inf,Inf,-Inf,Inf], kwargs...)
 
 Wrapper over the plot_trisurf function in matplotlib.
 """
@@ -635,7 +621,7 @@ end
 Wrapper over the plot_surface function in matplotlib.
 """
 function plot_surface(data::Data, var::AbstractString;
-   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, kwargs=Dict())
+   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, kwargs...)
 
    Xi, Yi, Wi = getdata(data, var, plotrange, plotinterval, 2)
 
@@ -658,7 +644,7 @@ function streamplot(data::Data, var::AbstractString;
    VarIndex1_ = findfirst(x->x==lowercase(VarStream[1]), wnames)
    VarIndex2_ = findfirst(x->x==lowercase(VarStream[2]), wnames)
 
-   if data.head.gencoord # Generalized coordinates
+   if data.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
diff --git a/src/plot/utility.jl b/src/plot/utility.jl
index 372a0bf1..a41daea7 100644
--- a/src/plot/utility.jl
+++ b/src/plot/utility.jl
@@ -1,8 +1,7 @@
 # Utility functions for plotting.
 
 "Prepare 2D data arrays for passing to plotting functions."
-function getdata(data::Data, var::AbstractString, plotrange, plotinterval,
-   griddim=1)
+function getdata(data::Data, var::AbstractString, plotrange, plotinterval, griddim=1)
    @assert data.head.ndim == 2 "data must be in 2D!"
 
    x, w = data.x, data.w
diff --git a/src/runReadData.jl b/src/runReadData.jl
deleted file mode 100644
index ce0151a8..00000000
--- a/src/runReadData.jl
+++ /dev/null
@@ -1,11 +0,0 @@
-# Script for using Batsrus loader in REPL
-#
-# Hongyang Zhou, hyzhou@umich.edu 08/02/2019
-
-# Check if filename
-
-# Check if present npict
-
-# If not, read in npict from command line
-
-filehead, data, filelist = readdata(filename);
diff --git a/src/unit/UnitfulBatsrus.jl b/src/unit/UnitfulBatsrus.jl
index 9c8e16fd..857d4835 100644
--- a/src/unit/UnitfulBatsrus.jl
+++ b/src/unit/UnitfulBatsrus.jl
@@ -8,7 +8,9 @@ export getunit, getunits
 
 
 # lengths
-@unit R_bu  "Rₑ"     Radii             (6378)*Unitful.km        false
+@unit R_bu  "Rₑ"     EarthRadii        (6378)*Unitful.km        false
+@unit Rg_bu  "Rg"    GanymedeRadii     (2634)*Unitful.km        false # Ganymede's radius
+@unit Rm_bu  "Rm"    MercuryRadii      (2444)*Unitful.km        false # Mercury's radius 
 
 # masses
 @unit me_bu "mₑ"     ElectronMass      me                       false
diff --git a/src/unit/batsrusmacro.jl b/src/unit/batsrusmacro.jl
index 5053d8b2..b718e68a 100644
--- a/src/unit/batsrusmacro.jl
+++ b/src/unit/batsrusmacro.jl
@@ -1,9 +1,9 @@
 """
     macro bu_str(unit)
 
-String macro to easily recall Batsrus units located in the `UnitfulBatsrus`
-package. Although all unit symbols in that package are suffixed with `_bu`,
-the suffix should not be used when using this macro.
+String macro to easily recall Batsrus units located in the `UnitfulBatsrus` package.
+Although all unit symbols in that package are suffixed with `_bu`, the suffix should not be
+used when using this macro.
 Note that what goes inside must be parsable as a valid Julia expression.
 Examples:
 ```
@@ -50,7 +50,7 @@ dottify(s) = s
 function replace_value(sym::Symbol)
    s = Symbol(sym, :_bu)
    if !(isdefined(UnitfulBatsrus, s) && typecheck_bool(getfield(UnitfulBatsrus, s)))
-       error("Symbol $s could not be found in UnitfulBatsrus.")
+      error("Symbol $s could not be found in UnitfulBatsrus.")
    end
 
    return getfield(UnitfulBatsrus, s)
diff --git a/test/runtests.jl b/test/runtests.jl
index 9e1f888c..0683d897 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,11 +1,12 @@
 # Tests of BATSRUS.jl
 
 using Batsrus, Test, SHA
+using Batsrus.UnitfulBatsrus, Unitful
 
 function filecmp(path1::AbstractString, path2::AbstractString)
    stat1, stat2 = stat(path1), stat(path2)
    if !(isfile(stat1) && isfile(stat2)) || filesize(stat1) != filesize(stat2)
-      return false # or should it throw if a file doesn't exist?
+      return false
    end
    stat1 == stat2 && return true # same file
    open(path1, "r") do file1
@@ -107,4 +108,23 @@ end
          @test isa(ax, PyPlot.PyObject)
       end
    end
+
+   @testset "Units" begin
+      @test 1.0bu"R" > 2.0bu"Rg"
+      filename = "y=0_var_1_t00000000_n00000000.out"
+      data = readdata(filename, dir="data")
+      varunit = getunit(data, "Rho")
+      @test varunit == bu"amucc"
+      varunit = getunit(data, "Ux")
+      @test varunit == u"km/s"
+      varunit = getunit(data, "Bx")
+      @test varunit == u"nT"
+      varunit = getunit(data, "P")
+      @test varunit == u"nPa"
+      varunit = getunit(data, "jx")
+      @test dimension(varunit) == dimension(Unitful.A/Unitful.m^2)
+      varunit = getunit(data, "ex")
+      @test varunit == u"mV/m"
+   end
+
 end
\ No newline at end of file