Simple energy balance models of planetary climate use the incoming solar radiation ({\em insolation}) as input. For these models, the mean annual insolation as a function of latitude is often the appropriate variable. This function also depends on the obliquity, or axial tilt, of the planet, which may change enough over long time spans to affect the climate. The insolation as a function of obliquity and latitude can be expressed as a definite integral, which can be approximated by a polynomial, cutting computing time significantly.

In this paper, we use a sixth order Legendre series expansion to approximate the mean annual insolation by latitude of a planet with obliquity angle $\beta$, leading to faster computations with little loss in the accuracy of results. We discuss differences between our method and selected computational results for insolation found in the literature.