Merge pull request #188 from alhom/currentReducer2

FHA-style jacobian reducer, incl. vg_J.

pyCalculations/variable.py

@@ -126,6 +126,10 @@ def get_scaled_units(self, vscale=None, env='EarthSpace', manualDict=None):
                  'eV/cm^3': {'defaultScale':1e-3,
                              1e-3:{'scaledUnits':'keV/cm^3', 'scaledLatexUnit':r'$\mathrm{keV}\,\mathrm{cm}^{-3}$'}
                             },
+                 'T/m': {'defaultScale':1e-12,
+                             1e-9:{'scaledUnits':'nT/m', 'scaledLatexUnit':r'$\mathrm{nT}\,\mathrm{m}^{-1}$'},
+                             1e-12:{'scaledUnits':'nT/km', 'scaledLatexUnit':r'$\mathrm{nT}\,\mathrm{km}^{-1}$'}
+                           },
                  'kg/m3'  : {'defaultScale':5.97863741e26,
                              5.97863741e26:{'scaledUnits':'amu/m^3', 'scaledLatexUnit':r'$\mathrm{amu}\,\mathrm{m}^{-3}$'},
                              5.97863741e20:{'scaledUnits':'amu/cm^3', 'scaledLatexUnit':r'$\mathrm{amu}\,\mathrm{cm}^{-3}$'}

pyVlsv/reduction.py

@@ -460,6 +460,44 @@ def PerpendicularTensorComponent( variables ):
    else:
       return 0.5*(RotatedTensor[:,0,0] + RotatedTensor[:,1,1])
 
+def J( variables ):
+   ''' Data reducer taking a jacobian (assume 9-component vector) and extracting the current
+   via curl from the components of the jacobian (background or perturbed, as long as it has 9
+    components)
+
+   '''
+   stack = True
+   if(variables[0].shape == (9,)):
+      stack = False
+      variables[0] = np.array([variables[0]]) # I want this to be a stack of tensors
+   jacobian = variables[0]#.reshape((variables[0].shape[0],3,3))
+   # jacobian = jacobian.transpose((0,2,1)) # The axes are flipped at some point, correcting for that
+   J = np.zeros((jacobian.shape[0],3))
+
+   J[:,0] = (jacobian[:,7]-jacobian[:,5])/mu_0
+   J[:,1] = (jacobian[:,2]-jacobian[:,6])/mu_0
+   J[:,2] = (jacobian[:,3]-jacobian[:,1])/mu_0
+   if stack:
+      return J
+   else:
+      return J[0,:]
+
+
+   print("Error in J")
+   return -1
+
+def TensorFromScalars(variables):
+   '''Construct a 9-element vector ("tensor") from nine scalar fields.
+   '''
+
+   return np.stack(np.array(
+                    [variables[0], variables[1], variables[2],
+                     variables[3], variables[4], variables[5],
+                     variables[6], variables[7], variables[8]]
+                   ),
+                   axis=-1)
+
+
 def Anisotropy( variables ):
    ''' Data reducer for finding the ratio of perpendicular to parallel components of a tensor
    '''
@@ -1028,6 +1066,11 @@ def makelambda(index):
 
 
 
+v5reducers["vg_jacobian_b"] =             DataReducerVariable(["vg_dbxvoldx","vg_dbxvoldy","vg_dbxvoldz","vg_dbyvoldx","vg_dbyvoldy","vg_dbyvoldz","vg_dbzvoldx","vg_dbzvoldy","vg_dbzvoldz"], TensorFromScalars, "T/m", 1, latex=r"$\vec{J}$",latexunits=r"$\mathrm{A}\,\mathrm{m}^{-2}$")
+v5reducers["vg_jacobian_bper"] =          DataReducerVariable(["vg_dperbxvoldx","vg_dperbxvoldy","vg_dperbxvoldz","vg_dperbyvoldx","vg_dperbyvoldy","vg_dperbyvoldz","vg_dperbzvoldx","vg_dperbzvoldy","vg_dperbzvoldz"], TensorFromScalars, "T/m", 1, latex=r"$\vec{J}$",latexunits=r"$\mathrm{A}\,\mathrm{m}^{-2}$")
+v5reducers["vg_j"] =                     DataReducerVariable(["vg_jacobian_bper"], J, "A/m^2", 1, latex=r"$\vec{J}$",latexunits=r"$\mathrm{A}\,\mathrm{m}^{-2}$")
+
+
 #multipopv5reducers
 multipopv5reducers = {}
 multipopv5reducers["pop/vg_rhom"] =                   DataReducerVariable(["pop/vg_rho"], rhom, "kg/m3", 1, latex=r"$\rho_{m,\mathrm{REPLACEPOP}}$",latexunits=r"$\mathrm{kg}\,\mathrm{m}^{-3}$")