The R Markdown code is available on GitHub

Dependencies

library(psych)      # For describe
library(lubridate)
library(data.table) # For fread
library(sphereplot)
library(tidyverse)
library(foreach)
library(plotly)

Functions

Function for aggregation

aggnrmean <- function(df, n = 5, FUN = mean) {
  aggregate(df,
            by = list(gl(ceiling(nrow(df)/n), n)[1:nrow(df)]),
            FUN = FUN)[-1]
}

Download Data

Downloading all 1 hour CPI data and combine into one

path <- 'https://spdf.gsfc.nasa.gov/pub/data/pioneer/pioneer10/particle/cpi/ip_1hour_ascii/p10cp_hr'
yseq <- seq(1972, 1992)
rcpi <- data.frame()      # Empty data frame for result
foreach(y = yseq) %do% {
  url <- paste(path, as.character(y), '.asc', sep = '')
  ycpi <- fread(url)
  rcpi <- rbind(rcpi, ycpi)
}

Download mission Trajectory

rtraj <- fread("https://spdf.gsfc.nasa.gov/pub/data/pioneer/pioneer10/traj/ip_project/p10tjall.asc")

Convert date time

cpi <- rcpi %>% mutate(ts = date_decimal(V1) + days(V2) + hours(V3))
traj <- rtraj %>% mutate(ts = date_decimal(V1) + as.duration(86400 * V2))

Select and rename

cpi <- cpi %>% select(
  ts,
  RID2P   = V4,   # ID-2 rate for 11-20 MeV protons [cps]
  RID2HE  = V5,   # ID-2 rate for 11-20 MeV/nucleon helium [cps]
  RID3P   = V6,   # ID-3 rate for 20-24 MeV protons [cps]
  RID3HE  = V7,   # ID-3 rate for 20-24 MeV/nucleon helium [cps]
  RID4P   = V8,   # ID-4 rate for 24-29 MeV protons [cps]
  RID4HE  = V9,   # ID-4 rate for 24-29 MeV/nucleon helium [cps]
  RID5P   = V10,  # ID-5 rate for 29-67 MeV protons [cps]
  RID5HE  = V11,  # ID-5 rate for 29-67 MeV/nucleon helium [cps]
  RID5E1  = V12,  # ID-5 rate for 7-17 MeV electrons [cps]
  RID5E2  = V13,  # ID-5 rate for 2 x minimum-ionizing [cps]
  RID7    = V14,  # ID-7 + ID-13 integral rate for ions at E > 67 MeV/nucleon [cps]
  RID7ZG5 = V15)  # ID-7 integral rate for Z > 5 ions at E > 67 MeV/nucleon [cps]
traj <- traj %>% select(
  ts,
  HRANGP = V3,  # Distance from Sun to spacecraft in km
  SECLAT = V4,  # Solar ecliptic latitude and
  SECLON = V5,  # longitude of spacecraft with respect to true-of-date Ecliptic
  HELLAT = V6,  # Heliographic latitude and
  HELLON = V7,  # longitude of spacecraft
  HILLON = V8,  # Heliographic inertial longitude of spacecraft with respect to
                # direction of zero heliographic longitude on 1 Jan. 1854 at 1200 UT
  REARSC = V9)  # Distance from Earth to spacecraft in AU
                # (used to calculate Earth Received Time - UT)

Replace fill in value with NA

cpi <- na_if(cpi, 1e31)

Time range

CPI data range

summary(cpi$ts)
##                  Min.               1st Qu.                Median 
## "1972-01-02 00:00:00" "1977-04-02 17:45:00" "1982-07-03 11:30:00" 
##                  Mean               3rd Qu.                  Max. 
## "1982-07-03 11:30:00" "1987-10-03 05:15:00" "1993-01-01 23:00:00"

Trajectory data range

summary(traj$ts)
##                  Min.               1st Qu.                Median 
## "1972-03-04 02:09:36" "1976-11-01 06:00:00" "1981-07-04 12:00:00" 
##                  Mean               3rd Qu.                  Max. 
## "1981-07-02 14:32:59" "1986-03-03 18:00:00" "2000-01-02 00:00:00"

Determine Date Range

cpimints <- min(cpi$ts)
cpimaxts <- max(cpi$ts)
trajmints <- min(traj$ts)
trajmaxts <- max(traj$ts)
mints <- max(cpimints, trajmints)
maxts <- min(cpimaxts, trajmaxts)

Filter all data by date range

cpi <- cpi %>% filter(ts >= mints & ts <= maxts)
traj <- traj %>% filter(ts >= mints & ts <= maxts)

Describe data

