Jump to content


aaplus - calculating planet altitudes, etc...

Recommended Posts

In November, I came across aaplus by PJ Naughter. 


 A class framework for Computational Astronomy

AA+ is a C++ implementation for the algorithms as presented in the book "Astronomical Algorithms" by Jean Meeus.


I was playing around and managed to use his library to calculate topocentric horizontal coordinates for the major planets, for a given local time - that means alt-az from my observing location. They were corrected for the speed of light and atmospheric refraction. I got it running on an Arduino for fun and was happy at that and left it.


Last night, in response to a post...

... I went back to my code and tweaked it to produce the following graph, amongst others:


It shows how high Mars gets in the sky between sunset and sunrise (yellow indicating meridian transit) over a period of five years and the error bars give an indication of relative apparent angular diameter, or size. I've read several times that Mars is "due to come back" to UK skies in 2018 but it's clear from this graph that we have to wait until late 2020 for decent views.


In this thread, I'll show how I put this together. It'll be heavy on C++ and templates but I'll try and give explanations within the limits of my understanding. Hopefully, through writing it all down I'll make some improvements over what I've done so far and identify any errors in the above.

Edited by furrysocks2
  • Like 3
Link to comment
Share on other sites

Quite by chance, there was a new release of code today... https://naughter.wordpress.com/2017/02/12/aa-v1-76/.

If you download the zip file, the documentation can be found in a file called AA+.htm which I will refer to as appropriate. It's nicely written and refers to specific formulae and variable names from the book, which I have found in PDF format but have not yet purchased a copy in print. I have not located an online version of the documentation but if you unzip the source to C:\aaplus then you could copy/paste links which I will include alongside bits of my code, for example file:///C:/aaplus/AA+.htm#APIReference.

If you've got a compiler (I used g++ for my first play and Visual Studio while I'm writing up here), then just make sure to include the aaplus/ folder in the include search path and add the cpp files to your build - compile AATest.cpp in the first instance to make sure it's all working then omit it or modify it to have a play yourself. The main include file is "AA+.h", but I'm not intending to write any further on the build process. If you'd like to know how to make it work on Arduino, drop me a PM.


Next post will contain code.

  • Like 1
Link to comment
Share on other sites

While the implementation behind the scenes is entirely aaplus, I've made a few different choices for my code.


I'm using simple pairs for my coordinates:

#include <tuple>

using LatLong = std::pair<double, double>;
using AltAz = std::pair<double, double>;
using RaDec = std::pair<double, double>;

static const LatLong GREENWICH{ 51.4769, -0.0005 }; // 51.4769° N, 0.0005° W

There's no type safety but their convenient and fairly explicit.


I'm using std::chrono for time:

#include <chrono>

namespace chrono = std::chrono;
using namespace std::chrono_literals;

const auto now = chrono::system_clock::now();
const auto one_hour_ago = now - 1h;


I've got a few helper conversion routines for coordinates and time, and some pretty print routines including:

format_ra_dec(x); // RA 22h 51m 23s, Dec -8° 12' 7"
format_alt_az(y); // Alt -9° 19' 3", Az 269° 9' 33"


After that, most of my work went into a class template, Planet:

template <CAAElliptical::EllipticalObject AAPlanetEnum>
class Planet;

using Venus = Planet<CAAElliptical::VENUS>;
using Mars = Planet<CAAElliptical::MARS>;


So basically, I use aaplus coordinates and dates in the implementation but prefer std::pair and std::chrono in my own code. I make use of the various planet-related types and values in aaplus, but get at them through a common implementation, Planet.


Edited by furrysocks2
Link to comment
Share on other sites

At the heart of it, aaplus does all the hard work... see file:///C:/aaplus/AA+.htm#CAAElliptical



This class provides for calculation of the position of an object in an elliptical orbit. This refers to Chapter 33 in the book.




Given a time, t, and a planet, AAPlanetEnum, to calculate for, a call to CAAElliptical::Calcuate returns, amongst other values, Geocentric RA and Dec - relative to the center of the earth.

const auto JD = aa_date(t).Julian();
const auto planetary_details = CAAElliptical::Calculate(JD, AAPlanetEnum, true);

See file:///C:/aaplus/AA+.htm#CAAElliptical_Calculate.


To convert geocentric to topocentric, you need to adjust for where on the globe you are, when, and how far away the object is that you're observing:

const auto topo_coords = CAAParallax::Equatorial2Topocentric(planetary_details.ApparentGeocentricRA,
                                                             obs_lat_long.second, obs_lat_long.first, 38, /* TODO height */

See  file:///C:/aaplus/AA+.htm#CAAParallax_Eq2TopCen.

It's then a simple make_pair() to turn topo_coords into an RaDec value.


To convert RaDec to AltAz, first calculate Greenwich sidereal time and adjust using observer longitude and RA of the planet to get the local hour angle:

const auto apparent_gst = CAASidereal::MeanGreenwichSiderealTime(JD);
const auto local_hour_angle = apparent_gst + obs_lat_long.second/15 - ra_dec.first;

See file:///C:/aaplus/AA+.htm#CAASideral_MeanGrennwichSiderealTime.

aaplus can then do the coordinate transformation for you:

const auto az_alt = CAACoordinateTransformation::Equatorial2Horizontal(local_hour_angle,

See file:///C:/aaplus/AA+.htm#CAACoordinateTransformation_2DCoordinateEquatorial2Horizontal.


And that's pretty much all my Planet template class does. Each planet has a value in the CAAElliptical::EllipticalObject enum:

  enum EllipticalObject


The only thing that differs in the topocentric RaDec and AltAz calculations above is the value of that enum, so the code is common. The class template lets you do this, share the code and parameterise, using the type system, on which enum value you're using - this is what makes the types "Mars", "Venus" and all the other major planets different, but share implementation.

template <CAAElliptical::EllipticalObject AAPlanetEnum>
class Planet
    Planet(chrono::system_clock::time_point t = chrono::system_clock::now())
    : t_(t) {}

    RaDec topocentric_ra_dec(LatLong obs_lat_long)
    { return topocentric_ra_dec(obs_lat_long, t_); }

    static RaDec topocentric_ra_dec(LatLong obs_lat_long, chrono::system_clock::time_point t)
    { ... } 
    AltAz topocentric_alt_az(LatLong obs_lat_long)
    { return topocentric_alt_az(obs_lat_long, t_); }

    static AltAz topocentric_alt_az(LatLong obs_lat_long, chrono::system_clock::time_point t)
    { ... }
    const chrono::system_clock::time_point t_;


Then each planet just becomes a specialisation of this class template, using the aaplus enumeration values:

using Sun = Planet<CAAElliptical::SUN>;
using Mercury = Planet<CAAElliptical::MERCURY>;
using Venus = Planet<CAAElliptical::VENUS>;
using Mars = Planet<CAAElliptical::MARS>;
using Jupiter = Planet<CAAElliptical::JUPITER>;
using Saturn = Planet<CAAElliptical::SATURN>;
using Uranus = Planet<CAAElliptical::URANUS>;
using Neptune = Planet<CAAElliptical::NEPTUNE>;


I can then do something relatively simple like this:

const auto now = chrono::system_clock::now();
for (auto t = now; t < now + 1h; t+=1min)
    const auto alt_az = Mars::topocentric_alt_az(GREENWICH, t);
    std::cout << format_time(t) << ": " << format_alt_az(alt_az) << "\n";

This should give me apparent alt-az coordinates for Mars, as viewed from Greenwich, every minute for the next hour.


2017/2/13 10:39:15: Alt 17° 21' 27", Az 104° 57' 43"
2017/2/13 10:40:15: Alt 17° 30' 29", Az 105° 10' 13"
2017/2/13 10:41:14: Alt 17° 39' 31", Az 105° 22' 44"
2017/2/13 10:42:14: Alt 17° 48' 32", Az 105° 35' 16"
2017/2/13 10:43:14: Alt 17° 57' 33", Az 105° 47' 50"
2017/2/13 10:44:14: Alt 18° 6' 33", Az 106° 0' 24"
2017/2/13 10:45:14: Alt 18° 15' 33", Az 106° 12' 59"


Edited by furrysocks2
Link to comment
Share on other sites

aaplus also does the hard work in calculating rise, meridian transit and set times:



This class provides for calculation of the time of rise, transit and set of a celestial body. This refers to Chapter 15 in the book.

See file:///C:/aaplus/AA+.htm#CAARiseTransitSet.

It requires three RA/Dec values for your object, a day either side of the time you're calculating for:

static CAARiseTransitSetDetails rise_transit_set(LatLong obs,
		chrono::system_clock::time_point t)
  const auto ad1 = topocentric_ra_dec(obs, t - 24h);
  const auto ad2 = topocentric_ra_dec(obs, t);
  const auto ad3 = topocentric_ra_dec(obs, t + 24h);

  const auto JD = aa_date(t).Julian();
  return CAARiseTransitSet::Calculate(
	ad1.first, ad1.second,
	ad2.first, ad2.second,
	ad3.first, ad3.second,
	obs.second, obs.first,


The structure returned contains several values, with times in decimal hour offset from the time for which you're calculating. I convert these back to chrono time_points for my convenience:

	chrono::system_clock::time_point transit()
	{ return t_ + (static_cast<uint32_t>(rise_transit_set().Transit * 3600) * 1s); }


So I can now do this:

	const auto now = chrono::system_clock::now();
	for (auto t = now; t < now + 5*364*24h; t += 7*24h)
		Mars mars(GREENWICH, t);
		const auto transit_time = mars.transit();
		const auto transit_alt_az = mars.topocentric_alt_az(transit_time);
		std::cout << formatDateOnly(transit_time) << "," << transit_alt_az.first << "\n";

Which outputs the following:




The following graph can then be produced:

mars transit.jpg


In summary, this graph shows the altitude at which Mars will transit the meridian as viewed from Greenwich, night and day, for the next five years.

  • Like 1
Link to comment
Share on other sites

Meridian transits occuring in the daytime are no use, so we can do the following:

Sun sun(GREENWICH, t);
Mars mars(GREENWICH, t);
const auto transit_time = mars.transit();
const bool sun_up = sun.set() < sun.rise();
const bool transit_visible_at_night =
	sun_up ? (transit_time > sun.set() && transit_time < sun.rise())
		   : (transit_time < sun.rise() || transit_time > sun.set());

This lets us tell whether a meridian transit occurs between sunset and sunrise.

By adding another column of data for when these occur, we can draw the following graph:

mars transit.jpg


Now, in red is the periods during which a meridian transit occurs between sunset and sunrise. If a meridian transit is visible at night, then this is as high as an object will rise during that night. We're not interested in meridian transits during the day, but it would be interesting to know how high an object rises during the night.

We can do that by selecting the meridian transit altitude if such is visible, otherwise, taking the maximum altitude of the object at either sunrise or sunset.

const auto transit_alt_az = mars.topocentric_alt_az(transit_time);

const auto sunrise_alt_az = mars.topocentric_alt_az(sun.rise());
const auto sunset_alt_az = mars.topocentric_alt_az(sun.set());

const auto max_alt = transit_visible_at_night ? transit_alt_az.first
	: std::max(0.0, std::max(sunrise_alt_az.first, sunset_alt_az.first));


This gives the following graph:

mars transit.jpg


Summary - by calculating sunrise and sunset times, we can differentiate meridian transits that occur during day and night. By replacing daytime meridian transits with the maximum altitude of an object at either sunrise or sunset, whichever is larger, we get a graph of maximum altitude as shown above. This graph now has the same shape as the graph back in the first post, minus the error bars.

Edited by furrysocks2
  • Like 2
Link to comment
Share on other sites

We can get the geocentric distance to any planet (in AU) quite easily:

	static double distance(chrono::system_clock::time_point t)
		const auto JD = aa_date(t).Julian();

		const auto planetary_details = CAAElliptical::Calculate(JD, AAPlanetEnum, true);
		return planetary_details.ApparentGeocentricDistance;

But there is no generic way that I have found to get either planetary radius/diameter or apparent size. If there was, then we could do our own trigonometric calculations for size.


aaplus has two classes (note that there are similar classes for the moon and the sun, but not for the other major planets):



This class provides for calculation of various physical parameters related to the Jupiter. This refers to Chapter 43 in the book.



This class provides for calculation of various physical parameters related to the Mars. This refers to Chapter 42 in the book.

See file:///C:/aaplus/AA+.htm#CAAPhysicalJupiter and file:///C:/aaplus/AA+.htm#CAAPhysicalMars.

Calculating the apparent diameter of Mars is very simple:

CAAPhysicalMars::Calculate(aa_date(t).Julian(), true).d; // d = apparent diameter in arc seconds

Unfortunately, there is not the same calculation performed for Jupiter.


So we're stuck... we can use aaplus to calculate angular size of Mars for us and we can use aaplus to calculate distance from Earth for any planet. If we had physical diameter for all planets, we could do our own calcs for angular size for any planet.


Ignoring the introduction of a generic solution just now, we can use aaplus to get sizes for Mars, and complete the example of how to arrive at the graph in the first post:




Graphing directly gives us:

mars transit.jpg


I created a further column dividing apparent size by the maximum value in the apparent size column, to give a value in the range 0-1, and multiplied this by 4, giving a value in the range 0-4. I'd never plotted error bars on a graph before, but using that 0-4 scale value of relative size gives you the original graph:

mars transit.jpg

Summary - aaplus has a class to calculate the apparent diameter of Mars, but not the other planets. Error bars can be added to the "maximum altitude" graph to give a reasonable indication of relative size.


  • Like 1
Link to comment
Share on other sites

You can interpret approximate opposition dates from the graph above by estimating the maximum apparent size from the error bars, or better finding the maximum value in the data table produced by the code above.

aaplus lets you calculate oppositionss, and other events:



This class provides for the calculation of several configurations involving planets Mercury to Neptune; oppositions and conjunctions with the Sun, greatest elongations, and stations. This refers to Chapter 36 in the book.

See file:///C:/aaplus/AA+.htm#CAAPlanetaryPhenomena.


For Mars, here's how I generated the next few oppositions:

    std::set<double> mars_oppositions;
    for (int year=2017; year<2030; ++year)
        const auto k = CAAPlanetaryPhenomena::K(year, CAAPlanetaryPhenomena::MARS, CAAPlanetaryPhenomena::OPPOSITION);
        const auto t = CAAPlanetaryPhenomena::True(k, CAAPlanetaryPhenomena::MARS, CAAPlanetaryPhenomena::OPPOSITION);
    for (auto opposition_date : mars_oppositions)
        std::cout << formatTime(opposition_date) << "\n";

This gives the following output:


2016-5-22 11:15:56
2018-7-27 05:32:30
2020-10-13 23:11:24
2022-12-8 04:24:46
2025-1-16 01:18:03
2027-2-19 15:13:14
2029-3-25 07:34:42


The enumeration CAAPlanetaryPhenomena::PlanetaryObject is different from CAAElliptical::EllipticalObject, but there is a value for each of the major planets, so I could add a further template parameter to my Planet class template, with a dummy value perhaps for the Sun, which I shoe-horned into being a Planet... this would let me do something like:

std::cout << "Next mars opposition on " << formatDateOnly(Mars::next_opposition());

But I don't currently feel the need.


I could, however, annotate my previous chart with opposition dates... a quick and dirty way might be to add another column with a single max-alt value on the first calculated date after such an event. I loop the set of calculated dates and write a value in the next column when I reach/exceed that date, and remove it from the set. The green dots show opposition.

mars transit.jpg

Edited by furrysocks2
  • Like 1
Link to comment
Share on other sites

That was all I was intending on posting, really.

One final graph, for my location rather than Greenwich I suffer a few degrees of altitude. I've limited my graph to altitudes above 30 degrees and the planets Mars, Jupiter and Saturn.


While the maximum altitude works for the Mars example further up in the thread, for a wide date range like this, it doesn't really add anything and clutters the graph somewhat. Just graphing the transits results in an improved image, in my opinion, which was the first graph I made a couple of nights ago.


I'll try to get Jupiter's angular size plotted up and give a next-few-years graph for it, rather like the Mars ones above.


All this sort of information is available with far less effort than it takes to code, but the aaplus library does so much more than just this, and it's always fun to play around with these things.

Hope you've learned something, please correct any mistakes so that I can learn something too. Thanks for reading.

  • Like 2
Link to comment
Share on other sites

Jupiter - apparent size error bars calculated as follows:

Same value for +ve and -ve, calculated as "((size / min(size)) - 1", ie minimum size represented by error bars having value 0, anything else can be read as "total height of error bars X (in degrees)" means "disc is (X/2)+1 times larger than minimum size". I think. Obscure perhaps, but visually ok... unless they're misleading. Whatever, never mind.

jupiter 5y.jpg

jupiter 10y.jpg

Edited by furrysocks2
Link to comment
Share on other sites

  • 1 month later...

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now
  • Create New...

Important Information

We have placed cookies on your device to help make this website better. You can adjust your cookie settings, otherwise we'll assume you're okay to continue. By using this site, you agree to our Terms of Use.