Dependencies

library(tidyverse)
library(lubridate)
library(data.table) # For fread
library(foreach)
library(ggtext)

Download Data

Voyager 1 and 2 daily merged data

Format information see:

pathv1 <- 'https://spdf.gsfc.nasa.gov/pub/data/voyager/voyager1/merged/vy1_'
pathv2 <- 'https://spdf.gsfc.nasa.gov/pub/data/voyager/voyager2/merged/vy2_'
yseqv1 <- seq(from = 2012, to = 2013)
yseqv2 <- seq(from = 2017, to = 2018)
rdmv1 <- data.frame()      # Empty data frame for result
rdmv2 <- data.frame()      # Empty data frame for result
foreach(y = yseqv1) %do% {
  urlv1 <- paste(pathv1, as.character(y), '.asc', sep = '')
  rdmv1 <- rbind(rdmv1, fread(urlv1)) 
}
foreach(y = yseqv2) %do% {
  urlv2 <- paste(pathv2, as.character(y), '.asc', sep = '')
  rdmv2 <- rbind(rdmv2, fread(urlv2))
}

Convert date time

dmv1 <- rdmv1 %>% mutate(ts = date_decimal(V1) + days(V2) + hours(V3))
dmv2 <- rdmv2 %>% mutate(ts = date_decimal(V1) + days(V2) + hours(V3))

Select and rename

dmv1 <- dmv1 %>% select(
  ts,             # Date time
  V1HRANG  = V4,  # Spacecraft Heliographic distance (AU)
  V1HELLAT = V5,  # Heliographic Inertial latitude of the Spacecraft (Degrees, +/-90.)
  V1HELLON = V6,  # Heliographic Inertial longitude of the Spacecraft (Degrees, 0-360)
  V1FMAB   = V7,  # Field Magnitude Average |B| (1/N SUM |B|, nT)
  V1MAF    = V8,  # Magnitude of Average Field (sqrt(Bx^2+By^2+Bz^2), nT)
  V1RTNBR  = V9,  # BR RTN-Coordinate System (nanoteslas)
  V1RTNBT  = V10, # BT RTN-Coordinate System (nanoteslas)
  V1RTNBN  = V11, # BN RTN-Coordinate System (nanoteslas)
  V1RTNPFS = V12, # Proton flow speed, RTN (km/s)
  V1RTNTHE = V13, # THETA-elevation angle of flow velocity vector (degrees in RTN-cordinate system)
  V1RTNPHI = V14, # PHI- azimuth angle of flow velocity vector (degrees in RTN-cordinate system)
  V1PD     = V15, # Proton density (n/cc)
  V1PT     = V16, # Proton Temperature (K) calculated from thermal speed width T=60.5*Vth*Vth
  # LECP H in 1/(cm^2 sec ster MeV)
  V1LECP1  = V17, # 0.57-1.78 H MeV
  V1LECP2  = V18, # 3.40-17.6 H Mev
  V1LECP3  = V19, # 22.0-31.0 H MeV
  # CRS H in 1/(cm^2 sec ster MeV)
  V1CRS01  = V20, # 3.000-4.600 MeV
  V1CRS02  = V21, # 4.600-6.200 MeV
  V1CRS03  = V22, # 6.200-7.700 MeV
  V1CRS04  = V23, # 7.700-12.800 MeV
  V1CRS05  = V24, # 12.800-17.900 MeV
  V1CRS06  = V25, # 17.900-30.000 MeV
  V1CRS07  = V26, # 30.000-48.000 MeV
  V1CRS08  = V27, # 48.000-56.000 MeV
  V1CRS09  = V28, # 74.471-83.661 MeV
  V1CRS10  = V29, # 132.834-154.911 MeV
  V1CRS11  = V30, # 154.911-174.866 MeV
  V1CRS12  = V31, # 174.866-187.713 MeV
  V1CRS13  = V32, # 187.713-220.475 MeV
  V1CRS14  = V33, # 220.475-270.050 MeV
  V1CRS16  = V34) # 270.050- 346.034 MeV
