Python numpy.tan() Examples

def set_camera(self, fov_vertical, z_near, z_far, aspect):
    width = 2*np.tan(np.deg2rad(fov_vertical)/2.0)*z_near*aspect;
    height = 2*np.tan(np.deg2rad(fov_vertical)/2.0)*z_near;
    egl_program = self.egl_program
    c = np.eye(4, dtype=np.float32)
    c[3,3] = 0
    c[3,2] = -1
    c[2,2] = -(z_near+z_far)/(z_far-z_near)
    c[2,3] = -2.0*(z_near*z_far)/(z_far-z_near)
    c[0,0] = 2.0*z_near/width
    c[1,1] = 2.0*z_near/height
    c = c.T
    
    projection_matrix_o = glGetUniformLocation(egl_program, 'uProjectionMatrix')
    projection_matrix = np.eye(4, dtype=np.float32)
    projection_matrix[...] = c
    projection_matrix = np.reshape(projection_matrix, (-1))
    glUniformMatrix4fv(projection_matrix_o, 1, GL_FALSE, projection_matrix) 

def cor2xybound(cor):
    ''' Helper function to clip max/min stretch factor '''
    corU = cor[0::2]
    corB = cor[1::2]
    zU = -50
    u = panostretch.coorx2u(corU[:, 0])
    vU = panostretch.coory2v(corU[:, 1])
    vB = panostretch.coory2v(corB[:, 1])

    x, y = panostretch.uv2xy(u, vU, z=zU)
    c = np.sqrt(x**2 + y**2)
    zB = c * np.tan(vB)
    xmin, xmax = x.min(), x.max()
    ymin, ymax = y.min(), y.max()

    S = 3 / abs(zB.mean() - zU)
    dx = [abs(xmin * S), abs(xmax * S)]
    dy = [abs(ymin * S), abs(ymax * S)]

    return min(dx), min(dy), max(dx), max(dy) 

def test_triangle_angles(self):
        msh = mesh_io.Msh()
        # 2 regular tetrahedron of edge size 1
        msh.elm = mesh_io.Elements(
            triangles=np.array(
                [[1, 2, 3],
                 [1, 2, 3],
                 [1, 4, 2],
                 [1, 4, 2]], dtype=int))
        msh.nodes = mesh_io.Nodes(np.array(
            [[0, 0, 0],
             [1, 0, 0],
             [0, np.tan(np.pi/6), 0],
             [0, 0, np.tan(np.pi/6)]], dtype=float))
        angles = msh.triangle_angles()
        assert np.allclose(angles[:3], [90, 30, 60])
        assert np.allclose(angles[3:], [90, 60, 30]) 

def vehicle_flat_reverse(zflag, params={}):
    # Get the parameter values
    b = params.get('wheelbase', 3.)

    # Create a vector to store the state and inputs
    x = np.zeros(3)
    u = np.zeros(2)

    # Given the flat variables, solve for the state
    x[0] = zflag[0][0]  # x position
    x[1] = zflag[1][0]  # y position
    x[2] = np.arctan2(zflag[1][1], zflag[0][1])  # tan(theta) = ydot/xdot

    # And next solve for the inputs
    u[0] = zflag[0][1] * np.cos(x[2]) + zflag[1][1] * np.sin(x[2])
    thdot_v = zflag[1][2] * np.cos(x[2]) - zflag[0][2] * np.sin(x[2])
    u[1] = np.arctan2(thdot_v, u[0]**2 / b)

    return x, u


# Function to compute the RHS of the system dynamics 

def vehicle_flat_reverse(zflag, params={}):
    # Get the parameter values
    b = params.get('wheelbase', 3.)

    # Create a vector to store the state and inputs
    x = np.zeros(3)
    u = np.zeros(2)

    # Given the flat variables, solve for the state
    x[0] = zflag[0][0]  # x position
    x[1] = zflag[1][0]  # y position
    x[2] = np.arctan2(zflag[1][1], zflag[0][1])  # tan(theta) = ydot/xdot

    # And next solve for the inputs
    u[0] = zflag[0][1] * np.cos(x[2]) + zflag[1][1] * np.sin(x[2])
    thdot_v = zflag[1][2] * np.cos(x[2]) - zflag[0][2] * np.sin(x[2])
    u[1] = np.arctan2(thdot_v, u[0]**2 / b)

    return x, u


# Function to compute the RHS of the system dynamics 

