def verify_roots_of_generated_polynomial_on_unit_circle(num_of_lamda, lamda_lst):

    # we need to flip the coefficient left and right, because
    # numpy.roots function expect the coefficient in the following format:
    # p[0] * x**n + p[1] * x**(n-1) + ... + p[n-1]*x + p[n]

    # numpy.fliplr require 2D array, for 1D array, we need the follow trick
    lamda_new = np.fliplr([lamda_lst])[0]

    roots = np.roots(lamda_new)

    # get the absolute of the roots to see if they are on unit circle
    roots_abs = np.absolute(roots)

    # since leeya polynomial requires one symmetry, that is p[0] = p[n], p[1] = p[n-1], etc...
    # above polynomial will not meet this symmetry, let's just create a symmetric version

    lamda_symmetric = np.concatenate([lamda_new, lamda_lst])
    roots_symmetric = np.roots(lamda_symmetric)

    # get the absolute of the roots to see if they are on unit circle
    roots_symmetric_abs = np.absolute(roots_symmetric)

    # print "num_of_lamda", num_of_lamda,  'roots_abs =', roots_abs
    print "num_of_lamda", num_of_lamda,  'roots_symmetric_abs =', roots_symmetric_abs
    # print "num_of_lamda", num_of_lamda,  'roots_symmetric =', roots_symmetric

    return roots, roots_abs, roots_symmetric, roots_symmetric_abs

def test_calculateTF():
	"""Test function for calculateTF()"""
	from ._utils import cplxpair
	ABCD = [[1.000000000000000, 0., 0., 0.044408783846879, -0.044408783846879],
	        [0.999036450096481, 0.997109907515262, -0.005777399147297, 0., 0.499759089304780],
	        [0.499759089304780, 0.999036450096481, 0.997109907515262,  0., -0.260002096136488],
	        [0,                 0,                 1.000000000000000,  0, -0.796730400347216]]
	ABCD = np.array(ABCD)
	ntf, stf = calculateTF(ABCD)
	ntf_zeros, ntf_poles = np.roots(ntf.num), np.roots(ntf.den)
	stf_zeros, stf_poles = np.roots(stf.num), np.roots(stf.den)
	mntf_poles = np.array((1.498975311463384, 1.102565142679772, 0.132677264750882))
	mntf_zeros = np.array((0.997109907515262 + 0.075972576202904j, 
	                       0.997109907515262 - 0.075972576202904j, 
	                       1.000000000000000 + 0.000000000000000j)
	                      )
	mstf_zeros = np.array((-0.999999999999996,))
	mstf_poles = np.array((1.498975311463384, 1.102565142679772, 0.132677264750882))
	# for some reason, sometimes the zeros are in different order.
	ntf_zeros, mntf_zeros = cplxpair(ntf_zeros), cplxpair(mntf_zeros)
	stf_zeros, mstf_zeros = cplxpair(stf_zeros), cplxpair(mstf_zeros)
	ntf_poles, mntf_poles = cplxpair(ntf_poles), cplxpair(mntf_poles)
	stf_poles, mstf_poles = cplxpair(stf_poles), cplxpair(mstf_poles)
	assert np.allclose(ntf_zeros, mntf_zeros, rtol=1e-5, atol=1e-8)
	assert np.allclose(ntf_poles, mntf_poles, rtol=1e-5, atol=1e-8)
	assert np.allclose(stf_zeros, mstf_zeros, rtol=1e-5, atol=1e-8)
	assert np.allclose(stf_poles, mstf_poles, rtol=1e-5, atol=1e-8)