dmv2 <- dmv2 %>% select(
  ts,             # Date time
  V2HRANG  = V4,  # Spacecraft Heliographic distance (AU)
  V2HELLAT = V5,  # Heliographic Inertial latitude of the Spacecraft (Degrees, +/-90.)
  V2HELLON = V6,  # Heliographic Inertial longitude of the Spacecraft (Degrees, 0-360)
  V2FMAB   = V7,  # Field Magnitude Average |B| (1/N SUM |B|, nT)
  V2MAF    = V8,  # Magnitude of Average Field (sqrt(Bx^2+By^2+Bz^2), nT)
  V2RTNBR  = V9,  # BR RTN-Coordinate System (nanoteslas)
  V2RTNBT  = V10, # BT RTN-Coordinate System (nanoteslas)
  V2RTNBN  = V11, # BN RTN-Coordinate System (nanoteslas)
  V2RTNBFS = V12, # Proton flow speed, RTN (km/s)
  V2RTNTHE = V13, # THETA-elevation angle of flow velocity vector (degrees in RTN-cordinate system)
  V2RTNPHI = V14, # PHI- azimuth angle of flow velocity vector (degrees in RTN-cordinate system)
  V2PPD    = V15, # Plasma Proton density (n/cc)
  V2PPT    = V16, # Plasma Proton Temperature (K) calculated from thermal speed width T=60.5*Vth*Vth
  # LECP H in 1/(cm^2 sec ster MeV)
  V2LECP1  = V17, # 0.52-1.45 MeV
  V2LECP2  = V18, # 3.04-17.3 Mev
  V2LECP3  = V19, # 22.0-30.0 MeV
  # CRS H in 1/(cm^2 sec ster MeV)
  V2CRS01  = V20, # 3.000-4.600 MeV
  V2CRS02  = V21, # 4.600-6.200 MeV
  V2CRS03  = V22, # 6.200-7.700 MeV
  V2CRS04  = V23, # 7.700-12.800 MeV
  V2CRS05  = V24, # 12.800-17.900 MeV
  V2CRS06  = V25, # 17.900-30.000 MeV
  V2CRS07  = V26, # 30.000-48.000 MeV
  V2CRS08  = V27, # 48.000-56.000 MeV
  V2CRS09  = V28, # 75.861-82.562 MeV
  V2CRS10  = V29, # 130.339 - 154.217 MeV
  V2CRS11  = V30, # 154.217 - 171.338 MeV
  V2CRS12  = V31, # 171.338 - 193.643 MeV
  V2CRS13  = V32, # 193.643 - 208.152 MeV
  V2CRS14  = V33, # 208.152 - 245.690 MeV
  V2CRS15  = V34, # 245.690 - 272.300 MeV
  V2CRS16  = V35, # 272.300 - 344.010 MeV
  V2CRS17  = V36, # 344.010 - 478.623 MeV
  V2CRS18  = V37) # 478.623 - 598.667 MeV

Replace fill in value with NA

dmv1 <- na_if(dmv1, 999.99)
dmv1 <- na_if(dmv1, 9999.9)
dmv1 <- na_if(dmv1, 999.999)
dmv1 <- na_if(dmv1, 99.99999)
dmv1 <- na_if(dmv1, 9999999)
dmv1 <- na_if(dmv1, 9.999E5)
dmv2 <- na_if(dmv2, 999.99)
dmv2 <- na_if(dmv2, 9999.9)
dmv2 <- na_if(dmv2, 999.999)
dmv2 <- na_if(dmv2, 99.99999)
dmv2 <- na_if(dmv2, 9999999)
dmv2 <- na_if(dmv2, 9.999E5)

Estimate the speed

From Trajectory data the speed can be estimated

modelv1 <- lm(formula = V1HRANG ~ ts, data = dmv1)
modelv2 <- lm(formula = V2HRANG ~ ts, data = dmv2)
coefficientsv1 <- coef(modelv1)
coefficientsv2 <- coef(modelv2)
speedv1 <- coefficientsv1[2] # AU/μs
speedv2 <- coefficientsv2[2] # AU/μs
speedv1 <- speedv1 * 149597870.7 # Convert to km/s
speedv2 <- speedv2 * 149597870.7 # Convert to km/s
speedv1

