I’ve now included a Keplerian orbital model in the Glaciation model code, and used it to generate numerically mean daily solar power at latitude 65º N on summer solstice over 310 kyears using Earth orbit eccentricity, obliquity, and perihelion longitude at 1000 year intervals provided by J Laskar, and a solar constant of 1367 W/m² (from M Kharseh) at 1 AU radius from Sun..
The values calculated compare quite well with Laskar’s own figures (below, using Milankovitch Orbital Data Viewer).
The period covers the current interglacial period (present day shown as 0), and the two preceding ones. Laskar has orbital elements that extend back one million years.
So I’m now able to generate Milankovitch cycle solar power at all Earth latitudes over this period. I currently only use daily mean values of solar power, but I could also generate hourly values should I need to do so.