def perspectiveprojectionnp(fovy, ratio=1.0, near=0.01, far=10.0):
    
    tanfov = np.tan(fovy / 2.0)
    # top = near * tanfov
    # right = ratio * top
    # mtx = [near / right, 0, 0, 0, \
    #          0, near / top, 0, 0, \
    #          0, 0, -(far+near)/(far-near), -2*far*near/(far-near), \
    #          0, 0, -1, 0]
    mtx = [[1.0 / (ratio * tanfov), 0, 0, 0], \
                [0, 1.0 / tanfov, 0, 0], \
                [0, 0, -(far + near) / (far - near), -2 * far * near / (far - near)], \
                [0, 0, -1.0, 0]]
    # return np.array(mtx, dtype=np.float32)
    return np.array([[1.0 / (ratio * tanfov)], [1.0 / tanfov], [-1]], dtype=np.float32)


##################################################### 

def makeArrowPath(headLen=20, tipAngle=20, tailLen=20, tailWidth=3, baseAngle=0):
    """
    Construct a path outlining an arrow with the given dimensions.
    The arrow points in the -x direction with tip positioned at 0,0.
    If *tipAngle* is supplied (in degrees), it overrides *headWidth*.
    If *tailLen* is None, no tail will be drawn.
    """
    headWidth = headLen * np.tan(tipAngle * 0.5 * np.pi/180.)
    path = QtGui.QPainterPath()
    path.moveTo(0,0)
    path.lineTo(headLen, -headWidth)
    if tailLen is None:
        innerY = headLen - headWidth * np.tan(baseAngle*np.pi/180.)
        path.lineTo(innerY, 0)
    else:
        tailWidth *= 0.5
        innerY = headLen - (headWidth-tailWidth) * np.tan(baseAngle*np.pi/180.)
        path.lineTo(innerY, -tailWidth)
        path.lineTo(headLen + tailLen, -tailWidth)
        path.lineTo(headLen + tailLen, tailWidth)
        path.lineTo(innerY, tailWidth)
    path.lineTo(headLen, headWidth)
    path.lineTo(0,0)
    return path 

def projectionMatrix(self, region=None):
        # Xw = (Xnd + 1) * width/2 + X
        if region is None:
            region = (0, 0, self.width(), self.height())
        
        x0, y0, w, h = self.getViewport()
        dist = self.opts['distance']
        fov = self.opts['fov']
        nearClip = dist * 0.001
        farClip = dist * 1000.

        r = nearClip * np.tan(fov * 0.5 * np.pi / 180.)
        t = r * h / w

        # convert screen coordinates (region) to normalized device coordinates
        # Xnd = (Xw - X0) * 2/width - 1
        ## Note that X0 and width in these equations must be the values used in viewport
        left  = r * ((region[0]-x0) * (2.0/w) - 1)
        right = r * ((region[0]+region[2]-x0) * (2.0/w) - 1)
        bottom = t * ((region[1]-y0) * (2.0/h) - 1)
        top    = t * ((region[1]+region[3]-y0) * (2.0/h) - 1)

        tr = QtGui.QMatrix4x4()
        tr.frustum(left, right, bottom, top, nearClip, farClip)
        return tr 

def pan(self, dx, dy, dz, relative=False):
        """
        Moves the center (look-at) position while holding the camera in place. 
        
        If relative=True, then the coordinates are interpreted such that x
        if in the global xy plane and points to the right side of the view, y is
        in the global xy plane and orthogonal to x, and z points in the global z
        direction. Distances are scaled roughly such that a value of 1.0 moves
        by one pixel on screen.
        
        """
        if not relative:
            self.opts['center'] += QtGui.QVector3D(dx, dy, dz)
        else:
            cPos = self.cameraPosition()
            cVec = self.opts['center'] - cPos
            dist = cVec.length()  ## distance from camera to center
            xDist = dist * 2. * np.tan(0.5 * self.opts['fov'] * np.pi / 180.)  ## approx. width of view at distance of center point
            xScale = xDist / self.width()
            zVec = QtGui.QVector3D(0,0,1)
            xVec = QtGui.QVector3D.crossProduct(zVec, cVec).normalized()
            yVec = QtGui.QVector3D.crossProduct(xVec, zVec).normalized()
            self.opts['center'] = self.opts['center'] + xVec * xScale * dx + yVec * xScale * dy + zVec * xScale * dz
        self.update() 

