I'm using ArcGIS Pro to try to identify permanently shadowed craters on the North Pole of the Moon. Pictured below is a screenshot of the DEM I'm using which was collected from the laser altimeter on the Lunar Reconnaissance Orbiter. It is available for download from USGS here:
I can use the Hillshade tool with a low altitude angle and model_shadows turned on to create a raster like the one below that shows the shadowed regions of the craters for those particular illumination conditions. The picture below shows the results of a Hillshade with Azimuth 315 deg and Altitude 5 deg.
In a simplified world, I could just run this hillshade for a variety of Azimuth inputs and sum the resulting rasters. The pixels that still had a value of 0 would be permanently shadowed regions. There is one problem with this plan though. The moon is not flat, especially over such a large region. This DEM covers the whole surface of the Moon north of 75 deg latitude and even lower at the corners. Because of this, the elevation of the sun would vary quite a lot for various pixels in the image.
To summarize: Is there a way to calculate shadows caused by terrain on a large, spherical surface?
The Hillshade tool operates only on planar maps.
EDIT: Esri has a help page where they basically give you the hillshade formula and explain all the terms. I was able to modify this formula by calculating the unique azimuth and zenith angle for each pixel. The result below is almost what I'm looking for. In this image, the sun is off the top of the page and it is easy to see how the curvature of the moon makes the other side darker due to the higher zenith angles. The only thing this image is missing is shadows, which is why some of the crater walls on the bottom image are still illuminated despite being on the far side from the light. Unfortunately, I can't find any documentation on how shadows are calculated.