Add E field computations; Add axis extraction

src/Batsrus.jl

@@ -12,7 +12,7 @@ using StaticArrays: SVector, @SMatrix, SA
 
 export BATS,
    load, readlogdata, readtecdata, showhead, # io
-   getvar, cutdata, subvolume, subsurface, # select
+   getvar, cutdata, subvolume, subsurface, get_convection_E, get_hall_E, # select
    Batl, convertTECtoVTU, convertIDLtoVTK, readhead, readtree, getConnectivity, # vtk
    interp1d, interp2d, slice1d, get_var_range, squeeze, get_range # plot/utility
 

src/select.jl

@@ -289,4 +289,43 @@ const variables_predefined = Dict{String, Function}(
    "U" => (bd -> _fill_vector_from_scalars(bd, :U)),
    "Anisotropy0" => (bd -> _get_anisotropy(bd, 0)),
    "Anisotropy1" => (bd -> _get_anisotropy(bd, 1)),
-)
+)
+
+"Return the convection electric field from PIC outputs."
+function get_convection_E(bd::BATS)
+   Bx, By, Bz = _get_vectors(bd, :B) # [nT]
+   # Let us use H+ velocities as the ion bulk velocity and ignore heavy ions
+   uix, uiy, uiz = _get_vectors(bd, :U1)
+
+   Econvx = similar(Bx)
+   Econvy = similar(By)
+   Econvz = similar(Bz)
+   # -Ui × B
+   @simd for i in eachindex(Econvx)
+      Econvx[i] = -uiy[i]*Bz[i] + uiz[i]*By[i]
+      Econvy[i] = -uiz[i]*Bx[i] + uix[i]*Bz[i]
+      Econvz[i] = -uix[i]*By[i] + uiy[i]*Bx[i]
+   end
+
+   Econvx, Econvy, Econvz
+end
+
+"Return the Hall electric field from PIC outputs."
+function get_hall_E(bd::BATS)
+   Bx, By, Bz = _get_vectors(bd, :B) # [nT]
+   uex, uey, uez = _get_vectors(bd, :U0) # [km/s]
+   # Let us use H+ velocities as the ion bulk velocity and ignore heavy ions
+   uix, uiy, uiz = _get_vectors(bd, :U1)
+
+   Ehallx = similar(Bx)
+   Ehally = similar(By)
+   Ehallz = similar(Bz)
+   # (Ui - Ue) × B
+   for i in eachindex(Ehallx)
+      Ehallx[i] = (uiy[i] - uey[i])*Bz[i] - (uiz[i] - uez[i])*By[i]
+      Ehally[i] = (uiz[i] - uez[i])*Bx[i] - (uix[i] - uex[i])*Bz[i]
+      Ehallz[i] = (uix[i] - uex[i])*By[i] - (uiy[i] - uey[i])*Bx[i]
+   end
+
+   Ehallx, Ehally, Ehallz
+end

src/utility.jl

@@ -89,6 +89,43 @@ function interp2d(bd::BATS{2, 3, T}, var::AbstractString,
    xi, yi, Wi
 end
 
+"Return the axis range for 2D outputs."
+function meshgrid(bd::BATS,
+   plotrange::Vector=[-Inf32, Inf32, -Inf32, Inf32], plotinterval::Real=Inf32)
+   x = bd.x
+
+   if bd.head.gencoord # Generalized coordinates
+      X, Y = eachslice(x, dims=3)
+      X, Y = vec(X), vec(Y)
+
+      adjust_plotrange!(plotrange, extrema(X), extrema(Y))
+      # Set a heuristic value if not set
+      if isinf(plotinterval)
+         plotinterval = (plotrange[2] - plotrange[1]) / size(X, 1)
+      end
+      xi = range(plotrange[1], stop=plotrange[2], step=plotinterval)
+      yi = range(plotrange[3], stop=plotrange[4], step=plotinterval)
+   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
+      else
+         adjust_plotrange!(plotrange, (xrange[1], xrange[end]), (yrange[1], yrange[end]))
+
+         if isinf(plotinterval)
+            xi = range(plotrange[1], stop=plotrange[2], step=xrange[2] - xrange[1])
+            yi = range(plotrange[3], stop=plotrange[4], step=yrange[2] - yrange[1])
+         else
+            xi = range(plotrange[1], stop=plotrange[2], step=plotinterval)
+            yi = range(plotrange[3], stop=plotrange[4], step=plotinterval)
+         end
+      end
+   end
+
+   xi, yi
+end
+
 "Find variable index in the BATSRUS data."
 function findindex(bd::BATS, var::AbstractString)
    varIndex_ = findfirst(x->lowercase(x)==lowercase(var), bd.head.wnames)

test/runtests.jl

@@ -69,13 +69,19 @@ end
 
       file = "z=0_fluid_region0_0_t00001640_n00010142.out"
       bd = load(joinpath(datapath, file))
+      x, y = Batsrus.meshgrid(bd)
+      @test length(x) == 601 && y[2] == 0.0f0
       @test bd["Emag"][2,1] == 2655.4805f0
       @test bd["E2"][2,1] == 7.051577f6
       @test bd["E"][:,2,1] == Float32[-241.05942, -2644.2058, -40.53219]
       @test bd["Ue2"][2,1] == 33784.973f0
       @test bd["Ui2"][2,1] == bd["Ui2"][2,1]
       @test bd["Anisotropy0"][1:2,1] ≈ Float32[1.2630985, 2.4700143]
       @test bd["Anisotropy1"][1:2,1] ≈ Float32[1.2906302, 2.6070855]
+      w = get_convection_E(bd)
+      @test w[2][2,1] ≈ -2454.3933f0
+      w = get_hall_E(bd)
+      @test w[2][2,1] ≈ -782.2945f0
    end
 
    @testset "Reading 2D unstructured ascii" begin