diff --git a/source_code/conductor.py b/source_code/conductor.py
index 90973c3..0f46f07 100644
--- a/source_code/conductor.py
+++ b/source_code/conductor.py
@@ -14,6 +14,7 @@
 # import classes
 from component_collection import ComponentCollection
 from conductor_flags import (
+    CONSTANT_INDUCTANCE,
     ANALYTICAL_INDUCTANCE,
     APPROXIMATE_INDUCTANCE,
     ELECTRIC_CONDUCTANCE_UNIT_LENGTH,
@@ -22,6 +23,7 @@
     IOP_CONSTANT,
     IOP_FROM_EXT_FUNCTION,
     IOP_FROM_FILE,
+    SELF_INDUCTANCE_MODE_0,
     SELF_INDUCTANCE_MODE_1,
     SELF_INDUCTANCE_MODE_2,
     STATIC_ELECTRIC_SOLVER,
@@ -182,7 +184,7 @@ def __init__(self: Self, simulation: object, sheetConductorsList: list, ICOND: i
             "INDUCTANCE_MODE",
             "ELECTRIC_SOLVER",
         ]
-        self.operations.update({key: 1 for key in keys if self.operations[key]})
+        self.operations.update({key: 1 for key in keys if self.operations[key]==True})
         self.operations.update(
             {key: 0 for key in keys if self.operations[key] == False}
         )
@@ -3688,6 +3690,62 @@ def build_electric_known_term_vector(self):
             self.total_elements_current_carriers :
         ] = self.dict_node_pt["op_current"]
 
+    # CONSTANT INDUCTANCE
+    def __constant_inductance(self, mode: int):
+        """Private method that assigns a constant value to the mutual inductance as defined by the user in sheet CONDUCTOR_operation of the input file conductor_definition.xlsx
+
+        Args:
+            mode (int): flag to select the equation for the analytical evaluation of self inductance. 0:constan value from sheet CONDUCTOR_operation of the input file conductor_definition.xlsx; 1: from method __self_inductance_mode1; 2: from method __self_inductance_mode2.
+        """
+
+        lmod = (
+            (
+                (
+                    self.nodal_coordinates.iloc[
+                        self.connectivity_matrix.loc[
+                            "StrandComponent",
+                            "end",
+                        ],
+                        :,
+                    ]
+                    - self.nodal_coordinates.iloc[
+                        self.connectivity_matrix.loc[
+                            "StrandComponent",
+                            "start",
+                        ],
+                        :,
+                    ]
+                )
+                ** 2
+            )
+            .sum(axis=1)
+            .apply(np.sqrt)
+        )
+        mutual_inductance = self.operations["MUTUAL_INDUCTANCE"] * np.ones(self.inductance_matrix.shape) 
+       
+       # The principal diagonal is set to 0
+        for ii in range(mutual_inductance.shape[0]):
+           mutual_inductance[ii, ii] = 0
+      
+        self_inductance_switch = {
+            SELF_INDUCTANCE_MODE_0: self.__constant_self_inductance_evaluation,
+            SELF_INDUCTANCE_MODE_1: self.__self_inductance_mode1,
+            SELF_INDUCTANCE_MODE_2: self.__self_inductance_mode2,
+        }
+        self_inductance = self_inductance_switch[mode](lmod)
+
+        # Evaluate internal inductance
+        internal_inductance = lmod.to_numpy() / 2.0
+
+        self.inductance_matrix = (
+            constants.mu_0
+            / (4.0 * constants.pi)
+            * (
+                np.diag(self_inductance + internal_inductance)
+                + mutual_inductance
+                + mutual_inductance.T
+            )
+        )
     # START: INDUCTANCE ANALYTICAL EVALUATION
 
     def __inductance_analytical_calculation(self, mode: int = 2):
@@ -3697,9 +3755,9 @@ def __inductance_analytical_calculation(self, mode: int = 2):
             mode (int,optional): flag to select the equation for the analytical evaluation of self inductance. 1: mode 1; 2: mode 2. Defaults to 2.
         """
 
-        if mode != 1 and mode != 2:
+        if mode!= SELF_INDUCTANCE_MODE_0 and mode != SELF_INDUCTANCE_MODE_1 and mode != SELF_INDUCTANCE_MODE_2:
             raise ValueError(
-                f"{self.identifier}\nArgument 'mode' must be equal to {SELF_INDUCTANCE_MODE_1 = } or to {SELF_INDUCTANCE_MODE_2 = }. Current value {mode = } is not allowed. Please check sheet {self.workbook_sheet_name[2]} in file {self.workbook_name}.\n"
+                f"{self.identifier}\nArgument 'mode' must be equal to {SELF_INDUCTANCE_MODE_0 = } or to {SELF_INDUCTANCE_MODE_1 = } or to {SELF_INDUCTANCE_MODE_2 = }. Current value {mode = } is not allowed. Please check sheet {self.workbook_sheet_name[2]} in file {self.workbook_name}.\n"
             )
         ABSTOL = 1e-6
         lmod = (
@@ -3743,6 +3801,7 @@ def __inductance_analytical_calculation(self, mode: int = 2):
 
         # Switch to evalutae self inductance.
         self_inductance_switch = {
+            SELF_INDUCTANCE_MODE_0: self.__constant_self_inductance_evaluation,
             SELF_INDUCTANCE_MODE_1: self.__self_inductance_mode1,
             SELF_INDUCTANCE_MODE_2: self.__self_inductance_mode2,
         }
@@ -3907,6 +3966,21 @@ def __vertex_to_vertex_distance(
             .apply(np.sqrt)
         )
 
+    #  CONSTANT SELF INDUCTANCE 
+    def __constant_self_inductance_evaluation(self, lmod: np.array) -> np.ndarray:
+        """Private method that assigns a constant value to the self inductance that is defined by the user in sheet CONDICTOR_operation in the input file conductor_definition.xlsx
+
+        Args:
+            lmod (np.ndarray): array with the distance between strand component nodal nodes.
+
+        Returns:
+            np.ndarray: self inductances.
+        """
+        
+        self_inductance = np.ones(lmod.shape) * self.operations["SELF_INDUCTANCE"]
+
+        return self_inductance
+
     def __self_inductance_mode1(self, lmod: np.ndarray) -> np.ndarray:
         """Private method that analytically evaluates self inductances according to mode 1.
 
@@ -3920,23 +3994,23 @@ def __self_inductance_mode1(self, lmod: np.ndarray) -> np.ndarray:
 
         for ii, obj in enumerate(self.inventory["StrandComponent"].collection):
             self_inductance[ii :: self.inventory["StrandComponent"].number] = (
-                2
-                * lmod[ii :: self.inventory["StrandComponent"].number]
-                * (
-                    np.arcsinh(
-                        lmod[ii :: self.inventory["StrandComponent"].number]
-                        / obj.radius
-                    )
-                    - np.sqrt(
-                        1.0
-                        + (
-                            obj.radius
-                            / lmod[ii :: self.inventory["StrandComponent"].number]
-                        )
-                        ** 2
-                    )
-                    + obj.radius / lmod[ii :: self.inventory["StrandComponent"].number]
-                )
+               2
+               * lmod[ii :: self.inventory["StrandComponent"].number]
+               * (
+                   np.arcsinh(
+                       lmod[ii :: self.inventory["StrandComponent"].number]
+                       / obj.radius
+                   )
+                   - np.sqrt(
+                       1.0
+                       + (
+                           obj.radius
+                           / lmod[ii :: self.inventory["StrandComponent"].number]
+                       )
+                       ** 2
+                   )
+                   + obj.radius / lmod[ii :: self.inventory["StrandComponent"].number]
+               )
             )
         return self_inductance
 
@@ -3952,25 +4026,26 @@ def __self_inductance_mode2(self, lmod: np.ndarray) -> np.ndarray:
 
         for ii, obj in enumerate(self.inventory["StrandComponent"].collection):
 