##       ts 
## 16.95751

speedv2

##       ts 
## 15.03117

Plots

Heliopause

Two Voyager program spacecraft explored the outer reaches of the heliosphere, passing through the termination shock and the heliosheath. Voyager 1 encountered the heliopause on 25 August 2012, when the spacecraft measured a forty-fold sudden increase in plasma density. Voyager 2 traversed the heliopause on 5 November 2018. Because the heliopause marks the boundary between matter originating from the Sun and matter originating from the rest of the galaxy, spacecraft that depart the heliosphere (such as the two Voyagers) are in interstellar space.

See Heliosphere

Heliosphere upwind Voyager 1 121.6 AU: dashed blue
Heliosphere upwind Voyager 2 119.0 AU: dashed green

hsauv1 <- 121.6
hsauv2 <- 119.0
# According to https://en.wikipedia.org/wiki/Heliosphere
hsdv1 <- as_datetime('2012-08-25')
hsdv2 <- as_datetime('2018-11-05')

Distance from sun

ggplot(data = dmv1, mapping = aes(x = ts, y = V1HRANG)) +
  geom_point(color = 'blue') + geom_smooth(method = 'lm') +
  labs(title = "Voyager 1 Distance from Sun", x = "Date", y = "Heliographic distance (AU)") +
  geom_hline(yintercept = hsauv1, linetype = 'dashed', color = 'blue') +
  geom_vline(xintercept = as.numeric(hsdv1), linetype = 'dashed', color = 'blue') +
  theme_light()

ggplot(data = dmv2, mapping = aes(x = ts, y = V2HRANG)) +
  geom_point(color = 'green') + geom_smooth(method = 'lm') +
  labs(title = "Voyager 2 Distance from Sun", x = "Date", y = "Heliographic distance (AU)") +
  geom_hline(yintercept = hsauv2, linetype = 'dashed', color = 'green') +
  geom_vline(xintercept = as.numeric(hsdv2), linetype = 'dashed', color = 'green') +
  theme_light()

Plot Voyager 2 Proton density

dmv2 %>% filter(!is.na(V2PPD)) %>%
  ggplot(mapping = aes(x = ts, y = V2PPD)) + geom_point(color = 'green') +
  scale_y_continuous(trans = 'log10') +
  labs(title = "Voyager 2 Proton Density", x = "Date", y = "Density (n/cc)") +
  geom_vline(xintercept = as.numeric(hsdv2), linetype = 'dashed', color = 'green') +
  theme_light()

dmv2 %>% filter(!is.na(V2PPD)) %>%
  ggplot(mapping = aes(x = V2HRANG, y = V2PPD)) + geom_point(color = 'green') +
  scale_y_continuous(trans = 'log10') +
  labs(title = "Voyager 2 Proton Density", x = "Distance (AU)", y = "Density (n/cc)") +
  geom_vline(xintercept = hsauv2, linetype = 'dashed', color = 'green') +
  theme_light()

Plot Voyager 1 and 2 LECP H

yat <- 'LECP H cm<sup>-2</sup> sec<sup>-1</sup> ster<sup>-1</sup> MeV<sup>-1</sup>'
yet <- element_textbox_simple(width = NULL, orientation = "left-rotated")
melt(dmv1 %>% select(date = ts,
  "0.57-1.78 MeV" = V1LECP1,
  "3.40-17.6 Mev" = V1LECP2,
  "22.0-31.0 MeV" = V1LECP3), id = "date") %>% filter(!is.na(value)) %>%
  ggplot(mapping = aes(x = date, y = value, color = variable)) + geom_point() +
  scale_y_continuous(trans = 'log10') +
  labs(title = "Voyager 1", x = "Date", y = yat) +
  geom_vline(xintercept = as.numeric(hsdv1), linetype = 'dashed', color = 'blue') +
  theme_light() + theme(axis.title.y = yet)

