diff --git a/_gridFunctions/magnetization.py b/_gridFunctions/magnetization.py
index a04ca33c382e51d7cbe228b7e1d813375d677879..e3e79e6ffe1d457bc6e8006d966a5d0396d0d719 100644
--- a/_gridFunctions/magnetization.py
+++ b/_gridFunctions/magnetization.py
@@ -156,19 +156,6 @@ class Magnetization:
     self.hs.Compute(solvetol)
   
 
-#------------------------------------------------------------------------------#
-
-
-  def InitializeEddyCurrent(self):
-    import sys
-    if ( self.maxwellCoupling != "eddyCurrent" ):
-      print("Error in _Mag.ComputeEddyCurrent():")
-      print("maxwellCoupling != <eddyCurrent> !! --> Terminate.")
-      sys.exit()
-    
-    self.H.SetFromMagnetization(self.m)
-  
-
 #------------------------------------------------------------------------------#
 
 
@@ -343,7 +330,6 @@ class Magnetization:
   def EnergyStrayField(self):
     import sys
     from ngsolve import Integrate
-  # TOECDO zu energy maxwell zusammenlegen
     if ( self.maxwellCoupling != "strayField" \
       and self.maxwellCoupling != "magnetostaticField" ):
       print("Error in _Mag.EnergyStrayField():")
@@ -354,21 +340,6 @@ class Magnetization:
     return energy
   
   
-#------------------------------------------------------------------------------#
-
-
-  def EnergyEddyCurrent(self):
-    import sys
-    from ngsolve import Integrate
-    if ( self.maxwellCoupling != "eddyCurrent" ):
-      print("Error in _Mag.EnergyEddyCurrent():")
-      print("maxwellCoupling != <eddyCurrent> !! --> Terminate.")
-      sys.exit()
-    
-    energy = Integrate( self.H.gf * self.H.gf, self.H.GetMesh())
-    return energy
-  
-  
 #------------------------------------------------------------------------------#
 
 
@@ -387,8 +358,6 @@ class Magnetization:
   def Save(self, filepath):
     from commics._tools import SaveMagnetization
     SaveMagnetization(filepath, self._mesh, self.m)
-#    if ( self.maxwellCoupling == "eddyCurrent" ): # TOECDO
-#      self.H.Save(filepath)
   
   
 #------------------------------------------------------------------------------#
diff --git a/_setup/geometries/__init__.py b/_setup/geometries/__init__.py
index 88c6dce3496964c8a5465e77a7fe5552bd7f4a0d..b8c9d9567b9261dc125ebad1219cf18a34d0df21 100644
--- a/_setup/geometries/__init__.py
+++ b/_setup/geometries/__init__.py
@@ -1,10 +1,8 @@
 print("Importing module", __name__)
 
-#from commics._setup.geometries.boundingBox import BoundingBox TOECDO
 from commics._setup.geometries.geometry import Geometry
 from commics._setup.geometries.standardGeometries import *
 
-#__all__ = ["BoundingBox", "Geometry"] TOECDO
 __all__ = ["Geometry"]
 
 for item in standardGeometries.__all__:
