def _solve_type_2(symo, X, Y, Z, th):
    """Solution for the equation:
    X*S + Y*C = Z
    """
    symo.write_line("# X*sin({0}) + Y*cos({0}) = Z".format(th))
    X = symo.replace(trigsimp(X), 'X', th)
    Y = symo.replace(trigsimp(Y), 'Y', th)
    Z = symo.replace(trigsimp(Z), 'Z', th)
    YPS = var('YPS'+str(th))
    if X == tools.ZERO and Y != tools.ZERO:
        C = symo.replace(Z/Y, 'C', th)
        symo.add_to_dict(YPS, (tools.ONE, - tools.ONE))
        symo.add_to_dict(th, atan2(YPS*sqrt(1-C**2), C))
    elif X != tools.ZERO and Y == tools.ZERO:
        S = symo.replace(Z/X, 'S', th)
        symo.add_to_dict(YPS, (tools.ONE, - tools.ONE))
        symo.add_to_dict(th, atan2(S, YPS*sqrt(1-S**2)))
    elif Z == tools.ZERO:
        symo.add_to_dict(YPS, (tools.ONE, tools.ZERO))
        symo.add_to_dict(th, atan2(-Y, X) + YPS*pi)
    else:
        B = symo.replace(X**2 + Y**2, 'B', th)
        D = symo.replace(B - Z**2, 'D', th)
        symo.add_to_dict(YPS, (tools.ONE, - tools.ONE))
        S = symo.replace((X*Z + YPS * Y * sqrt(D))/B, 'S', th)
        C = symo.replace((Y*Z - YPS * X * sqrt(D))/B, 'C', th)
        symo.add_to_dict(th, atan2(S, C))

    def _set_inv_trans_equations(curv_coord_name):
        """
        Store information about inverse transformation equations for
        pre-defined coordinate systems.

        Parameters
        ==========

        curv_coord_name : str
            Name of coordinate system

        """
        if curv_coord_name == 'cartesian':
            return lambda x, y, z: (x, y, z)

        if curv_coord_name == 'spherical':
            return lambda x, y, z: (
                sqrt(x**2 + y**2 + z**2),
                acos(z/sqrt(x**2 + y**2 + z**2)),
                atan2(y, x)
            )
        if curv_coord_name == 'cylindrical':
            return lambda x, y, z: (
                sqrt(x**2 + y**2),
                atan2(y, x),
                z
            )
        raise ValueError('Wrong set of parameters.'
                         'Type of coordinate system is defined')

    def computeSymbolicModel(self):
        """
        Symbollically computes G(X,t) and stores the models and lambda functions.
        :return:
        """
        # satellite position in ECI
        x = self._stateSymb[0]
        y = self._stateSymb[1]
        z = self._stateSymb[2]
        x_dot = self._stateSymb[3]
        y_dot = self._stateSymb[4]
        z_dot = self._stateSymb[5]

       # x, y, z = sp.symbols('x y z')
        r_xy = sp.sqrt(x**2+y**2)
        #x_dot, y_dot, z_dot = sp.symbols('x_dot y_dot z_dot')

        right_ascension = sp.atan2(y,x)

        declination = sp.atan2(z,r_xy)

        g1 = sp.lambdify((x, y, z, x_dot, y_dot, z_dot), right_ascension, "numpy")
        g2 = sp.lambdify((x, y, z, x_dot, y_dot, z_dot), declination, "numpy")

        self._modelSymb = [right_ascension, declination]
        self._modelLambda = [g1, g2]

        return self._modelSymb

