Python numpy.arcsin() Examples

def arcsine(x, null=(-np.inf, np.inf)):
    '''
    arcsine(x) is equivalent to asin(x) except that it also works on sparse arrays.

    The optional argument null (default, (-numpy.inf, numpy.inf)) may be specified to indicate what
    value(s) should be assigned when x < -1 or x > 1. If only one number is given, then it is used
    for both values; otherwise the first value corresponds to <-1 and the second to >1.  If null is
    None, then an error is raised when invalid values are encountered.
    '''
    if sps.issparse(x):
        x = x.copy()
        x.data = arcsine(x.data, null=null, rtol=rtol, atol=atol)
        return x
    else: x = np.asarray(x)
    try:    (nln,nlp) = null
    except Exception: (nln,nlp) = (null,null)
    ii = None if nln is None else np.where(x < -1)
    jj = None if nlp is None else np.where(x > 1)
    if ii: x[ii] = 0
    if jj: x[jj] = 0
    x = np.arcsin(x)
    if ii: x[ii] = nln
    if jj: x[jj] = nlp
    return x 

def inverse_method(self,N,d):
        
        t = np.linspace(1e-3,0.999,N)
        f = np.log( t / (1 - t) )
        f = f/f[0]
        
        psi= np.pi*f
        cosPsi = np.cos(psi)
        sinTheta = ( np.abs(cosPsi) + (1-np.abs(cosPsi))*np.random.rand(len(cosPsi)))
        
        theta = np.arcsin(sinTheta)
        theta = np.pi-theta + (2*theta - np.pi)*np.round(np.random.rand(len(t)))
        cosPhi = cosPsi/sinTheta
        phi = np.arccos(cosPhi)*(-1)**np.round(np.random.rand(len(t)))
        
        coords = SkyCoord(phi*u.rad,(np.pi/2-theta)*u.rad,d*np.ones(len(phi))*u.pc)

        return coords 

def mindmag(self, s):
        """Calculates the minimum value of dMag for projected separation
        
        Args:
            s (float):
                Projected separations (AU)
        
        Returns:
            mindmag (float):
                Minimum planet delta magnitude
        """
        if s == 0.0:
            mindmag = self.cdmin1
        elif s < self.rmin*np.sin(self.bstar):
            mindmag = self.cdmin1-2.5*np.log10(self.Phi(np.arcsin(s/self.rmin)))
        elif s < self.rmax*np.sin(self.bstar):
            mindmag = self.cdmin2+5.0*np.log10(s)
        elif s <= self.rmax:
            mindmag = self.cdmin3-2.5*np.log10(self.Phi(np.arcsin(s/self.rmax)))
        else:
            mindmag = np.inf
        
        return mindmag 

def maxdmag(self, s):
        """Calculates the maximum value of dMag for projected separation
        
        Args:
            s (float):
                Projected separation (AU)
        
        Returns:
            maxdmag (float):
                Maximum planet delta magnitude
        
        """
        
        if s == 0.0:
            maxdmag = self.cdmax - 2.5*np.log10(self.Phi(np.pi))
        elif s < self.rmax:
            maxdmag = self.cdmax - 2.5*np.log10(np.abs(self.Phi(np.pi-np.arcsin(s/self.rmax))))
        else:
            maxdmag = self.cdmax - 2.5*np.log10(self.Phi(np.pi/2.0))

        return maxdmag 

def xyz2uvN(xyz, planeID=1):
    ID1 = (int(planeID) - 1 + 0) % 3
    ID2 = (int(planeID) - 1 + 1) % 3
    ID3 = (int(planeID) - 1 + 2) % 3
    normXY = np.sqrt(xyz[:, [ID1]] ** 2 + xyz[:, [ID2]] ** 2)
    normXY[normXY < 0.000001] = 0.000001
    normXYZ = np.sqrt(xyz[:, [ID1]] ** 2 + xyz[:, [ID2]] ** 2 + xyz[:, [ID3]] ** 2)
    v = np.arcsin(xyz[:, [ID3]] / normXYZ)
    u = np.arcsin(xyz[:, [ID1]] / normXY)
    valid = (xyz[:, [ID2]] < 0) & (u >= 0)
    u[valid] = np.pi - u[valid]
    valid = (xyz[:, [ID2]] < 0) & (u <= 0)
    u[valid] = -np.pi - u[valid]
    uv = np.hstack([u, v])
    uv[np.isnan(uv[:, 0]), 0] = 0
    return uv 

