I actually very briefly looked into how get this basic data about sunrise/sunset/moonrise/moonset times and headings for any location on Earth (because I wanted to make a simple website to visualize the data as a fun side project). There must be some free/open source software models of the Sun-Earth-Moon system that could be combined with some basic GIS model of Earth to provide reasonably accurate data for any location on Earth, but I'm just so far away from these fields I probably don't know the right things to search for.
It basically depends on how accurate you want to be. I looked into this so I could generate lunar calendars for friends. I got the phases of the moon predictions from the BSD Games program pom by Keith E. Brandt[0], which is itself based on 'Practical Astronomy with Your Calculator'[1] by Peter Duffett-Smith. Lunar eclipse predictions came from the Basic Lunar Eclipse and Phase Calculator by Xavier Jubier[2]. These work well for the near future and with a reasonable amount of precision for regular people. You can get much more accurate using Astronomy Engine[3], but that was already a lot more complex than I wanted to get into
I'd agree about accuracy but it is also worth mentioning the time reach. Plus or minus half a century from now (or J2000.0 epoch) is fine but if you wanted to explore historical eclipses then you will need to look at higher precision series. This whole thing is about truncating trigonometric series with amplitudes and angles that change over time.
For sunrise/set and moonrise/set calculations you also need to allow for refraction in the Earth's atmosphere (around half a degree or so). Also for moonrise/set you need to allow for the parallax of the Moon's position because you are standing on the surface of the Earth not the centre (up to 1 degree difference in position). This all assumes a spherical Earth and an Ocean horizon. And the rise/set time can changes a little bit with your altitude above the horizon.
For latitudes near the Arctic/Antarctic circles any method based on iteration (as in Duffett-Smith's lovely book) is going to have convergence issues because there might not actually be a rise/set event near that day so if you really want an algorithm that will work anywhere on Earth you'll need to use a method based on fitting a quadratic curve to (say) 24 positions calculated for the day. Have a look at [1] in a library. The code is in Pascal, but it is the search method using a quadratic fit to sets of 3 consecutive positions you want.
Just for low-precision lunar positions, an interesting alternative that also summarises the various methods was produced by NASA [2] for space craft use.
Last time I looked into this I implemented some stuff from Astronomical Algorithms https://www.amazon.com/dp/0943396352 which was more fun for me than just grabbing the data.
