Code for solar simulations in PNAS paper "Ancient inhabitants of the Basin of Mexico kept an accurate agricultural calendar using sunrise observatories and mountain alignments."

To calculate the position of the sun at any specific day and time, all dates were first converted to Julian Day numbers, a continuous count of days from the beginning of year −4712. The solar coordinates for each specific time were calculated using standard astronomical algorithms, including (a) the equation of Kepler (eccentric anomaly), (b) solar longitude (corrected by the gravitational influence of the moon and planets), (c) eccentricity of the Earth’s orbit, (d) relative distance to the sun, (e) obliquity of the elliptic, (f) equation of time, and (g) solar declination. Solar angular elevation values were corrected for atmospheric refraction, assuming a mean atmospheric pressure of 76 kPa and an air temperature at sunrise of 10 °C. Finally, using spherical trigonometry, solar coordinates were converted to local azimuth and elevation coordinates for a specific site on the Earth’s surface. Simulating the position of the sun every second of the day and knowing the elevation of the landscape horizon above the celestial horizon, the time and azimuth of sunrise were calculated for each day of the year. The calculations were programmed in 64-bit Quick-Basic language (QB64; https://qb64.com/) and compiled as a stand-alone executable program file for faster processing times. The results of our simulations were checked against NOAA’s Solar Position Calculator (https://gml.noaa.gov/grad/solcalc/) to ensure that our results are robust and validate the model.