From 26ce83539106ec07199fa5c519ab198d6c13d6f9 Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Fri, 28 May 2021 12:54:34 +0300
Subject: [PATCH 01/25] Initial commit of divS helpers. Introduced fsgrid
 reader that crudely accounts for variable centering for edge and face
 vectors.

---
 pyPlots/plot_colormap.py |  2 +-
 pyPlots/plot_helpers.py  |  8 +++++++-
 pyVlsv/vlsvreader.py     | 32 ++++++++++++++++++++++++++++++++
 3 files changed, 40 insertions(+), 2 deletions(-)

diff --git a/pyPlots/plot_colormap.py b/pyPlots/plot_colormap.py
index 5eb876bd..3be73b69 100644
--- a/pyPlots/plot_colormap.py
+++ b/pyPlots/plot_colormap.py
@@ -816,7 +816,7 @@ def exprMA_cust(exprmaps, requestvariables=False):
     if lin is None:
         # Special SymLogNorm case
         if symlog is not None:
-            if LooseVersion(matplotlib.__version__) < LooseVersion("3.3.0"):
+            if LooseVersion(matplotlib.__version__) < LooseVersion("3.2.0"):
                 norm = SymLogNorm(linthresh=linthresh, linscale = 1.0, vmin=vminuse, vmax=vmaxuse, clip=True)
                 print("WARNING: colormap SymLogNorm uses base-e but ticks are calculated with base-10.")
                 #TODO: copy over matplotlib 3.3.0 implementation of SymLogNorm into pytools/analysator
diff --git a/pyPlots/plot_helpers.py b/pyPlots/plot_helpers.py
index 2a22b4c4..c092e82a 100644
--- a/pyPlots/plot_helpers.py
+++ b/pyPlots/plot_helpers.py
@@ -26,7 +26,7 @@
 import sys
 from rotation import rotateTensorToVector
 
-PLANE = 'XY'
+PLANE = 'XZ'
 # or alternatively, 'XZ'
 CELLSIZE = 300000.0 # cell size
 DT = 0.5 # time step
@@ -661,6 +661,12 @@ def expr_flowcompression(pass_maps, requestvariables=False):
         return ['V']
     Vmap = TransposeVectorArray(pass_maps['V']) # Bulk flow
     return numdiv(Vmap).T
+    
+def expr_divPoynting(pass_maps, requestvariables=False):
+    if requestvariables==True:
+        return ['poynting']
+    dSmap = TransposeVectorArray(pass_maps['poynting']) # Bulk flow
+    return numdiv(dSmap).T
 
 # def expr_gradB_aniso(pass_maps):
 #     Bmap = TransposeVectorArray(pass_maps['B']) # Magnetic field
diff --git a/pyVlsv/vlsvreader.py b/pyVlsv/vlsvreader.py
index 48e910ef..0c8a0f15 100644
--- a/pyVlsv/vlsvreader.py
+++ b/pyVlsv/vlsvreader.py
@@ -1027,6 +1027,38 @@ def calcLocalSize(globalCells, ntasks, my_n):
 
        return np.squeeze(orderedData)
 
+   def read_fg_variable_as_volumetric(self, name, centering="default", operator="pass"):
+      fgdata = self.read_fsgrid_variable(name, operator)
+
+      fssize=list(self.get_fsgrid_mesh_size())
+      fgdata=np.expand_dims(fgdata, fssize.index(1)) #expand to have a singleton dimension for a reduced dim - lets roll happen with ease
+      print('read in fgdata with shape', fgdata.shape, name)
+      celldata = np.zeros_like(fgdata)
+      known_centerings = {"fg_b":"face", "fg_e":"edge"}
+      if name.casefold() in known_centerings.keys() and centering == "default":
+         centering = known_centerings[name.casefold()]
+      else:
+         print("A variable with unknown centering! Aborting.")
+         return False
+      print('fgdata.shape', fgdata.shape)
+      #vector variable
+      if fgdata.shape[-1] == 3:
+         if centering=="face":
+            celldata[:,:,:,0] = (fgdata[:,:,:,0] + np.roll(fgdata[:,:,:,0],-1, 0))/2.0
+            celldata[:,:,:,1] = (fgdata[:,:,:,1] + np.roll(fgdata[:,:,:,1],-1, 1))/2.0
+            celldata[:,:,:,2] = (fgdata[:,:,:,2] + np.roll(fgdata[:,:,:,2],-1, 2))/2.0
+         elif centering=="edge":
+            celldata[:,:,:,0] = (fgdata[:,:,:,0] + np.roll(fgdata[:,:,:,0],-1, 1) + np.roll(fgdata[:,:,:,0],-1, 2) + np.roll(fgdata[:,:,:,0],-1, (1,2)))/4.0
+            celldata[:,:,:,1] = (fgdata[:,:,:,1] + np.roll(fgdata[:,:,:,1],-1, 0) + np.roll(fgdata[:,:,:,1],-1, 2) + np.roll(fgdata[:,:,:,1],-1, (0,2)))/4.0
+            celldata[:,:,:,2] = (fgdata[:,:,:,2] + np.roll(fgdata[:,:,:,2],-1, 0) + np.roll(fgdata[:,:,:,2],-1, 1) + np.roll(fgdata[:,:,:,2],-1, (0,1)))/4.0
+         else:
+            print("Unknown centering ('" +centering+ "')! Aborting.")
+            return False
+      else:
+         print("A scalar variable! I don't know what to do with this! Aborting.")
+         return False
+      return celldata
+
    def read_variable(self, name, cellids=-1,operator="pass"):
       ''' Read variables from the open vlsv file. 
       Arguments:

From 2370fc2899120089eb16f546a67ea34ce72091b7 Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Tue, 1 Jun 2021 13:20:05 +0300
Subject: [PATCH 02/25] Reduced dimension handling

---
 pyVlsv/vlsvreader.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/pyVlsv/vlsvreader.py b/pyVlsv/vlsvreader.py
index 0c8a0f15..bc996e6a 100644
--- a/pyVlsv/vlsvreader.py
+++ b/pyVlsv/vlsvreader.py
@@ -1031,7 +1031,8 @@ def read_fg_variable_as_volumetric(self, name, centering="default", operator="pa
       fgdata = self.read_fsgrid_variable(name, operator)
 
       fssize=list(self.get_fsgrid_mesh_size())
-      fgdata=np.expand_dims(fgdata, fssize.index(1)) #expand to have a singleton dimension for a reduced dim - lets roll happen with ease
+      if 1 in fssize:
+         fgdata=np.expand_dims(fgdata, fssize.index(1)) #expand to have a singleton dimension for a reduced dim - lets roll happen with ease
       print('read in fgdata with shape', fgdata.shape, name)
       celldata = np.zeros_like(fgdata)
       known_centerings = {"fg_b":"face", "fg_e":"edge"}

From 893d0c2ffab350df9f2136f0f29a07a44b33b4fa Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Fri, 4 Jun 2021 11:05:36 +0300
Subject: [PATCH 03/25] Note for using reconstructions

---
 pyVlsv/vlsvreader.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/pyVlsv/vlsvreader.py b/pyVlsv/vlsvreader.py
index bc996e6a..c27bf666 100644
--- a/pyVlsv/vlsvreader.py
+++ b/pyVlsv/vlsvreader.py
@@ -1048,6 +1048,7 @@ def read_fg_variable_as_volumetric(self, name, centering="default", operator="pa
             celldata[:,:,:,0] = (fgdata[:,:,:,0] + np.roll(fgdata[:,:,:,0],-1, 0))/2.0
             celldata[:,:,:,1] = (fgdata[:,:,:,1] + np.roll(fgdata[:,:,:,1],-1, 1))/2.0
             celldata[:,:,:,2] = (fgdata[:,:,:,2] + np.roll(fgdata[:,:,:,2],-1, 2))/2.0
+            # Use Leo's reconstuction for fg_b instead
          elif centering=="edge":
             celldata[:,:,:,0] = (fgdata[:,:,:,0] + np.roll(fgdata[:,:,:,0],-1, 1) + np.roll(fgdata[:,:,:,0],-1, 2) + np.roll(fgdata[:,:,:,0],-1, (1,2)))/4.0
             celldata[:,:,:,1] = (fgdata[:,:,:,1] + np.roll(fgdata[:,:,:,1],-1, 0) + np.roll(fgdata[:,:,:,1],-1, 2) + np.roll(fgdata[:,:,:,1],-1, (0,2)))/4.0

From 617573b11cd648a7fa21b56849ec1c2b60993ae2 Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Tue, 22 Jun 2021 13:54:39 +0300
Subject: [PATCH 04/25] an initial, dysfunctional divS function and
 derivatives.py file for future use

---
 pyCalculations/derivatives.py | 52 +++++++++++++++++++++++++++++++++++
 1 file changed, 52 insertions(+)
 create mode 100644 pyCalculations/derivatives.py

diff --git a/pyCalculations/derivatives.py b/pyCalculations/derivatives.py
new file mode 100644
index 00000000..3d91e6b7
--- /dev/null
+++ b/pyCalculations/derivatives.py
@@ -0,0 +1,52 @@
+# 
+# This file is part of Analysator.
+# Copyright 2013-2016 Finnish Meteorological Institute
+# Copyright 2017-2018 University of Helsinki
+# 
+# For details of usage, see the COPYING file and read the "Rules of the Road"
+# at http://www.physics.helsinki.fi/vlasiator/
+# 
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+# 
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+# 
+# You should have received a copy of the GNU General Public License along
+# with this program; if not, write to the Free Software Foundation, Inc.,
+# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+# 
+
+import numpy as np
+import scipy as sp
+import pytools as pt
+import warnings
+from scipy import interpolate
+
+def get_fg_divPoynting:
+   b = bulkReader.read_fg_variable_as_volumetric('fg_b')
+   print('b.shape=',b.shape)
+   e = bulkReader.read_fg_variable_as_volumetric('fg_e')
+   print('e.shape=',e.shape)
+
+   mu_0 = 1.25663706144e-6
+   S = np.cross(e,b)/mu_0
+   print('S.shape', S.shape)
+
+   #divS = pt.plot.plot_helpers.numdiv(pt.plot.plot_helpers.TransposeVectorArray(np.squeeze(S))).T
+   divS = (np.roll(S[:,:,:,0],-1, 0) - np.roll(S[:,:,:,0], 1, 0) +
+         np.roll(S[:,:,:,1],-1, 1) - np.roll(S[:,:,:,1], 1, 1) +
+         np.roll(S[:,:,:,2],-1, 2) - np.roll(S[:,:,:,2], 1, 2))/(2*1e6)
+   print('divS.shape', divS.shape)
+   #divS = np.div
+
+   fig=plt.pyplot.figure()
+   ax=plt.pyplot.axes()
+   im=ax.imshow(divS[:,int(divS.shape[1]/2),:], origin='lower', norm=plt.colors.SymLogNorm(base=10, linthresh=1e-12, linscale = 1.0, vmin=-1e-7, vmax=1e-7, clip=True), cmap='vik')
+   plt.pyplot.colorbar(im, ax=ax)
+   plt.pyplot.savefig('./divS_centered.png')
+   print(divS)
\ No newline at end of file

From cdb3ad299f556c1058232fcda8a72b4b0c9cf654 Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Thu, 19 Aug 2021 14:34:33 +0300
Subject: [PATCH 05/25] ElementTree.SubElement not accepting kwords, only
 arguments. Added vlsvWriter functions to copy given variables instead of all
 variables to a new file.

---
 pyVlsv/vlsvwriter.py | 49 ++++++++++++++++++++++++++++++++++++++++----
 1 file changed, 45 insertions(+), 4 deletions(-)

diff --git a/pyVlsv/vlsvwriter.py b/pyVlsv/vlsvwriter.py
index 222d7103..8e0a53e9 100644
--- a/pyVlsv/vlsvwriter.py
+++ b/pyVlsv/vlsvwriter.py
@@ -90,10 +90,15 @@ def __initialize( self, vlsvReader ):
             self.write( data=data, name=name, tag=tag, mesh=mesh, extra_attribs=extra_attribs )
 
 
-   def copy_variables( self, vlsvReader ):
-      ''' Copies all variables from vlsv reader to the file
-
+   def copy_variables( self, vlsvReader, varlist=None ):
+      ''' Copies variables from vlsv reader to the file.
+           varlist = None: list of variables to copy; 
       '''
