def _check_orthogonality(equations):
        """
        Helper method for _connect_to_cartesian. It checks if
        set of transformation equations create orthogonal curvilinear
        coordinate system

        Parameters
        ==========

        equations : Lambda
            Lambda of transformation equations

        """

        x1, x2, x3 = symbols("x1, x2, x3", cls=Dummy)
        equations = equations(x1, x2, x3)
        v1 = Matrix([diff(equations[0], x1),
                     diff(equations[1], x1), diff(equations[2], x1)])

        v2 = Matrix([diff(equations[0], x2),
                     diff(equations[1], x2), diff(equations[2], x2)])

        v3 = Matrix([diff(equations[0], x3),
                     diff(equations[1], x3), diff(equations[2], x3)])

        if any(simplify(i[0] + i[1] + i[2]) == 0 for i in (v1, v2, v3)):
            return False
        else:
            if simplify(v1.dot(v2)) == 0 and simplify(v2.dot(v3)) == 0 \
                and simplify(v3.dot(v1)) == 0:
                return True
            else:
                return False

def test_atan2():
    assert atan2(0, 0) == S.NaN
    assert atan2(0, 1) == 0
    assert atan2(1, 1) == pi/4
    assert atan2(1, 0) == pi/2
    assert atan2(1, -1) == 3*pi/4
    assert atan2(0, -1) == pi
    assert atan2(-1, -1) == -3*pi/4
    assert atan2(-1, 0) == -pi/2
    assert atan2(-1, 1) == -pi/4

    u = Symbol("u", positive=True)
    assert atan2(0, u) == 0
    u = Symbol("u", negative=True)
    assert atan2(0, u) == pi

    assert atan2(y, oo) ==  0
    assert atan2(y, -oo)==  2*pi*Heaviside(re(y)) - pi

    assert atan2(y, x).rewrite(log) == -I*log((x + I*y)/sqrt(x**2 + y**2))
    assert atan2(y, x).rewrite(atan) == 2*atan(y/(x + sqrt(x**2 + y**2)))

    assert diff(atan2(y, x), x) == -y/(x**2 + y**2)
    assert diff(atan2(y, x), y) == x/(x**2 + y**2)

    assert simplify(diff(atan2(y, x).rewrite(log), x)) == -y/(x**2 + y**2)
    assert simplify(diff(atan2(y, x).rewrite(log), y)) ==  x/(x**2 + y**2)

    assert isinstance(atan2(2, 3*I).n(), atan2)

def check_sign_algebraically(expr, param_names, species_id, initial_values):

    import sympy
    from sympy import Symbol

    # create variables for each param and species
    for parameter_name in param_names:
        exec("%s = Symbol('%s', positive=True)" % (parameter_name, parameter_name))

    # create variables for time and avogadro's constant: N_A
    exec("t = Symbol('t', positive=True)")
    exec("N_A = Symbol('N_A', positive=True)")

    exec("query_var = %s" % species_id)

    converted_rate_law = convert_rate_law(expr, output_type="sympy")
    if converted_rate_law == "piecewise":
        return "?"

    exec("rate = %s" % converted_rate_law)
    if sympy.simplify(rate.diff(query_var)).is_positive:
        return "monotonic_increasing"
    elif sympy.simplify(rate.diff(query_var)).is_negative:
        return "monotonic_decreasing"
    else:
        return "?"

def test_symbolic_twoport():
    circuit.default_toolkit = symbolic
    cir = SubCircuit()

    k = symbolic.kboltzmann
    var('R1 R0 C1 w T', real=True, positive=True)
    s = 1j*w

    cir['R0'] = R(1, gnd, r=R0)
    cir['R1'] = R(1, 2, r=R1)
#    cir['C1'] = C(2, gnd, c=C1)

    ## Add an AC source to verify that the source will not affect results