def tangent(x, null=(-np.inf, np.inf), rtol=default_rtol, atol=default_atol):
    '''
    tangent(x) is equivalent to tan(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 == -pi/2 or -pi/2. If only one number is given, then it is
    used for both values; otherwise the first value corresponds to -pi/2 and the second to pi/2.
    A value of x is considered to be equal to one of these valids based on numpy.isclose. The
    optional arguments rtol and atol are passed along to isclose. If null is None, then no
    replacement is performed.
    '''
    if sps.issparse(x):
        x = x.copy()
        x.data = tangent(x.data, null=null, rtol=rtol, atol=atol)
        return x
    else: x = np.asarray(x)
    if rtol is None: rtol = default_rtol
    if atol is None: atol = default_atol
    try:    (nln,nlp) = null
    except Exception: (nln,nlp) = (null,null)
    x = np.mod(x + pi, tau) - pi
    ii = None if nln is None else np.where(np.isclose(x, neghpi, rtol=rtol, atol=atol))
    jj = None if nlp is None else np.where(np.isclose(x, hpi,    rtol=rtol, atol=atol))
    x = np.tan(x)
    if ii: x[ii] = nln
    if jj: x[jj] = nlp
    return x 

def cotangent(x, null=(-np.inf, np.inf), rtol=default_rtol, atol=default_atol):
    '''
    cotangent(x) is equivalent to cot(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 == 0 or pi. If only one number is given, then it is used for
    both values; otherwise the first value corresponds to 0 and the second to pi.  A value of x is
    considered to be equal to one of these valids based on numpy.isclose. The optional arguments
    rtol and atol are passed along to isclose. If null is None, then no replacement is performed.
    '''
    if sps.issparse(x): x = x.toarray()
    else:               x = np.asarray(x)
    if rtol is None: rtol = default_rtol
    if atol is None: atol = default_atol
    try:    (nln,nlp) = null
    except Exception: (nln,nlp) = (null,null)
    x = np.mod(x + hpi, tau) - hpi
    ii = None if nln is None else np.where(np.isclose(x, 0,  rtol=rtol, atol=atol))
    jj = None if nlp is None else np.where(np.isclose(x, pi, rtol=rtol, atol=atol))
    x = np.tan(x)
    if ii: x[ii] = 1
    if jj: x[jj] = 1
    x = 1.0 / x
    if ii: x[ii] = nln
    if jj: x[jj] = nlp
    return x 

def get_camera_matrix(width, height, fov):
  """Returns a camera matrix from image size and fov."""
  xc = (width-1.) / 2.
  zc = (height-1.) / 2.
  f = (width / 2.) / np.tan(np.deg2rad(fov / 2.))
  camera_matrix = utils.Foo(xc=xc, zc=zc, f=f)
  return camera_matrix 

def render_nodes(self, nodes, perturb=None, aux_delta_theta=0.):
    self.set_building_visibility(True)
    if perturb is None:
      perturb = np.zeros((len(nodes), 4))

    imgs = []
    r = 2
    elevation_z = r * np.tan(np.deg2rad(self.robot.camera_elevation_degree))

    for i in range(len(nodes)):
      xyt = self.to_actual_xyt(nodes[i])
      lookat_theta = 3.0 * np.pi / 2.0 - (xyt[2]+perturb[i,2]+aux_delta_theta) * (self.task.delta_theta)
      nxy = np.array([xyt[0]+perturb[i,0], xyt[1]+perturb[i,1]]).reshape(1, -1)
      nxy = nxy * self.map.resolution
      nxy = nxy + self.map.origin
      camera_xyz = np.zeros((1, 3))
      camera_xyz[...] = [nxy[0, 0], nxy[0, 1], self.robot.sensor_height]
      camera_xyz = camera_xyz / 100.
      lookat_xyz = np.array([-r * np.sin(lookat_theta),
                             -r * np.cos(lookat_theta), elevation_z])
      lookat_xyz = lookat_xyz + camera_xyz[0, :]
      self.r_obj.position_camera(camera_xyz[0, :].tolist(),
                                 lookat_xyz.tolist(), [0.0, 0.0, 1.0])
      img = self.r_obj.render(take_screenshot=True, output_type=0)
      img = [x for x in img if x is not None]
      img = np.concatenate(img, axis=2).astype(np.float32)
      if perturb[i,3]>0:
        img = img[:,::-1,:]
      imgs.append(img)

    self.set_building_visibility(False)
    return imgs 