+
+      # Delegate to the variable list handler
+      if (varlist is not None):
+         self.copy_variables_list(vlsvReader, varlist)
+
       # Get the xml sheet:
       xml_root = vlsvReader._VlsvReader__xml_root
 
@@ -123,6 +128,41 @@ def copy_variables( self, vlsvReader ):
             self.write( data=data, name=name, tag=tag, mesh=mesh, extra_attribs=extra_attribs )
       return
 
+   def copy_variables_list( self, vlsvReader, vars ):
+      ''' Copies variables in the list vars from vlsv reader to the file
+
+      '''
+      # Get the xml sheet:
+      xml_root = vlsvReader._VlsvReader__xml_root
+
+      # Get list of tags to write:
+      tags = {}
+      tags['VARIABLE'] = ''
+
+      # Copy the xml root and write variables
+      for child in xml_root:
+         if child.tag in tags:
+            if 'name' in child.attrib:
+                name = child.attrib['name']
+                if not name in vars:
+                   continue
+            else:
+                continue
+            if 'mesh' in child.attrib:
+                mesh = child.attrib['mesh']
+            else:
+                mesh = ''
+            tag = child.tag
+            # Copy extra attributes:
+            extra_attribs = {}
+            for i in child.attrib.items():
+               if i[0] != 'name' and i[0] != 'mesh':
+                  extra_attribs[i[0]] = i[1]
+            data = vlsvReader.read( name=name, tag=tag, mesh=mesh )
+            # Write the data:
+            self.write( data=data, name=name, tag=tag, mesh=mesh, extra_attribs=extra_attribs )
+      return
+
    def write_velocity_space( self, vlsvReader, cellid, blocks_and_values ):
       ''' Writes given velocity space into vlsv file
 
@@ -162,7 +202,8 @@ def write(self, data, name, tag, mesh, extra_attribs={}):
       datatype = ''
 
       # Add the data into the xml data:
-      child = ET.SubElement(parent=self.__xml_root, tag=tag)
+      print(self.__xml_root, tag)
+      child = ET.SubElement(self.__xml_root, tag)
       child.attrib["name"] = name
       child.attrib["mesh"] = mesh
       child.attrib["arraysize"] = len(np.atleast_1d(data))

From 16ef53717bc70a8fc811067e3254d902a6f9e83f Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Mon, 23 Aug 2021 20:06:50 +0300
Subject: [PATCH 06/25] Functional fg to vg upsampling, todo: diffs beyond
 eyeballing, downsampling routines, writer examples

---
 pyVlsv/vlsvreader.py | 128 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 128 insertions(+)

diff --git a/pyVlsv/vlsvreader.py b/pyVlsv/vlsvreader.py
index a12d4bab..a071a484 100644
--- a/pyVlsv/vlsvreader.py
+++ b/pyVlsv/vlsvreader.py
@@ -1236,6 +1236,91 @@ def get_amr_level(self,cellid):
          AMR_count += 1
       return AMR_count - 1 
 
+   def get_cell_dx(self, cellid):
+      '''Returns the dx of a given cell defined by its cellid
+      
+      :param cellid:        The cell's cellid
+      :returns:             The cell's size [dx, dy, dz]
+      '''
+      return np.array([self.__dx,self.__dy,self.__dz])/2**self.get_amr_level(cellid)
+
+   def get_cell_bbox(self, cellid):
+      '''Returns the bounding box of a given cell defined by its cellid
+      
+      :param cellid:        The cell's cellid
+      :returns:             The cell's bbox [xmin,ymin,zmin],[xmax,ymax,zmax]
+      '''
+      
+      hdx = self.get_cell_dx(cellid)*0.5
+      mid = self.get_cell_coordinates(cellid)
+      #print('halfdx:', hdx, 'mid:', mid, 'low:', mid-hdx, 'hi:', mid+hdx)
+      return mid-hdx, mid+hdx
+
+   def get_cell_fsgrid_slicemap(self, cellid):
+      '''Returns a slice tuple of fsgrid indices that are contained in the SpatialGrid
+      cell.
+      '''
+      low, up = self.get_cell_bbox(cellid)
+      lowi, upi = self.get_fsgrid_slice_indices(low, up)
+      return lowi, upi
+
+   def get_fsgrid_subarray(self, cellid, array):
+      '''Returns a subarray of the fsgrid array, corresponding to the fsgrid
+      covered by the SpatialGrid cellid. [untested]
+      '''
+      lowi, upi = self.get_cell_fsgrid_slicemap(cellid)
+      #print('subarray:',lowi, upi)
+      if array.ndim == 4:
+         return array[lowi[0]:upi[0]+1, lowi[1]:upi[1]+1, lowi[2]:upi[2]+1, :]
+      else:
+         return array[lowi[0]:upi[0]+1, lowi[1]:upi[1]+1, lowi[2]:upi[2]+1]
+
+
+   def downsample_fsgrid_subarray(self, cellid, array):
+      '''Returns a mean value of fsgrid values underlying the SpatialGrid cellid.
+      [untested]
+      '''
+      fsarr = self.get_fsgrid_subarray(cellid, array)
+      n = fsarr.size
+      if fsarr.ndim == 4:
+         n = n/3
+      ncells = 8**(self.get_max_refinement_level()-self.get_amr_level(cellid))
+      if(n != ncells):
+         print("Warning: weird fs subarray size", n, 'for amrlevel', self.get_amr_level(cellid), 'expect', ncells)
+      return np.mean(fsarr,axis=(0,1,2))
+
+   def upsample_fsgrid_subarray(self, cellid, var, array):
+      '''Set the elements of the fsgrid array to the value of corresponding SpatialGrid
+      cellid. Mutator.
+      [untested]
+      '''
+      lowi, upi = self.get_fsgrid_slicemap(cellid)
+      value = self.read_variable(var, cellids=[cellid])
+      if array.ndim == 4:
+         array[lowi[0]:upi[0]+1,lowi[1]:upi[1]+1,lowi[2]:upi[2]+1,:] = value
+      else:
+         array[lowi[0]:upi[0]+1,lowi[1]:upi[1]+1,lowi[2]:upi[2]+1] = value
+      return
+
+   def get_cell_fsgrid(self, cellid):
+      '''Returns a slice tuple of fsgrid indices that are contained in the SpatialGrid
+      cell.
+      '''
+      low, up = self.get_cell_bbox(cellid)
+      lowi, upi = self.get_fsgrid_slice_indices(low, up)
+      return lowi, upi
+
+   def get_fsgrid_coordinates(self, ri):
+      '''Returns real-space center coordinates of the fsgrid 3-index.
+      '''
+      lowerlimit = self.get_fsgrid_mesh_extent()[0:3]
+      upperlimit = self.get_fsgrid_mesh_extent()[3:6]
+      delta = upperlimit-lowerlimit
+      sz=self.get_fsgrid_mesh_size()
+      dxs = delta/sz
+
+      return lowerlimit+dxs*(np.array(ri)+0.5)
+
    def get_unique_cellids(self, coords):
       ''' Returns a list of cellids containing all the coordinates in coords,
           with no duplicate cellids. Relative order of elements is conserved.
@@ -1833,6 +1918,49 @@ def get_fsgrid_mesh_extent(self):
       '''
       return np.array([self.__xmin, self.__ymin, self.__zmin, self.__xmax, self.__ymax, self.__zmax])
 
+   def get_fsgrid_indices(self, coords):
+      ''' Convert spatial coordinates coords to an index array [xi, yi, zi] for fsgrid
+
+      :returns 3-tuple of integers [xi, yi, zi] corresponding to fsgrid cell containing coords (low-inclusive)
+      Example:
+      ii = f.get_fsgrid_mesh_extent(coords)
+      fsvar_at_coords = fsvar_array.item(ii)
+      '''
+      lower = self.get_fsgrid_mesh_extent()[0:3]
+      upper = self.get_fsgrid_mesh_extent()[3:6]
+      delta = upper-lower
+      sz=self.get_fsgrid_mesh_size()
+      dx = delta/sz
+      r0 = coords-lower
+      ri = np.floor(r0/dx).astype(int)
+      if (ri < 0).any() or (ri>sz-1).any():
+         print("get_fsgrid_indices: Resulting index out of bounds, returning None")
+         return None
+      return tuple(ri)
+
+   def get_fsgrid_slice_indices(self, lower, upper, eps=1e-3):
+      ''' Get indices for a subarray of an fsgrid variable, in the cuboid from lower to upper.
+      This is meant for mapping a set of fsgrid cells to a given SpatialGrid cell.
+      Shifts the corners (lower, upper) by dx_fsgrid*eps inward, if direct low-inclusive behaviour
+      is required, set kword eps = 0.
+
+
+      :returns two 3-tuples of integers.
+      Example:
+      ii = f.get_fsgrid_mesh_extent(coords)
+      fsvar_at_coords = fsvar_array.item(ii)
+      '''
+      lowerlimit = self.get_fsgrid_mesh_extent()[0:3]
+      upperlimit = self.get_fsgrid_mesh_extent()[3:6]
+      delta = upperlimit-lowerlimit
+      sz=self.get_fsgrid_mesh_size()
+      dx = delta/sz
+      eps = dx*eps
+      loweri = self.get_fsgrid_indices(lower+eps)
+      upperi = self.get_fsgrid_indices(upper-eps)
+      return loweri, upperi
+      
+
    def get_velocity_mesh_size(self, pop="proton"):
       ''' Read velocity mesh size
       