-            self_inductance[ii :: self.inventory["StrandComponent"].number] = 2 * (
-                lmod[ii :: self.inventory["StrandComponent"].number]
-                * np.log(
-                    (
-                        lmod[ii :: self.inventory["StrandComponent"].number]
-                        + np.sqrt(
-                            lmod[ii :: self.inventory["StrandComponent"].number] ** 2
-                            + obj.radius ** 2
-                        )
-                    )
-                    / obj.radius
-                )
-                - np.sqrt(
-                    lmod[ii :: self.inventory["StrandComponent"].number] ** 2
-                    + obj.radius ** 2
-                )
-                + lmod[ii :: self.inventory["StrandComponent"].number] / 4
-                + obj.radius
-            )
+            self_inductance[ii :: self.inventory["StrandComponent"].number] = (
+           2 * (
+               lmod[ii :: self.inventory["StrandComponent"].number]
+               * np.log(
+                   (
+                       lmod[ii :: self.inventory["StrandComponent"].number]
+                       + np.sqrt(
+                           lmod[ii :: self.inventory["StrandComponent"].number] ** 2
+                           + obj.radius ** 2
+                       )
+                   )
+                   / obj.radius
+               )
+               - np.sqrt(
+                   lmod[ii :: self.inventory["StrandComponent"].number] ** 2
+                   + obj.radius ** 2
+               )
+               + lmod[ii :: self.inventory["StrandComponent"].number] / 4
+               + obj.radius
+            ))
 
         return self_inductance
 
@@ -3978,8 +4053,13 @@ def __self_inductance_mode2(self, lmod: np.ndarray) -> np.ndarray:
 
     # START: INDUCTANCE APPROXIMATE EVALUATION
 
-    def __inductance_approximate_calculation(self):
+    def __inductance_approximate_calculation(self, mode : int = 2):
         """Private method that approximate the inductance of the system. For an analytical evaluation of the inductance use private method __inductance_analytical_calculation."""