def xyz_from_rotm(R):
  R = R.reshape(-1,3,3)
  xyz = np.empty((R.shape[0],3), dtype=R.dtype)
  for bidx in range(R.shape[0]):
    if R[bidx,0,2] < 1:
      if R[bidx,0,2] > -1:
        xyz[bidx,1] = np.arcsin(R[bidx,0,2])
        xyz[bidx,0] = np.arctan2(-R[bidx,1,2], R[bidx,2,2])
        xyz[bidx,2] = np.arctan2(-R[bidx,0,1], R[bidx,0,0])
      else:
        xyz[bidx,1] = -np.pi/2
        xyz[bidx,0] = -np.arctan2(R[bidx,1,0],R[bidx,1,1])
        xyz[bidx,2] = 0
    else:
      xyz[bidx,1] = np.pi/2
      xyz[bidx,0] = np.arctan2(R[bidx,1,0], R[bidx,1,1])
      xyz[bidx,2] = 0
  return xyz.squeeze() 

def zyx_from_rotm(R):
  R = R.reshape(-1,3,3)
  zyx = np.empty((R.shape[0],3), dtype=R.dtype)
  for bidx in range(R.shape[0]):
    if R[bidx,2,0] < 1:
      if R[bidx,2,0] > -1:
        zyx[bidx,1] = np.arcsin(-R[bidx,2,0])
        zyx[bidx,0] = np.arctan2(R[bidx,1,0], R[bidx,0,0])
        zyx[bidx,2] = np.arctan2(R[bidx,2,1], R[bidx,2,2])
      else:
        zyx[bidx,1] = np.pi / 2
        zyx[bidx,0] = -np.arctan2(-R[bidx,1,2], R[bidx,1,1])
        zyx[bidx,2] = 0
    else:
      zyx[bidx,1] = -np.pi / 2
      zyx[bidx,0] = np.arctan2(-R[bidx,1,2], R[bidx,1,1])
      zyx[bidx,2] = 0
  return zyx.squeeze() 

def test_branch_cuts(self):
        # check branch cuts and continuity on them
        _check_branch_cut(np.log,   -0.5, 1j, 1, -1, True)
        _check_branch_cut(np.log2,  -0.5, 1j, 1, -1, True)
        _check_branch_cut(np.log10, -0.5, 1j, 1, -1, True)
        _check_branch_cut(np.log1p, -1.5, 1j, 1, -1, True)
        _check_branch_cut(np.sqrt,  -0.5, 1j, 1, -1, True)

        _check_branch_cut(np.arcsin, [ -2, 2],   [1j, 1j], 1, -1, True)
        _check_branch_cut(np.arccos, [ -2, 2],   [1j, 1j], 1, -1, True)
        _check_branch_cut(np.arctan, [0-2j, 2j],  [1,  1], -1, 1, True)

        _check_branch_cut(np.arcsinh, [0-2j,  2j], [1,   1], -1, 1, True)
        _check_branch_cut(np.arccosh, [ -1, 0.5], [1j,  1j], 1, -1, True)
        _check_branch_cut(np.arctanh, [ -2,   2], [1j, 1j], 1, -1, True)

        # check against bogus branch cuts: assert continuity between quadrants
        _check_branch_cut(np.arcsin, [0-2j, 2j], [ 1,  1], 1, 1)
        _check_branch_cut(np.arccos, [0-2j, 2j], [ 1,  1], 1, 1)
        _check_branch_cut(np.arctan, [ -2,  2], [1j, 1j], 1, 1)

        _check_branch_cut(np.arcsinh, [ -2,  2, 0], [1j, 1j, 1], 1, 1)
        _check_branch_cut(np.arccosh, [0-2j, 2j, 2], [1,  1,  1j], 1, 1)
        _check_branch_cut(np.arctanh, [0-2j, 2j, 0], [1,  1,  1j], 1, 1) 

