rm(list=ls())
library(chron)

url1 <- "ftp://ulysses.sr.unh.edu/NeutronMonitor/DailyAverages.1951-.txt"
url2 <- "ftp://ftp.ngdc.noaa.gov/STP/GEOMAGNETIC_DATA/AASTAR/aaindex"

print(paste("Reading GCR data from",url1,"..."))
a <- readLines(url1)
print("Saving GCR data in local file 'neutron2.dat'; skipping the header:")
#print(a[1:43])
a <- a[44:length(a)]
writeLines(a,"neutron2.dat")


gcr <- read.table("neutron2.dat",header=T, 
                  col.names=c("Date","SecsOf1904","Climax","Huancayo.NoGC",
                             "Huancayo.GC","Haleakala.IGY","Haleakala.S/M"))
climax <- filter(as.numeric(gcr$Climax),rep(1,30)/30)
nt.gcr <- length(climax)
climax <- climax[seq(15,nt.gcr,by=30)]
gcr.yy <- as.numeric(substr(as.character(gcr$Date),7,10))
gcr.mm <- as.numeric(substr(as.character(gcr$Date),1,2))
gcr.dd <- as.numeric(substr(as.character(gcr$Date),4,5))
jday.gcr  <- gcr.yy + (gcr.mm-1)/12 + (gcr.dd-0.5)/365.25
jday.gcr <- jday.gcr[seq(15,nt.gcr,by=30)]


x11()
par(cex.sub=0.8)
plot(range(gcr.yy),range(climax),type="n",
     main=paste("Climax Galactic Cosmic Ray (GCR) counts",min(gcr.yy),"-",max(gcr.yy)),lwd=4,
     sub="URL http://ulysses.uchicago.edu/NeutronMonitor/neutron_mon.html (gcr.R)",
     ylab="30-day mean Counts",xlab="Time")
grid()
points(jday.gcr,climax,pch=19,cex=0.7,col="grey50")
lines(jday.gcr,climax,lty=3)
grid()


dev2bitmap(file="gcr.png")