Smoky Cogs

As explained in my previous tutorial, equatorial coordinates are based on a projection of the Earth’s equator to demarcate position in the celestial sphere. Now we move on to ecliptic coordinates.

The Earth’s axis of rotation is not quite the same as the plane in which it revolves around the Sun. The Earth’s axis is tilted by approximately 23.4 degrees to this plane, called the ecliptic plane. This tilt is what causes our seasons and other related phenomena.

Now, ecliptic coordinates uses the ecliptic plane as its base plane, much like the equatorial coordinate system uses the Earth’s equator as it’s base. This gives the position of the sky relative to this plane, which can be very useful when doing calculations which include the Sun, since the Earth and the Sun would both be on the base plane.

The ecliptic latitude is denoted by the β symbol, and the ecliptic longitude is denoted by λ.

To convert between equatorial and ecliptic coordinates, we only need the obliquity of the ecliptic (the angle of the Earth’s tilt) as additional information to calculate the conversion.

		public static void ConvEquToEcl(double fOblique, double fRA, double fDecl, ref double fELong, ref double fELat)
		{
			double fX;
			double fY;
			double fSinELat;

			fDecl = Trig.DegToRad(fDecl);
			fRA = Trig.DegToRad(fRA * 15);
			fOblique = Trig.DegToRad(fOblique);
			fSinELat = (Math.Sin(fDecl) * Math.Cos(fOblique)) - (Math.Cos(fDecl) * Math.Sin(fOblique) * Math.Sin(fRA));
			fELat = Math.Asin(fSinELat);
			fY = (Math.Sin(fRA) * Math.Cos(fOblique)) + (Math.Tan(fDecl) * Math.Sin(fOblique));
			fX = Math.Cos(fRA);
			fELong = Math.Atan(fY / fX);
			fELong = Trig.RadToDeg(fELong);
			fELat = Trig.RadToDeg(fELat);
			fELong = Trig.TanQuadrant(fX, fY, fELong);
		}

		public static void ConvEclToEqu(double fOblique, double fELong, double fELat, ref double fRA, ref double fDecl)
		{
			double fX;
			double fY;
			double fSinDecl;

			fELong = Trig.DegToRad(fELong);
			fELat = Trig.DegToRad(fELat);
			fOblique = Trig.DegToRad(fOblique);
			fSinDecl = (Math.Sin(fELat) * Math.Cos(fOblique)) + (Math.Cos(fELat) * Math.Sin(fOblique) * Math.Sin(fELong));
			fDecl = Math.Asin(fSinDecl);
			fY = (Math.Sin(fELong) * Math.Cos(fOblique)) - (Math.Tan(fELat) * Math.Sin(fOblique));
			fX = Math.Cos(fELong);
			fRA = Math.Atan(fY / fX);
			fRA = Trig.RadToDeg(fRA);
			fDecl = Trig.RadToDeg(fDecl);
			fRA = Trig.TanQuadrant(fX, fY, fRA);
			fRA = fRA / 15.0;
		}

We now begin with a vital part of astronomical calculations – how to convert between different coordinate systems.

Each coordinate system has advantages and disadvantages for different areas of work, so each is used where applicable. The main coodinate systems in use in astronomy are the Equatorial, Horizontal, Ecliptic and Galactic coordinate systems.

The most important of these for astronomical work are Equatorial coordinates, so I tend to use this as a basis to convert all the other coordinates to and from each other. In this tutorial I will handle the conversion between equatorial and horizontal coodinates, and then cover the others in later tutorials, but first some definitions.

The equatorial coordinate system is the projection of the normal geographic coordinates on the earth. So, what this means, is that the lines of longitude and latitude on the earth are projected to an imaginary celestial sphere, with the celestial equator above the earth’s equator, and the celestial poles above the earth’s poles.

The sphere is split up into longitude and latitude, just like the Earth is. The longitude is known as the right ascension, and is measured in hours, with the zero point being the First Point of Aries, which is the point where the Sun is in at the Vernal Equinox (on March 21st). This value will change over time due to precession and nutation, but that is out of the scope of this tutorial.

The latitude is known as declination (denoted by δ), and is the distance in degrees that an object is away from the celestial equator, with the poles being 90 degrees, just like latitude on the earth’s surface.

