Astronomical Calculations: Solar Coordinates

In the last post I showed how to calculate the Julian Day (JD) and the Ephemeris Time (T). Now I want to build on that to calculate the geocentric (or apparent) coordinates of the Sun for any given moment in time. This algorithm is taken from Astronomical Algorithms by Jean Meeus (2nd Edition) with modifications from the U.S. Naval Observatory. All page number references are to Meeus’ book.

The goal here is to input any date and time (say the launch of Sputnik 1 on 4 Oct 1957 at 19:29 UTC) and get back the coordinates of the Sun as it appears to us on the celestial sphere at that moment. I recommend reading an introductory article Position of the Sun before continuing. For the complete software solution implemented in C# see my GitHub repo. Note that you will need to download and install the free Microsoft Visual Studio Community Edition 2017 to build and run the code. 

Continue reading “Astronomical Calculations: Solar Coordinates”

Astronomical Calculations: The Julian Day

Later I want to post an implementation of how to calculate geocentric solar coordinates for a given date and time based on algorithms published in Astronomical Algorithms by Jean Meeus. But that is very involved and going to be a long post. So I thought it would be wise first to talk about dates and times used by astronomers. This will “prepare the ground” so to speak.

In the West we use the Gregorian calendar which began on 4 Oct 1582 and replaced the Julian calendar. To fix the drift from the Julian calendar that had accumulated over the centuries,  4 Oct was followed by 15 Oct. As you can imagine this led to a lot of confusion since some people still used the old calendar and others referred to the new one. (Think of the Imperial and the metric system in our own time.) Also, not all cultures adopted the new calendar right away. An astronomer in Munich might still be using the Julian calendar while another in Venice used dates that referred to the new Gregorian calendar. So in 1583 a scholar named Joseph Scaliger proposed an abstraction which he called the Julian Period.  The period is 7,980 years and runs from 1 Jan 4713 BCE (which is Year 1) to 31 Dec 3268 CE.  We have a similar system today with POSIX (or Unix) time in which the epoch begins on 1 Jan 1970 UTC and counts the number of elapsed seconds since that moment.

From the Julian Period astronomers refer to the Julian Day (JD) when making calculations in celestial mechanics. A JD value is the number of days since the beginning of the Julian Period. Meeus writes:

The Julian Day number or, more simply, the Julian Day … is a continuous count of days and fractions thereof from the beginning of the year -4712.

Why -4712? Because dates always start at 0 so if 4713 BCE is year one then -4712 is the very beginning of the period for calculation purposes. Civil dates begin at midnight UTC. This is inconvenient for astronomers since celestial events might overlap from one day to the next. For example if you observe a comet at 23:00 hours and report on its movement until 02:00 hours you’d have to use two different dates. So the Julian Day begins at apparent solar noon UTC so that there’s a 12 hour difference from the beginning of the civil date.

So with all of that esoterica out of the way how do we calculate the JD? Here I’ll modify Meeus to show the pseudocode:

Let Y == Year
Let M == Month Number // 1 is Jan, 2 is Feb, etc.
Let D == Day Number // where D is double (64-bit floating-point value)

If M == 1 OR M == 2 Then:
    
    Y == Y - 1
    M == M + 12

// That is if the date is Jan or Feb then the date is the 13th or
// 14th month of the preceding year for calculation purposes.

If the date is on the Gregorian calendar (after 14 Oct 1582) then:

Let A = (int)(Y / 100)
Let B = 2 - A + (int)(A / 4)

If the date is on the Julian calendar (prior to 4 Oct 1582) then:

Let B = 0

So far we’ve set up the variables for the main calculation. D is a double because it could represent a fraction of a day. For example, 4.812 is the 4th day of the month at 19:29 (7:29 PM). Here’s the algorithm for the Julian Day itself:

JD = (int)(365.25 * (Y + 4716)) + (int)(30.6001 * (M + 1)) + D + B - 1524.5

I can implement this in C# as a struct:

/// <summary>
/// A moment is a specific point in time down to the millisecond.
/// </summary>
public struct Moment
{
    public int Year { get; }
    public int Month { get; }
    public int Day { get; }
    public int Hour { get; }
    public int Minute { get; }
    public int Second { get; }
    public int Millisecond { get; }

    /// <summary>
    /// Creates a Moment with known values down to the millisecond.
    /// </summary>
    public Moment(int y, int m, int d, int h, int m, int s = 0, int ms = 0)
    {
        Year = y;
        Month = m;
        Day = d;
        Hour = h;
        Minute = m;
        Second = s;
        Millisecond = ms;
    }

    /// <summary>
    /// A Julian Day (JD) is a continuous count of days and fractions thereof
    /// starting at 1 Jan -4712 at noon UTC to a given point in time thereafter.
    /// </summary>
    public double JulianDay
    {
        get
        {
            int Y = Year;
            int M = Month;
            int B = 0; // Julian calendar default

            // if the date is Jan or Feb then it is considered to be in the 
            // 13th or 14th month of the preceding year.
            switch (M)
            {
                case 1:
                case 2:
                    Y = Y - 1;
                    M = M + 12;
                    break;

                default:
                    break;
            }

            if (!IsJulianDate()) // convert to Gregorian calendar
            {
                var A = Y / 100;
                B = 2 - A + (A / 4);
            }

            return 
                (int)(365.25 * (Y + 4716)) + 
                (int)(30.6001 * (M + 1)) + DayOfMonth + B - 1524.5;
        }
    }


    /// <summary>
    /// Pope Gregory introduced the Gregorian calendar in October 1582 when the 
    /// calendar had drifted 10 days. Dates prior to 4 Oct 1582 are Julian dates
    /// and dates after 15 Oct 1582 are Gregorian dates. Any date in the gap is
    /// invalid on the Gregorian calendar.
    /// </summary>
    /// <returns></returns>
    public bool IsJulianDate()
    {
        if (Year > 1582)
            return false;

        if (Year < 1582)
            return true;

        // year is 1582 so check month
        if (Month > 10)
            return false;

        if (Month < 10)
            return true;

        // month is 10 so check days
        if (Day > 14)
            return false;

        return true;
    }

    public double DayOfMonth
    {
        get
        {
            return 
                Day + 
                (Hour / 24.0) + 
                (Minute / 1440.0) + 
                (Second + Millisecond / 1000.0) / 86400.0;
        }
    }
}

Now i can write a console app to test this implementation:

class Program
{
    static void Main(string[] args)
    {
        Console.WriteLine("Calculation of the Julian Day");
        Console.WriteLine(Environment.NewLine);

        var moment = new Moment(1957, 10, 4, 19, 29, 0);  // Sputnik 1 launched (UTC)

        Console.WriteLine("Sputnik 1 launched on 4 Oct 1957 at 19:29 UTC");
        Console.WriteLine("Sputnik 1 JD: " + moment.JulianDay);
        Console.ReadLine();
    }
}

And the output is:

Calculation of the Julian Day

Sputnik 1 launched on 4 Oct 1957 at 19:29 UTC
Sputnik 1 JD: 2436116.31180556

I will be using this struct in a future blog post when I dive headlong into calculating geocentric solar coordinates. That is to say, given any date and time (like the moment Sputnik 1 launched) I can calculate the position of the sun in the sky.