diff --git a/_gridFunctions/_zeroField.py b/_gridFunctions/_zeroField.py
new file mode 100644
index 0000000000000000000000000000000000000000..4200f1974461eee648d22914cb486a52dfa18c93
--- /dev/null
+++ b/_gridFunctions/_zeroField.py
@@ -0,0 +1,46 @@
+# class ZeroField;
+# used for OerstedField(v) OerstedField(m_old) if (d/dt)Je = 0
+
+
+################################################################################
+
+
+class ZeroField:
+ 
+
+#------------------------------------------------------------------------------#
+
+  
+  def __init__(self):
+    from ngsolve import CoefficientFunction
+    self.cf = CoefficientFunction((0.0, 0.0, 0.0))
+
+
+#------------------------------------------------------------------------------#
+
+
+  def SetZero(self): pass
+
+
+#------------------------------------------------------------------------------#
+
+
+  def SetFromOther(self, other=None): pass
+
+
+#------------------------------------------------------------------------------#
+
+
+  def Compute(self, solvetol=None): pass
+
+
+#------------------------------------------------------------------------------#
+
+
+  def Save(self, path=None): pass
+
+
+#------------------------------------------------------------------------------#
+
+
+################################################################################
diff --git a/_gridFunctions/derivative.py b/_gridFunctions/derivative.py
index 26a106b455a442ceaabfa00d8c97f1f18697e9bd..fdfd70e55eacfc95f741be1e9832881bd34299ad 100644
--- a/_gridFunctions/derivative.py
+++ b/_gridFunctions/derivative.py
@@ -24,11 +24,17 @@ class Derivative:
     
     if( self.maxwellCoupling == "magnetostaticField" ):
       self.hs = strayField.StrayField(self._mesh, self.cf, solvetol)
-      self.hc = oerstedField.OerstedField(self._mesh, j_cf, solvetol)
+      if j_cf != None:
+        self.hc = oerstedField.OerstedField(self._mesh, j_cf, solvetol)
+      else:
+        self.hc = _zeroField.ZeroField()
     if( self.maxwellCoupling == "strayField" ):
       self.hs = strayField.StrayField(self._mesh, self.cf, solvetol)
     elif( self.maxwellCoupling == "oerstedField" ):
-      self.hc = oerstedField.OerstedField(self._mesh, j_cf, solvetol)
+      if j_cf != None:
+        self.hc = oerstedField.OerstedField(self._mesh, j_cf, solvetol)
+      else:
+        self.hc = _zeroField.ZeroField()
   
 
 #-----------------------------------------------------------------------------#
@@ -43,11 +49,17 @@ class Derivative:
 
     if( self.maxwellCoupling == "magnetostaticField" ):
       self.hs = strayField.StrayField(self._mesh, self.cf, solvetol)
-      self.hc = oerstedField.OerstedField(self._mesh, j_cf, solvetol)
+      if j_cf != None:
+        self.hc = oerstedField.OerstedField(self._mesh, j_cf, solvetol)
+      else:
+        self.hc = _zeroField.ZeroField()
     elif( self.maxwellCoupling == "strayField" ):
       self.hs = strayField.StrayField(self._mesh, self.cf, solvetol)
     elif( self.maxwellCoupling == "oerstedField" ):
-      self.hc = oerstedField.OerstedField(self._mesh, j_cf, solvetol)
+      if j_cf != None:
+        self.hc = oerstedField.OerstedField(self._mesh, j_cf, solvetol)
+      else:
+        self.hc = _zeroField.ZeroField()
     elif( self.maxwellCoupling == "eddyCurrent" ):
       pass
     elif( self.maxwellCoupling != "none" ): 
diff --git a/_gridFunctions/magnetization.py b/_gridFunctions/magnetization.py
index 7ee8e9d7b489d4901aedfe3b7a791bc7bb1662a8..17d23ff04cab098c9f9dd309054f1418e55489e2 100644
--- a/_gridFunctions/magnetization.py
+++ b/_gridFunctions/magnetization.py
@@ -12,7 +12,7 @@ class Magnetization:
   def __init__(self, geometry, V, X, m0, maxwellCoupling, solvetol, j_cf):
     import sys
     from ngsolve import GridFunction