def _solve_type_2(symo, X, Y, Z, th):
    """Solution for the equation:
    X*S + Y*C = Z
    """
    symo.write_line("# X*sin({0}) + Y*cos({0}) = Z".format(th))
    X = symo.replace(symo.CS12_simp(X), "X", th)
    Y = symo.replace(symo.CS12_simp(Y), "Y", th)
    Z = symo.replace(symo.CS12_simp(Z), "Z", th)
    YPS = var("YPS" + str(th))
    if X == ZERO and Y != ZERO:
        C = symo.replace(Z / Y, "C", th)
        symo.add_to_dict(YPS, (ONE, -ONE))
        symo.add_to_dict(th, atan2(YPS * sqrt(1 - C ** 2), C))
    elif X != ZERO and Y == ZERO:
        S = symo.replace(Z / X, "S", th)
        symo.add_to_dict(YPS, (ONE, -ONE))
        symo.add_to_dict(th, atan2(S, YPS * sqrt(1 - S ** 2)))
    elif Z == ZERO:
        symo.add_to_dict(YPS, (ONE, ZERO))
        symo.add_to_dict(th, atan2(-Y, X) + YPS * pi)
    else:
        B = symo.replace(X ** 2 + Y ** 2, "B", th)
        D = symo.replace(B - Z ** 2, "D", th)
        symo.add_to_dict(YPS, (ONE, -ONE))
        S = symo.replace((X * Z + YPS * Y * sqrt(D)) / B, "S", th)
        C = symo.replace((Y * Z - YPS * X * sqrt(D)) / B, "C", th)
        symo.add_to_dict(th, atan2(S, C))

def test_im():
    x, y = symbols('x,y')
    a, b = symbols('a,b', real=True)

    r = Symbol('r', real=True)
    i = Symbol('i', imaginary=True)

    assert im(nan) == nan

    assert im(oo*I) == oo
    assert im(-oo*I) == -oo

    assert im(0) == 0

    assert im(1) == 0
    assert im(-1) == 0

    assert im(E*I) == E
    assert im(-E*I) == -E

    assert im(x) == im(x)
    assert im(x*I) == re(x)
    assert im(r*I) == r
    assert im(r) == 0
    assert im(i*I) == 0
    assert im(i) == -I * i

    assert im(x + y) == im(x + y)
    assert im(x + r) == im(x)
    assert im(x + r*I) == im(x) + r

    assert im(im(x)*I) == im(x)

    assert im(2 + I) == 1
    assert im(x + I) == im(x) + 1

    assert im(x + y*I) == im(x) + re(y)
    assert im(x + r*I) == im(x) + r

    assert im(log(2*I)) == pi/2

    assert im((2 + I)**2).expand(complex=True) == 4

    assert im(conjugate(x)) == -im(x)
    assert conjugate(im(x)) == im(x)

    assert im(x).as_real_imag() == (im(x), 0)

    assert im(i*r*x).diff(r) == im(i*x)
    assert im(i*r*x).diff(i) == -I * re(r*x)

    assert im(
        sqrt(a + b*I)) == (a**2 + b**2)**Rational(1, 4)*sin(atan2(b, a)/2)
    assert im(a * (2 + b*I)) == a*b

    assert im((1 + sqrt(a + b*I))/2) == \
        (a**2 + b**2)**Rational(1, 4)*sin(atan2(b, a)/2)/2

    assert im(x).rewrite(re) == x - re(x)
    assert (x + im(y)).rewrite(im, re) == x + y - re(y)