def test_branch_cuts_complex64(self):
        # check branch cuts and continuity on them
        _check_branch_cut(np.log,   -0.5, 1j, 1, -1, True, np.complex64)
        _check_branch_cut(np.log2,  -0.5, 1j, 1, -1, True, np.complex64)
        _check_branch_cut(np.log10, -0.5, 1j, 1, -1, True, np.complex64)
        _check_branch_cut(np.log1p, -1.5, 1j, 1, -1, True, np.complex64)
        _check_branch_cut(np.sqrt,  -0.5, 1j, 1, -1, True, np.complex64)

        _check_branch_cut(np.arcsin, [ -2, 2],   [1j, 1j], 1, -1, True, np.complex64)
        _check_branch_cut(np.arccos, [ -2, 2],   [1j, 1j], 1, -1, True, np.complex64)
        _check_branch_cut(np.arctan, [0-2j, 2j],  [1,  1], -1, 1, True, np.complex64)

        _check_branch_cut(np.arcsinh, [0-2j,  2j], [1,   1], -1, 1, True, np.complex64)
        _check_branch_cut(np.arccosh, [ -1, 0.5], [1j,  1j], 1, -1, True, np.complex64)
        _check_branch_cut(np.arctanh, [ -2,   2], [1j, 1j], 1, -1, True, np.complex64)

        # check against bogus branch cuts: assert continuity between quadrants
        _check_branch_cut(np.arcsin, [0-2j, 2j], [ 1,  1], 1, 1, False, np.complex64)
        _check_branch_cut(np.arccos, [0-2j, 2j], [ 1,  1], 1, 1, False, np.complex64)
        _check_branch_cut(np.arctan, [ -2,  2], [1j, 1j], 1, 1, False, np.complex64)

        _check_branch_cut(np.arcsinh, [ -2,  2, 0], [1j, 1j, 1], 1, 1, False, np.complex64)
        _check_branch_cut(np.arccosh, [0-2j, 2j, 2], [1,  1,  1j], 1, 1, False, np.complex64)
        _check_branch_cut(np.arctanh, [0-2j, 2j, 0], [1,  1,  1j], 1, 1, False, np.complex64) 

def test_against_cmath(self):
        import cmath

        points = [-1-1j, -1+1j, +1-1j, +1+1j]
        name_map = {'arcsin': 'asin', 'arccos': 'acos', 'arctan': 'atan',
                    'arcsinh': 'asinh', 'arccosh': 'acosh', 'arctanh': 'atanh'}
        atol = 4*np.finfo(complex).eps
        for func in self.funcs:
            fname = func.__name__.split('.')[-1]
            cname = name_map.get(fname, fname)
            try:
                cfunc = getattr(cmath, cname)
            except AttributeError:
                continue
            for p in points:
                a = complex(func(np.complex_(p)))
                b = cfunc(p)
                assert_(abs(a - b) < atol, "%s %s: %s; cmath: %s" % (fname, p, a, b)) 

def HaversineDistance(lon1,lat1,lon2,lat2):
    """
    Function to calculate the great circle distance between two points
    using the Haversine formula
    """
    R = 6371. #Mean radius of the Earth

    # convert decimal degrees to radians
    lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])

    # haversine formula
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = np.sin(dlat/2.)**2. + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.)**2.
    c = 2.*np.arcsin(np.sqrt(a))
    distance = R * c

    return distance 