-    from . import strayField, oerstedField
+    from . import strayField, oerstedField, _zeroField
 
     # X is the compound FESpace \subset H1(\Omega)^3
     self._mesh = geometry.GetMesh()
@@ -25,11 +25,17 @@ class Magnetization:
    
     if( self.maxwellCoupling == "magnetostaticField" ):
       self.hs = strayField.StrayField(self._mesh, self.cf, solvetol)
-      self.hc = oerstedField.OerstedField(self._mesh, j_cf, solvetol)
+      if j_cf != None:
+        self.hc = oerstedField.OerstedField(self._mesh, j_cf, solvetol)
+      else:
+        self.hc = _zeroField.ZeroField()
     elif( self.maxwellCoupling == "strayField" ):
       self.hs = strayField.StrayField(self._mesh, self.cf, solvetol)
     elif( self.maxwellCoupling == "oerstedField" ):
-      self.hc = oerstedField.OerstedField(self._mesh, j_cf, solvetol)
+      if j_cf != None:
+        self.hc = oerstedField.OerstedField(self._mesh, j_cf, solvetol)
+      else:
+        self.hc = _zeroField.ZeroField()
     elif( self.maxwellCoupling == "eddyCurrent" ):
       pass
     elif( self.maxwellCoupling != "none" ): 
@@ -109,7 +115,7 @@ class Magnetization:
 
 
   def UpdateMaxwellCoupling(self, newMaxwellCoupling, solvetol, j_cf):
-    from . import strayField, oerstedField
+    from . import strayField, oerstedField, _zeroField
 
     self.maxwellCoupling = newMaxwellCoupling
     self.hs = None
@@ -117,11 +123,17 @@ class Magnetization:
 
     if( self.maxwellCoupling == "magnetostaticField" ):
       self.hs = strayField.StrayField(self._mesh, self.cf, solvetol)
-      self.hc = oerstedField.OerstedField(self._mesh, j_cf, solvetol)
+      if j_cf != None:
+        self.hc = oerstedField.OerstedField(self._mesh, j_cf, solvetol)
+      else:
+        self.hc = _zeroField.ZeroField()
     elif( self.maxwellCoupling == "strayField" ):
       self.hs = strayField.StrayField(self._mesh, self.cf, solvetol)
     elif( self.maxwellCoupling == "oerstedField" ):
-      self.hc = oerstedField.OerstedField(self._mesh, j_cf, solvetol)
+      if j_cf != None:
+        self.hc = oerstedField.OerstedField(self._mesh, j_cf, solvetol)
+      else:
+        self.hc = _zeroField.ZeroField()
     elif( self.maxwellCoupling == "eddyCurrent" ):
       pass
     elif( self.maxwellCoupling != "none" ): 
diff --git a/integrators/_details/physicalToMathematical.py b/integrators/_details/physicalToMathematical.py
index 7e0ba3f78f5e895474f11d144412b182d0508efd..65988062b257641124451b7d5a18e3eb84cb596f 100644
--- a/integrators/_details/physicalToMathematical.py
+++ b/integrators/_details/physicalToMathematical.py
@@ -198,7 +198,7 @@ def ConvertCurrentDensity_ZhangLi(self):
 
     self._jprime_cf = (jnew_cf - self._j_cf) / self._k
   else:
-    self._jprime_cf = CoefficientFunction((0,0,0))
+    self._jprime_cf = None
   
   
 #------------------------------------------------------------------------------#
@@ -222,7 +222,7 @@ def ConvertCurrentDensity_Slonczewski(self):
     / self.parameters.Ms
     self._jprime_cf = (jnew_cf - self._j_cf) / self._k
   else:
-    self._jprime_cf = CoefficientFunction(0)
+    self._jprime_cf = None
 
   
 #------------------------------------------------------------------------------#
diff --git a/integrators/_integrator.py b/integrators/_integrator.py
index f0808c9f40565c3d8110830be9c6a7e960471fdb..bb10abcdacf61e0080c25a6fa9f4ef92f770f60c 100644
--- a/integrators/_integrator.py
+++ b/integrators/_integrator.py
@@ -167,6 +167,12 @@ class _Integrator:
     import sys
     ## TODO search for all the prints and sys.exits --> replace with raise(...)
 
