For my EQ3 mount, I had to add planet position calculation to my code... Well, it was not that easy finding useable code online, so I am providing you with the code that I am using here if you ever needed it!

I actually have 2 versions of the code. One uses floating points. And the other one is optimized for smaller memory and devices without FPU (like an arduino) and uses fixed point arithmetics on integers.

Version that uses floating points

//  N = longitude of the ascending node

//  i = inclination to the ecliptic (plane of the Earth's orbit)

//  w = argument of perihelion

//  a = semi-major axis, or mean distance from Sun

//  e = eccentricity (0=circle, 0-1=ellipse, 1=parabola)

//  M = mean anomaly (0 at perihelion; increases uniformly with time)

struct TOrbitalElements { float N1, N2, i1, i2, w1, w2, a1, a2, e1, e2, M1, M2; };

class TOrbitalElementsAt { public:

    float N, i, w, a, e, M; 

    float fmodf360(float a)

    {

        a= fmodf(a, 360.0f); if (a<0.0f) a+= 360.0f;

        return a;

    }

    TOrbitalElementsAt(TOrbitalElements const &o, float d)

    { 

        N= fmodf360(o.N1+o.N2*d); i=fmodf360(o.i1+o.i2*d); w= o.w1+o.w2*d;

        a= o.a1+o.a2*d; e= o.e1+o.e2*d; M= fmodf360(o.M1+o.M2*d);

    }

};

static TOrbitalElements const oe[9]={

    // N                        i                      w                         a                       e                        M

    { 0.0f,       0.0f,         0.0f,     0.0f,        282.9404f,  4.70935E-5f,  1.0f,       0.0f,       0.016709f,  -1.151E-9f,  356.0470f,  0.9856002585f}, // sun

    { 48.3313f,   3.24587E-5f,  7.0047f,  5.00E-8f,    29.1241f,   1.01444E-5f,  0.387098f,  0.0f,       0.205635f,  5.59E-10f,   168.6562f,  4.0923344368f}, // mercury

    { 76.6799f,   2.46590E-5f,  3.3946f,  2.75E-8f,    54.8910f,   1.38374E-5f,  0.723330f,  0.0f,       0.006773f,  1.302E-9f,   48.0052f,   1.6021302244f}, // Venus

    { 49.5574f,   2.11081E-5f,  1.8497f,  -1.78E-8f,   286.5016f,  2.92961E-5f,  1.523688f,  0.0f,       0.093405f,  2.516E-9f,   18.6021f,   0.5240207766f}, // Mars

    { 100.4542f,  2.76854E-5f,  1.3030f,  -1.557E-7f,  273.8777f,  1.64505E-5f,  5.20256f,   0.0f,       0.048498f,  4.469E-9f,   19.8950f,   0.0830853001f}, // Jupiter

    { 113.6634f,  2.38980E-5f,  2.4886f,  -1.081E-7f,  339.3939f,  2.97661E-5f,  9.55475f,   0.0f,       0.055546f,  -9.499E-9f,  316.9670f,  0.0334442282f}, // Saturne

    { 74.0005f,   1.3978E-5f,   0.7733f,  1.9E-8f,     96.6612f,   3.0565E-5f,   19.18171f,  -1.55E-8f,  0.047318f,  7.45E-9f,    142.5905f,  0.011725806f}, // Uranus

    { 131.7806f,  3.0173E-5f,   1.7700f,  -2.55E-7f,   272.8461f,  -6.027E-6f,   30.05826f,  3.313E-8f,  0.008606f,  2.15E-9f,    260.2471f,  0.005995147f} // Nepture

};

static float const MPI= 3.14159265358979323846264338327f;

static float const degToRad = MPI/180.0f;

static float inline dtr(float a) { return float(a)*degToRad; }

// Get sun position compared with earth. unit = AU

static void sunPos(float d, float &sx, float &sy, float &sz, float &ex, float &ey)