From 0b110fed12a2a3bbfbd32a9e24504ca2942ae2e6 Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Wed, 25 Aug 2021 17:33:16 +0300
Subject: [PATCH 07/25] fg to vg writer function in VlsvWriter, testing and
 upsampling-downsampling routines (nb. no filtering yet)

---
 pyVlsv/vlsvreader.py | 21 +++++++++++++++++++--
 pyVlsv/vlsvwriter.py | 11 ++++++++++-
 2 files changed, 29 insertions(+), 3 deletions(-)

diff --git a/pyVlsv/vlsvreader.py b/pyVlsv/vlsvreader.py
index a071a484..02784830 100644
--- a/pyVlsv/vlsvreader.py
+++ b/pyVlsv/vlsvreader.py
@@ -1278,7 +1278,6 @@ def get_fsgrid_subarray(self, cellid, array):
 
    def downsample_fsgrid_subarray(self, cellid, array):
       '''Returns a mean value of fsgrid values underlying the SpatialGrid cellid.
-      [untested]
       '''
       fsarr = self.get_fsgrid_subarray(cellid, array)
       n = fsarr.size
@@ -1294,14 +1293,32 @@ def upsample_fsgrid_subarray(self, cellid, var, array):
       cellid. Mutator.
       [untested]
       '''
-      lowi, upi = self.get_fsgrid_slicemap(cellid)
+      lowi, upi = self.get_cell_fsgrid_slicemap(cellid)
+      #print(lowi,upi)
       value = self.read_variable(var, cellids=[cellid])
       if array.ndim == 4:
+         #print(value)
          array[lowi[0]:upi[0]+1,lowi[1]:upi[1]+1,lowi[2]:upi[2]+1,:] = value
+         #print(array[lowi[0]:upi[0]+1,lowi[1]:upi[1]+1,lowi[2]:upi[2]+1,:])
       else:
          array[lowi[0]:upi[0]+1,lowi[1]:upi[1]+1,lowi[2]:upi[2]+1] = value
       return
 
+   def read_variable_as_fg(self, var):
+      cellids = self.read_variable('CellID')
+      sz = self.get_fsgrid_mesh_size()
+      testvar = self.read_variable(var, [cellids[0]])
+      varsize = testvar.size
+      if(varsize > 1):
+         fgarr = np.zeros([sz[0], sz[1], sz[2], varsize], dtype=testvar.dtype)
+      else:
+         fgarr = np.zeros(sz, dtype=testvar.dtype)
+      for c in cellids:
+         self.upsample_fsgrid_subarray(c, var, fgarr)
+         #print
+      return fgarr
+
+
    def get_cell_fsgrid(self, cellid):
       '''Returns a slice tuple of fsgrid indices that are contained in the SpatialGrid
       cell.
diff --git a/pyVlsv/vlsvwriter.py b/pyVlsv/vlsvwriter.py
index 8e0a53e9..e8eab964 100644
--- a/pyVlsv/vlsvwriter.py
+++ b/pyVlsv/vlsvwriter.py
@@ -202,7 +202,6 @@ def write(self, data, name, tag, mesh, extra_attribs={}):
       datatype = ''
 
       # Add the data into the xml data:
-      print(self.__xml_root, tag)
       child = ET.SubElement(self.__xml_root, tag)
       child.attrib["name"] = name
       child.attrib["mesh"] = mesh
@@ -248,6 +247,16 @@ def write(self, data, name, tag, mesh, extra_attribs={}):
       # write the xml footer:
       self.__write_xml_footer()
 
+   def write_fgarray_to_SpatialGrid(self, reader, data, name, extra_attribs={}):
+      # get a reader for the target file
+      #print(data.shape[0:3], reader.get_fsgrid_mesh_size(), (data.shape[0:3] == reader.get_fsgrid_mesh_size()))
+      if not (data.shape[0:3] == reader.get_fsgrid_mesh_size()).all():
+         print("Data shape does not match target fsgrid mesh")
+         return
+      cids = reader.read_variable("CellID")
+      vgdata = [reader.downsample_fsgrid_subarray(cid, data) for cid in cids]
+      self.write(vgdata, name, "VARIABLE", "SpatialGrid",extra_attribs)
+
    def __write_xml_footer( self ):
       # Write the xml footer:
       max_xml_size = 1000000

From 9323843583fcf1d885a06ce7ba9eadf3a1568c73 Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Tue, 7 Sep 2021 14:55:25 +0300
Subject: [PATCH 08/25] Derivatives who may not be correct yet, bounding box
 reader functions, etc

---
 pyCalculations/calculations.py |  1 +
 pyCalculations/derivatives.py  | 34 +++++++++++++++++++++++++---------
 pyVlsv/vlsvreader.py           | 27 +++++++++++++++++++++------
 3 files changed, 47 insertions(+), 15 deletions(-)

diff --git a/pyCalculations/calculations.py b/pyCalculations/calculations.py
index 483f8bc7..10a06270 100644
--- a/pyCalculations/calculations.py
+++ b/pyCalculations/calculations.py
@@ -56,4 +56,5 @@
 import fit
 from fieldtracer import static_field_tracer
 from fieldtracer import dynamic_field_tracer
+from derivatives import fg_divPoynting, fg_PoyntingFlux, fg_vol_curl
 
diff --git a/pyCalculations/derivatives.py b/pyCalculations/derivatives.py
index 3d91e6b7..d887267b 100644
--- a/pyCalculations/derivatives.py
+++ b/pyCalculations/derivatives.py
@@ -27,7 +27,7 @@
 import warnings
 from scipy import interpolate
 
-def get_fg_divPoynting:
+def fg_PoyntingFlux(bulkReader):
    b = bulkReader.read_fg_variable_as_volumetric('fg_b')
    print('b.shape=',b.shape)
    e = bulkReader.read_fg_variable_as_volumetric('fg_e')
@@ -35,18 +35,34 @@ def get_fg_divPoynting:
 
    mu_0 = 1.25663706144e-6
    S = np.cross(e,b)/mu_0
-   print('S.shape', S.shape)
+   return S
+
+
+def fg_divPoynting(bulkReader, dx=1e6):
+   S = fg_PoyntingFlux(bulkReader)
 
    #divS = pt.plot.plot_helpers.numdiv(pt.plot.plot_helpers.TransposeVectorArray(np.squeeze(S))).T
    divS = (np.roll(S[:,:,:,0],-1, 0) - np.roll(S[:,:,:,0], 1, 0) +
          np.roll(S[:,:,:,1],-1, 1) - np.roll(S[:,:,:,1], 1, 1) +
-         np.roll(S[:,:,:,2],-1, 2) - np.roll(S[:,:,:,2], 1, 2))/(2*1e6)
+         np.roll(S[:,:,:,2],-1, 2) - np.roll(S[:,:,:,2], 1, 2))/(2*dx)
    print('divS.shape', divS.shape)
    #divS = np.div
 