def test_branch_cuts_complex64(self):
        # check branch cuts and continuity on them
        yield _check_branch_cut, np.log,   -0.5, 1j, 1, -1, True, np.complex64
        yield _check_branch_cut, np.log2,  -0.5, 1j, 1, -1, True, np.complex64
        yield _check_branch_cut, np.log10, -0.5, 1j, 1, -1, True, np.complex64
        yield _check_branch_cut, np.log1p, -1.5, 1j, 1, -1, True, np.complex64
        yield _check_branch_cut, np.sqrt,  -0.5, 1j, 1, -1, True, np.complex64

        yield _check_branch_cut, np.arcsin, [ -2, 2],   [1j, 1j], 1, -1, True, np.complex64
        yield _check_branch_cut, np.arccos, [ -2, 2],   [1j, 1j], 1, -1, True, np.complex64
        yield _check_branch_cut, np.arctan, [0-2j, 2j],  [1,  1], -1, 1, True, np.complex64

        yield _check_branch_cut, np.arcsinh, [0-2j,  2j], [1,   1], -1, 1, True, np.complex64
        yield _check_branch_cut, np.arccosh, [ -1, 0.5], [1j,  1j], 1, -1, True, np.complex64
        yield _check_branch_cut, np.arctanh, [ -2,   2], [1j, 1j], 1, -1, True, np.complex64

        # check against bogus branch cuts: assert continuity between quadrants
        yield _check_branch_cut, np.arcsin, [0-2j, 2j], [ 1,  1], 1, 1, False, np.complex64
        yield _check_branch_cut, np.arccos, [0-2j, 2j], [ 1,  1], 1, 1, False, np.complex64
        yield _check_branch_cut, np.arctan, [ -2,  2], [1j, 1j], 1, 1, False, np.complex64

        yield _check_branch_cut, np.arcsinh, [ -2,  2, 0], [1j, 1j, 1], 1, 1, False, np.complex64
        yield _check_branch_cut, np.arccosh, [0-2j, 2j, 2], [1,  1,  1j], 1, 1, False, np.complex64
        yield _check_branch_cut, np.arctanh, [0-2j, 2j, 0], [1,  1,  1j], 1, 1, False, np.complex64 

def test_against_cmath(self):
        import cmath

        points = [-1-1j, -1+1j, +1-1j, +1+1j]
        name_map = {'arcsin': 'asin', 'arccos': 'acos', 'arctan': 'atan',
                    'arcsinh': 'asinh', 'arccosh': 'acosh', 'arctanh': 'atanh'}
        atol = 4*np.finfo(np.complex).eps
        for func in self.funcs:
            fname = func.__name__.split('.')[-1]
            cname = name_map.get(fname, fname)
            try:
                cfunc = getattr(cmath, cname)
            except AttributeError:
                continue
            for p in points:
                a = complex(func(np.complex_(p)))
                b = cfunc(p)
                assert_(abs(a - b) < atol, "%s %s: %s; cmath: %s" % (fname, p, a, b)) 

def get_direction(self):
        tau = [0, 0, 0]
        theta = [0, 0, 0]

        buf = b''.join(self.queue)
        buf = np.fromstring(buf, dtype='int16')
        for i, v in enumerate(self.pair):
            tau[i], _ = gcc_phat(buf[v[0]::8], buf[v[1]::8], fs=self.sample_rate, max_tau=MAX_TDOA_6P1,
                                 interp=1)
            theta[i] = np.arcsin(tau[i] / MAX_TDOA_6P1) * 180 / np.pi

        min_index = np.argmin(np.abs(tau))
        if (min_index != 0 and theta[min_index - 1] >= 0) or (min_index == 0 and theta[len(self.pair) - 1] < 0):
            best_guess = (theta[min_index] + 360) % 360
        else:
            best_guess = (180 - theta[min_index])

        best_guess = (best_guess + 120 + min_index * 60) % 360

        return best_guess 

def get_direction(self):
        tau = [0, 0, 0]
        theta = [0, 0, 0]

        buf = b''.join(self.queue)
        buf = np.fromstring(buf, dtype='int16')
        for i, v in enumerate(self.pair):
            tau[i], _ = gcc_phat(buf[v[0]::8], buf[v[1]::8], fs=self.sample_rate, max_tau=MAX_TDOA_6,
                                 interp=1)
            theta[i] = np.arcsin(tau[i] / MAX_TDOA_6) * 180 / np.pi

        min_index = np.argmin(np.abs(tau))
        if (min_index != 0 and theta[min_index - 1] >= 0) or (min_index == 0 and theta[len(self.pair) - 1] < 0):
            best_guess = (theta[min_index] + 360) % 360
        else:
            best_guess = (180 - theta[min_index])

        best_guess = (best_guess + 30 + min_index * 60) % 360

        return best_guess 