Now, let us look at horizontal coordinates, also known as alt-azimuth coordinates. Here, the coodinates in the sky are demarcated by the height about the horizon (altitude) expressed in degrees, with 90 degrees being vertically upwards, and the azimuth, which is the direction to the object in relation to north, expressed in degrees in the range 0-360 degrees.

Before we can convert between the two coordinate systems, we need to find the hour angle from the right ascension. The right ascension is based on the star’s position in the celestial sphere, whereas the hour angle converts that value to a value based on the earth’s longitude and latitude.

This is done by calculating the Local Sidereal Time, using the Universal Time and longitude of the observer’s position, and then subtracting the right ascension from that to get the hour angle.
continue reading…

Before I move on to other stuff, let us talk about right ascension. All angular representations in astronomy use degrees, except for one. This is the right ascension, the longitude component of the equatorial coordinate system, the full detail of which I will cover in more detail later. This value is stored in hours rather than in degrees, so needs a bit of special consideration.

Since hours, minutes and seconds are very similar to degrees, minutes and seconds, the class to handle it is very similar to what we discussed in the previous tutorial. The major difference between the two is that the hours are in the range 0-24 hours, while degrees are in the range 0-360.

To convert hours to degrees, 1 degree = 4 minutes (of time), and was handled in the tutorial covering the basics of time.
continue reading…

A lot of astronomy is based on mapping positions of objects in either the celestial sphere, or ground based positions (ie longitude and latitude). Because of this, almost every astronomical calculation uses the coordinates – expressed in degrees – of objects to calculate their various attributes.

Fortunately, most people are used to working with degrees, although, the problem is that degrees are not exactly that easy to work with in calculations.

Degrees are made up of 60 minutes, which in turn are made up of 60 seconds per minute. With all the calculations in the application, we need a decimal format of the degree value so that we can easily work with it.

Before I move on to further topics covering coodinate system, let me introduce you to a class that works with degrees, and allows easy conversion between DMS and decimal formats.

The UraniaDMS class encapsulates all the functionality required to handle degrees.

Most of the class is pretty much self-explanatory. There are functions to set and set the degrees, minutes and seconds individually, as well as a whole. The class stores the values for degree, minute and second separately, as well as a separate value for the sign.

The most complicated function is the function that sets the value as a decimal – setDec(). It sets the degrees, minutes and seconds in turn, each time taking the whole part nad then multiplying the fractional part by 60.

There is also a getString() function that outputs a nice friendly representation of the value.
continue reading…

Easter is a rather tricky beast to calculate. It is usually the first Sunday after the fourteenth day after the first new moon after the Vernal Equinox (March 21st). Quite a mouthful that.

There have been numerous ways worked out to find the date. The method listed below is an adaptation of the method appearing in Butcher’s Ecclesiastical Calendar in 1876. It works for the Gregorian calendar, so cannot be used to find dates before 1583.

		public static DateTime CalcEasterDate(int piYear)
		{
			DateTime dEaster = new DateTime();

			double iA;
			double iB;
			double iC;
			double iD;
			double iE;
			double iFA;
			double iG;
			double iH;
			double iK;
			double iI;
			double iL;
			double iM;
			double iNA;
			double iP;

			iA = (piYear % 19);
			iB = Math.Floor(piYear / 100.0);
			iC = (piYear % 100);
			iD = Math.Floor(iB / 4.0);
			iE = (iB % 4);
			iFA = Math.Floor((iB +8) / 25.0);
			iG = Math.Floor((iB - iFA + 1) / 3);
			iH = ((19 * iA) + iB - iD - iG + 15) % 30;
			iI = Math.Floor(iC / 4.0);
			iK = (iC % 4);
			iL = (32 + (2 * iE) + (2 * iI) - iH - iK) % 7;
			iM = Math.Floor((iA + (11 * iH) + (22 * iL)) / 451.0);
			iNA = Math.Floor((iH + iL - (7 * iM) + 114) / 31.0);
			iP = (iH + iL - (7 * iM) + 114) % 31;

			dEaster = new DateTime(piYear, Convert.ToInt16(iNA), Convert.ToInt16(iP + 1));
			return dEaster;
		}

Sidereal time is the time in relation to the stars. Normal time is calculated as the time it takes for the earth to spin once about its axis to the same point in relation to the Sun, so for example from when the Sun is directly overhead to the next time the Sun is directly overhead. The total time for this to happen is 24 hours.