def daub(p):
    """
    The coefficients for the FIR low-pass filter producing Daubechies wavelets.

    p>=1 gives the order of the zero at f=1/2.
    There are 2p filter coefficients.

    Parameters
    ----------
    p : int
        Order of the zero at f=1/2, can have values from 1 to 34.

    """
    sqrt = np.sqrt
    if p < 1:
        raise ValueError("p must be at least 1.")
    if p==1:
        c = 1/sqrt(2)
        return np.array([c,c])
    elif p==2:
        f = sqrt(2)/8
        c = sqrt(3)
        return f*np.array([1+c,3+c,3-c,1-c])
    elif p==3:
        tmp  = 12*sqrt(10)
        z1 = 1.5 + sqrt(15+tmp)/6 - 1j*(sqrt(15)+sqrt(tmp-15))/6
        z1c = np.conj(z1)
        f = sqrt(2)/8
        d0 = np.real((1-z1)*(1-z1c))
        a0 = np.real(z1*z1c)
        a1 = 2*np.real(z1)
        return f/d0*np.array([a0, 3*a0-a1, 3*a0-3*a1+1, a0-3*a1+3, 3-a1, 1])
    elif p<35:
        # construct polynomial and factor it
        if p<35:
            P = [comb(p-1+k,k,exact=1) for k in range(p)][::-1]
            yj = np.roots(P)
        else:  # try different polynomial --- needs work
            P = [comb(p-1+k,k,exact=1)/4.0**k for k in range(p)][::-1]
            yj = np.roots(P) / 4
        # for each root, compute two z roots, select the one with |z|>1
        # Build up final polynomial
        c = np.poly1d([1,1])**p
        q = np.poly1d([1])
        for k in range(p-1):
            yval = yj[k]
            part = 2*sqrt(yval*(yval-1))
            const = 1-2*yval
            z1 = const + part
            if (abs(z1)) < 1:
                z1 = const - part
            q = q * [1,-z1]

        q = c * np.real(q)
        # Normalize result
        q = q / np.sum(q) * sqrt(2)
        return q.c[::-1]
    else:
        raise ValueError("Polynomial factorization does not work "
              "well for p too large.")

    def splitBimodal(self, x, y, largepoly=30):
        p = np.polyfit(x, y, largepoly) # polynomial coefficients for fit

        extrema = np.roots(np.polyder(p))
        extrema = extrema[np.isreal(extrema)]
        extrema = extrema[(extrema - x[1]) * (x[-2] - extrema) > 0] # exclude the endpoints due false maxima during fitting
        try:
            root_vals = [sum([p[::-1][i]*(root**i) for i in range(len(p))]) for root in extrema]
            peaks = extrema[np.argpartition(root_vals, -2)][-2:] # find two peaks of bimodal distribution

            mid, = np.where((x - peaks[0])* (peaks[1] - x) > 0)
             # want data points between the peaks
        except:
            warnings.warn("Peak finding failed!")
            return None

        try:
            p_mid = np.polyfit(x[mid], y[mid], 2) # fit middle section to a parabola
            midpoint = np.roots(np.polyder(p_mid))[0]
        except:
            warnings.warn("Polynomial fit between peaks of distribution poorly conditioned. Falling back on using the minimum! May result in inaccurate split determination.")
            if len(mid) == 0:
                return None

            midx = np.argmin(y[mid])
            midpoint = x[mid][midx]

        return midpoint

    def setUp(self):
        ABCD = [[1.0, 0.0, 0.0, 0.044408783846879, -0.044408783846879],
                [0.999036450096481, 0.997109907515262, -0.005777399147297,
                 0.0, 0.499759089304780],
                [0.499759089304780, 0.999036450096481, 0.997109907515262,
                 0.0, -0.260002096136488],
                [0.0, 0.0, 1.0,  0.0, -0.796730400347216]]
        ABCD = np.array(ABCD)
        (ntf, stf) = ds.calculateTF(ABCD)
        (ntf_zeros, ntf_poles) = (np.roots(ntf.num), np.roots(ntf.den))
        (stf_zeros, stf_poles) = (np.roots(stf.num), np.roots(stf.den))
        mntf_poles = np.array((1.498975311463384, 1.102565142679772,
                               0.132677264750882))
        mntf_zeros = np.array((0.997109907515262 + 0.075972576202904j,
                               0.997109907515262 - 0.075972576202904j,
                               1.000000000000000 + 0.000000000000000j))
        mstf_zeros = np.array((-0.999999999999996,))
        mstf_poles = np.array((1.498975311463384, 1.102565142679772,
                               0.132677264750882))

        # for some reason, sometimes the zeros are in different order.
        (self.ntf_zeros, self.mntf_zeros) = (cplxpair(ntf_zeros),
                                             cplxpair(mntf_zeros))
        (self.stf_zeros, self.mstf_zeros) = (cplxpair(stf_zeros),
                                             cplxpair(mstf_zeros))
        (self.ntf_poles, self.mntf_poles) = (cplxpair(ntf_poles),
                                             cplxpair(mntf_poles))
        (self.stf_poles, self.mstf_poles) = (cplxpair(stf_poles),
                                             cplxpair(mstf_poles))

