Smoky Cogs

Like with what we did the Moon phase, the method of calculating the Moon age is the same, for the most part, as finding the Moon’s position.

The only difference is in finding the formula for finding the distance
Distance = (Semi-major axis) * (1 – (ecc * ecc)) / (1 + (ecc * Cos(MM1 + MEC)))
where MM1 and MEC come from the calculation

		public static void CalcMoonDistance(DateTime dDate, DateTime dEpoch, double fMEpochLong, double fMPeriLong, double fMAscNode, double fMIncl, double fMEcc, double fSEpochEclLong, double fSPeriEclLong, double fSEcc, double fMSMA, ref double fMDistance)
		{
			double fN, fSM, fSE, fSLambda;
			double fL, fMM, fMN, fME, fAE, fMEC, fA3, fMM1;
			double fJD1, fJD2, fDays;

			fJD1 = UraniaTime.GetJulianDay(dDate, 0);
			fJD2 = UraniaTime.GetJulianDay(dEpoch, 0);
			fDays = (fJD1 - fJD2);
			fDays += 1;

			fN = (360.0/365.242191) * fDays;
			fN = Trig.PutIn360Deg(fN);
			fSM = fN + fSEpochEclLong - fSPeriEclLong;
			fSM = Trig.PutIn360Deg(fSM);

			fSE = (360.0 / Math.PI) * fSEcc * Math.Sin(Trig.DegToRad(fSM));
			fSLambda = fN + fSE + fSEpochEclLong;

			fL = (13.176396 * fDays) + fMEpochLong;
			fL = Trig.PutIn360Deg(fL);
			
			fMM = fL - (0.111404 * fDays) - fMPeriLong;
			fMM = Trig.PutIn360Deg(fMM);

			fMN = fMAscNode - (0.0529539 * fDays);
			fMN = Trig.PutIn360Deg(fMN);

			fME = 1.2739 * Trig.Sin((2.0 * (fL - fSLambda)) - fMM);
			fAE = 0.1858 * Trig.Sin(fSM);
			fA3 = 0.37 * Trig.Sin(fSM);

			fMM1 = fMM + fME - fAE + fA3;
			
			fMEC = 6.2886 * Trig.Sin(fMM1);

			fMDistance = fMSMA * ((1.0 - (fMEcc * fMEcc))/(1.0 + (fMEcc * Trig.Cos(fMM1 + fMEC))));
		}

Once we have Moon’s position, as we did in the last post, it is very easy to calculate the phase of the Moon.

The calculation is identical, except for the last two lines, where we get the age of the Moon, and then get the phase with the formula
Phase = 0.5 * (1 – Cos(Age))

The resulting answer is in the range 0 – 1.

		public static void CalcMoonPhase(DateTime dDate, DateTime dEpoch, double fMEpochLong, double fMPeriLong, double fMAscNode, double fMIncl, double fMEcc, double fSEpochEclLong, double fSPeriEclLong, double fSEcc, ref double fMPhase)
		{
			double fN, fSM, fSE, fSLambda;
			double fL, fMM, fMN, fME, fAE, fMEC, fA3, fA4, fMV, fMM1, fL1, fL2;
			double fJD1, fJD2, fDays, fMD;

			fJD1 = UraniaTime.GetJulianDay(dDate, 0);
			fJD2 = UraniaTime.GetJulianDay(dEpoch, 0);
			fDays = (fJD1 - fJD2);
			fDays += 1;

			fN = (360.0/365.242191) * fDays;
			fN = Trig.PutIn360Deg(fN);
			fSM = fN + fSEpochEclLong - fSPeriEclLong;
			fSM = Trig.PutIn360Deg(fSM);

			fSE = (360.0 / Math.PI) * fSEcc * Math.Sin(Trig.DegToRad(fSM));
			fSLambda = fN + fSE + fSEpochEclLong;

			fL = (13.176396 * fDays) + fMEpochLong;
			fL = Trig.PutIn360Deg(fL);
			
			fMM = fL - (0.111404 * fDays) - fMPeriLong;
			fMM = Trig.PutIn360Deg(fMM);

			fMN = fMAscNode - (0.0529539 * fDays);
			fMN = Trig.PutIn360Deg(fMN);

			fME = 1.2739 * Trig.Sin((2.0 * (fL - fSLambda)) - fMM);
			fAE = 0.1858 * Trig.Sin(fSM);
			fA3 = 0.37 * Trig.Sin(fSM);

			fMM1 = fMM + fME - fAE + fA3;
			
			fMEC = 6.2886 * Trig.Sin(fMM1);
			fA4 = 0.214 * Trig.Sin(2.0 * fMM1);
			fL1 = fL + fME + fMEC - fAE + fA4;

			fMV = 0.6583 * Trig.Sin(2.0 * (fL1 - fSLambda));
			fL2 = fL1 + fMV;

			fMD = fL2 - fSLambda;
			fMPhase = 0.5 * (1.0 - Trig.Cos(fMD));
		}