#    cir['IS'] = IS(1, gnd, iac=1) 

    ## Run symbolic 2-port analysis
    twoport_ana = TwoPortAnalysis(cir, Node('1'), gnd, Node('2'), gnd,
                                  noise = True, toolkit=symbolic,
                                  noise_outquantity = 'v')
    result = twoport_ana.solve(freqs=s, complexfreq=True)
    
    ABCD = Matrix(result['twoport'].A)
    ABCD.simplify()

    assert_array_equal(ABCD, np.array([[1 + 0*R1*C1*s, R1],
                                    [(1 + 0*R0*C1*s + 0*R1*C1*s) / R0,  (R0 + R1)/R0]]))

    assert_array_equal(simplify(result['Sin'] - (4*k*T/R0 + 4*R1*k*T/R0**2)), 0)
    assert_array_equal(simplify(result['Svn']), 4*k*T*R1)

def test_gosper_sum_AeqB_part3():
    f3a = 1/n**4
    f3b = (6*n + 3)/(4*n**4 + 8*n**3 + 8*n**2 + 4*n + 3)
    f3c = 2**n*(n**2 - 2*n - 1)/(n**2*(n + 1)**2)
    f3d = n**2*4**n/((n + 1)*(n + 2))
    f3e = 2**n/(n + 1)
    f3f = 4*(n - 1)*(n**2 - 2*n - 1)/(n**2*(n + 1)**2*(n - 2)**2*(n - 3)**2)
    f3g = (n**4 - 14*n**2 - 24*n - 9)*2**n/(n**2*(n + 1)**2*(n + 2)**2*
           (n + 3)**2)

    # g3a -> no closed form
    g3b = m*(m + 2)/(2*m**2 + 4*m + 3)
    g3c = 2**m/m**2 - 2
    g3d = S(2)/3 + 4**(m + 1)*(m - 1)/(m + 2)/3
    # g3e -> no closed form
    g3f = -(-S(1)/16 + 1/((m - 2)**2*(m + 1)**2))  # the AeqB key is wrong
    g3g = -S(2)/9 + 2**(m + 1)/((m + 1)**2*(m + 3)**2)

    g = gosper_sum(f3a, (n, 1, m))
    assert g is None
    g = gosper_sum(f3b, (n, 1, m))
    assert g is not None and simplify(g - g3b) == 0
    g = gosper_sum(f3c, (n, 1, m - 1))
    assert g is not None and simplify(g - g3c) == 0
    g = gosper_sum(f3d, (n, 1, m))
    assert g is not None and simplify(g - g3d) == 0
    g = gosper_sum(f3e, (n, 0, m - 1))
    assert g is None
    g = gosper_sum(f3f, (n, 4, m))
    assert g is not None and simplify(g - g3f) == 0
    g = gosper_sum(f3g, (n, 1, m))
    assert g is not None and simplify(g - g3g) == 0

def reciprocal_frame_test():
    Print_Function()

    metric = "1 # #," + "# 1 #," + "# # 1,"

    (e1, e2, e3) = MV.setup("e1 e2 e3", metric)

    print("g_{ij} =\n", MV.metric)

    E = e1 ^ e2 ^ e3
    Esq = (E * E).scalar()
    print("E =", E)
    print("E**2 =", Esq)
    Esq_inv = 1 / Esq

    E1 = (e2 ^ e3) * E
    E2 = (-1) * (e1 ^ e3) * E
    E3 = (e1 ^ e2) * E

    print("E1 = (e2^e3)*E =", E1)
    print("E2 =-(e1^e3)*E =", E2)
    print("E3 = (e1^e2)*E =", E3)

    w = E1 | e2
    w = w.expand()
    print("E1|e2 =", w)

    w = E1 | e3
    w = w.expand()
    print("E1|e3 =", w)

    w = E2 | e1
    w = w.expand()
    print("E2|e1 =", w)

    w = E2 | e3
    w = w.expand()
    print("E2|e3 =", w)

    w = E3 | e1
    w = w.expand()
    print("E3|e1 =", w)

    w = E3 | e2
    w = w.expand()
    print("E3|e2 =", w)

    w = E1 | e1
    w = (w.expand()).scalar()
    Esq = expand(Esq)
    print("(E1|e1)/E**2 =", simplify(w / Esq))

    w = E2 | e2
    w = (w.expand()).scalar()
    print("(E2|e2)/E**2 =", simplify(w / Esq))

    w = E3 | e3
    w = (w.expand()).scalar()
    print("(E3|e3)/E**2 =", simplify(w / Esq))
    return