diff --git a/_setup/geometries/_details.py b/_setup/geometries/_details.py
deleted file mode 100644
index d8e9fbfdb4e3222da3a9b025d260e846a5bf5878..0000000000000000000000000000000000000000
--- a/_setup/geometries/_details.py
+++ /dev/null
@@ -1,110 +0,0 @@
-# TOECDO
-#
-#################################################################################
-## suppose innerMesh is a true interior of  outerMesh,
-## in the sense that it is subset and that the bounderies are disjoint
-##        _______________________________
-##       |                               |
-##       |               1               |
-##       |        _______________        |
-##       |       |               |       |
-##       |       |               |       |
-##       |       |       2       |       |
-##       |       |               |       |
-##       |       |_______________|       |
-##       |                               |
-##       |_______________________________|
-##
-#################################################################################
-#
-#
-#def ExtractInnerMesh(combinedMesh, nr_outside=1, nr_inside=2, nr_void=0):
-#    from netgen.meshing import Mesh, FaceDescriptor, Element3D, Element2D
-#  
-#  mesh = Mesh()
-#  
-#  nr_outer_dom = nr_outside
-#  nr_outer_bc = nr_outside
-#  nr_outer_surf = nr_outside
-#  
-#  nr_inner_dom = nr_inside
-#  nr_inner_bc = nr_inside
-#  nr_inner_surf = nr_inside
-#  
-#  restrictionMap = []
-#  extensionMap = {}
-#  
-#  for K in combinedMesh.Elements3D():
-#    if ( K.index == nr_inner_dom ):
-#      for v in K.vertices:
-#        if( v not in extensionMap ):
-#          extensionMap[v] = mesh.Add( combinedMesh[v] )
-#          restrictionMap.append( v.nr-1 )
-#      
-#      mesh.Add( Element3D(index=nr_inner_dom \
-#        , vertices=[extensionMap[v] for v in K.vertices] ) )
-#  
-#  
-#  ## boundary conditions
-#  fd_inner_boundary = mesh.Add( FaceDescriptor(surfnr=nr_inner_surf \
-#    , domin=nr_inner_dom, domout=nr_void, bc=nr_inner_bc) )
-#  
-#  for E in combinedMesh.Elements2D():
-#    if ( E.index == nr_inner_bc ):
-#      mesh.Add( Element2D(index=fd_inner_boundary \
-#        , vertices=[extensionMap[v] for v in E.vertices] ) )
-#  
-#  return mesh, restrictionMap
-#
-#
-#################################################################################
-#
-#
-#def CombineMeshes(innerMesh, outerMesh, nr_outside=1, nr_inside=2, nr_void=0):
-#    from netgen.meshing import Mesh, FaceDescriptor, Element2D
-#  # adapted from 
-#  # https://gitlab.asc.tuwien.ac.at/jschoeberl/ngsolve-docu/wikis/ngpy-meshes 
-#  # on 2017/03/13
-#  
-#  mesh = Mesh()
-#  
-#  nr_outer_dom = nr_outside
-#  nr_outer_bc = nr_outside
-#  nr_outer_surf = nr_outside
-#  
-#  nr_inner_dom = nr_inside
-#  nr_inner_bc = nr_inside
-#  nr_inner_surf = nr_inside
-#  fd_outside = mesh.Add( FaceDescriptor(surfnr=nr_outer_surf \
-#    , domin=nr_outer_dom,  bc=nr_outer_bc) )
-#  fd_inside = mesh.Add( FaceDescriptor(surfnr=nr_inner_surf \
-#    , domin=nr_inner_dom, domout=nr_outer_dom, bc=nr_inner_bc) )
-#  
-#  
-#  pmap1 = { }
-#  for e in outerMesh.Elements2D():
-#    for v in e.vertices:
-#      if (v not in pmap1):
-#        pmap1[v] = mesh.Add( outerMesh[v] )
-#  
-#  
-#  for e in outerMesh.Elements2D():
-#    mesh.Add( Element2D(fd_outside, [pmap1[v] for v in e.vertices]) )
-#  
-#  
-#  pmap2 = { }
-#  for e in innerMesh.Elements2D():
-#    for v in e.vertices:
-#      if (v not in pmap2):
-#        pmap2[v] = mesh.Add( innerMesh[v] )
-#  
-#  
-#  for e in innerMesh.Elements2D():
-#    mesh.Add( Element2D(fd_inside, [pmap2[v] for v in e.vertices]) )
-#  
-#  mesh.GenerateVolumeMesh()
-#  
-#  return mesh
-#
-#
-#################################################################################
diff --git a/_setup/geometries/boundingBox.py b/_setup/geometries/boundingBox.py
deleted file mode 100644
index 7c3d8a215b5dc7369435d7b7e3b09c5ed8c88e69..0000000000000000000000000000000000000000
--- a/_setup/geometries/boundingBox.py
+++ /dev/null
@@ -1,116 +0,0 @@
-#### class BoundingBox ### TOECDO
-#
-#
-#################################################################################
-#
-#
-#class BoundingBox:
-#  
-#  def __init__(self, geometry, h_max):
-#    
-#    self.geometry = geometry
-#    self.h_max = h_max
-#    
-#    self._mesh = None
-#    
-#    self._outer_domain_nr = 1
-#    self._inner_domain_nr = 2
-#    self._restrictionMap = None
-#
-#    self._scaling = None
-#  
-#  
-##------------------------------------------------------------------------------#
-#
-#
-#  def GetReady(self, geometry_in, h_max_in, scaling_in):
-#    from commics._setup.geometries._details import ExtractInnerMesh
-#    from commics._setup.geometries._details import CombineMeshes
-#    from ngsolve import Mesh as ngsolveMesh
-#    from netgen.csg import CSGeometry
-#    self._scaling = scaling_in
-#
-#    geo_in = CSGeometry()
-#    geo_in.Add( geometry_in )
-#    ngmesh_in = geo_in.GenerateMesh( maxh = h_max_in )
-#    
-#    geo_out = CSGeometry()
-#    geo_out.Add( self.geometry )
-#    ngmesh_out = geo_out.GenerateMesh( maxh = self.h_max )
-#    
-#    ngmesh = CombineMeshes(ngmesh_in, ngmesh_out \
-#      , self._outer_domain_nr, self._inner_domain_nr)
-#    ngmesh_restricted, self._restrictionMap = ExtractInnerMesh(ngmesh \
-#      , self._outer_domain_nr, self._inner_domain_nr)
-#    
-#    self._mesh = ngsolveMesh( ngmesh )
-#    return ngsolveMesh( ngmesh_restricted )
-#  
-#  
-##------------------------------------------------------------------------------#
-#
-#
-#  def IsBoundingBox(self):
-#    return True
-#  
-#  
-##------------------------------------------------------------------------------#
-#
-#
-#  def Get_restrictionMap(self):
-#    return self._restrictionMap
-#  
-#  
-##------------------------------------------------------------------------------#
-#
-#
-#  def Get_inner_domain_nr(self):
-#    return self._inner_domain_nr
-#  
-#  
-##------------------------------------------------------------------------------#
-#
-#
-#  def Get_outer_domain_nr(self):
-#    return self._outer_domain_nr
-#  
-#  
-##------------------------------------------------------------------------------#
-#
-#
-#  def __str__(self):
-#    string = "\n\n ----- " + self.__class__.__name__ + " ----- \n\n"
-#    if not(self._mesh == None):
-#      string += "Vertices: " + str(self._mesh.nv) + "\n"
-#      string += "Elements: " + str(self._mesh.ne) + "\n"
-#      string += "h_max = " + str(self.h_max) + "\n"
-#      string += "scaling = " + str(self._scaling) + "\n"
-#    else:
-#      string = "From BoundingBox.__str__(): \
-#        Mesh not ready: Call GetReady() first!"
-#  
-#    return string
-#  
-##------------------------------------------------------------------------------#
-#
-#
-#  def GetMesh(self):
-#    if not(self._mesh == None):
-#      return self._mesh
-#    else:
-#      print("From BoundingBox.GetMesh() : \
-#        Mesh not ready: Call GetReady() first!")
-#  
-#  
-##------------------------------------------------------------------------------#
-#
-#
-#  def Save(self, filename):
-#    if not(self._mesh == None):
-#      self._mesh.ngmesh.Save( filename )
-#    else:
-#      print("From BoundingBox.Save(): \
-#        Mesh not ready: Call GetReady() first!")
-#
-#
-#################################################################################
diff --git a/_setup/geometries/geometry.py b/_setup/geometries/geometry.py
index 7e593e526c1f313233d1273211201d35f0ee1822..8f20786a319b053381d59f02e0cd1c3eb36478d5 100644
--- a/_setup/geometries/geometry.py
+++ b/_setup/geometries/geometry.py
@@ -1,4 +1,3 @@
-#from . import boundingBox # TOECDO
 
 
 ################################################################################
