def _findRootAttempt(self, startingEne, multipler, type):
     try:
         fun = lambda e: self._getDet(e, multipler)
         if type == 'muller':
             return complex(mpmath.findroot(fun, startingEne, solver='muller', maxsteps=10000, tol=2.16840434497100886801e-19))
         else:
             return complex(mpmath.findroot(fun, startingEne[0], solver='secant'))
     except ValueError:
         return None

def getZero(Func, init=(1.+1.j, 1.1+1.j, 1.+1.1j), solver='muller', ftol=1.e-3, tol=1.e-3, maxsteps=600):
  import mpmath as mp
  w0 = complex(float('nan'), float('nan'))
  try:
     w0 = mp.findroot(lambda w : Func(complex(w)), init, solver=solver, maxsteps=maxsteps, ftol=ftol, tol=tol)
  except:
   try:
     w0 = mp.findroot(lambda w : Func(complex(w)), init[0], solver='halley', maxsteps=maxsteps, ftol=ftol, tol=tol)
   except:   
     print "Could not find zero"

  return complex(w0)

 def draw_waiting_time(self, avail_devices, failed_devices):
     self.curr_avail_devices = avail_devices
     self.curr_failed_devices = failed_devices
     self.curr_uniform_variate = random.uniform(0,1)
     
     while self.curr_uniform_variate == 0:
             self.curr_uniform_variate = random.uniform(0,1)
     
     try:
         wait_time = findroot(self.func, [0,100], solver='secant')
     except:
         wait_time = findroot(self.func, [0,100], solver='secant', maxsteps=1000)
         #wait_time = findroot(self.func, [0, 87600] , maxsteps=1000, solver='secant', verify=False)
         
     return abs(wait_time)

def nodes(n):
    left  = mp.mpf(0)
    right = mp.mpf(n+(n-1)*sqrt(n)*1.01)

    i = 2
    factor = 2

    while True:
        l = [ (x,mp.laguerre(n,0,x)) for x in mp.linspace(left,right,n*factor**i)]

        intervals = []
        for j in range(len(l)-1):
            prod = l[j][1]*l[j+1][1]
            if prod < 0:
                intervals.append([l[j][0], l[j+1][0]])

        if len(intervals) == n:
            break
        i += 1

    roots = []
    f = lambda x: mp.laguerre(n,0,x)
    for ab in intervals:
        a,b = ab
        try:
            z = mp.findroot(f, (a, b), tol=1e-50, solver='bisect')
        except:
            z = bisect(f, a, b, tol=1e-50)
        roots.append( z )

    return roots

    def update_mix(self):
        if self.esolute.shape != self.esolvent.shape:
            return
        
        eeff=np.empty(self.esolute.shape, dtype=complex)
        for i in range(len(eeff)):
            e1 = self.esolute[i]#.real    
            em = self.esolvent[i]#.real   #THIS IS GENERALLY THE SOLVENT
            
            partialfunc = partial(self.fill_func_dielectric, em=em, e1=e1, v=self.Vfrac, w=self.w)
            
            
            eeff[i] = findroot(partialfunc, 
                               (0.5, 1, 2),  #This is a root that we are guessing for 3 params
                               solver='muller',
                               )            
            
            #eeff[i] = scipy.optimize.newton(self.fill_func_dielectric, 
                                            #self.Root,
                                            #fprime=None, 
                                            #args=(em, e1, self.Vfrac, self.w),
                                            #tol=.0001, 
                                            #maxiter=1000)
#			self.Root=eeff[i]  #Reset the root
        self.mixedarray=eeff

 def long_interval(w, n, t):
     """
     Return the interval of integration for longitudinal impedance, which
     is of the form [0, root, inf], where root satisfies
     d_l(root, w, n, t) = 0.
     
     """
     if w < mp.sqrt(1+n):
         return [0, mp.inf]
     if n == 0:
         guesses = [4, 6, 8, 10]
     else:
         guesses = [4, 6, 8, 10, 12, 14]
     int_range = [0, mp.inf]    
     for guess in guesses:
         try:
             root = mp.findroot(lambda z: BiMax.d_l(z, w, n, t), guess)
         except ValueError:
             continue
         unique = True
         for z in int_range:
             if mp.fabs(z - mp.fabs(root)) < 1e-4:
                 unique = False
         if unique:
             int_range += [mp.fabs(root)]
     int_range = np.sort(int_range)
     return int_range

    def generate_zeta(self):
        import mpmath
        import sys
        root_list		= {}
        for s in [complex(0.5,t) for t in range(10,100)]:
            try:
                root_complex	= mpmath.findroot(mpmath.zeta, s)
                if self.debug:
                    sys.stderr.write("[NOTICE] Found nontrivial zeta zero at %s.\n" % str(root_complex))
                root_imag	= float(root_complex.imag)
                try:
                    root_list[root_imag]
                    if self.debug:
                        sys.stderr.write("[NOTICE] Already found this zero, skipping.\n")
                    next
                except KeyError:
                    if root_imag > 1:
                        yield root_imag
                        root_list[root_imag] = True
                    else:
                        if self.debug:
                            sys.stderr.write("[NOTICE] Irrelevant zero, skipping.\n")
            except ValueError, e:
