Once you have used pyephem or an R script to calculate the rise/set azimuths of stars, the Sun and the Moon at a particular date for a particular latitude, you need to correct the calculation for horizon features.

In order to determine the horizon features, you need to obtain the topographic information for the surrounding region. The National Geographic Data Center of the National Oceanic and Atmospheric Administration has a webpage where you can download topographic information in 1′ grids of latitude and longitude.  To determine how big a grid you need to download you can use the rough rule of thumb that at most you would likely be able to see ~100km from a very tall mountain on a very clear day.  You can use this online calculator to determine the range in latitude and longitude from your central position in order to get a distance of 100km.  Note that the range in latitude needed to achieve this is pretty constant (about +/- 1 degree), but the range in longitude needed depends on the latitude because the longitude points get closer together at the poles.  At the latitude of the Merry Maidens stone circle site in Cornwall (50N) a range of +/-1.5 degrees of longitude produces a range of ~+/-100km.

On the NGDC website choose the lon,lat,depth format with 4-byte floating precision (under advanced options), and space delimited format.