def genCheby(self):
		# find order first
		n = np.ceil((np.arccosh(np.sqrt(self.sat-1)/self.epsilon))/(np.arccosh(self.ws/self.w1)));

		#find coefficients
		beta = np.sinh(np.arctanh(1/np.sqrt(1+self.epsilon**2))/n)

		cbk = [1];
		for i in range(1,len(str(n))+1):
			cbk.append(2*np.sin(((2*i-1)*np.pi)/(2*n)));

		ck = np.ones(len(str(n))+1);

		for i in range(1,len(cbk)):
			if (i == 1):
				ck[i] = cbk[i]/beta;
			else: 
				ck[i] = (cbk[i]*cbk[i-1])/(ck[i-1]*(beta**2+np.sin((i-1)*np.pi/n)**2));

		if (self.firsty == "l"):
			for i in range(1, len(str(n))+1):
				if (i%2 == 1):
					tmp = "L" + str(i)
					self.components.append(Component(tmp, "inductor", self.zin*ck[i]/self.w1));
				elif (i%2 == 0):
					tmp = "C" + str(i)
					self.components.append(Component(tmp, "capacitor", ck[i]/(self.w1*self.zin)));
		elif (self.firsty == "c"):
			for i in range(1, len(str(n))+1):
				if (i%2 == 1):
					tmp = "C" + str(i)
					self.components.append(Component(tmp, "capacitor", ck[i]/(self.w1*self.zin)));
				elif (i%2 == 0):
					tmp = "L" + str(i)
					self.components.append(Component(tmp, "inductor", self.zin*ck[i]/self.w1));

def generate_dolph_chebyshev(lobe_fraction, tolerance):
    def cheb(m, x):
        if np.abs(x) <= 1:
            return np.cos(m * np.arccos(x))
        else:
            return np.cosh(m * np.arccosh(np.abs(x))).real

    w = int((1 / np.pi) * (1 / lobe_fraction) * np.arccosh(1 / tolerance))

    if w % 2 == 0:
        w -= 1

    beta = np.cosh(np.arccosh(1 / tolerance) / float(w - 1))

    print("lobe_fraction = {0}, tolerance = {1}, w = {2}, beta = {3}".format(lobe_fraction, tolerance, w, beta))

    x = np.empty(w, dtype=np.complex128)

    for i in xrange(w):
        x[i] = cheb(w - 1, beta * np.cos(np.pi * i / w)) * tolerance

    x = fft(x, n=w)
    x = fftshift(x)
    x = np.real(x)

    return {"x": x, "w": w}

 def test_arccosh(self):
     q = [1, 2, 3, 4, 6] * self.ureg.dimensionless
     x = np.ones((1,1))  * self.ureg.rad
     self.assertEqual(
         np.arccosh(q),
         np.arccosh(q.magnitude * self.ureg.rad)
     )

