Skip to content

Archive

Tag: Astronomy

This calculation is based on the position calculation, but also includes finding the phase and distance, which are required for finding the magnitude.

The formula for the magnitude is quite straightforward once we have these values, as shown below. The calculation needs the magnitude of the object at 1AU as the basis to calculate the magnitude at the specified distance.

		public static void CalcPlanetMag(DateTime dEpoch, DateTime dDate, double fPOrbPeriod, double fPEcc, double fPEpochLong, double fPPeriLong, double fPSMA, double fPAscNode, double fPIncl, double fEOrbPeriod, double fEEcc, double fEEpochLong, double fEPeriLong, double fESMA, bool bInnerPlanet,  double fMag1AU, ref double fMag)
		{
			double fD, fPN, fPM, fPL, fPV, fPR;
			double fEN, fEM, fEL, fEV, fER;
			double fPsi, fX, fY, fA, fL1, fR1, fLambda, fT, fPhase, fDistance;

			fD = UraniaTime.GetDaysBetween(dDate, dEpoch);
			fPN = (360.0/365.242191) * (fD / fPOrbPeriod);
			fPN = Trig.PutIn360Deg(fPN);
			fPM = fPN + fPEpochLong - fPPeriLong;
			fPL = fPN + ((360.0 / Math.PI) * fPEcc * Math.Sin(Trig.DegToRad(fPM))) + fPEpochLong;
			fPL = Trig.PutIn360Deg(fPL);
			fPV = fPL - fPPeriLong;
			fPR = (fPSMA * (1 - (fPEcc * fPEcc))) / (1 + (fPEcc * Math.Cos(Trig.DegToRad(fPV))));

			fEN = (360.0/365.242191) * (fD / fEOrbPeriod);
			fEN = Trig.PutIn360Deg(fEN);
			fEM = fEN + fEEpochLong - fEPeriLong;
			fEL = fEN + ((360.0 / Math.PI) * fEEcc * Math.Sin(Trig.DegToRad(fEM))) + fEEpochLong;
			fEL = Trig.PutIn360Deg(fEL);
			fEV = fEL - fEPeriLong;
			fER = ((1 - (fEEcc * fEEcc))) / (1 + (fEEcc * Math.Cos(Trig.DegToRad(fEV))));

			fPsi = Math.Asin(Math.Sin(Trig.DegToRad(fPL - fPAscNode)) * Math.Sin(Trig.DegToRad(fPIncl)));
			fPsi = Trig.RadToDeg(fPsi);

			fY = Math.Sin(Trig.DegToRad(fPL - fPAscNode));
			fX = Math.Cos(Trig.DegToRad(fPL - fPAscNode));
			fT = Math.Atan((fY/fX) * Math.Sin(Trig.DegToRad(fPIncl)));
			fT = Math.Atan((fY/fX));
			fT = Trig.RadToDeg(fT) ;
			fT = Trig.TanQuadrant(fX, fY, fT);
			fL1 = fT + fPAscNode;
			fL1 = Trig.PutIn360Deg(fL1);
			fR1 = fPR * Math.Cos(Trig.DegToRad(fPsi));

			if (bInnerPlanet == true)
			{
				fA = Math.Atan((fR1 * Math.Sin(Trig.DegToRad(fEL - fL1)))/(fER - (fR1 * Math.Cos(Trig.DegToRad(fEL - fL1)))));
				fA = Trig.RadToDeg(fA);
				fLambda = 180 + fEL + fA;
			}
			else
			{
				fLambda = Math.Atan((fER * Math.Sin(Trig.DegToRad(fL1 - fEL)))/(fR1 - (fER * Math.Cos(Trig.DegToRad(fL1 - fEL)))));
				fLambda = Trig.RadToDeg(fLambda) + fL1;
				fLambda = Trig.PutIn360Deg(fLambda);
			}

			fPhase = 0.5 * (1 + Math.Abs(Math.Cos(Trig.DegToRad(fLambda - fPL))));
			fDistance = Math.Sqrt((fER * fER) + (fPR * fPR) - (2.0 * fER * fPR * Math.Cos(Trig.DegToRad(fPL - fEL))));
	
			fMag = (5 * Math.Log10((fPR * fDistance)/(Math.Sqrt(fPhase)))) + fMag1AU;

		}
Share

The calculation for the distance of the Earth to a planet is also based on the position calculation, but is a little bit simpler. To get all the information we need, we do not need the full position calculation, and it is enough to find the heliographical coordinates of the planet and the Earth.

