Here is a little programme I wrote to calculate the date and time of equinoxes and solstices in a given year (not earlier than 4713 BC). Actually I am no astronomer: I have only collected information on calendars and astronomical calculations from various sources on the Internet and assembled what I needed.
Years BC are given as negative numbers. Note that there is no year 0 BC; the day before 1 Jan 1 AD is 31 Dec 1 BC. Note also that 1 BC, 5 BC, etc. are leap years.
Pope Gregory XIII decreed that the Julian calendar would end on 4 Oct 1582 AD and that the next day would be 15 Oct 1582 in the Gregorian Calendar. The only other change is that centesimal years (years ending in 00) would no longer be leap years unless divisible by 400. Britain and its possessions and colonies continued to use the Julian calendar up until 2 Sep 1752, when the next day became 14 Sep 1752 in the Gregorian Calendar. This routine can be compiled to use either convention. By default, as I am Italian and Roman Catholic, Pope Gregory's convention is used. Simply comment out
#define PAPALto use British convention.
This programme uses formulae taken from Jean Meeus's "Astronomical Algorithms" (1991). CalculateDate() function is based on formulae originally posted by Tom Van Flandern Washington, DC email@example.com in the UseNet newsgroup sci.astro. Reposted 14 May 1991 in FidoNet C Echo conference by Paul Schlyter (Stockholm). Minor corrections and additional information on JDNs and calendars by Raymond Gardner Englewood, Colorado.