@@ -57,8 +56,6 @@ class Geometry:
     self.structuredMesh = structuredMesh
     self.structuredObtuseMesh = structuredObtuseMesh
     
-#TOECDO    self.boundingBox = None
-
     self._mesh = None
     
 
@@ -79,15 +76,6 @@ class Geometry:
       print("From Geometry.GetReady(): Nothing to do, mesh is ready.")
       return
 
-    
-    # TOECDO
-#    if (self.boundingBox == None):
-#      self.boundingBox = self
-#    else:
-#      self._mesh = self.boundingBox.GetReady(self._geometry, self.meshSize \
-#        , self.scaling)
-#      return
-
     if(isinstance(self._geometry, str)):
       if(self.meshSize != None):
         print("Error in Geometry.GetReady():")
@@ -155,45 +143,6 @@ class Geometry:
       return StructuredObtuseMesh_Cuboid(p1, p2, self.meshSize)
 
 
-
-#------------------------------------------------------------------------------#
-# TOECDO
-#  
-#  def IsBoundingBox(self):
-#    return False
-#  
-#  
-#------------------------------------------------------------------------------#
-# TOECDO
-#  
-#  def Get_restrictionMap(self):
-#    return range(0,self._mesh.nv)
-#  
-#  
-#------------------------------------------------------------------------------#
-# TOECDO
-#  
-#  def Get_inner_domain_nr(self):
-#    import sys
-#    print("Error in Geometry.Get_inner_domain_nr():")
-#    print("self.IsBoundingBox() = ", self.IsBoundingBox())
-#    print("Something went wrong with the meshes:")
-#    print("This method should never be called!")
-#    sys.exit()
-#  
-#  
-#------------------------------------------------------------------------------#
-# TOECDO
-#  
-#  def Get_outer_domain_nr(self):
-#    import sys
-#    print("Error in Geometry.Get_outer_domain_nr():")
-#    print("self.IsBoundingBox() = ", self.IsBoundingBox())
-#    print("Something went wrong with the meshes:")
-#    print("This method should never be called!")
-#    sys.exit()
-#  
-#  
 #------------------------------------------------------------------------------#
 
   