+    if(self.parameters.maxwellCoupling == "oerstedField" \
+      or self.parameters.maxwellCoupling == "magnetostaticField"):
+      print("maxwellCoupling set to", self.parameters.maxwellCoupling)
+      print("Oersted field implementation is Work In Progress; not finished.")
+      sys.exit()
+
     if(False): # TODO check bad combinations
       print("message TODO")
       sys.exit()
diff --git a/integrators/_tpsx.py b/integrators/_tpsx.py
index 1288e5c4b208a7564c0f3fef91ef580738ca5fa4..f38061710e4b92254ab0b811305ee4972fa2ae6d 100644
--- a/integrators/_tpsx.py
+++ b/integrators/_tpsx.py
@@ -345,9 +345,10 @@ class TPSX(_integrator._Integrator):
   def _SetUpAuxiliaryFunctions(self):
     from commix._gridFunctions import Derivative
     if self.parameters.spintronicsCoupling == "slonczewski":
+      jprime_cf = self._jprime_cf*self._p_cf if self._jprime_cf != None else \
+        self._jprime_cf
       self._Mag.v = Derivative(self._GetMesh(), self._V, self._X \
-        , self.parameters.maxwellCoupling, self._solvetol \
-        , self._jprime_cf * self._p_cf)
+        , self.parameters.maxwellCoupling, self._solvetol, jprime_cf)
     else:
       self._Mag.v = Derivative(self._GetMesh(), self._V, self._X \
         , self.parameters.maxwellCoupling, self._solvetol, self._jprime_cf)
diff --git a/integrators/tps2.py b/integrators/tps2.py
index daecf470fe9fd7e9cef14759a06fa6e2d22d9443..9cade8a734491100b80cd2dcc0ae35d8d87d2472 100644
--- a/integrators/tps2.py
+++ b/integrators/tps2.py
@@ -43,6 +43,8 @@ class TPS2(_tpsx.TPSX):
     if (self.parameters.theta != 0.5):
       raise ValueError("theta!=0.5 not available for", type(self))
 
+    _tpsx.TPSX._CheckForBadParameterCombinations(self)
+
 
 #------------------------------------------------------------------------------#
 
diff --git a/integrators/tps2ab.py b/integrators/tps2ab.py
index 47f85373b65061a9ad4309777f4aecdfbd0fe35e..57eb4453fcdf9b9efbd9aeca283c0c105c43dfb1 100644
--- a/integrators/tps2ab.py
+++ b/integrators/tps2ab.py
@@ -54,14 +54,20 @@ class TPS2AB(tps2.TPS2):
   
   def _SetUpOldMagnetization(self):
     from commix._gridFunctions import Magnetization
-    if self.parameters.spintronicsCoupling == "slonczewski":
-      self._Mag_old = Magnetization(self._geometry, self._V, self._X \
-        , self._Mag \
-        , self.parameters.maxwellCoupling, self._solvetol, self._j_cf*self._p_cf)
-    else:
-      self._Mag_old = Magnetization(self._geometry, self._V, self._X \
-        , self._Mag \
-        , self.parameters.maxwellCoupling, self._solvetol, self._j_cf )
+#    if self.parameters.spintronicsCoupling == "slonczewski":
+      # for _Mag_old, the OerstedField is never computed.
+      # of j_prime = None, it's also not needed
+      # TODO: OERSTED take care of this in the future, when Oersted is supported
+      # jprime_cf = self._jprime_cf*self._p_cf if self._jprime_cf != None else \
+      # self._jprime_cf
+#    else:
+#      self._Mag_old = Magnetization(self._geometry, self._V, self._X \
+#        , self._Mag \
+#        , self.parameters.maxwellCoupling, self._solvetol, self._j_cf )
+    jprime_cf = None
+    self._Mag_old = Magnetization(self._geometry, self._V, self._X \
+      , self._Mag \
+      , self.parameters.maxwellCoupling, self._solvetol, jprime_cf)
 
 
 #------------------------------------------------------------------------------#