Unsure about the cause of the DiracDelta error with my code using symfit

Question

I have been carrying out QMSA analysis of PdSn4 and when trying to implement a multicarrier fitting I keep getting an error

NameError: name 'DiracDelta' is not defined

My output includes only one instance where the DiracDelta function is present when looking at the second band fitting (n2, mu2). From what I understand there shouldn't be any DiracDelta functions within the output based on my model design.

import sympy as sp
from sympy.functions.special.delta_functions import DiracDelta
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt 
import math
import symfit
from symfit import DiracDelta
from symfit import Variable, Parameter, parameters, Fit, Model, cosh
from symfit.contrib.interactive_guess import InteractiveGuess
import numpy as np
plt.ion()

data = np.loadtxt('PdSn4_010K_rhoData.txt')
xdata = np.array(data[:, 0])
rhoxx_exp = np.array(data[:, 1])
rhoxy_exp = np.array(data[:, 2])

B = Variable('B')
y1 = Variable('y1')
y2 = Variable('y2')
n_1 = Parameter('n1', value=n1, min=1.5e27, max=1.55e27) 
n_2 = Parameter('n2', value=n2, min=-1.52e27, max=-1.51e27)
n_3 = Parameter('n3', value=n3, min=0.91e27, max=0.92e27)
n_4 = Parameter('n4', value=n4, min=-1.05e27, max=-1.045e27)
mu_1 = Parameter('mu1', value=mu1, min=1.4, max=1.42)
mu_2 = Parameter('mu2', value=mu2, min=0.63, max=0.65)
mu_3 = Parameter('mu3', value=mu3, min=0.085, max=0.09)
mu_4 = Parameter('mu4', value=mu4, min=0.088, max=0.092)

a = 1.1513
f = 0.992
q = 1.6e-19

R1 = -1 / (q*n_1)
R2 = -1 / (q*n_2)
R3 = -1 / (q*n_3)
R4 = -1 / (q*n_4)

rho1 = 1 / (abs(n_1) * q * mu_1)
rho2 = 1 / (abs(n_2) * q * mu_2)
rho3 = 1 / (abs(n_3) * q * mu_3)
rho4 = 1 / (abs(n_4) * q * mu_4)
    
dictionary = {

    y1: (rho1 / (rho1**2 + (R1*B)**2) + rho2 / (rho2**2 + (R2*B)**2) 
    
    + rho3 / (rho3**2 + (R3*B)**2) + rho4 / (rho4**2 + (R4*B)**2)) 
    
    / ((rho1 / (rho1**2 + (R1*B)**2) + rho2 / (rho2**2 + (R2*B)**2) 
    
    + rho3 / (rho3**2 + (R3*B)**2) + rho4 / (rho4**2 + (R4*B)**2))**2
       
       + ((-R1*B)/(rho1**2+(R1*B)**2) + (-R2*B)/(rho2**2+(R2*B)**2)
                                                
       + (-R3*B)/(rho3**2+(R3*B)**2) + (-R4*B)/(rho4**2+(R4*B)**2))**2)
    
    + (a/cosh(B/f))*((rho1 / (rho1**2 + (R1*B)**2) + rho2 / (rho2**2 + (R2*B)**2) 
    
    + rho3 / (rho3**2 + (R3*B)**2) + rho4 / (rho4**2 + (R4*B)**2)) 
    
    / ((rho1 / (rho1**2 + (R1*B)**2) + rho2 / (rho2**2 + (R2*B)**2) 
    
    + rho3 / (rho3**2 + (R3*B)**2) + rho4 / (rho4**2 + (R4*B)**2))**2
       
       + ((-R1*B)/(rho1**2+(R1*B)**2) + (-R2*B)/(rho2**2+(R2*B)**2)
                                                
       + (-R3*B)/(rho3**2+(R3*B)**2) + (-R4*B)/(rho4**2+(R4*B)**2))**2))
    ,
    
   
       
    y2: - ((-R1*B)/(rho1**2+(R1*B)**2) + (-R2*B)/(rho2**2+(R2*B)**2)
                                                
       + (-R3*B)/(rho3**2+(R3*B)**2) + (-R4*B)/(rho4**2+(R4*B)**2)) / 
    
  ((rho1 / (rho1**2 + (R1*B)**2) + rho2 / (rho2**2 + (R2*B)**2) 
    
    + rho3 / (rho3**2 + (R3*B)**2) + rho4 / (rho4**2 + (R4*B)**2))**2
       
       + ((-R1*B)/(rho1**2+(R1*B)**2) + (-R2*B)/(rho2**2+(R2*B)**2)
                                                
       + (-R3*B)/(rho3**2+(R3*B)**2) + (-R4*B)/(rho4**2+(R4*B)**2))**2)  
    

  
}