+        
+        if mode!= SELF_INDUCTANCE_MODE_0 and mode != SELF_INDUCTANCE_MODE_1 and mode != SELF_INDUCTANCE_MODE_2:
+            raise ValueError(
+                f"{self.identifier}\nArgument 'mode' must be equal to {SELF_INDUCTANCE_MODE_0 = } or to {SELF_INDUCTANCE_MODE_1 = } or to {SELF_INDUCTANCE_MODE_2 = }. Current value {mode = } is not allowed. Please check sheet {self.workbook_sheet_name[2]} in file {self.workbook_name}.\n"
+            )
 
         ll = (
             self.nodal_coordinates.iloc[
@@ -4013,8 +4093,13 @@ def __inductance_approximate_calculation(self):
                 ll, mutual_inductance
             )
 
-        # Evaluate self inductance
-        self_inductance = self.__self_inductance_approximate(lmod)
+        self_inductance_switch = {
+            SELF_INDUCTANCE_MODE_0: self.__constant_self_inductance_evaluation,
+            SELF_INDUCTANCE_MODE_1: self.__self_inductance_mode1,
+            SELF_INDUCTANCE_MODE_2: self.__self_inductance_mode2,
+        }
+        self_inductance = self_inductance_switch[mode](lmod)
+
         # Evaluate internal inductance
         internal_inductance = lmod.to_numpy() / 2.0
 
@@ -4139,19 +4224,21 @@ def __build_electric_mass_matrix(self):
         """
 
         if (
-            self.operations["INDUCTANCE_MODE"] != 0
-            and self.operations["INDUCTANCE_MODE"] != 1
+            self.operations["INDUCTANCE_MODE"] != CONSTANT_INDUCTANCE
+            and self.operations["INDUCTANCE_MODE"] != ANALYTICAL_INDUCTANCE
+            and self.operations["INDUCTANCE_MODE"] != APPROXIMATE_INDUCTANCE
         ):
             raise ValueError(
-                f"{self.identifier = }\nArgument self.operations['INDUCTANCE_MODE'] should be equal to {APPROXIMATE_INDUCTANCE = } or {ANALYTICAL_INDUCTANCE = }. Current value ({self.operations['INDUCTANCE_MODE'] = }) is not allowed. Please check {self.workbook_sheet_name[2]} in file {self.workbook_name}.\n"
+                f"{self.identifier = }\nArgument self.operations['INDUCTANCE_MODE'] should be equal to {CONSTANT_INDUCTANCE = } or {APPROXIMATE_INDUCTANCE = } or {ANALYTICAL_INDUCTANCE = }. Current value ({self.operations['INDUCTANCE_MODE'] = }) is not allowed. Please check {self.workbook_sheet_name[2]} in file {self.workbook_name}.\n"
             )
 
         inductance_switch = {
+            CONSTANT_INDUCTANCE: self.__constant_inductance,
             ANALYTICAL_INDUCTANCE: self.__inductance_analytical_calculation,
             APPROXIMATE_INDUCTANCE: self.__inductance_approximate_calculation,
         }
 
-        inductance_switch[self.operations["INDUCTANCE_MODE"]]()
+        inductance_switch[self.operations["INDUCTANCE_MODE"]](self.operations["SELF_INDUCTANCE_MODE"])
 
         self.electric_mass_matrix[
             : self.total_elements_current_carriers,
@@ -4234,6 +4321,12 @@ def electric_solution_reorganization(self):
                 ii :: self.inventory["StrandComponent"].number
             ]
 
+            obj.dict_node_pt["current_along"] = np.interp(
+                        self.grid_features["zcoord"],
+                        self.grid_features["zcoord_gauss"],
+                        obj.dict_Gauss_pt["current_along"]
+                        )
+
             obj.dict_Gauss_pt["delta_voltage_along"] = delta_voltage_along[
                 ii :: self.inventory["StrandComponent"].number
             ]
@@ -4368,8 +4461,7 @@ def __electric_method_steady(self):
         # conductor_initialization);
         # 2) all the Joule power contribution due to the electric solution 
         # are evaluated calling method self.__get_heat_source_em_steady.
-        # N.B. arrays obj.dict_Gauss_pt["integral_power_el_res_mod1"]
-        # obj.dict_Gauss_pt["integral_power_el_res_mod2"] and 
+        # N.B. arrays obj.dict_Gauss_pt["integral_power_el_res"] and 
         # obj.dict_node_pt["integral_power_el_cond"] are set to zero for 
         # the next evaluation inside method __build_heat_source_gauss_pt
         self.__build_heat_source_gauss_pt()
@@ -4627,8 +4719,7 @@ def __build_heat_source_gauss_pt(self):
         """Private method that builds heat source therms in Gauss points for 
         strand and jacket objects.
         Sets to zeros the following arrays as preliminary step for the next evaluation:
-            * strand.dict_Gauss_pt["integral_power_el_res_mod1"]
-            * strand.dict_Gauss_pt["integral_power_el_res_mod2"]
+            * strand.dict_Gauss_pt["integral_power_el_res"]
             * strand.dict_node_pt["integral_power_el_cond"]
         """
 
@@ -4649,8 +4740,7 @@ def __build_heat_source_gauss_pt(self):
                 + strand.dict_Gauss_pt["linear_power_el_resistance"]
             )
 
-            # Set arrays strand.dict_Gauss_pt["integral_power_el_res_mod1"], 
-            # strand.dict_Gauss_pt["integral_power_el_res_mod2"] and 
+            # Set arrays strand.dict_Gauss_pt["integral_power_el_res"] and 
             # strand.dict_node_pt["integral_power_el_cond"] to zero for the 
             # next evaluation.
             strand.set_power_array_to_zeros(self)
@@ -4700,13 +4790,12 @@ def __build_heat_source_gauss_pt(self):
     
     def eval_integral_joule_power(self):
         
-        """Method that evaluates the numerator of the expression used to evaluate the integral value of the Joule power. For the Joule power due to the electric resistance along the current carriers two different but equivalent approaches are used:
-            * mode1 -> P_Joule = R * I^2 (R electric resistance in Ohm, I electric current)
-            * mode2 -> P_Joule = \Delta_Phi * I
+        """Method that evaluates the numerator of the expression used to evaluate the integral value of the Joule power. The Joule power along the current carriers is computed as:
+            P_Joule = \Delta_Phi * I
         with
-            * R electric resistance in Ohm
             * I electric current in A
             * \Delta_Ph electric voltage potential difference along current carriers in V
+        This approach was discussed with prof. Zach Hartwig and Dr. Nicolò Riva and is more general and coservative that the evaluation that takes into account only the contribution of the electric resistance.
         For the Joule power due to the electric conductances between current carriers it is exploited the power computed in method get_total_joule_power_electric_conductance
         Regardless of the kind of Joule power, the integration can be performed as follows:
             P_Joule = 1/Delta_t_TH * int_0^Delta_t_TH (dt P_{Joule,i})
@@ -4730,27 +4819,17 @@ def eval_integral_joule_power(self):
         # Loop on StrandComponent objects.
         for strand in self.inventory["StrandComponent"].collection:
             
-            # Compute the numerator of the integral Joule power due to the 
-            # electric resistance along the current carrier, mode1.
-            strand.dict_Gauss_pt["integral_power_el_res_mod1"] += (
-                strand.dict_Gauss_pt["current_along"] ** 2
-                * strand.dict_Gauss_pt["electric_resistance"]
-                * self.electric_time_step
-            )
-
-            # Compute the numerator of the integral Joule power due to the 
-            # electric resistance along the current carrier, mode2.
-            strand.dict_Gauss_pt["integral_power_el_res_mod2"] += (
+            # Compute the numerator of the integral Joule power along the 
+            # current carrier.
+            # N.B. this evaluation accounts aslo for the voltage due to the 
+            # inductance and is a conservative an more general approach.
+            # Discussed with prof. Zach Hartwig and Dr. Nicolò Riva.
+            strand.dict_Gauss_pt["integral_power_el_res"] += (
                 strand.dict_Gauss_pt["current_along"]
                 * strand.dict_Gauss_pt["delta_voltage_along"]
                 * self.electric_time_step
             )
 
-            # Check equivalence of Mode 1 and Mode 2 (they should be equivalent 
-            # but Mode 2 may give numerical cancellation).
-            if not np.allclose(strand.dict_Gauss_pt["integral_power_el_res_mod1"],strand.dict_Gauss_pt["integral_power_el_res_mod2"]):
-                warnings.warn("P_Joule = R I^2 dt_em != Delta_Phi I dt_em.Possible violation of the energy conservation!")
-
             # Compute the numerator of the integral Joule power due to the 
             # electric conductance between current carriers.
             strand.dict_node_pt["integral_power_el_cond"] += (
@@ -4812,6 +4891,10 @@ def operating_conditions_em(self):
                     strand.dict_node_pt["electrical_resistivity_stabilizer"] = strand.strand_electrical_resistivity_not_sc(
                             strand.dict_node_pt
                         )
+                elif isinstance(strand, StackComponent):
+                    strand.dict_Gauss_pt["electrical_resistivity_stabilizer"] = strand.stack_electrical_resistivity_not_sc(
+                        strand.dict_Gauss_pt
+                        )
                 elif isinstance(strand, StrandStabilizerComponent):
                     strand.dict_node_pt["electrical_resistivity_stabilizer"] = strand.strand_electrical_resistivity(
                             strand.dict_node_pt
@@ -4900,7 +4983,11 @@ def __eval_gauss_point_em(self):
                 # electric time step.
                 if isinstance(strand,StrandMixedComponent):
                     strand.dict_Gauss_pt["electrical_resistivity_stabilizer"] = strand.strand_electrical_resistivity_not_sc(
-                            strand.dict_Gauss_pt
+                        strand.dict_Gauss_pt
+                        )
+                elif isinstance(strand, StackComponent):
+                    strand.dict_Gauss_pt["electrical_resistivity_stabilizer"] = strand.stack_electrical_resistivity_not_sc(
+                        strand.dict_Gauss_pt
                         )
                 elif isinstance(strand, StrandStabilizerComponent):
                     strand.dict_Gauss_pt["electrical_resistivity_stabilizer"] = strand.strand_electrical_resistivity(
diff --git a/source_code/conductor_flags.py b/source_code/conductor_flags.py
index 7c61647..4d87e10 100644
--- a/source_code/conductor_flags.py
+++ b/source_code/conductor_flags.py
@@ -16,15 +16,19 @@
 # Electric conductance is not defined per unit length
 ELECTRIC_CONDUCTANCE_NOT_UNIT_LENGTH = 2
 
+# Self inductance readed from input file conductor_definition.xlsx
+SELF_INDUCTANCE_MODE_0 = 0
 # Analytical self inductance is evaluated according to mode 1.
 SELF_INDUCTANCE_MODE_1 = 1
 # Analytical self inductance is evaluated according to mode 2.
 SELF_INDUCTANCE_MODE_2 = 2
 
+# Constant mutual inductance readed from input file conductor_definition.xlsx
+CONSTANT_INDUCTANCE = 0
 # Flag to evaluate inductance analytically
-ANALYTICAL_INDUCTANCE = 0
+ANALYTICAL_INDUCTANCE = 1
 # Flag to evaluate inductance using an approximation
-APPROXIMATE_INDUCTANCE = 1
+APPROXIMATE_INDUCTANCE = 2
 
 # Flag to solve the electric problem in steady state conditions.
 STATIC_ELECTRIC_SOLVER = 0
diff --git a/source_code/input_files/input_file_template/template_conductor_definition.xlsx b/source_code/input_files/input_file_template/template_conductor_definition.xlsx
index fbe93ab..2f9a33f 100644
Binary files a/source_code/input_files/input_file_template/template_conductor_definition.xlsx and b/source_code/input_files/input_file_template/template_conductor_definition.xlsx differ
diff --git a/source_code/properties_of_materials/rare_earth_123.py b/source_code/properties_of_materials/rare_earth_123.py
index 9722c3b..21ccf45 100644
--- a/source_code/properties_of_materials/rare_earth_123.py
+++ b/source_code/properties_of_materials/rare_earth_123.py
@@ -393,7 +393,7 @@ def critical_current_density_bisection_re123(TT, BB, JOP, TC0M, BC20M, C):
 
     for _, vv in enumerate(JC_ind):
 
-        T_lower = np.array([3.5])
+        T_lower = np.array([0])
         # T_upper = np.array([40.0])
         T_upper = np.array([TC0M])
 
diff --git a/source_code/solid_component.py b/source_code/solid_component.py
index e3d98f1..f33a086 100644
--- a/source_code/solid_component.py
+++ b/source_code/solid_component.py
@@ -267,6 +267,10 @@ def get_current(self, conductor: object):
                     self.dict_Gauss_pt["op_current_sc"] = self.dict_Gauss_pt[
                         "op_current"
                     ]
+                if conductor.cond_time[-1] == 0:
+                    self.dict_node_pt["current_along"] = self.dict_node_pt["op_current"]
+                    self.dict_Gauss_pt["current_along"] = self.dict_Gauss_pt["op_current"]
+
 
                 if self.flagSpecfield_current == 2:
                     # Add also a logger
@@ -297,6 +301,9 @@ def get_current(self, conductor: object):
                         self.dict_node_pt["op_current_sc"][:-1]
                         + self.dict_node_pt["op_current_sc"][1:]
                     ) / 2.0
+                if conductor.cond_time[-1] == 0:
+                    self.dict_node_pt["current_along"] = self.dict_node_pt["op_current"]
+                    self.dict_Gauss_pt["current_along"] = self.dict_Gauss_pt["op_current"]
 
             elif conductor.inputs["I0_OP_MODE"] == IOP_NOT_DEFINED:
                 # User does not specify a current: set current carrient 
@@ -659,8 +666,7 @@ def initialize_electric_quantities(self, conductor):
         * self.dict_Gauss_pt["current_along"]
         * self.dict_Gauss_pt["delta_voltage_along"]
         * self.dict_Gauss_pt["delta_voltage_along_sum"]
-        * self.dict_Gauss_pt["integral_power_el_res_mod1"]
-        * self.dict_Gauss_pt["integral_power_el_res_mod2"]
+        * self.dict_Gauss_pt["integral_power_el_res"]
         * self.dict_Gauss_pt["linear_power_el_resistance"]
         * self.dict_node_pt["total_linear_power_el_cond"]
         * self.dict_node_pt["total_power_el_cond"]
@@ -675,12 +681,10 @@ def initialize_electric_quantities(self, conductor):
         self.dict_Gauss_pt["current_along"] = np.zeros(n_elem)
         self.dict_Gauss_pt["delta_voltage_along"] = np.zeros(n_elem)
         self.dict_Gauss_pt["delta_voltage_along_sum"] = np.zeros(n_elem)
-        # Introduced to evaluate the integral value Joule power due to electric 
-        # resistance along the cable in the electric module as R * I^2.
-        self.dict_Gauss_pt["integral_power_el_res_mod1"] = np.zeros(n_elem)
+
         # Introduced to evaluate the integral value Joule power due to electric 
         # resistance along the cable in the electric module as Delta V * I.
-        self.dict_Gauss_pt["integral_power_el_res_mod2"] = np.zeros(n_elem)
+        self.dict_Gauss_pt["integral_power_el_res"] = np.zeros(n_elem)
 
         self.dict_node_pt["total_power_el_cond"] = np.zeros(n_nod)
         # Introduced to evaluate the integral value Joule power due to electric 
@@ -716,12 +720,11 @@ def set_power_array_to_zeros(self, conductor):
         n_elem = conductor.grid_input["NELEMS"]
 
         # Set the following arrays to zeros
-        self.dict_Gauss_pt["integral_power_el_res_mod1"] = np.zeros(n_elem)
-        self.dict_Gauss_pt["integral_power_el_res_mod2"] = np.zeros(n_elem)
+        self.dict_Gauss_pt["integral_power_el_res"] = np.zeros(n_elem)
         self.dict_node_pt["integral_power_el_cond"] = np.zeros(n_nod)
 
     def get_joule_power_along(self, conductor: object):
-        """Method that evaluate the contribution to the integral power in the element of Joule power (in W/m) due to the electic resistances along the SolidComponent objects.
+        """Method that evaluate the contribution to the integral power in the element of Joule power (in W/m) due to the voltage difference (due to electric resistance and inductance) along the SolidComponent objects. This approach was discussed with prof. Zach Hartwig and Dr. Nicolò Riva and is more general and coservative that the evaluation that takes into account only the contribution of the electric resistance.
         This method should be called in the electric method, when the transient solution is used and only for current carriers.
 
         Args:
@@ -733,14 +736,14 @@ def get_joule_power_along(self, conductor: object):
 
         # Finalize the evaluation of the integral value of the Joule power due 
         # to electric resistance along the current carrier started in method 
-        # conductor.eval_integral_joule_power. Array integral_power_el_res_mod1 
+        # conductor.eval_integral_joule_power. Array integral_power_el_res 
         # stores the numerator (energy in J) that here is divided by the 
         # thermal hydraulic to get again a power (W), which is further divided 
         # by the length of the discreticazion element corrected by cos(theta) 
         # to get a linear power density (W/m):
         # P_Joule = sum_1^N_em P_{Joule,i} / (Delta_t_TH * Delta_z * costheta)
         integral_j_pow_along = (
-            self.dict_Gauss_pt["integral_power_el_res_mod1"]
+            self.dict_Gauss_pt["integral_power_el_res"]
             / (conductor.grid_features["delta_z"] * self.inputs["COSTETA"]
             * conductor.time_step)
         )
@@ -823,7 +826,8 @@ def get_joule_power_across(self, conductor: object):
             )
 
     def get_joule_power_along_steady(self, conductor: object):
-        """Method that evaluate the contribution to the total power in the element of Joule power (in W/m) due to the electic resistances along the SolidComponent objects.
+        """Method that evaluate the contribution to the total power in the element of Joule power (in W/m) due to the voltage difference (due to electric resistance and inductance) along the SolidComponent objects.
+        This approach was discussed with prof. Zach Hartwig and Dr. Nicolò Riva and is more general and coservative that the evaluation that takes into account only the contribution of the electric resistance.
         This method should be called in the electric method, when the transient solution is used and for the solid component that actually carries a current. It works when the steady state solution for the electric module is computed.
 
         Args:
@@ -831,31 +835,24 @@ def get_joule_power_along_steady(self, conductor: object):
         """
 
         # Alias
-        el_res = self.dict_Gauss_pt["electric_resistance"]
         current = self.dict_Gauss_pt["current_along"]
         voltage = self.dict_Gauss_pt["delta_voltage_along"]
         d_z_tilde = conductor.grid_features["delta_z"] * self.inputs["COSTETA"]
 
-        # Mode 1: evaluate Joule linear power along the strand in W, due
-        # to electric resistances only for current carriers:
-        # P_along = R_along * I_along ^2
-        self.dict_Gauss_pt["integral_power_el_res_mod1"] = current ** 2 * el_res
-
-        # Mode 2: evaluate Joule linear power along the strand in W, due
+        # Evaluate Joule linear power along the strand in W, due
         # to electric resistances only for current carriers:
         # P_along = Delta_Phi_along * I_along
-        self.dict_Gauss_pt["integral_power_el_res_mod2"] = voltage * current
+        # N.B. this evaluation accounts aslo for the voltage due to the 
+        # inductance and is a conservative an more general approach. Discussed 
+        # with prof. Zach Hartwig and Dr. Nicolò Riva.
+        self.dict_Gauss_pt["integral_power_el_res"] = voltage * current
 
-        # Check equivalence of Mode 1 and Mode 2 (they should be equivalent 
-        # but Mode 2 may give numerical cancellation).
-        if not np.allclose(self.dict_Gauss_pt["integral_power_el_res_mod1"],self.dict_Gauss_pt["integral_power_el_res_mod2"]):
-            warnings.warn("P_Joule = R I^2 != Delta_Phi I. Possible violation of the energy conservation!")
         # Convert W in W/m keping into account the cos(theta).
         # This is independent of the time integration method since at the 
         # initialization (when this function is used) BE, CN and AM4 should 
         # fill only the 0 index column.
         self.dict_Gauss_pt["linear_power_el_resistance"][:, 0] = (
-            self.dict_Gauss_pt["integral_power_el_res_mod1"] / d_z_tilde
+            self.dict_Gauss_pt["integral_power_el_res"] / d_z_tilde
         )
 
     def get_joule_power_across_steady(self, conductor: object):
diff --git a/source_code/stack_component.py b/source_code/stack_component.py
index 4de5212..1913b16 100644
--- a/source_code/stack_component.py
+++ b/source_code/stack_component.py
@@ -78,7 +78,7 @@
     cu=density_cu,
     ge=density_ge,
     hc276=density_hc276,
-    re123=density_re123,
+    ybco=density_re123,
     sn60pb40=density_sn60pb40,
     ss=density_ss,
 )