melt(dmv1 %>% select(distance = V1HRANG,
  "0.57-1.78 MeV" = V1LECP1,
  "3.40-17.6 Mev" = V1LECP2,
  "22.0-31.0 MeV" = V1LECP3), id = "distance") %>% filter(!is.na(value)) %>%
  ggplot(mapping = aes(x = distance, y = value, color = variable)) + geom_point() +
  scale_y_continuous(trans = 'log10') +
  labs(title = "Voyager 1", x = "Distance (AU)", y = yat) +
  geom_vline(xintercept = hsauv1, linetype = 'dashed', color = 'blue') +
  theme_light() + theme(axis.title.y = yet)

melt(dmv2 %>% select(date = ts,
  "0.52-1.45 MeV" = V2LECP1,
  "3.04-17.3 Mev" = V2LECP2,
  "22.0-30.0 MeV" = V2LECP3), id = "date") %>% filter(!is.na(value)) %>%
  ggplot(mapping = aes(x = date, y = value, color = variable)) + geom_point() +
  scale_y_continuous(trans = 'log10') +
  labs(title = "Voyager 2", x = "Date", y = yat) +
  geom_vline(xintercept = as.numeric(hsdv2), linetype = 'dashed', color = 'green') +
  theme_light() + theme(axis.title.y = yet)

melt(dmv2 %>% select(distance = V2HRANG,
  "0.52-1.45 MeV" = V2LECP1,
  "3.04-17.3 Mev" = V2LECP2,
  "22.0-30.0 MeV" = V2LECP3), id = "distance") %>% filter(!is.na(value)) %>%
  ggplot(mapping = aes(x = distance, y = value, color = variable)) + geom_point() +
  scale_y_continuous(trans = 'log10') +
  labs(title = "Voyager 2", x = "Distance (AU)", y = yat) +
  geom_vline(xintercept = hsauv2, linetype = 'dashed', color = 'green') +
  theme_light() + theme(axis.title.y = yet)

Plot Voyager 1 and 2 CRS H

yat <- 'CRS H cm<sup>-2</sup> sec<sup>-1</sup> ster<sup>-1</sup> MeV<sup>-1</sup>'
yet <- element_textbox_simple(width = NULL, orientation = "left-rotated")
melt(dmv1 %>% select(date = ts,
  "3.000-4.600 MeV" = V1CRS01,
  "4.600-6.200 MeV" = V1CRS02,
  "6.200-7.700 MeV" = V1CRS03,
  "7.700-12.800 MeV" = V1CRS04,
  "12.800-17.900 MeV" = V1CRS05,
  "17.900-30.000 MeV" = V1CRS06,
  "30.000-48.000 MeV" = V1CRS07,
  "48.000-56.000 MeV" = V1CRS08,
  "74.471-83.661 MeV" = V1CRS09,
  "132.834-154.911 MeV" = V1CRS10,
  "154.911-174.866 MeV" = V1CRS11,
  "174.866-187.713 MeV" = V1CRS12,
  "187.713-220.475 MeV" = V1CRS13,
  "220.475-270.050 MeV" = V1CRS14,
  "270.050-346.034 MeV" = V1CRS16), id = "date") %>% filter(!is.na(value)) %>%
  ggplot(mapping = aes(x = date, y = value, color = variable)) + geom_point() +
  scale_y_continuous(trans = 'log10') +
  labs(title = "Voyager 1", x = "Date", y = yat) +
  geom_vline(xintercept = as.numeric(hsdv1), linetype = 'dashed', color = 'blue') +
  theme_light() + theme(axis.title.y = yet)