def AF_zeros(a, M, R, dist_type, nbar=False, alpha=0):
    r"""
    This function gives array-factor zeros corresponding to different
    types of array distributions. Unless you know what you are doing exactly,
    do not use this function directly. Instead, user can use the function :func:`dist`.
    
    :param a:        separation between the elements along the x-axis in wavelengths
    :param M:        number of elements along the x-axis
    :param R:        side-lobe ratio in linear scale
    :param dist_type:type of the distribution, e.g., 'Dolph-Chebyshev', 'Riblet', etc.
    :param nbar:     transition index for dilation
    :param alpha:    Taylor's asymptotic tapering parameter
    
    :rtype:          U0, a Numpy array of size (*,1)
    """
    k = 2 * np.pi  # (angular) wave-number, which is 2*pi when lambda = 1
    m = np.ceil((M - 2) / 2)
    n = np.arange(1, 1 + m, 1)  # number of zeros for symmetric array-factors
    na = np.arange(1, M, 1)  # number of zeros for 'asymmetric' array-factors
    
    if(dist_type == "Dolph-Chebyshev"):  # Dolph zeros
        c = np.cosh(np.arccosh(R) / (M - 1))
        print c
        U0 = (2 / (a * k)) * np.arccos((np.cos(np.pi * (2 * n - 1) / (2 * M - 2))) / c)        
    elif(dist_type == "Riblet"):  # Riblet zeros
        c1 = np.cosh(np.arccosh(R) / m)
        c = np.sqrt((1 + c1) / (2 + (c1 - 1) * np.cos(k * a / 2) ** 2))
        alph = c * np.cos(k * a / 2)
        xi = (1 / c) * np.sqrt(((1 + alph ** 2) / 2) + ((1 - alph ** 2) / 2) * 
                               np.cos(((2 * n - 1) * np.pi) / (2 * m)))
        U0 = (2 / (a * k)) * np.arccos(xi)
    elif(dist_type == "Duhamel-b"):  # Duhamel bi-directional end-fire array zeros
        if(a < 0.5): c = np.cosh(np.arccosh(R) / (M - 1)) / np.sin((k * a) / 2)
        else: c = np.cosh(np.arccosh(R) / (M - 1))
        U0 = (2 / (a * k)) * np.arcsin((np.cos(np.pi * (2 * n - 1) / (2 * M - 2))) / c)
    elif(dist_type == "Duhamel-u"):  # Duhamel uni-directional end-fire array zeros
        Lamb = np.cosh(np.arccosh(R) / (M - 1))
        xi = (2 / a) * (0.5 * np.pi - (np.arctan(np.tan(k * a / 2) * ((Lamb + 1) / (Lamb - 1)))))
        c = 1 / (np.sin((xi - k) * a / 2))
        U0 = -(xi / k) + (2 / (a * k)) * np.arcsin((
                            np.cos(np.pi * (2 * na - 1) / (2 * M - 2))) / c)
    elif(dist_type == "McNamara-s"):  # McNamara-Zolotarev sum-pattern zeros
        U0 = "Yet to be done"
    elif(dist_type == "McNamara-d"):  # McNamara-Zolotarev difference-pattern zeros
        if(a < 0.5): c = 1 / np.sin((k * a) / 2)
        else: c = 1
        m1 = Zol.z_m_frm_R(M - 1, R)
        xn = Zol.z_Zolotarev_poly(N=M - 1, m=m1)[1][m + 1:]
        U0 = (2 / (a * k)) * np.arcsin(xn / c)
    if(nbar):  # Taylor's Dilation procedure
        if((dist_type == "Dolph-Chebyshev") or (dist_type == "Riblet") or (dist_type == "McNamara-s")):
            n_gen = np.arange(nbar, 1 + m, 1)  # indices of the generic zeros
            U0_gen = (n_gen + alpha / 2) * (1 / (M * a))  # generic sum zeros
        elif(dist_type == "McNamara-d"):
            n_gen = np.arange(nbar, 1 + m, 1)  # indices of the generic zeros
            U0_gen = (n_gen + (alpha + 1) / 2) * (1 / (M * a))  # generic difference zeros
        sigma = U0_gen[0] / U0[nbar - 1]  # Dilation factor
        U0 = np.hstack((sigma * U0[0:nbar - 1], U0_gen))  # Dilated zeros
    U0 = np.reshape(U0, (len(U0), -1))
    return U0

 def geodesic(self, P1, P2, theta):
     if np.ndim(P1)==1:
          arg =np.array([np.arccosh(2*P1[0]*P2[0]-np.inner(P1, P2))])
     else:
         arg = np.arccosh(2*P1[:,0]*P2[:,0]-np.einsum('ij...,ij...->i...',P1, P2))
         arg = arg[:, np.newaxis,...]
     return ((1-theta)*sinhc((1-theta)*arg)*P1 + theta*sinhc(theta*arg)*P2)/sinhc(arg)

    def test_arccosh(self):
        import math
        from numpy import arccosh

        for v in [1.0, 1.1, 2]:
            assert math.acosh(v) == arccosh(v)
        for v in [-1.0, 0, 0.99]:
            assert math.isnan(arccosh(v))

def localized(T, a, xi):
    '''
    total transmission distribution for 1D localized system;
    effective localization length also fit parameter (in units
    of L)
    '''
    return (a * np.sqrt(np.arccosh(T**(-1./2.))) / (T**(3./2.) * (1.-T)**(1./4.)) \
        * np.exp(-0.5 * xi * np.arccosh(T**(-1./2.))**2.)).real