@@ -89,7 +89,7 @@
     cu=thermal_conductivity_cu_nist,
     ge=thermal_conductivity_ge,
     hc276=thermal_conductivity_hc276,
-    re123=thermal_conductivity_re123,
+    ybco=thermal_conductivity_re123,
     sn60pb40=thermal_conductivity_sn60pb40,
     ss=thermal_conductivity_ss,
 )
@@ -100,7 +100,7 @@
     cu=isobaric_specific_heat_cu_nist,
     ge=isobaric_specific_heat_ge,
     hc276=isobaric_specific_heat_hc276,
-    re123=isobaric_specific_heat_re123,
+    ybco=isobaric_specific_heat_re123,
     sn60pb40=isobaric_specific_heat_sn60pb40,
     ss=isobaric_specific_heat_ss,
 )
@@ -115,6 +115,8 @@
     ss=electrical_resistivity_ss,
 )
 
+INGEGNERISTIC_MODE = 0
+PHYSICAL_MODE = 1
 
 class StackComponent(StrandComponent):
     """Class that defines StackComponents objects to model HTS stacks of tapes.
@@ -167,6 +169,7 @@ def __init__(
             current_along=dict(),
             delta_voltage_along=dict(),
             linear_power_el_resistance=dict(),
+            delta_voltage_along_sum=dict(),
         )
         self.dict_scaling_input = dict()
         # Dictionary initialization: inputs.
@@ -190,10 +193,16 @@ def __init__(
 
         # Check that costheta is in the range (0,1].
         check_costheta(self,dict_file_path["input"],sheet)
+        self.__compute_cross_section(sheet, dict_file_path["input"])
+        self.__get_current_density_cross_section(sheet, dict_file_path["input"])
 
         # Call SolidComponent class constructor to deal with StrandMixedComponent time \
         # steps for current, external heating and so on
+        
         SolidComponent(simulation, self)
+        if self.inputs["Stabilizer_material"] != "Cu":
+            # remove key RRR from inputs if stabilizer is not Cu (cdp, 07/2020)
+            self.inputs.pop("RRR")
         if self.operations["IBIFUN"] != -1:
             # Remove key B_field_units.
             del self.operations["B_field_units"]
@@ -203,39 +212,34 @@ def __init__(
         # of flag self.operations["IOP_MODE"].
         self.deal_with_flag_IOP_MODE()
 
-        self.__reorganize_input()
-        self.__check_consistency(conductor)
-
         # Flag to check if evaluation of homogenized isobaric specific heat can
         # be done or not (depends on homogenized density evaluation).
         self.__stack_density_flag = False
+        self.__reorganize_input()
+        self.__check_consistency(conductor)
 
         # Call method deal_with_fixed_potential to manipulate input about fixed
         # potential values.
         self.deal_with_fixed_potential(conductor.inputs["ZLENGTH"])
+       
+    def __repr__(self):
+        return f"{self.__class__.__name__}(Type: {self.name}, identifier: {self.identifier})"
 
+    def __str__(self):
+        pass
+    
+    def __compute_cross_section(self,sheet, file_path:str):
+        self.cross_section = dict()
+        self.cross_section["total"] = self.inputs["CROSSECTION"]
         # Superconductor total cross section in m^2
-        self.sc = (
+        self.cross_section["sc"] = (
             self.inputs["HTS_thickness"]
             * self.inputs["Stack_width"]
-            * self.inputs["Tape_number"]
+            * self.inputs["N_tape"]
         )
+        self.cross_section["sc_sloped"] = self.cross_section["sc"] / self.inputs["COSTETA"]
+        
         # Stabilizer (not sc) total cross section in m^2
-        self.stabilizer_cross_section = (
-            self.tape_thickness_not_sc
-            * self.inputs["Stack_width"]
-            * self.inputs["Tape_number"]
-        )
-
-    def __repr__(self):
-        return f"{self.__class__.__name__}(Type: {self.name}, identifier: {self.identifier})"
-
-    def __str__(self):
-        pass
-
-    def __reorganize_input(self):
-        """Private method that reorganizes input data stored in dictionary self.inputs to simplify the procedure of properties homogenization."""
-
         # Create numpy array of string with the identifier of tape material
         self.tape_material = np.array(
             [
@@ -258,7 +262,7 @@ def __reorganize_input(self):
                 value.lower()
                 for key, value in self.inputs.items()
                 if key.endswith("material")
-                and key != "HTS_material"
+                and key != "superconducting_material"
                 and value != "none"
             ],
             dtype=str,
@@ -291,6 +295,51 @@ def __reorganize_input(self):
         )
         # Total not superconducting tape thickness in m.
         self.tape_thickness_not_sc = self.material_thickness_not_sc.sum()
+        
+        self.cross_section["stab"] = (
+            self.tape_thickness_not_sc
+            * self.inputs["Stack_width"]
+            * self.inputs["N_tape"]
+        )
+
+        self.cross_section["stab_sloped"] = self.cross_section["stab"] / self.inputs["COSTETA"]
+
+        # Evaluate alpha_0 (stabilizer non stabilizer ratio defined wrt the 
+        # not segregated stabilizer cross section only).
+        self.cross_section["alpha_0"] = self.cross_section["stab"] / self.cross_section["sc"]
+
+
+    def __get_current_density_cross_section(self, sheet, file_path:str):
+        """Private method that evalutates cross section used to compute the current density and the current sharing temperature according to the definiton of the critical current density scaling parameter c0.
+
+        Args:
+            sheet (Chartsheet | Worksheet | ReadOnlyWorksheet): sheet with input data for StrandMixedComponent definition.
+            file_path (str): path to the input file.
+
+        Raises:
+            ValueError: if self.inputs["C0_MODE"] =! 0 or self.inputs["C0_MODE"] != 1
+        """
+                
+        if self.inputs["C0_MODE"] == INGEGNERISTIC_MODE:
+            # Ingegneristic definition for c0 is used, i.e., the value is 
+            # normalized with respect to the total perpendicular cross section 
+            # of superconducting strand.
+
+            # Convert the ingegneristic definition to the physical definition 
+            # (normalized with respect to the total perpendicular cross section 
+            # of the superconductor) in order to use the superconducting cross 
+            # section in the electric model.
+            self.inputs["c0"] = self.inputs["c0"] * (1. + self.cross_section["alpha_0"])
+        elif (self.inputs["C0_MODE"] != PHYSICAL_MODE
+                and self.inputs["C0_MODE"] != INGEGNERISTIC_MODE
+            ):
+            raise ValueError(
+                f"Not valid value for flag self.inputs['C0_MODE']. Flag value should be {PHYSICAL_MODE = } or {INGEGNERISTIC_MODE = }, current value is {self.inputs['C0_MODE'] = }.\nPlease check in sheet {sheet} of input file {file_path}.\n"
+            )
+
+
+    def __reorganize_input(self):
+        """Private method that reorganizes input data stored in dictionary self.inputs to simplify the procedure of properties homogenization."""
 
         # Create numpy array with density functions according to the tape
         # material; order is consistent with values in self.tape_material.
@@ -436,15 +485,14 @@ def stack_isobaric_specific_heat(self, property: dict) -> np.ndarray:
         # Set flag to false to trigger error in the next homogenized isobaric
         # specific heat evaluation if not done properly.
         self.__stack_density_flag = False
-        isobaric_specific_heat = np.array(
-            [
-                func(property["temperature"])
-                for func in self.isobaric_specific_heat_function
-            ]
+        isobaric_specific_heat = np.zeros(
+            (property["temperature"].size, self.inputs["Material_number"])
         )
+        for ii, func in enumerate(self.isobaric_specific_heat_function):
+            isobaric_specific_heat[:, ii] = func(property["temperature"])
         # Evaluate homogenized isobaric specific heat of the stack:
         # cp_eq = sum(s_i*rho_i*cp_i)/sum(s_i*rho_i)
-        return (isobaric_specific_heat.T * self.__density_numerator).sum(
+        return (isobaric_specific_heat * self.__density_numerator).sum(
             axis=1
         ) / self.__density_numerator_sum
 
@@ -498,10 +546,11 @@ def stack_electrical_resistivity_not_sc(self, property: dict) -> np.ndarray:
                 electrical_resistivity[:, ii] = func(property["temperature"])
         if self.inputs["Material_number"] - 1 > 1:
             # Evaluate homogenized electrical resistivity of the stack:
-            # rho_el_eq = s_not_sc * (sum(s_i/rho_el_i))^-1 for any i not sc
-            return self.tape_thickness_not_sc * np.reciprocal(
-                (self.material_thickness_not_sc / electrical_resistivity).sum(axis=1)
-            )
+            # rho_el_eq = A_not_sc * (sum(A_i/rho_el_i))^-1 for any i not sc
+            # In this case becomes rho_el_eq = t_not_sc * (sum(t_i/rho_el_i))^-1 where t is the thickness
+            
+            return self.tape_thickness_not_sc * np.reciprocal((
+                self.material_thickness_not_sc / electrical_resistivity).sum(axis=1))
         elif self.inputs["Material_number"] - 1 == 1:
             return electrical_resistivity.reshape(property["temperature"].size)
 
@@ -575,14 +624,22 @@ def solve_current_divider(
             raise ValueError(
                 f"Arrays rho_el_stabilizer and critical_current must have the same shape.\n {rho_el_stabilizer.shape = };\n{critical_current.shape}.\n"
             )
-
+        if rho_el_stabilizer.shape != current.shape:
+            raise ValueError(
+                f"Arrays rho_el_stabilizer and current must have the same shape.\n {rho_el_stabilizer.shape = };\n{current.shape}.\n"
+            )
+        if critical_current.shape != current.shape:
+            raise ValueError(
+                f"Arrays critical_current and current must have the same shape.\n {critical_current.shape = };\n{current.shape}.\n"
+            )
+        
         # Evaluate constant value:
         # psi = rho_el_stab*I_c^n/(E_0*A_stab)
         psi = (
             rho_el_stabilizer
             * critical_current ** self.inputs["nn"]
             / self.inputs["E0"]
-            / self.stabilizer_cross_section
+            / self.cross_section["stab"]
         )
 
         # Initialize guess.
@@ -600,6 +657,19 @@ def solve_current_divider(
                 disp=False,
             )
 
+        # Tolerance on newton halley increased in case I_critical is very small,
+        # to avoid inaccuracies on the divider that could lead to voltage 
+        # differences between sc and stab that are not expected in the parallel
+        # of electric resistances.
+        if min(critical_current)>1e-6:
+            # Default tollerance in optimize.newton method
+            tollerance = 1.48e-8
+        else:
+            # Value found trial and error iteration
+            # tollerance = 1e-12
+            # Other possible solution for the correct tollerance
+            tollerance = min(critical_current)/1e3
+        
         # Evaluate superconducting with Halley's method
         sc_current = optimize.newton(
             self.__sc_current_residual,
@@ -607,6 +677,7 @@ def solve_current_divider(
             args=(psi, current),
             fprime=self.__d_sc_current_residual,
             fprime2=self.__d2_sc_current_residual,
+            tol = tollerance,
             maxiter=1000,
         )
 
@@ -714,174 +785,81 @@ def get_electric_resistance(self, conductor: object) -> np.ndarray:
             np.ndarray: array of electrical resistance in Ohm of length {conductor.grid_input["NELEMS"] = }.
         """
 