describe(cpi[,-1], skew = FALSE)
##         vars     n mean   sd min     max   range   se
## RID2P      1 78980 0.28 7.38   0  610.00  610.00 0.03
## RID2HE     2 78980 0.20 5.44   0  391.00  391.00 0.02
## RID3P      3 78980 0.03 0.55   0   26.30   26.30 0.00
## RID3HE     4 78980 0.00 0.22   0   49.50   49.50 0.00
## RID4P      5 78980 0.01 0.25   0   25.90   25.90 0.00
## RID4HE     6 78980 0.00 0.08   0   14.30   14.30 0.00
## RID5P      7 78980 0.04 2.52   0  326.00  326.00 0.01
## RID5HE     8 78980 0.01 0.83   0  146.00  146.00 0.00
## RID5E1     9 78980 0.02 0.38   0   24.80   24.80 0.00
## RID5E2    10 78980 0.03 1.48   0  133.00  133.00 0.01
## RID7      11 74484 0.65 8.11   0 1980.00 1980.00 0.03
## RID7ZG5   12 74484 0.00 0.01   0    1.17    1.17 0.00

Calculate averages

cpi <- cpi %>% mutate(
  RIDMP = (RID2P + RID3P + RID4P + RID5P) / 4,
  RIDMHE = (RID2HE + RID3HE + RID4HE + RID5HE) / 4,
  RIDME = (RID5E1 + RID5E2) / 2,
  RIDM7 = (RID7 + RID7ZG5) / 2)

Aggregate by Mean

mcpi <- aggnrmean(cpi, 10)

Describe aggregated data

describe(mcpi[,-1], skew = FALSE)
##         vars    n mean    sd min    max  range   se
## RID2P      1 2511 0.76 10.87   0 283.60 283.60 0.22
## RID2HE     2 2511 0.49  7.80   0 262.90 262.90 0.16
## RID3P      3 2511 0.06  0.80   0  18.11  18.11 0.02
## RID3HE     4 2511 0.01  0.11   0   2.96   2.96 0.00
## RID4P      5 2511 0.02  0.26   0   6.58   6.58 0.01
## RID4HE     6 2511 0.00  0.06   0   3.04   3.04 0.00
## RID5P      7 2511 0.06  1.81   0  84.28  84.28 0.04
## RID5HE     8 2511 0.02  0.77   0  38.50  38.50 0.02
## RID5E1     9 2511 0.05  0.57   0  11.90  11.90 0.01
## RID5E2    10 2511 0.06  1.45   0  57.09  57.09 0.03
## RID7      11 2353 0.56  0.22   0   7.10   7.10 0.00
## RID7ZG5   12 2353 0.00  0.00   0   0.05   0.05 0.00
## RIDMP     13 2511 0.22  3.15   0  77.67  77.67 0.06
## RIDMHE    14 2511 0.13  2.03   0  66.51  66.51 0.04
## RIDME     15 2511 0.05  0.86   0  30.30  30.30 0.02
## RIDM7     16 2353 0.28  0.11   0   3.55   3.55 0.00

Transforms 3D spherical coordinates to Cartesian coordinates

hcart <- sph2car(traj$SECLON, traj$SECLAT, traj$HRANGP)
traj <- traj %>% mutate(x = hcart[,1], y = hcart[,2], z = hcart[,3])

Convert to data tables

setDT(mcpi)
setDT(traj)

Combine data

Combine CPI and Trajectory data into one set

setkey(mcpi, ts)
setkey(traj, ts)
cpitraj <- traj[mcpi, roll = "nearest"]

Plots

Distance from sun

ggplot(data = traj) +
  geom_point(aes(x = ts, y = HRANGP)) +
  labs(title = "Pioneer 10 Trajectory", x = "Time", y = "Distance from Sun (AU)") +
  theme_light()

Protons

pcpi <- reshape2::melt(mcpi %>% select(time = ts, "11-20 MeV" = RID2P, "20-24 MeV" = RID3P, "24-29 MeV" = RID4P, "29-67 MeV" = RID5P), id = "time") %>% filter(!is.na(value))
ggplot(data = pcpi) +
# geom_line(aes(x = time, y = value, color = variable)) +
  geom_point(aes(x = time, y = value, color = variable)) +
  scale_y_continuous(trans = 'log10') +
  labs(title = "Pioneer 10 CPI Protons", x = "Time", y = "Rate (cps)") +
  theme_light()

pcpi <- reshape2::melt(cpitraj %>% select(distance = HRANGP, "11-20 MeV" = RID2P, "20-24 MeV" = RID3P, "24-29 MeV" = RID4P, "29-67 MeV" = RID5P), id = "distance") %>% filter(!is.na(value))
ggplot(data = pcpi) +
# geom_line(aes(x = distance, y = value, color = variable)) +
  geom_point(aes(x = distance, y = value, color = variable)) +
  scale_y_continuous(trans = 'log10') +
  labs(title = "Pioneer 10 CPI Protons", x = "Distance from Sun (AU)", y = "Rate (cps)") +
  theme_light()

Nucleon helium

pcpi <- reshape2::melt(mcpi %>% select(time = ts, "11-20 MeV" = RID2HE, "20-24 MeV" = RID3HE, "24-29 MeV" = RID4HE, "29-67 MeV" = RID5HE), id = "time") %>% filter(!is.na(value))
ggplot(data = pcpi) +
# geom_line(aes(x = time, y = value, color = variable)) +
  geom_point(aes(x = time, y = value, color = variable)) +
  scale_y_continuous(trans = 'log10') +
  labs(title = "Pioneer 10 CPI nucleon He", x = "Time", y = "Rate (cps)") +
  theme_light()

