How to use the pvlib.pvl_tools.cosd function in pvlib

plane of array irradiance and the horizontal irradiance. 
    
    Input all angles in degrees.
    
    :param surf_tilt: float or Series. Panel tilt from horizontal.
    :param surf_az: float or Series. Panel azimuth from north.
    :param sun_zen: float or Series. Solar zenith angle.
    :param sun_az: float or Series. Solar azimuth angle.
    
    :returns: float or Series. Ratio of the plane of array irradiance to the
              horizontal plane irradiance
    """
    
    cos_poa_zen = aoi_projection(surf_tilt, surf_az, sun_zen, sun_az)
    
    cos_sun_zen = pvl_tools.cosd(sun_zen)
    
    # ratio of titled and horizontal beam irradiance
    ratio = cos_poa_zen / cos_sun_zen
    
    try:
        ratio.name = 'poa_ratio'
    except AttributeError:
        pass
    
    return ratio

pvl_ephemeris   
    pvl_extraradiation   
    pvl_isotropicsky
    pvl_reindl1990   
    pvl_perez 
    pvl_klucher1979   
    pvl_kingdiffuse
    pvl_spa

    '''
    
    pvl_logger.debug('diffuse_sky.haydavies()')
    
    cos_tt = pvlib.planeofarray.aoi_projection(surf_tilt, surf_az, sun_zen, sun_az)
    
    cos_sun_zen = pvl_tools.cosd(sun_zen)
    
    # ratio of titled and horizontal beam irradiance
    # Will H: is this to avoid / 0? Shouldn't it be element-wise?
    #Rb = np.max(cos_tt, 0) / np.max(pvl_tools.cosd(sun_zen), 0.01745) 
    Rb = cos_tt / cos_sun_zen
    
    # Anisotropy Index
    AI = DNI / DNI_ET

    # these are actually the () and [] sub-terms of the second term of eqn 7
    term1 = 1 - AI
    term2 = 0.5 * (1 + pvl_tools.cosd(surf_tilt))

    sky_diffuse = DHI * ( AI*Rb + term1 * term2 )
    sky_diffuse[sky_diffuse < 0] = 0

F1 = F1.astype(float)

    F2 = F2c[ebin,0] + F2c[ebin,1]*delta + F2c[ebin,2]*z
    F2[F2 < 0] = 0
    F2 = F2.astype(float)

    A = pvlib.planeofarray.aoi_projection(surf_tilt, surf_az, sun_zen, sun_az)
    A[A < 0] = 0

    B = pvl_tools.cosd(sun_zen);
    B[B < pvl_tools.cosd(85)] = pvl_tools.cosd(85)


    #Calculate Diffuse POA from sky dome
    
    term1 = 0.5 * (1 - F1) * (1 + pvl_tools.cosd(surf_tilt))
    term2 = F1 * A[ebin.index] / B[ebin.index]
    term3 = F2*pvl_tools.sind(surf_tilt)
    
    sky_diffuse = DHI[ebin.index] * (term1 + term2 + term3)
    sky_diffuse[sky_diffuse < 0] = 0

    return sky_diffuse