def hamiltonian(Lagrangian, t = Symbol('t'), delta = False):
    """
    Returns the Hamiltonian of the Lagrangian.

    Examples
    ========

    >>> from sympy import *
    >>> t, k = symbols('t k')
    >>> x = symbols('x', cls=Function)
    >>> hamiltonian(diff(x(t),t)**2/2 - k*x(t)**2/2)
    k*x**2/2 + v_x**2/2
    """
    Lagrangian = simplify(Lagrangian)
    var_list = [list(x.atoms(Function))[0] for x in Lagrangian.atoms(Derivative)]
    nvar = len(var_list)
    # New variables.
    str_list = [ str(variable).replace("("+str(t)+")","") for variable in var_list ]
    v_subs = {diff(var_list[i],t): Symbol('v_' + str_list[i]) for i in range(nvar)}
    x_subs = {var_list[i]: Symbol(str_list[i]) for i in range(nvar)}
    # Hamiltonian calculus.
    dxdLv = 0
    for variable in var_list:
        dxdLv += diff(variable,t)*diff(Lagrangian, diff(variable,t))
    result = simplify((dxdLv - Lagrangian).subs(v_subs).subs(x_subs))
    if delta:
        v0_subs = {Symbol('v_' + str_list[i]): Symbol('v_' + str_list[i] + "0") for i in range(nvar)}
        x0_subs = {Symbol(str_list[i]): Symbol(str_list[i] + "0") for i in range(nvar)}
        return result - result.subs(v0_subs).subs(x0_subs)
    else:
        return result

def test_medium():
    m1 = Medium('m1')
    assert m1.intrinsic_impedance == sqrt(u0/e0)
    assert m1.speed == 1/sqrt(e0*u0)
    assert m1.refractive_index == c*sqrt(e0*u0)
    assert m1.permittivity == e0
    assert m1.permeability == u0
    m2 = Medium('m2', epsilon, mu)
    assert m2.intrinsic_impedance == sqrt(mu/epsilon)
    assert m2.speed == 1/sqrt(epsilon*mu)
    assert m2.refractive_index == c*sqrt(epsilon*mu)
    assert m2.permittivity == epsilon
    assert m2.permeability == mu
    # Increasing electric permittivity and magnetic permeability
    # by small amount from its value in vacuum.
    m3 = Medium('m3', 9.0*10**(-12)*s**4*A**2/(m**3*kg), 1.45*10**(-6)*kg*m/(A**2*s**2))
    assert m3.refractive_index > m1.refractive_index
    assert m3 > m1
    # Decreasing electric permittivity and magnetic permeability
    # by small amount from its value in vacuum.
    m4 = Medium('m4', 7.0*10**(-12)*s**4*A**2/(m**3*kg), 1.15*10**(-6)*kg*m/(A**2*s**2))
    assert m4.refractive_index < m1.refractive_index
    assert m4 < m1
    m5 = Medium('m5', permittivity=710*10**(-12)*s**4*A**2/(m**3*kg), n=1.33)
    assert simplify(m5.intrinsic_impedance - 6.24845417765552*kg*m**2/(A**2*s**3)) == 0
    assert abs(m5.speed - 225407863.157895*m/s) < 1e-6*m/s
    assert simplify(m5.refractive_index - 1.33000000000000) == 0
    assert simplify(m5.permittivity - 7.1e-10*A**2*s**4/(kg*m**3)) == 0
    assert simplify(m5.permeability - 2.77206575232851e-8*kg*m/(A**2*s**2)) == 0

    def get_S(self, theta=0., strain_type="Engineering"):
        """
        Calculates the compliance matrix for the ply oriented at angle theta
        in the specified strain units
        Parameters
        ----------
        theta : 'Int' or 'Float'
            Angle in degrees of ply orientation.
        strain_type : 'String'
            Specifies 'Engineering' or 'Tensorial' strain.
        """

        if strain_type.lower().startswith('e'):
            strain_multiplier = 1
        else:
            strain_multiplier = 2

        compliance = sp.Matrix([[1 / self.E11, -self.v21 / self.E22, 0],
                               [-self.v12 / self.E11, 1 / self.E22, 0],
                               [0, 0, 1 / (strain_multiplier * self.G12)]])

        if theta == 0.:
            return compliance
        else:
            T = T_Matrix(theta)
            if strain_type.lower().startswith('e'):
                R = R_Matrix()
                TI = sp.simplify(R * T.inv() * R.inv())
            else:
                TI = T.inv()
            return sp.simplify(sp.N(TI * compliance * T, chop=1e-10))