pcpi <- reshape2::melt(cpitraj %>% select(distance = HRANGP, "11-20 MeV" = RID2HE, "20-24 MeV" = RID3HE, "24-29 MeV" = RID4HE, "29-67 MeV" = RID5HE), id = "distance") %>% filter(!is.na(value))
ggplot(data = pcpi) +
# geom_line(aes(x = distance, y = value, color = variable)) +
  geom_point(aes(x = distance, y = value, color = variable)) +
  scale_y_continuous(trans = 'log10') +
  labs(title = "Pioneer 10 CPI nucleon He", x = "Distance from Sun (AU)", y = "Rate (cps)") +
  theme_light()

Electrons and minimum-ionizing rate

pcpi <- reshape2::melt(mcpi %>% select(time = ts, "7-17 MeV Electrons" = RID5E1, "Double minimum-ionizing" = RID5E2), id = "time") %>% filter(!is.na(value))
ggplot(data = pcpi) +
# geom_line(aes(x = time, y = value, color = variable)) +
  geom_point(aes(x = time, y = value, color = variable)) +
  scale_y_continuous(trans = 'log10') +
  labs(title = "Pioneer 10 CPI Electrons and minimum ionizing", x = "Time", y = "Rate (cps)") +
  theme_light()

pcpi <- reshape2::melt(cpitraj %>% select(distance = HRANGP, "7-17 MeV Electrons" = RID5E1, "Double minimum-ionizing" = RID5E2), id = "distance") %>% filter(!is.na(value))
ggplot(data = pcpi) +
# geom_line(aes(x = distance, y = value, color = variable)) +
  geom_point(aes(x = distance, y = value, color = variable)) +
  scale_y_continuous(trans = 'log10') +
  labs(title = "Pioneer 10 CPI Electrons and minimum ionizing", x = "Distance from Sun (AU)", y = "Rate (cps)") +
  theme_light()

Integral rate

pcpi <- reshape2::melt(mcpi %>% select(time = ts, "Integral rate" = RID7, "Integral rate Z > 5" = RID7ZG5), id = "time") %>% filter(!is.na(value))
ggplot(data = pcpi) +
# geom_line(aes(x = time, y = value, color = variable)) +
  geom_point(aes(x = time, y = value, color = variable)) +
  scale_y_continuous(trans = 'log10') +
  labs(title = "Pioneer 10 CPI ions at E > 67 MeV/nucleon", x = "Time", y = "Rate (cps)") +
  theme_light()

pcpi <- reshape2::melt(cpitraj %>% select(distance = HRANGP, "Integral rate" = RID7, "Integral rate Z > 5" = RID7ZG5), id = "distance") %>% filter(!is.na(value))
ggplot(data = pcpi) +
# geom_line(aes(x = distance, y = value, color = variable)) +
  geom_point(aes(x = distance, y = value, color = variable)) +
  scale_y_continuous(trans = 'log10') +
  labs(title = "Pioneer 10 CPI ions at E > 67 MeV/nucleon", x = "Distance from Sun (AU)", y = "Rate (cps)") +
  theme_light()

3D Scatter

Average protons rate in Trajectory

p3d <- cpitraj%>% select(x, y, z, Rate = RIDMP)
plot_ly(p3d, x = ~x, y = ~y, z = ~z, color = ~Rate) %>%
  add_markers() %>%
  layout(scene = list(
    xaxis = list(title = 'AU'),
    yaxis = list(title = 'AU'),
    zaxis = list(title = 'AU')))

Average rate of nucleon helium in Trajectory

p3d <- cpitraj%>% select(x, y, z, Rate = RIDMHE)
plot_ly(p3d, x = ~x, y = ~y, z = ~z, color = ~Rate) %>%
  add_markers() %>%
  layout(scene = list(
    xaxis = list(title = 'AU'),
    yaxis = list(title = 'AU'),
    zaxis = list(title = 'AU')))

Average rate of electrons and ionizing in Trajectory

p3d <- cpitraj%>% select(x, y, z, Rate = RIDME)
plot_ly(p3d, x = ~x, y = ~y, z = ~z, color = ~Rate) %>%
  add_markers() %>%
  layout(scene = list(
    xaxis = list(title = 'AU'),
    yaxis = list(title = 'AU'),
    zaxis = list(title = 'AU')))

Average integral rate in Trajectory

p3d <- cpitraj%>% select(x, y, z, Rate = RIDM7)
plot_ly(p3d, x = ~x, y = ~y, z = ~z, color = ~Rate) %>%
  add_markers() %>%
  layout(scene = list(
    xaxis = list(title = 'AU'),
    yaxis = list(title = 'AU'),
    zaxis = list(title = 'AU')))