-        critical_current_node = self.sc_cross_section * self.dict_node_pt["J_critical"]
         critical_current_gauss = (
-            self.sc_cross_section * self.dict_Gauss_pt["J_critical"]
+            self.cross_section["sc"] * self.dict_Gauss_pt["J_critical"]
         )
 
-        # Get index that correspond to superconducting regime.
-        ind_sc_node = np.nonzero(
-            self.dict_node_pt["op_current_sc"] / critical_current_node < 0.95
-        )[0]
-        ind_sc_gauss = np.nonzero(
-            self.dict_Gauss_pt["op_current_sc"] / critical_current_gauss < 0.95
-        )[0]
-        # Get index that correspond to current sharing regime.
-        ind_sh_node = np.nonzero(
-            self.dict_node_pt["op_current_sc"] / critical_current_node >= 0.95
-        )[0]
-        ind_sh_gauss = np.nonzero(
-            self.dict_Gauss_pt["op_current_sc"] / critical_current_gauss >= 0.95
-        )[0]
-
-        # Initialize electric resistance arrays in both nodal and Gauss points;
-        # this is the equivalent electrical resistance, thus it is defined in
-        # this way:
-        # R_eq = R_sc if superconducting regime
-        # R_eq = R_sc * R_stab/(R_sc + R_stab) is sharing or normal regime
-        self.dict_node_pt["electric_resistance"] = 10.0 * np.ones(
-            self.dict_node_pt["temperature"].shape
-        )
-        self.dict_Gauss_pt["electric_resistance"] = 10.0 * np.ones(
-            self.dict_Gauss_pt["temperature"].shape
-        )
-
-        ## SUPERCONDUCTING REGIME ##
-
-        # Strand current in superconducting regime is the one carriend by the
-        # superconducting material only.
-        self.dict_node_pt["op_current"][ind_sc_node] = self.dict_node_pt[
-            "op_current_sc"
-        ][ind_sc_node]
-        self.dict_Gauss_pt["op_current"][ind_sc_gauss] = self.dict_Gauss_pt[
-            "op_current_sc"
-        ][ind_sc_gauss]
-
-        # Initialize array of superconducting electrical resistivit in nodal and Gauss points to None.
-        self.dict_node_pt["electrical_resistivity_superconductor"] = np.full_like(
-            self.dict_node_pt["temperature"], None
-        )
-        self.dict_Gauss_pt["electrical_resistivity_superconductor"] = np.full_like(
-            self.dict_Gauss_pt["temperature"], None
-        )
-
-        # Compute superconducting electrical resistivity only in index for
-        # which the superconducting regime is guaranteed, using the power low.
-        self.dict_node_pt["electrical_resistivity_superconductor"][
-            ind_sc_node
-        ] = self.superconductor_power_law(
-            self.dict_node_pt["op_current_sc"][ind_sc_node],
-            critical_current_node[ind_sc_node],
-            self.dict_node_pt["J_critical"][ind_sc_node],
-        )
-        self.dict_Gauss_pt["electrical_resistivity_superconductor"][
-            ind_sc_gauss
-        ] = self.superconductor_power_law(
-            self.dict_Gauss_pt["op_current_sc"][ind_sc_gauss],
-            critical_current_gauss[ind_sc_gauss],
-            self.dict_Gauss_pt["J_critical"][ind_sc_gauss],
-        )
-        # Evaluate electic resistance in superconducting region (superconductor
-        # only).
-        self.dict_Gauss_pt["electric_resistance"][
-            ind_sc_gauss
-        ] = self.electric_resistance(
-            conductor, "electrical_resistivity_superconductor", ind_sc_gauss
-        )
-
-        ## SHARING OR NORMAL REGIME ##
+        # Make initialization only once for each conductor object.
+        if conductor.cond_num_step == 0:
+            # Initialize electric resistance arrays in Gauss point; this is the 
+            # equivalent electrical resistance, thus it is defined in this way:
+            # R_eq = R_sc if superconducting regime
+            # R_eq = R_sc * R_stab/(R_sc + R_stab) is normal regime.
+            self.dict_Gauss_pt["electric_resistance"] = np.full_like(
+                self.dict_Gauss_pt["temperature"], None
+            )
 
