def METHOD_NAME(self, beta, phiIndex=None):
    """Calculate the phase function. Prototype method uses the Lambert phase
    function from Sobolev 1975.
    Args:
        beta (~astropy.units.Quantity(~numpy.ndarray(float))):
            Planet phase angles at which the phase function is to be calculated,
            in units of rad
    Returns:
        ~numpy.ndarray(float):
            Planet phase function
    """
    beta = beta.to("rad").value
    Phi = (np.sin(beta) + (np.pi - beta) * np.cos(beta)) / np.pi
    return Phi