def real_roots(poly_coeffs, x1=None, x2=None, warn=False):

# Supress warnings (default)
     if (warn==False):
          warnings.simplefilter('ignore', np.RankWarning)

# Evaluate roots, keeping only the real parts or those without an imaginary component
     re_roots = \
         np.roots(poly_coeffs)[np.roots(poly_coeffs).imag == 0.].real

# Go through limit possibilities, returning the appropriate values
# If no limits were given then return all real roots
     if (x1==None and x2==None):
          return re_roots
# The following are cases where either or both limits are given
     elif (x2==None):  # If only lower limit was given
          return re_roots[(re_roots >= x1)]
     elif (x1==None):  # If only upper limit was given
          return re_roots[(re_roots <= x2)]
     else:             # If both limits are given
          # Check that x1 < x2 and fix if necessary
          if (x1 > x2):
               temp_x = x1
               x1 = x2
               x2 = temp_x
          return re_roots[(re_roots >= x1) & (re_roots <= x2)]

def tf2zpk(b, a):
    """Return zero, pole, gain (z,p,k) representation from a numerator,
    denominator representation of a linear filter.

    Parameters
    ----------
    b : ndarray
        Numerator polynomial.
    a : ndarray
        Denominator polynomial.

    Returns
    -------
    z : ndarray
        Zeros of the transfer function.
    p : ndarray
        Poles of the transfer function.
    k : float
        System gain.

    Notes
    -----
    If some values of `b` are too close to 0, they are removed. In that case,
    a BadCoefficients warning is emitted.

    """
    b, a = normalize(b, a)
    b = (b + 0.0) / a[0]
    a = (a + 0.0) / a[0]
    k = b[0]
    b /= b[0]
    z = roots(b)
    p = roots(a)
    return z, p, k

def isstable(b,a,ftype='digital'):
    """Determine whether IIR filter (b,a) is stable
    
    Parameters
    ----------
        b: ndarray
            filter numerator coefficients
        a: ndarray
            filter denominator coefficients
        ftype: string
            type of filter (`digital` or `analog`)
    Returns
    -------
        stable: bool
            whether filter is stable or not
                
    """

    if ftype=='digital':
        v = np.roots(a)
        if np.any(np.abs(v)>1.0):
            return False
        else:
            return True
    elif ftype=='analog':
        v = np.roots(a)
        if np.any(np.real(v)<0):
            return False
        else:
            return True

def lpc_to_lsf(all_lpc):
    if len(all_lpc.shape) < 2:
        all_lpc = all_lpc[None]
    order = all_lpc.shape[1] - 1
    all_lsf = np.zeros((len(all_lpc), order))
    for i in range(len(all_lpc)):
        lpc = all_lpc[i]
        lpc1 = np.append(lpc, 0)
        lpc2 = lpc1[::-1]
        sum_filt = lpc1 + lpc2
        diff_filt = lpc1 - lpc2

        if order % 2 != 0:
            deconv_diff, _ = sg.deconvolve(diff_filt, [1, 0, -1])
            deconv_sum = sum_filt
        else:
            deconv_diff, _ = sg.deconvolve(diff_filt, [1, -1])
            deconv_sum, _ = sg.deconvolve(sum_filt, [1, 1])

        roots_diff = np.roots(deconv_diff)
        roots_sum = np.roots(deconv_sum)
        angle_diff = np.angle(roots_diff[::2])
        angle_sum = np.angle(roots_sum[::2])
        lsf = np.sort(np.hstack((angle_diff, angle_sum)))
        if len(lsf) != 0:
            all_lsf[i] = lsf
    return np.squeeze(all_lsf)