model = Model(dictionary)
#print(model)
model_sim = model(B=xdata, n1=n1, n2=n2, n3=n3, n4=n4, mu1=mu1, mu2=mu2, mu3=mu3, mu4=mu4)
#print(data)
rho_xx_sim = model_sim.y1
rho_xy_sim = model_sim.y2
f = plt.figure(figsize=(15,15))
ax = f.add_subplot(211)
plt.xlabel("Magnetic Field [T]")
plt.ylabel("Rho$_{xx}$")
plt.plot(xdata, rhoxx_exp, 'k.')
plt.plot(xdata, rho_xx_sim, 'r-')
plt.xlim([0,8])


ax2 = f.add_subplot(212)
plt.xlabel("Magnetic Field [T]")
plt.ylabel("Rho$_{xy}$")
plt.plot(xdata, rhoxy_exp, 'k.')
plt.plot(xdata, rho_xy_sim, 'r-')
plt.xlim([0,8])

plt.show()

guess = InteractiveGuess(dictionary, B=xdata, y1=rho_xx_sim, y2=rho_xy_sim, n_points=250)
guess.execute()
print(guess)

fit = Fit(model, B=xdata, y1=rho_xx_sim, y2=rho_xy_sim)
fit_result = fit.execute()
print(fit_result)

Link to the data: Multicarrier Fitting Data It should be saved as PdSn4_010K_rhoData.txt in the same directory as the file.

If anyone could figure out why there is this DiracDelta function present in my output or how to overcome this error that would be great!

Error traceback goes as follows:

NameError                                 Traceback (most recent call last)
<ipython-input-19-68c235066887> in <module>
    197 
    198 fit = Fit(model, B=xdata, y1=rho_xx_sim, y2=rho_xy_sim)
--> 199 fit_result = fit.execute()
    200 print(fit_result)
    201 #y1result, y2result = fit.model(B=xdata, **fit_result.params)

~\anaconda5\lib\site-packages\symfit\core\fit.py in execute(self, **minimize_options)
    579         """
    580         minimizer_ans = self.minimizer.execute(**minimize_options)
--> 581         minimizer_ans.covariance_matrix = self.covariance_matrix(
    582             dict(zip(self.model.params, minimizer_ans._popt))
    583         )

~\anaconda5\lib\site-packages\symfit\core\fit.py in covariance_matrix(self, best_fit_params)
    278         :return: covariance matrix.
    279         """
--> 280         cov_matrix = self._covariance_matrix(best_fit_params,
    281                                              objective=self.objective)
    282         if cov_matrix is None:

~\anaconda5\lib\site-packages\symfit\core\fit.py in _covariance_matrix(self, best_fit_params, objective)
    236         # Helper function for self.covariance_matrix.
    237         try:
--> 238             hess = objective.eval_hessian(**key2str(best_fit_params))
    239         except AttributeError:
    240             # Some models do not have an eval_hessian, in which case we give up

~\anaconda5\lib\site-packages\symfit\core\objectives.py in eval_hessian(self, ordered_parameters, **parameters)
    367             ordered_parameters, **parameters
    368         )
