Python cmath.phase() Examples

def angle_average(angles: np.ndarray) -> float:
    """
    Function to calculate the average value of a list of angles.

    Parameters
    ----------
    angles : numpy.ndarray
        Parallactic angles (deg).

    Returns
    -------
    float
        Average angle (deg).
    """

    cmath_rect = sum(cmath.rect(1, math.radians(ang)) for ang in angles)
    cmath_phase = cmath.phase(cmath_rect/len(angles))

    return math.degrees(cmath_phase) 

def isBlocked(self, id, allcells):
        for i in range(len(allcells)):
            if self.id != i and allcells[i].radius + 0.5 > allcells[self.id].radius:
                selfTheta = cmath.phase(complex(allcells[id].pos[0] - allcells[self.id].pos[0],
                                                allcells[self.id].pos[1] - allcells[id].pos[1])) % (2 * math.pi)
                deltaTheta = math.fabs(cmath.phase(complex(allcells[i].pos[0] - allcells[self.id].pos[0],
                                                           allcells[self.id].pos[1] - allcells[i].pos[1])) % (
                                               2 * math.pi) - selfTheta)
                if deltaTheta < math.pi / 2 and allcells[self.id].distance_from(allcells[i]) < allcells[
                    self.id].distance_from(allcells[id]):
                    k = (allcells[self.id].pos[1] - allcells[id].pos[1]) / (
                            allcells[id].pos[0] - allcells[self.id].pos[0])

                    d = math.fabs(k * (allcells[i].pos[0] - allcells[self.id].pos[0]) - allcells[self.id].pos[1] -
                                  allcells[i].pos[1]) / (1 + k * k) ** 0.5 - allcells[self.id].radius - allcells[
                            i].radius
                    if d < 0.5:
                        return True
        return False 

def calculate_E(self, id, allcells):
        theta = cmath.phase(
            complex(self.simplify_pos(id, allcells)[0], self.simplify_pos(id, allcells)[1])) % (
                        2 * math.pi)
        cos = math.cos(theta)
        sin = math.sin(theta)
        veloc_r = (allcells[id].veloc[0] - allcells[self.id].veloc[0]) * cos + (
                allcells[self.id].veloc[1] - allcells[id].veloc[1]) * sin
        d = allcells[self.id].distance_from(allcells[id]) - allcells[self.id].radius - allcells[
            id].radius
        if allcells[id].radius + 0.5 > allcells[self.id].radius:
            E = self.switchVeloc(veloc_r) * self.switchdistance(d)
        else:
            q = math.exp(-veloc_r) / (1.00001 - allcells[id].radius / allcells[self.id].radius)
            E = 1000 * q / d
        return [E * cos, E * sin, E, veloc_r] 

def arc_to(self, x1, y1):
        x0, y0 = self.get_current_point()
        dx, dy = x1 - x0, y1 - y0

        if abs(dx) < 1e-8:
            self.line_to(x1, y1)

        else:
            center = 0.5 * (x0 + x1) + 0.5 * (y0 + y1) * dy / dx
            theta0 = cmath.phase(x0 - center + y0*1j)
            theta1 = cmath.phase(x1 - center + y1*1j)
            r = abs(x0 - center + y0*1j)

            # we must ensure that the arc ends at (x1, y1)
            if x0 < x1:
                self.arc_negative(center, 0, r, theta0, theta1)
            else:
                self.arc(center, 0, r, theta0, theta1) 