def chebwin(M, at, sym=True):
    """Return a Dolph-Chebyshev window.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero or less, an empty
        array is returned.
    at : float
        Attenuation (in dB).
    sym : bool, optional
        When True, generates a symmetric window, for use in filter design.
        When False, generates a periodic window, for use in spectral analysis.

    Returns
    -------
    w : ndarray
        The window, with the maximum value always normalized to 1

    """
    if M < 1:
        return np.array([])
    if M == 1:
        return np.ones(1, 'd')

    odd = M % 2
    if not sym and not odd:
        M = M + 1

    # compute the parameter beta
    order = M - 1.0
    beta = np.cosh(1.0 / order * np.arccosh(10 ** (np.abs(at) / 20.)))
    k = np.r_[0:M] * 1.0
    x = beta * np.cos(np.pi * k / M)
    # Find the window's DFT coefficients
    # Use analytic definition of Chebyshev polynomial instead of expansion
    # from scipy.special. Using the expansion in scipy.special leads to errors.
    p = np.zeros(x.shape)
    p[x > 1] = np.cosh(order * np.arccosh(x[x > 1]))
    p[x < -1] = (1 - 2 * (order % 2)) * np.cosh(order * np.arccosh(-x[x < -1]))
    p[np.abs(x) <= 1] = np.cos(order * np.arccos(x[np.abs(x) <= 1]))

    # Appropriate IDFT and filling up
    # depending on even/odd M
    if M % 2:
        w = np.real(fft(p))
        n = (M + 1) / 2
        w = w[:n] / w[0]
        w = np.concatenate((w[n - 1:0:-1], w))
    else:
        p = p * np.exp(1.j * np.pi / M * np.r_[0:M])
        w = np.real(fft(p))
        n = M / 2 + 1
        w = w / w[1]
        w = np.concatenate((w[n - 1:0:-1], w[1:n]))
    if not sym and not odd:
        w = w[:-1]
    return w

def acosh(x):
    """
    Inverse hyperbolic cosine
    """
    if isinstance(x, UncertainFunction):
        mcpts = np.arccosh(x._mcpts)
        return UncertainFunction(mcpts)
    else:
        return np.arccosh(x)

def test_arccosh():
    # Domain for arccosh starts at 1
    a = afnumpy.random.random((2,3))+1
    b = numpy.array(a)
    fassert(afnumpy.arccosh(a), numpy.arccosh(b))
    c = afnumpy.random.random((2,3))
    d = numpy.array(a)
    fassert(afnumpy.arccosh(a, out=c), numpy.arccosh(b, out=d))
    fassert(c, d)

def localized_log(T, a, xi):
    '''
    total transmission distribution for 1D localized system;
    distribution for lnT ( P(lnT) = T * P(T); prefactor ln(10), which can
    be considered as contained in amplitude a
    effective localization length also fit parameter (in units
    of L)
    '''
    return (a * T * np.sqrt(np.arccosh(T**(-1./2.))) / (T**(3./2.) * (1.-T)**(1./4.)) \
        * np.exp(-0.5 * xi * np.arccosh(T**(-1./2.))**2.)).real

    def test_der_arccosh():
        x = np.linspace(1.2, 5, 5)
        h = 1e-8
        der1 = np.arccosh(bicomplex(x + h * 1j, 0)).imag1 / h
        np.testing.assert_allclose(der1, 1. / np.sqrt(x**2 - 1))

        h = (_default_base_step(x, scale=2.5) + 1) - 1
        der2 = np.arccosh(bicomplex(x + h * 1j, h)).imag12 / h**2
        true_der2 = -x / (x**2-1)**(3. / 2)
        np.testing.assert_allclose(der2, true_der2, atol=1e-5)

def arccosh(x, out=None):
    """
    Raises a ValueError if input cannot be rescaled to a dimensionless
    quantity.
    """
    if not isinstance(x, Quantity):
        return np.arccosh(x, out)

    return Quantity(
        np.arccosh(x.rescale(dimensionless).magnitude, out),
        dimensionless,
        copy=False
    )