def test_target_completeness_def(self):
        """
        Compare calculated completenesses for multiple targets under default population
        settings.
        """
            
        with RedirectStreams(stdout=self.dev_null):
            TL = TargetList(ntargs=100,**copy.deepcopy(self.spec))

            mode = list(filter(lambda mode: mode['detectionMode'] == True, TL.OpticalSystem.observingModes))[0]
            IWA = mode['IWA']
            OWA = mode['OWA']
            rrange = TL.PlanetPopulation.rrange
            maxd = (rrange[1]/np.tan(IWA)).to(u.pc).value
            mind = (rrange[0]/np.tan(OWA)).to(u.pc).value

            #want distances to span from outer edge below IWA to inner edge above OWA
            TL.dist = np.logspace(np.log10(mind/10.),np.log10(maxd*10.),TL.nStars)*u.pc


        Brown = EXOSIMS.Completeness.BrownCompleteness.BrownCompleteness(**copy.deepcopy(self.spec))
        Garrett = EXOSIMS.Completeness.GarrettCompleteness.GarrettCompleteness(**copy.deepcopy(self.spec))

        cBrown = Brown.target_completeness(TL)
        cGarrett = Garrett.target_completeness(TL)
        
        np.testing.assert_allclose(cGarrett,cBrown,rtol=0.1,atol=1e-6)
        
        # test when scaleOrbits == True
        TL.L = np.exp(np.random.uniform(low=np.log(0.1), high=np.log(10.), size=TL.nStars))
        Brown.PlanetPopulation.scaleOrbits = True
        Garrett.PlanetPopulation.scaleOrbits = True

        cBrown = Brown.target_completeness(TL)
        cGarrett = Garrett.target_completeness(TL)

        cGarrett = cGarrett[cBrown != 0 ]
        cBrown = cBrown[cBrown != 0]
        meandiff = np.mean(np.abs(cGarrett - cBrown)/cBrown)

        self.assertLessEqual(meandiff,0.1) 

def test_target_completeness_constrainOrbits(self):
        """
        Compare calculated completenesses for multiple targets with constrain orbits set to true
        """
            
        
        with RedirectStreams(stdout=self.dev_null):
            TL = TargetList(ntargs=100,constrainOrbits=True,**copy.deepcopy(self.spec))

            mode = list(filter(lambda mode: mode['detectionMode'] == True, TL.OpticalSystem.observingModes))[0]
            IWA = mode['IWA']
            OWA = mode['OWA']
            rrange = TL.PlanetPopulation.rrange
            maxd = (rrange[1]/np.tan(IWA)).to(u.pc).value
            mind = (rrange[0]/np.tan(OWA)).to(u.pc).value

            #want distances to span from outer edge below IWA to inner edge above OWA
            TL.dist = np.logspace(np.log10(mind/10.),np.log10(maxd*10.),TL.nStars)*u.pc


        Brown = EXOSIMS.Completeness.BrownCompleteness.BrownCompleteness(constrainOrbits=True,**copy.deepcopy(self.spec))
        Garrett = EXOSIMS.Completeness.GarrettCompleteness.GarrettCompleteness(constrainOrbits=True,**copy.deepcopy(self.spec))

        cBrown = Brown.target_completeness(TL)
        cGarrett = Garrett.target_completeness(TL)

        np.testing.assert_allclose(cGarrett,cBrown,rtol=0.1,atol=1e-6) 

def outside_IWA_filter(self):
        """Includes stars with planets with orbits outside of the IWA 
        
        """
        
        PPop = self.PlanetPopulation
        OS = self.OpticalSystem
        
        s = np.tan(OS.IWA)*self.dist
        L = np.sqrt(self.L) if PPop.scaleOrbits else 1.
        i = np.where(s < L*np.max(PPop.rrange))[0]
        self.revise_lists(i) 

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 xyzpers(h_fov, v_fov, u, v, out_hw, in_rot):
    out = np.ones((*out_hw, 3), np.float32)

    x_max = np.tan(h_fov / 2)
    y_max = np.tan(v_fov / 2)
    x_rng = np.linspace(-x_max, x_max, num=out_hw[1], dtype=np.float32)
    y_rng = np.linspace(-y_max, y_max, num=out_hw[0], dtype=np.float32)
    out[..., :2] = np.stack(np.meshgrid(x_rng, -y_rng), -1)
    Rx = rotation_matrix(v, [1, 0, 0])
    Ry = rotation_matrix(u, [0, 1, 0])
    Ri = rotation_matrix(in_rot, np.array([0, 0, 1.0]).dot(Rx).dot(Ry))

    return out.dot(Rx).dot(Ry).dot(Ri) 

