GOES-13 proton flux for particles >= 10 Mev

goes.R

# Plots GOES-13 Proton Flux for particles >= 10 Mev.
#
# Data source: ftp://ftp.swpc.noaa.gov/pub/lists/particle/
#
# Check out http://stuffin.space to see where GOES-13 flies.

library(dplyr)
library(tidyr)
library(lubridate)
library(ggplot2)
library(scales)

options(scipen=10000)
theme_set(theme_bw())

ingest_particle <- function(filename) {
  cols <- c("Time", "M", "D", "HHMM", "JDay", "JSec", "P1", "P5", "P10", "P30",
            "P50", "P100", "E0.8", "E2.0", "E4.0")
  read.table(filename, col.names=cols, skip=2) %>%
    filter(P1 >= 0) %>%
    mutate(Time=ymd(Time * 10000 + M * 100 + D) + hours(round(HHMM / 100)) +
           minutes(HHMM %% 100)) %>%
    select(-(2:6)) 
}

goes13_files <- list.files("particle", ".*_Gp_part_5m\\.txt", full.names=TRUE)
goes13 <- do.call(rbind, lapply(goes13_files, ingest_particle))

p <- goes13 %>%
  group_by(y=year(Time), m=month(Time), d=day(Time)) %>%
  summarize(time=min(Time), min=min(P10), median=median(P10), max=max(P10)) %>%
  ggplot(aes(x=time, y=median, ymin=min, ymax=max)) +
    geom_line(color=I("red")) +
    geom_ribbon(alpha=.25) +
    scale_y_log10(breaks=10^(-1:5)) +
    annotation_logticks(sides="l") +
    labs(x="Time", y=expression("Particles at" >= "10 Mev"))

ggsave("goes-13.pdf", p)
ggsave("goes-13.png", p)