What Sidereal Time Is It?

It was one of those rare clear nights we sometimes get this time of year in the Pacific Northwest. My wife had been outside looking up at the night sky. Soon she came inside the house and asked, “what’s the sidereal time right now?” It was good timing because I had just added a sidereal time object to the SquareWidget.Astronomy.Core code library. While she looked over my shoulder, I started to write some code to answer the question. I looked up our local longitude (L) and in four lines of code I had the answer:

using SquareWidget.Astronomy.Core.UnitsOfMeasure;

Moment moment = new(DateTime.UtcNow);
SiderealTime st = new(moment);
SexigesimalAngle L = new(-123, 15, 43.34);
Console.WriteLine(st.ToLocalMean(L));

It returned 7 hours and some odd arcminutes and arcseconds. She nodded approvingly because she already knew the answer having been outside and saw that the constellation Orion had recently crossed our local meridian.

It’s been about four weeks since I released the first version of the code library. It took a lot out of me, so my plan was to chill out and relax for a good long while. And I did. For about 24 hours. I kept identifying and making a list of new features that I wanted to add at some point. Soon I found myself peeling them off one at a time and building releases around each one. So now the library is up to v1.5.1 and I think it’s time to do “Chill Out 2.0” for real this time. But first, let me walk through some of the new features to illustrate them. See my previous blog post on the subject as well as the README documentation on GitHub for more details.

Solar Coordinates

The first version had a solar longitude calculator. The calculator has been deprecated and the functionality is now in the Sun object. The coordinates are returned as apparent (geocentric) equatorial coordinates at a given moment. Suppose I want to know the Sun’s coordinates on 4 Jul 2028:

using SquareWidget.Astronomy.Core.UnitsOfMeasure;
using SquareWidget.Astronomy.Core.CelestialObjects.Stars;

DateTime d = new(2028, 7, 4);
Moment moment = new(d);
Sun sun = new(moment);

EquatorialCoordinates eqc = sun.GetGeocentricPosition();

Console.WriteLine($"RA {eqc.α.ToString()}");
Console.WriteLine($"Dec {eqc.δ.ToString()}.");

/* Displays:
 
RA 6h 54m 36.194953s
Dec +22° 50' 36.22891852".

*/

For a deep dive into the algorithm on solar coordinates see my earlier blog post Astronomical Calculations: Solar Coordinates.

Moon Phases

There is now a calculator to return moon phases within a given date range. Some phases fall outside the range in order to get a baseline for the previous or next phase:

using SquareWidget.Astronomy.Core.Models;
using SquareWidget.Astronomy.Core.Calculators;
using SquareWidget.Astronomy.Core.CelestialObjects.Moons;

DateOnly startDate = new(2025, 1, 1);
DateOnly endDate = new(2025, 3, 31);
DateRange dateRange = new(startDate, endDate);

List<MoonPhase> list = MoonPhaseDatesCalculator.Calculate(dateRange);

foreach (var item in list)
{
    Console.WriteLine(item.PhaseName + " on " +
        item.Moment
            .ToDateTime()
            .ToString("MM/dd/yyyy hh:mm:ss.fff tt"));
}



/* Displays:
 
NewMoon on 12/30/2024 10:26:46.440 PM
FirstQuarter on 01/06/2025 11:56:49.386 PM
LastQuarter on 01/21/2025 08:30:00.587 PM
FullMoon on 02/12/2025 01:53:20.083 PM
NewMoon on 01/29/2025 12:35:53.996 PM
FirstQuarter on 02/05/2025 08:01:06.996 AM
LastQuarter on 02/20/2025 05:31:34.392 PM
FullMoon on 03/14/2025 06:54:32.963 AM
NewMoon on 02/28/2025 12:44:39.649 AM
FirstQuarter on 03/06/2025 04:30:40.312 PM
LastQuarter on 03/22/2025 11:29:33.015 AM
FullMoon on 04/13/2025 12:22:13.127 AM

*/

Solar Eclipses

As with moon phases above you build a date range and pass that into the calculator to get a list of solar eclipses within the range. Here I’m entering the current year 2024 and I see the total eclipse this April that will swing through the U.S. Midwest.

using SquareWidget.Astronomy.Core.Calculators;
using SquareWidget.Astronomy.Core.Models;

DateOnly startDate = new(2024, 1, 1);
DateOnly endDate = new(2024, 12, 31);
DateRange dateRange = new(startDate, endDate);