def fit_edge_hist(bins, counts, fwhm_guess=10.0):
    if len(bins) == len(counts)+1: bins = bins[:-1]+0.5*(bins[1]-bins[0]) # convert bin edge to bin centers if neccesary
    pfit = np.polyfit(bins, counts, 3)
    edgeGuess = np.roots(np.polyder(pfit,2))
    try:
        preGuessX, postGuessX = np.sort(np.roots(np.polyder(pfit,1)))
    except:
        raise ValueError("failed to generate guesses")
    use = bins>(edgeGuess+2*fwhm_guess)
    if np.sum(use)>4:
        pfit2 = np.polyfit(bins[use], counts[use],1)
        slope_guess = pfit2[0]
    else:
        slope_guess=1
    pGuess = np.array([edgeGuess, np.polyval(pfit,preGuessX), np.polyval(pfit,postGuessX),fwhm_guess,slope_guess],dtype='float64')

    try:
        pOut = curve_fit(edge_model, bins, counts, pGuess)
    except:
        return (0,0,0,0,0,0)
    (edgeCenter, preHeight, postHeight, fwhm, bgSlope) = pOut[0]
    model_counts = edge_model(bins, edgeCenter, preHeight, postHeight, fwhm, bgSlope)
    num_degree_of_freedom = float(len(bins)-1-5) # num points - 1 - number of fitted parameters
    chi2 = np.sum(((counts - model_counts)**2)/model_counts)/num_degree_of_freedom
    return (edgeCenter, preHeight, postHeight, fwhm, bgSlope, chi2)

def lsf(fir_filt):
  """
  Find the Line Spectral Frequencies (LSF) from a given FIR filter.

  Parameters
  ----------
  filt :
    A LTI FIR filter as a LinearFilter object.

  Returns
  -------
  A tuple with all LSFs in rad/sample, alternating from the forward prediction
  and backward prediction filters, starting with the lowest LSF value.

  """
  den = fir_filt.denominator
  if len(den) != 1:
    raise ValueError("Filter has feedback")
  elif den[0] != 1: # So we don't have to worry with the denominator anymore
    fir_filt /= den[0]

  from numpy import roots
  rev_filt = ZFilter(fir_filt.numerator[::-1]) * z ** -1
  P = fir_filt + rev_filt
  Q = fir_filt - rev_filt
  roots_p = roots(P.numerator[::-1])
  roots_q = roots(Q.numerator[::-1])
  lsf_p = sorted(phase(roots_p))
  lsf_q = sorted(phase(roots_q))
  return reduce(operator.concat, xzip(*sorted([lsf_p, lsf_q])), tuple())

def polyxval(poly_coeffs, y, x1=None, x2=None, warn=False):
     if (x1==None or x2==None):
          print "Must assign range [x1, x2] in which to search for x values."
          return None

     if(x1==x2):
          print "x1 must not equal x2."
          return None

     if (x1 > x2):
          temp_x = x1
          x1 = x2
          x2 = temp_x

# Supress warnings (default)
     if (warn==False):
          warnings.simplefilter('ignore', np.RankWarning)

     poly_coeffs_y = poly_coeffs
     # subtract y-value from zeroth order coefficient (i.e. no x's)
     poly_coeffs_y[len(poly_coeffs_y)-1] -= y
     re_roots = \
         np.roots(poly_coeffs_y)[np.roots(poly_coeffs_y).imag == 0.].real

# restrict solution to range [x1, x2]
     x_val = re_roots[(re_roots >= x1) & (re_roots <= x2)]

     return x_val