def test_SVCVS_laplace_d3_n1():
    """Test VCCS with a laplace defined transfer function with second order
    numerator and third order denominator
    """

    pycircuit.circuit.circuit.default_toolkit = symbolic
    cir = SubCircuit()

    n1,n2 = cir.add_nodes('1','2')

    b0,a0,a1,a2,a3,Gdc = [sympy.Symbol(symname, real=True) for
                                symname in 'b0,a0,a1,a2,a3,Gdc'
                                .split(',')]

    s = sympy.Symbol('s', complex=True)

    cir['VS']   = VS( n1, gnd, vac=1)
    cir['VCVS'] = SVCVS( n1, gnd, n2, gnd,
                        denominator = [a0, a1, a2, a3],
                        numerator   = [b0, 0, 0])

    res = AC(cir, toolkit=symbolic).solve(s, complexfreq=True)

    assert_equal(sympy.simplify(res.v(n2,gnd)),
                 sympy.simplify((b0*s*s)/(a0*s*s*s+a1*s*s+a2*s+a3)))

def test_nullor_vva():
    """Test nullor element by building a V-V amplifier"""
    pycircuit.circuit.circuit.default_toolkit = symbolic

    c = SubCircuit()

    Vin = Symbol('Vin')
    R1 =Symbol('R1')
    R2 = Symbol('R2')
    
    nin = c.add_node('in')
    n1 = c.add_node('n1')
    nout = c.add_node('out')
     
    c['vin'] = VS(nin, gnd, vac=Vin)
    c['R1'] = R(n1, gnd, r=R1)
    c['R2'] = R(nout, n1, r=R2)
    c['nullor'] = Nullor(n1, nin, gnd, nout)
    
    result = AC(c, toolkit=symbolic).solve(Symbol('s'))
    
    vout = result.v(nout)

    assert simplify(vout - Vin * (R1 + R2) / R1) == 0, \
        'Did not get the expected result, %s != 0'% \
        str(simplify(vout - Vin * (R1 + R2) / R1))

    def calculate_eigenvalues(self):
        r"""Calculate the two eigenvalues :math:`\lambda_i(x)` of the potential :math:`V(x)`.
        We can do this by symbolic calculations. The multiplicities are taken into account.
        Note: This function is idempotent and the eigenvalues are memoized for later reuse.
        """
        if self._eigenvalues_s is not None:
            return

        # Symbolic formula for the eigenvalues of a general 2x2 matrix
        T = self._potential_s.trace()
        D = self._potential_s.det()

        l1 = (T + sympy.sqrt(T**2 - 4*D)) * sympy.Rational(1,2)
        l2 = (T - sympy.sqrt(T**2 - 4*D)) * sympy.Rational(1,2)

        # Symbolic simplification may fail
        if self._try_simplify:
            try:
                l1 = sympy.simplify(l1)
                l2 = sympy.simplify(l2)
            except:
                pass

        # The symbolic expressions for the eigenvalues
        self._eigenvalues_s = (l1, l2)

        # The numerical functions for the eigenvalues
        self._eigenvalues_n = tuple([ sympy.lambdify(self._variables, item, "numpy") for item in self._eigenvalues_s ])