def get_border_coord(xy, other_xy, wh, n_type):
    (x, y) = xy
    (other_x, other_y) = other_xy
    (w, h) = wh
    if n_type == TYPE_REACTION:
        edge_angle = degrees(atan2(other_y - y, other_x - x)) if other_y != y or other_x != x else 0
        diag_angle = degrees(atan2(h, w))
        abs_edge_angle = abs(edge_angle)
        if diag_angle < abs_edge_angle < 180 - diag_angle:
            y += h if edge_angle > 0 else -h
        else:
            x += w if abs_edge_angle <= 90 else -w
        return x, y
    elif n_type == TYPE_COMPARTMENT:
        c_bottom_x, c_bottom_y, c_top_x, c_top_y = x - w, y - h, x + w, y + h
        inside_y = c_bottom_y <= other_y <= c_top_y
        inside_x = c_bottom_x <= other_x <= c_top_x

        if inside_x:
            return other_x, (c_bottom_y if abs(other_y - c_bottom_y) < abs(other_y - c_top_y) else c_top_y)
        elif inside_y:
            return (c_bottom_x if abs(other_x - c_bottom_x) < abs(other_x - c_top_x) else c_top_x), other_y
        else:
            return max(c_bottom_x, min(other_x, c_top_x)), max(c_bottom_y, min(other_y, c_top_y))
    else:
        diag = pow(pow(x - other_x, 2) + pow(y - other_y, 2), 0.5)
        transformation = lambda z, other_z: (w * (((other_z - z) / diag) if diag else 1)) + z
        return transformation(x, other_x), transformation(y, other_y)

def test_functional_differential_geometry_ch2():
    # From "Functional Differential Geometry" as of 2011
    # by Sussman and Wisdom.
    x0, y0, r0, theta0 = symbols('x0, y0, r0, theta0', real=True)
    x, y, r, theta = symbols('x, y, r, theta', real=True)
    f = Function('f')

    assert (R2_p.point_to_coords(R2_r.point([x0, y0])) ==
                Matrix([sqrt(x0**2+y0**2), atan2(y0, x0)]))
    assert (R2_r.point_to_coords(R2_p.point([r0, theta0])) ==
                Matrix([r0*cos(theta0), r0*sin(theta0)]))
    #TODO jacobian page 12 - 32

    field = ScalarField(R2_r, [x, y], f(x, y))
    p1_in_rect = R2_r.point([x0, y0])
    p1_in_polar = R2_p.point([sqrt(x0**2 + y0**2), atan2(y0,x0)])
    assert field(p1_in_rect) == f(x0, y0)
    # TODO better simplification for the next one
    #print simplify(field(p1_in_polar))
    #assert simplify(field(p1_in_polar)) == f(x0, y0)

    p_r = R2_r.point([x0, y0])
    p_p = R2_p.point([r0, theta0])
    assert R2.x(p_r) == x0
    assert R2.x(p_p) == r0*cos(theta0)
    assert R2.r(p_p) == r0
    assert R2.r(p_r) == sqrt(x0**2 + y0**2)
    assert R2.theta(p_r) == atan2(y0, x0)

    h = R2.x*R2.r**2 + R2.y**3
    assert h(p_r) == x0*(x0**2 + y0**2) + y0**3
    assert h(p_p) == r0**3*sin(theta0)**3 + r0**3*cos(theta0)

def toeuler(quat):
    """Convert quaternion rotation to roll-pitch-yaw Euler angles."""
    q0, q1, q2, q3 = quat
    roll = sympy.atan2(2*(q2*q3 + q0*q1), q0**2 - q1**2 - q2**2 + q3**2)
    pitch = -sympy.asin(2*(q1*q3 - q0*q2))
    yaw = sympy.atan2(2*(q1*q2 + q0*q3), q0**2 + q1**2 - q2**2 - q3**2)
    return np.array([roll, pitch, yaw])