Once we have these values, we can then find the distance between the two.

		public static void CalcPlanetDistance(DateTime dEpoch, DateTime dDate, double fPOrbPeriod, double fPEcc, double fPEpochLong, double fPPeriLong, double fPSMA, double fPAscNode, double fPIncl, double fEOrbPeriod, double fEEcc, double fEEpochLong, double fEPeriLong, double fESMA,  ref double fDistance)
		{
			double fD, fPN, fPM, fPL, fPV, fPR;
			double fEN, fEM, fEL, fEV, fER, fRho;
			
			fD = UraniaTime.GetDaysBetween(dDate, dEpoch);
			fPN = (360.0/365.242191) * (fD / fPOrbPeriod);
			fPN = Trig.PutIn360Deg(fPN);
			fPM = fPN + fPEpochLong - fPPeriLong;
			fPL = fPN + ((360.0 / Math.PI) * fPEcc * Math.Sin(Trig.DegToRad(fPM))) + fPEpochLong;
			fPL = Trig.PutIn360Deg(fPL);
			fPV = fPL - fPPeriLong;
			fPR = (fPSMA * (1 - (fPEcc * fPEcc))) / (1 + (fPEcc * Math.Cos(Trig.DegToRad(fPV))));

			fEN = (360.0/365.242191) * (fD / fEOrbPeriod);
			fEN = Trig.PutIn360Deg(fEN);
			fEM = fEN + fEEpochLong - fEPeriLong;
			fEL = fEN + ((360.0 / Math.PI) * fEEcc * Math.Sin(Trig.DegToRad(fEM))) + fEEpochLong;
			fEL = Trig.PutIn360Deg(fEL);
			fEV = fEL - fEPeriLong;
			fER = ((1 - (fEEcc * fEEcc))) / (1 + (fEEcc * Math.Cos(Trig.DegToRad(fEV))));
			
			fRho = Math.Sqrt((fER * fER) + (fPR * fPR) - (2.0 * fER * fPR * Math.Cos(Trig.DegToRad(fPL - fEL))));
			fDistance = fRho;
		}
Share

Once we know how to find the position of the planet, it is very easy to find the phase of the planet.

The first part of the calculation is identical to the calculation to find the planet position.

The actual formula to find the phase is
Phase = 0.5 * (1 + Abs(Cos(Ecliptic longitude – heliographical longitude)))

This is the only addition to the calculation we did previously

		public static void CalcPlanetPhase(DateTime dEpoch, DateTime dDate, double fPOrbPeriod, double fPEcc, double fPEpochLong, double fPPeriLong, double fPSMA, double fPAscNode, double fPIncl, double fEOrbPeriod, double fEEcc, double fEEpochLong, double fEPeriLong, double fESMA, bool bInnerPlanet,  ref double fPhase)
		{
			double fD, fPN, fPM, fPL, fPV, fPR;
			double fEN, fEM, fEL, fEV, fER;
			double fPsi, fX, fY, fA, fL1, fR1, fLambda, fT;

			fD = UraniaTime.GetDaysBetween(dDate, dEpoch);
			fPN = (360.0/365.242191) * (fD / fPOrbPeriod);
			fPN = Trig.PutIn360Deg(fPN);
			fPM = fPN + fPEpochLong - fPPeriLong;
			fPL = fPN + ((360.0 / Math.PI) * fPEcc * Math.Sin(Trig.DegToRad(fPM))) + fPEpochLong;
			fPL = Trig.PutIn360Deg(fPL);
			fPV = fPL - fPPeriLong;
			fPR = (fPSMA * (1 - (fPEcc * fPEcc))) / (1 + (fPEcc * Math.Cos(Trig.DegToRad(fPV))));

			fEN = (360.0/365.242191) * (fD / fEOrbPeriod);
			fEN = Trig.PutIn360Deg(fEN);
			fEM = fEN + fEEpochLong - fEPeriLong;
			fEL = fEN + ((360.0 / Math.PI) * fEEcc * Math.Sin(Trig.DegToRad(fEM))) + fEEpochLong;
			fEL = Trig.PutIn360Deg(fEL);
			fEV = fEL - fEPeriLong;
			fER = ((1 - (fEEcc * fEEcc))) / (1 + (fEEcc * Math.Cos(Trig.DegToRad(fEV))));

			fPsi = Math.Asin(Math.Sin(Trig.DegToRad(fPL - fPAscNode)) * Math.Sin(Trig.DegToRad(fPIncl)));
			fPsi = Trig.RadToDeg(fPsi);

			fY = Math.Sin(Trig.DegToRad(fPL - fPAscNode));
			fX = Math.Cos(Trig.DegToRad(fPL - fPAscNode));
			fT = Math.Atan((fY/fX) * Math.Sin(Trig.DegToRad(fPIncl)));
			fT = Math.Atan((fY/fX));
			fT = Trig.RadToDeg(fT) ;
			fT = Trig.TanQuadrant(fX, fY, fT);
			fL1 = fT + fPAscNode;
			fL1 = Trig.PutIn360Deg(fL1);
			fR1 = fPR * Math.Cos(Trig.DegToRad(fPsi));

			if (bInnerPlanet == true)
			{
				fA = Math.Atan((fR1 * Math.Sin(Trig.DegToRad(fEL - fL1)))/(fER - (fR1 * Math.Cos(Trig.DegToRad(fEL - fL1)))));
				fA = Trig.RadToDeg(fA);
				fLambda = 180 + fEL + fA;
			}
			else
			{
				fLambda = Math.Atan((fER * Math.Sin(Trig.DegToRad(fL1 - fEL)))/(fR1 - (fER * Math.Cos(Trig.DegToRad(fL1 - fEL)))));
				fLambda = Trig.RadToDeg(fLambda) + fL1;
				fLambda = Trig.PutIn360Deg(fLambda);
			}

			fPhase = 0.5 * (1 + Math.Abs(Math.Cos(Trig.DegToRad(fLambda - fPL))));

		}