def calc_rho(beta):
    return 1 / (1 - np.tan(beta / 360 * 2 * np.pi)) - 1 

def layout_2_depth(cor_id, h, w, return_mask=False):
    # Convert corners to per-column boundary first
    # Up -pi/2,  Down pi/2
    vc, vf = cor_2_1d(cor_id, h, w)
    vc = vc[None, :]  # [1, w]
    vf = vf[None, :]  # [1, w]
    assert (vc > 0).sum() == 0
    assert (vf < 0).sum() == 0

    # Per-pixel v coordinate (vertical angle)
    vs = ((np.arange(h) + 0.5) / h - 0.5) * np.pi
    vs = np.repeat(vs[:, None], w, axis=1)  # [h, w]

    # Floor-plane to depth
    floor_h = 1.6
    floor_d = np.abs(floor_h / np.sin(vs))

    # wall to camera distance on horizontal plane at cross camera center
    cs = floor_h / np.tan(vf)

    # Ceiling-plane to depth
    ceil_h = np.abs(cs * np.tan(vc))      # [1, w]
    ceil_d = np.abs(ceil_h / np.sin(vs))  # [h, w]

    # Wall to depth
    wall_d = np.abs(cs / np.cos(vs))  # [h, w]

    # Recover layout depth
    floor_mask = (vs > vf)
    ceil_mask = (vs < vc)
    wall_mask = (~floor_mask) & (~ceil_mask)
    depth = np.zeros([h, w], np.float32)    # [h, w]
    depth[floor_mask] = floor_d[floor_mask]
    depth[ceil_mask] = ceil_d[ceil_mask]
    depth[wall_mask] = wall_d[wall_mask]

    assert (depth == 0).sum() == 0
    if return_mask:
        return depth, floor_mask, ceil_mask, wall_mask
    return depth 

def edgeFromImg2Pano(edge):
    edgeList = edge['edgeLst']
    if len(edgeList) == 0:
        return np.array([])

    vx = edge['vx']
    vy = edge['vy']
    fov = edge['fov']
    imH, imW = edge['img'].shape

    R = (imW/2) / np.tan(fov/2)

    # im is the tangent plane, contacting with ball at [x0 y0 z0]
    x0 = R * np.cos(vy) * np.sin(vx)
    y0 = R * np.cos(vy) * np.cos(vx)
    z0 = R * np.sin(vy)
    vecposX = np.array([np.cos(vx), -np.sin(vx), 0])
    vecposY = np.cross(np.array([x0, y0, z0]), vecposX)
    vecposY = vecposY / np.sqrt(vecposY @ vecposY.T)
    vecposX = vecposX.reshape(1, -1)
    vecposY = vecposY.reshape(1, -1)
    Xc = (0 + imW-1) / 2
    Yc = (0 + imH-1) / 2

    vecx1 = edgeList[:, [0]] - Xc
    vecy1 = edgeList[:, [1]] - Yc
    vecx2 = edgeList[:, [2]] - Xc
    vecy2 = edgeList[:, [3]] - Yc

    vec1 = np.tile(vecx1, [1, 3]) * vecposX + np.tile(vecy1, [1, 3]) * vecposY
    vec2 = np.tile(vecx2, [1, 3]) * vecposX + np.tile(vecy2, [1, 3]) * vecposY
    coord1 = [[x0, y0, z0]] + vec1
    coord2 = [[x0, y0, z0]] + vec2

    normal = np.cross(coord1, coord2, axis=1)
    normal = normal / np.linalg.norm(normal, axis=1, keepdims=True)

    panoList = np.hstack([normal, coord1, coord2, edgeList[:, [-1]]])

    return panoList 

def np_coor2xy(coor, z=50, coorW=1024, coorH=512, floorW=1024, floorH=512):
    '''
    coor: N x 2, index of array in (col, row) format
    '''
    coor = np.array(coor)
    u = np_coorx2u(coor[:, 0], coorW)
    v = np_coory2v(coor[:, 1], coorH)
    c = z / np.tan(v)
    x = c * np.sin(u) + floorW / 2 - 0.5
    y = -c * np.cos(u) + floorH / 2 - 0.5
    return np.hstack([x[:, None], y[:, None]]) 