#DEBUG#                sys.stderr.write("Exception at s=%s: %s.\n" % (str(s), e.args[0])) #DEBUG#
                pass
            except OverflowError, e:
#DEBUG#                sys.stderr.write("Exception at s=%s: %s.\n" % (str(s), e.args[0])) #DEBUG#
                pass

 def int_interval(wrel, n, t, k):
     """
     find out the almost-singular point and
     divide the integration integrals.
     
     """
     if wrel < 1 or wrel > 1.3:
         return [0, mp.inf]
     else:
         guesses = [4, 6, 8, 10, 12, 14]
     wc = wrel * mp.sqrt(1+n)
     int_range = [0, mp.inf]
     for guess in guesses:
         try:
             root = mp.findroot(lambda zc: MaxKappa.e_l(zc, wc, n, t, k), guess)
         except ValueError:
             continue
         unique = True
         for z in int_range:
             if mp.fabs(z - mp.fabs(root)) < 1e-4:
                 unique = False
         if unique:
             int_range += [mp.fabs(root)]
     int_range = np.sort(int_range)
     return int_range

	def mini(r):
		nonlocal txi, tv
		xi = [r]
		y = f(r)
		v = r*f(r) + mpmath.quad(f, [r, mpmath.inf])
		if verbose:
			print('Trying r={0} (v={1})'.format(r, v), file=sys.stderr)
		for i in itertools.count():
			xm1 = xi[i]
			h = v / xm1
			y += h
			if y >= maximum or mpmath.almosteq(y, maximum, abs_eps=mpmath.mp.eps * 2**10):
				break
			# We solve for x via secant method instead of using f's inverse.
			x = mpmath.findroot(lambda x: f(x) - y, xm1)
			xi.append(x)
		xi.append(mpmath.mpf())
		if len(xi) == nseg:
			if mpmath.almosteq(y, maximum, abs_eps=mpmath.mp.eps * 2**10):
				txi, tv = xi[::-1], v
				return 0
			# If y > maximum, then v is too large, which means r is too far
			# left, so we want to return a negative value. The opposite holds
			# true when y < maximum.
			return maximum - y
		return len(xi) - nseg + h*mpmath.sign(len(xi) - nseg)

