From e4898a5e05ff07a4970c121834d36daa34dd2c44 Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Mon, 15 Jul 2024 16:59:53 +0800
Subject: [PATCH] Add point & line interpolations; Refactor (#57)

* Bake dimensionality into BATLData

* Bump minor version
---
 Project.toml                           |   2 +-
 docs/src/man/manual.md                 |  21 ++
 ext/BatsrusMakieExt/BatsrusMakieExt.jl |   2 +-
 ext/BatsrusMakieExt/typerecipe.jl      |   2 +-
 src/Batsrus.jl                         |   7 +-
 src/io.jl                              |  82 +++--
 src/plot/plots.jl                      |  49 +--
 src/plot/pyplot.jl                     | 486 +++++++++++++------------
 src/plot/utility.jl                    |  89 ++++-
 src/select.jl                          |  23 +-
 test/runtests.jl                       |  23 +-
 11 files changed, 457 insertions(+), 329 deletions(-)

diff --git a/Project.toml b/Project.toml
index 82c9d396..2186a5b8 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.5.11"
+version = "0.6.0"
 
 [deps]
 FortranFiles = "c58ffaec-ab22-586d-bfc5-781a99fd0b10"
diff --git a/docs/src/man/manual.md b/docs/src/man/manual.md
index d1c8edfc..065a54fa 100644
--- a/docs/src/man/manual.md
+++ b/docs/src/man/manual.md
@@ -27,6 +27,12 @@ head, data = readlogdata(logfilename)
 
 ### Data Extraction
 
+- Checking variable range
+
+```julia
+get_var_range(bd, "rho")
+```
+
 - Raw variables
 
 ```julia
@@ -41,6 +47,21 @@ v = getvars(bd, ["Bx", "By", "Bz"])
 Bmag = bd["Bmag"]
 ```
 
+- Extracting data at a given location
+
+```julia
+loc = Float32[0.0, 0.0] # The type determines the output type
+d = interp1d(bd, "rho", loc)
+```
+
+- Extracting data along a given line
+
+```julia
+point1 = Float32[-10.0, -1.0]
+point2 = Float32[10.0, 1.0]
+w = slice1d(bd, "rho", point1, point2)
+```
+
 Here is a full list of predefined derived quantities:
 
 | Derived variable name | Meaning                          | Required variable |
diff --git a/ext/BatsrusMakieExt/BatsrusMakieExt.jl b/ext/BatsrusMakieExt/BatsrusMakieExt.jl
index e3e93914..851a1e59 100644
--- a/ext/BatsrusMakieExt/BatsrusMakieExt.jl
+++ b/ext/BatsrusMakieExt/BatsrusMakieExt.jl
@@ -1,7 +1,7 @@
 module BatsrusMakieExt
 
 using Batsrus, Printf
-import Batsrus: findindex, hasunit, getunit, getdata2d
+import Batsrus: findindex, hasunit, getunit, slice2d
 import Batsrus.UnitfulBatsrus
 import Makie
 
diff --git a/ext/BatsrusMakieExt/typerecipe.jl b/ext/BatsrusMakieExt/typerecipe.jl
index fa1c78eb..a795d512 100644
--- a/ext/BatsrusMakieExt/typerecipe.jl
+++ b/ext/BatsrusMakieExt/typerecipe.jl
@@ -19,7 +19,7 @@ end
 "Conversion for 2D plots."
 function Makie.convert_arguments(P::Makie.GridBased, bd::BATLData, var::String;
    plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1)
-   x, y, w = getdata2d(bd, var, plotrange, plotinterval)
+   x, y, w = slice2d(bd, var, plotrange, plotinterval)
 
    unitx = getunit(bd, bd.head.variables[1])
    unity = getunit(bd, bd.head.variables[2])
diff --git a/src/Batsrus.jl b/src/Batsrus.jl
index 039a56a8..eb802354 100644
--- a/src/Batsrus.jl
+++ b/src/Batsrus.jl
@@ -5,13 +5,14 @@ module Batsrus
 
 using Printf, Reexport, Requires
 using Parsers
-using NaturalNeighbours: interpolate, Triangle
+using Interpolations: cubic_spline_interpolation, BSpline, Linear, scale, interpolate
+import NaturalNeighbours as NN
 
 export BATLData,
    load, readlogdata, readtecdata, showhead, # io
    getvars, getvar, cutdata, subvolume, subsurface, # select
    Batl, convertTECtoVTU, convertIDLtoVTK, readhead, readtree, getConnectivity, # vtk
-   getdata2d # plot/utility
+   interp1d, slice1d, slice2d, get_var_range # plot/utility
 
 "Type for the file information."
 struct FileList
@@ -30,7 +31,7 @@ struct FileList
 end
 
 "Primary Batsrus data storage type."
-struct BATLData{T<:AbstractFloat}
+struct BATLData{dim, T<:AbstractFloat}
    "header information"
    head::NamedTuple
    "grid"
diff --git a/src/io.jl b/src/io.jl
index 5f185082..29926444 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -21,31 +21,26 @@ function load(file::AbstractString; npict::Int=1, verbose::Bool=false)
    seekstart(fileID) # Rewind to start
 
    ## Read data from files
-   # Skip npict-1 snapshots (because we only want the npict-th snapshot)
+   # Skip npict-1 snapshots (since we only want the npict-th snapshot)
    skip(fileID, pictsize*(npict-1))
 
    filehead = getfilehead(fileID, filelist)
-
    # Read data
    if filelist.type == :ascii
       x, w = allocateBuffer(filehead, Float64) # why Float64?
-      getascii!(x, w, fileID, filehead)
+      getascii!(x, w, fileID)
    else
       skip(fileID, TAG) # skip record start tag
       T = filelist.type == :real4 ? Float32 : Float64
       x, w = allocateBuffer(filehead, T)
-      getbinary!(x, w, fileID, filehead)
+      getbinary!(x, w, fileID)
    end
 
    close(fileID)
 
    #setunits(filehead,"")
 
-	data = BATLData(filehead, x, w, filelist)
-
-   verbose && @info "Finished reading $(filelist.name)"
-
-   data
+	data = BATLData{Int(filehead.ndim), eltype(w)}(filehead, x, w, filelist)
 end
 
 "Read information from log file."
@@ -172,22 +167,22 @@ function readtecdata(file::AbstractString; verbose::Bool=false)
    data = Array{Float32,2}(undef, length(VARS), nNode)
 
    if nDim == 3
-	   connectivity = Array{Int32,2}(undef,8,nCell)
+	   connectivity = Array{Int32,2}(undef, 8, nCell)
    elseif nDim == 2
-	   connectivity = Array{Int32,2}(undef,4,nCell)
+	   connectivity = Array{Int32,2}(undef, 4, nCell)
    end
 
-   IsBinary = false
+   isBinary = false
    try
       Parsers.parse.(Float32, split(readline(f)))
    catch
-      IsBinary = true
+      isBinary = true
       verbose && @info "reading binary file"
    end
 
    seek(f, pt0)
 
-   if IsBinary
+   if isBinary
 	   @inbounds for i in 1:nNode
 		   read!(f, @view data[:,i])
 	   end
@@ -271,8 +266,7 @@ function getfilehead(fileID::IOStream, filelist::FileList)
 
    if type == :ascii
       headline = readline(fileID)
-      line = readline(fileID)
-      line = split(line)
+      line = readline(fileID) |> split
       it = Parsers.parse(Int, line[1])
       t = Parsers.parse(Float64, line[2])
       ndim = Parsers.parse(Int8, line[3])
@@ -326,8 +320,8 @@ end
 
 function skipline(s::IO)
    while !eof(s)
-       c = read(s, Char)
-       c == '\n' && break
+      c = read(s, Char)
+      c == '\n' && break
    end
 
    return
@@ -367,11 +361,9 @@ function getfilesize(fileID::IOStream, type::Symbol, lenstr::Int32)
       read(fileID, lenstr)
       skip(fileID, TAG)
    end
-
    # Header length
    pointer1 = position(fileID)
    headlen = pointer1 - pointer0
-
    # Calculate the snapshot size = header + data + recordmarks
    nxs = prod(nx)
    pictsize =
@@ -408,29 +400,59 @@ function allocateBuffer(filehead::NamedTuple, T::DataType)
 end
 
 "Read ascii format coordinates and data values."
-function getascii!(x, w, fileID::IOStream, filehead::NamedTuple)
-   ndim = filehead.ndim
-   Ids = CartesianIndices(size(x)[1:ndim])
-   @inbounds @views for ids in Ids
+function getascii!(x::Array{T, 2}, w, fileID::IOStream) where T
+   @inbounds @views for id in axes(x, 1)
+      temp = Parsers.parse.(Float64, split(readline(fileID)))
+      x[id] = temp[1]
+      w[id,:] = temp[2:end]
+   end
+end
+
+function getascii!(x::Array{T, 3}, w, fileID::IOStream) where T
+   @inbounds @views for id in CartesianIndices(size(x)[1:2])
       temp = Parsers.parse.(Float64, split(readline(fileID)))
-      x[ids,:] = temp[1:ndim]
-      w[ids,:] = temp[ndim+1:end]
+      x[id,:] = temp[1:2]
+      w[id,:] = temp[3:end]
    end
+end
 
-   return
+function getascii!(x::Array{T, 4}, w, fileID::IOStream) where T
+   @inbounds @views for id in CartesianIndices(size(x)[1:3])
+      temp = Parsers.parse.(Float64, split(readline(fileID)))
+      x[id,:] = temp[1:3]
+      w[id,:] = temp[4:end]
+   end
 end
 
 "Read binary format coordinates and data values."
-function getbinary!(x, w, fileID::IOStream, filehead::NamedTuple)
-   dimlast = filehead.ndim + 1
+function getbinary!(x::Array{T, 2}, w, fileID::IOStream) where T
    read!(fileID, x)
    skip(fileID, 2*TAG)
+   dimlast = 2
    @inbounds for iw in axes(w, dimlast)
       read!(fileID, selectdim(w, dimlast, iw))
       skip(fileID, 2*TAG)
    end
+end
 
-   return
+function getbinary!(x::Array{T, 3}, w, fileID::IOStream) where T
+   read!(fileID, x)
+   skip(fileID, 2*TAG)
+   dimlast = 3
+   @inbounds for iw in axes(w, dimlast)
+      read!(fileID, selectdim(w, dimlast, iw))
+      skip(fileID, 2*TAG)
+   end
+end
+
+function getbinary!(x::Array{T, 4}, w, fileID::IOStream) where T
+   read!(fileID, x)
+   skip(fileID, 2*TAG)
+   dimlast = 4
+   @inbounds for iw in axes(w, dimlast)
+      read!(fileID, selectdim(w, dimlast, iw))
+      skip(fileID, 2*TAG)
+   end
 end
 
 """
diff --git a/src/plot/plots.jl b/src/plot/plots.jl
index b2e07e2d..550ecdc5 100644
--- a/src/plot/plots.jl
+++ b/src/plot/plots.jl
@@ -3,30 +3,31 @@
 using RecipesBase
 
 # Build a recipe which acts on a custom type.
-@recipe function f(bd::BATLData, var::AbstractString;
-   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1)
-   ndim = bd.head.ndim
+@recipe function f(bd::BATLData{1, T}, var::AbstractString) where {T}
    hasunits = hasunit(bd)
 
-   if ndim == 1
-      VarIndex_ = findindex(bd, var)
-      if hasunits
-         unitx = getunit(bd, bd.head.variables[1])
-         unitw = getunit(bd, var)
-         x = bd.x .* unitx
-         y = bd.w[:,VarIndex_] .* unitw
-      else
-         x = bd.x
-         y = @view bd.w[:,VarIndex_]
-      end
+   if hasunits
+      unitx = getunit(bd, bd.head.variables[1])
+      unitw = getunit(bd, var)
+      x = bd.x .* unitx
+      y = getview(bd, var) .* unitw
+   else
+      x = bd.x
+      y = getview(bd, var)
+   end
 
-      @series begin
-         seriestype --> :path
-         x, y
-      end
-   elseif ndim == 2
-      x, y, w = getdata2d(bd, var, plotrange, plotinterval)
+   @series begin
+      seriestype --> :path
+      x, y
+   end
+end
+
+@recipe function f(bd::BATLData{2, T}, var::AbstractString;
+   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1) where {T}
+   hasunits = hasunit(bd)
 
+   x, y, w = slice2d(bd, var, plotrange, plotinterval)
+   if hasunits
       unitx = getunit(bd, bd.head.variables[1])
       unity = getunit(bd, bd.head.variables[2])
       unitw = getunit(bd, var)
@@ -36,10 +37,10 @@ using RecipesBase
          y *= unity
          w *= unitw
       end
+   end
 
-      @series begin
-         seriestype --> :contourf  # use := if you want to force it
-         x, y, w'
-      end
+   @series begin
+      seriestype --> :contourf  # use := if you want to force it
+      x, y, w'
    end
 end
\ No newline at end of file
diff --git a/src/plot/pyplot.jl b/src/plot/pyplot.jl
index 6a02b2e9..cc14177e 100644
--- a/src/plot/pyplot.jl
+++ b/src/plot/pyplot.jl
@@ -1,7 +1,6 @@
 # Plotting functionalities.
 
 using PyPlot
-using Interpolations: cubic_spline_interpolation
 
 export plotdata, plotlogdata, plot, scatter, contour, contourf, plot_surface, tripcolor,
    tricontourf, plot_trisurf, streamplot, streamslice, quiver, cutplot, pcolormesh
@@ -55,7 +54,7 @@ end
 
 
 """
-    plotdata(bd, func, args, kwargs...)
+    plotdata(bd::BATLData, func, args, kwargs...)
 
 Plot the variable from SWMF output.
 
@@ -78,252 +77,279 @@ Plot the variable from SWMF output.
 - `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.
+- `multifigure::Bool`: whether to use multifigure display or 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.
 """
-function plotdata(bd::BATLData, func::AbstractString; dir="x", plotmode="contbar",
-   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, sequence=1, multifigure=true,
-   getrangeonly::Bool=false, levels::Int=0, innermask::Bool=false, verbose::Bool=false,
-   density=1.0, stride::Int=10, kwargs...)
-
+function plotdata(bd::BATLData{1, T}, func::AbstractString; verbose::Bool=false,
+   plotmode="contbar", multifigure::Bool=true, kwargs...) where T
    x, w = bd.x, bd.w
    plotmode = split(plotmode)
    vars     = split(func)
-   ndim     = bd.head.ndim
    nvar     = length(vars)
 
-   if verbose || getrangeonly
-      @info "============ PLOTTING PARAMETERS ==============="
-      @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(bd, var)
-         if ndim == 1
-            wmin[ivar], wmax[ivar] = extrema(@view w[:,varIndex_])
-         elseif ndim == 2
-            wmin[ivar], wmax[ivar] = extrema(@view w[:,:,varIndex_])
-         end
-         println("Min & Max value for $(var) :$(wmin[ivar])",", $(wmax[ivar])")
+   verbose && show_ranges(bd, vars)
+
+   ## plot multiple variables with same plotmode
+   if length(plotmode) < nvar
+      [push!(plotmode, plotmode[i]) for i in 1:nvar-length(plotmode)]
+   end
+
+   ## Plot
+   for (ivar, var) in enumerate(vars)
+      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...)
+      else
+         scatter(x, w[:,varIndex_]; kwargs...)
+      end
+      if occursin("grid", plotmode[ivar])
+         grid(true)
       end
-      if getrangeonly return wmin, wmax end
+      xlabel("x"); ylabel("$(var)")
+      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)
+      at.patch.set_boxstyle("round,pad=0.,rounding_size=0.2")
+      ax.add_artist(at)
    end
 
+   return
+end
+
+
+function plotdata(bd::BATLData{2, T}, func::AbstractString; verbose::Bool=false,
+   plotmode="contbar", plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, density=1.0,
+   multifigure::Bool=true, levels::Int=0, innermask::Bool=false, stride::Int=10, kwargs...) where T
+   x, w = bd.x, bd.w
+   plotmode = split(plotmode)
+   vars     = split(func)
+   nvar     = length(vars)
+
+   verbose && show_ranges(bd, vars)
+
    ## plot multiple variables with same plotmode
    if length(plotmode) < nvar
-      [push!(plotmode, plotmode[i]) for i = 1:nvar-length(plotmode)]
+      [push!(plotmode, plotmode[i]) for i in 1:nvar-length(plotmode)]
    end
 
-   ## Plot
-   if ndim == 1
-      for (ivar,var) in enumerate(vars)
+   for (ivar, var) in enumerate(vars)
+      occursin("over", plotmode[ivar]) && (multifigure = false)
+      if ivar == 1 || multifigure fig, ax = subplots() else ax = gca() end
+      if !occursin(";", 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...)
-         else
-            scatter(x, w[:,varIndex_]; kwargs...)
-         end
-         if occursin("grid", plotmode[ivar])
-            grid(true)
-         end
-         xlabel("x"); ylabel("$(var)")
-         dim = [0.125, 0.013, 0.2, 0.045]
-         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)
-         at.patch.set_boxstyle("round,pad=0.,rounding_size=0.2")
-         ax.add_artist(at)
       end
-   elseif ndim == 2
-      for (ivar,var) in enumerate(vars)
-         occursin("over", plotmode[ivar]) && (multifigure = false)
-         if ivar == 1 || multifigure fig, ax = subplots() else ax = gca() end
-         if !occursin(";",var)
-            varIndex_ = findindex(bd, var)
-         end
 
-         if plotmode[ivar] ∈ ("surf","surfbar","surfbarlog","cont","contbar",
-            "contlog","contbarlog")
-
-            if plotmode[ivar] ∈ ("contbar", "contbarlog")
-               c = contourf(bd, var; levels, plotrange, plotinterval, innermask, kwargs...)
-            elseif plotmode[ivar] ∈ ("cont", "contlog")
-               c = contour(bd, var; levels, plotrange, plotinterval, innermask, kwargs...)
-            elseif plotmode[ivar] ∈ ("surfbar", "surfbarlog")
-               c = plot_surface(bd, var; plotrange, plotinterval, innermask, kwargs...)
-            end
-
-            occursin("bar", plotmode[ivar]) && colorbar(c; fraction=0.04, pad=0.02)
-            occursin("log", plotmode[ivar]) && (c.locator = matplotlib.ticker.LogLocator())
-            title(bd.head.wnames[varIndex_])
-
-         elseif plotmode[ivar] ∈ ("trimesh","trisurf","tricont","tristream")
-            X = vec(x[:,:,1])
-            Y = vec(x[:,:,2])
-            W = vec(w[:,:,varIndex_])
-
-            #TODO This needs to be modified!!!
-            if !all(isinf.(plotrange))
-               xyIndex = X .> plotrange[1] .& X .< plotrange[2] .&
-                  Y .> plotrange[3] .& Y .< plotrange[4]
-               X = X[xyIndex]
-               Y = Y[xyIndex]
-               W = W[xyIndex]
-            end
-
-            if plotmode[ivar] == "trimesh"
-               triang = matplotlib.tri.Triangulation(X, Y)
-               c = ax.triplot(triang, kwargs...)
-            elseif plotmode[ivar] == "trisurf"
-               c = ax.plot_trisurf(X, Y, W'; kwargs...)
-            elseif plotmode[ivar] == "tricont"
-               c = ax.tricontourf(X, Y, W'; levels, kwargs...)
-               fig.colorbar(c; ax)
-            elseif plotmode[ivar] == "tristream"
-               throw(ArgumentError("tristream not yet implemented!"))
-            end
-
-            title(bd.head.wnames[varIndex_])
-
-         elseif plotmode[ivar] ∈ ("stream","streamover")
-            s = streamplot(bd, var; plotrange, plotinterval, color="w", linewidth=1.0,
-               density, kwargs...)
-         elseif occursin("quiver", plotmode[ivar])
-            q = quiver(bd, var; stride, color="w", kwargs...)
-         elseif occursin("grid", plotmode[ivar])
-            # This does not take subdomain plot into account!
-            X, Y = eachslice(x, dims=3)
-            scatter(X, Y, marker=".", alpha=0.6)
-            title("Grid illustration")
-         else
-            throw(ArgumentError("unknown plot mode: $(plotmode[ivar])"))
+      if plotmode[ivar] ∈ ("surf", "surfbar", "surfbarlog", "cont", "contbar",
+         "contlog", "contbarlog")
+
+         if plotmode[ivar] ∈ ("contbar", "contbarlog")
+            c = contourf(bd, var; levels, plotrange, plotinterval, innermask, kwargs...)
+         elseif plotmode[ivar] ∈ ("cont", "contlog")
+            c = contour(bd, var; levels, plotrange, plotinterval, innermask, kwargs...)
+         elseif plotmode[ivar] ∈ ("surfbar", "surfbarlog")
+            c = plot_surface(bd, var; plotrange, plotinterval, innermask, kwargs...)
          end
 
-         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" 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)
-         at.patch.set_boxstyle("round,pad=0.,rounding_size=0.2")
-         ax.add_artist(at)
-         # recover status
-         occursin("over", plotmode[ivar]) && (multifigure = true)
-      end
+         occursin("bar", plotmode[ivar]) && colorbar(c; fraction=0.04, pad=0.02)
+         occursin("log", plotmode[ivar]) && (c.locator = matplotlib.ticker.LogLocator())
+         title(bd.head.wnames[varIndex_])
 
-   else # 2D cut from 3D output; now only for Cartesian output
-      X, Y, Z = eachslice(x, dims=4)
-      for (ivar,var) in enumerate(vars)
-         if plotmode[ivar] ∈ ("surf","surfbar","surfbarlog","cont","contbar", "contlog",
-               "contbarlog")
-
-            varIndex_ = findindex(bd, var)
-
-            if ivar == 1 || multifigure fig, ax = subplots() else ax = gca() end
-
-            W = w[:,:,:,varIndex_]
-
-            if dir == "x"
-               cut1 = @view X[sequence,:,:]
-               cut2 = @view Y[sequence,:,:]
-               W    = @view W[sequence,:,:]
-            elseif dir == "y"
-               cut1 = @view X[:,sequence,:]
-               cut2 = @view Z[:,sequence,:]
-               W    = @view W[:,sequence,:]
-            elseif dir == "z"
-               cut1 = @view X[:,:,sequence]
-               cut2 = @view Y[:,:,sequence]
-               W    = @view W[:,:,sequence]
-            end
-         elseif plotmode[ivar] ∈ ("stream","streamover")
-            varstream  = split(var,";")
-            var1_ = findindex(bd, varstream[1])
-            var2_ = findindex(bd, varstream[2])
-
-            v1 = @view w[:,:,:,var1_]
-            v2 = @view w[:,:,:,var2_]
-
-            if dir == "x"
-               cut1 = @view Y[sequence,:,:]
-               cut2 = @view Z[sequence,:,:]
-               v1   = v1[sequence,:,:]'
-               v2   = v2[sequence,:,:]'
-            elseif dir == "y"
-               cut1 = @view X[:,sequence,:]
-               cut2 = @view Z[:,sequence,:]
-               v1   = v1[:,sequence,:]'
-               v2   = v2[:,sequence,:]'
-            elseif dir == "z"
-               cut1 = @view X[:,:,sequence]
-               cut2 = @view Y[:,:,sequence]
-               v1   = v1[:,:,sequence]'
-               v2   = v2[:,:,sequence]'
-            end
-            cut1, cut2 = cut1', cut2'
-         end
+      elseif plotmode[ivar] ∈ ("trimesh","trisurf","tricont","tristream")
+         X = vec(x[:,:,1])
+         Y = vec(x[:,:,2])
+         W = vec(w[:,:,varIndex_])
 
+         #TODO This needs to be modified!!!
          if !all(isinf.(plotrange))
-            cut1, cut2, v1, v2 = subsurface(cut1, cut2, v1, v2, plotrange)
+            xyIndex = X .> plotrange[1] .& X .< plotrange[2] .&
+               Y .> plotrange[3] .& Y .< plotrange[4]
+            X = X[xyIndex]
+            Y = Y[xyIndex]
+            W = W[xyIndex]
          end
 
-         if plotmode[ivar] ∈ ("surf", "surfbar", "surfbarlog", "cont", "contbar", "contlog",
-               "contbarlog")
-            if levels == 0
-               c = ax.contourf(cut1, cut2, W'; kwargs...)
-            else
-               c = ax.contourf(cut1, cut2, W', levels; kwargs...)
-            end
-            fig.colorbar(c, ax=ax)
-            title(bd.head.wnames[varIndex_])
-
-         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))
-            yi = range(cut2[1,1], stop=cut2[end,1],
-               step=(cut2[end,1] - cut2[1,1]) / (size(cut2,1) - 1))
-
-            s = streamplot(xi, yi, v1, v2; color="w", linewidth=1.0, density, kwargs...)
+         if plotmode[ivar] == "trimesh"
+            triang = matplotlib.tri.Triangulation(X, Y)
+            c = ax.triplot(triang, kwargs...)
+         elseif plotmode[ivar] == "trisurf"
+            c = ax.plot_trisurf(X, Y, W'; kwargs...)
+         elseif plotmode[ivar] == "tricont"
+            c = ax.tricontourf(X, Y, W'; levels, kwargs...)
+            fig.colorbar(c; ax)
+         elseif plotmode[ivar] == "tristream"
+            throw(ArgumentError("tristream not yet implemented!"))
          end
 
+         title(bd.head.wnames[varIndex_])
+
+      elseif plotmode[ivar] ∈ ("stream","streamover")
+         s = streamplot(bd, var; plotrange, plotinterval, color="w", linewidth=1.0,
+            density, kwargs...)
+      elseif occursin("quiver", plotmode[ivar])
+         q = quiver(bd, var; stride, color="w", kwargs...)
+      elseif occursin("grid", plotmode[ivar])
+         # This does not take subdomain plot into account!
+         X, Y = eachslice(x, dims=3)
+         scatter(X, Y, marker=".", alpha=0.6)
+         title("Grid illustration")
+      else
+         throw(ArgumentError("unknown plot mode: $(plotmode[ivar])"))
+      end
+
+      xlabel(bd.head.variables[1]); ylabel(bd.head.variables[2])
+      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)
+      at.patch.set_boxstyle("round,pad=0.,rounding_size=0.2")
+      ax.add_artist(at)
+      # recover status
+      occursin("over", plotmode[ivar]) && (multifigure = true)
+   end
+
+   return
+end
+
+function plotdata(bd::BATLData{3, T}, func::AbstractString; dir="x", plotmode="contbar",
+   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, sequence=1, multifigure=true,
+   levels::Int=0, innermask::Bool=false, verbose::Bool=false, density=1.0, kwargs...) where T
+   x, w = bd.x, bd.w
+   plotmode = split(plotmode)
+   vars     = split(func)
+   nvar     = length(vars)
+
+   verbose && show_ranges(bd, vars)
+
+   ## plot multiple variables with same plotmode
+   if length(plotmode) < nvar
+      [push!(plotmode, plotmode[i]) for i = 1:nvar-length(plotmode)]
+   end
+
+   # 2D cut from 3D output; now only for Cartesian output
+   X, Y, Z = eachslice(x, dims=4)
+   for (ivar, var) in enumerate(vars)
+      if plotmode[ivar] ∈ ("surf","surfbar","surfbarlog","cont","contbar", "contlog",
+            "contbarlog")
+         varIndex_ = findindex(bd, var)
+
+         if ivar == 1 || multifigure fig, ax = subplots() else ax = gca() end
+
+         W = w[:,:,:,varIndex_]
+
+         if dir == "x"
+            cut1 = @view X[sequence,:,:]
+            cut2 = @view Y[sequence,:,:]
+            W    = @view W[sequence,:,:]
+         elseif dir == "y"
+            cut1 = @view X[:,sequence,:]
+            cut2 = @view Z[:,sequence,:]
+            W    = @view W[:,sequence,:]
+         elseif dir == "z"
+            cut1 = @view X[:,:,sequence]
+            cut2 = @view Y[:,:,sequence]
+            W    = @view W[:,:,sequence]
+         end
+      elseif plotmode[ivar] ∈ ("stream","streamover")
+         varstream  = split(var,";")
+         var1_ = findindex(bd, varstream[1])
+         var2_ = findindex(bd, varstream[2])
+
+         v1 = @view w[:,:,:,var1_]
+         v2 = @view w[:,:,:,var2_]
+
          if dir == "x"
-            xlabel("y"); ylabel("z")
+            cut1 = @view Y[sequence,:,:]
+            cut2 = @view Z[sequence,:,:]
+            v1   = v1[sequence,:,:]'
+            v2   = v2[sequence,:,:]'
          elseif dir == "y"
-            xlabel("x"); ylabel("z")
+            cut1 = @view X[:,sequence,:]
+            cut2 = @view Z[:,sequence,:]
+            v1   = v1[:,sequence,:]'
+            v2   = v2[:,sequence,:]'
          elseif dir == "z"
-            xlabel("x"); ylabel("y")
+            cut1 = @view X[:,:,sequence]
+            cut2 = @view Y[:,:,sequence]
+            v1   = v1[:,:,sequence]'
+            v2   = v2[:,:,sequence]'
+         end
+         cut1, cut2 = cut1', cut2'
+      end
+
+      if !all(isinf.(plotrange))
+         cut1, cut2, v1, v2 = subsurface(cut1, cut2, v1, v2, plotrange)
+      end
+
+      if plotmode[ivar] ∈ ("surf", "surfbar", "surfbarlog", "cont", "contbar", "contlog",
+            "contbarlog")
+         if levels == 0
+            c = ax.contourf(cut1, cut2, W'; kwargs...)
+         else
+            c = ax.contourf(cut1, cut2, W', levels; kwargs...)
          end
+         fig.colorbar(c, ax=ax)
+         title(bd.head.wnames[varIndex_])
+
+      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))
+         yi = range(cut2[1,1], stop=cut2[end,1],
+            step=(cut2[end,1] - cut2[1,1]) / (size(cut2,1) - 1))
+         s = streamplot(xi, yi, v1, v2; color="w", linewidth=1.0, density, kwargs...)
+      end
 
-         ax = gca()
-         dim = [0.125, 0.013, 0.2, 0.045]
-         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)
-         at.patch.set_boxstyle("round,pad=0.,rounding_size=0.2")
-         ax.add_artist(at)
+      if dir == "x"
+         xlabel("y"); ylabel("z")
+      elseif dir == "y"
+         xlabel("x"); ylabel("z")
+      elseif dir == "z"
+         xlabel("x"); ylabel("y")
       end
+
+      ax = gca()
+      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)
+      at.patch.set_boxstyle("round,pad=0.,rounding_size=0.2")
+      ax.add_artist(at)
    end
 
    return
 end
 
+function show_ranges(bd::BATLData, vars)
+   @info "============ PLOTTING PARAMETERS ==============="
+   @info "wnames = $(bd.head.wnames)"
+   nvar = length(vars)
+   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(bd, var)
+      if ndim == 1
+         wmin[ivar], wmax[ivar] = extrema(@view w[:,varIndex_])
+      elseif ndim == 2
+         wmin[ivar], wmax[ivar] = extrema(@view w[:,:,varIndex_])
+      end
+      println("Min & Max value for $(var) :$(wmin[ivar])",", $(wmax[ivar])")
+   end
+end
 
 """
     cutplot(data, var, ax=nothing; plotrange=[-Inf,Inf,-Inf,Inf], dir="x", sequence=1)
 
 2D plane cut pcolormesh of 3D box data. `sequence` is the index along `dir`.
 """
-function cutplot(bd::BATLData, var::AbstractString, ax=nothing;
-   plotrange=[-Inf,Inf,-Inf,Inf], dir="x", sequence=1)
+function cutplot(bd::BATLData{3, T}, var::AbstractString, ax=nothing;
+   plotrange=[-Inf,Inf,-Inf,Inf], dir="x", sequence=1) where {T}
    x, w = bd.x, bd.w
    varIndex_ = findindex(bd, var)
 
@@ -372,8 +398,8 @@ end
 Plot streamlines on 2D slices of 3D box data. Variable names in `var` string must be
 separated with `;`.
 """
-function streamslice(bd::BATLData, var::AbstractString, ax=nothing;
-   plotrange=[-Inf,Inf,-Inf,Inf], dir="x", sequence=1, kwargs...)
+function streamslice(bd::BATLData{3, T}, var::AbstractString, ax=nothing;
+   plotrange=[-Inf,Inf,-Inf,Inf], dir="x", sequence=1, kwargs...) where {T}
    x,w = bd.x, bd.w
    varstream  = split(var, ";")
    var1_ = findindex(bd, varstream[1])
@@ -430,7 +456,7 @@ end
 
 Wrapper over `plot` in matplotlib.
 """
-function PyPlot.plot(bd::BATLData, var::AbstractString, ax=nothing; kwargs...)
+function PyPlot.plot(bd::BATLData{1, T}, var::AbstractString, ax=nothing; kwargs...) where T
    x, w = bd.x, bd.w
    varIndex_ = findindex(bd, var)
    if isnothing(ax) ax = plt.gca() end
@@ -443,7 +469,7 @@ end
 
 Wrapper over `scatter` in matplotlib.
 """
-function PyPlot.scatter(bd::BATLData, var::AbstractString, ax=nothing; kwargs...)
+function PyPlot.scatter(bd::BATLData{1, T}, var::AbstractString, ax=nothing; kwargs...) where T
    x, w = bd.x, bd.w
    varIndex_ = findindex(bd, var)
    if isnothing(ax) ax = plt.gca() end
@@ -457,9 +483,9 @@ end
 
 Wrapper over `contour` in matplotlib.
 """
-function PyPlot.contour(bd::BATLData, var::AbstractString, ax=nothing; levels=0,
-   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, innermask=false, kwargs...)
-   Xi, Yi, Wi = getdata2d(bd, var, plotrange, plotinterval; innermask)
+function PyPlot.contour(bd::BATLData{2, T}, var::AbstractString, ax=nothing; levels=0,
+   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, innermask=false, kwargs...) where T
+   Xi, Yi, Wi = slice2d(bd, var, plotrange, plotinterval; innermask)
    if isnothing(ax) ax = plt.gca() end
 
    if levels != 0
@@ -473,11 +499,11 @@ end
     contourf(data, var, ax=nothing; levels=0, plotrange=[-Inf,Inf,-Inf,Inf],
        plotinterval=0.1, innermask=false, kwargs...)
 
-Wrapper over `contourf` in matplotlib. See [`getdata2d`](@ref) for some related keywords.
+Wrapper over `contourf` in matplotlib. See [`slice2d`](@ref) for some related keywords.
 """
-function PyPlot.contourf(bd::BATLData, var::AbstractString, ax=nothing; levels::Int=0,
-   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, innermask=false, kwargs...)
-   Xi, Yi, Wi = getdata2d(bd, var, plotrange, plotinterval; innermask)
+function PyPlot.contourf(bd::BATLData{2, T}, var::AbstractString, ax=nothing; levels::Int=0,
+   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, innermask=false, kwargs...) where T
+   Xi, Yi, Wi = slice2d(bd, var, plotrange, plotinterval; innermask)
    if isnothing(ax) ax = plt.gca() end
 
    if levels != 0
@@ -492,8 +518,8 @@ end
 
 Wrapper over `tricontourf` in matplotlib.
 """
-function PyPlot.tricontourf(bd::BATLData, var::AbstractString, ax=nothing;
-   plotrange=[-Inf,Inf,-Inf,Inf], kwargs...)
+function PyPlot.tricontourf(bd::BATLData{2, T}, var::AbstractString, ax=nothing;
+   plotrange=[-Inf,Inf,-Inf,Inf], kwargs...) where T
    x, w = bd.x, bd.w
    varIndex_ = findindex(bd, var)
 
@@ -520,8 +546,8 @@ end
 
 Wrapper over `plot_trisurf` in matplotlib.
 """
-function PyPlot.plot_trisurf(bd::BATLData, var::AbstractString, ax=nothing;
-   plotrange=[-Inf,Inf,-Inf,Inf], kwargs...)
+function PyPlot.plot_trisurf(bd::BATLData{2, T}, var::AbstractString, ax=nothing;
+   plotrange=[-Inf,Inf,-Inf,Inf], kwargs...) where T
    x, w = bd.x, bd.w
    varIndex_ = findindex(bd, var)
 
@@ -549,9 +575,9 @@ end
 
 Wrapper over `plot_surface` in matplotlib.
 """
-function PyPlot.plot_surface(bd::BATLData, var::AbstractString;
-   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, innermask=false, kwargs...)
-   xi, yi, Wi = getdata2d(bd, var, plotrange, plotinterval; innermask)
+function PyPlot.plot_surface(bd::BATLData{2, T}, var::AbstractString;
+   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, innermask=false, kwargs...) where T
+   xi, yi, Wi = slice2d(bd, var, plotrange, plotinterval; innermask)
    Xi, Yi = meshgrid(xi, yi)
 
    plot_surface(Xi, Yi, Wi; kwargs...)
@@ -564,9 +590,9 @@ end
 
 Wrapper over `pcolormesh` in matplotlib.
 """
-function PyPlot.pcolormesh(bd::BATLData, var::AbstractString, ax=nothing;
-   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, innermask=false, kwargs...)
-   xi, yi, Wi = getdata2d(bd, var, plotrange, plotinterval; innermask)
+function PyPlot.pcolormesh(bd::BATLData{2, T}, var::AbstractString, ax=nothing;
+   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, innermask=false, kwargs...) where T
+   xi, yi, Wi = slice2d(bd, var, plotrange, plotinterval; innermask)
 
    if isnothing(ax) ax = plt.gca() end
 
@@ -580,9 +606,9 @@ end
 
 Wrapper over `tripcolor` in matplotlib.
 """
-function PyPlot.tripcolor(bd::BATLData, var::AbstractString, ax=nothing;
-   plotrange=[-Inf,Inf,-Inf,Inf], innermask=false, kwargs...)
-   x, w, ndim = bd.x, bd.w, bd.head.ndim
+function PyPlot.tripcolor(bd::BATLData{2, T}, var::AbstractString, ax=nothing;
+   plotrange=[-Inf,Inf,-Inf,Inf], innermask=false, kwargs...) where {T}
+   x, w = bd.x, bd.w
 
    varIndex_ = findindex(bd, var)
 
@@ -597,7 +623,7 @@ function PyPlot.tripcolor(bd::BATLData, var::AbstractString, ax=nothing;
    if innermask
       varIndex_ = findlast(x->x=="rbody", bd.head.variables)
       isnothing(varIndex_) && error("rbody not found in file header parameters!")
-      ParamIndex_ = varIndex_ - ndim - bd.head.nw
+      ParamIndex_ = varIndex_ - 2 - bd.head.nw
       r2 = bd.head.eqpar[ParamIndex_]^2
 
       ids = triang.triangles .+ Int32(1)
@@ -627,8 +653,8 @@ end
 
 Wrapper over `streamplot` in matplotlib .
 """
-function PyPlot.streamplot(bd::BATLData, var::AbstractString, ax=nothing;
-   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, kwargs...)
+function PyPlot.streamplot(bd::BATLData{2, T}, var::AbstractString, ax=nothing;
+   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, kwargs...) where T
    x, w = bd.x, bd.w
    varstream = split(var, ";")
    wnames = lowercase.(bd.head.wnames)
@@ -685,8 +711,8 @@ end
 
 Wrapper over `quiver` in matplotlib. Only supports Cartesian grid for now.
 """
-function PyPlot.quiver(bd::BATLData, var::AbstractString, ax=nothing;
-   stride::Integer=10, kwargs...)
+function PyPlot.quiver(bd::BATLData{2, T}, var::AbstractString, ax=nothing;
+   stride::Integer=10, kwargs...) where T
    x, w = bd.x, bd.w
    VarQuiver  = split(var, ";")
    var1_ = findindex(bd, VarQuiver[1])
diff --git a/src/plot/utility.jl b/src/plot/utility.jl
index 780977b9..7f3921c8 100644
--- a/src/plot/utility.jl
+++ b/src/plot/utility.jl
@@ -1,7 +1,7 @@
 # Utility functions for plotting.
 
 """
-    getdata2d(bd::BATLData, var::AbstractString, plotrange=[-Inf, Inf, -Inf, Inf],
+    slice2d(bd::BATLData, var::AbstractString, plotrange=[-Inf, Inf, -Inf, Inf],
        plotinterval=Inf; kwargs...)
 
 Return 2D slices of data `var` from `bd`. If `plotrange` is not set, output data resolution
@@ -12,12 +12,10 @@ is the same as the original.
 - `rbody=1.0`: Radius of the inner mask. Used when the rbody parameter is not found in the header.
 - `useMatplotlib=true`: Whether to Matplotlib (somehow faster) or NaturalNeighbours for scattered interpolation.
 """
-function getdata2d(bd::BATLData, var::AbstractString,
+function slice2d(bd::BATLData{2, T}, var::AbstractString,
    plotrange::Vector=[-Inf, Inf, -Inf, Inf], plotinterval::Real=Inf;
-   innermask::Bool=false, rbody::Real=1.0, useMatplotlib::Bool=true)
-   x, w, ndim = bd.x, bd.w, bd.head.ndim
-   @assert ndim == 2 "data must be in 2D!"
-
+   innermask::Bool=false, rbody::Real=1.0, useMatplotlib::Bool=true) where {T}
+   x, w = bd.x, bd.w
    varIndex_ = findindex(bd, var)
 
    if bd.head.gencoord # Generalized coordinates
@@ -58,9 +56,9 @@ function getdata2d(bd::BATLData, var::AbstractString,
             xi = range(plotrange[1], stop=plotrange[2], step=plotinterval)
             yi = range(plotrange[3], stop=plotrange[4], step=plotinterval)
          end
-
-         interp = @views cubic_spline_interpolation((xrange, yrange), w[:,:,varIndex_])
-         Wi = [interp(i, j) for j in yi, i in xi]
+         itp = @views scale(interpolate(w[:,:,varIndex_], BSpline(Linear())),
+            (xrange, yrange))
+         Wi = [itp(i, j) for j in yi, i in xi]
       end
    end
 
@@ -75,6 +73,7 @@ function getdata2d(bd::BATLData, var::AbstractString,
             end
          end
       else
+         ndim = 2
          ParamIndex_ = varIndex_ - ndim - bd.head.nw
          @inbounds @simd for i in CartesianIndices(Wi)
             if xi[i[1]]^2 + yi[i[2]]^2 < bd.head.eqpar[ParamIndex_]^2
@@ -82,7 +81,6 @@ function getdata2d(bd::BATLData, var::AbstractString,
             end
          end
       end
-
    end
 
    xi, yi, Wi
@@ -105,11 +103,7 @@ function meshgrid(x, y)
 end
 
 @inline function hasunit(bd::BATLData)
-   if startswith(bd.head.headline, "normalized")
-      return false
-   else
-      return true
-   end
+   startswith(bd.head.headline, "normalized") ? false : true
 end
 
 "Adjust 2D plot ranges."
@@ -124,12 +118,67 @@ end
 
 "Perform Triangle interpolation of 2D data `W` on grid `X`, `Y`."
 function interpolate2d_generalized_coords(X::T, Y::T, W::T,
-   plotrange::Vector{<:AbstractFloat}, plotinterval::Real) where
-   {T<:AbstractVector}
+   plotrange::Vector{<:AbstractFloat}, plotinterval::Real) where {T<:AbstractVector}
    xi = range(plotrange[1], stop=plotrange[2], step=plotinterval)
    yi = range(plotrange[3], stop=plotrange[4], step=plotinterval)
-   itp = interpolate(X, Y, W)
-   Wi = [itp(x, y; method=Triangle()) for y in yi, x in xi]::Matrix{eltype(W)}
+   itp = NN.interpolate(X, Y, W)
+   Wi = [itp(x, y; method=NN.Triangle()) for y in yi, x in xi]::Matrix{eltype(W)}
 
    xi, yi, Wi
-end
\ No newline at end of file
+end
+
+"""
+    interp1d(bd::BATLData, var::AbstractString, loc::AbstractVector{<:AbstractFloat})
+
+Interpolate `var` at spatial point `loc` in `bd`.
+"""
+function interp1d(bd::BATLData{2, T}, var::AbstractString, loc::AbstractVector{<:AbstractFloat}) where {T}
+   @assert !bd.head.gencoord "Only accept structured grids!"
+
+   x = bd.x
+   v = getview(bd, var)
+   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))
+   itp = scale(interpolate(v, BSpline(Linear())), (xrange, yrange))
+
+   Wi = itp(loc...)
+end
+
+"""
+    slice1d(bd::BATLData, var::AbstractString, point1::Vector, point2::Vecto)
+
+Interpolate `var` along a line from `point1` to `point2` in `bd`.
+"""
+function slice1d(bd::BATLData{2, T}, var::AbstractString, point1::Vector, point2::Vector) where {T}
+   @assert !bd.head.gencoord "Only accept structured grids!"
+
+   x = bd.x
+   v = getview(bd, var)
+   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))
+   itp = scale(interpolate(v, BSpline(Linear())), (xrange, yrange))
+   Δestimate = √(xrange.step^2 + yrange.step^2)
+   l = √(sum((point2 .- point1).^2))
+   ns = Int(l ÷ Δestimate)
+   dx = (point2[1] - point1[1]) / ns
+   dy = (point2[2] - point1[2]) / ns
+   points = [(point1[1] + i*dx, point1[2] + i*dy) for i in 0:ns]
+
+   Wi = [itp(loc...) for loc in points]
+end
+
+"Return view of variable `var` in `bd`."
+function getview(bd::BATLData{1, T}, var) where T
+   varIndex_ = findindex(bd, var)
+
+   v = @view bd.w[:,varIndex_]
+end
+
+function getview(bd::BATLData{2, T}, var) where T
+   varIndex_ = findindex(bd, var)
+
+   v = @view bd.w[:,:,varIndex_]
+end
+
+"Return value range of `var` in `bd`."
+get_var_range(bd::BATLData, var) = getview(bd, var) |> extrema
\ No newline at end of file
diff --git a/src/select.jl b/src/select.jl
index ffd4a58c..28f7a7d0 100644
--- a/src/select.jl
+++ b/src/select.jl
@@ -8,7 +8,6 @@ The returned 2D data lies in the `sequence` plane from - to + in `dir`.
 """
 function cutdata(bd::BATLData, var::AbstractString;
    plotrange=[-Inf,Inf,-Inf,Inf], dir::String="x", sequence::Int=1)
-   x, w = bd.x, bd.w
    var_ = findfirst(x->x==lowercase(var), lowercase.(bd.head.wnames))
    isempty(var_) && error("$(var) not found in header variables!")
 
@@ -20,9 +19,9 @@ function cutdata(bd::BATLData, var::AbstractString;
       dim, d1, d2 = 3, 1, 2
    end
 
-   cut1 = selectdim(view(x,:,:,:,d1), dim, sequence)
-   cut2 = selectdim(view(x,:,:,:,d2), dim, sequence)
-   W = selectdim(view(w,:,:,:,var_), dim, sequence)
+   cut1 = selectdim(view(bd.x,:,:,:,d1), dim, sequence)
+   cut2 = selectdim(view(bd.x,:,:,:,d2), dim, sequence)
+   W = selectdim(view(bd.w,:,:,:,var_), dim, sequence)
 
    if !all(isinf.(plotrange))
       cut1, cut2, W = subsurface(cut1, cut2, W, plotrange)
@@ -175,13 +174,13 @@ bd["rho"]
 
 See also: [`getvars`](@ref).
 """
-function getvar(bd::BATLData, var::AbstractString)
+function getvar(bd::BATLData{ndim, T}, var::AbstractString) where {ndim, T}
    if var in keys(variables_predefined)
       w = variables_predefined[var](bd)
    else
       var_ = findfirst(x->x==lowercase(var), lowercase.(bd.head.wnames))
       isnothing(var_) && error("$var not found in file header variables!")
-      w = selectdim(bd.w, bd.head.ndim+1, var_)
+      w = selectdim(bd.w, ndim+1, var_)
    end
 
    w
@@ -191,14 +190,14 @@ end
    getvar(bd, var)
 
 """
-    getvars(bd::BATLData, Names::Vector) -> Dict
+    getvars(bd::BATLData, names::Vector) -> Dict
 
-Return variables' data as a dictionary from string vector.
+Return variables' data as a dictionary from vector of `names`.
 See also: [`getvar`](@ref).
 """
-function getvars(bd::BATLData{U}, Names::Vector{T}) where {U, T<:AbstractString}
+function getvars(bd::BATLData{ndim, U}, names::Vector{T}) where {ndim, U, T<:AbstractString}
    dict = Dict{T, Array{U}}()
-   for name in Names
+   for name in names
       dict[name] = getvar(bd, name)
    end
 
@@ -206,11 +205,11 @@ function getvars(bd::BATLData{U}, Names::Vector{T}) where {U, T<:AbstractString}
 end
 
 "Construct vectors from scalar components."
-function _fill_vector_from_scalars(bd::BATLData, vstr1, vstr2, vstr3)
+function _fill_vector_from_scalars(bd::BATLData{ndim, T}, vstr1, vstr2, vstr3) where {ndim, T}
    v1 = getvar(bd, vstr1)
    v2 = getvar(bd, vstr2)
    v3 = getvar(bd, vstr3)
-   v = Array{eltype(v1), ndims(v1)+1}(undef, 3, size(v1)...)
+   v = Array{T, ndims(v1)+1}(undef, 3, size(v1)...)
 
    Rpost = CartesianIndices(size(v1))
    for Ipost in Rpost
diff --git a/test/runtests.jl b/test/runtests.jl
index ddd7fa90..22aec918 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -50,16 +50,26 @@ end
       @test extrema(bd.x) == (-127.5f0, 127.5f0)
       @test extrema(bd.w) == (-0.79985905f0, 1.9399388f0)
       plotrange = [-10.0, 10.0, -Inf, Inf]
-      x, y, w = Batsrus.getdata2d(bd, "rho", plotrange)
-      @test w[1,end] == 0.6848978549242021
+      x, y, w = slice2d(bd, "rho", plotrange)
+      @test w[1,end] == 0.6848635077476501
       @test bd["B"][:,end,end] == Float32[1.118034, -0.559017, 0.0]
+      # Linear interpolation at a given point
+      d = interp1d(bd, "rho", Float32[0.0, 0.0])
+      @test d == 0.6936918f0
+      # Linear interpolation along a line
+      point1 = Float32[-10.0, -1.0]
+      point2 = Float32[10.0, 1.0]
+      w = slice1d(bd, "rho", point1, point2)
+      @test sum(w) == 10.676028f0
+
+      @test get_var_range(bd, "rho") == (0.11626893f0, 1.0f0)
    end
 
    @testset "Reading 2D unstructured" begin
       file = "bx0_mhd_6_t00000100_n00000352.out"
       bd = load(joinpath(datapath, file))
       plotrange = [-Inf, Inf, -Inf, Inf]
-      x, y, w = Batsrus.getdata2d(bd, "rho", plotrange, useMatplotlib=false)
+      x, y, w = slice2d(bd, "rho", plotrange, useMatplotlib=false)
       @test w[1,2] == 5.000018304080387
    end
 
@@ -68,8 +78,7 @@ end
       bd = load(joinpath(datapath, file))
       plotrange = [-50.0, 50.0, -0.5, 0.5]
       X, Z, p = cutdata(bd, "p"; dir="y", sequence=1, plotrange)
-      @test p[1] ≈ 0.560976f0
-      @test p[2] ≈ 0.53704995f0
+      @test p[1] ≈ 0.560976f0 && p[2] ≈ 0.53704995f0
       @test size(bd["p"]) == (8,8,8)
       vars = getvars(bd, ["p"])
       @test size(vars["p"]) == (8,8,8)
@@ -96,8 +105,8 @@ end
       filetag = joinpath(datapath, "3d_mhd_amr/3d__mhd_1_t00000000_n00000000")
       batl = Batl(readhead(filetag*".info"), readtree(filetag)...)
       # local block index check
-      @test Batsrus.find_grid_block(batl, [1.0, 0.0, 0.0]) == 4
-      @test Batsrus.find_grid_block(batl, [100.0, 0.0, 0.0]) == -100
+      @test Batsrus.find_grid_block(batl, [1.0, 0.0, 0.0]) == 4 &&
+         Batsrus.find_grid_block(batl, [100.0, 0.0, 0.0]) == -100
       connectivity = getConnectivity(batl)
       sha_str = bytes2hex(sha256(string(connectivity)))
       @test sha_str == "c6c5a65a46d86a9ba4096228c1516f89275e45e295cd305eb70c281a770ede74"