-        # Strand current in sharing regime is the one carried by the both the
-        # superconducting and the stabilizer materials.
-        # self.dict_node_pt["op_current"][ind_sh_node] = self.dict_node_pt["op_current"][ind_sh_node]
-        # self.dict_Gauss_pt["op_current"][ind_sh_node] = self.dict_node_pt["op_current"][ind_sh_gauss]
+            # Initialize array of superconducting electrical resistivit in 
+            # Gauss point to None.
+            self.dict_Gauss_pt["electrical_resistivity_superconductor"] = np.full_like(
+                self.dict_Gauss_pt["temperature"], None
+            )
 
-        # Evaluate how the current is distributed solving the current divider
-        # problem in both nodal and Gauss points.
-        sc_current_node, stab_current_node = self.solve_current_divider(
-            self.dict_node_pt["electrical_resistivity_stabilizer"][ind_sh_node],
-            critical_current_node[ind_sh_node],
-            self.dict_node_pt["op_current"][ind_sh_node],
-        )
-        sc_current_gauss, stab_current_gauss = self.solve_current_divider(
-            self.dict_Gauss_pt["electrical_resistivity_stabilizer"][ind_sh_gauss],
-            critical_current_gauss[ind_sh_gauss],
-            self.dict_Gauss_pt["op_current"][ind_sh_gauss],
-        )
+        # Get index for which abs(critical_current_gauss) == 0 (inside normal 
+        # zone by definition).
+        ind_zero = np.nonzero(abs(critical_current_gauss) == 0)[0]
+        # Get index for which abs(critical_current_gauss) > 0 (outside normal 
+        # zone by definition).
+        ind_not_zero = np.nonzero(abs(critical_current_gauss) > 0)[0]
 