{

    float ecl= 23.4393f - 3.563E-7f * d; // ecl angle

    TOrbitalElementsAt sun(oe[0], d);

    float E= sun.M+sun.e*(180.0f/MPI)*sinf(sun.M*degToRad)*(1.0f+sun.e*cosf(sun.M*degToRad)); //ccentric anomaly E from the mean anomaly M and from the eccentricity e (E and M in degrees):

    // compute the Sun's distance r and its true anomaly v from:

    float xv = cosf(E*degToRad)-sun.e;

    float yv = sqrt(1.0f-sun.e*sun.e) * sinf(E*degToRad);

    float v = atan2f( yv, xv );

    float rs = sqrt( xv*xv + yv*yv );

    //compute the Sun's true longitude

    float lonsun = v + sun.w*degToRad;

    // Convert lonsun,r to ecliptic rectangular geocentric coordinates xs,ys:

    ex= sx = rs * cosf(-lonsun);

    ey = -rs * sinf(-lonsun);

    // (since the Sun always is in the ecliptic plane, zs is of course zero). xs,ys is the Sun's position in a coordinate system in the plane of the ecliptic. To convert this to equatorial, rectangular, geocentric coordinates, compute:

    sy= ey * cosf(ecl*degToRad);

    sz= ey * sinf(ecl*degToRad);

}

// get an object position in solar system

static void objPos(float d, int ob, float &x, float &y, float &z)

{

    TOrbitalElementsAt o(oe[ob], d);

    float E= o.M+o.e*(180.0f/MPI)*sinf(o.M*degToRad)*(1.0f+o.e*cosf(o.M*degToRad)); //ccentric anomaly E from the mean anomaly M and from the eccentricity e (E and M in degrees):

    float E1= E;

    for (int i=0; i<5; i++) E= E - ( E - o.e*(180.0f/MPI) * sinf(E*degToRad) - o.M) / ( 1 - o.e * cosf(E*degToRad) );

    float xv = o.a * ( cosf(E*degToRad) - o.e );

    float yv = o.a * ( sqrt(1.0f - o.e*o.e) * sinf(E*degToRad) );

    float v = atan2f( yv, xv );

    float r = sqrt( xv*xv + yv*yv );

    x = r * ( cosf(o.N*degToRad) * cosf(v+o.w*degToRad) - sinf(o.N*degToRad) * sinf(v+o.w*degToRad) * cosf(o.i*degToRad) );

    y = r * ( sinf(o.N*degToRad) * cosf(v+o.w*degToRad) + cosf(o.N*degToRad) * sinf(v+o.w*degToRad) * cosf(o.i*degToRad) );

    z = r * ( sinf(v+o.w*degToRad) * sinf(o.i*degToRad) );

}

static void objPos(float d, int ob, float &x, float &y, float &z, float sex, float sey)

{

    // Sun centered, ecliptic centered coordinates...

    // Need to rotate around earth center by ecliptic!

    // so, add earth coordinate, rotate by ecliptic and stay there...

    float tx, ty, tz;

    objPos(d, ob, tx, ty, tz);

    x= tx+sex, ty= ty+sey;

    float ecl= 23.4393f - 3.563E-7f * d; // ecl angle

    y= ty * cosf(ecl*degToRad) - tz*sinf(ecl*degToRad);

    z= ty * sinf(ecl*degToRad) + tz*cosf(ecl*degToRad);

}

static char const planetNames[8][8]= { "Soleil", "Mercure", "Venus", "Mars", "Jupiter", "Saturne", "Uranus", "Neptune" };

void planetPos(int y, int m, int d, int h, int M, int s, int planet, float ra, float dec) // return planet position. 0 is sun, 1, 2 are mercury, venus, 3-7 are Mars to Neptune

{

    float d= 367*y - (7*(y+((m+9)/12)))/4 + (275*m)/9 + d - 730530L;

    d+= (h+M/60.0f+s/3600.0f)/24.0f;

    sunPos(d, sx, sy, sz, sex, sey);

    float ox, oy, oz; objPos(d, planet, ox, oy, oz, sex, sey);

    ra= atan2f(oy, ox)/degToRad/15.0f, dec= atan2f(oz, sqrtf(ox*ox+oy*oy))/degToRad;

}

