New, accurate numerical solutions of the radiative transfer equation are compared with the Hapke <1981, 1984, 1986> analytic approximation, which is widely used in planetary data analyses. The numerical solutions use the Ambartsumian invariance principle as do the well-known Chandrasekhar <1960> H function solutions. The invariance principle has been reexpressed in a form which allows high order-accurate numerical integrations without any required interpolations. The new numerical solutions reproduce the Chandrasekhar H function solutions for Legendre phase functions but also allow single-scattering phase functions of arbitrary form. Accurate numerical solutions of the radiative transfer equation for Henyey-Greenstein phase functions reveal that the errors in the Hapke model for estimating bidirectional reflectance values range from <2% rms for dark surfaces like the Moon to <7% rms for bright surfaces such as Europa. Further comparisons demonstrate that Hapke model fitting procedures estimate single particle scattering albedo values to within <3% for both dark and bright surfaces. ¿ 2000 American Geophysical Union