def test_as_real_imag():
    n = pi**1000
    # the special code for working out the real
    # and complex parts of a power with Integer exponent
    # should not run if there is no imaginary part, hence
    # this should not hang
    assert n.as_real_imag() == (n, 0)

    # issue 6261
    x = Symbol('x')
    assert sqrt(x).as_real_imag() == \
        ((re(x)**2 + im(x)**2)**(S(1)/4)*cos(atan2(im(x), re(x))/2),
     (re(x)**2 + im(x)**2)**(S(1)/4)*sin(atan2(im(x), re(x))/2))

    # issue 3853
    a, b = symbols('a,b', real=True)
    assert ((1 + sqrt(a + b*I))/2).as_real_imag() == \
           (
               (a**2 + b**2)**Rational(
                   1, 4)*cos(atan2(b, a)/2)/2 + Rational(1, 2),
               (a**2 + b**2)**Rational(1, 4)*sin(atan2(b, a)/2)/2)

    assert sqrt(a**2).as_real_imag() == (sqrt(a**2), 0)
    i = symbols('i', imaginary=True)
    assert sqrt(i**2).as_real_imag() == (0, abs(i))

    def _set_inv_trans_equations(self, curv_coord_name):
        """
        Store information about some default, pre-defined inverse
        transformation equations.

        Parameters
        ==========

        curv_coord_name : str
            The type of the new coordinate system.

        """

        equations_mapping = {
            'cartesian': (self.x, self.y, self.z),
            'spherical': (sqrt(self.x**2 + self.y**2 + self.z**2),
                          acos((self.z) / sqrt(self.x**2 + self.y**2 + self.z**2)),
                          atan2(self.y, self.x)),
            'cylindrical': (sqrt(self.x**2 + self.y**2),
                            atan2(self.y, self.x),
                            self.z)
        }
        if curv_coord_name not in equations_mapping:
            raise ValueError('Wrong set of parameters.'
                             'Type of coordinate system is defined')
        self._inv_transformation_eqs = equations_mapping[curv_coord_name]

def test_twave():
    A1, phi1, A2, phi2, f = symbols('A1, phi1, A2, phi2, f')
    n = Symbol('n')  # Refractive index
    t = Symbol('t')  # Time
    x = Symbol('x')  # Spatial varaible
    k = Symbol('k')  # Wave number
    E = Function('E')
    w1 = TWave(A1, f, phi1)
    w2 = TWave(A2, f, phi2)
    assert w1.amplitude == A1
    assert w1.frequency == f
    assert w1.phase == phi1
    assert w1.wavelength == c/(f*n)
    assert w1.time_period == 1/f
    w3 = w1 + w2
    assert w3.amplitude == sqrt(A1**2 + 2*A1*A2*cos(phi1 - phi2) + A2**2)
    assert w3.frequency == f
    assert w3.wavelength == c/(f*n)
    assert w3.time_period == 1/f
    assert w3.angular_velocity == 2*pi*f
    assert w3.wavenumber == 2*pi*f*n/c
    assert simplify(w3.rewrite('sin') - sqrt(A1**2 + 2*A1*A2*cos(phi1 - phi2)
    + A2**2)*sin(pi*f*n*x*s/(149896229*m) - 2*pi*f*t + atan2(A1*cos(phi1)
    + A2*cos(phi2), A1*sin(phi1) + A2*sin(phi2)) + pi/2)) == 0
    assert w3.rewrite('pde') == epsilon*mu*Derivative(E(x, t), t, t) + Derivative(E(x, t), x, x)
    assert w3.rewrite(cos) == sqrt(A1**2 + 2*A1*A2*cos(phi1 - phi2)
    + A2**2)*cos(pi*f*n*x*s/(149896229*m) - 2*pi*f*t + atan2(A1*cos(phi1)
    + A2*cos(phi2), A1*sin(phi1) + A2*sin(phi2)))
    assert w3.rewrite('exp') == sqrt(A1**2 + 2*A1*A2*cos(phi1 - phi2)
    + A2**2)*exp(I*(pi*f*n*x*s/(149896229*m) - 2*pi*f*t
    + atan2(A1*cos(phi1) + A2*cos(phi2), A1*sin(phi1) + A2*sin(phi2))))