def np_refine_by_fix_z(coory0, coory1, z0=50, coorH=512):
    '''
    Refine coory1 by coory0
    coory0 are assumed on given plane z
    '''
    v0 = np_coory2v(coory0, coorH)
    v1 = np_coory2v(coory1, coorH)

    c0 = z0 / np.tan(v0)
    z1 = c0 * np.tan(v1)
    z1_mean = mean_percentile(z1)
    v1_refine = np.arctan2(z1_mean, c0)
    coory1_refine = (-v1_refine / PI + 0.5) * coorH - 0.5

    return coory1_refine, z1_mean 

def infer_coory(coory0, h, z0=50, coorH=512):
    v0 = np_coory2v(coory0, coorH)
    c0 = z0 / np.tan(v0)
    z1 = z0 + h
    v1 = np.arctan2(z1, c0)
    return (-v1 / PI + 0.5) * coorH - 0.5 

def _uv_tri(w, h):
    uv = uv_meshgrid(w, h)
    sin_u = np.sin(uv[..., 0])
    cos_u = np.cos(uv[..., 0])
    tan_v = np.tan(uv[..., 1])
    return sin_u, cos_u, tan_v 

def uv2xy(u, v, z=-50):
    c = z / np.tan(v)
    x = c * np.cos(u)
    y = c * np.sin(u)
    return x, y 

def pano_connect_points(p1, p2, z=-50, w=1024, h=512):
    if p1[0] == p2[0]:
        return np.array([p1, p2], np.float32)

    u1 = coorx2u(p1[0], w)
    v1 = coory2v(p1[1], h)
    u2 = coorx2u(p2[0], w)
    v2 = coory2v(p2[1], h)

    x1, y1 = uv2xy(u1, v1, z)
    x2, y2 = uv2xy(u2, v2, z)

    if abs(p1[0] - p2[0]) < w / 2:
        pstart = np.ceil(min(p1[0], p2[0]))
        pend = np.floor(max(p1[0], p2[0]))
    else:
        pstart = np.ceil(max(p1[0], p2[0]))
        pend = np.floor(min(p1[0], p2[0]) + w)
    coorxs = (np.arange(pstart, pend + 1) % w).astype(np.float64)
    vx = x2 - x1
    vy = y2 - y1
    us = coorx2u(coorxs, w)
    ps = (np.tan(us) * x1 - y1) / (vy - np.tan(us) * vx)
    cs = np.sqrt((x1 + ps * vx) ** 2 + (y1 + ps * vy) ** 2)
    vs = np.arctan2(z, cs)
    coorys = v2coory(vs, h)

    return np.stack([coorxs, coorys], axis=-1) 

def eval(self, coords, grad=False):
        v1 = (coords[self.i] - coords[self.j]) * angstrom
        v2 = (coords[self.k] - coords[self.j]) * angstrom
        dot_product = np.dot(v1, v2) / (norm(v1) * norm(v2))
        if dot_product < -1:
            dot_product = -1
        elif dot_product > 1:
            dot_product = 1
        phi = np.arccos(dot_product)
        if not grad:
            return phi
        if abs(phi) > pi - 1e-6:
            grad = [
                (pi - phi) / (2 * norm(v1) ** 2) * v1,
                (1 / norm(v1) - 1 / norm(v2)) * (pi - phi) / (2 * norm(v1)) * v1,
                (pi - phi) / (2 * norm(v2) ** 2) * v2,
            ]
        else:
            grad = [
                1 / np.tan(phi) * v1 / norm(v1) ** 2
                - v2 / (norm(v1) * norm(v2) * np.sin(phi)),
                (v1 + v2) / (norm(v1) * norm(v2) * np.sin(phi))
                - 1 / np.tan(phi) * (v1 / norm(v1) ** 2 + v2 / norm(v2) ** 2),
                1 / np.tan(phi) * v2 / norm(v2) ** 2
                - v1 / (norm(v1) * norm(v2) * np.sin(phi)),
            ]
        return phi, grad 

def numpy_tan(a):
  return np.tan(a) 

def tan(y, x):
  cx = numpy.cos(x)
  d[x] = d[y] / (cx * cx)