def plot_range(num, den):

    # The corner frequencies
    zero = sort(abs(roots(num)))    
    pole = sort(abs(roots(den)))
    
    # Calculate the minimum and maximum corner frequencies needed
    if len(pole) == 0:
        corner_min = zero[0]
        corner_max = zero[-1]
    
    elif len(zero) == 0:
        corner_min = pole[0]
        corner_max = pole[-1]
        
    elif len(zero) > 0 and len(pole) > 0:
        corner_min = min(zero[0], pole[0])
        corner_max = max(zero[-1], pole[-1]) 
    
    else:
        corner_min, corner_max = 0.1, 10
    
    # start from 2 decades lower than the lowest corner 
    # and end at 2 decades above the highest corner
    freq_range = [10 ** (floor(log10(corner_min)) - 1), 
                   10 ** (floor(log10(corner_max)) + 2)]
    
    return freq_range

def ARLineSpectra(ar):
    """
    Convert AR coeffs to LSPs

    From wikipedia:
    A palindromic polynomial (i.e., P) of odd degree has -1 as a root.
    An antipalindromic polynomial (i.e., Q) has 1 as a root.
    An antipalindromic polynomial of even degree has -1 and 1 as roots
    """
    order = ar.shape[-1]
    ret = np.zeros(ar.shape)
    for a, o in core.refiter([ar, ret], core.newshape(ar.shape)):
        p = np.ones((order+2))
        q = np.ones((order+2))
        q[-1] = -1.0
        for i in range(order):
            p[i+1] = -a[i] - a[order-i-1]
            q[i+1] = -a[i] + a[order-i-1]
        pr = np.roots(p)
        qr = np.roots(q)

        j = 0
        an = np.ndarray((order+2))
        for i in range(len(pr)):
            if np.imag(pr[i]) >= 0.0:
                an[j] = np.angle(pr[i])
                j += 1
            if np.imag(qr[i]) >= 0.0:
                an[j] = np.angle(qr[i])
                j += 1
        # The angle list (an) will always contain both 0 and pi; they
        # will move to the ends after the sort
        o[...] = np.sort(an)[1:-1]
    return ret;

def asymptote(num, den):

    # create a Python list for the zeros and the poles of the system
    zero = list(sort(abs(roots(num))))    
    pole = list(sort(abs(roots(den))))
    
    #calculate the low frequency gain -- type 0 system
    lf_gain = 20 * log10(abs((num[-1] / den[-1])))
    
    # create an empty matrix to contain the corner frequencies and
    # the corresponding slope indicator (+1 or -1)
    corners = zeros((len(zero) + len(pole) + 2, 2))
    
    starting_freq, end_freq = plot_range(num, den)
    corners[0] = [starting_freq, 0]
    corners[-1] = [end_freq, 0]
    
    # take the first elements from the list of poles and zeros
    # compare them and assign the slope indicator
    # delete the corresponding valuefrom the original list of poles and zeros
    for count in range(len(zero) + len(pole)): 
            
        if len(zero) > 0:
            a = zero[0]
        else:
            a = inf
            
        if len(pole) > 0:
            b = pole[0]
        else:
            b = inf
        
        c = min(a, b)
                
        if c == a:
            corners[count + 1] = [c, 1]
            if len(zero) > 0:
                zero.pop(0)

        if c == b:
            corners[count + 1] = [c, -1]
            if len(pole) > 0:
                pole.pop(0)
         
    # now calculate the gains at the corners using 
    # gain = +/- 20log10(upper_corner / lower_corner) 
    asymptotic_gain = zeros_like(corners)   
    asymptotic_gain[0, 1] = lf_gain
    asymptotic_gain[1, 1] = lf_gain
    
    gain = lf_gain
    multiplier = cumsum(corners[:, 1])
    for k in range(2, len(corners)):    
        gain += multiplier[k-1] * 20 * log10(corners[k, 0] / corners[k-1, 0])
        asymptotic_gain[k, 1] = gain
        
    asymptotic_gain[:, 0] = corners[:, 0]
                   
    return asymptotic_gain

	def free_ene(self, trange):
		self.get_data()
		#self.get_E0()
		#ax = plt.subplot(111)
		ndic = collections.OrderedDict(sorted(self.__dic.items()))	#sort dictionary
		minF = []
		minl = []
		for temp in trange:
			xdata = []
			ydata = []
			i=0
			for out in ndic:
				xdata.append(out)
				ind = self.__dic[out]['T'].index(temp)
				
				ydata.append(self.__dic[out]['F'][ind]/self.conv + self.__E0[i] + 13.5)
				i+=1
			#polyfit:
			coeff = np.polyfit(xdata,ydata,3)
			p = np.poly1d(coeff)
			polyx = np.linspace(min(xdata),max(xdata),1000)
			
			#ax.plot(xdata,ydata,'+')
			#ax.plot(polyx,p(polyx))
			minl.append(np.roots(p.deriv())[1])
			minF.append(p(np.roots(p.deriv())[1]))
			
			#polyfit F-T
			coeff = np.polyfit(minl,minF,21)
			p = np.poly1d(coeff)
			polyx = np.linspace(min(minl),max(minl),1000)
			
			#ax.plot(polyx,p(polyx))
			#ax.plot(minl,minF,'o')
		
		
		
		
		#polyfit thermal expansion:
		
		coeff = np.polyfit(trange,minl,4)
		p = np.poly1d(coeff)
		polyx = np.linspace(min(trange),max(trange),1000)
		
		#ax1 = plt.plot(polyx,p(polyx))
		#ax1 = plt.plot(trange,minl,'o')
		
		#thermal expansion/T
		
		alpha = []
		for i in range(len(trange)-1):
			alpha.append((minl[i+1]-minl[i])/(trange[i+1]-trange[i])/minl[0]*10**6.)
		plt.plot(trange[:-1], alpha, label = self._lname, lw=2., ls=self.style)
		self.__alpha = alpha
		self.__minl = minl
		self.__minF = minF