def test_re():
    x, y = symbols('x,y')
    a, b = symbols('a,b', real=True)

    r = Symbol('r', real=True)
    i = Symbol('i', imaginary=True)

    assert re(nan) == nan

    assert re(oo) == oo
    assert re(-oo) == -oo

    assert re(0) == 0

    assert re(1) == 1
    assert re(-1) == -1

    assert re(E) == E
    assert re(-E) == -E

    assert re(x) == re(x)
    assert re(x*I) == -im(x)
    assert re(r*I) == 0
    assert re(r) == r
    assert re(i*I) == I * i
    assert re(i) == 0

    assert re(x + y) == re(x + y)
    assert re(x + r) == re(x) + r

    assert re(re(x)) == re(x)

    assert re(2 + I) == 2
    assert re(x + I) == re(x)

    assert re(x + y*I) == re(x) - im(y)
    assert re(x + r*I) == re(x)

    assert re(log(2*I)) == log(2)

    assert re((2 + I)**2).expand(complex=True) == 3

    assert re(conjugate(x)) == re(x)
    assert conjugate(re(x)) == re(x)

    assert re(x).as_real_imag() == (re(x), 0)

    assert re(i*r*x).diff(r) == re(i*x)
    assert re(i*r*x).diff(i) == I*r*im(x)

    assert re(
        sqrt(a + b*I)) == (a**2 + b**2)**Rational(1, 4)*cos(atan2(b, a)/2)
    assert re(a * (2 + b*I)) == 2*a

    assert re((1 + sqrt(a + b*I))/2) == \
        (a**2 + b**2)**Rational(1, 4)*cos(atan2(b, a)/2)/2 + Rational(1, 2)

    assert re(x).rewrite(im) == x - im(x)
    assert (x + re(y)).rewrite(re, im) == x + y - im(y)

 def points(self):
     from sympy import atan2
     res = []
     for p in self.boundary_points:
         x,y = p.point.x, p.point.y
         rotation = atan2(p.next_inner.y - y, p.next_inner.x - x)
         angle = atan2(p.prev_inner.y - y, p.prev_inner.x - x) - rotation
         while angle < 0:
             angle += 2*pi
         if angle >= pi:
             p.prev_inner, p.next_inner = p.next_inner, p.prev_inner

def test_as_real_imag():
    n = pi**1000
    # the special code for working out the real
    # and complex parts of a power with Integer exponent
    # should not run if there is no imaginary part, hence
    # this should not hang
    assert n.as_real_imag() == (n, 0)

    # issue 3162
    x = Symbol('x')
    assert sqrt(x).as_real_imag() == \
    ((re(x)**2 + im(x)**2)**(S(1)/4)*cos(atan2(im(x), re(x))/2), \
     (re(x)**2 + im(x)**2)**(S(1)/4)*sin(atan2(im(x), re(x))/2))

    def assertPairAngle(self, dyn):
        x1 = dyn.obja["tr.x"]
        y1 = dyn.obja["tr.y"]
        x2 = dyn.objb["tr.x"]
        y2 = dyn.objb["tr.y"]
        a1 = dyn.obja["rt.angle"]
        a2 = dyn.objb["rt.angle"]
        thetaa = dyn.thetaa
        name = dyn.name

        if not isTimeConstant(thetaa, self.symbols):
            return True

        # Convert the attachments to polar, add the body angle
        # and then reconvert to rectangular

        att1 = dyn.getAttachment(dyn.obja, "p")
        if att1[0] != 0:
            att1 = (att1[0], a1 + att1[1], att1[2])
            i1, j1, m1 = Globals.convertAttachment(att1, "r")
        else:
            i1 = 0
            j1 = 0

        att2 = dyn.getAttachment(dyn.objb, "p")
        if att2[0] != 0:
            att2 = (att2[0], a2 + att2[1], att2[2])
            i2, j2, m2 = Globals.convertAttachment(att2, "r")
        else:
            i2 = 0
            j2 = 0

        dx = sympy.sympify((x2 + i2) - (x1 + i1))
        dy = sympy.sympify((y2 + j2) - (y1 + j1))

        if dx == 0 and sympy.simplify(sympy.Mod(thetaa, sympy.pi)) != 0:
            self.printer.print_diagnostic(
                2,
                "thetaa assertion failed for dynamic %s, set=%s, inferred=%s."
                % (name, str(thetaa), str(sympy.sympify(sympy.atan2(dx, dy)))),
            )
            return False
        elif dy == 0 and sympy.simplify(sympy.Mod(thetaa + sympy.pi / 2, sympy.pi)) != 0:
            self.printer.print_diagnostic(
                2,
                "thetaa assertion failed for dynamic %s, set=%s, inferred=%s."
                % (name, str(thetaa), str(sympy.sympify(sympy.atan2(dx, dy)))),
            )
            return False

        return True

