aaplus - calculating planet altitudes, etc...

[furrysocks2]

In November, I came across aaplus by PJ Naughter. 

Quote

 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.

[furrysocks2]

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.

[johnfosteruk]

Good work sir. Downloading to have a play as we speak.

[furrysocks2]

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.

 

[furrysocks2]

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

Quote

CAAElliptical

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

...

CAAElliptical::Calculate

 

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,
                                                             planetary_details.ApparentGeocentricDeclination,
                                                             planetary_details.ApparentGeocentricDistance,
                                                             obs_lat_long.second, obs_lat_long.first, 38, /* TODO height */
                                                             JD);

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,
                                                                       ra_dec.second,
                                                                       obs_lat_long.first);
  

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
  {
    SUN,
    MERCURY,
    VENUS,
    MARS,
    JUPITER,
    SATURN,
    URANUS,
    NEPTUNE,
    PLUTO,
    UNKNOWN
  };

 

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
{
public:
    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)
    { ... }
  
private:
    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.

Quote

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"

etc...

[furrysocks2]

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

Quote

CAARiseTransitSet

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(
	JD,
	ad1.first, ad1.second,
	ad2.first, ad2.second,
	ad3.first, ad3.second,
	obs.second, obs.first,
	-0.5667);
}

 

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:

Quote

2017/2/13,43.1639
2017/2/20,45.2431
2017/2/27,47.2611
2017/3/6,49.2042
2017/3/13,51.059
2017/3/20,52.8133
2017/3/27,54.4562
2017/4/3,55.9769
2017/4/10,57.3648
2017/4/17,58.6108
2017/4/24,59.7071
2017/5/1,60.6471
2017/5/8,61.4247
2017/5/15,62.0356
2017/5/22,62.4774
2017/5/29,62.7493
2017/6/5,62.8516
2017/6/12,62.7863
2017/6/19,62.5571
2017/6/26,62.1692
2017/7/3,61.6288
2017/7/10,60.9431
etc...

 

The following graph can then be produced:

 

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.

[furrysocks2]

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:

 

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:

 

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.

[furrysocks2]

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):

Quote

CAAPhysicalJupiter

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

...

CAAPhysicalMars

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:

Quote

2017/2/13,36.9284,,4.84026
2017/2/20,36.5782,,4.72615
2017/2/27,35.8596,,4.6183
2017/3/6,34.8268,,4.5167
2017/3/13,33.5336,,4.42108
2017/3/20,32.0227,,4.33109
2017/3/27,30.3296,,4.24659
2017/4/3,28.4908,,4.16757
2017/4/10,26.5343,,4.09382
2017/4/17,24.4751,,4.02505
2017/4/24,22.3384,,3.96111
2017/5/1,20.1495,,3.90198
2017/5/9,17.9379,,3.84756
2017/5/16,15.7365,,3.79758

etc...

Graphing directly gives us:

 

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:

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.

 

[furrysocks2]

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:

Quote

CAAPlanetaryPhenomena

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);
        mars_oppositions.insert(t);
    }
    for (auto opposition_date : mars_oppositions)
    {
        std::cout << formatTime(opposition_date) << "\n";
    }

This gives the following output:

Quote

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.

[furrysocks2]

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.

[furrysocks2]

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.

[furrysocks2]

And Mars in the same way:

[furrysocks2]

Something to incorporate. Not sure how yet, though.

Source: https://www.timeanddate.com/astronomy/different-types-twilight.html

[furrysocks2]

Another thread using aaplus - for a given date, correcting J2000 RA/Dec coordinates for proper motion, precession, nutation and aberration: