Solar angle calculations taking into account elevation of viewer, taken from Prayer Guardian

SunPosition.cs

using System;

namespace PrayerGuardian
{
    internal static partial class JulianDateExtensions
    {
        private static readonly DateTimeOffset JulianEpoch = new DateTimeOffset(new DateTime(1858, 11, 16, 12, 0, 0, DateTimeKind.Utc));
        public static double ToJulianDate(this DateTimeOffset date)
        {
            // Based off of https://en.wikipedia.org/wiki/Julian_day
            // 12h Nov 16, 1858 is 2400000JD
            // We can calculate days since then easily.

            var diff = date - JulianEpoch;
            return diff.TotalDays + 2400000;
        }
    }

    public class SunPosition
    {
        // Helper math functions that take and return degrees directly
        private static double DegreesToRadians(double degrees) => (Math.PI / 180) * degrees;
        private static double RadiansToDegrees(double radians) => radians * 180 / Math.PI;
        public static double sin(double d) => Math.Sin(DegreesToRadians(d));
        public static double asin(double v) => RadiansToDegrees(Math.Asin(v));
        public static double cos(double d) => Math.Cos(DegreesToRadians(d));
        public static double acos(double v) => RadiansToDegrees(Math.Acos(v));

        private DateTimeOffset _date;
        private double _latitude;
        private double _longitude;
        private double? _elevation;

        public SunPosition(DateTimeOffset date, double latitude, double longitude, double? elevationInMeters)
        {
            // Zero out the time portion to remove any offsets in calculations
            _date = new DateTimeOffset(new DateTime(date.Year, date.Month, date.Day, 0, 0, 0, DateTimeKind.Utc)); // .NET is stupid and doesn't account for `kind` in add/subtract operations
            _latitude = latitude;
            _longitude = longitude;
            _elevation = elevationInMeters.HasValue && !double.IsNaN(elevationInMeters.Value) ? elevationInMeters : null;
        }

        static double AdjustAngleForElevation(double alpha, double elevation)
        {
            return alpha - 2.076 * Math.Sqrt(elevation) / 60;
        }

        /// <summary>
        /// Converts a Julian date into a C# <c>DateTimeOffset</c> value.
        /// </summary>
        /// <param name="jd">The Julian date</param>
        /// <returns>The <c>DateTimeOffset</c> corresponding to <paramref name="jd"/>.</returns>
        private static DateTimeOffset FromJulianDate(double jd)
        {
            var epoch = new DateTimeOffset(new DateTime(1858, 11, 16, 12, 0, 0, DateTimeKind.Utc));
            return epoch.AddDays(jd - 2400000);
        }

        /// The azimuth and elevation angle (aka altitude) give a 3D location of a celestial object (their intersection in the 3rd plane)
        /// This image best illustrates it: https://en.wikipedia.org/wiki/File:Azimuth-Altitude_schematic.svg
        /// This is not an absolute measure, it is relative to the observer's location
        static AzEl CalculateAzEl(double omega, double declination, double latitude)
        {
            var elevation = AltitudeAtOmega(omega, declination, latitude);
            return new AzEl
            {
                Azimuth = CalculateAzimuth(elevation, omega, declination),
                Elevation = elevation
            };
        }

        /// <summary>
        /// Calculates the altitude (or elevation) of the sun above the horizon at a given hour angle.
        /// The altitude angle of the sun is an important measure in solar energy applications,
        /// and it represents the angle between the rays of the sun and the plane of the observer's horizon.
        /// </summary>
        /// <param name="omega">The hour angle, which is the time offset from solar noon, expressed in degrees.
        /// Positive values indicate times before solar noon, and negative values indicate times after solar noon.</param>
        /// <param name="declination">The solar declination, which is the angle between the rays of the sun
        /// and the plane of the Earth's equator, expressed in degrees. It varies with the time of year.</param>
        /// <param name="latitude">The latitude of the observer's location, expressed in degrees.
        /// Positive values are for the Northern Hemisphere, and negative values are for the Southern Hemisphere.</param>
        /// <returns>The altitude (aka alpha or elevation) of the sun in degrees. This is the angle between sun's rays
        /// and the plane of the observer's horizon at the specified time, location, and date.
        /// A positive angle means the sun is above the horizon, while a negative angle indicates that the sun is below the horizon.</returns>
        static double AltitudeAtOmega(double omega, double declination, double latitude)
        {
            return asin(sin(declination) * sin(latitude) + cos(declination) * cos(latitude) * cos(omega));
        }

        /// <summary>
        /// Calculates the azimuth angle in an astronomical context.
        /// </summary>
        /// <param name="altitude">The altitude of the celestial object, in degrees.</param>
        /// <param name="omega">The hour angle of the celestial object, in degrees.</param>
        /// <param name="declination">The declination of the celestial object, in degrees.</param>
        /// <returns>
        /// The azimuth of the celestial object, in degrees. The azimuth is typically measured clockwise from the north.
        /// </returns>
        /// <remarks>
        /// The azimuth is a horizontal angle measured at the observer's location.
        /// This function computes the azimuth based on the given altitude, hour angle, and declination.
        /// Note that this formula is a simplification and might not represent a complete calculation
        /// of the azimuth angle, as it only considers a part of what is typically required.
        /// </remarks>
        static double CalculateAzimuth(double altitude, double omega, double declination)
        {
            return sin(omega) * cos(declination) * cos(altitude);
        }

        /// <summary>
        /// Represents astronomical data regarding the sun's position and events.
        /// </summary>
        public class SunInfo
        {
            /// <summary>
            /// The angle of the sun at solar noon in degrees.
            /// Solar noon is the time of the day when the sun is at its highest point in the sky.
            /// </summary>
            public double AngleAtSolarNoon { get; init; }

            /// <summary>
            /// The declination of the sun in degrees.
            /// Declination is the angle between the rays of the sun and the plane of the Earth's equator.
            /// </summary>
            public required double Declination { get; init; }

            /// <summary>
            /// The time of sunrise.
            /// Sunrise is defined as the moment when the upper limb of the sun just appears above the horizon.
            /// </summary>
            public required DateTimeOffset Sunrise {get; init;}

            /// <summary>
            /// The time of sunset.
            /// Sunset is defined as the moment when the upper limb of the sun just disappears below the horizon.
            /// </summary>
            public required DateTimeOffset Sunset {get; init;}

            /// <summary>
            /// The time of solar noon.
            /// Solar noon is the time during the day when the sun crosses the meridian and is at its highest point above the horizon.
            /// </summary>
            public required DateTimeOffset SolarNoon {get; init;}

            /// <summary>
            /// A function that calculates the <c>DateTimeOffset</c> at which the sun will be at a specified altitude above or below the horizon.
            /// </summary>
            /// <param name="altitude">The altitude angle of the sun in degrees for which the time needs to be calculated.</param>
            /// <param name="angleReference">Determines whether calculations ar for the matching angle before or after solar noon</param>
            /// <returns>The <c>DateTimeOffset</c> at which the sun reaches the specified altitude.</returns>
            public required Func<double, AngleReference, DateTimeOffset> TimeAtSunAltitude {get; init;}

            /// <summary>
            /// A function that calculates the altitude of the sun, in degrees, at a given <c>DateTimeOffset</c>.
            /// </summary>
            /// <param name="time">The specific time for which the sun's altitude needs to be calculated.</param>
            /// <returns>The altitude angle of the sun in degrees at the specified time.</returns>
            public required Func<DateTimeOffset, double> AltitudeAtTime {get; init;}
        }

        /// <summary>
        /// Represents the solar position in the sky using azimuth and elevation angles.
        /// </summary>
        public struct AzEl
        {
            /// <summary>
            /// The azimuth angle of the sun.
            /// </summary>
            /// <value>
            /// The azimuth angle, in degrees, measured clockwise from the North.
            /// </value>
            public double Azimuth;

            /// <summary>
            /// The elevation angle of the sun.
            /// </summary>
            /// <value>
            /// The elevation angle, in degrees, above the horizon.
            /// A positive value indicates the sun is above the horizon, while a negative value indicates it is below.
            /// </value>
            public double Elevation;
        }

        /// <summary>
        /// Computes the hour angle at sunrise or sunset for a specified location,
        /// based on the given apparent altitude, solar declination, and observer's latitude.
        /// The hour angle is the angular distance of the sun from the observer's meridian at
        /// sunrise or sunset and it helps determine the exact time of these events.
        /// </summary>
        /// <param name="alpha">
        /// The apparent altitude for sunrise or sunset. This value accounts for adjustments like atmospheric refraction
        /// and the solar disc's elevation.
        /// </param>
        /// <param name="declination">
        /// The solar declination, which is the angular distance of the sun north or south of the celestial equator.
        /// </param>
        /// <param name="latitude">
        /// The observer's latitude on the Earth's surface.
        /// </param>
        /// <returns>
        /// The hour angle at which sunrise (or sunset) occurs, expressed in degrees.
        /// </returns>
        static double SunriseEquation(double alpha, double declination, double latitude)
        {
            /// <summary>
            /// Normalizes an astronomical angle to a range between -1 and 1.
            /// </summary>
            double normalize(double l)
            {
                if (l > 0)
                {
                    return (l + 1) % 2 - 1;
                }
                return 1 - Math.Abs((l - 1) % 2);
            }

            // This accounts for atmospheric refraction and the non-zero angle subtended by the solar disc
            // elevation is accounted for by adjusting alpha separately
            var step1 = (sin(alpha) - sin(latitude) * sin(declination)) / (cos(latitude) * cos(declination));
            var normalized = normalize(step1);
            //var step2 = acos(normalized);
            var step2 = acos(step1);

            return step2;
        }

        public enum AngleReference
        {
            East,
            West,
            BeforeSolarNoon = East,
            AfterSolarNoon = West
        }

        static public SunInfo Calculate(DateTimeOffset date, double latitude, double longitude, double? elevation)
        {
            return new SunPosition(date, latitude, longitude, elevation).Calculate();
        }