def solve(f, x0=None, nseg=128, verbose=False):
	"""Find r, v, and the x coordinates for f."""
	r = x0
	if r is None:
		if verbose:
			print('Calculating initial guess... ', end='', file=sys.stderr)
		# The area we seek is nseg equal-area rectangles surrounding f(x), not
		# f(x) itself, but we can get a good approximation from it.
		v = mpmath.quad(f, [0, mpmath.inf]) / nseg
		r = mpmath.findroot(lambda x: x*f(x) + mpmath.quad(f, [x, mpmath.inf]) - v, x0=1, x1=100, maxsteps=100)
		if verbose:
			print(r, file=sys.stderr)
	# We know that f(0) is the maximum because f must decrease monotonically.
	maximum = f(0)
	txi = []
	tv = mpmath.mpf()
	def mini(r):
		nonlocal txi, tv
		xi = [r]
		y = f(r)
		v = r*f(r) + mpmath.quad(f, [r, mpmath.inf])
		if verbose:
			print('Trying r={0} (v={1})'.format(r, v), file=sys.stderr)
		for i in itertools.count():
			xm1 = xi[i]
			h = v / xm1
			y += h
			if y >= maximum or mpmath.almosteq(y, maximum, abs_eps=mpmath.mp.eps * 2**10):
				break
			# We solve for x via secant method instead of using f's inverse.
			x = mpmath.findroot(lambda x: f(x) - y, xm1)
			xi.append(x)
		xi.append(mpmath.mpf())
		if len(xi) == nseg:
			if mpmath.almosteq(y, maximum, abs_eps=mpmath.mp.eps * 2**10):
				txi, tv = xi[::-1], v
				return 0
			# If y > maximum, then v is too large, which means r is too far
			# left, so we want to return a negative value. The opposite holds
			# true when y < maximum.
			return maximum - y
		return len(xi) - nseg + h*mpmath.sign(len(xi) - nseg)
	r = mpmath.findroot(mini, r)
	assert len(txi) == nseg
	if verbose:
		print('Done calculating r, v, x[i].', file=sys.stderr)
	return r, tv, txi

    def _eval_evalf(self, prec):
        """Evaluate this complex root to the given precision."""
        with workprec(prec):
            g = self.poly.gen
            if not g.is_Symbol:
                d = Dummy('x')
                func = lambdify(d, self.expr.subs({g: d}), "mpmath")
            else:
                func = lambdify(g, self.expr, "mpmath")

            try:
                interval = self.interval
            except DomainError:
                return super()._eval_evalf(prec)

            while True:
                if self.is_extended_real:
                    a = mpf(str(interval.a))
                    b = mpf(str(interval.b))
                    if a == b:
                        root = a
                        break
                    x0 = mpf(str(interval.center))
                else:
                    ax = mpf(str(interval.ax))
                    bx = mpf(str(interval.bx))
                    ay = mpf(str(interval.ay))
                    by = mpf(str(interval.by))
                    x0 = mpc(*map(str, interval.center))
                    if ax == bx and ay == by:
                        root = x0
                        break

                try:
                    root = findroot(func, x0)
                    # If the (real or complex) root is not in the 'interval',
                    # then keep refining the interval. This happens if findroot
                    # accidentally finds a different root outside of this
                    # interval because our initial estimate 'x0' was not close
                    # enough. It is also possible that the secant method will
                    # get trapped by a max/min in the interval; the root
                    # verification by findroot will raise a ValueError in this
                    # case and the interval will then be tightened -- and
                    # eventually the root will be found.
                    if self.is_extended_real:
                        if (a <= root <= b):
                            break
                    elif (ax <= root.real <= bx and ay <= root.imag <= by
                          and (interval.ay > 0 or interval.by < 0)):
                        break
                except (ValueError, UnboundLocalError):
                    pass
                self.refine()
                interval = self.interval

        return ((Float._new(root.real._mpf_, prec) if not self.is_imaginary else 0) +
                I*Float._new(root.imag._mpf_, prec))

    def _eval_evalf(self, prec):
        """Evaluate this complex root to the given precision. """
        with workprec(prec):
            g = self.poly.gen
            if not g.is_Symbol:
                d = Dummy('x')
                func = lambdify(d, self.expr.subs(g, d))
            else:
                func = lambdify(g, self.expr)

            interval = self._get_interval()
            if not self.is_real:
                # For complex intervals, we need to keep refining until the
                # imaginary interval is disjunct with other roots, that is,
                # until both ends get refined.
                ay = interval.ay
                by = interval.by
                while interval.ay == ay or interval.by == by:
                    interval = interval.refine()

            while True:
                if self.is_real:
                    x0 = mpf(str(interval.center))
                else:
                    x0 = mpc(*map(str, interval.center))
                try:
                    root = findroot(func, x0, verify=False)
                    # If the (real or complex) root is not in the 'interval',
                    # then keep refining the interval. This happens if findroot
                    # accidentally finds a different root outside of this
                    # interval because our initial estimate 'x0' was not close
                    # enough.
                    if self.is_real:
                        a = mpf(str(interval.a))
                        b = mpf(str(interval.b))
                        if a == b:
                            root = a
                            break
                        if not (a < root < b):
                            raise ValueError("Root not in the interval.")
                    else:
                        ax = mpf(str(interval.ax))
                        bx = mpf(str(interval.bx))
                        ay = mpf(str(interval.ay))
                        by = mpf(str(interval.by))
                        if ax == bx and ay == by:
                            root = ax + S.ImaginaryUnit*by
                            break
                        if not (ax < root.real < bx and ay < root.imag < by):
                            raise ValueError("Root not in the interval.")
                except ValueError:
                    interval = interval.refine()
                    continue
                else:
                    break

        return Float._new(root.real._mpf_, prec) + I*Float._new(root.imag._mpf_, prec)

 def solver(args):
     equation, lower_boundary, upper_boundary, initial_guess, f, use_mpmath = args
     if not use_mpmath:
         warnings.filterwarnings('ignore')
         solution = root(equation, initial_guess, args=f, method='hybr')
         solution = solution.x[0]
         warnings.resetwarnings()
     else:
         equation_wrap = partial(equation, f=f)
         try:
             solution = mp.findroot(equation_wrap, (lower_boundary, upper_boundary),
                                    maxsteps=1000, solver='anderson', tol=5e-16)
         except ValueError as err:
             print err
             print 'Lowering tolerance to 5e-6'
             solution = mp.findroot(equation_wrap, (lower_boundary, upper_boundary),
                                    maxsteps=1000, solver='anderson', tol=5e-6)
         solution = np.float(solution)
     return solution

