def METHOD_NAME( nwno, cos_theta, xint_at_top, gweight, tweight,F0PI): 
    """
    Last step in albedo code. Integrates over phase angle based on the 
    Gaussian-Chebychev weights in geometry.json 
    Parameters
    ----------
    nwno : int 
        Number of wavenumbers 
    cos_theta : float 
        Cosine of phase angle 
    xint_at_top : ndarray of floats 
        Planetary intensity at the top of the atmosphere with dimensions (ng, nt, nwno)
    gweight : ndarray of floats 
        Gaussian weights for integration 
    tweight : ndarray of floats 
        Chebychev weights for integration
    F0PI : ndarray of floats 
        Stellar flux 
    """
    ng = len(gweight)
    nt = len(tweight)
    if nt==1 : sym_fac = 2*pi #azimuthal symmetry  
    else: sym_fac = 1
    albedo=zeros(nwno)
    #for w in range(nwno):
    #   albedo[w] = 0.5 * sum((xint_at_top[:,:,w]*tweight).T*gweight)
    for ig in range(ng): 
        for it in range(nt): 
            albedo = albedo + xint_at_top[ig,it,:] * gweight[ig] * tweight[it]
    albedo = sym_fac * 0.5 * albedo /F0PI * (cos_theta + 1.0)
    return albedo