Version that uses only integers.

It looks bigger than the other one, but only because it does not rely on an external FPU to do the work. Compile the version above on a CPU with no FPU and it will take more memory that this version.

It is not as precise, but should be ok enough to have the planet in the EP assuming that you are not zooming like crazy!

typedef int32_t Cartesian; // Coordinates in AU in 1+7+24 format (1 bit sign, 7 bit IP, 24 bit FP) for a max value of 127AU which is OK here...

Cartesian Cmul(Cartesian v1, Cartesian v2) { return Cartesian((int64_t(v1)*v2)>>24); }

Cartesian Cmul(Cartesian v1, Cartesian v2, Cartesian v3) { return Cartesian((((int64_t(v1)*v2)>>24)*v3)>>24); }

Cartesian Cdiv(Cartesian v1, Cartesian v2) { return Cartesian((int64_t(v1)<<24)/v2); }

Cartesian sqrt(Cartesian a, Cartesian b, bool sign=false) // return sqrt a*a+/-b*b (substract b if sign = true)

{

    int64_t v= int64_t(b)*b; // do the multiplication in 64 bit mode...

    if (sign) v= -v;

    v+= int64_t(a)*a;

    if (v<=0) return 0;

    // Here, v is a^2+/-b^2 with 48 bit precision...

    uint32_t t= uint32_t(v>>48); uint8_t bits= 0; while (t!=0) { bits++; t>>=1; } // count the number of bits of ip to generate the initial value close to the sqrt(v)

    int64_t x= v>>((bits>>1)+24); // remove the extra 24 bits of precision from the mul and shift by bits/2 bits as first guess on sqrt..

    for (int8_t i=10; --i>=0; ) x= (x+v/x)>>1; // 10 iterations of newton for sqrt

    return int32_t(x);

}

typedef uint16_t Angle; // [0, 65536[ correspond to [0,360deg[

// sinTable[i]=65536*sin(i in 256th of circle. One full circle is 256)

static uint16_t const sinTable[65]={0, 1608, 3215, 4821, 6423, 8022, 9616, 11204, 12785, 14359, 15923, 17479, 19024, 20557, 22078, 23586, 25079, 26557, 28020, 29465, 30893, 32302, 33692, 35061, 36409, 37736, 39039, 40319, 41575, 42806, 44011, 45189, 46340, 47464, 48558, 49624, 50660, 51665, 52639, 53581, 54491, 55368, 56212, 57022, 57797, 58538, 59243, 59913, 60547, 61144, 61705, 62228, 62714, 63162, 63571, 63943, 64276, 64571, 64826, 65043, 65220, 65358, 65457, 65516, 65535};

// sin (a+b) = sin(a)*cos(b) + cos(a)*sin(b)

// if b is small cos(b)=1 and  sin(b)=b

// sin (a+b) = sin(a) + b*cos(a)

// here we choose b to be between 0 and 2pi/256 or 0 and 0.024rad

// we use a table for sin/cos for values for 0 to 2pi with 8 bit accuracy...

// The max error is 0.114% with an average error of 0.018%

static Cartesian sin(Angle a) // a: [0, 65536[ correspond to [0,360deg[. Returns sin(a)*(1<<24)

{

    int32_t s, c;

    uint8_t a2= a>>8; // separte a into a (high part) and a2 (low part)

    // find sin/cos of a2 in table (with 16 bit precision)

    if (a2<0x40) s= sinTable[a2], c= sinTable[64-a2];

    else if (a2<0x80) s= sinTable[128-a2], c= -int32_t(sinTable[a2-64]);

    else if (a2<0xC0) s= -int32_t(sinTable[a2-128]), c= -int32_t(sinTable[192-a2]);

    else s= -int32_t(sinTable[uint8_t(0-a2)]), c= sinTable[a2-192];

    int32_t b= (a&0xff)*25; // b has to be in radiant. a is in 1/2^16 part of 0..2pi interval.

                            // convert a&0xff to a "real number", by multiply by 2pi and divide by 2^16..

                            // Here we will take 18 bit precision. a having 16 bits, and 2*pi*4=25 adding an extra 2 bits

    b= (b*c)>>10;           // multiply by cos with 16 bits for 34 bit precision (with room to spare as our numbers are <<1)... and remove 10 bits to get to our final 24 bits

    return (s<<8)+b;        // return sin(a) with 24 bit precsion + b*cos(a) with also the right precsion..

}