-        # Get index of the normal region, to avoid division by 0 in evaluation
-        # of sc electrical resistivity with the power law.
-        ind_normal_node = np.nonzero(
-            (
-                stab_current_node / self.dict_node_pt["op_current"][ind_sh_node]
-                > 0.999999
-            )
-            | (sc_current_node < 1.0)
-        )[0]
-        ind_normal_gauss = np.nonzero(
-            (
-                stab_current_gauss / self.dict_Gauss_pt["op_current"][ind_sh_gauss]
-                > 0.999999
-            )
-            | (sc_current_gauss < 1.0)
-        )[0]
-
-        ## NORMAL REGIME ONLY ##
-        if ind_normal_node.any():
-            # Get the index of location of true current sharing region;
-            # overwrite ind_sh_node.
-            ind_sh_node = np.nonzero(
-                (
-                    stab_current_node / self.dict_node_pt["op_current"][ind_sh_node]
-                    <= 0.999999
-                )
-                | (sc_current_node >= 1.0)
-            )[0]
-        if ind_normal_gauss.any():
-            # Get the index of location of true current sharing region;
-            # overwrite ind_sh_gauss.
-            ind_sh_gauss = np.nonzero(
-                (
-                    stab_current_gauss / self.dict_Gauss_pt["op_current"][ind_sh_gauss]
-                    <= 0.999999
-                )
-                | (sc_current_gauss >= 1.0)
-            )[0]
+        # Check if np array ind_zero is not empty: NORMAL REGION BY DEFINITION
+        if ind_zero.any():
             # Evaluate electic resistance in normal region (stabilizer only).
             self.dict_Gauss_pt["electric_resistance"][
-                ind_normal_gauss
+                ind_zero
             ] = self.electric_resistance(
-                conductor, "electrical_resistivity_stabilizer", ind_normal_gauss
+                conductor, "electrical_resistivity_stabilizer", "stab", ind_zero
             )
 
-        ## SHARING REGIME ONLY ##
-        # Evaluate the electrical resistivity of the superconductor according
-        # to the power low in both nodal and Gauss points in Ohm*m.
-        self.dict_node_pt["electrical_resistivity_superconductor"][
-            ind_sh_node
-        ] = self.superconductor_power_law(
-            sc_current_node[ind_sh_node],
-            critical_current_node[ind_sh_node],
-            self.dict_node_pt["J_critical"][ind_sh_node],
-        )
-        self.dict_Gauss_pt["electrical_resistivity_superconductor"][
-            ind_sh_gauss
-        ] = self.superconductor_power_law(
-            sc_current_gauss[ind_sh_gauss],
-            critical_current_gauss[ind_sh_gauss],
-            self.dict_Gauss_pt["J_critical"][ind_sh_gauss],
-        )
+        # Check if np array ind_not_zero is not empty: deal with index that 
+        # are outside normal zone by definition; however some of them could 
+        # still identify a normal region.
+        if ind_not_zero.any():
 
-        # Evaluate the equivalent electric resistance in Ohm.
-        self.dict_Gauss_pt["electric_resistance"][
-            ind_sh_gauss
-        ] = self.parallel_electric_resistance(
-            conductor,
-            [
-                "electrical_resistivity_superconductor",
-                "electrical_resistivity_stabilizer",
-            ],
-            ind_sh_gauss,
-        )
+            sc_current_gauss, stab_current_gauss = self.solve_current_divider(
+                self.dict_Gauss_pt["electrical_resistivity_stabilizer"][ind_not_zero],
+                critical_current_gauss[ind_not_zero],
+                self.dict_Gauss_pt["op_current"][ind_not_zero]
+                )
+            
+            self.dict_Gauss_pt["electrical_resistivity_superconductor"][
+                ind_not_zero
+            ] = self.superconductor_power_law(
+                sc_current_gauss,
+                critical_current_gauss[ind_not_zero],
+                self.dict_Gauss_pt["J_critical"][ind_not_zero],
+                )
+            # Evaluate the equivalent electric resistance in Ohm.
+            self.dict_Gauss_pt["electric_resistance"][
+                ind_not_zero] = self.parallel_electric_resistance(
+                conductor,
+                [
+                    "electrical_resistivity_superconductor",
+                    "electrical_resistivity_stabilizer",
+                ],["sc","stab"],
+                ind_not_zero,
+            )
 
-        return self.dict_Gauss_pt["electric_resistance"]
+            # Compute voltage along stabilizer.
+            v_stab = self.dict_Gauss_pt["electrical_resistivity_stabilizer"][
+                ind_not_zero
+            ] * stab_current_gauss / self.cross_section["stab"]
+            # Compute voltage along superconductor
+            v_sc = self.inputs["E0"] * (sc_current_gauss / critical_current_gauss[ind_not_zero]) ** self.inputs["nn"]
+            # Check that the voltage along stabilizer is equal to the 
+            # voltage along superconductor (i.e, check the reliability 
+            # of the current divider).
+            if all(np.isclose(v_stab,v_sc)) == False:
+                raise ValueError(f"Voltage difference along superconductor and stabilizer must be the same.")
+
+        return self.dict_Gauss_pt["electric_resistance"]
\ No newline at end of file
diff --git a/source_code/strand_component.py b/source_code/strand_component.py
index f1314a2..db4e685 100644
--- a/source_code/strand_component.py
+++ b/source_code/strand_component.py
@@ -270,7 +270,7 @@ def eval_tcs(self, dict_dummy):
         # evaluation of the current densities the total superconducting 
         # perpendicular cross section is always used.
         jop = (
-            np.abs(dict_dummy["op_current"])
+            np.abs(dict_dummy["current_along"])
             / (self.cross_section["sc"])
         )
 
diff --git a/source_code/strand_mixed_component.py b/source_code/strand_mixed_component.py
index 76aeb11..218d64c 100644
--- a/source_code/strand_mixed_component.py
+++ b/source_code/strand_mixed_component.py
@@ -641,6 +641,19 @@ def solve_current_divider(
                 maxiter=10,
                 disp=False,
             )
