Home

Photography

Silver

Computers

Contact

Barry Thomas

Planetary Motion Source Code

This code is from my demo Android app offering a tutorial in coding Keppler's first two laws in Java for Android. The snippet here is intended to be read with the truly excellent work of Jean Meeus in his book, which I strongly recommend to anyone with an interest in astronomy and computing.

Feel free to use this code for any purpose. I don't guarantee it works. The code isn't written for speed or efficiency, just so it works and you can see what is going on. Feel free to improve it.

Keppler

Mr Keppler. What a guy. Also, Jean Meeus. I thank you both. In Meeus's book, chapter 30, he describes the process of "calculating the position of a body (planet, minor planet or periodic comet) on its elliptical orbit around the Sun at a given instant ... from the orbital elements of the body".

My notes here do not teach you all about this topic - you go and do the work yourself, following Mr Meeus. My code here is just one way of calculating the position of the body, according to Meeus's description.

I'm assuming you know how to create an Android app, know Java and have a copy of his book.

This is my method for drawing the planet at the 'next' position in the orbit. Think of the steps like stepping round a circle, one degree at a time, on a circle which has the same period as the planet you are trying to track. Outside of this method I use a global double as the step counter - called dTime, which contains a number of degrees of rotation.

The key parameters passed to the method are, dEccentricty, dScalar (a scaling factor so the orbit all fits on the display), dYear (the duration of the orbit in Earth years) and to orient the orbit so that perihelion is at the right place on the dial, so to speak, dLongPeri - the Longitude of Perihelion.

Read the comments.

public void drawPlanet (double dEccentricity, double dScalar, double dYear, Canvas canvas, Paint paint, 
    		String sName, Bitmap bmp, double dLongPeri)
        {
        double dE, dr, dv, dSatX, dSatY, dSatXCorrected, dSatYCorrected;
        float fX, fY;
        int iSunXOffset = getWidth() / 2;
        int iSunYOffset = getHeight() / 2;
        
        // get the value of E from the angle travelled in this 'tick'
	    
	    dE = getE (dTime * (1 / dYear), dEccentricity);
	    
	    // get r: the length of 'radius' vector
	    
	    dr = getRfromE (dE, dEccentricity, dScalar);
	    
	    // calculate v - the true anomaly
	    
	    dv = 2 * Math.atan (
	    		Math.sqrt((1 + dEccentricity) / (1 - dEccentricity))
	    		*
	    		Math.tan(dE / 2)
	    		); 
	    
	    // get X and Y coords based on the origin
	    
	    dSatX = dr / Math.sin(Math.PI / 2) * Math.sin(dv);
	    dSatY = Math.sin((Math.PI / 2) - dv) * (dSatX / Math.sin(dv));
	    
	    // now correct for Longitude of Perihelion for this planet
	    
	    dSatXCorrected = dSatX * (float)Math.cos (Math.toRadians(dLongPeri)) - 
            dSatY * (float)Math.sin(Math.toRadians(dLongPeri));
	    dSatYCorrected = dSatX * (float)Math.sin (Math.toRadians(dLongPeri)) + 
            dSatY * (float)Math.cos(Math.toRadians(dLongPeri));
	    
	    // offset the origin to nearer the centre of the display
	    
	    fX = (float)dSatXCorrected + (float)iSunXOffset;
	    fY = (float)dSatYCorrected + (float)iSunYOffset;
	    
	    if (bDrawOrbits)
            {
            // draw the path of the orbit travelled
            paint.setColor(Color.WHITE);
            paint.setStyle(Paint.Style.STROKE);
            paint.setAntiAlias(true);
        
            // get the size of the rect which encloses the elliptical orbit
        
            dE = getE (0.0, dEccentricity);
            dr = getRfromE (dE, dEccentricity, dScalar);
            rectOval.bottom = (float)dr;
            dE = getE (180.0, dEccentricity);
            dr = getRfromE (dE, dEccentricity, dScalar);
            rectOval.top = (float)(0 - dr);
        
            // calculate minor axis from major axis and eccentricity
            // http://www.1728.org/ellipse.htm

            double dMajor = rectOval.bottom - rectOval.top;
            double dMinor = Math.sqrt(1 - (dEccentricity * dEccentricity)) * dMajor;
        
            rectOval.left = 0 - (float)(dMinor / 2);
            rectOval.right = (float)(dMinor / 2);

            rectOval.left += (float)iSunXOffset;
            rectOval.right += (float)iSunXOffset;
            rectOval.top += (float)iSunYOffset;
            rectOval.bottom += (float)iSunYOffset;
       
            // now correct for Longitude of Perihelion for this orbit's path

            canvas.save();
            	canvas.rotate((float)dLongPeri, (float)iSunXOffset, (float)iSunYOffset);
            	canvas.drawOval(rectOval, paint);
            canvas.restore();
            }
	    
	    int iBitmapHeight = bmp.getHeight();
		
	    canvas.drawBitmap(bmp, fX - (iBitmapHeight / 2), fY - (iBitmapHeight / 2), null);

	    // draw planet label
	    
		myPaint.setColor(Color.WHITE);
		paint.setTextSize(30);
		canvas.drawText(sName, fX+20, fY-20, paint);
        }

The method above calls two further methods which provide values of E (the mean anomaly) and r, the length of the vector at the end of which the planet is found. Check out the equation numbers and line numbers in Meeus's Basic program on page 206:

    public double getE (double dTime, double dEccentricity)
        {
    	// we are passed the degree count in degrees (duh) 
        // and the eccentricity value
	    // the method returns E

    	double dM1, dD, dE0, dE = 0; // return value E = the mean anomaly
    	double dM; // local value of M in radians
    	      
    	dM = Math.toRadians (dTime);
    	
    	int iSign = 1;
    	
    	if (dM > 0) iSign = 1; else iSign = -1;
    	
    	dM = Math.abs(dM) / (2 * Math.PI); // Meeus, p 206, line 110
    	dM = (dM - (long)dM) * (2 * Math.PI) * iSign; // line 120
    	if (dM < 0)
    		dM = dM + (2 * Math.PI); // line 130
    	iSign = 1;
    	if (dM > Math.PI) iSign = -1; // line 150
    	if (dM > Math.PI) dM = 2 * Math.PI - dM; // line 160
    	
    	dE0 = Math.PI / 2; // line 170
    	dD = Math.PI / 4; // line 170
    	
        for (int i = 0; i < 33; i++) // line 180 
            {
        	dM1 = dE0 - dEccentricity * Math.sin(dE0); // line 190
        	dE0 = dE0 + dD * Math.signum((float)(dM - dM1));
        	dD = dD / 2; 
            }
        
        dE = dE0 * iSign;
        
        return dE;
        }

And once you have E, you need the length of the vector r:

public double getRfromE (double dE, double dEccentricty, double dScalar)
        {
    	return Math.min(getWidth(), getHeight()) / 2 * dScalar * (1 - (dEccentricty * Math.cos(dE)));
        }

If you would like the entire source code of my app, go here.