Cartesian cos(Angle v) { return sin(v+0x4000); }

Angle atan2(Cartesian y, Cartesian x)

{ 

    if (x==0) // x=0, can be +/-90 depending on y

    {

        if (y==0) return 0; // error!!!

        if (y>0) return 0x4000;

        return 0xC000;

    }

    // We will calculate atan(x in [0..1]).

    // But it this only correspond to an eight of the circle... so we will need to correct the final answer by potentially negating it and adding a value...

    // the sign of x and y as well as the order of their magnitude indicates the octan that we are in

    // The table changes encodes the 8 possibilities with low order bit indicating an inversion and high order bits a number to add

    uint8_t const changes[8]= { 0x00, 0x41, 0x81, 0x40, 0x01, 0xC0, 0x80, 0xC1 };

    uint8_t addAtEnd= changes[((y>>29)&0b100)+((x>>30)&0b10)+(abs(x)>abs(y)?0:1)]; // Get the right number in the table...

    // now calculate arctan(y/x) or x/y so that v=y/x is in [0..1]

    Cartesian v; if (abs(x)<abs(y)) v= Cdiv(x,y); else v= Cdiv(y,x); if (v<0) v=-v;

    // approximate arctan(x) by a 4th order polynomial centered at x=0.5 (taylor series)...

    // 0.463647609001+0.8*(x-0.5)-0.32*(x-0.5)^2-4.26666666667e-2*(x-0.5)^3+0.1536*(x-0.5)^4

    // But the coefitians here are for a result in radians.. and we want a result in our angle unit, so multiply by 0x8000/pi=10430

    // so the coefs becomes: 4836, 8344, -3338, -445, 1602

    // We will use a ((ax)*x+b)*x... type polynomial evaluation to optimize...

    // so, when it comes to the coefs for x3, and lower, we need a version which is shifted by 24 bits...

    // so we precalculate it for higher precision as X3:−7466355650, X2:−55997667375, X:139994168438, X0: 81134951838

    int64_t X= v-(1L<<23); // X=v-0.5

    int64_t r= ((1602*X-7466355650LL)*X)>>24;

    r= ((r-55997667375LL)*X)>>24;

    r= ((r+139994168438LL)*X)>>24;

    r= r+81134951838LL;

    if (r<0) r=-r;

    uint16_t r2= uint16_t(r>>24);

    if ((addAtEnd&1)!=0) r2= -r2;

    r2+= uint16_t(addAtEnd&0xf0)<<8;

    return r2;

}

// Table of orbital elements

//  N = longitude of the ascending node

//  i = inclination to the ecliptic (plane of the Earth's orbit)

//  w = argument of perihelion

//  a = semi-major axis, or mean distance from Sun

//  e = eccentricity (0=circle, 0-1=ellipse, 1=parabola)

//  M = mean anomaly (0 at perihelion; increases uniformly with time)

struct { uint16_t N; uint16_t N2; uint16_t i; uint16_t w; int16_t w2; uint32_t a; uint32_t e; uint16_t M; int32_t M2; } const oe[8]={

    //N      N2   i     w      w2   a           e         M      M2

    { 0,     0,   0,    51507, 561, 16777216L,  280330L,  64816, 11758669L }, // sun. or the earth with negative coordinates...

    { 8798,  387, 1275, 5301,  121, 6494427L,   3449982L, 30702, 48823452L }, // mercury

    { 13959, 294, 617,  9992,  165, 12135464L,  113632L,  8739,  19114158L }, // venus

