#!/usr/bin/env Rscript
# Read File with Data Points
# Data obtained from National Snow & Ice Data Center
# --------------------------------------------------
rm(list=ls())
filename <- "sea-ice.txt"
tt <- read.table(file=filename,header=TRUE)
tt
x <- tt$Year
y <- tt$Ice
x2 <- x^2

# Polynomial Regression
f1 <- lm(formula = y ~ x + x2)
summary(f1)
anova(f1)
xarr  <- seq(from=min(x),to=max(x),by=0.05)
beta0 <- f1$coefficients[1]
beta1 <- f1$coefficients[2]
beta2 <- f1$coefficients[3]
yarr  <- beta0 + beta1*xarr + beta2 *xarr^2

# Make the plots
outfile <- "sea-ice.pdf"
pdf(file=outfile)
plot(x,y,main="Extent Artic Sea Ice in month of September")
lines(xarr,yarr,lty=3,col="red",lwd=2)
dev.off()