Finding the position of a planet is a rather involved calculation. mybloggingplanet.com
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.
public static void CalcPlanetCoord(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 fRA, ref double fDecl)
{
double fD, fPN, fPM, fPL, fPV, fPR;
double fEN, fEM, fEL, fEV, fER;
double fPsi, fX, fY, fA, fL1, fR1, fLambda, fBeta, 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;
fBeta = Math.Atan((fR1 * Math.Tan(Trig.DegToRad(fPsi))* Math.Sin(Trig.DegToRad(fLambda - fL1)))/(fER * Math.Sin(Trig.DegToRad(fL1 - fEL))));
fBeta = Trig.RadToDeg(fBeta);
}
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);
fBeta = Math.Atan((fR1 * Math.Tan(Trig.DegToRad(fPsi))* Math.Sin(Trig.DegToRad(fLambda - fL1)))/(fER * Math.Sin(Trig.DegToRad(fL1 - fEL))));
fBeta = Trig.RadToDeg(fBeta);
}
UraniaCoord.ConvEclToEqu(23.441884,fLambda, fBeta, ref fRA, ref fDecl);
}
Now, this calculation requires quite a bit of planetary data to find the positions, which I put into a class as constants.
public class UraniaPlanetData
{
public int MERCURY = 0;
public int VENUS = 1;
public int EARTH = 2;
public int MARS = 3;
public int JUPITER = 4;
public int SATURN = 5;
public int URANUS = 6;
public int NEPTUNE = 7;
public int PLUTO = 8;
public DateTime dEpoch = new DateTime(1990, 1, 1, 0, 0, 0, 0);
public double[] fPeriod = new double[9];
public double[] fEpochLong = new double[9];
public double[] fPerihelionLong = new double[9];
public double[] fEcc = new double[9];
public double[] fSMA = new double[9];
public double[] fIncl = new double[9];
public double[] fAscNode = new double[9];
public double[] fAngDiam = new double[9];
public double[] fVMag = new double[9];
public string[] sName = new string[9];
public double[] fRadius = new double[9];
public double[] fMass = new double[9];
public double[] fEscapeVelocity = new double[9];
public double[] fAlbedo = new double[9];
public double[] fOblate = new double[9];
public double[] fAvgOrbitVel = new double[9];
public double[] fRotationPeriod = new double[9];
public double[] fAtmPressure = new double[9];
public double[] fTilt = new double[9];
public double[] fGravity = new double[9];
public double[] fDensity = new double[9];
public UraniaPlanetData()
{
sName[0] = "Mercury";
sName[1] = "Venus";
sName[2] = "Earth";
sName[3] = "Mars";
sName[4] = "Jupiter";
sName[5] = "Saturn";
sName[6] = "Uranus";
sName[7] = "Neptune";
sName[8] = "Pluto";
//Tropical years
fPeriod[0] = 0.240852;
fPeriod[1] = 0.615211;
fPeriod[2] = 1.00004;
fPeriod[3] = 1.880932;
fPeriod[4] = 11.863075;
fPeriod[5] = 29.471362;
fPeriod[6] = 84.039492;
fPeriod[7] = 164.79246;
fPeriod[8] = 246.77027;
//Degrees
fEpochLong[0] = 60.750646;
fEpochLong[1] = 88.455855;
fEpochLong[2] = 99.403308;
fEpochLong[3] = 240.739474;
fEpochLong[4] = 90.638185;
fEpochLong[5] = 287.690033;
fEpochLong[6] = 271.063148;
fEpochLong[7] = 282.349556;
fEpochLong[8] = 221.4127;
//Degrees
fPerihelionLong[0] = 77.299833;
fPerihelionLong[1] = 131.430236;
fPerihelionLong[2] = 102.768413;
fPerihelionLong[3] = 335.874939;
fPerihelionLong[4] = 14.170747;
fPerihelionLong[5] = 92.861407;
fPerihelionLong[6] = 172.884833;
fPerihelionLong[7] = 48.009758;
fPerihelionLong[8] = 224.133;
fEcc[0] = 0.205633;
fEcc[1] = 0.006778;
fEcc[2] = 0.016713;
fEcc[3] = 0.093396;
fEcc[4] = 0.048482;
fEcc[5] = 0.055581;
fEcc[6] = 0.046321;
fEcc[7] = 0.009003;
fEcc[8] = 0.24624;
//AU
fSMA[0] = 0.387099;
fSMA[1] = 0.723332;
fSMA[2] = 1.000000;
fSMA[3] = 1.523688;
fSMA[4] = 5.202561;
fSMA[5] = 9.554747;
fSMA[6] = 19.21814;
fSMA[7] = 30.109570;
fSMA[8] = 39.3414;
//Degrees
fIncl[0] = 7.004540;
fIncl[1] = 3.394535;
fIncl[2] = 0;
fIncl[3] = 1.849736;
fIncl[4] = 1.303613;
fIncl[5] = 2.488980;
fIncl[6] = 0.773059;
fIncl[7] = 1.770646;
fIncl[8] = 17.1420;
//Degrees
fAscNode[0] = 48.212740;
fAscNode[1] = 76.589820;
fAscNode[2] = 0;
fAscNode[3] = 49.480308;
fAscNode[4] = 100.353142;
fAscNode[5] = 113.576139;
fAscNode[6] = 73.926961;
fAscNode[7] = 131.670599;
fAscNode[8] = 110.144;
//Arcsec at 1AU
fAngDiam[0] = 6.74 / 3600.0;
fAngDiam[1] = 16.92 / 3600.0;
fAngDiam[2] = 17.0 / 3600.0;
fAngDiam[3] = 9.36 / 3600.0;
fAngDiam[4] = 196.74 / 3600.0;
fAngDiam[5] = 165.60 / 3600.0;
fAngDiam[6] = 65.80 / 3600.0;
fAngDiam[7] = 62.20 / 3600.0;
fAngDiam[8] = 8.20 / 3600.0;
//at 1AU
fVMag[0] = -0.42;
fVMag[1] = -4.40;
fVMag[2] = 0;
fVMag[3] = -1.52;
fVMag[4] = -9.40;
fVMag[5] = -8.88;
fVMag[6] = -7.19;
fVMag[7] = -6.87;
fVMag[8] = -1.00;
//km
fRadius[0] = 2439.7;
fRadius[1] = 6051.8;
fRadius[2] = 6378.14;
fRadius[3] = 3397.2;
fRadius[4] = 71492.0;
fRadius[5] = 60268.0;
fRadius[6] = 25559.0;
fRadius[7] = 24746.0;
fRadius[8] = 1137.0;
//degrees
fTilt[0] = 0.0;
fTilt[1] = 177.36;
fTilt[2] = 23.441884;
fTilt[3] = 25.19;
fTilt[4] = 3.13;
fTilt[5] = 25.33;
fTilt[6] = 97.86;
fTilt[7] = 28.31;
fTilt[8] = 122.52;
//hours
fRotationPeriod[0] = 58.6462;
fRotationPeriod[1] = 243.0187;
fRotationPeriod[2] = 23.9345;
fRotationPeriod[3] = 24.6229;
fRotationPeriod[4] = 9.925;
fRotationPeriod[5] = 10.23306;
fRotationPeriod[6] = 17.9;
fRotationPeriod[7] = 16.10833;
fRotationPeriod[8] = 6.154861;
//km/s
fAvgOrbitVel[0] = 47.88;
fAvgOrbitVel[1] = 35.02;
fAvgOrbitVel[2] = 29.79;
fAvgOrbitVel[3] = 24.13;
fAvgOrbitVel[4] = 13.07;
fAvgOrbitVel[5] = 9.67;
fAvgOrbitVel[6] = 6.81;
fAvgOrbitVel[7] = 5.45;
fAvgOrbitVel[8] = 4.74;
//kg
fMass[0] = 3.303E+23;
fMass[1] = 4.869E+24;
fMass[2] = 5.976E+24;
fMass[3] = 6.421E+23;
fMass[4] = 1.900E+27;
fMass[5] = 5.688E+26;
fMass[6] = 8.69E+25;
fMass[7] = 1.024E+26;
fMass[8] = 1.27E+22;
//km/s
fEscapeVelocity[0] = 4.25;
fEscapeVelocity[1] = 10.36;
fEscapeVelocity[2] = 11.18;
fEscapeVelocity[3] = 5.02;
fEscapeVelocity[4] = 59.56;
fEscapeVelocity[5] = 35.49;
fEscapeVelocity[6] = 21.3;
fEscapeVelocity[7] = 23.5;
fEscapeVelocity[8] = 1.22;
//kg/m^3
fDensity[0] = 5427.0;
fDensity[1] = 5243.0;
fDensity[2] = 5515.0;
fDensity[3] = 3933.0;
fDensity[4] = 1326.0;
fDensity[5] = 687.0;
fDensity[6] = 1290.0;
fDensity[7] = 1640.0;
fDensity[8] = 2000.0;
//g
fGravity[0] = 0.378;
fGravity[1] = 0.907;
fGravity[2] = 1.000;
fGravity[3] = 0.377;
fGravity[4] = 23.12;
fGravity[5] = 0.916;
fGravity[6] = 0.79;
fGravity[7] = 0.378;
fGravity[8] = 0.04;
fAlbedo[0] = 0.1;
fAlbedo[1] = 0.65;
fAlbedo[2] = 0.37;
fAlbedo[3] = 0.15;
fAlbedo[4] = 0.52;
fAlbedo[5] = 0.47;
fAlbedo[6] = 0.51;
fAlbedo[7] = 0.41;
fAlbedo[8] = 0.3;
//atm
fAtmPressure[0] = 0.0;
fAtmPressure[1] = 90.82;
fAtmPressure[2] = 1.0;
fAtmPressure[3] = 0.007;
fAtmPressure[4] = 0.7;
fAtmPressure[5] = 1.4;
fAtmPressure[6] = 1.2;
fAtmPressure[7] = 1.0;
fAtmPressure[8] = 0.00001;
fOblate[0] = 0.0;
fOblate[1] = 0.0;
fOblate[2] = 0.0034;
fOblate[3] = 0.009;
fOblate[4] = 0.0637;
fOblate[5] = 0.102;
fOblate[6] = 0.0;
fOblate[7] = 0.027;
fOblate[8] = 0.0;
}
}
Comments