def find_bragg(lambd, lattice, ord_max):
    H = {}
    d = {}
    phi = {}
    
    for h in range(0, ord_max+1):
        for k in range(0, ord_max+1):
            for l in range(0, ord_max+1):
                
                if h == k == l == 0:
                    continue
                            
                hkl = (h,k,l)
    
                H_hkl = h*lattice.b1 + k*lattice.b2 + l*lattice.b3
                d_hkl = 1. / np.linalg.norm(H_hkl)
                                        
                sin_phi = lambd / (2. * d_hkl)
                
                if np.abs(sin_phi) <= 1.0:
                    H[hkl] = H_hkl
                    d[hkl] = d_hkl
                    phi[hkl] = np.arcsin(sin_phi) * 180.0 / np.pi

    return H, d, phi 

def sun_topo_azimuth_zenith(latitude, delta_prime, H_prime, temperature=14.6, pressure=1013):
        """Compute the sun's topocentric azimuth and zenith angles
        azimuth is measured eastward from north, zenith from vertical
        temperature = average temperature in C (default is 14.6 = global average in 2013)
        pressure = average pressure in mBar (default 1013 = global average)
        """
        phi = np.deg2rad(latitude)
        dr, Hr = map(np.deg2rad,(delta_prime, H_prime))
        P, T = pressure, temperature
        e0 = np.rad2deg(np.arcsin(np.sin(phi)*np.sin(dr) + np.cos(phi)*np.cos(dr)*np.cos(Hr)))
        tmp = np.deg2rad(e0 + 10.3/(e0+5.11))
        delta_e = (P/1010.0)*(283.0/(273+T))*(1.02/(60*np.tan(tmp)))
        e = e0 + delta_e
        zenith = 90 - e

        gamma = np.rad2deg(np.arctan2(np.sin(Hr), np.cos(Hr)*np.sin(phi) - np.tan(dr)*np.cos(phi))) % 360
        Phi = (gamma + 180) % 360 #azimuth from north
        return Phi, zenith 

def comp_alpha(self):
    """The opening angle with a W3 teeth width and Rbo - H0 radius

    Parameters
    ----------
    self : HoleM52
        A HoleM52 object

    Returns
    -------
    alpha: float
        Angle between P1 and P9 (cf schematics) [rad]

    """
    Rbo = self.get_Rbo()

    alpha_tooth = 2 * arcsin(self.W3 / (2 * (Rbo - self.H0)))
    slot_pitch = 2 * pi / self.Zh

    return slot_pitch - alpha_tooth 

def comp_angle_opening(self):
    """Compute the average opening angle of the Slot

    Parameters
    ----------
    self : SlotCirc
        A SlotCirc object

    Returns
    -------
    alpha: float
        Average opening angle of the slot [rad]
    """

    Rbo = self.get_Rbo()

    return float(2 * arcsin(self.W0 / (2 * Rbo))) 

def comp_angle_opening(self):
    """Compute the average opening angle of the Slot

    Parameters
    ----------
    self : SlotW21
        A SlotW21 object

    Returns
    -------
    alpha: float
        Average opening angle of the slot [rad]

    """

    Rbo = self.get_Rbo()

    return float(2 * arcsin(self.W0 / (2 * Rbo))) 

def spherical_distance(pt0, pt1):
    '''
    spherical_distance(a, b) yields the angular distance between points a and b, both of which
      should be expressed in spherical coordinates as (longitude, latitude).
    If a and/or b are (2 x n) matrices, then the calculation is performed over all columns.
    The spherical_distance function uses the Haversine formula; accordingly it may suffer from
    rounding errors in the case of nearly antipodal points.
    '''
    dtheta = pt1[0] - pt0[0]
    dphi   = pt1[1] - pt0[1]
    a = np.sin(dphi/2)**2 + np.cos(pt0[1]) * np.cos(pt1[1]) * np.sin(dtheta/2)**2
    return 2 * np.arcsin(np.sqrt(a)) 

