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

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));
		}
Share