Finding the position of the Moon is in many ways similar to finding the position of the Sun.

The Moon moves along a path that is near the ecliptic but is slightly offset – by 5.14 degrees. If this were not so, then we would get eclipses every new and full moon.

The other information we need to calculate the Moon’s position, is the ascending node, which is the point where the orbit of the moon crosses the ecliptic, which is 318.351648 degrees, and the eccentricity of the orbit, which is 0.054900. The perigee of the Moon is also needed which is 36.340410 degrees.

We then need to know the position of the Sun too, so we also need the solar ecliptic longitude, perigee, and eccentricity.

The last thing we need now is a starting point for our calculation, so for the epoch given, we need to supply the ecliptic longitude for that time, from which the current position can be calculated.

		public static void CalcMoonPos(DateTime dDate, DateTime dEpoch, double fMEpochLong, double fMPeriLong, double fMAscNode, double fMIncl, double fMEcc, double fSEpochEclLong, double fSPeriEclLong, double fSEcc, ref double fMRA, ref double fMDecl)
		{
			double fN, fSM, fSE, fSLambda;
			double fL, fMM, fMN, fME, fAE, fMEC, fA3, fA4, fMV, fMM1, fL1, fL2, fN1, fX, fY;
			double fT, fMLambda, fMBeta, fJD1, fJD2, fDays;

			fJD1 = UraniaTime.GetJulianDay(dDate, 0);
			fJD2 = UraniaTime.GetJulianDay(dEpoch, 0);
			fDays = (fJD1 - fJD2);
			fDays += 1;

			fN = (360.0/365.242191) * fDays;
			fN = Trig.PutIn360Deg(fN);
			fSM = fN + fSEpochEclLong - fSPeriEclLong;
			fSM = Trig.PutIn360Deg(fSM);

			fSE = (360.0 / Math.PI) * fSEcc * Math.Sin(Trig.DegToRad(fSM));
			fSLambda = fN + fSE + fSEpochEclLong;

			fL = (13.176396 * fDays) + fMEpochLong;
			fL = Trig.PutIn360Deg(fL);
			
			fMM = fL - (0.111404 * fDays) - fMPeriLong;
			fMM = Trig.PutIn360Deg(fMM);

			fMN = fMAscNode - (0.0529539 * fDays);
			fMN = Trig.PutIn360Deg(fMN);

			fME = 1.2739 * Trig.Sin((2.0 * (fL - fSLambda)) - fMM);
			fAE = 0.1858 * Trig.Sin(fSM);
			fA3 = 0.37 * Trig.Sin(fSM);

			fMM1 = fMM + fME - fAE + fA3;
			
			fMEC = 6.2886 * Trig.Sin(fMM1);
			fA4 = 0.214 * Trig.Sin(2.0 * fMM1);
			fL1 = fL + fME + fMEC - fAE + fA4;

			fMV = 0.6583 * Trig.Sin(2.0 * (fL1 - fSLambda));
			fL2 = fL1 + fMV;

			fN1 = fMN - (0.16 * Trig.Sin(fSM));
			fY = Trig.Sin(fL2 - fN1) * Trig.Cos(fMIncl);
			fX = Trig.Cos(fL2 - fN1);

			fT = Trig.Atan(fY / fX);

			fT = Trig.TanQuadrant(fX, fY, fT);
			
			fMLambda = fT + fN1;
			fMBeta = Trig.Asin(Trig.Sin(fL2 - fN1) * Trig.Sin(fMIncl));
			UraniaCoord.ConvEclToEqu(23.441884,fMLambda, fMBeta, ref fMRA, ref fMDecl);
		}

Ok. So maybe this is not quite strictly and astronomical calculation, but the calculation is valid for any spherical object, so is equally applicable for finding the distance to the horizon on Earth, as it is on the Moon, or Mars, so have decided to include it with my astronomical calculations.

The distance to the horizon is dependent on the height above the surface of the sphere, which, in most cases which we care about would be the surface of the Earth.

So, given a height, h, in metres, and a radius, r (which for the Earth is 6371km), we can use the following formula to find the distance (in km) to the horizon:

D = √(2 * r) * √(h / 1000).

The inverse of this can be used to find the height required to see a horizon of particular distance.

	public class UraniaHorizon
	{

		public static double fEarthRadius = 6371.0;
		
		public static double HorizonDistance(double pdHeight, double pdRadius)
		{
			if (pdRadius == 0)
			{
				pdRadius = UraniaHorizon.fEarthRadius;
			}
			return System.Math.Sqrt(pdRadius * 2.0) * System.Math.Sqrt(pdHeight / 1000.0);
		}
		
		public static double HorizonHeight(double pdDistance, double pdRadius)
		{
			if (pdRadius == 0)
			{
				pdRadius = UraniaHorizon.fEarthRadius;
			}
			return ((pdDistance / System.Math.Sqrt(pdRadius * 2.0)) *(pdDistance / 112.88)) * 1000.0;
		}
	}