def chebwin(M, at, sym=1):
    """Dolph-Chebyshev window.

    INPUTS:

      M : int
        Window size
      at : float
        Attenuation (in dB)
      sym : bool
        Generates symmetric window if True.

    """
    if M < 1:
        return array([])
    if M == 1:
        return ones(1,'d')

    odd = M % 2
    if not sym and not odd:
        M = M+1

    # compute the parameter beta
    order = M - 1.0
    beta = cosh(1.0/order * arccosh(10**(abs(at)/20.)))
    k = r_[0:M]*1.0
    x = beta*cos(pi*k/M)
    #find the window's DFT coefficients
    # Use analytic definition of Chebyshev polynomial instead of expansion
    # from scipy.special. Using the expansion in scipy.special leads to errors.
    p = zeros(x.shape)
    p[x > 1] = cosh(order * arccosh(x[x > 1]))
    p[x < -1] = (1 - 2*(order%2)) * cosh(order * arccosh(-x[x < -1]))
    p[np.abs(x) <=1 ] = cos(order * arccos(x[np.abs(x) <= 1]))

    # Appropriate IDFT and filling up
    # depending on even/odd M
    if M % 2:
        w = real(fft(p))
        n = (M + 1) / 2
        w = w[:n] / w[0]
        w = concatenate((w[n - 1:0:-1], w))
    else:
        p = p * exp(1.j*pi / M * r_[0:M])
        w = real(fft(p))
        n = M / 2 + 1
        w = w / w[1]
        w = concatenate((w[n - 1:0:-1], w[1:n]))
    if not sym and not odd:
        w = w[:-1]
    return w

def chebwin(M, at, sym=1):
    """Dolph-Chebyshev window.

    INPUTS:

      M : int
        Window size
      at : float
        Attenuation (in dB)
      sym : bool
        Generates symmetric window if True.

    """
    if M < 1:
        return array([])
    if M == 1:
        return ones(1,'d')

    odd = M % 2
    if not sym and not odd:
        M = M+1

    # compute the parameter beta
    beta = cosh(1.0/(M-1.0)*arccosh(10**(at/20.)))
    k = r_[0:M]*1.0
    x = beta*cos(pi*k/M)
    #find the window's DFT coefficients
    p = zeros(x.shape) * 1.0
    for i in range(len(x)):
        if x[i] < 1:
            p[i] = cos((M - 1) * arccos(x[i]))
        else:
            p[i] = cosh((M - 1) * arccosh(x[i]))

    # Appropriate IDFT and filling up
    # depending on even/odd M
    if M % 2:
        w = real(fft(p));
        n = (M + 1) / 2;
        w = w[:n] / w[0];
        w = concatenate((w[n - 1:0:-1], w))
    else:
        p = p * exp(1.j*pi / M * r_[0:M])
        w = real(fft(p));
        n = M / 2 + 1;
        w = w / w[1];
        w = concatenate((w[n - 1:0:-1], w[1:n]));
    if not sym and not odd:
        w = w[:-1]
    return w

    def pinv(self, v):
        # first transform coordinates in inner elliptic coordinates
        # with (-a, 0) (a,0)
        v_in = dot(self.__rotmat.T, v - self.__centrum) / self.__a

        # m1 = cosh(mu), n1 = cos(nu)
        # solve for
        # v[0] = m1* n1
        # v[1]**2 = (m1**2 -1)(1 - n1**2)
        a = dot(v_in, v_in) + 1.
        m1 = sqrt(a /2. + sqrt(a**2 - 4 * v_in[0]**2) / 2.)

        # 1 <= m1, only rounding errors could say something else
        if abs(m1)  < 1.:
            mu = 0.
            m1 = 1.
        else:
            mu = arccosh(m1)

        n1 = v_in[0] / m1
        nu = acos(n1)

        # one has to consider that acos gives only values for the upper half
        if v_in[1] < 0. and abs(v_in[1]) > finfo(float).eps:
            nu = 2. * pi - nu

        return array([mu, nu])

