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