        /// <summary>
        /// Calculates information about the location of the sun for the date, latitude, longitude, and elevation
        /// that the class was initialized with.
        /// </summary>
        public SunInfo Calculate()
        {
            // the date to calculate, expressed as a Julian date
            var jd = _date.ToJulianDate();

            // days since Jan 1 2000 12:00
            var n = (_date - new DateTime(2000, 1, 1, 12, 0, 0, DateTimeKind.Utc)).TotalDays;
            var J = n - _longitude / 360; // mean solar noon
            var M = (357.5291 + 0.98560028 * J) % 360; // mean solar anomaly
            var C = 1.9148 * sin(M) + 0.200 * sin(2 * M) + 0.0003 * sin(3 * M); // equation of the center
            var argumentOfPerihelion = 102.9372; // an approximation...
            {
                // ...further refined
                // https://ssd.jpl.nasa.gov/txt/aprx_pos_planets.pdf
                var J2000 = new DateTimeOffset(new DateTime(2000, 1, 1, 12, 0, 0, DateTimeKind.Utc)).ToJulianDate();
                var daysInCentury = 36525; // according to JPL: 36524.22 is more accurate
                var centuriesSinceJ2000 = (_date.ToJulianDate() - J2000) / daysInCentury;

                var degreesAtJ2000 = 102.9300589;
                var degreesPerCentury = 0.3179526; // valid until 3000AD
                argumentOfPerihelion = degreesAtJ2000 + degreesPerCentury * centuriesSinceJ2000;
            }

            var lambda = (M + C + 180 + argumentOfPerihelion) % 360; // ecliptic longitude

            var jTransit = 2451545.5 + J + 0.0053 * sin(M) - 0.0069 * sin(2 * lambda); // solar noon

            var declination = asin(sin(lambda) * sin(23.44)); // the angle between the rays of the sun and the plane of the earth's equator

            // Sunrise/sunset is defined as the time when the sun is at ±0.83 degrees
            // However, the apparent solar altitude angle differs from the actual depending on elevation. If we know the elevation, we can adjust for that.
            var alpha = -50 / 60.0;
            var angle = _elevation is null ? alpha : AdjustAngleForElevation(alpha, _elevation.Value); //
            var omegaS = SunriseEquation(angle, declination, _latitude); // sunset and sunrise angle

            var jSunrise = jTransit - omegaS / 360;
            var jSunset = jTransit + omegaS / 360;

            // Calculate sun altitude angle at solar noon.
            // While one would expect this to be zero, that is not the case for many reasons (observer not at equator so sun will not be directly overhead, elliptical earth, axis tilt, etc)
            var azEl = CalculateAzEl(0, declination, _latitude); // omega (the hour angle) is zero at solar noon.

            /// <summary>
            /// Calculates the time of day when the Sun reaches a specified altitude angle from the observer's horizon.
            /// </summary>
            /// <param name="alpha">The altitude angle of the Sun in degrees that is desired to be calculated.</param>
            /// <param name="reference">Specifies the angular reference direction (East or West) for calculation. Determines whether the calculation is for sunrise or sunset.</param>
            /// <returns>
            /// A <see cref="DateTimeOffset"/> that represents the time of day (UTC) when the Sun reaches the specified altitude angle.
            /// </returns>
            /// <remarks>
            /// The function uses the sunrise equation to determine the time when the Sun reaches a given altitude. If the direct calculation using the sunrise equation is not applicable, a numerical method (Newton's Method) is used for more general cases. The calculation depends on the observer's geographical latitude and the Sun's declination on the given day.
            /// </remarks>
            /// <exception cref="ArgumentOutOfRangeException">Thrown when the altitude angle is outside the feasible range for sunrise or sunset calculations.</exception>
            /// <exception cref="InvalidOperationException">Thrown when the calculations fail to converge to a solution using the Newton Method.</exception>
            DateTimeOffset TimeAtSunAltitude(double alpha, AngleReference reference)
            {
                // alpha = -1 * Math.Abs(alpha); // regardless of sign, angle is always same
                // alpha = AdjustAngleForElevation(alpha);

                // This is what we used to use..
                // However, the sunrise equation seems to only be usable within a small subset of values
                var omega = SunriseEquation(alpha, declination, _latitude);
                // System.Diagnostics.Debug.Assert(!double.IsNaN(omega));
                if (!double.IsNaN(omega))
                {
                    return reference switch
                    {
                        AngleReference.East => FromJulianDate(jTransit - omega / 360),
                        AngleReference.West => FromJulianDate(jTransit + omega / 360)
                    };
                }

                (DateTimeOffset start, DateTimeOffset end) = reference switch
                {
                    AngleReference.East => (FromJulianDate(jSunrise), FromJulianDate(jTransit)),
                    AngleReference.West => (FromJulianDate(jTransit), FromJulianDate(jSunset)),
                };

                return new NewtonMethod<DateTimeOffset, double>(start, end, Math.Abs(alpha), (t1, t2) => new DateTimeOffset(new DateTime(t1.Ticks + (t2.Ticks - t1.Ticks) / 2, DateTimeKind.Utc)), t => AltitudeAtTime(t), 0.001).Result;
            }

            /// <summary>
            /// Calculates the altitude (or elevation) of the sun above the horizon at a given <c>DateTimeOffset</c>.
            /// </summary>
            double AltitudeAtTime(DateTimeOffset time)
            {
                time = time.ToUniversalTime();
                var solarNoon = jTransit;
                var diff = Math.Abs(time.ToJulianDate() - solarNoon); // in days
                var omega2 = diff * 360; // a day is 360 degrees, so days to degrees is to multiply by 360
                return AltitudeAtOmega(omega2, declination, _latitude);
            }

            var sunInfo = new SunInfo()
            {
                AngleAtSolarNoon = azEl.Elevation,
                Declination = declination,
                Sunrise = FromJulianDate(jSunrise),
                SolarNoon = FromJulianDate(jTransit),
                Sunset = FromJulianDate(jSunset),
                TimeAtSunAltitude = TimeAtSunAltitude,
                AltitudeAtTime = AltitudeAtTime
            };

            return sunInfo;
        }
    }
}