def interpolation(f, psi, points):
    """
    Given a function f(x), return the approximation to
    f(x) in the space V, spanned by psi, that interpolates
    f at the given points. Must have len(points) = len(psi)
    """
    N = len(psi) - 1
    A = sym.zeros((N+1, N+1))
    b = sym.zeros((N+1, 1))
    # Wrap psi and f in Python functions rather than expressions
    # so that we can evaluate psi at points[i] (alternative to subs?)
    x = sym.Symbol('x')
    psi = [sym.lambdify([x], psi[i]) for i in range(N+1)]
    f = sym.lambdify([x], f)
    print '...evaluating matrix...'
    for i in range(N+1):
        for j in range(N+1):
            print '(%d,%d)' % (i, j)
            A[i,j] = psi[j](points[i])
        b[i,0] = f(points[i])
    print
    print 'A:\n', A, '\nb:\n', b
    c = A.LUsolve(b)
    # c is a sympy Matrix object, turn to list
    c = [sym.simplify(c[i,0]) for i in range(c.shape[0])]
    print 'coeff:', c
    u = 0
    for i in range(len(psi)):
        u += c[i]*psi[i](x)
    # Alternative:
    # u = sum(c[i,0]*psi[i] for i in range(len(psi)))
    print 'approximation:', sym.simplify(u)
    return u, c

    def __pow__(self, e, z=None):
        """
        Perform multinomial expansion and return terms of order less than
        `self.order`.
        """
        data = {}
        exps,coeffs = zip(*(self.d.items()))
        m = len(exps)

        print "EXPONENT:",e

        multinoms = sympy.multinomial.multinomial_coefficients_iterator(m,int(e))
        for k,mcoeff in multinoms:
            exp = sum( expi*ki for expi,ki in zip(exps,k) )
            
            if exp < self.order:
                try:
                    coeff = sympy.simplify(
                        mcoeff * sympy.prod(ci**ki for ci,ki in zip(coeffs,k))
                        ) 
                    coeff = sympy.simplify(coeff + data[exp])                
                except KeyError:
                    data[exp] = coeff
            
        return PuiseuxSeries(data, self.alpha, self.order, self.var)

def test_equation():
    p1 = Point(0, 0)
    p2 = Point(1, 1)
    l1 = Line(p1, p2)
    l3 = Line(Point(x1, x1), Point(x1, 1 + x1))

    assert simplify(l1.equation()) in (x - y, y - x)
    assert simplify(l3.equation()) in (x - x1, x1 - x)
    assert simplify(l1.equation()) in (x - y, y - x)
    assert simplify(l3.equation()) in (x - x1, x1 - x)

    assert Line(p1, Point(1, 0)).equation(x=x, y=y) == y
    assert Line(p1, Point(0, 1)).equation() == x
    assert Line(Point(2, 0), Point(2, 1)).equation() == x - 2
    assert Line(p2, Point(2, 1)).equation() == y - 1

    assert Line3D(Point(x1, x1, x1), Point(y1, y1, y1)
        ).equation() == (-x + y, -x + z)
    assert Line3D(Point(1, 2, 3), Point(2, 3, 4)
        ).equation() == (-x + y - 1, -x + z - 2)
    assert Line3D(Point(1, 2, 3), Point(1, 3, 4)
        ).equation() == (x - 1, -y + z - 1)
    assert Line3D(Point(1, 2, 3), Point(2, 2, 4)
        ).equation() == (y - 2, -x + z - 2)
    assert Line3D(Point(1, 2, 3), Point(2, 3, 3)
        ).equation() == (-x + y - 1, z - 3)
    assert Line3D(Point(1, 2, 3), Point(1, 2, 4)
        ).equation() == (x - 1, y - 2)
    assert Line3D(Point(1, 2, 3), Point(1, 3, 3)
        ).equation() == (x - 1, z - 3)
    assert Line3D(Point(1, 2, 3), Point(2, 2, 3)
        ).equation() == (y - 2, z - 3)