def InitializeCalculations(output):
    """This functions calculates parameters from the input"""
    # SYSTEM
    if output['series'] is False:	
    	output['system'] = 'F' + "{0:.2f}".format(output['BASEfraction']).replace('.','_') + 'D1_8V' + \
                           '0'*(len(str(output['version']))%2) + str(output['version'])
        output['systemfolder'] = output['system']
    # BASE PARTICLES	
    output['BASEparticlevolume'] = (4.0/3.0) * np.pi * (output['BASEsigmacorrected']/2.0)**3	
    output['BASEdiffusion'] = (output['ReferenceT']*output['Boltzmann']) / (6.0*np.pi*output['viscosity']*(0.5*output['BASEsigma']))
    output['BASEtau'] = (output['BASEsigma']**2) / output['BASEdiffusion']	
    output['eta'] = output['viscosity']*(output['BASEsigma']**3)/(output['Boltzmann']*output['ReferenceT']*output['BASEtau'])
    output['BASEfriction'] = 0.5*6.0*output['BASEsigmacorrected']*np.pi*output['eta'];
    if output['series'] is False:
        if output['basegenerator'] == 'Lbox' or output['InitialConfiguration'] == 'random':
            output['BOXvolume'] = output['BOXsize_x']*output['BOXsize_y']*output['BOXsize_z']
            TotalParticleVolume = output['BOXvolume'] * output['BASEfraction']
            output['BASEparticles'] = int(np.round(TotalParticleVolume / output['BASEparticlevolume'],0))
        elif output['InitialConfiguration'] == 'BCC' or output['InitialConfiguration'] == 'FCC':	
            output['BOXsize'] = (output['BASEparticlevolume']*output['BASEparticles']/output['BASEfraction'])**(1/3.0)
            output['BOXsize_x'] = output['BOXsize']
            output['BOXsize_y'] = output['BOXsize']
            output['BOXsize_z'] = output['BOXsize']
            output['BOXvolume'] = output['BOXsize_x']*output['BOXsize_y']*output['BOXsize_z']			
    if output['InitialConfiguration'] == 'BCC' and output['basegenerator'] == 'Nparticles':
        roots = np.roots([2,0,0,-1.0*output['BASEparticles']])[::-1] # Nparticles of the basis crystal along one axis
        m = roots[0]
        output['m'] = int(round(m.real,0))
        output['n'] = output['m']-1 # number of unit cells
        output['a'] = output['BOXsize'] / output['m']
    elif output['InitialConfiguration'] == 'FCC': 
        roots = np.roots([4,-6,3,-output['BASEparticles']]) # number of particles of the basis crystal along one axis
        m = roots[0]
        output['m'] = int(round(m.real,0))
        output['n'] = output['m']-1 # number of unit cells
        output['a'] = output['BOXsize'] / output['m']
    elif output['InitialConfiguration'] == 'BCC' and output['basegenerator'] == 'Lbox':
        # for a box with ratio's x,y = 1 and z = 0.5
        output['m'] = int(np.round(np.power(output['BASEparticles'],(1/3.0)),0))
        output['n'] = output['m']-1 # number of unit cells
        output['a'] = output['BOXsize'] / output['m']
    if output['method'] == 'wigner':
        output['ALLparticles'] = output['BASEparticles'] 
    # INTERSTITIALS
    elif output['method'] == 'interstitial':		
        output['INTERsigmacorrected'] = output['INTERsigma'] / output['BASEsigma']	
        output['INTERfriction'] = 0.5*6.0*output['INTERsigmacorrected']*np.pi*output['eta']
        output['INTERdiffusion'] = (output['ReferenceT']*output['Boltzmann']) / (6.0*np.pi*output['viscosity']*(0.5*output['INTERsigma']))
        output['INTERtau'] = (output['INTERsigma']**2) / output['INTERdiffusion']
        output['IBepsilon'] = 0.5*(output['BASEepsilon'] + output['INTERepsilon'])
        output['IBkappa'] = 0.5*(output['BASEkappa'] + output['INTERkappa'])
        if output['InitialConfiguration'] == 'BCC':
            output['INTERsites'] = 12*output['m']*output['n']*output['n']
            output['INTERparticles'] = int(output['INTERsites']*output['INTERfraction'])
        output['ALLparticles'] = output['BASEparticles'] + output['INTERparticles']
    return output

