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