If we instead use the stars as our reference point, the length of a day is slightly shorter, by 4 minutes, thus making a sidereal day 23 hours 54 minutes.

The reason for this is that during the day, the Earth has moved 1/365th of the way further in its orbit around the sun, so once the Earth has made one full rotation, it has to rotate just that little bit more to have the Sun in the same position.

What this means though as well is that a year is made up of 366 sidereal days as opposed to 365 days in the normal year (rounded off here).

We also talk about Greenwich Sidereal Time (GST) and Local Sideral Time (LST) which are analogous to Universal Time and Local Mean Time.

Finding the GST, we work with the Julian day, which was covered in a previous tutorial. Working from an epoch (a date where we know where everything is) Julian Day of 2451545.0, where we know the Sidereal time was the same as the UT. First we need work out the time difference from the current date to that date, and once we have this we can calculate the amount of precession to take into account, which I will talk more about in a future tutorial, which is stored in fT0. We then convert the UT to a decimal format, so that we can multiply by the multiplication factor to put the UT into GST, and then add in the precessional amount.

We then need to put the final answer into the 24 hour range before returning the result.

		public static DateTime CalcGSTFromUT(DateTime dDate)
		{
			double fJD;
			double fS;
			double fT;
			double fT0;
			DateTime dGST;
			double fUT;
			double fGST;

			fJD = GetJulianDay(dDate.Date, 0);
			fS = fJD - 2451545.0;
			fT = fS / 36525.0;
			fT0 = 6.697374558 + (2400.051336 * fT) + (0.000025862 * fT * fT);
			fT0 = Trig.PutIn24Hour(fT0);
			fUT = ConvTimeToDec(dDate);
			fUT = fUT * 1.002737909;
			fGST = fUT + fT0;
			while (fGST > 24)
			{
				fGST = fGST - 24;
			}
			dGST = ConvDecTUraniaTime(fGST);
			return dGST;
		}

		public static DateTime ConvDecTUraniaTime(double fTime)
		{
			DateTime dDate;

			dDate = new DateTime();
			dDate = dDate.AddHours(fTime);
			return (dDate);
		}

		public static double ConvTimeToDec(DateTime dDate)
		{
			double fHour;

			fHour = dDate.Hour + (dDate.Minute / 60.0) + (dDate.Second / 60.0 / 60.0) + (dDate.Millisecond / 60.0 / 60.0 / 1000.0);
			return fHour;
		}

Now converting the other way round, we just do the reverse. We work out fT0 as before, and then subtract it from the GST. We then put it into 24 hours, before multiplying by the inverse of the scaling factor we used in when going in the other direction.

We then have the UT.

		public static DateTime CalcUTFromGST(DateTime dGSTDate, DateTime dCalDate)
		{
			double fJD;
			double fS;
			double fT;
			double fT0;
			double fGST;
			double fUT;
			DateTime dUT;

			bool bPrevDay;
			bool bNextDay;

			bPrevDay = false;
			bNextDay = false;

			fJD = GetJulianDay(dCalDate.Date, 0);
			fS = fJD - 2451545.0;
			fT = fS / 36525.0;
			fT0 = 6.697374558 + (2400.051336 * fT) + (0.000025862 * fT * fT);

			fT0 = Trig.PutIn24Hour(fT0);

			fGST = ConvTimeToDec(dGSTDate);
			fGST = fGST - fT0;
		
			while (fGST > 24)
			{
				fGST = fGST - 24;
				bNextDay = true;
			}
			while (fGST < 0)
			{
				fGST = fGST + 24;
				bPrevDay = true;
			}

			fUT = fGST * 0.9972695663;
			// fUT = fGST
		
			dUT = dCalDate.Date;
			dUT = dUT.AddHours(fUT);
			fUT = dUT.Millisecond;
			
			if (bNextDay == true)
			{
				dUT = dUT.AddDays(1);
			}

			if (bPrevDay == true)
			{
				dUT = dUT.Subtract(new TimeSpan(1, 0, 0, 0, 0));
			}
			return dUT;
		}

The Julian day is the standard form of reporting the date in astronomy. One Julian day is equal to one normal day, with the fractional part being the time of day expressed as a fraction of a day.

The Julian day starts on 1st January 4713BC, at Greenwich noon. This means that the Julian day is offset from a normal day by half, with the day starting at noon.