IEnumerable<SolarEclipse> eclipses = SolarEclipseCalculator.Calculate(dateRange);

foreach (var eclipse in eclipses)
{
    Console.WriteLine(eclipse.ToString());
}

/* Displays:
 
04-08-2024 18:18 UTC  g:  0.3438   Total
10-02-2024 18:46 UTC  g: -0.3515   Annular

*/

Check out NASA’s Five Millennium Canon of Solar Eclipses -1999 to +3000 to see a visual diagram that you can compare with the results.

Galilean Moons of Jupiter

A few years ago, I showed how to calculate, for a given date, the apparent rectangular coordinates (X, Y) of the four Galilean moons of Jupiter (here). You can see a visual representation of the data — a sort of corkscrew diagram — that accompanies that article here. This past month, I put the feature into the code library. You use the same date range pattern as with the other calculators above:

using SquareWidget.Astronomy.Core.Calculators;
using SquareWidget.Astronomy.Core.Models;

DateTime startDate = DateTime.Today;
DateTime endDate = startDate.AddDays(1);

for (var d = startDate; d <= endDate; d = d.AddDays(1))
{
    SatellitePositions positions = JupiterSatellitePositionCalculator.Calculate(d);
    Console.WriteLine(positions.ToString());
}

/* Displays:
 
2/24/2024 12:01 AM
      Io - X:  -4.06   Y:   0.22
  Europa - X:  -0.58   Y:   0.48
Ganymede - X: -14.48   Y:   0.20
Callisto - X:  10.73   Y:  -1.24

2/25/2024 12:01 AM
      Io - X:   5.38   Y:  -0.12
  Europa - X:  -9.02   Y:  -0.12
Ganymede - X: -12.32   Y:  -0.44
Callisto - X:  18.74   Y:  -0.95

*/

The SatellitePositions object holds the date, all four moons, and their respective X, Y coordinates with respect to Jupiter.

Horizontal Coordinates

I’ve added horizontal coordinates as a third system for Alt-Az support. See the code sample for a detailed implementation. If you already have horizontal coordinates, you can convert them to equatorial:

// convert horizontal coordinates to equatorial coordinates
EquatorialCoordinates eqc = hc.ToΕquatorialCoordinates(moment, φ, L);

You need a Moment instance, and the observer’s latitude (φ) and longitude (L) for the conversion.

Sidereal Time

The code library uses sidereal time for conversion between coordinates. You can use it to find the Greenwich Mean Sidereal Time (GMST) or the Greenwich Apparent Sidereal Time (GAST) for any given moment in time. Always use UTC when initializing a SiderealTime object:

using SquareWidget.Astronomy.Core.Planets;
using SquareWidget.Astronomy.Core.UnitsOfMeasure;

DateTime datetime = new(1987, 4, 10, 19, 21, 0);
Moment moment = new Moment(datetime);
SiderealTime st = new(moment);

RightAscension gmst = new(st.GreenwichMean);
RightAscension gast = new(st.GreenwichApparent);

Console.WriteLine(gmst.ToString());
Console.WriteLine(gast.ToString());


/* Displays:
 
8h 34m 57.089579s
8h 34m 57.073455s

*/

If the observer’s longitude (L) is known, GMST can be converted to LMST:

using SquareWidget.Astronomy.Core.Planets;
using SquareWidget.Astronomy.Core.UnitsOfMeasure;

DateTime datetime = new(1987, 4, 10, 19, 21, 0);
Moment moment = new Moment(datetime);
SiderealTime st = new(moment);

RightAscension gmst = new(st.GreenwichMean);

// longitude in Corvallis, OR
SexigesimalAngle L = new(-123, 15, 43.34); 

// local mean sidereal time on the date (PST -8 offset)
RightAscension lmst = st.ToLocalMean(L);   

Console.WriteLine(gmst.ToString());
Console.WriteLine(lmst.ToString());


/* Displays:
 
8h 34m 57.089579s
0h 21m 54.200245s

*/

Finally, the hour angle (H) for an observer’s location and an object’s equatorial right ascension (α) is available by calling ToHourAngle:

using SquareWidget.Astronomy.Core.Planets;
using SquareWidget.Astronomy.Core.UnitsOfMeasure;

DateTime datetime = new(2024, 7, 4, 17, 0, 0); // UTC
Moment moment = new(datetime);
SiderealTime st = new(moment);