def dual_band(F, P, E, eps, eps_R=1, x1=0.5, map_type=1):
    r"""
    Function to give modified F, P, and E polynomials after lowpass 
    prototype dual band transformation.
    
    :param F:       Polynomial F, i.e., numerator of S11 (in s-domain)
    :param P:       Polynomial P, i.e., numerator of S21 (in s-domain)
    :param E:       polynomial E is the denominator of S11 and S21 (in s-domain)
    :param eps:     Constant term associated with S21
    :param eps_R:   Constant term associated with S11
    :param x1:
    :param map_type:
                  
    :rtype:    
    """

    if (map_type == 1) or (map_type == 2):

        s = sp.Symbol("s")
        if map_type == 1:
            a = -2j / (1 - x1 ** 2)
            b = -1j * (1 + x1 ** 2) / (1 - x1 ** 2)
        elif map_type == 2:
            a = 2j / (1 - x1 ** 2)
            b = 1j * (1 + x1 ** 2) / (1 - x1 ** 2)
        s1 = a * s ** 2 + b

        F = sp.Poly(F.ravel().tolist(), s)
        F1 = sp.simplify(F.subs(s, s1))
        F1 = sp.Poly(F1, s).all_coeffs()
        F1 = I_to_i(F1)

        E = sp.Poly(E.ravel().tolist(), s)
        E1 = sp.simplify(E.subs(s, s1))
        E1 = sp.Poly(E1, s).all_coeffs()
        E1 = I_to_i(E1)

        P = sp.Poly(P.ravel().tolist(), s)
        P1 = sp.simplify(P.subs(s, s1))
        P1 = sp.Poly(P1, s).all_coeffs()
        P1 = I_to_i(P1)

    elif map_type == 3:

        F_roots = np.roots(F.ravel().tolist())
        P_roots = np.roots(P.ravel().tolist())
        F1_roots = Lee_roots_map(F_roots, x1)
        P1_roots = Lee_roots_map(P_roots, x1)
        P1_roots = np.concatenate((P1_roots, np.array([0, 0])))
        F1 = np.poly(F1_roots)
        P1 = np.poly(P1_roots)
        F1 = np.reshape(F1, (len(F1), -1))
        P1 = np.reshape(P1, (len(P1), -1))
        E1 = poly_E(eps, eps_R, F1, P1)[0]

    return F1, P1, E1