@@ -209,9 +158,6 @@ class Geometry:
 
       if isinstance(self._geometry, str):
         string += "file = " + str(self._geometry) + "\n"
-# TOECDO
-#      if ( self.boundingBox.IsBoundingBox() ):
-#        string += self.boundingBox.__str__()
 
       return string
     else:
@@ -240,11 +186,6 @@ class Geometry:
       if(filename[-4:] == ext):
         ext = ""
       self._mesh.ngmesh.Save( filename + ext )
-
-# TOECDO
-#      if ( self.boundingBox.IsBoundingBox() ):
-#        tmp = filename.rfind("/")+1
-#        self.boundingBox.Save( filename[:tmp]+"boundingBox_"+filename[tmp:] )
     else:
       print("From Geometry.Save(): \
         Mesh not ready: Call GetReady() first!")
diff --git a/_setup/outputOptions.py b/_setup/outputOptions.py
index f3ced67c1f25b36cc017da06606f15c9069008f2..2847b6bdd3357a6942a09eeb493f851afcb9b1a0 100644
--- a/_setup/outputOptions.py
+++ b/_setup/outputOptions.py
@@ -65,7 +65,6 @@ class OutputOptions:
         + "energy_strayField,".ljust(self.columnWidth) + " " \
         + "energy_applied,".ljust(self.columnWidth) + " " \
         + "energy_total,".ljust(self.columnWidth)
-#TOECDO        + "energy_eddyCurrent,".ljust(self.columnWidth) + " " \
       
       self.headerUnprojectedExchangeEnergy \
         = "columns are: t,".ljust(self.columnWidth-2) + " " \
diff --git a/_setup/parameters.py b/_setup/parameters.py
index 2e405199a87e29e930f3cbd191658a6bcd12af99..ec6691f172e6abf8210698a9e64a6463d37f8d3f 100644
--- a/_setup/parameters.py
+++ b/_setup/parameters.py
@@ -65,8 +65,6 @@ class Parameters:
     self.dmCoupling = dmCoupling
     self.dmImplicitTreatment = dmImplicitTreatment
     self.maxwellCoupling = maxwellCoupling
-#    self.sigma_in = 1.0 #TOECDO
-#    self.sigma_out = None #TOECDO
     
     # spintronics par
     self.spintronicsCoupling = spintronicsCoupling
diff --git a/integrators/_details/interfaces/writeInterface.py b/integrators/_details/interfaces/writeInterface.py
index 3e16d036b352067bcef76102e55edd3f0e3f88da..830af11db0bcc5a8480c9f7e107604766aee3da6 100644
--- a/integrators/_details/interfaces/writeInterface.py
+++ b/integrators/_details/interfaces/writeInterface.py
@@ -224,22 +224,13 @@ class WriteInterface(_interface._Interface):
     else:
       E_stray = 0.0
     