// longitude in Los Angeles, USA
SexigesimalAngle L = new(-118, 14, 38); 

// An object's equatorial right ascension
RightAscension α = new(5.838);

Degrees H = st.ToHourAngle(L, α);

Console.WriteLine(H.ToString());

/* Displays:
 
290°.6014494764589

*/

I’d be lying if I said I wasn’t thinking of new features and yes, I’ve started another list of enhancements. But for now, I’m going to take a break from this project for a bit and recharge my batteries. When I pick it up again it will be to do more testing. Currently, there are about 60 unit tests but there’s always room for more coverage and more tests. I also want to focus on useability. Then over time I’ll look to add more enhancements: equinoxes, solstices, equation of time, transits, conjunctions, perihelion and aphelion… the list goes on!

 

 

 

 

 

 

 

 

 

SquareWidget.Astronomy.Core: A C# Code Library

I’m happy to announce the release of SquareWidget.Astronomy.Core, a new code library targeting .NET 8.0 and written in C#, that supports common astronomical calculations and algorithms. The current version is 1.1.0. You can reference the NuGet package at the link above. It has a baseline of functionality and I have plans to add more over time.

Why a library? I’ve been writing code (and snippets of code) over several years to solve various problems in astronomy. When something looked interesting enough to share, I published some of that code on my GitHub repo and on this blog. Recently I was ploughing through the position angle of Saturn’s north pole (Jean Meeus,  Astronomical Algorithms, Chapter 49). After some long debugging sessions, I realized I was using the double data type for everything. Several times I would forget whether a value was in degrees, radians, large angles, reduced angles, or something else entirely. So, I decided I needed a basic code library to bring a little order out of the chaos. The first thing I added were the units of measure. Then planets and calculators came later. I have two design goals in mind for the library: useability and robustness. Certainly, the library should return results accurate to within a few minutes of degrees. But it should also be intuitive and easy to use.

I’ve written more detailed documentation on the GitHub README. Here I’ll just introduce a few of the features. At the heart of the library is a Moment struct. You can pass in date and time parameters, a DateTime instance, or a Julian Day (JD) value:

Moment m1 = new(2024, 1, 31, 9, 36, 0);
Moment m2 = new(2024, 1, 31.4);
Moment m3 = new(2460340.90001);

string s1 = m1.ToDateTime().ToString("yyyy MMM dd HH:mm");
string s2 = m2.ToDateTime().ToString("yyyy MMM dd HH:mm");
string s3 = m3.ToDateTime().ToString("yyyy MMM dd HH:mm");

Console.WriteLine(s1); // 2024 Jan 31 09:36
Console.WriteLine(s2); // 2024 Jan 31 09:36
Console.WriteLine(s3); // 2024 Jan 31 09:36

Once you have a Moment you can get its JD, modified JDE, Day D, or Time T (measure of Julian centuries from epoch J2000.0).

The code library also includes several units of measure with the ability to convert between them. For example, suppose I want to model the Earth’s obliquity of the ecliptic.  I can create a Degrees object by passing in decimal degrees:

Degrees d = new(23.41);
Console.WriteLine(d.ToString()); // 23°.41

Units can be converted to other units. Here I’ll convert to Radians, then a SexesimalAngle, and back to Degrees again:

Degrees d = new(23.41);
Console.WriteLine(d.ToString()); // 23°.41

Radians r = d.ToRadians();
Console.WriteLine(r.ToString()); // 0.4085816

SexigesimalAngle sa = new(d);
Console.WriteLine(sa.ToString()); // +23° 24' 36"

Degrees d1 = new(sa);
Console.WriteLine(d1.ToString()); // 23°.41

Quite often in astronomical calculations we deal with very large angles. Those can be reduced to within the range [0, 360] with the ToReducedAngle function:

Degrees d = new(13512.45);
Console.WriteLine(d.ToString()); // 13512°.45
Console.WriteLine(d.ToReducedAngle().ToString()); // 192°.45

A RightAscension can be constructed with the standard hours, minutes, and seconds (HMS) or decimal degrees:

// passing in HMS
RightAscension ra = new(12, 37, 27);
Console.WriteLine(ra.ToString());   // 12h 37m 27s

// passing decimal degrees
RightAscension ra = new(189.3625);
Console.WriteLine(ra.ToString()); // 12h 37m 27s