-   fig=plt.pyplot.figure()
-   ax=plt.pyplot.axes()
-   im=ax.imshow(divS[:,int(divS.shape[1]/2),:], origin='lower', norm=plt.colors.SymLogNorm(base=10, linthresh=1e-12, linscale = 1.0, vmin=-1e-7, vmax=1e-7, clip=True), cmap='vik')
-   plt.pyplot.colorbar(im, ax=ax)
-   plt.pyplot.savefig('./divS_centered.png')
-   print(divS)
\ No newline at end of file
+#   print(divS)
+   return divS
+
+def fg_vol_jacobian(array, dx=1e6):
+   dFx_dx, dFx_dy, dFx_dz = np.gradient(array[:,:,:,0], dx)
+   dFy_dx, dFy_dy, dFy_dz = np.gradient(array[:,:,:,1], dx)
+   dFz_dx, dFz_dy, dFz_dz = np.gradient(array[:,:,:,2], dx)
+
+def fg_vol_curl(array, dx=1e6):
+   dummy,  dFx_dy, dFx_dz = np.gradient(array[:,:,:,0], dx)
+   dFy_dx, dummy,  dFy_dz = np.gradient(array[:,:,:,1], dx)
+   dFz_dx, dFz_dy, dummy  = np.gradient(array[:,:,:,2], dx)
+
+   rotx = dFz_dy - dFy_dz
+   roty = dFx_dz - dFz_dx
+   rotz = dFy_dx - dFx_dy
+
+   return np.stack([rotx, roty, rotz], axis=-1)
\ No newline at end of file
diff --git a/pyVlsv/vlsvreader.py b/pyVlsv/vlsvreader.py
index 02784830..2880ca34 100644
--- a/pyVlsv/vlsvreader.py
+++ b/pyVlsv/vlsvreader.py
@@ -1082,8 +1082,8 @@ def read_fg_variable_as_volumetric(self, name, centering="default", operator="pa
       print('read in fgdata with shape', fgdata.shape, name)
       celldata = np.zeros_like(fgdata)
       known_centerings = {"fg_b":"face", "fg_e":"edge"}
-      if name.casefold() in known_centerings.keys() and centering == "default":
-         centering = known_centerings[name.casefold()]
+      if name.lower() in known_centerings.keys() and centering == "default":
+         centering = known_centerings[name.lower()]
       else:
          print("A variable with unknown centering! Aborting.")
          return False
@@ -1264,7 +1264,13 @@ def get_cell_fsgrid_slicemap(self, cellid):
       lowi, upi = self.get_fsgrid_slice_indices(low, up)
       return lowi, upi
 
-   def get_fsgrid_subarray(self, cellid, array):
+   def get_bbox_fsgrid_slicemap(self, low, up):
+      '''Returns a slice tuple of fsgrid indices that are contained in the (low, up) bounding box.
+      '''
+      lowi, upi = self.get_fsgrid_slice_indices(low, up)
+      return lowi, upi
+
+   def get_cell_fsgrid_subarray(self, cellid, array):
       '''Returns a subarray of the fsgrid array, corresponding to the fsgrid
       covered by the SpatialGrid cellid. [untested]
       '''
@@ -1275,11 +1281,21 @@ def get_fsgrid_subarray(self, cellid, array):
       else:
          return array[lowi[0]:upi[0]+1, lowi[1]:upi[1]+1, lowi[2]:upi[2]+1]
 
+   def get_bbox_fsgrid_subarray(self, low, up, array):
+      '''Returns a subarray of the fsgrid array, corresponding to the (low, up) bounding box.
+      '''
+      lowi, upi = self.get_bbox_fsgrid_slicemap(low,up)
+      #print('subarray:',lowi, upi)
+      if array.ndim == 4:
+         return array[lowi[0]:upi[0]+1, lowi[1]:upi[1]+1, lowi[2]:upi[2]+1, :]
+      else:
+         return array[lowi[0]:upi[0]+1, lowi[1]:upi[1]+1, lowi[2]:upi[2]+1]
+
 
    def downsample_fsgrid_subarray(self, cellid, array):
       '''Returns a mean value of fsgrid values underlying the SpatialGrid cellid.
       '''
-      fsarr = self.get_fsgrid_subarray(cellid, array)
+      fsarr = self.get_cell_fsgrid_subarray(cellid, array)
       n = fsarr.size
       if fsarr.ndim == 4:
          n = n/3
@@ -1290,8 +1306,7 @@ def downsample_fsgrid_subarray(self, cellid, array):
 
    def upsample_fsgrid_subarray(self, cellid, var, array):
       '''Set the elements of the fsgrid array to the value of corresponding SpatialGrid
-      cellid. Mutator.
-      [untested]
+      cellid. Mutator for array.
       '''
       lowi, upi = self.get_cell_fsgrid_slicemap(cellid)
       #print(lowi,upi)

From cc788a35147a4a20be63f71cc2910c953e402c6a Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Tue, 7 Sep 2021 17:23:45 +0300
Subject: [PATCH 09/25] pyVlsv.VlsvWriter updates: no empty mesh attributes,
 pretty printed footer as with Vlasiator output so that the vlsv VisIt plugin
 doesnt die.

---
 pyVlsv/vlsvwriter.py | 37 +++++++++++++++++++++++++++++++++----
 1 file changed, 33 insertions(+), 4 deletions(-)

diff --git a/pyVlsv/vlsvwriter.py b/pyVlsv/vlsvwriter.py
index e8eab964..9cacb609 100644
--- a/pyVlsv/vlsvwriter.py
+++ b/pyVlsv/vlsvwriter.py
@@ -69,6 +69,8 @@ def __initialize( self, vlsvReader ):
       tags['MESH_NODE_CRDS_Z'] = ''
       tags['MESH'] = ''
       tags['MESH_DOMAIN_SIZES'] = ''
+      tags['MESH_GHOST_DOMAINS'] = ''
+      tags['MESH_GHOST_LOCALIDS'] = ''
       tags['CellID'] = ''
       tags['MESH_BBOX'] = ''
       tags['COORDS'] = ''
@@ -79,7 +81,7 @@ def __initialize( self, vlsvReader ):
             if 'name' in child.attrib: name = child.attrib['name']
             else: name = ''
             if 'mesh' in child.attrib: mesh = child.attrib['mesh']
-            else: mesh = ''
+            else: mesh = None
             tag = child.tag
             extra_attribs = {}
             for i in child.attrib.items():
@@ -87,6 +89,7 @@ def __initialize( self, vlsvReader ):
                   extra_attribs[i[0]] = i[1]
             data = vlsvReader.read( name=name, tag=tag, mesh=mesh )
             # Write the data:
+            print('Writing tag', tag)
             self.write( data=data, name=name, tag=tag, mesh=mesh, extra_attribs=extra_attribs )
 
 
@@ -151,7 +154,10 @@ def copy_variables_list( self, vlsvReader, vars ):
             if 'mesh' in child.attrib:
                 mesh = child.attrib['mesh']
             else:
-                mesh = ''
+               if tag in ['VARIABLE']:
+                  print('MESH required')
+                  return
+               mesh = None
             tag = child.tag
             # Copy extra attributes:
             extra_attribs = {}
@@ -204,7 +210,8 @@ def write(self, data, name, tag, mesh, extra_attribs={}):
       # Add the data into the xml data:
       child = ET.SubElement(self.__xml_root, tag)
       child.attrib["name"] = name
-      child.attrib["mesh"] = mesh
+      if mesh is not None:
+         child.attrib["mesh"] = mesh
       child.attrib["arraysize"] = len(np.atleast_1d(data))
       if extra_attribs != '':
          for i in extra_attribs.items():
@@ -217,7 +224,14 @@ def write(self, data, name, tag, mesh, extra_attribs={}):
          return False
       else:
          child.attrib["vectorsize"] = 1
-         datatype = str(type(data[0]))
+         if(len(data) == 0):
+            if tag=="MESH_GHOST_DOMAINS" or tag=="MESH_GHOST_LOCALIDS":
+               datatype="int32"
+            else:
+               print("Trying to extract datatype from an empty array. I will fail as usual, since this is not the special case that is guarded against!")
+               datatype = str(type(data[0]))
+         else:
+            datatype = str(type(data[0]))
 
       # Parse the data types:
       if 'uint' in datatype:
@@ -271,6 +285,7 @@ def __write_xml_footer( self ):
          for i in child.attrib.items():
             child.attrib[i[0]] = str(child.attrib[i[0]])
       tree = ET.ElementTree( self.__xml_root)
+      xml_footer_indent( self.__xml_root)
       tree.write(fptr)
       # Write the offset (first 8 bytes = endianness):
       offset_endianness = 8
@@ -284,3 +299,17 @@ def close( self ):
       self.__write_xml_footer()
       self.__fptr.close()
 
+def xml_footer_indent(elem, level=0):
+   i = "\n" + level*"   "
+   if len(elem):
+      if not elem.text or not elem.text.strip():
+            elem.text = i + "   "
+      if not elem.tail or not elem.tail.strip():
+            elem.tail = i
+      for elem in elem:
+            xml_footer_indent(elem, level+1)
+      if not elem.tail or not elem.tail.strip():
+            elem.tail = i
+   else:
+      if level and (not elem.tail or not elem.tail.strip()):
+            elem.tail = i
\ No newline at end of file

From 9386f1d5f8db67f7231a1f3586e176e07284ec98 Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Wed, 8 Sep 2021 13:07:04 +0300
Subject: [PATCH 10/25] Added some of Jonas' helpers and ballooning criterion.

---
 pyCalculations/derivatives.py | 78 +++++++++++++++++++++++++++++++----
 1 file changed, 70 insertions(+), 8 deletions(-)

diff --git a/pyCalculations/derivatives.py b/pyCalculations/derivatives.py
index d887267b..b9aa14d2 100644
--- a/pyCalculations/derivatives.py
+++ b/pyCalculations/derivatives.py
@@ -41,20 +41,19 @@ def fg_PoyntingFlux(bulkReader):
 def fg_divPoynting(bulkReader, dx=1e6):
    S = fg_PoyntingFlux(bulkReader)
 
-   #divS = pt.plot.plot_helpers.numdiv(pt.plot.plot_helpers.TransposeVectorArray(np.squeeze(S))).T
    divS = (np.roll(S[:,:,:,0],-1, 0) - np.roll(S[:,:,:,0], 1, 0) +
          np.roll(S[:,:,:,1],-1, 1) - np.roll(S[:,:,:,1], 1, 1) +
          np.roll(S[:,:,:,2],-1, 2) - np.roll(S[:,:,:,2], 1, 2))/(2*dx)
    print('divS.shape', divS.shape)
-   #divS = np.div
 
-#   print(divS)
    return divS
 
-def fg_vol_jacobian(array, dx=1e6):
-   dFx_dx, dFx_dy, dFx_dz = np.gradient(array[:,:,:,0], dx)
-   dFy_dx, dFy_dy, dFy_dz = np.gradient(array[:,:,:,1], dx)
-   dFz_dx, dFz_dy, dFz_dz = np.gradient(array[:,:,:,2], dx)
+def fg_vol_jacobian(b, dx=1e6):
+   dFx_dx, dFx_dy, dFx_dz = np.gradient(b[:,:,:,0], dx)
+   dFy_dx, dFy_dy, dFy_dz = np.gradient(b[:,:,:,1], dx)
+   dFz_dx, dFz_dy, dFz_dz = np.gradient(b[:,:,:,2], dx)
+
+   return np.stack()
 
 def fg_vol_curl(array, dx=1e6):
    dummy,  dFx_dy, dFx_dz = np.gradient(array[:,:,:,0], dx)
@@ -65,4 +64,67 @@ def fg_vol_curl(array, dx=1e6):
    roty = dFx_dz - dFz_dx
    rotz = dFy_dx - dFx_dy
 
-   return np.stack([rotx, roty, rotz], axis=-1)
\ No newline at end of file
+   return np.stack([rotx, roty, rotz], axis=-1)
+
+def fg_vol_div(array, dx=1e6):
+   dFx_dx, dummy, dummy = np.gradient(array[:,:,:,0], dx)
+   dummy, dFy_dy, dummy = np.gradient(array[:,:,:,1], dx)
+   dummy, dummy, dFz_dz = np.gradient(array[:,:,:,2], dx)
+
+   return dFx_dx+dFy_dy+dFz_dz
+
+def vfield3_dot(a, b):
+    """Calculates dot product of vectors a and b in 3D vector field"""
+
+    return (
+        a[:, :, :, 0] * b[:, :, :, 0]
+        + a[:, :, :, 1] * b[:, :, :, 1]
+        + a[:, :, :, 2] * b[:, :, :, 2]
+    )
+
+def vfield3_matder(a, b, dr):
+    """Calculates material derivative of 3D vector fields a and b"""
+
+    bx = b[:, :, :, 0]
+    by = b[:, :, :, 1]
+    bz = b[:, :, :, 2]
+
+    grad_bx = np.gradient(bx, dr)
+    grad_by = np.gradient(by, dr)
+    grad_bz = np.gradient(bz, dr)
+
+    resx = vfield3_dot(a, grad_bx)
+    resy = vfield3_dot(a, grad_by)
+    resz = vfield3_dot(a, grad_bz)
+
+    return np.stack((resx, resy, resz), axis=-1)
+
+def vfield3_normalise(a):
+
+    amag = np.linalg.norm(a, axis=-1)
+
+    resx = a[:, :, :, 0] / amag
+    resy = a[:, :, :, 1] / amag
+    resz = a[:, :, :, 2] / amag
+
+    return np.stack((resx, resy, resz), axis=-1)
+
+def ballooning_crit(B, P, beta, dr=1000e3):
+
+    # Bmag = np.linalg.norm(B, axis=-1)
+
+    b = vfield3_normalise(B)
+
+    n = vfield3_matder(b, b, dr)
+
+    nnorm = vfield3_normalise(n)
+
+    gradP = np.gradient(P,dr)
+
+    kappaP = vfield3_dot(nnorm, gradP) / P
+    
+    kappaC = vfield3_dot(nnorm, n)
+
+    balloon = (2 + beta) / 4.0 * kappaP / (kappaC + 1e-27)
+
+    return (balloon, nnorm, kappaC, gradP)
\ No newline at end of file

From c927cae8a195534d0cd21175ac52a62cbe3e9b28 Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Thu, 9 Sep 2021 12:40:43 +0300
Subject: [PATCH 11/25] a merge

---
 pyCalculations/calculations.py |  3 +--
 pyCalculations/derivatives.py  | 33 ++++++++++++++++++---------------
 2 files changed, 19 insertions(+), 17 deletions(-)

diff --git a/pyCalculations/calculations.py b/pyCalculations/calculations.py
index 10a06270..c7bb589c 100644
--- a/pyCalculations/calculations.py
+++ b/pyCalculations/calculations.py
@@ -56,5 +56,4 @@
 import fit
 from fieldtracer import static_field_tracer
 from fieldtracer import dynamic_field_tracer
-from derivatives import fg_divPoynting, fg_PoyntingFlux, fg_vol_curl
-
+from derivatives import fg_divPoynting, fg_PoyntingFlux, fg_vol_curl, ballooning_crit,vfield3_curvature
diff --git a/pyCalculations/derivatives.py b/pyCalculations/derivatives.py
index b9aa14d2..a401485a 100644
--- a/pyCalculations/derivatives.py
+++ b/pyCalculations/derivatives.py
@@ -76,11 +76,12 @@ def fg_vol_div(array, dx=1e6):
 def vfield3_dot(a, b):
     """Calculates dot product of vectors a and b in 3D vector field"""
 
-    return (
-        a[:, :, :, 0] * b[:, :, :, 0]
-        + a[:, :, :, 1] * b[:, :, :, 1]
-        + a[:, :, :, 2] * b[:, :, :, 2]
-    )
+    return np.sum(np.multiply(a,b), axis=-1)
+    #return (
+    #    a[:, :, :, 0] * b[:, :, :, 0]
+    #    + a[:, :, :, 1] * b[:, :, :, 1]
+    #    + a[:, :, :, 2] * b[:, :, :, 2]
+    #)
 
 def vfield3_matder(a, b, dr):
     """Calculates material derivative of 3D vector fields a and b"""
@@ -103,19 +104,21 @@ def vfield3_normalise(a):
 
     amag = np.linalg.norm(a, axis=-1)
 
-    resx = a[:, :, :, 0] / amag
-    resy = a[:, :, :, 1] / amag
-    resz = a[:, :, :, 2] / amag
+    res=a/amag
+    return res
+    #resx = a[:, :, :, 0] / amag
+    #resy = a[:, :, :, 1] / amag
+    #resz = a[:, :, :, 2] / amag
 
-    return np.stack((resx, resy, resz), axis=-1)
-
-def ballooning_crit(B, P, beta, dr=1000e3):
+    #return np.stack((resx, resy, resz), axis=-1)
 
-    # Bmag = np.linalg.norm(B, axis=-1)
+def vfield3_curvature(a, dr=1000e3):
+   a = vfield3_normalise(a)
+   return vfield3_matder(a, a, dr)
 
-    b = vfield3_normalise(B)
+def ballooning_crit(B, P, beta, dr=1000e3):
 
-    n = vfield3_matder(b, b, dr)
+    n = vfield3_curvature(B, dr)
 
     nnorm = vfield3_normalise(n)
 
@@ -127,4 +130,4 @@ def ballooning_crit(B, P, beta, dr=1000e3):
 
     balloon = (2 + beta) / 4.0 * kappaP / (kappaC + 1e-27)
 
-    return (balloon, nnorm, kappaC, gradP)
\ No newline at end of file
+    return (balloon, kappaC, gradP)
\ No newline at end of file

From 345ccbfc3722545e3cdcb8ea467aae3e2070242b Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Thu, 23 Sep 2021 09:10:12 +0300
Subject: [PATCH 12/25] curvature exposed, streamlined

---
 pyCalculations/calculations.py |  3 +--
 pyCalculations/derivatives.py  | 33 ++++++++++++++++++---------------
 2 files changed, 19 insertions(+), 17 deletions(-)

diff --git a/pyCalculations/calculations.py b/pyCalculations/calculations.py
index 10a06270..c7bb589c 100644
--- a/pyCalculations/calculations.py
+++ b/pyCalculations/calculations.py
@@ -56,5 +56,4 @@
 import fit
 from fieldtracer import static_field_tracer
 from fieldtracer import dynamic_field_tracer
-from derivatives import fg_divPoynting, fg_PoyntingFlux, fg_vol_curl
-
+from derivatives import fg_divPoynting, fg_PoyntingFlux, fg_vol_curl, ballooning_crit,vfield3_curvature
diff --git a/pyCalculations/derivatives.py b/pyCalculations/derivatives.py
index b9aa14d2..a401485a 100644
--- a/pyCalculations/derivatives.py
+++ b/pyCalculations/derivatives.py
@@ -76,11 +76,12 @@ def fg_vol_div(array, dx=1e6):
 def vfield3_dot(a, b):
     """Calculates dot product of vectors a and b in 3D vector field"""
 
-    return (
-        a[:, :, :, 0] * b[:, :, :, 0]
-        + a[:, :, :, 1] * b[:, :, :, 1]
-        + a[:, :, :, 2] * b[:, :, :, 2]
-    )
+    return np.sum(np.multiply(a,b), axis=-1)
+    #return (
+    #    a[:, :, :, 0] * b[:, :, :, 0]
+    #    + a[:, :, :, 1] * b[:, :, :, 1]
+    #    + a[:, :, :, 2] * b[:, :, :, 2]
+    #)
 
 def vfield3_matder(a, b, dr):
     """Calculates material derivative of 3D vector fields a and b"""
@@ -103,19 +104,21 @@ def vfield3_normalise(a):
 
     amag = np.linalg.norm(a, axis=-1)
 
-    resx = a[:, :, :, 0] / amag
-    resy = a[:, :, :, 1] / amag
-    resz = a[:, :, :, 2] / amag
+    res=a/amag
+    return res
+    #resx = a[:, :, :, 0] / amag
+    #resy = a[:, :, :, 1] / amag
+    #resz = a[:, :, :, 2] / amag
 
-    return np.stack((resx, resy, resz), axis=-1)
-
-def ballooning_crit(B, P, beta, dr=1000e3):
+    #return np.stack((resx, resy, resz), axis=-1)
 
-    # Bmag = np.linalg.norm(B, axis=-1)
+def vfield3_curvature(a, dr=1000e3):
+   a = vfield3_normalise(a)
+   return vfield3_matder(a, a, dr)
 
-    b = vfield3_normalise(B)
+def ballooning_crit(B, P, beta, dr=1000e3):
 
-    n = vfield3_matder(b, b, dr)
+    n = vfield3_curvature(B, dr)
 
     nnorm = vfield3_normalise(n)
 
@@ -127,4 +130,4 @@ def ballooning_crit(B, P, beta, dr=1000e3):
 
     balloon = (2 + beta) / 4.0 * kappaP / (kappaC + 1e-27)
 
-    return (balloon, nnorm, kappaC, gradP)
\ No newline at end of file
+    return (balloon, kappaC, gradP)
\ No newline at end of file

From 049e6e46b2706616ca32437f8e770a3acf8a5e65 Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Tue, 28 Sep 2021 12:26:10 +0300
Subject: [PATCH 13/25] Requested changes

---
 pyCalculations/derivatives.py | 73 +++++++++++++++++++++--------------
 pyVlsv/vlsvreader.py          | 36 +++++++++--------
 pyVlsv/vlsvwriter.py          |  1 -
 3 files changed, 62 insertions(+), 48 deletions(-)

diff --git a/pyCalculations/derivatives.py b/pyCalculations/derivatives.py
index a401485a..531cfdff 100644
--- a/pyCalculations/derivatives.py
+++ b/pyCalculations/derivatives.py
@@ -27,10 +27,10 @@
 import warnings
 from scipy import interpolate
 
-def fg_PoyntingFlux(bulkReader):
-   b = bulkReader.read_fg_variable_as_volumetric('fg_b')
+def fg_PoyntingFlux(reader):
+   b = reader.read_fg_variable_as_volumetric('fg_b')
    print('b.shape=',b.shape)
-   e = bulkReader.read_fg_variable_as_volumetric('fg_e')
+   e = reader.read_fg_variable_as_volumetric('fg_e')
    print('e.shape=',e.shape)
 
    mu_0 = 1.25663706144e-6
@@ -38,24 +38,35 @@ def fg_PoyntingFlux(bulkReader):
    return S
 
 
-def fg_divPoynting(bulkReader, dx=1e6):
-   S = fg_PoyntingFlux(bulkReader)
-
-   divS = (np.roll(S[:,:,:,0],-1, 0) - np.roll(S[:,:,:,0], 1, 0) +
-         np.roll(S[:,:,:,1],-1, 1) - np.roll(S[:,:,:,1], 1, 1) +
-         np.roll(S[:,:,:,2],-1, 2) - np.roll(S[:,:,:,2], 1, 2))/(2*dx)
+def fg_divPoynting(reader):
+   S = fg_PoyntingFlux(reader)
+   dx = reader.get_fsgrid_cell_size()
+   divS = ((np.roll(S[:,:,:,0],-1, 0) - np.roll(S[:,:,:,0], 1, 0))/(2*dx[0]) +
+           (np.roll(S[:,:,:,1],-1, 1) - np.roll(S[:,:,:,1], 1, 1))/(2*dx[1]) +
+           (np.roll(S[:,:,:,2],-1, 2) - np.roll(S[:,:,:,2], 1, 2))/(2*dx[2])
+          )
    print('divS.shape', divS.shape)
 
    return divS
 
-def fg_vol_jacobian(b, dx=1e6):
-   dFx_dx, dFx_dy, dFx_dz = np.gradient(b[:,:,:,0], dx)
-   dFy_dx, dFy_dy, dFy_dz = np.gradient(b[:,:,:,1], dx)
-   dFz_dx, dFz_dy, dFz_dz = np.gradient(b[:,:,:,2], dx)
-
-   return np.stack()
-
-def fg_vol_curl(array, dx=1e6):
+def fg_vol_jacobian(reader, b):
+   # Return the jacobian of an fsgrid volumetric variable
+   # result is a 3x3 array for each cell
+   # Last index of output gives the vector component,
+   # second-to-last the direction of derivative
+   dx = reader.get_fsgrid_cell_size()
+   dFx_dx, dFx_dy, dFx_dz = np.gradient(b[:,:,:,0], dx[0])
+   dFy_dx, dFy_dy, dFy_dz = np.gradient(b[:,:,:,1], dx[1])
+   dFz_dx, dFz_dy, dFz_dz = np.gradient(b[:,:,:,2], dx[2])
+
+   dFx = np.stack([dFx_dx, dFx_dy, dFx_dz], axis=-1)
+   dFy = np.stack([dFy_dx, dFy_dy, dFy_dz], axis=-1)
+   dFz = np.stack([dFz_dx, dFz_dy, dFz_dz], axis=-1)
+   
+   return np.stack([dFx,dFy,dFz], axis=-1)
+
+def fg_vol_curl(reader, array):
+   dx = reader.get_fsgrid_cell_size()
    dummy,  dFx_dy, dFx_dz = np.gradient(array[:,:,:,0], dx)
    dFy_dx, dummy,  dFy_dz = np.gradient(array[:,:,:,1], dx)
    dFz_dx, dFz_dy, dummy  = np.gradient(array[:,:,:,2], dx)
@@ -66,33 +77,34 @@ def fg_vol_curl(array, dx=1e6):
 
    return np.stack([rotx, roty, rotz], axis=-1)
 
-def fg_vol_div(array, dx=1e6):
-   dFx_dx, dummy, dummy = np.gradient(array[:,:,:,0], dx)
-   dummy, dFy_dy, dummy = np.gradient(array[:,:,:,1], dx)
-   dummy, dummy, dFz_dz = np.gradient(array[:,:,:,2], dx)
+def fg_vol_div(reader, array):
+   dx = reader.get_fsgrid_cell_size()
+   dFx_dx = np.gradient(array[:,:,:,0], dx[0], axis=0)
+   dFy_dy = np.gradient(array[:,:,:,1], dx[1], axis=1)
+   dFz_dz = np.gradient(array[:,:,:,2], dx[2], axis=2)
 
    return dFx_dx+dFy_dy+dFz_dz
 
 def vfield3_dot(a, b):
     """Calculates dot product of vectors a and b in 3D vector field"""
 
-    return np.sum(np.multiply(a,b), axis=-1)
+    return (a*b).sum(-1)
     #return (
     #    a[:, :, :, 0] * b[:, :, :, 0]
     #    + a[:, :, :, 1] * b[:, :, :, 1]
     #    + a[:, :, :, 2] * b[:, :, :, 2]
     #)
 
-def vfield3_matder(a, b, dr):
+def vfield3_matder(reader, a, b):
     """Calculates material derivative of 3D vector fields a and b"""
-
+    dr = reader.get_fsgrid_cell_size()
     bx = b[:, :, :, 0]
     by = b[:, :, :, 1]
     bz = b[:, :, :, 2]
 
-    grad_bx = np.gradient(bx, dr)
-    grad_by = np.gradient(by, dr)
-    grad_bz = np.gradient(bz, dr)
+    grad_bx = np.gradient(bx, dr[0])
+    grad_by = np.gradient(by, dr[1])
+    grad_bz = np.gradient(bz, dr[2])
 
     resx = vfield3_dot(a, grad_bx)
     resy = vfield3_dot(a, grad_by)
@@ -112,12 +124,13 @@ def vfield3_normalise(a):
 
     #return np.stack((resx, resy, resz), axis=-1)
 
-def vfield3_curvature(a, dr=1000e3):
+def vfield3_curvature(reader, a):
+   dr = reader.get_fsgrid_cell_size()
    a = vfield3_normalise(a)
    return vfield3_matder(a, a, dr)
 
-def ballooning_crit(B, P, beta, dr=1000e3):
-
+def ballooning_crit(reader, B, P, beta):
+    dr = reader.get_fsgrid_cell_size()
     n = vfield3_curvature(B, dr)
 
     nnorm = vfield3_normalise(n)
diff --git a/pyVlsv/vlsvreader.py b/pyVlsv/vlsvreader.py
index 2880ca34..20ce7950 100644
--- a/pyVlsv/vlsvreader.py
+++ b/pyVlsv/vlsvreader.py
@@ -1078,16 +1078,18 @@ def read_fg_variable_as_volumetric(self, name, centering="default", operator="pa
 
       fssize=list(self.get_fsgrid_mesh_size())
       if 1 in fssize:
-         fgdata=np.expand_dims(fgdata, fssize.index(1)) #expand to have a singleton dimension for a reduced dim - lets roll happen with ease
-      print('read in fgdata with shape', fgdata.shape, name)
+         #expand to have a singleton dimension for a reduced dim - lets roll happen with ease
+         singletons = [i for i, sz in enumerate(fssize) if sz == 1]
+         for dim in singletons:
+            fgdata=np.expand_dims(fgdata, dim)
+      #print('read in fgdata with shape', fgdata.shape, name)
       celldata = np.zeros_like(fgdata)
       known_centerings = {"fg_b":"face", "fg_e":"edge"}
       if name.lower() in known_centerings.keys() and centering == "default":
          centering = known_centerings[name.lower()]
       else:
-         print("A variable with unknown centering! Aborting.")
+         print("A variable ("+name+") with unknown centering! Aborting.")
          return False
-      print('fgdata.shape', fgdata.shape)
       #vector variable
       if fgdata.shape[-1] == 3:
          if centering=="face":
@@ -1346,10 +1348,7 @@ def get_fsgrid_coordinates(self, ri):
       '''Returns real-space center coordinates of the fsgrid 3-index.
       '''
       lowerlimit = self.get_fsgrid_mesh_extent()[0:3]
-      upperlimit = self.get_fsgrid_mesh_extent()[3:6]
-      delta = upperlimit-lowerlimit
-      sz=self.get_fsgrid_mesh_size()
-      dxs = delta/sz
+      dxs = self.get_fsgrid_cell_size()
 
       return lowerlimit+dxs*(np.array(ri)+0.5)
 
@@ -1950,6 +1949,16 @@ def get_fsgrid_mesh_extent(self):
       '''
       return np.array([self.__xmin, self.__ymin, self.__zmin, self.__xmax, self.__ymax, self.__zmax])
 
+   def get_fsgrid_cell_size(self):
+      ''' Read fsgrid cell size
+      
+      :returns: Maximum and minimum coordinates of the mesh, [dx, dy, dz]
+      '''
+      size = self.get_fsgrid_mesh_size()
+      ext = self.get_fsgrid_mesh_extent()
+      ext = ext[3:6]-ext[0:3]
+      return ext/size
+
    def get_fsgrid_indices(self, coords):
       ''' Convert spatial coordinates coords to an index array [xi, yi, zi] for fsgrid
 
@@ -1959,10 +1968,7 @@ def get_fsgrid_indices(self, coords):
       fsvar_at_coords = fsvar_array.item(ii)
       '''
       lower = self.get_fsgrid_mesh_extent()[0:3]
-      upper = self.get_fsgrid_mesh_extent()[3:6]
-      delta = upper-lower
-      sz=self.get_fsgrid_mesh_size()
-      dx = delta/sz
+      dx = self.get_fsgrid_cell_size()
       r0 = coords-lower
       ri = np.floor(r0/dx).astype(int)
       if (ri < 0).any() or (ri>sz-1).any():
@@ -1982,11 +1988,7 @@ def get_fsgrid_slice_indices(self, lower, upper, eps=1e-3):
       ii = f.get_fsgrid_mesh_extent(coords)
       fsvar_at_coords = fsvar_array.item(ii)
       '''
-      lowerlimit = self.get_fsgrid_mesh_extent()[0:3]
-      upperlimit = self.get_fsgrid_mesh_extent()[3:6]
-      delta = upperlimit-lowerlimit
-      sz=self.get_fsgrid_mesh_size()
-      dx = delta/sz
+      dx = self.get_fsgrid_cell_size()
       eps = dx*eps
       loweri = self.get_fsgrid_indices(lower+eps)
       upperi = self.get_fsgrid_indices(upper-eps)
diff --git a/pyVlsv/vlsvwriter.py b/pyVlsv/vlsvwriter.py
index 9cacb609..40379a4f 100644
--- a/pyVlsv/vlsvwriter.py
+++ b/pyVlsv/vlsvwriter.py
@@ -89,7 +89,6 @@ def __initialize( self, vlsvReader ):
                   extra_attribs[i[0]] = i[1]
             data = vlsvReader.read( name=name, tag=tag, mesh=mesh )
             # Write the data:
-            print('Writing tag', tag)
             self.write( data=data, name=name, tag=tag, mesh=mesh, extra_attribs=extra_attribs )
 
 

From 31c21d70a1261cf89a7cb4716e45c9b3283893e5 Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Tue, 28 Sep 2021 12:31:26 +0300
Subject: [PATCH 14/25] Requested changes 2: errant plot_helpers.PLANE change

---
 pyPlots/plot_helpers.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyPlots/plot_helpers.py b/pyPlots/plot_helpers.py
index c092e82a..39edf51c 100644
--- a/pyPlots/plot_helpers.py
+++ b/pyPlots/plot_helpers.py
@@ -26,7 +26,7 @@
 import sys
 from rotation import rotateTensorToVector
 
-PLANE = 'XZ'
+PLANE = 'XY'
 # or alternatively, 'XZ'
 CELLSIZE = 300000.0 # cell size
 DT = 0.5 # time step

From d7979c0658adfd5d65bfa55f18a918cf77a9eeac Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Tue, 28 Sep 2021 13:30:30 +0300
Subject: [PATCH 15/25] Default centering to None, allowed manual centering

---
 pyVlsv/vlsvreader.py | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/pyVlsv/vlsvreader.py b/pyVlsv/vlsvreader.py
index 20ce7950..ae7506ae 100644
--- a/pyVlsv/vlsvreader.py
+++ b/pyVlsv/vlsvreader.py
@@ -1073,7 +1073,7 @@ def calcLocalSize(globalCells, ntasks, my_n):
 
        return np.squeeze(orderedData)
 
-   def read_fg_variable_as_volumetric(self, name, centering="default", operator="pass"):
+   def read_fg_variable_as_volumetric(self, name, centering=None, operator="pass"):
       fgdata = self.read_fsgrid_variable(name, operator)
 
       fssize=list(self.get_fsgrid_mesh_size())
@@ -1085,11 +1085,13 @@ def read_fg_variable_as_volumetric(self, name, centering="default", operator="pa
       #print('read in fgdata with shape', fgdata.shape, name)
       celldata = np.zeros_like(fgdata)
       known_centerings = {"fg_b":"face", "fg_e":"edge"}
-      if name.lower() in known_centerings.keys() and centering == "default":
-         centering = known_centerings[name.lower()]
-      else:
-         print("A variable ("+name+") with unknown centering! Aborting.")
-         return False
+      if centering is None:
+         try:
+            centering = known_centerings[name.lower()]
+         except KeyError:
+            print("A variable ("+name+") with unknown centering! Aborting.")
+            return False
+         
       #vector variable
       if fgdata.shape[-1] == 3:
          if centering=="face":

From 9f64b6ff3403f91de78687c2bf593d83893e70d0 Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Tue, 28 Sep 2021 15:49:20 +0300
Subject: [PATCH 16/25] Last missing requests

---
 pyVlsv/vlsvreader.py | 22 +++++++++++++++++++++-
 1 file changed, 21 insertions(+), 1 deletion(-)

diff --git a/pyVlsv/vlsvreader.py b/pyVlsv/vlsvreader.py
index ae7506ae..6cf84afa 100644
--- a/pyVlsv/vlsvreader.py
+++ b/pyVlsv/vlsvreader.py
@@ -975,6 +975,22 @@ def read_interpolated_variable(self, name, coordinates, operator="pass",periodic
          # Done.
          return ret_array
 
+   def read_fsgrid_variable_cellid(self, name, cellids=-1, operator="pass"):
+      ''' Reads fsgrid variables from the open vlsv file.
+       Arguments:
+       :param name:     Name of the variable
+       :param cellids:  SpatialGrid cellids for which to fetch data. Default: return full data
+       :param operator: Datareduction operator. "pass" does no operation on data
+       :returns: *ordered* list of numpy arrays with the data
+
+       ... seealso:: :func:`read_fsgrid_variable`
+       '''
+      var = self.read_fsgrid_variable(name, operator=operator)
+      if cellids == -1:
+         return var
+      else:
+         return [self.downsample_fsgrid_subarray(cid, var) for cid in cellids]
+
    def read_fsgrid_variable(self, name, operator="pass"):
        ''' Reads fsgrid variables from the open vlsv file.
        Arguments:
@@ -1276,7 +1292,7 @@ def get_bbox_fsgrid_slicemap(self, low, up):
 
    def get_cell_fsgrid_subarray(self, cellid, array):
       '''Returns a subarray of the fsgrid array, corresponding to the fsgrid
-      covered by the SpatialGrid cellid. [untested]
+      covered by the SpatialGrid cellid.
       '''
       lowi, upi = self.get_cell_fsgrid_slicemap(cellid)
       #print('subarray:',lowi, upi)
@@ -1308,6 +1324,10 @@ def downsample_fsgrid_subarray(self, cellid, array):
          print("Warning: weird fs subarray size", n, 'for amrlevel', self.get_amr_level(cellid), 'expect', ncells)
       return np.mean(fsarr,axis=(0,1,2))
 
+   def get_cellid_at_fsgrid_index(i,j,k):
+      coords = self.get_fsgrid_coordinates([i,j,k])
+      return self.get_cellid(coord)
+
    def upsample_fsgrid_subarray(self, cellid, var, array):
       '''Set the elements of the fsgrid array to the value of corresponding SpatialGrid
       cellid. Mutator for array.

From 9b957faaf38f406a19c0297e905a44419413c5c7 Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Wed, 29 Sep 2021 16:34:56 +0300
Subject: [PATCH 17/25] A vlsvwriter.write wrapper to explicitly require
 necessary variable metadata.

---
 pyVlsv/vlsvwriter.py | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/pyVlsv/vlsvwriter.py b/pyVlsv/vlsvwriter.py
index 40379a4f..3e252d70 100644
--- a/pyVlsv/vlsvwriter.py
+++ b/pyVlsv/vlsvwriter.py
@@ -260,6 +260,24 @@ def write(self, data, name, tag, mesh, extra_attribs={}):
       # write the xml footer:
       self.__write_xml_footer()
 
+   def write_variable(self, data, name, mesh, variableLaTex, unit, unitLaTeX, unitConversion, extra_attribs={}):
+      ''' Writes an array into the vlsv file as a variable; requires input of metadata required by VlsvReader
+
+      :param name: Name of the data array
+      :param mesh: Mesh for the data array
+      :param variableLaTeX: LaTeX string representation of variable
+      :param unit: plaintext string represenation of the unit
+      :param unitLaTeX: LaTeX string represenation of the unit
+      :param unitConversion: string represenation of the unit conversion to get to SI
+      :param extra_attribs: Dictionary with whatever xml attributes that should be defined in the array that aren't name, tag, or mesh
+
+      :returns: True if the data was written successfully
+
+      '''
+
+      return self.write(data, name, 'VARIABLE', mesh, extra_attribs={'variableLaTeX':variableLaTex, 'unit':unit, 'unitLaTeX':unitLaTeX, 'unitConversion':unitConversion}.update(extra_attribs))
+
+
    def write_fgarray_to_SpatialGrid(self, reader, data, name, extra_attribs={}):
       # get a reader for the target file
       #print(data.shape[0:3], reader.get_fsgrid_mesh_size(), (data.shape[0:3] == reader.get_fsgrid_mesh_size()))

From c9b026320f01aefa74c1dbd7150dbf09497f63bf Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Wed, 26 Jan 2022 11:20:07 +0200
Subject: [PATCH 18/25] Utilities

---
 pyCalculations/derivatives.py |  6 +++---
 pyPlots/plot_colormap.py      |  6 +++++-
 pyVlsv/vlsvreader.py          | 25 ++++++++++++++++++++++---
 pyVlsv/vlsvwriter.py          |  4 ++--
 4 files changed, 32 insertions(+), 9 deletions(-)

diff --git a/pyCalculations/derivatives.py b/pyCalculations/derivatives.py
index 531cfdff..d80a6720 100644
--- a/pyCalculations/derivatives.py
+++ b/pyCalculations/derivatives.py
@@ -67,9 +67,9 @@ def fg_vol_jacobian(reader, b):
 
 def fg_vol_curl(reader, array):
    dx = reader.get_fsgrid_cell_size()
-   dummy,  dFx_dy, dFx_dz = np.gradient(array[:,:,:,0], dx)
-   dFy_dx, dummy,  dFy_dz = np.gradient(array[:,:,:,1], dx)
-   dFz_dx, dFz_dy, dummy  = np.gradient(array[:,:,:,2], dx)
+   dummy,  dFx_dy, dFx_dz = np.gradient(array[:,:,:,0], *dx)
+   dFy_dx, dummy,  dFy_dz = np.gradient(array[:,:,:,1], *dx)
+   dFz_dx, dFz_dy, dummy  = np.gradient(array[:,:,:,2], *dx)
 
    rotx = dFz_dy - dFy_dz
    roty = dFx_dz - dFz_dx
diff --git a/pyPlots/plot_colormap.py b/pyPlots/plot_colormap.py
index b705e917..8f92e936 100644
--- a/pyPlots/plot_colormap.py
+++ b/pyPlots/plot_colormap.py
@@ -249,7 +249,11 @@ def exprMA_cust(exprmaps, requestvariables=False):
     if filename:
         f=pt.vlsvfile.VlsvReader(filename)
     elif (filedir and step is not None):
-        filename = glob.glob(filedir+'bulk*'+str(step).rjust(7,'0')+'.vlsv')[0]
+        try:
+            filename = glob.glob(filedir+'bulk*'+str(step).rjust(7,'0')+'.vlsv')[0]
+        except IndexError as err:
+            print("IndexError for filedir =", filedir, ", step = ", step)
+            return
         #filename = filedir+'bulk.'+str(step).rjust(7,'0')+'.vlsv'
         f=pt.vlsvfile.VlsvReader(filename)
     elif vlsvobj:
diff --git a/pyVlsv/vlsvreader.py b/pyVlsv/vlsvreader.py
index 6cf84afa..dda3d4c6 100644
--- a/pyVlsv/vlsvreader.py
+++ b/pyVlsv/vlsvreader.py
@@ -979,7 +979,7 @@ def read_fsgrid_variable_cellid(self, name, cellids=-1, operator="pass"):
       ''' Reads fsgrid variables from the open vlsv file.
        Arguments:
        :param name:     Name of the variable
-       :param cellids:  SpatialGrid cellids for which to fetch data. Default: return full data
+       :param cellids:  SpatialGrid cellids for which to fetch data. Default: return full fsgrid data
        :param operator: Datareduction operator. "pass" does no operation on data
        :returns: *ordered* list of numpy arrays with the data
 
@@ -1324,9 +1324,27 @@ def downsample_fsgrid_subarray(self, cellid, array):
          print("Warning: weird fs subarray size", n, 'for amrlevel', self.get_amr_level(cellid), 'expect', ncells)
       return np.mean(fsarr,axis=(0,1,2))
 
-   def get_cellid_at_fsgrid_index(i,j,k):
+   def fsgrid_array_to_vg(self, array):
+      cellIds=self.read_variable("CellID")
+
+      #if not hasattr(self, 'fsCellIdTargets'):
+      #   self.fsCellIdTargets = np.zeros(self.get_fsgrid_mesh_size())
+      #   for cid in cellIds:
+      #      lowi, upi = self.get_cell_fsgrid_slicemap(cid)
+      #      self.fsCellIdTargets[lowi[0]:upi[0]+1, lowi[1]:upi[1]+1, lowi[2]:upi[2]+1] = cid
+      
+      vgarr = [self.downsample_fsgrid_subarray(cid, array) for cid in cellIds]
+      return vgarr
+
+   def vg_uniform_grid_process(self, variable, expr, exprtuple):
+      cellIds=self.read_variable("CellID")
+      array = self.read_variable_as_fg(variable)
+      array = expr(*exprtuple)
+      return self.fsgrid_array_to_vg(array)
+
+   def get_cellid_at_fsgrid_index(self, i,j,k):
       coords = self.get_fsgrid_coordinates([i,j,k])
-      return self.get_cellid(coord)
+      return self.get_cellid(coords)
 
    def upsample_fsgrid_subarray(self, cellid, var, array):
       '''Set the elements of the fsgrid array to the value of corresponding SpatialGrid
@@ -1993,6 +2011,7 @@ def get_fsgrid_indices(self, coords):
       dx = self.get_fsgrid_cell_size()
       r0 = coords-lower
       ri = np.floor(r0/dx).astype(int)
+      sz = self.get_fsgrid_mesh_size()
       if (ri < 0).any() or (ri>sz-1).any():
          print("get_fsgrid_indices: Resulting index out of bounds, returning None")
          return None
diff --git a/pyVlsv/vlsvwriter.py b/pyVlsv/vlsvwriter.py
index 3e252d70..9f18d9ed 100644
--- a/pyVlsv/vlsvwriter.py
+++ b/pyVlsv/vlsvwriter.py
@@ -240,7 +240,7 @@ def write(self, data, name, tag, mesh, extra_attribs={}):
       elif 'float' in datatype:
          child.attrib["datatype"] = "float"
       else:
-         print("BAD DATATYPE")
+         print("BAD DATATYPE: " + datatype)
          return False
 
       if '64' in datatype:
@@ -275,7 +275,7 @@ def write_variable(self, data, name, mesh, variableLaTex, unit, unitLaTeX, unitC
 
       '''
 
-      return self.write(data, name, 'VARIABLE', mesh, extra_attribs={'variableLaTeX':variableLaTex, 'unit':unit, 'unitLaTeX':unitLaTeX, 'unitConversion':unitConversion}.update(extra_attribs))
+      return self.write(data, name, 'VARIABLE', mesh, extra_attribs={'variableLaTeX':variableLaTex, 'unit':unit, 'unitLaTeX':unitLaTeX, 'unitConversion':unitConversion})
 
 
    def write_fgarray_to_SpatialGrid(self, reader, data, name, extra_attribs={}):

From 38b5f173d622dfc7badbc96f8acd0c8360808861 Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Thu, 19 May 2022 13:18:06 +0300
Subject: [PATCH 19/25] Guarding against missing paths

---
 pyVlsv/vlsvreader.py | 8 ++++++--
 pyVlsv/vlsvwriter.py | 7 +++++--
 2 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/pyVlsv/vlsvreader.py b/pyVlsv/vlsvreader.py
index dda3d4c6..890f9e38 100644
--- a/pyVlsv/vlsvreader.py
+++ b/pyVlsv/vlsvreader.py
@@ -25,7 +25,7 @@
 import xml.etree.ElementTree as ET
 import ast
 import numpy as np
-import os
+import os,sys
 import numbers
 import vlsvvariables
 from reduction import datareducers,multipopdatareducers,data_operators,v5reducers,multipopv5reducers
@@ -54,7 +54,11 @@ def __init__(self, file_name):
       file_name = os.path.abspath(file_name)
 
       self.file_name = file_name
-      self.__fptr = open(self.file_name,"rb")
+      try:
+         self.__fptr = open(self.file_name,"rb")
+      except FileNotFoundError as e:
+         print("File not found: ", self.file_name)
+         raise e
       self.__xml_root = ET.fromstring("<VLSV></VLSV>")
       self.__fileindex_for_cellid={}
       self.__max_spatial_amr_level = -1
diff --git a/pyVlsv/vlsvwriter.py b/pyVlsv/vlsvwriter.py
index 9f18d9ed..7f2c2fb5 100644
--- a/pyVlsv/vlsvwriter.py
+++ b/pyVlsv/vlsvwriter.py
@@ -41,8 +41,11 @@ def __init__(self, vlsvReader, file_name ):
           :param file_name:     Name of the vlsv file where to input data
       '''
       self.file_name = os.path.abspath(file_name)
-      self.__fptr = open(self.file_name,"wb")
-
+      try:
+         self.__fptr = open(self.file_name,"wb")
+      except FileNotFoundError as e:
+         print("No such path: ", self.file_name)
+         raise e
       self.__xml_root = ET.fromstring("<VLSV></VLSV>")
       self.__fileindex_for_cellid={}
 

From 10e463673a85988dc701bd5f46bf62bd968fd329 Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Thu, 19 May 2022 13:28:02 +0300
Subject: [PATCH 20/25] Commented out unused import, essentially dummy commit

---
 pyCalculations/derivatives.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyCalculations/derivatives.py b/pyCalculations/derivatives.py
index d80a6720..2e3d5bce 100644
--- a/pyCalculations/derivatives.py
+++ b/pyCalculations/derivatives.py
@@ -25,7 +25,7 @@
 import scipy as sp
 import pytools as pt
 import warnings
-from scipy import interpolate
+#from scipy import interpolate
 
 def fg_PoyntingFlux(reader):
    b = reader.read_fg_variable_as_volumetric('fg_b')

From 659da7689305d1b62dca6c812e3db93cdf1b7e35 Mon Sep 17 00:00:00 2001
From: ykempf <yann.kempf@helsinki.fi>
Date: Mon, 30 May 2022 16:53:04 +0300
Subject: [PATCH 21/25] Much faster implementation of read_variable_as_fg

---
 pyVlsv/vlsvreader.py | 32 ++++++++++++++++++++++----------
 1 file changed, 22 insertions(+), 10 deletions(-)

diff --git a/pyVlsv/vlsvreader.py b/pyVlsv/vlsvreader.py
index 890f9e38..06ee6490 100644
--- a/pyVlsv/vlsvreader.py
+++ b/pyVlsv/vlsvreader.py
@@ -1366,19 +1366,31 @@ def upsample_fsgrid_subarray(self, cellid, var, array):
       return
 
    def read_variable_as_fg(self, var):
-      cellids = self.read_variable('CellID')
+      vg_cellids = self.read_variable('CellID')
       sz = self.get_fsgrid_mesh_size()
-      testvar = self.read_variable(var, [cellids[0]])
-      varsize = testvar.size
+      sz_amr = self.get_spatial_mesh_size()
+      vg_var = self.read_variable(var)
+      varsize = vg_var[0].size
       if(varsize > 1):
-         fgarr = np.zeros([sz[0], sz[1], sz[2], varsize], dtype=testvar.dtype)
+         fg_var = np.zeros([sz[0], sz[1], sz[2], varsize], dtype=vg_var.dtype)
       else:
-         fgarr = np.zeros(sz, dtype=testvar.dtype)
-      for c in cellids:
-         self.upsample_fsgrid_subarray(c, var, fgarr)
-         #print
-      return fgarr
-
+         fg_var = np.zeros(sz, dtype=vg_var.dtype)
+      vg_cellids_on_fg = np.zeros(sz, dtype=np.int64) + 1000000000 # big number to catch errors in the latter code, 0 is not good for that
+      current_amr_level = self.get_amr_level(np.min(vg_cellids))
+      max_amr_level = int(np.log2(sz[0] / sz_amr[0]))
+      current_max_cellid = 0
+      for level in range(current_amr_level+1):
+         current_max_cellid += 2**(3*(level))*(self.__xcells*self.__ycells*self.__zcells)
+      for id in np.argsort(vg_cellids):
+         if vg_cellids[id] > current_max_cellid:
+            current_amr_level += 1
+            current_max_cellid += 2**(3*(current_amr_level))*(self.__xcells*self.__ycells*self.__zcells)
+         this_cell_indices = np.array(self.get_cell_indices(vg_cellids[id], current_amr_level), dtype=np.int64)
+         refined_ids_start = this_cell_indices * 2**(max_amr_level-current_amr_level)
+         refined_ids_end = refined_ids_start + 2**(max_amr_level-current_amr_level)
+         vg_cellids_on_fg[refined_ids_start[0]:refined_ids_end[0],refined_ids_start[1]:refined_ids_end[1],refined_ids_start[2]:refined_ids_e>
+      fg_var = vg_var[vg_cellids_on_fg]
+      return fg_var
 
    def get_cell_fsgrid(self, cellid):
       '''Returns a slice tuple of fsgrid indices that are contained in the SpatialGrid

From 3600ce19d1d0517b3d280f340239acad9c81172b Mon Sep 17 00:00:00 2001
From: ykempf <yann.kempf@helsinki.fi>
Date: Mon, 30 May 2022 17:32:57 +0300
Subject: [PATCH 22/25] Fix copy-paste error.

---
 pyVlsv/vlsvreader.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyVlsv/vlsvreader.py b/pyVlsv/vlsvreader.py
index 06ee6490..0fd346b6 100644
--- a/pyVlsv/vlsvreader.py
+++ b/pyVlsv/vlsvreader.py
@@ -1388,7 +1388,7 @@ def read_variable_as_fg(self, var):
          this_cell_indices = np.array(self.get_cell_indices(vg_cellids[id], current_amr_level), dtype=np.int64)
          refined_ids_start = this_cell_indices * 2**(max_amr_level-current_amr_level)
          refined_ids_end = refined_ids_start + 2**(max_amr_level-current_amr_level)
-         vg_cellids_on_fg[refined_ids_start[0]:refined_ids_end[0],refined_ids_start[1]:refined_ids_end[1],refined_ids_start[2]:refined_ids_e>
+         vg_cellids_on_fg[refined_ids_start[0]:refined_ids_end[0],refined_ids_start[1]:refined_ids_end[1],refined_ids_start[2]:refined_ids_end[2]] = id
       fg_var = vg_var[vg_cellids_on_fg]
       return fg_var
 

From f493473d5e3a2e569896a5c17a8d14a2570543e7 Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Wed, 1 Jun 2022 15:07:03 +0300
Subject: [PATCH 23/25] renamed derivatives.py

---
 pyCalculations/calculations.py                        | 2 +-
 pyCalculations/{derivatives.py => differentialops.py} | 0
 2 files changed, 1 insertion(+), 1 deletion(-)
 rename pyCalculations/{derivatives.py => differentialops.py} (100%)

diff --git a/pyCalculations/calculations.py b/pyCalculations/calculations.py
index c7bb589c..cb114147 100644
--- a/pyCalculations/calculations.py
+++ b/pyCalculations/calculations.py
@@ -56,4 +56,4 @@
 import fit
 from fieldtracer import static_field_tracer
 from fieldtracer import dynamic_field_tracer
-from derivatives import fg_divPoynting, fg_PoyntingFlux, fg_vol_curl, ballooning_crit,vfield3_curvature
+from differentialops import fg_divPoynting, fg_PoyntingFlux, fg_vol_curl, ballooning_crit,vfield3_curvature
diff --git a/pyCalculations/derivatives.py b/pyCalculations/differentialops.py
similarity index 100%
rename from pyCalculations/derivatives.py
rename to pyCalculations/differentialops.py

From 9ecc815e08c9c46a10af156399517f6804fea67e Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Tue, 27 Sep 2022 11:19:35 +0300
Subject: [PATCH 24/25] A guard against files with zero populations

---
 pyVlsv/vlsvreader.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyVlsv/vlsvreader.py b/pyVlsv/vlsvreader.py
index 890f9e38..a0550a4e 100644
--- a/pyVlsv/vlsvreader.py
+++ b/pyVlsv/vlsvreader.py
@@ -675,7 +675,7 @@ def read(self, name="", tag="", mesh="", operator="pass", cellids=-1):
          reducer_multipop = multipopdatareducers
             
       # If this is a variable that can be summed over the populations (Ex. rho, PTensorDiagonal, ...)
-      if self.check_variable(self.active_populations[0]+'/'+name): 
+      if len(self.active_populations) > 0 and self.check_variable(self.active_populations[0]+'/'+name): 
          tmp_vars = []
          for pname in self.active_populations:
             vlsvvariables.activepopulation = pname

From 9c765551f9b606527cff8ce5258daa6f53b6cb8e Mon Sep 17 00:00:00 2001
From: Markku Alho <markku.alho@helsinki.fi>
Date: Sun, 26 Nov 2023 14:51:30 +0200
Subject: [PATCH 25/25] cleaning up a bit

---
 pyCalculations/differentialops.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/pyCalculations/differentialops.py b/pyCalculations/differentialops.py
index 2e3d5bce..93f8dc1a 100644
--- a/pyCalculations/differentialops.py
+++ b/pyCalculations/differentialops.py
@@ -141,6 +141,4 @@ def ballooning_crit(reader, B, P, beta):
     
     kappaC = vfield3_dot(nnorm, n)
 
-    balloon = (2 + beta) / 4.0 * kappaP / (kappaC + 1e-27)
-
-    return (balloon, kappaC, gradP)
\ No newline at end of file
+    return (kappaC, gradP)
\ No newline at end of file