How can I calculate the positions of Jupiter and Saturn?
An answer to this question on Stack Overflow.
Question
Jupiter and Saturn have just met up in the sky. Where can I look for them? When will they do it again?
PyEphem seems to be a standard programming tool for performing this kind of calculation, but how can it be employed in this instance?
Answer
Producing Planetary Data (Code)
First, I used pyephem to generate a lot of data about the relative positions of the planets:
#!/usr/bin/env python3
#Running this script produces a CSV file containing information about the
#positions of Jupiter and Saturn
import ephem
from datetime import datetime, timedelta
my_lat = 40.0150 # Boulder, Colorado is a space-y place
my_lng = -105.2705
ephem_tdfmt = '%Y-%m-%d %H:%M:%S' # Used for formatting PyEphem strings
p1 = ephem.Jupiter() # Planet of interest
p2 = ephem.Saturn() # Second planet of interest
observer = ephem.Observer()
observer.lat = str(my_lat) # PyEphem takes lat/lng in strings
observer.lon = str(my_lng)
observer.elev = 1624 # Elevation in metres
# We'll generate a large amount of data, one point for each day,
# starting at the time below. We choose the hour/minute to be a
# time when we're interested in looking at the sky
base_time = datetime(year=2018, month=12, day=21, hour=23, minute=0, second=0, microsecond=0)
# Let's print a CSV to stdout with info about the planets
print("date time Planet az alt sep mag")
for x in range(2000*365): # 365 days a year for 2000 years
calc_time = base_time + timedelta(days=x)
observer.date = ephem.date(calc_time.strftime(ephem_tdfmt))
p1.compute(observer)
p2.compute(observer)
sep = ephem.separation(p1, p2)
print(calc_time, "Jupiter", float(p1.az), float(p1.alt), float(sep), float(p1.mag))
print(calc_time, "Saturn", float(p2.az), float(p2.alt), float(sep), float(p2.mag))
Analyzing and Plotting the Data (Code)
Next, let's use R to crunch the numbers and make some plots
library(directlabels)
library(dplyr)
library(ggplot2)
library(ggrepel)
library(lubridate)
library(reshape2)
library(tidyr)
# Used for numbering the conjunctions
get.groups = function(x) with(rle(x),rep(1:length(lengths),lengths))
# Read the data
dat = read.table("output", header=TRUE)
dat = dat %>% mutate(date=lubridate:::date(date), sep=sep*180/2/pi)
# Plot angular separation for the next 2,000 years
# Since `sep` is common to both planets, we arbitrarily choose one
sep_over_time = dat %>% filter(Planet=="Jupiter")
ggplot(sep_over_time, aes(x=date, y=sep*180/2/pi)) +
geom_line() +
xlab("Date") +
ylab("Separation (degrees)") +
scale_y_log10() +
ggtitle("Angular Separation over time of Jupiter and Saturn (2018-4017)")
ggsave("/z/sep_over_time.png", width=10, height=6)
# Data used to plot the "shape" of a number of future conjunctions
conjunctions_time_normalized = dat %>%
group_by(Planet) %>%
mutate(cnum=get.groups(sep<1.5)) %>%
ungroup() %>%
filter(sep<1.5) %>%
group_by(cnum) %>%
mutate(days=date-first(date))
# Get information about the 2020 conjunction
this_year = conjunctions_time_normalized %>%
filter(year(date)==2020 | year(date)==2021) %>%
ungroup()
# Minimal angular separation of the 2020 conjunction
this_year_sep = min(this_year$sep)
# Plot the "shape" of future conjunctions by angular separation
ggplot(conjunctions_time_normalized %>% filter(Planet=="Jupiter"), aes(x=days, y=sep, group=as.factor(cnum))) +
geom_line(color="blue", alpha=0.3) +
geom_line(data=this_year, aes(x=days, y=sep), color="green") +
geom_hline(yintercept=this_year_sep) +
xlab("Days from first <1.5deg separation") +
ylab("Separation (deg)") +
ggtitle("Angular Separation and Timing of Jupiter-Saturn Conjunctions (2018-4017)")
ggsave("/z/all_conjunctions.png", width=10, height=6)
# Get positions of the planets for the current (2020) conjunction
temp = this_year %>% mutate(az=az*180/2/pi,alt=alt*180/2/pi) %>% mutate(date=ifelse(wday(date)==1,as.character(date),NA))
ggplot(temp, aes(x=az, y=alt, color=Planet, label=date)) +
geom_point() +
geom_text_repel(nudge_y=2) +
xlab("Azimuth") +
ylab("Zenith") +
theme(legend.position="bottom") +
ggtitle("2020 Conjunction of Jupiter and Saturn (16:00 local time)")
ggsave("/z/conjunction_timing.png")
# Get times of next conjunction
print("Times of next conjunctions")
print(dat %>% filter(sep<this_year_sep))
The Current Conjunction
So, here we have the path of Jupiter and Saturn through the sky. Note that if you missed the closest point of separation on the 21st, tonight (the 22nd), and the 23rd will both be good as well. After that, you can enjoy watching them pull away from each other.
When Will It Happen Again?
Although Jupiter and Saturn regularly get close to each other (the period is about 20 years), the chart below shows that it's rare for them to get this close.
In fact, over the next 2,000 years they'll only be this close a handful of times:
date time Planet az alt sep mag
1 2080-03-14 23:00:00 Jupiter 4.565937 -0.31129250 0.05045080 -1.86
2 2080-03-14 23:00:00 Saturn 4.564148 -0.31174162 0.05045080 0.76
3 2417-08-24 23:00:00 Jupiter 4.923358 0.27409384 0.04687376 -1.70
4 2417-08-24 23:00:00 Saturn 4.922710 0.27258733 0.04687376 0.24
5 2874-12-24 23:00:00 Jupiter 3.558560 0.49130434 0.04148634 -1.83
6 2874-12-24 23:00:00 Saturn 3.557792 0.49258307 0.04148634 0.63
9 3271-08-14 23:00:00 Jupiter 4.730337 0.46493667 0.03977666 -1.70
10 3271-08-14 23:00:00 Saturn 4.728929 0.46435314 0.03977666 0.13
13 3669-02-27 23:00:00 Jupiter 4.407206 -0.08765785 0.05079357 -1.73
14 3669-02-27 23:00:00 Saturn 4.408436 -0.08571938 0.05079357 0.58
15 3728-12-29 23:00:00 Jupiter 3.501876 0.55620575 0.04263970 -1.83
16 3728-12-29 23:00:00 Saturn 3.500159 0.55650318 0.04263970 0.61
Are Conjunctions Always Like This?
Nope! As the graph below shows, conjunctions in which Jupiter and Saturn get very close last only a brief time. However, conjunctions in which they get kinda close can last almost a year. Here I've highlighted the current conjunction in bright green and the minimal angular separation with a black line.
Happy Star/Planet-gazing!