def sym_matrix_to_euler(mat):
	'''(symbolic!) Converts a rotation matrix (YXZ order) to 3 euler angles XYZ. Note that
	for simplicity this ignores the singularities.'''
	#XYZ
	#x = sympy.atan2(-mat[1, 2], mat[2, 2])
	#y = sympy.asin(mat[0, 2])
	#z = sympy.atan2(-mat[0, 1], mat[0, 0])

	#YXZ
	x = sympy.asin(-mat[1, 2])
	y = sympy.atan2(mat[0, 2], mat[2, 2])
	z = sympy.atan2(mat[1, 0], mat[1, 1])

	return [x, y, z]

def test_transformation_equations():
    from sympy import symbols
    x, y, z = symbols('x y z')
    a = CoordSys3D('a')
    # Str
    a._connect_to_standard_cartesian('spherical')
    assert a._transformation_equations() == (a.x * sin(a.y) * cos(a.z),
                          a.x * sin(a.y) * sin(a.z),
                          a.x * cos(a.y))
    assert a.lame_coefficients() == (1, a.x, a.x * sin(a.y))
    assert a._inverse_transformation_equations() == (sqrt(a.x**2 + a.y**2 +a.z**2),
                          acos((a.z) / sqrt(a.x**2 + a.y**2 + a.z**2)),
                          atan2(a.y, a.x))
    a._connect_to_standard_cartesian('cylindrical')
    assert a._transformation_equations() == (a.x * cos(a.y), a.x * sin(a.y), a.z)
    assert a.lame_coefficients() == (1, a.y, 1)
    assert a._inverse_transformation_equations() == (sqrt(a.x**2 + a.y**2),
                            atan2(a.y, a.x), a.z)
    a._connect_to_standard_cartesian('cartesian')
    assert a._transformation_equations() == (a.x, a.y, a.z)
    assert a.lame_coefficients() == (1, 1, 1)
    assert a._inverse_transformation_equations() == (a.x, a.y, a.z)
    # Variables and expressions
    a._connect_to_standard_cartesian(((x, y, z), (x, y, z)))
    assert a._transformation_equations() == (a.x, a.y, a.z)
    assert a.lame_coefficients() == (1, 1, 1)
    assert a._inverse_transformation_equations() == (a.x, a.y, a.z)
    a._connect_to_standard_cartesian(((x, y, z), ((x * cos(y), x * sin(y), z))), inverse=False)
    assert a._transformation_equations() == (a.x * cos(a.y), a.x * sin(a.y), a.z)
    assert simplify(a.lame_coefficients()) == (1, sqrt(a.x**2), 1)
    a._connect_to_standard_cartesian(((x, y, z), (x * sin(y) * cos(z), x * sin(y) * sin(z), x * cos(y))), inverse=False)
    assert a._transformation_equations() == (a.x * sin(a.y) * cos(a.z),
                          a.x * sin(a.y) * sin(a.z),
                          a.x * cos(a.y))
    assert simplify(a.lame_coefficients()) == (1, sqrt(a.x**2), sqrt(sin(a.y)**2*a.x**2))
    # Equations
    a._connect_to_standard_cartesian((a.x*sin(a.y)*cos(a.z), a.x*sin(a.y)*sin(a.z), a.x*cos(a.y)), inverse=False)
    assert a._transformation_equations() == (a.x * sin(a.y) * cos(a.z),
                          a.x * sin(a.y) * sin(a.z),
                          a.x * cos(a.y))
    assert simplify(a.lame_coefficients()) == (1, sqrt(a.x**2), sqrt(sin(a.y)**2*a.x**2))
    a._connect_to_standard_cartesian((a.x, a.y, a.z))
    assert a._transformation_equations() == (a.x, a.y, a.z)
    assert simplify(a.lame_coefficients()) == (1, 1, 1)
    assert a._inverse_transformation_equations() == (a.x, a.y, a.z)
    a._connect_to_standard_cartesian((a.x * cos(a.y), a.x * sin(a.y), a.z), inverse=False)
    assert a._transformation_equations() == (a.x * cos(a.y), a.x * sin(a.y), a.z)
    assert simplify(a.lame_coefficients()) == (1, sqrt(a.x**2), 1)

    raises(ValueError, lambda: a._connect_to_standard_cartesian((x, y, z)))