Two coordinate systems have been implemented, equitorial and ecliptical,  and you can convert back and forth between them. Suppose we have the equatorial coordinates (α, δ) of Pollux from a star catalog or the wiki and we know Earth’s obliquity of the ecliptic (ε):

// Pollux (β Gem)
RightAscension α = new(7, 45, 18.946);
SexigesimalAngle δ = new(28, 01, 34.26);
Degrees ε = new(23.4392911);

EquitorialCoordinates eqc = new(δ, α, ε);
Console.WriteLine(eqc.α.ToString()); // 7h 45m 18.946s
Console.WriteLine(eqc.δ.ToString()); // +28° 1' 34.26"

We can take those coordinates and convert them to ecliptical coordinates:

EquitorialCoordinates eqc = new(δ, α, ε);
EclipticalCoordinates ec = eqc.ToΕclipticCoordinates();

Console.WriteLine(ec.λ.ToDegrees().ToReducedAngle().ToString()); // 113°.21
Console.WriteLine(ec.β.ToDegrees().ToString());                  //   6°.68

You can instantiate a planet with the new operator or use the factory method. You must pass in a Moment struct either way.

DateTime datetime = new(2028, 7, 1);
Moment moment = new(datetime);

// using new operator
Earth earth = new(moment);

// using factory method for polymorphism
Planet planet = PlanetFactory.Create(PlanetName.Earth, moment);

Once you have a planet you can get its properties which include orbital elements and heliocentric spherical coordinates.

The library has six calculators right now and I plan to add at least a couple more in the next few weeks.

    • Geocentric position calculator
    • Moon phase and dates calculator
    • Nutation calculator
    • Position Angle of Saturn’s semiminor axis calculator
    • Sundial calculator
    • Solar longitude calculator

Consult the README documentation and the unit tests for details. Let me just illustrate one here. Suppose I want to know the four main phases of the moon for the month of August 2024:

DateOnly startDate = new(2024, 7, 1);
DateOnly endDate = new(2024, 9, 30);
DateRange dateRange = new(startDate, endDate);

var list = MoonPhaseDatesCalculator.Calculate(dateRange);

// use a collection expression to get just August
list = [.. list
    .Where(i => i.Moment.Month == 8)
    .OrderBy(o => o.Moment.ToDateTime())];

foreach (var item in list)
{
    Console.WriteLine(item.PhaseName + " on " + 
        item.Moment.ToDateTime().ToString("MM/dd/yyyy hh:mm:ss.fff tt"));
}

This displays:

NewMoon on 08/04/2024 11:12:56.525 AM
FirstQuarter on 08/12/2024 03:19:31.125 PM
FullMoon on 08/19/2024 06:25:37.281 PM
LastQuarter on 08/26/2024 09:27:16.726 AM

What’s next? I’ve got two enhancements open right now. One is to add a feature to calculate solar eclipses and the other is to calculate the positions of the four Galilean moons of Jupiter. Sometime after that I’ll want to address magnitude, planetary perihelion and aphelion, and probably more fun stuff with the Moon. Drop me a line if you have any requests.

Calculate Future Solar Eclipses

According to Herodotus, who wrote The History around 430 BCE, a total solar eclipse interrupted a battle between the Lydians and the Medes:

[J]ust as the battle was growing warm, day was on a sudden changed into night. This event had been foretold by Thales, the Milesian, who forewarned the Ionians of it, fixing for it the very year in which it actually took place. The Medes and Lydians, when they observed the change, ceased fighting, and were alike anxious to have terms of peace agreed on.

(Source: Internet Classic Archives)

This “battle of the eclipse” took place in what is now Turkey, and the eclipse itself could have been the one that occurred on 28 May 585 BCE. Unfortunately, if Thales wrote anything it didn’t survive the centuries. So, we don’t know if he really did predict the eclipse or by what means he did it. All we have is Herodotus to go on. In any case, it’s the first mention in history of anyone predicting an eclipse so it’s interesting on that basis alone.

Wouldn’t it be cool if we could figure out when future solar eclipses occurred for us in the next decade or so? That’s what we’re going to do in this blog post so read on. Continue reading “Calculate Future Solar Eclipses”

Let’s Make a Sundial!

Let’s make a sundial that tells the hour of the day. I know what you’re thinking. Can’t I just go to my local garden store and buy one? Nope. Those are just for decoration only. A real sundial, one that accurately tells the time, must be built for your location (latitude). All sundials are local.