def besselJZeros(m, a, b):
  require_mpmath()
  if not hasattr(mpmath, 'besseljzero'):
    besseljn = lambda x: mpmath.besselj(m, x)
    results = [mpmath.findroot(besseljn, mpmath.pi*(kp - 1./4 + 0.5*m)) for kp in range(a, b+1)]
  else:
    results = [mpmath.besseljzero(m, i) for i in xrange(a, b+1)]
  # Check that we haven't found double roots or missed a root. All roots should be separated by approximately pi
  assert all([0.6*mpmath.pi < (b - a) < 1.4*mpmath.pi for a, b in zip(results[:-1], results[1:])]), "Separation of Bessel zeros was incorrect."
  return results

def printAnaRoots(V, fndStates):
  print "\n\nFor V=" + str(V)
  def denumWrap(k):
    return denum(k,a,V)
  for state in fndStates:
    try:
      root = mpm.findroot(denumWrap, state)
      diff = abs(root-state)
      print str(root.real) + "\t" + str(root.imag) + "\t" + str(diff)
    except ValueError:
      print "Not Fnd"

def getIslandWidth(x,y,Psi):
    from scipy import interpolate
    import mpmath as mp      
    Nx = len(x)
    Ny = len(y)
    XPointValue_1 = Psi[0,Nx/2]
    XPointValue_2 = Psi[-1,Nx/2]

    
    # Interpolate Psi
    Psi_inter = interpolate.splrep(x,Psi[Ny/2,:],s=0)

    # Find zero
    def W(x) : 
      return XPointValue_1 - interpolate.splev(float(x),Psi_inter,der=0)
    try:
        w1 = float(mp.findroot(W,  1))
        w2 = float(mp.findroot(W, -1))
    except:
        w1 = 0.
        w2 = 0.
    return (w1, w2) 

 def find_s_min(S_max, N, B):
     '''Kraft, Burrows and Nousek suggest to integrate from N-B in both
     directions at once, so that S_min and S_max move similarly (see the article
     for details). Here, this is implemented differently:
     Treat S_max as the optimization parameters in func and then calculate the
     matching s_min that has has eqn7(S_max) = eqn7(S_min) here.
     '''
     y_S_max = eqn7(S_max, N, B)
     if eqn7(0, N, B) >= y_S_max:
         return 0.
     else:
         def eqn7ysmax(x):
             return eqn7(x, N, B) - y_S_max
         return findroot(eqn7ysmax , (N - B) / 2.)

def diodeResistorIMPLICITfunction(x, b, c=None):

    """
    Возвращает ток, текущей по цепи

    коэффициенты модели: Ток утечки Is, коэффициент неидеальности N, омическое сопротивление, последовательное диоду R
    входные параметры: напряжение, приложенное источником к системе резистор-диод
    +-----------|||||---------->|--------- -
    Резистор подключен до диода

    ===ОГРАНИЧЕНИЯ===
    #Сочетание боьльших напряжений  (>1В) и нулевого сопротивления даёт идиотский результат (единицу, то есть метод не отрабатывает)
    Аналогичное нехорошее поведение может повстречаться и при прочих сочетаниях. Этот момент должно иметь в виду

    :return:
    """

    global numnone
    global resultsOfEstimation
    # V=x[0] #напряжение на диоде
    # Is=b[0]
    # N=b[1]
    # R=b[2]
    FT = 0.02586419  # подогнанное по pspice

    dfdy = lambda y, x, b, c=None: np.array(
        [[-1 - b[0] * b[2] * math.exp((-b[2] * y[0] + x[0]) / (FT * b[1])) / (FT * b[1])]]
    )
    # function=lambda y,x,b,c=None: [b[0]*(math.exp((x[0]-y[0]*b[2])/(FT*b[1])) -1)-y[0]] #в подготовленном для взятия производной виде

    function = lambda y: func(y, x, b, c)

    solvinit = [1]

    #    solx=optimize.root(function, solvinit, args=(x,b,c), jac=dfdy, method='lm').x

    try:
        solx = [mpm.findroot(function, solvinit, verify=False, solver="secant", verbose=False)]
    except BaseException as e:
        print("Error in findroot=", e)
        return None

    # if solx-solvinit==[0]*len(solx):
    #     numnone+=1
    #     return None

    return solx