def test_atan2_expansion():
    assert cancel(atan2(x + 1, x ** 2).diff(x) - atan((x + 1) / x ** 2).diff(x)) == 0
    assert cancel(atan(x / y).series(x, 0, 5) - atan2(x, y).series(x, 0, 5) + atan2(0, y) - atan(0)) == O(x ** 5)
    assert cancel(atan(x / y).series(y, 1, 4) - atan2(x, y).series(y, 1, 4) + atan2(x, 1) - atan(x)) == O(y ** 4)
    assert cancel(
        atan((x + y) / y).series(y, 1, 3) - atan2(x + y, y).series(y, 1, 3) + atan2(1 + x, 1) - atan(1 + x)
    ) == O(y ** 3)
    assert Matrix([atan2(x, y)]).jacobian([x, y]) == Matrix([[y / (x ** 2 + y ** 2), -x / (x ** 2 + y ** 2)]])

def test_atan2_expansion():
    assert cancel(atan2(x ** 2, x + 1).diff(x) - atan(x ** 2 / (x + 1)).diff(x)) == 0
    assert cancel(atan(y / x).series(y, 0, 5) - atan2(y, x).series(y, 0, 5) + atan2(0, x) - atan(0)) == O(y ** 5)
    assert cancel(atan(y / x).series(x, 1, 4) - atan2(y, x).series(x, 1, 4) + atan2(y, 1) - atan(y)) == O(x ** 4)
    assert cancel(
        atan((y + x) / x).series(x, 1, 3) - atan2(y + x, x).series(x, 1, 3) + atan2(1 + y, 1) - atan(1 + y)
    ) == O(x ** 3)
    assert Matrix([atan2(y, x)]).jacobian([y, x]) == Matrix([[x / (y ** 2 + x ** 2), -y / (y ** 2 + x ** 2)]])

    def inferPairAngle(self, obja, atta, objb, attb):
        x1 = obja["tr.x"]
        y1 = obja["tr.y"]
        x2 = objb["tr.x"]
        y2 = objb["tr.y"]
        a1 = obja["rt.angle"]
        a2 = objb["rt.angle"]

        # Convert the attachments to polar, add the body angle
        # and then reconvert to rectangular

        att1 = Globals.convertAttachment(atta, "p")
        if att1[0] != 0:
            att1 = (att1[0], a1 + att1[1], att1[2])
            i1, j1, m1 = Globals.convertAttachment(att1, "r")
        else:
            i1 = 0
            j1 = 0

        att2 = Globals.convertAttachment(attb, "p")
        if att2[0] != 0:
            att2 = (att2[0], a2 + att2[1], att2[2])
            i2, j2, m2 = Globals.convertAttachment(att2, "r")
        else:
            i2 = 0
            j2 = 0

        dx = sympy.sympify((x2 + i2) - (x1 + i1))
        dy = sympy.sympify((y2 + j2) - (y1 + j1))

        return sympy.sympify(sympy.atan2(dx, dy))