def test_gosper_sum_AeqB_part2():
    f2a = n**2*a**n
    f2b = (n - r/2)*binomial(r, n)
    f2c = factorial(n - 1)**2/(factorial(n - x)*factorial(n + x))

    g2a = -a*(a + 1)/(a - 1)**3 + a**(
        m + 1)*(a**2*m**2 - 2*a*m**2 + m**2 - 2*a*m + 2*m + a + 1)/(a - 1)**3
    g2b = (m - r)*binomial(r, m)/2
    ff = factorial(1 - x)*factorial(1 + x)
    g2c = 1/ff*(
        1 - 1/x**2) + factorial(m)**2/(x**2*factorial(m - x)*factorial(m + x))

    g = gosper_sum(f2a, (n, 0, m))
    assert g is not None and simplify(g - g2a) == 0
    g = gosper_sum(f2b, (n, 0, m))
    assert g is not None and simplify(g - g2b) == 0
    g = gosper_sum(f2c, (n, 1, m))
    assert g is not None and simplify(g - g2c) == 0

    # delete these lines and unXFAIL the nan test below when it passes
    f2d = n*(n + a + b)*a**n*b**n/(factorial(n + a)*factorial(n + b))
    g2d = 1/(factorial(a - 1)*factorial(
        b - 1)) - a**(m + 1)*b**(m + 1)/(factorial(a + m)*factorial(b + m))
    assert simplify(
        sum(f2d.subs(n, i) for i in range(3)) - g2d.subs(m, 2)) == 0

    def _check_orthogonality(self):
        """
        Helper method for _connect_to_cartesian. It checks if
        set of transformation equations create orthogonal curvilinear
        coordinate system

        Parameters
        ==========

        equations : tuple
            Tuple of transformation equations

        """

        eq = self._transformation_equations()

        v1 = Matrix([diff(eq[0], self.x), diff(eq[1], self.x), diff(eq[2], self.x)])
        v2 = Matrix([diff(eq[0], self.y), diff(eq[1], self.y), diff(eq[2], self.y)])
        v3 = Matrix([diff(eq[0], self.z), diff(eq[1], self.z), diff(eq[2], self.z)])

        if any(simplify(i[0] + i[1] + i[2]) == 0 for i in (v1, v2, v3)):
            return False
        else:
            if simplify(v1.dot(v2)) == 0 and simplify(v2.dot(v3)) == 0 and simplify(v3.dot(v1)) == 0:
                return True
            else:
                return False

def get_solid_harmonics(shell, xyz):
    x, y, z = xyz
    theta = Symbol("theta", real=True)
    phi = Symbol("phi", real=True)
    r = Symbol("r", real=True)
    result = []
    for m in range(shell+1):
        if m > 0:
            part_plus = (Ylm(shell,m,theta,phi)*r**shell).expand()
            part_plus = part_plus.subs(exp(I*phi),(x+I*y)/sin(theta)/r)
            part_min = (Ylm(shell,-m,theta,phi)*r**shell).expand()
            part_min = part_min.subs(exp(-I*phi),(x-I*y)/sin(theta)/r)
            sym = simplify(((-1)**m*part_plus+part_min)/sqrt(2))
            sym = sym.subs(r*cos(theta),z)
            sym = sym.subs(r**2,x**2+y**2+z**2)
            sym = sym.subs(cos(theta)**2,1-sin(theta)**2)
            sym = simplify(sym)
            asym = simplify(((-1)**m*part_plus-part_min)/I/sqrt(2))
            asym = asym.subs(r*cos(theta),z)
            asym = asym.subs(r**2,x**2+y**2+z**2)
            asym = asym.subs(cos(theta)**2,1-sin(theta)**2)
            asym = simplify(asym)
            result.append(sym)
            result.append(asym)
        else:
            part = (Ylm(shell,0,theta,phi)*r**shell).expand()
            part = part.subs(r*cos(theta),z)
            part = part.subs(r**2,x**2+y**2+z**2)
            part = simplify(part)
            result.append(part)
    return result

def test_M23():
    x = symbols('x', complex=True)

    assert solve(x - 1/sqrt(1 + x**2)) == [
        simplify(-I*sqrt((sqrt(5) + 1)/2)),
        simplify(   sqrt((sqrt(5) - 1)/2)),
    ]

def calculate(str_expr, return_float=False):
    """
    Evaluate a string STR_EXPR.
    if STR_EXPR is a mathematical expression, simplify it with sympy.simplify
    Otherwise, search WolframAlpha for answers
    """
    try:
        result = simplify(str_expr)
        return result
    except Exception:
        # use Wolfram Alpha
        
        payload = {'input': str_expr, 'appid': 'UAGAWR-3X6Y8W777Q'}
        r = requests.get('http://api.wolframalpha.com/v2/query', params=payload)

        soup = BeautifulSoup(r.content,"lxml")
        if soup.queryresult['success'] == 'true':
            pods = soup.queryresult.findAll('pod')
            assert len(pods) >= 2
            num_string =  pods[1].plaintext.string.split()[0]
            if return_float:
                return float(simplify(num_string.replace(u"\u00D7", '*')))
            else:
                return num_string
        else:
            return "I don't know the answer"