    { 9021,  251, 336,  52156, 349, 25563242L,  1567075L, 3386,  6251811L },  // mars

    { 18287, 330, 237,  49857, 196, 87284472L,  813661L,  3621,  991246L },   // jupiter

    { 20691, 285, 453,  61784, 355, 160302112L, 931907L,  57702, 399005L },   // saturne

    { 13471, 166, 140,  17596, 364, 321815680L, 793864L,  25957, 139894L },   // uranus

    { 23989, 359, 322,  49670, -71, 504293920L, 144384L,  47376, 71524L },    // neptune

};

static void objPos(int32_t d, int8_t obj, Cartesian &x, Cartesian &y, Cartesian &z)

{

    // get an object position in Sun centered, ecliptic centered coordinates (ecept for the earth where it is the reverse)

    Angle N= oe[obj].N+((oe[obj].N2*d)>>16), w= oe[obj].w+((oe[obj].w2*d)>>16), M= oe[obj].M+((oe[obj].M2*d)>>16);

    Angle E= M+Angle( (Cmul( oe[obj].e, sin(M), (1L<<24)+Cmul(oe[obj].e, cos(M)))*10430LL)>>24); // excentric anomaly E from the mean anomaly M and from the eccentricity e:

    if (obj!=0)  // not needed for the sun...

      for (int8_t i=5; --i>=0;) 

        E= E - Angle( ((int64_t(E-M)<<24) - Cmul(oe[obj].e, sin(E))*10430LL) / ((1L<<24) - Cmul(oe[obj].e, cos(E)) ) ); // not used for sun...

    Cartesian xv= Cmul(oe[obj].a, cos(E)-oe[obj].e);

    Cartesian yv= Cmul(oe[obj].a, sqrt(1L<<24, oe[obj].e, true), sin(E));

    Angle v= atan2(yv, xv);

    Cartesian r= sqrt(xv,yv);

    Angle lon = v + w; //compute the longitude. i and N are actually 0 for sun... but we are size optimizing here...

    Cartesian clon= cos(lon), slon= sin(lon);

    Cartesian sloncosi= Cmul(slon, cos(oe[obj].i));

    x = Cmul(r, (int64_t(cos(N))*clon - int64_t(sin(N))*sloncosi )>>24);

    y = Cmul(r, (int64_t(sin(N))*clon + int64_t(cos(N))*sloncosi )>>24);

    z = Cmul(r, slon, sin(oe[obj].i)); 

}

static char planetNames[8][8]= { "Soleil", "Mercure", "Venus", "Mars", "Jupiter", "Saturne", "Uranus", "Neptune" };

int8_t year=23, month=1, day=1; // Parameters for planetPos

void planetPos(int p, int32_t &ra, int32_t &dec) // return planet position. 0 is sun, 1, 2 are mercury, venus, 3-7 are Mars to Neptune

{

    int32_t d= 367*year - (7*((year+2000L)+((month+9)/12)))/4 + (275*month)/9 + day - (730530L-367*2000L);

    Cartesian sex, sey, ox, oy, oz, t;

    objPos(d, 0, sex, sey, t);  // earth position (sun is -earth)

    objPos(d, p, ox, oy, oz);   // object position

    ox= ox+sex, oy= oy+sey;    // object position compared with earth... sunz=0...

    // Need to rotate around earth center by ecliptic!

    int64_t const sinecl= 6673596LL, cosecl= 15392794LL; // sin and cos of ecliptic

    t= (oy*cosecl - oz*sinecl)>>24; oz= (oy*sinecl + oz*cosecl)>>24; oy= t;    // ecliptic rotation...

    int32_t raf= atan2(oy, ox), decf= atan2(oz, sqrt(ox,oy));  // cartesian to spherical transform...

    ra= (raf*5400)>>12; dec= (decf*5062)>>8; // to my fixed point system, which has ra in hours*3600 and dec in deg*3600

    if (dec>180*3600L) dec= dec-360*3600L; // modulo...

}