The Carrington Rotation Number is defined as the number of rotations which the Sun has done since 9 November 1853.

The time for one rotation of the Sun (in essence, one “day” on the Sun) is 27.2753 days.

The calculation is based on the fact that rotation 1690 started on 27 Dec 1979, so what it does is works out the number of days between that date and the date you want the CRN for, and divide by 27.2753, and then adding 1690.

		public static double CalcCRN(DateTime dDate)
		{
			double fCRN;
			double fJulian;

			fJulian = UraniaTime.GetJulianDay(dDate, 0);
			fCRN = 1690 + ((fJulian - 2444235.34)/27.2753);
			return fCRN;
		}

The method of find the angular diameter of the Sun is essentially the same as that of finding the distance to the Sun.

We need to find the true anomaly (v) and then using that we plug it into the formula
Angular Diameter = (angular Diameter at perigee) * ((1 + (eccentricy * cos(v)) / (1 – eccentricity²))

The resulting angular diameter is expressed in degrees, if the original angular diameter is also expressed in degrees.

		public static double CalcSunAngDiam(DateTime dDate, DateTime dEpoch)
		{
			double fD;
			double fN;
			double fM;
			double fE;
			double fTanV2;
			double fV;
			double fAngDiam;
			double fSolarMEL;
			double fSolarPL;
			double fSunEarthEcc;
			double fAcc;
			double fOblique;
			double fAngDiam0;
			

			fAngDiam0 = 0.533128;
			fAcc = 0.0000001;
			fD = UraniaTime.GetDaysBetween(dDate, dEpoch);
			fSolarMEL = GetSolarMEL(dEpoch, true);
			fSolarPL = GetSolarPerigeeLong(dEpoch, true);
			fSunEarthEcc = GetSunEarthEcc(dEpoch, true);
			fOblique = GetEarthObliquity(dEpoch, true);

			fN = (360.0 / 365.242191) * fD;
			fN = Trig.PutIn360Deg(fN);
			fM = fN + fSolarMEL - fSolarPL;
			fM = Trig.PutIn360Deg(fM);
			fM = Trig.DegToRad(fM);
			fE = CalcEccentricAnomaly(fM, fM, fSunEarthEcc, fAcc);
			fTanV2 = Math.Sqrt((1.0 + fSunEarthEcc) / (1.0 - fSunEarthEcc)) * Math.Tan(fE / 2.0);
			fV = Math.Atan(fTanV2) * 2.0;
			fAngDiam = fAngDiam0 * ((1.0 + (fSunEarthEcc * Math.Cos(fV))) / (1.0 - (fSunEarthEcc * fSunEarthEcc)));
			return fAngDiam;
		}

Finding the distance to the Sun is quite straightforward once we know the true anomaly of the Sun, which we did in the previous post in the quest to find the Sun’s position.

In fact the first half of the code is pretty much identical.

The one thing left to do is find the distance using the formula

Distance = (Semi-major axis * (1 – v * v)) / (1 + (eccentricity * Cos(v)))

which gives a distance in the same units as the semi-mahor axis, which is usually kilometres.

		public static double CalcSunDistance(DateTime dDate, DateTime dEpoch)
		{
			double fD;
			double fN;
			double fM;
			double fE;
			double fTanV2;
			double fV;
			double fDistance;
			double fSolarMEL;
			double fSolarPL;
			double fSunEarthEcc;
			double fAcc;
			double fOblique;
			double fSMA;

			fAcc = 0.0000001;
			fD = UraniaTime.GetDaysBetween(dDate, dEpoch);
			fSolarMEL = GetSolarMEL(dEpoch, true);
			fSolarPL = GetSolarPerigeeLong(dEpoch, true);
			fSunEarthEcc = GetSunEarthEcc(dEpoch, true);
			fOblique = GetEarthObliquity(dEpoch, true);
			fSMA = 149598500.0;

			fN = (360.0 / 365.242191) * fD;
			fN = Trig.PutIn360Deg(fN);
			fM = fN + fSolarMEL - fSolarPL;
			fM = Trig.PutIn360Deg(fM);
			fM = Trig.DegToRad(fM);
			fE = CalcEccentricAnomaly(fM, fM, fSunEarthEcc, fAcc);
			fTanV2 = Math.Sqrt((1.0 + fSunEarthEcc) / (1.0 - fSunEarthEcc)) * Math.Tan(fE / 2.0);
			fV = Math.Atan(fTanV2) * 2.0;
			fDistance = (fSMA * (1.0 - (fSunEarthEcc * fSunEarthEcc))) / (1.0 + (fSunEarthEcc * Math.Cos(fV)));

			return fDistance;
		}