def max_dmag_filter(self):
        """Includes stars if maximum delta mag is in the allowed orbital range
        
        Removed from prototype filters. Prototype is already calling the 
        int_cutoff_filter with OS.dMag0 and the completeness_filter with Comp.dMagLim
        
        """
        
        PPop = self.PlanetPopulation
        PPMod = self.PlanetPhysicalModel
        Comp = self.Completeness
        
        # s and beta arrays
        s = np.tan(self.OpticalSystem.WA0)*self.dist
        if PPop.scaleOrbits:
            s /= np.sqrt(self.L)
        beta = np.array([1.10472881476178]*len(s))*u.rad
        
        # fix out of range values
        below = np.where(s < np.min(PPop.rrange)*np.sin(beta))[0]
        above = np.where(s > np.max(PPop.rrange)*np.sin(beta))[0]
        s[below] = np.sin(beta[below])*np.min(PPop.rrange)
        beta[above] = np.arcsin(s[above]/np.max(PPop.rrange))
        
        # calculate delta mag
        p = np.max(PPop.prange)
        Rp = np.max(PPop.Rprange)
        d = s/np.sin(beta)
        Phi = PPMod.calc_Phi(beta)
        i = np.where(deltaMag(p, Rp, d, Phi) < Comp.dMagLim)[0]
        self.revise_lists(i) 

def _evaluate(self, X, out, *args, **kwargs):
        g = self.g2(X)
        f = g.reshape((-1, 1)) * np.ones((X.shape[0], self.n_obj))
        f[:, 1:] *= np.sin(0.5 * np.pi * X[:, (self.n_obj - 2)::-1])
        cos = np.cos(0.5 * np.pi * X[:, :(self.n_obj - 1)])
        f[:, 0:-1] *= np.flip(np.cumprod(cos, axis=1), axis=1)

        f_squared = (f ** 2).sum(axis=1)
        g0 = f_squared - (1.25 - self.LA2(0.5, 6.0, 1.0, 2.0, np.arcsin(f[:, -1] / np.sqrt(f_squared)))) * (
                1.25 - self.LA2(0.5, 6.0, 1.0, 2.0, np.arcsin(f[:, -1] / np.sqrt(f_squared))))
        out["F"] = f
        out["G"] = g0.reshape((-1, 1)) 

def _calc_pareto_front(self, ref_dirs):
        F = ref_dirs
        F = F / np.sqrt(np.sum(F ** 2, axis=1)).reshape((-1, 1))
        c = (1.25 - 0.5 * np.sin(6 * np.arcsin(F[:, -1])) ** 2) ** 2 - np.sum(F ** 2, axis=1)
        return F[c >= 0] 

def make20(b):
    theta1 = numpy.arccos(numpy.sqrt(5)/3)
    theta2 = numpy.arcsin(r2edge(theta1,1)/2/numpy.sin(numpy.pi/5))
    r = b/2./numpy.sin(theta1/2)
    rot72 = rotmatz(numpy.pi*2/5)
    s2 = numpy.sin(theta2)
    c2 = numpy.cos(theta2)
    s3 = numpy.sin(theta1+theta2)
    c3 = numpy.cos(theta1+theta2)
    p1 = numpy.array(( s2*r,  0,  c2*r))
    p2 = numpy.array(( s3*r,  0,  c3*r))
    p3 = numpy.array((-s3*r,  0, -c3*r))
    p4 = numpy.array((-s2*r,  0, -c2*r))
    coord = []
    for i in range(5):
        coord.append(p1)
        p1 = numpy.dot(p1, rot72)
    for i in range(5):
        coord.append(p2)
        p2 = numpy.dot(p2, rot72)
    for i in range(5):
        coord.append(p3)
        p3 = numpy.dot(p3, rot72)
    for i in range(5):
        coord.append(p4)
        p4 = numpy.dot(p4, rot72)
    return numpy.array(coord) 