def FS_Chebyshev(N, R, x_min, x_max, x_num, plot_far=False, dB_limit= -40):
    """Function to evaluate Fourier series coefficients of Chebyshev far-field
       pattern"""
       
    c = np.cosh(np.arccosh(R) / (N))
    c = str(c)     
           
    if(N % 2 == 0):
        m_start = int(-N / 2)
        m_stop = int(N / 2)
        N = str(N)       
        fun_str_re = 'special.eval_chebyt(' + N + ',' + c + '*np.cos(x/2))'
        m_index, zm = FS(fun_str_re, m_start=m_start, m_stop=m_stop, err_lim=1e-5)        
    else:
        m_start = -N # make this (2*P+1) ... and take fourier for only half period
        m_stop = N
        N = str(N)
        fun_str_re = 'special.eval_chebyt(' + N + ',' + c + '*np.cos(x))'
        m_index, zm = FS(fun_str_re, m_start=m_start, m_stop=m_stop, err_lim=1e-5)
        
    if(plot_far):
        x, AF = IFS(zm, 2 * np.pi, m_start, m_stop, x_min, x_max, x_num)        
        AF = 20 * np.log10(abs(AF))
        AF = pl.cutoff(AF, dB_limit)
        plt.plot(x * (180 / np.pi), AF); plt.axis('tight'); plt.grid(True)        
        plt.title('Far-field Pattern')
        plt.xlabel(r'$\phi$')
        plt.ylabel('AF')
        plt.show()        
                     
    return m_index, zm

def elliptic_hamiltonian(u,v, opts):
    
    '''
    
    Returns arrays of values for the first elliptic invariant (Hamiltonian) for a system with NL magnetic potential
    
    '''
    #beta = 0.6539
    #beta =  fixed['beta']
    #beta = 6.1246858201346921
    beta = opts.betae
    
    
    t = opts.t 
    
    c = opts.c
    #c = opts.c * np.sqrt(beta)
    
    f2u = u * np.sqrt(u**2 -1.) * np.arccosh(u)
    g2v = v * np.sqrt(1. - v**2) * (-np.pi/2 + np.arccos(v))
    
    elliptic = (f2u + g2v) / (u**2 - v**2)
    kfac = -1.*t*c*c
    
    return kfac*elliptic

    def computeInvariant(self, xHat, pxHat, yHat, pyHat):
        """Compute the 2nd invariant for the integrable elliptic potential"""

        xN = xHat / self.ellipticC
        yN = yHat / self.ellipticC

        # Elliptic coordinates
        u = (np.sqrt((xN + 1)**2 + yN**2) +\
                 np.sqrt((xN - 1)**2 + yN**2))/2.
        v = (np.sqrt((xN + 1)**2 + yN**2) -\
                 np.sqrt((xN - 1)**2 + yN**2))/2.

        # Elaborate x^2 + y^2
        f1u = self.ellipticC**2 * u**2 * (u**2 - 1)
        g1v = self.ellipticC**2 * v**2 * (1 - v**2)

        f2u = u * np.sqrt(u**2 - 1.) * np.arccosh(u)
        g2v = v * np.sqrt(1. - v**2) * (-np.pi/2 + np.arccos(v))

        kfac = self.ellipticT * self.ellipticC**2

        invariant = (x*py - y*px)**2 + self.ellipticC**2*(px**2) +\
            (2*self.ellipticC**2 ) * ((0.5*f1u - kfac*f2u)*v**2 \
            + (0.5 * g1v + kfac * g2v) * u**2 ) /(u**2 - v**2)

        return invariant

def ising_v(dtau, U, L=32, polar=0.5):
    """initialize the vector V of Ising fields

    .. math:: V = \\lambda (\\sigma_1, \\sigma_2, \\cdots, \\sigma_L)

    where the vector entries :math:`\\sigma_n=\\pm 1` are randomized subject
    to a threshold given by polar. And

    .. math:: \\cosh(\\lambda) = \\exp(\\Delta \\tau \\frac{U}{2})

    Parameters
    ----------
    dtau : float
        time spacing :math::`\\Delta\\Tau`
    U : float
        local Coulomb repulsion
    L : integer
        length of the array
    polar : float :math:`\\in (0, 1)`
        polarization threshold, probability of :math:`\\sigma_n=+ 1`

    Returns
    -------
    out : single dimension ndarray
    """
    lam = np.arccosh(np.exp(dtau*U/2))
    vis = np.ones(L)
    rand = np.random.rand(L)
    vis[rand > polar] = -1
    return vis*lam