The current value of the Julian day at the time of writing this tutorial is 2455117.64.

The algorithm for calculating the Julian day for Gregorian dates is
JD = (1461 × (Y + 4800 + (M − 14)/12))/4 +(367 × (M − 2 − 12 × ((M − 14)/12)))/12 − (3 × ((Y + 4900 + (M – 14)/12)/100))/4 + D − 32075

and converting Julian dates, the formula is
JD = 367 × Y − (7 × (Y + 5001 + (M − 9)/7))/4 + (275 × M)/9 + D + 1729777

We then have to add in the fractional time portion which is defined by
JD = JD + (hour – 12) / 24 + minute / 1440 + second / 86400.

The functions below also pass the time zone, thus converting the dates (in Local Zone Time) UT, and vice versa.

		public static double GetJulianDay(DateTime dDate, int iZone)
		{
			double fJD;
			double iYear;
			double iMonth;
			double iDay;
			double iHour;
			double iMinute;
			double iSecond;
			double iGreg;
			double fA;
			double fB;
			double fC;
			double fD;
			double fFrac;

			dDate = CalcUTFromZT(dDate, iZone);

			iYear = dDate.Year;
			iMonth = dDate.Month;
			iDay = dDate.Day;
			iHour = dDate.Hour;
			iMinute = dDate.Minute;
			iSecond = dDate.Second;
			fFrac = iDay + ((iHour + (iMinute / 60) + (iSecond / 60 / 60)) / 24);
			if (iYear < 1582)
			{
				iGreg = 0;
			}
			else
			{
				iGreg = 1;
			}
			if ((iMonth == 1) || (iMonth == 2))
			{
				iYear = iYear - 1;
				iMonth = iMonth + 12;
			}

			fA = (long)Math.Floor(iYear / 100);
			fB = (2 - fA + (long)Math.Floor(fA / 4)) * iGreg;
			if (iYear < 0)
			{
				fC = (int)Math.Floor((365.25 * iYear) - 0.75);
			}
			else
			{
				fC = (int)Math.Floor(365.25 * iYear);
			}
			fD = (int)Math.Floor(30.6001 * (iMonth + 1));
			fJD = fB + fC + fD + 1720994.5;
			fJD = fJD + fFrac;
			return fJD;
		}

Finding the date from the Julian day is found simply by inverting the calculations done above

		public static DateTime getDateFromJD(double fJD, int iZone)
		{
			DateTime dDate;
			int iYear;
			int fMonth;
			int iDay;
			int iHour;
			int iMinute;
			int iSecond;
			double fFrac;
			double fFracDay;
			int fI;
			int fA;
			int fB;
			int fC;
			int fD;
			int fE;
			int fG;

			fJD = fJD + 0.5;
			fI = (int)Math.Floor(fJD);
			fFrac = fJD - fI;

			if (fI > 2299160)
			{
				fA = (int)Math.Floor((fI - 1867216.25) / 36524.25);
				fB = fI + 1 + fA - (int)Math.Floor((double)(fA / 4));
			}
			else
			{
				fA = 0;
				fB = fI;
			}
			fC = fB + 1524;
			fD = (int)Math.Floor((fC - 122.1) / 365.25);
			fE = (int)Math.Floor(365.25 * fD);
			fG = (int)Math.Floor((fC - fE) / 30.6001);
			fFracDay = fC - fE + fFrac - (long)Math.Floor((double)(30.6001 * fG));
			iDay = (int)Math.Floor(fFracDay);
			fFracDay = (fFracDay - iDay) * 24;
			iHour = (int)Math.Floor(fFracDay);
			fFracDay = (fFracDay - iHour) * 60;
			iMinute = (int)Math.Floor(fFracDay);
			fFracDay = (fFracDay - iMinute) * 60;
			iSecond = (int)Math.Floor(fFracDay);

			if (fG < 13.5)
			{
				fMonth = fG - 1;
			}
			else
			{
				fMonth = fG - 13;
			}

			if (fMonth > 2.5)
			{
				iYear = (int)Math.Floor((double)(fD - 4716.0));
			}
			else
			{
				iYear = (int)Math.Floor((double)(fD - 4715.0));
			}


			dDate = new DateTime(iYear, (int)Math.Floor((double)fMonth), iDay, iHour, iMinute, iSecond);
			dDate = CalcZTFromUT(dDate, iZone);
			return dDate;
		}