-#    if (self.parameters.maxwellCoupling == "eddyCurrent"):
-#      E_eddy = 0.5 * self.parameters.mu0 \
-#        * self.parameters.Ms * self.parameters.Ms \
-#        * self._geometry.scaling * self._geometry.scaling \
-#        * self._geometry.scaling * self._Mag.EnergyEddyCurrent()
-#    else:
-#      E_eddy = 0.0
-    
     E_appl = -1.0 * self.parameters.mu0 \
       * self.parameters.Ms * self.parameters.Ms \
       * self._geometry.scaling * self._geometry.scaling \
       * self._geometry.scaling * self._Mag.EnergyApplied(self._hext_cf)
       
-    E_total = E_exch + E_ani + E_dm + E_stray + E_appl #TOECDO + E_eddy
+    E_total = E_exch + E_ani + E_dm + E_stray + E_appl
     
-#TOECDO    en = [E_exch, E_ani, E_dm, E_stray, E_eddy, E_appl, E_total]
     en = [E_exch, E_ani, E_dm, E_stray, E_appl, E_total]
     
     LiveSaveValues(self.output.foldername+"/"+self.output.filenameEnergies \
diff --git a/integrators/kw1.py b/integrators/kw1.py
index e380cc29b80cf49a6d28fd5f1313f1645369c31d..0b7ebf611fb9b08cb832a2d1a4bc6a2a81beba86 100644
--- a/integrators/kw1.py
+++ b/integrators/kw1.py
@@ -302,9 +302,6 @@ class KW1(_integrator._Integrator, \
     if( self.parameters.maxwellCoupling == "magnetostaticField" \
       or self.parameters.maxwellCoupling == "strayField" ):
       self._Mag.v.hs.SetZero()
-    elif( self.parameters.maxwellCoupling == "eddyCurrent" ):
-      pass #TOECDO
-#      self.H.SetEtaToH()
     
     # assemble blf and lf
     self._Compute_LHS_instat()
@@ -331,12 +328,6 @@ class KW1(_integrator._Integrator, \
       diff_gf.vec.FV().NumPy()[:] = self._Mag.v.gf.vec.FV().NumPy()-oldIterate
       itError = sqrt( Integrate(diff_cf*diff_cf, self._GetMesh()) )
     
-      if ( self.parameters.maxwellCoupling == "eddyCurrent" ):
-        pass # TOECDO
-#        self.H.SolveForEta(self.m, self.v.gf)
-#        itError += self.H.IterationError()
-    
-      # TOECDO consider different error for Eddy Currents?
       if (itError <= self._itertol):
         break
     
diff --git a/integrators/mps.py b/integrators/mps.py
index 32dbc93391cd898180e29e45264c34fe05a51fde..b4048494e03b79f059b85a3b42e28da5d15106c6 100644
--- a/integrators/mps.py
+++ b/integrators/mps.py
@@ -341,7 +341,6 @@ class MPS(_integrator._Integrator, \
       diff_gf.vec.FV().NumPy()[:] = self._Mag_eta.m.vec.FV().NumPy()-oldIterate
       itError = sqrt( Integrate(diff_cf*diff_cf, self._GetMesh()) )
     
-      # TOECDO consider different error for Eddy Currents?
       # TODO consider error in heff like in mpsPaper
       if (itError <= self._itertol):
         break
diff --git a/integrators/tps2.py b/integrators/tps2.py
index 5af8bd548caa3b5e08d67427600ade0471c72a37..23a31979fb5854b74583ab31bd01848ce4b161dc 100644
--- a/integrators/tps2.py
+++ b/integrators/tps2.py
@@ -242,12 +242,6 @@ class TPS2(_tpsx.TPSX, \
       diff_gf.vec.FV().NumPy()[:] = self._Mag.v.gf.vec.FV().NumPy()-oldIterate
       itError = sqrt( Integrate(diff_cf*diff_cf, self._GetMesh()) )
     
-      if ( self.parameters.maxwellCoupling == "eddyCurrent" ):
-        pass # TOECDO
-#        self.H.SolveForEta(self.m, self.v.gf)
-#        itError += self.H.IterationError()
-    
-      # TOECDO consider different error for Eddy Currents?
       if (itError <= self._itertol):
         break