def _dec_diag(self):
        """
        Call to create a circuit implementing the diagonal gate.
        """
        q = QuantumRegister(self.num_qubits)
        circuit = QuantumCircuit(q)
        # Since the diagonal is a unitary, all its entries have absolute value one and the diagonal
        # is fully specified by the phases of its entries
        diag_phases = [cmath.phase(z) for z in self.params]
        n = len(self.params)
        while n >= 2:
            angles_rz = []
            for i in range(0, n, 2):
                diag_phases[i // 2], rz_angle = _extract_rz(diag_phases[i], diag_phases[i + 1])
                angles_rz.append(rz_angle)
            num_act_qubits = int(np.log2(n))
            contr_qubits = q[self.num_qubits - num_act_qubits + 1:self.num_qubits]
            target_qubit = q[self.num_qubits - num_act_qubits]
            circuit.ucrz(angles_rz, contr_qubits, target_qubit)
            n //= 2
        return circuit


# extract a Rz rotation (angle given by first output) such that exp(j*phase)*Rz(z_angle)
# is equal to the diagonal matrix with entires exp(1j*ph1) and exp(1j*ph2) 

def phase2t(self, psi):
        """Given phase -pi < psi <= pi,
        returns the t value such that
        exp(1j*psi) = self.u1transform(self.point(t)).
        """
        def _deg(rads, domain_lower_limit):
            # Convert rads to degrees in [0, 360) domain
            degs = degrees(rads % (2*pi))

            # Convert to [domain_lower_limit, domain_lower_limit + 360) domain
            k = domain_lower_limit // 360
            degs += k * 360
            if degs < domain_lower_limit:
                degs += 360
            return degs

        if self.delta > 0:
            degs = _deg(psi, domain_lower_limit=self.theta)
        else:
            degs = _deg(psi, domain_lower_limit=self.theta)
        return (degs - self.theta)/self.delta 

def compute_g(derts__, Ave, fa=0):   # compute g from dx, dy

    for x0, derts_ in derts__:
        for derts in derts_:
            dy, dx = derts[-1][-1]

            if not fa:
                g = hypot(dy, dx)
            else:
                ga = hypot(phase(dy), phase(dx))
                if ga > pi: ga = two_pi - ga  # translate ga scope into (0, pi), unsigned
                g = int(ga * angle_coef)      # transform to fit in scope (-128, 127)

            derts[-1] = (g-Ave,) + derts[-1]  # return

    # ---------- compute_g() end ------------------------------------------------------------------------------------------- 

def standardize_phases(complex_array):
    """Uses `extract_phase` to obtain global phase from `array` and returns standardized array with global phase factor
    standardized.

    Parameters
    ----------
    complex_array: ndarray
        complex

    Returns
    -------
    ndarray (complex)
    """
    phase = extract_phase(complex_array)
    std_array = complex_array * np.exp(-1j * phase)
    return std_array 

def extract_phase(complex_array, position=None):
    """Extracts global phase from `complex_array` at given `position`. If position is not specified, the `position` is
    set to to an intermediate position to avoid machine-precision problems with tails of wavefunctions at beginning
    or end of the array.

    Parameters
    ----------
    complex_array: ndarray
        complex-valued array
    position: int, optional
        position where the phase is extracted (default value = None)
    """
    if position is None:
        flattened_position = np.argmax(
            np.abs(complex_array))  # extract phase from element with largest amplitude modulus
        position = np.unravel_index(flattened_position, complex_array.shape)
    return cmath.phase(complex_array[position]) 

def test_common(self):
        # random values
        M12 = 0.12716600+0.08765385j
        G12 = -0.34399429-0.1490931j
        aM12 = abs(M12)
        aG12 = abs(G12)
        phi12 = cmath.phase(-M12/G12)
        qp = flavio.physics.mesonmixing.common.q_over_p(M12, G12)
        DM = flavio.physics.mesonmixing.common.DeltaM(M12, G12)
        DG = flavio.physics.mesonmixing.common.DeltaGamma(M12, G12)
        # arXiv:0904.1869
        self.assertAlmostEqual(DM**2-1/4.*DG**2, 4*aM12**2-aG12**2, places=10) # (35)
        self.assertAlmostEqual(qp, -(DM+1j*DG/2)/(2*M12-1j*G12), places=10) # (37)
        self.assertAlmostEqual(qp, -(2*M12.conjugate()-1j*G12.conjugate())/(DM+1j*DG/2), places=10) # (37)
        self.assertAlmostEqual(DM*DG, -4*(M12*G12.conjugate()).real, places=10) # (36)
        self.assertAlmostEqual(DM*DG, 4*aM12*aG12*cos(phi12), places=10) # (39) 

def test_phase(self):
        self.assertAlmostEqual(phase(0), 0.)
        self.assertAlmostEqual(phase(1.), 0.)
        self.assertAlmostEqual(phase(-1.), pi)
        self.assertAlmostEqual(phase(-1.+1E-300j), pi)
        self.assertAlmostEqual(phase(-1.-1E-300j), -pi)
        self.assertAlmostEqual(phase(1j), pi/2)
        self.assertAlmostEqual(phase(-1j), -pi/2)

        # zeros
        self.assertEqual(phase(complex(0.0, 0.0)), 0.0)
        self.assertEqual(phase(complex(0.0, -0.0)), -0.0)
        self.assertEqual(phase(complex(-0.0, 0.0)), pi)
        self.assertEqual(phase(complex(-0.0, -0.0)), -pi)

        # infinities
        self.assertAlmostEqual(phase(complex(-INF, -0.0)), -pi)
        self.assertAlmostEqual(phase(complex(-INF, -2.3)), -pi)
        self.assertAlmostEqual(phase(complex(-INF, -INF)), -0.75*pi)
        self.assertAlmostEqual(phase(complex(-2.3, -INF)), -pi/2)
        self.assertAlmostEqual(phase(complex(-0.0, -INF)), -pi/2)
        self.assertAlmostEqual(phase(complex(0.0, -INF)), -pi/2)
        self.assertAlmostEqual(phase(complex(2.3, -INF)), -pi/2)
        self.assertAlmostEqual(phase(complex(INF, -INF)), -pi/4)
        self.assertEqual(phase(complex(INF, -2.3)), -0.0)
        self.assertEqual(phase(complex(INF, -0.0)), -0.0)
        self.assertEqual(phase(complex(INF, 0.0)), 0.0)
        self.assertEqual(phase(complex(INF, 2.3)), 0.0)
        self.assertAlmostEqual(phase(complex(INF, INF)), pi/4)
        self.assertAlmostEqual(phase(complex(2.3, INF)), pi/2)
        self.assertAlmostEqual(phase(complex(0.0, INF)), pi/2)
        self.assertAlmostEqual(phase(complex(-0.0, INF)), pi/2)
        self.assertAlmostEqual(phase(complex(-2.3, INF)), pi/2)
        self.assertAlmostEqual(phase(complex(-INF, INF)), 0.75*pi)
        self.assertAlmostEqual(phase(complex(-INF, 2.3)), pi)
        self.assertAlmostEqual(phase(complex(-INF, 0.0)), pi)

        # real or imaginary part NaN
        for z in complex_nans:
            self.assertTrue(math.isnan(phase(z))) 

def get_ckmangle_beta(par):
    r"""Returns the CKM angle $\beta$."""
    V = get_ckm(par)
    # see eq. (12.16) of http://pdg.lbl.gov/2015/reviews/rpp2014-rev-ckm-matrix.pdf
    return phase(-V[1,0]*V[1,2].conj()/V[2,0]/V[2,2].conj()) 

def phi12(wc_obj, par, meson):
    r"""$\phi_{12}=\text{arg}(-M_{12}/\Gamma_{12})"""
    M12, G12 = get_M12_G12(wc_obj, par, meson)
    return phase(M12 / G12) 

def a_fs(M12, G12):
    r"""Flavour-specific CP asymmetry in meson mixing as a function of
    $M_{12}$ and $\Gamma_{12}$."""
    aM12 = abs(M12)
    aG12 = abs(G12)
    phi12 = phase(-M12/G12)
    return aG12 / aM12 * sin(phi12) 

def test_bmixing(self):
        # just some trivial tests to see if calling the functions raises an error
        m12d = amplitude.M12(par, wc_B0, 'B0')
        m12s = amplitude.M12(par, wc_Bs, 'Bs')
        # check whether order of magnitudes of SM predictions are right
        ps = 1e-12*s
        self.assertAlmostEqual(observables.DeltaM_positive(wc_obj, par, 'B0')*ps, 0.55, places=0)
        self.assertAlmostEqual(observables.DeltaM_positive(wc_obj, par, 'Bs')*ps, 18, places=-1)
        self.assertAlmostEqual(observables.DeltaGamma_B(wc_obj, par, 'B0')/0.00261*ps, 1, places=0)
        self.assertAlmostEqual(observables.DeltaGamma_B(wc_obj, par, 'Bs')/0.088*ps, 1, places=0)
        self.assertAlmostEqual(observables.a_fs(wc_obj, par, 'B0')/-4.7e-4, 1, places=0)
        self.assertAlmostEqual(observables.a_fs(wc_obj, par, 'Bs')/2.22e-5, 1, places=0)
        self.assertAlmostEqual(observables.S_BJpsiK(wc_obj, par), 0.73, places=1)
        self.assertAlmostEqual(observables.S_Bspsiphi(wc_obj, par), asin(+0.038), places=2)
        # test classic formula: numerics of Wolfi's thesis
        w_par = par.copy()
        GF = w_par['GF']
        mW = w_par['m_W']
        mBs = w_par['m_Bs']
        fBs = 0.245
        w_par['f_Bs'] = fBs
        S0 = 2.31
        V = flavio.physics.ckm.ckm_wolfenstein(0.2254, 0.808, 0.177, 0.360)
        w_par['Vub'] = abs(V[0,2])
        w_par['Vcb'] = abs(V[1,2])
        w_par['Vus'] = abs(V[0,1])
        w_par['delta'] = -cmath.phase(V[0,2])
        etaB = 0.55
        BBsh = 1.22 # 0.952 * 1.517
        w_par['bag_Bs_1'] = BBsh/1.5173
        M12 = (GF**2 * mW**2/12/pi**2 * etaB * mBs * fBs**2 * BBsh
               * S0 * (V[2,2] * V[2,1].conj())**2)
        self.assertAlmostEqual(amplitude.M12(w_par, wc_Bs, 'Bs')/M12, 1, delta=0.01)
        self.assertAlmostEqual(observables.DeltaM_positive(wc_obj, w_par, 'Bs')*ps, 18.3, delta=0.2) 

def radial_average(L):
    # r = round(math.degrees(cmath.phase(sum(cmath.rect(1, math.radians(float(d))) for d in L)/len(L))),2)
    scipy = circmean(L, -180, 180)
    return scipy 

def predictTime(self, id, allcells):
        for t in range(10, 100):
            predictPos = [0, 0]
            predictPos[0] = allcells[id].pos[0] + allcells[id].veloc[0] * t * Consts["FRAME_DELTA"]
            predictPos[1] = allcells[id].pos[1] + allcells[id].veloc[1] * t * Consts["FRAME_DELTA"]
            mypos = [0, 0]
            mypos[0] = allcells[self.id].pos[0] - predictPos[0]
            mypos[1] = predictPos[1] - allcells[self.id].pos[1]
            myveloc = [0, 0]
            myveloc[0] = allcells[self.id].veloc[0]
            myveloc[1] = - allcells[self.id].veloc[1]
            theta = cmath.phase(complex(mypos[0], mypos[1])) % (2 * math.pi)
            cos = math.cos(theta)
            sin = math.sin(theta)
            veloc_r = myveloc[0] * cos + myveloc[1] * sin
            veloc_t = myveloc[1] * cos - myveloc[0] * sin
            v_max = 1.5 if self.time < 200 else 0.5
            d = (mypos[0] ** 2 + mypos[1] ** 2) ** 0.5 - allcells[id].radius - allcells[self.id].radius + 10
            a = Consts["EJECT_MASS_RATIO"] * Consts["DELTA_VELOC"] / Consts["FRAME_DELTA"]
            if d + veloc_r * math.fabs(veloc_t) / a + (veloc_r - v_max) * math.fabs(v_max - veloc_r) / (2 * a) < 0:
                t1 = math.fabs(veloc_t) / a
                delta = veloc_r ** 2 + 2 * a * (d + veloc_r * t1)
                # print(d, veloc_r * t1)
                if delta >= 0:
                    t2 = (delta ** 0.5 - veloc_r) / a
                    T = (t1 + t2) / Consts["FRAME_DELTA"]
                else:
                    continue
            else:
                t1 = math.fabs(veloc_t) / a
                t2 = math.fabs(v_max - veloc_r) / a
                t3 = (d + veloc_r * t1 + (veloc_r - v_max) * t2 / 2) / v_max
                T = (t1 + t2 + t3) / Consts["FRAME_DELTA"]
            leftarea = allcells[self.id].area() * (1 - Consts["EJECT_MASS_RATIO"]) ** (
                    (t1 + t2) / Consts["FRAME_DELTA"])
            if math.fabs(t - T) < 2 and leftarea > allcells[id].area() and t < self.deadTime(id,
                                                                                             allcells) and leftarea + \
                    allcells[id].area() > allcells[self.id].area():
                return t
        return None 

def mean_angle(self, deg):
        return phase(sum(rect(1, d) for d in deg) / len(deg)) 

def deconstruct_single_qubit_matrix_into_angles(
        mat: np.ndarray) -> Tuple[float, float, float]:
    """Breaks down a 2x2 unitary into more useful ZYZ angle parameters.

    Args:
        mat: The 2x2 unitary matrix to break down.

    Returns:
        A tuple containing the amount to phase around Z, then rotate around Y,
        then phase around Z (all in radians).
    """
    # Anti-cancel left-vs-right phase along top row.
    right_phase = cmath.phase(mat[0, 1] * np.conj(mat[0, 0])) + math.pi
    mat = np.dot(mat, _phase_matrix(-right_phase))

    # Cancel top-vs-bottom phase along left column.
    bottom_phase = cmath.phase(mat[1, 0] * np.conj(mat[0, 0]))
    mat = np.dot(_phase_matrix(-bottom_phase), mat)

    # Lined up for a rotation. Clear the off-diagonal cells with one.
    rotation = math.atan2(abs(mat[1, 0]), abs(mat[0, 0]))
    mat = np.dot(_rotation_matrix(-rotation), mat)

    # Cancel top-left-vs-bottom-right phase.
    diagonal_phase = cmath.phase(mat[1, 1] * np.conj(mat[0, 0]))

    # Note: Ignoring global phase.
    return right_phase + diagonal_phase, rotation * 2, bottom_phase 

def _sum_score_detail(reciprocal_space_vectors, solutions, granularity=None, amax=None):
    """Evaluates the probability that the trial value of (S0_vector | origin_offset) is correct,
    given the current estimate and the observations.  The trial value comes through the
    reciprocal space vectors, and the current estimate comes through the short list of
    DPS solutions. Actual return value is a sum of NH terms, one for each DPS solution, each ranging
    from -1.0 to 1.0"""

    nh = min(solutions.size(), 20)  # extended API
    sum_score = 0.0
    for t in range(nh):
        dfft = Directional_FFT(
            angle=Direction(solutions[t]),
            xyzdata=reciprocal_space_vectors,
            granularity=5.0,
            amax=amax,  # extended API XXX These values have to come from somewhere!
            F0_cutoff=11,
        )
        kval = dfft.kval()
        kmax = dfft.kmax()
        kval_cutoff = reciprocal_space_vectors.size() / 4.0
        if kval > kval_cutoff:
            ff = dfft.fft_result
            kbeam = ((-dfft.pmin) / dfft.delta_p) + 0.5
            Tkmax = cmath.phase(ff[kmax])
            backmax = math.cos(
                Tkmax + (2 * math.pi * kmax * kbeam / (2 * ff.size() - 1))
            )
            ### Here it should be possible to calculate a gradient.
            ### Then minimize with respect to two coordinates.  Use lbfgs?  Have second derivatives?
            ### can I do something local to model the cosine wave?
            ### direction of wave travel.  Period. phase.
            sum_score += backmax
    return sum_score 

def get_ckmangle_alpha(par):
    r"""Returns the CKM angle $\alpha$."""
    V = get_ckm(par)
    # see eq. (12.16) of http://pdg.lbl.gov/2015/reviews/rpp2014-rev-ckm-matrix.pdf
    return phase(-V[2,0]*V[2,2].conj()/V[0,0]/V[0,2].conj()) 

def get_ckmangle_gamma(par):
    r"""Returns the CKM angle $\gamma$."""
    V = get_ckm(par)
    # see eq. (12.16) of http://pdg.lbl.gov/2015/reviews/rpp2014-rev-ckm-matrix.pdf
    return phase(-V[0,0]*V[0,2].conj()/V[1,0]/V[1,2].conj())


# Some useful shorthands for CKM combinations appearing in FCNC amplitudes 

def arg(x):
    return cmath.phase(x) 

def test_phase(self):
        self.assertAlmostEqual(phase(0), 0.)
        self.assertAlmostEqual(phase(1.), 0.)
        self.assertAlmostEqual(phase(-1.), pi)
        self.assertAlmostEqual(phase(-1.+1E-300j), pi)
        self.assertAlmostEqual(phase(-1.-1E-300j), -pi)
        self.assertAlmostEqual(phase(1j), pi/2)
        self.assertAlmostEqual(phase(-1j), -pi/2)

        # zeros
        self.assertEqual(phase(complex(0.0, 0.0)), 0.0)
        self.assertEqual(phase(complex(0.0, -0.0)), -0.0)
        self.assertEqual(phase(complex(-0.0, 0.0)), pi)
        self.assertEqual(phase(complex(-0.0, -0.0)), -pi)

        # infinities
        self.assertAlmostEqual(phase(complex(-INF, -0.0)), -pi)
        self.assertAlmostEqual(phase(complex(-INF, -2.3)), -pi)
        self.assertAlmostEqual(phase(complex(-INF, -INF)), -0.75*pi)
        self.assertAlmostEqual(phase(complex(-2.3, -INF)), -pi/2)
        self.assertAlmostEqual(phase(complex(-0.0, -INF)), -pi/2)
        self.assertAlmostEqual(phase(complex(0.0, -INF)), -pi/2)
        self.assertAlmostEqual(phase(complex(2.3, -INF)), -pi/2)
        self.assertAlmostEqual(phase(complex(INF, -INF)), -pi/4)
        self.assertEqual(phase(complex(INF, -2.3)), -0.0)
        self.assertEqual(phase(complex(INF, -0.0)), -0.0)
        self.assertEqual(phase(complex(INF, 0.0)), 0.0)
        self.assertEqual(phase(complex(INF, 2.3)), 0.0)
        self.assertAlmostEqual(phase(complex(INF, INF)), pi/4)
        self.assertAlmostEqual(phase(complex(2.3, INF)), pi/2)
        self.assertAlmostEqual(phase(complex(0.0, INF)), pi/2)
        self.assertAlmostEqual(phase(complex(-0.0, INF)), pi/2)
        self.assertAlmostEqual(phase(complex(-2.3, INF)), pi/2)
        self.assertAlmostEqual(phase(complex(-INF, INF)), 0.75*pi)
        self.assertAlmostEqual(phase(complex(-INF, 2.3)), pi)
        self.assertAlmostEqual(phase(complex(-INF, 0.0)), pi)

        # real or imaginary part NaN
        for z in complex_nans:
            self.assertTrue(math.isnan(phase(z))) 

def test_phase(self):
        self.assertAlmostEqual(phase(0), 0.)
        self.assertAlmostEqual(phase(1.), 0.)
        self.assertAlmostEqual(phase(-1.), pi)
        self.assertAlmostEqual(phase(-1.+1E-300j), pi)
        self.assertAlmostEqual(phase(-1.-1E-300j), -pi)
        self.assertAlmostEqual(phase(1j), pi/2)
        self.assertAlmostEqual(phase(-1j), -pi/2)

        # zeros
        self.assertEqual(phase(complex(0.0, 0.0)), 0.0)
        self.assertEqual(phase(complex(0.0, -0.0)), -0.0)
        self.assertEqual(phase(complex(-0.0, 0.0)), pi)
        self.assertEqual(phase(complex(-0.0, -0.0)), -pi)

        # infinities
        self.assertAlmostEqual(phase(complex(-INF, -0.0)), -pi)
        self.assertAlmostEqual(phase(complex(-INF, -2.3)), -pi)
        self.assertAlmostEqual(phase(complex(-INF, -INF)), -0.75*pi)
        self.assertAlmostEqual(phase(complex(-2.3, -INF)), -pi/2)
        self.assertAlmostEqual(phase(complex(-0.0, -INF)), -pi/2)
        self.assertAlmostEqual(phase(complex(0.0, -INF)), -pi/2)
        self.assertAlmostEqual(phase(complex(2.3, -INF)), -pi/2)
        self.assertAlmostEqual(phase(complex(INF, -INF)), -pi/4)
        self.assertEqual(phase(complex(INF, -2.3)), -0.0)
        self.assertEqual(phase(complex(INF, -0.0)), -0.0)
        self.assertEqual(phase(complex(INF, 0.0)), 0.0)
        self.assertEqual(phase(complex(INF, 2.3)), 0.0)
        self.assertAlmostEqual(phase(complex(INF, INF)), pi/4)
        self.assertAlmostEqual(phase(complex(2.3, INF)), pi/2)
        self.assertAlmostEqual(phase(complex(0.0, INF)), pi/2)
        self.assertAlmostEqual(phase(complex(-0.0, INF)), pi/2)
        self.assertAlmostEqual(phase(complex(-2.3, INF)), pi/2)
        self.assertAlmostEqual(phase(complex(-INF, INF)), 0.75*pi)
        self.assertAlmostEqual(phase(complex(-INF, 2.3)), pi)
        self.assertAlmostEqual(phase(complex(-INF, 0.0)), pi)

        # real or imaginary part NaN
        for z in complex_nans:
            self.assertTrue(math.isnan(phase(z))) 

def arc(hrs, angle=60, mul=1.0):
    vs = rect(RADIUS * mul, PHI)  # Coords relative to arc origin
    ve = rect(RADIUS * mul, PHI)
    pe = PHI - angle * RV + TV
    rv = rect(1, -RV)  # Rotation vector for 1 minute (about OR)
    rot = rect(1, (3 - hrs) * pi / 6)  # hrs rotation (about [0,0])
    while phase(vs) > pe:
        ve *= rv
        # Translate to 0, 0
        yield ((vs + XLT) * rot, (ve + XLT) * rot)
        vs *= rv

# Currently unused. Draw the unit circle in a way which may readily be
# ported to other displays. 

def lateral_comp(P_):  # horizontal comparison between pixels at distance == rng

    derts__ = []

    for P in P_:
        x0 = P[1]
        derts_ = P[-1]
        new_derts_ = []     # new derts buffer

        _derts = derts_[0]
        gd, dx, dy, _ = _derts[-1]        # comp_angle always follows comp_gradient or hypot_g
        _a = complex(dx, dy)              # to complex number: _a = dx + dyj
        _a /= abs(_a)                     # normalize _a so that abs(_a) == 1 (hypot() from real and imaginary part of _a == 1)
        _a_radian = phase(_a)             # angular value of _a in radian: _a_radian in (-pi, pi)
        _dax, _ncomp = 0j, 0              # init ncomp, dx(complex) buffers

        for derts in derts_[1:]:
            # compute angle:
            g, dx, dy, _ = derts[-1]      # derts_ and new_derts_ are separate
            a = complex(dx, dy)           # to complex number: a = dx + dyj
            a /= abs(a)                   # normalize a so that abs(a) == 1 (hypot() from real and imaginary part of a == 1)
            a_radian = phase(a)           # angular value of a in radian: aa in (-pi, pi)

            da = rect(1, a_radian - _a_radian)   # convert bearing difference into complex form (rectangular coordinate)
            # complex d doesn't need to correct angle diff: if d > pi: d -= 255; elif d < -127: d += 255
            dx = da
            _dax += da      # bilateral accumulation
            _ncomp += 1     # bilateral accumulation

            new_derts_.append(_derts + [(_a, _a_radian), (0j, _dax, _ncomp)])  # return a and _dy = 0 + 0j in separate tuples
            _derts = derts                         # buffer derts
            _a, _aa, _dx, _ncomp = a, a_radian, dx, 1    # buffer last ncomp and dx

        new_derts_.append(_derts + [(_a, _a_radian), (0j, _dax, _ncomp)]) # return last derts

        derts__.append((x0, new_derts_))    # new line of P derts_ appended with new_derts_

    return derts__

    # ---------- lateral_comp() end ----------------------------------------------------------------------------------------- 

def ga_from_da(da_x, da_y, ncomp):
    " convert dx, dy to angular value then compute g"
    da_x = phase(da_x / ncomp)      # phase(da_x) is the same as phase(da_x / ncomp)
    da_y = phase(da_y / ncomp)      # phase(da_y) is the same as phase(da_y / ncomp)

    ga = hypot(da_x, da_y)
    if ga > pi: ga = two_pi - ga        # translate ga's scope into [0, pi) (g is unsigned)

    return int(ga * angle_coef) 

def compute_a(P_):  # compute angle from last dy, dx

    for P in P_:
        derts_ = P[-1]
        for derts in derts_:

            g, (dy, dx), _ = derts[-1]
            a = complex(dx, dy)  # to complex number: a = dx + dyj
            a /= abs(a)   # normalize a to make abs(a) == 1 (hypot() from real and imaginary part of a == 1)
            a_radian = phase(a)  # angular value of a in radians: a_radian in (-pi, pi)

            derts += [(a,), (a_radian,)]  # return

    # ---------- compute_a() end -------------------------------------------------------------------------------------------- 

def calc_angle(derts):  # compute a, return derts with angle

    g = derts[-1][0]
    dy, dx = derts[-1][-1]

    a = (dx + dy * 1j) / g
    a_radian = phase(a)

    derts[-1][1].insert(a, a_radian)  # a = dert[1], for i = dert[fia]

    return derts[-1]

    # ---------- calc_angle() end -----------------------------------------------------------------------------------------------