--> 369         evaluated_hess = super(LeastSquares, self).eval_hessian(
    370             ordered_parameters, **parameters
    371         )

~\anaconda5\lib\site-packages\symfit\core\objectives.py in eval_hessian(self, ordered_parameters, **parameters)
    221         parameters.update(dict(zip(self.model.free_params, ordered_parameters)))
    222         parameters.update(self._invariant_kwargs)
--> 223         result = self.model.eval_hessian(**key2str(parameters))._asdict()
    224         # Return only the components corresponding to the dependent data.
    225         return self._shape_of_dependent_data(

~\anaconda5\lib\site-packages\symfit\core\models.py in eval_hessian(self, *args, **kwargs)
    897         # Evaluate the hessian model and use the resulting Ans namedtuple as a
    898         # dict. From this, take the relevant components.
--> 899         eval_hess_dict = self.hessian_model(*args, **kwargs)._asdict()
    900         hess = [[[np.broadcast_to(eval_hess_dict.get(D(var, p1, p2), 0),
    901                                   eval_hess_dict[var].shape)

~\anaconda5\lib\site-packages\symfit\core\models.py in __call__(self, *args, **kwargs)
    665             even for scalar valued functions. This is done for consistency.
    666         """
--> 667         return ModelOutput(self.keys(), self.eval_components(*args, **kwargs))
    668 
    669 

~\anaconda5\lib\site-packages\symfit\core\models.py in eval_components(self, *args, **kwargs)
    613                 dependencies_kwargs = {d.name: kwargs[d.name]
    614                                        for d in dependencies}
--> 615                 kwargs[symbol.name] = components[symbol](**dependencies_kwargs)
    616 
    617         return [np.atleast_1d(kwargs[var.name]) for var in self]

<lambdifygenerated-1024> in _lambdifygenerated(B, mu1, mu2, mu3, mu4, n1, n2, n3, n4)
      1 def _lambdifygenerated(B, mu1, mu2, mu3, mu4, n1, n2, n3, n4):
----> 2     return (6.5536e-76*(-3.90625e+37*(1.31072e-75*(7.8125e+37*B**2/n2**3 + 7.8125e+37/(mu2**2*n2**3))/(mu2*(B**2/n2**2 + 1/(mu2**2*n2**2))**2*abs(n2)) - 2*sign(n2)/(mu2*n2**2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))))*(1/(mu4*(3.90625e+37*B**2/n4**2 + 3.90625e+37/(mu4**2*n4**2))*abs(n4)) + 1/(mu3*(3.90625e+37*B**2/n3**2 + 3.90625e+37/(mu3**2*n3**2))*abs(n3)) + 1/(mu2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))*abs(n2)) + 1/(mu1*(3.90625e+37*B**2/n1**2 + 3.90625e+37/(mu1**2*n1**2))*abs(n1))) - 3.90625e+37*(1.31072e-75*B*(7.8125e+37*B**2/n2**3 + 7.8125e+37/(mu2**2*n2**3))/(n2*(B**2/n2**2 + 1/(mu2**2*n2**2))**2) - 2*B/(n2**2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))))*(B/(n4*(3.90625e+37*B**2/n4**2 + 3.90625e+37/(mu4**2*n4**2))) + B/(n3*(3.90625e+37*B**2/n3**2 + 3.90625e+37/(mu3**2*n3**2))) + B/(n2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))) + B/(n1*(3.90625e+37*B**2/n1**2 + 3.90625e+37/(mu1**2*n1**2)))))*(-2*(1.31072e-75*(7.8125e+37*B**2/n2**3 + 7.8125e+37/(mu2**2*n2**3))/(mu2*(B**2/n2**2 + 1/(mu2**2*n2**2))**2*abs(n2)) - 2*sign(n2)/(mu2*n2**2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))))*(1/(mu4*(3.90625e+37*B**2/n4**2 + 3.90625e+37/(mu4**2*n4**2))*abs(n4)) + 1/(mu3*(3.90625e+37*B**2/n3**2 + 3.90625e+37/(mu3**2*n3**2))*abs(n3)) + 1/(mu2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))*abs(n2)) + 1/(mu1*(3.90625e+37*B**2/n1**2 + 3.90625e+37/(mu1**2*n1**2))*abs(n1))) - 2*(1.31072e-75*B*(7.8125e+37*B**2/n2**3 + 7.8125e+37/(mu2**2*n2**3))/(n2*(B**2/n2**2 + 1/(mu2**2*n2**2))**2) - 2*B/(n2**2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))))*(B/(n4*(3.90625e+37*B**2/n4**2 + 3.90625e+37/(mu4**2*n4**2))) + B/(n3*(3.90625e+37*B**2/n3**2 + 3.90625e+37/(mu3**2*n3**2))) + B/(n2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))) + B/(n1*(3.90625e+37*B**2/n1**2 + 3.90625e+37/(mu1**2*n1**2)))))*(-6.25e+18*B/(n4*(3.90625e+37*B**2/n4**2 + 3.90625e+37/(mu4**2*n4**2))) - 6.25e+18*B/(n3*(3.90625e+37*B**2/n3**2 + 3.90625e+37/(mu3**2*n3**2))) - 6.25e+18*B/(n2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))) - 6.25e+18*B/(n1*(3.90625e+37*B**2/n1**2 + 3.90625e+37/(mu1**2*n1**2))))/((1/(mu4*(3.90625e+37*B**2/n4**2 + 3.90625e+37/(mu4**2*n4**2))*abs(n4)) + 1/(mu3*(3.90625e+37*B**2/n3**2 + 3.90625e+37/(mu3**2*n3**2))*abs(n3)) + 1/(mu2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))*abs(n2)) + 1/(mu1*(3.90625e+37*B**2/n1**2 + 3.90625e+37/(mu1**2*n1**2))*abs(n1)))**2 + (B/(n4*(3.90625e+37*B**2/n4**2 + 3.90625e+37/(mu4**2*n4**2))) + B/(n3*(3.90625e+37*B**2/n3**2 + 3.90625e+37/(mu3**2*n3**2))) + B/(n2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))) + B/(n1*(3.90625e+37*B**2/n1**2 + 3.90625e+37/(mu1**2*n1**2))))**2)**3 + 1.31072e-75*(-3.90625e+37*(1.31072e-75*(7.8125e+37*B**2/n2**3 + 7.8125e+37/(mu2**2*n2**3))/(mu2*(B**2/n2**2 + 1/(mu2**2*n2**2))**2*abs(n2)) - 2*sign(n2)/(mu2*n2**2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))))*(1/(mu4*(3.90625e+37*B**2/n4**2 + 3.90625e+37/(mu4**2*n4**2))*abs(n4)) + 1/(mu3*(3.90625e+37*B**2/n3**2 + 3.90625e+37/(mu3**2*n3**2))*abs(n3)) + 1/(mu2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))*abs(n2)) + 1/(mu1*(3.90625e+37*B**2/n1**2 + 3.90625e+37/(mu1**2*n1**2))*abs(n1))) - 3.90625e+37*(1.31072e-75*B*(7.8125e+37*B**2/n2**3 + 7.8125e+37/(mu2**2*n2**3))/(n2*(B**2/n2**2 + 1/(mu2**2*n2**2))**2) - 2*B/(n2**2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))))*(B/(n4*(3.90625e+37*B**2/n4**2 + 3.90625e+37/(mu4**2*n4**2))) + B/(n3*(3.90625e+37*B**2/n3**2 + 3.90625e+37/(mu3**2*n3**2))) + B/(n2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))) + B/(n1*(3.90625e+37*B**2/n1**2 + 3.90625e+37/(mu1**2*n1**2)))))*(-4.096e-57*B*(7.8125e+37*B**2/n2**3 + 7.8125e+37/(mu2**2*n2**3))/(n2*(B**2/n2**2 + 1/(mu2**2*n2**2))**2) + 6.25e+18*B/(n2**2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))))/((1/(mu4*(3.90625e+37*B**2/n4**2 + 3.90625e+37/(mu4**2*n4**2))*abs(n4)) + 1/(mu3*(3.90625e+37*B**2/n3**2 + 3.90625e+37/(mu3**2*n3**2))*abs(n3)) + 1/(mu2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))*abs(n2)) + 1/(mu1*(3.90625e+37*B**2/n1**2 + 3.90625e+37/(mu1**2*n1**2))*abs(n1)))**2 + (B/(n4*(3.90625e+37*B**2/n4**2 + 3.90625e+37/(mu4**2*n4**2))) + B/(n3*(3.90625e+37*B**2/n3**2 + 3.90625e+37/(mu3**2*n3**2))) + B/(n2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))) + B/(n1*(3.90625e+37*B**2/n1**2 + 3.90625e+37/(mu1**2*n1**2))))**2)**2 + (-4.096e-57*B*(-2.34375e+38*B**2/n2**4 - 2.34375e+38/(mu2**2*n2**4))/(n2*(B**2/n2**2 + 1/(mu2**2*n2**2))**2) - 4.096e-57*B*(4*B**2/n2**3 + 4/(mu2**2*n2**3))*(7.8125e+37*B**2/n2**3 + 7.8125e+37/(mu2**2*n2**3))/(n2*(B**2/n2**2 + 1/(mu2**2*n2**2))**3) + 8.192e-57*B*(7.8125e+37*B**2/n2**3 + 7.8125e+37/(mu2**2*n2**3))/(n2**2*(B**2/n2**2 + 1/(mu2**2*n2**2))**2) - 1.25e+19*B/(n2**3*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))))/(3.90625e+37*(1/(mu4*(3.90625e+37*B**2/n4**2 + 3.90625e+37/(mu4**2*n4**2))*abs(n4)) + 1/(mu3*(3.90625e+37*B**2/n3**2 + 3.90625e+37/(mu3**2*n3**2))*abs(n3)) + 1/(mu2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))*abs(n2)) + 1/(mu1*(3.90625e+37*B**2/n1**2 + 3.90625e+37/(mu1**2*n1**2))*abs(n1)))**2 + 3.90625e+37*(B/(n4*(3.90625e+37*B**2/n4**2 + 3.90625e+37/(mu4**2*n4**2))) + B/(n3*(3.90625e+37*B**2/n3**2 + 3.90625e+37/(mu3**2*n3**2))) + B/(n2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))) + B/(n1*(3.90625e+37*B**2/n1**2 + 3.90625e+37/(mu1**2*n1**2))))**2) + 6.5536e-76*((-5.12e-38*(7.8125e+37*B**2/n2**3 + 7.8125e+37/(mu2**2*n2**3))/(mu2*(B**2/n2**2 + 1/(mu2**2*n2**2))**2*abs(n2)) + 7.8125e+37*sign(n2)/(mu2*n2**2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))))*(6.5536e-76*(7.8125e+37*B**2/n2**3 + 7.8125e+37/(mu2**2*n2**3))/(mu2*(B**2/n2**2 + 1/(mu2**2*n2**2))**2*abs(n2)) - sign(n2)/(mu2*n2**2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2)))) + (-5.12e-38*B*(7.8125e+37*B**2/n2**3 + 7.8125e+37/(mu2**2*n2**3))/(n2*(B**2/n2**2 + 1/(mu2**2*n2**2))**2) + 7.8125e+37*B/(n2**2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))))*(6.5536e-76*B*(7.8125e+37*B**2/n2**3 + 7.8125e+37/(mu2**2*n2**3))/(n2*(B**2/n2**2 + 1/(mu2**2*n2**2))**2) - B/(n2**2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2)))) + (1/(mu4*(3.90625e+37*B**2/n4**2 + 3.90625e+37/(mu4**2*n4**2))*abs(n4)) + 1/(mu3*(3.90625e+37*B**2/n3**2 + 3.90625e+37/(mu3**2*n3**2))*abs(n3)) + 1/(mu2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))*abs(n2)) + 1/(mu1*(3.90625e+37*B**2/n1**2 + 3.90625e+37/(mu1**2*n1**2))*abs(n1)))*(-5.12e-38*(-2.34375e+38*B**2/n2**4 - 2.34375e+38/(mu2**2*n2**4))/(mu2*(B**2/n2**2 + 1/(mu2**2*n2**2))**2*abs(n2)) - 5.12e-38*(4*B**2/n2**3 + 4/(mu2**2*n2**3))*(7.8125e+37*B**2/n2**3 + 7.8125e+37/(mu2**2*n2**3))/(mu2*(B**2/n2**2 + 1/(mu2**2*n2**2))**3*abs(n2)) + 1.024e-37*(7.8125e+37*B**2/n2**3 + 7.8125e+37/(mu2**2*n2**3))*sign(n2)/(mu2*n2**2*(B**2/n2**2 + 1/(mu2**2*n2**2))**2) + 1.5625e+38*DiracDelta(n2)/(mu2*n2**2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))) - 1.5625e+38*sign(n2)/(mu2*n2**3*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2)))) + (B/(n4*(3.90625e+37*B**2/n4**2 + 3.90625e+37/(mu4**2*n4**2))) + B/(n3*(3.90625e+37*B**2/n3**2 + 3.90625e+37/(mu3**2*n3**2))) + B/(n2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))) + B/(n1*(3.90625e+37*B**2/n1**2 + 3.90625e+37/(mu1**2*n1**2))))*(-5.12e-38*B*(-2.34375e+38*B**2/n2**4 - 2.34375e+38/(mu2**2*n2**4))/(n2*(B**2/n2**2 + 1/(mu2**2*n2**2))**2) - 5.12e-38*B*(4*B**2/n2**3 + 4/(mu2**2*n2**3))*(7.8125e+37*B**2/n2**3 + 7.8125e+37/(mu2**2*n2**3))/(n2*(B**2/n2**2 + 1/(mu2**2*n2**2))**3) + 1.024e-37*B*(7.8125e+37*B**2/n2**3 + 7.8125e+37/(mu2**2*n2**3))/(n2**2*(B**2/n2**2 + 1/(mu2**2*n2**2))**2) - 1.5625e+38*B/(n2**3*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2)))))*(-6.25e+18*B/(n4*(3.90625e+37*B**2/n4**2 + 3.90625e+37/(mu4**2*n4**2))) - 6.25e+18*B/(n3*(3.90625e+37*B**2/n3**2 + 3.90625e+37/(mu3**2*n3**2))) - 6.25e+18*B/(n2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))) - 6.25e+18*B/(n1*(3.90625e+37*B**2/n1**2 + 3.90625e+37/(mu1**2*n1**2))))/((1/(mu4*(3.90625e+37*B**2/n4**2 + 3.90625e+37/(mu4**2*n4**2))*abs(n4)) + 1/(mu3*(3.90625e+37*B**2/n3**2 + 3.90625e+37/(mu3**2*n3**2))*abs(n3)) + 1/(mu2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))*abs(n2)) + 1/(mu1*(3.90625e+37*B**2/n1**2 + 3.90625e+37/(mu1**2*n1**2))*abs(n1)))**2 + (B/(n4*(3.90625e+37*B**2/n4**2 + 3.90625e+37/(mu4**2*n4**2))) + B/(n3*(3.90625e+37*B**2/n3**2 + 3.90625e+37/(mu3**2*n3**2))) + B/(n2*(3.90625e+37*B**2/n2**2 + 3.90625e+37/(mu2**2*n2**2))) + B/(n1*(3.90625e+37*B**2/n1**2 + 3.90625e+37/(mu1**2*n1**2))))**2)**2)

NameError: name 'DiracDelta' is not defined