Share

Finding the position of a planet is a rather involved calculation.

The first step is to find the position of the planet in it’s orbit around the Sun.

We then find the heliocentric coordinates of the planet – which is essentially the longitude and latitude with reference to the Sun.

The last step involves converting these coordinates to ecliptic coordinates using Earth as the reference point.

There is a slight difference in the calculation between inner planets – planets closer to the Sun than Earth, and the outer planets, due to the positions of the orbits.

The parameters we need is a long list. We need the planet’s orbital period, eccentricity, longitude at the epoch, perigee, semi-major axis, ascending node and inclination. We also need the same values for Earth as well, and whether the planet is an inner planet or not.

continue reading…

Share

Here we are not going to try find the actual age of the Moon, which is several billion years old, but rather the age of the Moon since the last new Moon. This is just the amount of time in days, and fraction of days passed since the new Moon.

The Normalize() function just gets the fraction part of a number and if the resulting number is less than 1, returns 1 – number, otherwise it returns the resulting number itself.

The length of one lunar cycle is 29.53 days, so what the calculation does is works out how many lunar cycles have passed since a known new Moon. It then gets the fraction part, multiplied by the length of a lunar cycle, and we then have the age of the current lunar cycle in days.

		public static double CalcMoonAge(DateTime dDate, int iZone)
		{
			double fJD, fIP, fAge;

			fJD = UraniaTime.GetJulianDay(dDate, iZone);
			fIP = Normalize((fJD - 2451550.1) / 29.530588853);
			fAge = fIP*29.530588853;
            return fAge;
		}

		private static double Normalize(double fN)
		{
			fN = fN - Math.Floor(fN);
			if (fN < 0)
			{
				fN = fN + 1;
			}
			return fN;
		}
Share

Parallax, in the case of the Moon, is caused by the difference in location of the observer as the Earth turns, since over the course of an evening, the Earth has rotated halfway though its full rotation, and therefore the position of the observer is rather markedly moved.

The means that the apparent position of the Moon will shift slightly as the position of the observer shifts slightly.

The parallax is also determined by the distance of the Moon, so to find the angle of parallax, we need to know the parallax at perigee, and then we can use a similar formula to find the angular diameter
Parallax = (Parallax at perigee) * (Semi-major axis) / (Distance)

		public static void CalcMoonParallax(DateTime dDate, DateTime dEpoch, double fMEpochLong, double fMPeriLong, double fMAscNode, double fMIncl, double fMEcc, double fSEpochEclLong, double fSPeriEclLong, double fSEcc, double fMSMA, double fVParallax, ref double fMParallax)
		{
			double fRho;

			fRho = 0;
			CalcMoonDistance(dDate, dEpoch, fMEpochLong, fMPeriLong, fMAscNode, fMIncl, fMEcc, fSEpochEclLong, fSPeriEclLong, fSEcc, fMSMA, ref fRho);

			fMParallax = (fVParallax * fMSMA) / fRho;
		}
Share

Once we can find the distance to the Moon, we can use that to find the angular diameter.

Given that we know what the angular diameter is at perigee, we can find the angular diameter at any point in the orbit.

We plug the distance into the formula
Angular diameter = (Angular diameter at perigee) * (Semi-major axis) / (Distance)

		public static void CalcMoonDiam(DateTime dDate, DateTime dEpoch, double fMEpochLong, double fMPeriLong, double fMAscNode, double fMIncl, double fMEcc, double fSEpochEclLong, double fSPeriEclLong, double fSEcc, double fMSMA, double fVAngDiam, ref double fMAngDiam)
		{
			double fRho;

			fRho = 0;
			CalcMoonDistance(dDate, dEpoch, fMEpochLong, fMPeriLong, fMAscNode, fMIncl, fMEcc, fSEpochEclLong, fSPeriEclLong, fSEcc, fMSMA, ref fRho);

			fMAngDiam = (fVAngDiam * fMSMA) / fRho;
		}
Share