def make20(b):
    theta1 = numpy.arccos(numpy.sqrt(5)/3)
    theta2 = numpy.arcsin(r2edge(theta1,1)/2/numpy.sin(numpy.pi/5))
    r = b/2/numpy.sin(theta1/2)
    rot72 = rotmatz(numpy.pi*2/5)
    s2 = numpy.sin(theta2)
    c2 = numpy.cos(theta2)
    s3 = numpy.sin(theta1+theta2)
    c3 = numpy.cos(theta1+theta2)
    p1 = numpy.array(( s2*r,  0,  c2*r))
    p2 = numpy.array(( s3*r,  0,  c3*r))
    p3 = numpy.array((-s3*r,  0, -c3*r))
    p4 = numpy.array((-s2*r,  0, -c2*r))
    coord = []
    for i in range(5):
        coord.append(p1)
        p1 = numpy.dot(p1, rot72)
    for i in range(5):
        coord.append(p2)
        p2 = numpy.dot(p2, rot72)
    for i in range(5):
        coord.append(p3)
        p3 = numpy.dot(p3, rot72)
    for i in range(5):
        coord.append(p4)
        p4 = numpy.dot(p4, rot72)
    return numpy.array(coord) 

def assignVanishingType(lines, vp, tol, area=10):
    numLine = len(lines)
    numVP = len(vp)
    typeCost = np.zeros((numLine, numVP))
    # perpendicular
    for vid in range(numVP):
        cosint = (lines[:, :3] * vp[[vid]]).sum(1)
        typeCost[:, vid] = np.arcsin(np.abs(cosint).clip(-1, 1))

    # infinity
    u = np.stack([lines[:, 4], lines[:, 5]], -1)
    u = u.reshape(-1, 1) * 2 * np.pi - np.pi
    v = computeUVN_vec(lines[:, :3], u, lines[:, 3])
    xyz = uv2xyzN_vec(np.hstack([u, v]), np.repeat(lines[:, 3], 2))
    xyz = multi_linspace(xyz[0::2].reshape(-1), xyz[1::2].reshape(-1), 100)
    xyz = np.vstack([blk.T for blk in np.split(xyz, numLine)])
    xyz = xyz / np.linalg.norm(xyz, axis=1, keepdims=True)
    for vid in range(numVP):
        ang = np.arccos(np.abs((xyz * vp[[vid]]).sum(1)).clip(-1, 1))
        notok = (ang < area * np.pi / 180).reshape(numLine, 100).sum(1) != 0
        typeCost[notok, vid] = 100

    I = typeCost.min(1)
    tp = typeCost.argmin(1)
    tp[I > tol] = numVP + 1

    return tp, typeCost 

def numpy_arcsin(a):
  return np.arcsin(a) 

def arcsin(y, x):
  d[x] = d[y] / numpy.sqrt(1.0 - x * x) 

def beam_block_frac(Th, Bh, a):
    """Partial beam blockage fraction. Unitless.

    From Bech et al. (2003), Eqn 2 and Appendix

    Parameters
    ----------
    Th : float
        Terrain height [m]
    Bh : float
        Beam height [m]
    a : float
        Half power beam radius [m]

    Notes
    -----
    This procedure uses a simplified interception function where no vertical
      gradient of refractivity is considered.  Other algorithms treat this
      more thoroughly.  However, this is accurate in most cases other than
      the super-refractive case.

    See the the half_power_radius function to calculate variable a

    The heights must be the same units!
    """

    # First find the difference between the terrain and height of
    # radar beam (Bech et al. (2003), Fig.3)
    y = Th - Bh

    Numer = (y * np.sqrt(a**2 - y**2)) + (a**2 * np.arcsin(y/a)) + (np.pi * a**2 /2.)

    Denom = np.pi * a**2

    PBB = Numer / Denom

    return PBB