def poly2lsf(a):
    """Prediction polynomial to line spectral frequencies.

    converts the prediction polynomial specified by A,
    into the corresponding line spectral frequencies, LSF.
    normalizes the prediction polynomial by A(1).

    .. doctest::

        >>> from spectrum import poly2lsf
        >>> a = [1.0000,  0.6149, 0.9899, 0.0000 ,0.0031, -0.0082]
        >>> lsf = poly2lsf(a)
        >>> lsf =  array([0.7842, 1.5605, 1.8776, 1.8984, 2.3593])

    .. seealso:: lsf2poly, poly2rc, poly2qc, rc2is
    """

    #Line spectral frequencies are not defined for complex polynomials.

    # Normalize the polynomial

    a = numpy.array(a)
    if a[0] != 1:
        a/=a[0]

    if max(numpy.abs(numpy.roots(a))) >= 1.0:
        error('The polynomial must have all roots inside of the unit circle.');


    # Form the sum and differnce filters

    p  = len(a)-1   # The leading one in the polynomial is not used
    a1 = numpy.concatenate((a, numpy.array([0])))
    a2 = a1[-1::-1]
    P1 = a1 - a2        # Difference filter
    Q1 = a1 + a2        # Sum Filter

    # If order is even, remove the known root at z = 1 for P1 and z = -1 for Q1
    # If odd, remove both the roots from P1

    if p%2: # Odd order
        P, r = deconvolve(P1,[1, 0 ,-1])
        Q = Q1
    else:          # Even order
        P, r = deconvolve(P1, [1, -1])
        Q, r = deconvolve(Q1, [1,  1])

    rP  = numpy.roots(P)
    rQ  = numpy.roots(Q)

    aP  = numpy.angle(rP[1::2])
    aQ  = numpy.angle(rQ[1::2])

    lsf = sorted(numpy.concatenate((-aP,-aQ)))

    return lsf

def zplane(b,a,filename=None):
    """Plot the complex z-plane given a transfer function.
    """

    # get a figure/plot
    ax = plt.subplot(111)

    # create the unit circle
    uc = patches.Circle((0,0), radius=1, fill=False,
                        color='black', ls='dashed')
    ax.add_patch(uc)

    # The coefficients are less than 1, normalize the coeficients
    if np.max(b) > 1:
        kn = np.max(b)
        b = b/float(kn)
    else:
        kn = 1

    if np.max(a) > 1:
        kd = np.max(a)
        a = a/float(kd)
    else:
        kd = 1
        
    # Get the poles and zeros
    p = np.roots(a)
    z = np.roots(b)
    k = kn/float(kd)
    
    # Plot the zeros and set marker properties    
    t1 = plt.plot(z.real, z.imag, 'go', ms=10)
    plt.setp( t1, markersize=10.0, markeredgewidth=1.0,
              markeredgecolor='k', markerfacecolor='g')

    # Plot the poles and set marker properties
    t2 = plt.plot(p.real, p.imag, 'rx', ms=10)
    plt.setp( t2, markersize=12.0, markeredgewidth=3.0,
              markeredgecolor='r', markerfacecolor='r')

    ax.spines['left'].set_position('center')
    ax.spines['bottom'].set_position('center')
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)

    # set the ticks
    r = 1.5; plt.axis('scaled'); plt.axis([-r, r, -r, r])
    ticks = [-1, -.5, .5, 1]; plt.xticks(ticks); plt.yticks(ticks)

    if filename is None:
        plt.show()
    else:
        plt.savefig(filename)
    

    return z, p, k