melt(dmv1 %>% select(distance = V1HRANG,
  "3.000-4.600 MeV" = V1CRS01,
  "4.600-6.200 MeV" = V1CRS02,
  "6.200-7.700 MeV" = V1CRS03,
  "7.700-12.800 MeV" = V1CRS04,
  "12.800-17.900 MeV" = V1CRS05,
  "17.900-30.000 MeV" = V1CRS06,
  "30.000-48.000 MeV" = V1CRS07,
  "48.000-56.000 MeV" = V1CRS08,
  "74.471-83.661 MeV" = V1CRS09,
  "132.834-154.911 MeV" = V1CRS10,
  "154.911-174.866 MeV" = V1CRS11,
  "174.866-187.713 MeV" = V1CRS12,
  "187.713-220.475 MeV" = V1CRS13,
  "220.475-270.050 MeV" = V1CRS14,
  "270.050-346.034 MeV" = V1CRS16), id = "distance") %>% filter(!is.na(value)) %>%
  ggplot(mapping = aes(x = distance, y = value, color = variable)) + geom_point() +
  scale_y_continuous(trans = 'log10') +
  labs(title = "Voyager 1", x = "Distance (AU)", y = yat) +
  geom_vline(xintercept = hsauv1, linetype = 'dashed', color = 'blue') +
  theme_light() + theme(axis.title.y = yet)

melt(dmv2 %>% select(date = ts,
  "3.000-4.600 MeV" = V2CRS01,
  "4.600-6.200 MeV" = V2CRS02,
  "6.200-7.700 MeV" = V2CRS03,
  "7.700-12.800 MeV" = V2CRS04,
  "12.800-17.900 MeV" = V2CRS05,
  "17.900-30.000 MeV" = V2CRS06,
  "30.000-48.000 MeV" = V2CRS07,
  "48.000-56.000 MeV" = V2CRS08,
  "75.861-82.562 MeV" = V2CRS09,
  "130.339-154.217 MeV" = V2CRS10,
  "154.217-171.338 MeV" = V2CRS11,
  "171.338-193.643 MeV" = V2CRS12,
  "193.643-208.152 MeV" = V2CRS13,
  "208.152-245.690 MeV" = V2CRS14,
  "245.690-272.300 MeV" = V2CRS15,
  "272.300-344.010 MeV" = V2CRS16,
  "344.010-478.623 MeV" = V2CRS17,
  "478.623-598.667 MeV" = V2CRS18), id = "date") %>% filter(!is.na(value)) %>%
  ggplot(mapping = aes(x = date, y = value, color = variable)) + geom_point() +
  scale_y_continuous(trans = 'log10') +
  labs(title = "Voyager 2", x = "Date", y = yat) +
  geom_vline(xintercept = as.numeric(hsdv2), linetype = 'dashed', color = 'green') +
  theme_light() + theme(axis.title.y = yet)

melt(dmv2 %>% filter(!is.na(V2HRANG)) %>% select(distance = V2HRANG,
  "3.000-4.600 MeV" = V2LECP1,
  "4.600-6.200 MeV" = V2LECP2,
  "6.200-7.700 MeV" = V2LECP3,
  "7.700-12.800 MeV" = V2CRS04,
  "12.800-17.900 MeV" = V2CRS05,
  "17.900-30.000 MeV" = V2CRS06,
  "30.000-48.000 MeV" = V2CRS07,
  "48.000-56.000 MeV" = V2CRS08,
  "75.861-82.562 MeV" = V2CRS09,
  "130.339-154.217 MeV" = V2CRS10,
  "154.217-171.338 MeV" = V2CRS11,
  "171.338-193.643 MeV" = V2CRS12,
  "193.643-208.152 MeV" = V2CRS13,
  "208.152-245.690 MeV" = V2CRS14,
  "245.690-272.300 MeV" = V2CRS15,
  "272.300-344.010 MeV" = V2CRS16,
  "344.010-478.623 MeV" = V2CRS17,
  "478.623-598.667 MeV" = V2CRS18), id = "distance") %>% filter(!is.na(value)) %>%
  ggplot(mapping = aes(x = distance, y = value, color = variable)) + geom_point() +
  scale_y_continuous(trans = 'log10') +
  labs(title = "Voyager 2", x = "Distance (AU)", y = yat) +
  geom_vline(xintercept = hsauv2, linetype = 'dashed', color = 'green') +
  theme_light() + theme(axis.title.y = yet)

Investigating NASA Voyager Merged Data for Heliopause

Geirmundur Orri Sigurdsson

4/26/2021