+        
+        # Tolerance on newton halley increased in case I_critical is very small,
+        # to avoid inaccuracies on the divider that could lead to potential 
+        # differences between sc and stab
+        if min(critical_current)>1e-6:
+            # Default tollerance in optimize.newton method
+            tollerance = 1.48e-8
+        else:
+            # Value found trial and error iteration
+            tollerance = 1e-12
+            # Other possible solution for the correct tollerance
+            # tollerance = min(critical_current)/1000
+
         # Evaluate superconducting with Halley's method
         sc_current = optimize.newton(
             self.__sc_current_residual,
@@ -648,6 +661,7 @@ def solve_current_divider(
             args=(psi, current),
             fprime=self.__d_sc_current_residual,
             fprime2=self.__d2_sc_current_residual,
+            tol = tollerance,
             maxiter=1000,
         )
 
@@ -797,149 +811,40 @@ def get_electric_resistance(self, conductor: object) -> np.ndarray:
         # still identify a normal region.
         if ind_not_zero.any():
 
-            # Get index that correspond to superconducting regime.
-            ind_sc_gauss = np.nonzero(
-                self.dict_Gauss_pt["op_current_sc"][ind_not_zero] / critical_current_gauss[ind_not_zero] < 0.95
-            )[0]
-            # Get index that correspond to the normal regime.
-            ind_normal_gauss = np.nonzero(
-                self.dict_Gauss_pt["op_current_sc"][ind_not_zero] / critical_current_gauss[ind_not_zero] >= 0.95
-            )[0]
-
-            ## SUPERCONDUCTING REGIME ##
-
-            # Check if np array ind_sc_gauss is not empty.
-            if ind_sc_gauss.any():
-                # Current in superconducting regime is the carried by
-                # superconducting material only.
-                self.dict_Gauss_pt["op_current"][ind_not_zero[ind_sc_gauss]] = self.dict_Gauss_pt[
-                    "op_current_sc"
-                ][ind_not_zero[ind_sc_gauss]]
-
-                # Compute superconducting electrical resistivity only in index 
-                # for which the superconducting regime is guaranteed, using the 
-                # power low.
-                self.dict_Gauss_pt["electrical_resistivity_superconductor"][ind_not_zero[ind_sc_gauss]] = self.superconductor_power_law(
-                    self.dict_Gauss_pt["op_current"][ind_not_zero[ind_sc_gauss]],
-                    critical_current_gauss[ind_not_zero[ind_sc_gauss]],
-                    self.dict_Gauss_pt["J_critical"][ind_not_zero[ind_sc_gauss]]
+            sc_current_gauss, stab_current_gauss = self.solve_current_divider(
+                self.dict_Gauss_pt["electrical_resistivity_stabilizer"][ind_not_zero],
+                critical_current_gauss[ind_not_zero],
+                self.dict_Gauss_pt["op_current"][ind_not_zero]
                 )
-
-                # Evaluate electic resistance in superconducting region 
-                # (superconductor only).
-                self.dict_Gauss_pt["electric_resistance"][ind_not_zero[
-                    ind_sc_gauss
-                ]] = self.electric_resistance(
-                    conductor, "electrical_resistivity_superconductor", "sc", ind_not_zero[ind_sc_gauss]
-                )
-
-            ## NORMAL REGIME ##
-
-            # Check if np array ind_normal_gauss is not empty.
-            if ind_normal_gauss.any():
-
-                # Evaluate how the current is distributed solving the current
-                # divider problem in Gauss point.
-                sc_current_gauss, stab_current_gauss = self.solve_current_divider(
-                    self.dict_Gauss_pt["electrical_resistivity_stabilizer"][ind_not_zero[ind_normal_gauss]],
-                    critical_current_gauss[ind_not_zero[ind_normal_gauss]],
-                    self.dict_Gauss_pt["op_current"][ind_not_zero[ind_normal_gauss]]
+            
+            self.dict_Gauss_pt["electrical_resistivity_superconductor"][
+                ind_not_zero
+            ] = self.superconductor_power_law(
+                sc_current_gauss,
+                critical_current_gauss[ind_not_zero],
+                self.dict_Gauss_pt["J_critical"][ind_not_zero],
                 )
+            # Evaluate the equivalent electric resistance in Ohm.
+            self.dict_Gauss_pt["electric_resistance"][
+                ind_not_zero] = self.parallel_electric_resistance(
+                conductor,
+                [
+                    "electrical_resistivity_superconductor",
+                    "electrical_resistivity_stabilizer",
+                ],["sc","stab"],
+                ind_not_zero,
+            )
 
-                # Get index of the normal region where all the current is
-                # carried by the stabilizer, to avoid division by 0 in 
-                # evaluation of superconducting electrical resistivity with the 
-                # power law.
-                ind_stab_gauss = np.nonzero(
-                    (
-                        stab_current_gauss / self.dict_Gauss_pt["op_current"][ind_not_zero[ind_normal_gauss]]
-                        > 0.999999
-                    )
-                    | (sc_current_gauss < 1.0)
-                )[0]
-
-                ## CURRENT CARRIED BY THE STABILIZER ##
-                # Check if np array ind_stab_gauss is not empty.
-                if ind_stab_gauss.any():
-                    # Get the index of location of current sharing region, if 
-                    # any.
-                    ind_sh_gauss = np.nonzero(
-                        (
-                            stab_current_gauss
-                            / self.dict_Gauss_pt["op_current"][ind_not_zero[ind_normal_gauss]]
-                            <= 0.999999
-                        )
-                        | (sc_current_gauss >= 1.0)
-                    )[0]
-                    
-                    # Check if nparray ind_sh_gauss is not empty.
-                    if ind_sh_gauss.any():
-                        # ind_sh_gauss is not empty.
-                        # Get final index of the location of the current 
-                        # sharing zone, keeping into account that 
-                        # sc_current_gauss is already filtered on 
-                        # ind_not_zero[ind_normal_gauss].
-                        ind_shf_gauss = {1:ind_sh_gauss,2:ind_not_zero[ind_normal_gauss[ind_sh_gauss]]}
-                    else:
-                        # ind_sh_gauss is empty.
-                        # Get final index of the location of the current 
-                        # sharing zone, keeping into account that 
-                        # sc_current_gauss is already filtered on 
-                        # ind_not_zero[ind_normal_gauss] and that ind_sh_gauss 
-                        # is empty.
-                        ind_shf_gauss = {1:ind_sh_gauss,2:ind_not_zero[ind_normal_gauss]}
-
-                    # Evaluate electic resistance in normal region (stabilizer 
-                    # only).
-                    self.dict_Gauss_pt["electric_resistance"][
-                        ind_shf_gauss[2]
-                    ] = self.electric_resistance(
-                        conductor, "electrical_resistivity_stabilizer", "stab", ind_shf_gauss[2]
-                    )
-                else:
-                    # Get final index of the location of the current sharing 
-                    # zone, keeping into account that sc_current_gauss is 
-                    # already filtered on ind_not_zero[ind_normal_gauss]. In 
-                    # this case the current is shared between the 
-                    # superconductor and the stabilizer, so we need to exploit 
-                    # all the index in array ind_normal_gauss for the array 
-                    # sc_current_gauss.
-                    ind_shf_gauss = {1:np.nonzero(ind_normal_gauss>=0)[0],2:ind_not_zero[ind_normal_gauss]}
-
-                ## CURRENT SHARED BY THE SUPERCONDUCTOR AND THE STABILIZER ##
-                if ind_shf_gauss[1].any():
-                    # Evaluate the electrical resistivity of the superconductor
-                    # according to the power low in Gauss point in Ohm*m.
-                    self.dict_Gauss_pt["electrical_resistivity_superconductor"][
-                        ind_shf_gauss[2]
-                    ] = self.superconductor_power_law(
-                        sc_current_gauss[ind_shf_gauss[1]],
-                        critical_current_gauss[ind_shf_gauss[2]],
-                        self.dict_Gauss_pt["J_critical"][ind_shf_gauss[2]],
-                    )
-
-                    # Evaluate the equivalent electric resistance in Ohm.
-                    self.dict_Gauss_pt["electric_resistance"][
-                        ind_shf_gauss[2]
-                    ] = self.parallel_electric_resistance(
-                        conductor,
-                        [
-                            "electrical_resistivity_superconductor",
-                            "electrical_resistivity_stabilizer",
-                        ],["sc","stab"],
-                        ind_shf_gauss[2],
-                    )
-                    
-                    # Compute voltage along stabilizer.
-                    v_stab = self.dict_Gauss_pt["electrical_resistivity_stabilizer"][
-                        ind_shf_gauss[2]
-                    ] * stab_current_gauss[ind_shf_gauss[1]] / self.cross_section["stab"]
-                    # Compute voltage along superconductor
-                    v_sc = self.inputs["E0"] * (sc_current_gauss[ind_shf_gauss[1]] / critical_current_gauss[ind_shf_gauss[2]]) ** self.inputs["nn"]
-                    # Check that the voltage along stabilizer is equal to the 
-                    # voltage along superconductor (i.e, check the reliability 
-                    # of the current divider).
-                    if all(np.isclose(v_stab,v_sc)) == False:
-                        raise ValueError(f"Voltage difference along superconductor and stabilizer must be the same.")
+            # Compute voltage along stabilizer.
+            v_stab = self.dict_Gauss_pt["electrical_resistivity_stabilizer"][
+                ind_not_zero
+            ] * stab_current_gauss / self.cross_section["stab"]
+            # Compute voltage along superconductor
+            v_sc = self.inputs["E0"] * (sc_current_gauss / critical_current_gauss[ind_not_zero]) ** self.inputs["nn"]
+            # Check that the voltage along stabilizer is equal to the 
+            # voltage along superconductor (i.e, check the reliability 
+            # of the current divider).
+            if all(np.isclose(v_stab,v_sc)) == False:
+                raise ValueError(f"Voltage difference along superconductor and stabilizer must be the same.")
 
         return self.dict_Gauss_pt["electric_resistance"]
\ No newline at end of file