I recommend reading my previous post on building a solar calendar. Also, if you’ve never given much thought to the apparent motion of the sun across the sky, then check out my post on sidereal time. If you know all about sundials and you just want to build one, then read on!

Continue reading “Let’s Make a Sundial!”

Astronomical Calculations: Sidereal Time

What is Sidereal Time?

The vernal equinox for 2020 was last week (as of this writing) and I had planned on posting this on that day. But there’s a global pandemic spreading through the world right now and, like so many others in North America, I was preoccupied with other things. In any case, in honor of the vernal equinox this post is going to talk about sidereal time. If you’re an amateur astronomer interested in the algorithm for local sidereal time then skip down to the next section. If you’re coming in cold then read on!

Continue reading “Astronomical Calculations: Sidereal Time”

Astronomical Calculations: Moons of Jupiter

Introduction

In this post I’m going to walk through an algorithm for calculating the positions of the four Galilean moons of Jupiter for any given moment in time. This method was developed by Jean Meeus and published in his book Astronomical Algorithms (Willmann-Bell, 2nd Ed., 1998).

All source code is checked into my GitHub repo. You can see the working software in action here: https://moonsofjupiter.azurewebsites.net. If you’re a raw beginner and want to know what you’ve just stumbled upon before diving into this post here are some good places to start:

(1) Ernie Wright provides an introduction with some cool animations that show what Galileo saw at the eyepiece and what he sketched in his own notebook.

(2) Dominic Ford has all kinds of astronomy goodness, including a nice looking corkscrew diagram for the four Galilean moons.

(3) I published an example web site  to  accompany this post. It would be a good idea to play around with it a bit before proceeding so you get an idea of what we’re trying to do.

Continue reading “Astronomical Calculations: Moons of Jupiter”

Astronomical Calculations: Delta T

Historically, Greenwich Mean Time (GMT) — solar noon at the Royal Observatory in Greenwich, England — was used as a universal reference point for trains and ships across the British empire. It is a timekeeping system based on the Earth’s rotation.  And it is a “mean” (average) because the exact moment the Sun crosses the Observatory’s meridian varies throughout the year. (This is a function of the Earth’s elliptical orbit and  angular velocity that I won’t go into here; see the equation of time for a full explanation.)

Eventually we needed something more  universal than a time standard kept in sync with one spot on Earth. Enter the atomic clock. An atomic clock can keep an exact 24 hours, day after day, century after century regardless of when the Sun crosses the meridian of any spot on Earth. Using cesium-133 atoms as a reference point, an atomic clock has an error deviation of about one second in 1.4 million years! We refer to this time as International Atomic Time (TAI).

We now have a solid, universal fixed-length period. Is that it? Can we now just use TAI for our civil time? No because the Earth and the natural forces that affect it are not precise. Things speed up and slow down due to tidal friction and other smaller perturbations. In a nutshell, the Earth cannot keep in sync with the precision of an atomic clock. It’s rotation is slowing down. It’s just like a toy top that a child spins up. The top spins very fast in the first few seconds and then over time loses its angular momentum and gradually spins slower. Our Earth is like a top that is slowing down. But not to worry. Long before it wobbles to a stop it will have been swallowed up by an expanding Sun!

So on the one hand we have a very precise TAI. And on the other hand we have an uncooperative Earth that is slowing down gradually over time. So in 1972, an international standards body introduced Coordinated Universal Time (UTC), a time standard like GMT that is still based on the Earth’s rotation, but with adjustments based on the precision of the underlying TAI.

The idea is to allow UTC to follow the irregularity of the Earth’s rotational period but also to adjust it at regular intervals according to TAI. These adjustments are made by an international committee charged with tracking this deviation. Twice a year they decide whether to add leap seconds to UTC. It’s all incredibly nerdy and I’m going to leave that right there. But as I write this on 1 Feb 2020, I took note of the time as it reached 16:39:00 TAI. At that exact moment it was also 16:38:23 UTC. So TAI is currently 37 seconds ahead of UTC. (This is because TAI was 10 seconds ahead of UTC when it was adopted in 1972 and there have been 27 leap seconds since then.)

That’s all well and good for civil time. But in astronomy, in order to calculate eclipses or transits you need to know the exact orbital positions of the Sun, Moon and planets at a given moment in time. So we need yet another time scale called Terrestrial Time (TT). This is another uniform time scale that works over long distances and long periods of time. Along with TT is a current baseline called the J2000.0 epoch. Both TT and the current epoch will come into play shortly.