def getGrowthNakata(ky_list, Setup, init = -0.07 -.015j):
  
  results = []
  
  eta_e     = Setup['eta_e']
  kx        = Setup['kx']
  v_te      = Setup['v_te']
  rho_te2   = Setup['rho_te2']
  tau       = Setup['tau']
  theta     = Setup['theta']

  for ky in ky_list:
    
    kp  = mp.sqrt(2.) * theta * ky 
    # Dispersion Relation Equation (9)
    if ky > 2.5 : init = 0.01 - 0.01j
    def DispersionRelation(w):
        ko2 = kx**2 + ky**2
        def Lambda(b): return mp.exp(b) * (1. + tau - Gamma0(b))
        #def Lambda(b): return mp.exp(b) * (tau + b/(1.+b))
      
        #zeta = w / (kp * v_te)
        zeta = w / (kp * v_te)
        w_star_e = ky * pylab.sqrt(rho_te2) * v_te
 
        # Take care of extra pie due to normalization
        #return  1. + Lambda(ko2 * rho_te2) + zeta * Z(zeta) - ky/kp *  eta_e * zeta - ky/kp * ( eta_e * zeta**2 +\
        return  1. + Lambda(ko2 * rho_te2) + zeta * Z(zeta) - ky/kp *  eta_e * zeta - ky/kp * ( eta_e * zeta**2 +\
                    (1. - eta_e/2. * (1. + ko2 * rho_te2)))*Z(zeta) * mp.sqrt(mp.pi)
                    #(1. - eta_e/2. * (1. + ko2 * rho_te2)))*Z(zeta) 

    try:
        omega = complex(mp.findroot(DispersionRelation, init, solver='muller', maxsteps=1000))
        #omega = complex(PT.zermuller(DispersionRelation, 0., -0.2 - 0.05j, dx=.001)[0])

    except:
        omega = .0
        print "Not found : ", ky,  "  Theta : ", theta
    results.append(float(pylab.real(omega))  + 1.j * pylab.imag(omega))

  return (pylab.array(ky_list), pylab.array(results))

    def _eval_evalf(self, prec):
        """Evaluate this complex root to the given precision. """
        with workprec(prec):
            g = self.poly.gen
            if not g.is_Symbol:
                d = Dummy('x')
                func = lambdify(d, self.expr.subs(g, d))
            else:
                func = lambdify(g, self.expr)

            interval = self._get_interval()
            if not self.is_real:
                # For complex intervals, we need to keep refining until the
                # imaginary interval is disjunct with other roots, that is,
                # until both ends get refined.
                ay = interval.ay
                by = interval.by
                while interval.ay == ay or interval.by == by:
                    interval = interval.refine()

            while True:
                if self.is_real:
                    a = mpf(str(interval.a))
                    b = mpf(str(interval.b))
                    if a == b:
                        root = a
                        break
                    x0 = mpf(str(interval.center))
                else:
                    ax = mpf(str(interval.ax))
                    bx = mpf(str(interval.bx))
                    ay = mpf(str(interval.ay))
                    by = mpf(str(interval.by))
                    if ax == bx and ay == by:
                        # the sign of the imaginary part will be assigned
                        # according to the desired index using the fact that
                        # roots are sorted with negative imag parts coming
                        # before positive (and all imag roots coming after real
                        # roots)
                        deg = self.poly.degree()
                        i = self.index  # a positive attribute after creation
                        if (deg - i) % 2:
                            if ay < 0:
                                ay = -ay
                        else:
                            if ay > 0:
                                ay = -ay
                        root = mpc(ax, ay)
                        break
                    x0 = mpc(*map(str, interval.center))

                try:
                    root = findroot(func, x0)
                    # If the (real or complex) root is not in the 'interval',
                    # then keep refining the interval. This happens if findroot
                    # accidentally finds a different root outside of this
                    # interval because our initial estimate 'x0' was not close
                    # enough. It is also possible that the secant method will
                    # get trapped by a max/min in the interval; the root
                    # verification by findroot will raise a ValueError in this
                    # case and the interval will then be tightened -- and
                    # eventually the root will be found.
                    #
                    # It is also possible that findroot will not have any
                    # successful iterations to process (in which case it
                    # will fail to initialize a variable that is tested
                    # after the iterations and raise an UnboundLocalError).
                    if self.is_real:
                        if (a <= root <= b):
                            break
                    elif (ax <= root.real <= bx and ay <= root.imag <= by):
                        break
                except (UnboundLocalError, ValueError):
                    pass
                interval = interval.refine()

        return (Float._new(root.real._mpf_, prec)
                + I*Float._new(root.imag._mpf_, prec))