Realistic sun position based on time of year and latitude/longitude in Unity 3D. Making it possible to simulate phenomena like the midnight sun. If this came in handy, say hi over at @SnutiHQ on Twitter.

SunPosition.cs

// This code is licensed under the terms of the MIT license
using UnityEngine;

public class SunPosition : MonoBehaviour
{
    [SerializeField, Range(-90f, 90f)]
    public float latitude = 69.6546f;
    [SerializeField, Range(-180f, 180f)]
    public float longitude = 18.9637f;
    
    [SerializeField, Range(0, 365)]
    public int dayOfYear = 180;
    [SerializeField, Range(0f, 24f)]
    public float hourOfDay = 22f;

    private float zenith = 41f;
    private float azimuth = 160f;
    private float solarHourAngle = 13f;

    private void OnValidate()
    {
        UpdateSun();
    }

    public void UpdateSun()
    {
        Calculate();
        SetRotation();
    }
    
    private void Calculate()
    {
        // Time of year
        var time = Mathf.Deg2Rad * (360f / 365.25f) * (90f - dayOfYear  + hourOfDay / 24f);

        // Solar Declination
        var declination = 0.396372f - 22.91327f * Mathf.Sin(time) + 4.02543f * Mathf.Sin(time) 
            - 0.387205f * Mathf.Cos(2f * time) + 0.051967f * Mathf.Sin(2f * time) 
            - 0.154527f * Mathf.Cos(3f * time) + 0.084798f * Mathf.Sin(3f * time);

        // Time correction for solar angle
        var timeCorrection = 0.004297f + 0.107029f * Mathf.Cos(time) 
                               - 1.837877f * Mathf.Sin(time) 
                               - 0.837378f * Mathf.Cos(2f * time) 
                               - 2.340475f * Mathf.Sin(2f * time);

        // Solar Hour Angle
        solarHourAngle = (hourOfDay - 12f) * 15f + longitude + timeCorrection;
        if (solarHourAngle > 180f) solarHourAngle -= 360f;
        if (solarHourAngle < -180f) solarHourAngle += 360f;

        // Sun Zenith Angle
        var cosZenith = Mathf.Sin(Mathf.Deg2Rad * latitude) * Mathf.Sin(Mathf.Deg2Rad * declination) +
                          Mathf.Cos(Mathf.Deg2Rad * latitude) * Mathf.Cos(Mathf.Deg2Rad * declination) *
                          Mathf.Cos(Mathf.Deg2Rad * solarHourAngle);
        zenith = Mathf.Rad2Deg * Mathf.Acos(cosZenith);

        // Azimuth Angle
        var cosAzimuth = (Mathf.Sin(Mathf.Deg2Rad * declination) -
                            Mathf.Sin(Mathf.Deg2Rad * latitude) * Mathf.Cos(Mathf.Deg2Rad * zenith))
                           / (Mathf.Cos(Mathf.Deg2Rad * latitude) * Mathf.Sin(Mathf.Deg2Rad * zenith));
        azimuth = Mathf.Rad2Deg * Mathf.Acos(cosAzimuth);
    }

    private void SetRotation()
    {
        var x = 90f - zenith;
        var y = 180f + (solarHourAngle > 0f ? -azimuth : azimuth); // Z+ is north
        transform.eulerAngles = new Vector3(x, y, 0f);
    }
}