Finally that brings me to Delta T. From the wiki:

In precise timekeepingΔT (Delta Tdelta-TdeltaT, or DT) is a measure of the cumulative effect of the departure of the Earth’s rotation period from the fixed-length day of atomic time. Formally it is the time difference obtained by subtracting Universal Time (UT, defined by the Earth’s rotation) from Terrestrial Time (TT, independent of the Earth’s rotation): ΔT = TT − UT. The value of ΔT for the start of 1902 is approximately zero; for 2002 it is about 64 seconds. So the Earth’s rotations over that century took about 64 seconds longer than would be required for days of atomic time.

In his book Astronomical Algorithms, Meeus lists a general algorithm you can follow to calculate an approximate value of ΔT. Here we let t be the time (in centuries which is why we divide by 100) measured from the current J2000.0 epoch:

$$t = {year -2000 \over 100}$$

For years after 2000 Meeus introduces another calculation published in a paper in Paris written by Chapront, Chapront-Touzé & Francou (1997):

$$ΔT = 102 + 102 t + 25.3 t^2$$

I do not have access to their paper so I can’t explain the use of these “magic numbers” but you can see them in use on this Delta T calculator on a retro 90s-era (but awesome) web site maintained by Professor van Gent of the University of Utrecht.

Next we add the following correction for years between 2000 and 2100:

$$0.37 * (year – 2100)$$

That’s it. Let’s implement this algorithm for the year 2020. I’m using C# language syntax and Visual Studio Code as my editor:

class Program
{
    static void Main(string[] args)
    {
        var year = 2020;
        var t = (year - 2000) / 100;
        var deltaT = 102 + (102 * t) + (25.3 * Math.Pow(t, 2));
        var correction = 0.37 * (year - 2100);
        var correctedDeltaT = deltaT + correction;

        Console.WriteLine("Corrected DT for the year {0} is: {1}", 
            year, 
            correctedDeltaT);
    }
}

And my output is +72.4 seconds. (Delta T is always expressed in positive terms.) So at the start of 2020 the Earth’s rotation was roughly 72.4 seconds behind UTC.

Now you’ll recall in my discussion of TAI and UTC above that when TAI was at 16:39:00 UTC was at 16:38:23, a difference of 37 seconds. What’s the relationship here?  Well, TT is always 32.184 seconds ahead of TAI. Why 32.184 seconds? Just take my word for it. It involves an offset from the start of TAI that takes into account an earlier time standard. With that bit of hand-waving out of the way here are our relationships when TAI was at 16:39:00:

TT == 16:39:32.184
TAI == 16:39:00
UTC == 16:38:23

So TT is about 69 seconds ahead of UTC. But if ΔT = TT − UT results in 69 what about our algorithm above that produced 72.4? Remember the leap second additions? Every time a leap second is added to UTC the gap between TT and UTC widens. Three leap seconds were added since the J2000.0 epoch began (2005, 2008, and 2016). This is why our algorithm uses J2000.0 epoch as its baseline.

NASA has published another algorithm for deriving ΔT that is specific for years between 2005 and 2050. If I implement their algorithm into C# code we get this:

private static double ApproximateDeltaT(int year)
{
    var t = year - 2000;
    return 62.92 + 0.32217 * t + 0.005589 * t * t;
}

This functions returns 71.599. These are approximations so we mustn’t worry about why the  two algorithms don’t produce the exact same results. The point to all of this is to know exactly when in TT an event like the next solar eclipse will occur so that you can convert it to UTC in order to observe it!

Kepler’s Equation

Introduction

Since at least the time of Hipparchus and Eudoxus the ancient world believed that the Sun, moon, planets, and stars moved around the earth in circular orbits. Aristotle’s Physics put the heavenly bodies on perfect crystal spheres. This theory was further formalized in the second century by Ptolemy in his Almagest which served the basic needs of astronomers for the next 1,500 years. Over that time, Aristotelian physics became an article of faith not to be questioned.

But there were unexplained irregularities that never quite fit the theory. The ancients knew that the seasons were of different lengths. For example, Winter is about 89 days in length while Summer is about 94 days. Why did the Sun sometimes speed up and slow down like that? They also noticed that Mars would move East across the night sky for a few years, slow down, and then reverse course and move West for a few months before looping back to its original course. (What we now call apparent retrograde motion.)

Continue reading “Kepler’s Equation”

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.