library(Cairo)

#t <- read.table("http://processtrends.com/Files/RClimate_consol_temp_anom_latest.csv", sep=",", header=T)
t <- read.table("RClimate_consol_temp_anom_latest.csv", sep=",", header=T)

# note: this does not yet handle NAs in the trend well
# be careful with the bounderies - 1979 for sats, pre-1900 for ground
# need boundary tests

y1 <- 2000 # end of trend
t0 <- 100 # years of trend
y0 <- y1 - t0 # start of trend
t1 <- t0 # number of years to display 
y2 <- y0 + t1 # max of x-axis
ymin <- -0.5 # min of y-axis
ymax <- 0.8 # max of y-axis
dstamp=Sys.Date()

# column idx
gis <- 3
cru <- 4
noa <- 5
rss <- 6
uah <- 7
idx <- cru

src = c("GISTEMP", "HADCRU", "NOAA", "RSS", "UAH")

# annualize monthly data
t$year <- as.integer(t[,1])

# get the record of interest, monthly to annual
Tann <- aggregate(t[,idx],list(t$year), FUN=mean, na.rm=T)
target <- Tann[Tann[,1]> y0 & Tann[,1]<= y1,]
target2 <- Tann[Tann[,1]<= y1,]

# a series for x-axis
tyear <- 1880:2010

# some set setup
x <- target[,1]
y <- target[,2]
new <- data.frame(x = seq(y0+1,y1))

# BEST EXPONENTIAL TREND
s <- NULL
for (a in seq(130,150,1)) {
for (bb in seq(5,15,1)) {
for (k in seq(30,60,1)) {
    y1= -a/100 + (bb/10) * exp((k/10000)*(target[,1]-1880))
    s <- rbind(s, c(a,bb,k,sum((target[,2]-y1)^2)))
}}}

s[s[,4]==min(s[,4],na.rm=T),]
a <- s[s[,4]==min(s[,4],na.rm=T),1]
bb <- s[s[,4]==min(s[,4],na.rm=T),2]
k <- s[s[,4]==min(s[,4],na.rm=T),3]
a
bb
k
#a = 149
#bb = 11
#k = 40

    y1= -a/100 + (bb/10) * exp((k/10000)*(target[,1]-1880))
    y2= -a/100 + (bb/10) * exp((k/10000)*(target2[,1]-1880))
    y3= -a/100 + (bb/10) * exp((k/10000)*(Tann[,1]-1880))
    y4= -a/100 + (bb/10) * exp((k/10000)*(1750:2100-1880))

# residuals
res <- target[,2] - y1
sum(res^2)
cor(target[,2],y1)

# 20th
CairoPNG(paste("img/tt-e-exp-",src[idx-2],"-1900-2000.png",sep=""), bg="white")
plot(target[,1],y1,ylim=c(-1,1),xlim=c(1880,2010),type="l")
title(main="20th Century Exponential Trend", sub=paste("http://rhinohide.wordpress.com       ",dstamp), cex.lab=0.75)
points(target[,1],target[,2]) # 
dev.off()

# long view
CairoPNG(paste("img/tt-e-exp-",src[idx-2],"-1750-2100.png",sep=""), bg="white")
plot(target[,1],y1,ylim=c(-1,1),xlim=c(1750,2100),type="l", main=)
lines(1750:2100,y4,type="l")
points(Tann[,1],Tann[,2],col="green") # 
points(target[,1],target[,2]) # 
dev.off()

# plot resid
CairoPNG(paste("img/tt-e-exp-residuals-",src[idx-2],"-1900-2000.png",sep=""), bg="white")
plot(target[,1],res,ylim=c(ymin,ymax),xlim=c(1880,2010))
title(main="20th Century Exponential Trend Residuals", sub=paste("http://rhinohide.wordpress.com       ",dstamp), cex.lab=0.75)
dev.off()

# BEST SINE WAVE FOR RESIDUALS

s <- NULL
for (A in 0:20) {
for (b in 20:30) {
for (T in 50:60) {
	sina <- (A/100) * sin(((x-b)/T)*(2*pi))
	s <- rbind(s,cbind(A,b,T,sum((res-sina)^2)))
}}}
#s[s[,4]==min(s[,4],na.rm=T),]
A <- s[s[,4]==min(s[,4],na.rm=T),1]
b <- s[s[,4]==min(s[,4],na.rm=T),2]
T <- s[s[,4]==min(s[,4],na.rm=T),3]
A
b
T
#A = 8
#b = 24
#T = 56

# file name to write plot
CairoPNG(paste("img/tt-e-exp-with-sine-residuals-",src[idx-2],"-1900-2000.png",sep=""),bg="white")
plot(tyear,rep(0,length(tyear)), ylim=c(ymin,ymax), col="white", xlab="Year", ylab="Anomaly (deg C)") 
title(main="20th Century Sinusoidal Trend of Exponential Residuals", sub=paste("http://rhinohide.wordpress.com       ",dstamp), cex.lab=0.75)
lines(target[,1],sina) # 
points(target[,1],res) # 
dev.off()

sina <- (A/100) * sin(((x-b)/T)*(2*pi))
res2 <- res - sina
s2 <- sum(res2^2)
s2 # 0.857 

# COMBINE LINEAR AND SINE
CairoPNG(paste("img/tt-e-exp-sine-",src[idx-2],"-1900-2000.png",sep=""),bg="white")
plot(tyear,rep(0,length(tyear)), ylim=c(ymin,ymax), col="white", xlab="Year", ylab="Anomaly (deg C)") 
title(main="20th Century Exponential Trend with Sine Fit to Residuals\n1900-2000", sub=paste("http://rhinohide.wordpress.com       ",dstamp), cex.lab=0.75)
lines(target[,1],y1 + sina) # 
points(target[,1],target[,2]) # 
dev.off()

res3 <- target[,2]-(y1 + sina)
s3 <- sum(res3^2)
s3
cor(target[,2],y1+sina) # 0.889

# add leading and trailing
CairoPNG(paste("img/tt-e-exp-sine-all-",src[idx-2],"-1880-2000.png",sep=""),col="white")
plot(tyear,rep(0,length(tyear)), ylim=c(ymin,ymax), col="white", xlab="Year", ylab="Anomaly (deg C)") 
title(main="20th Century Exponential Trend with Sine Fit of Residuals\1880-2000", sub=paste("http://rhinohide.wordpress.com       ",dstamp), cex.lab=0.75)
sinb <- (A/100) * sin(((target2[,1]-b)/T)*(2*pi))
lines(target2[,1],y2+sinb,col="red") # 
points(target2[,1],target2[,2],col="green") # 
lines(target[,1],y1 + sina) # 
points(target[,1],target[,2]) # 

res3 <- target2[,2]-(y2 + sinb)
s3 <- sum(res3^2)
s3
cor(target2[,2],y2+sinb) # 0.889

# add leading and trailing
CairoPNG(paste("img/tt-e-exp-sine-all-",src[idx-2],"-1880-2010.png",sep=""),col="white")
plot(tyear,rep(0,length(tyear)), ylim=c(ymin,ymax), col="white", xlab="Year", ylab="Anomaly (deg C)") 
title(main="20th Century Exponential Trend with Sine Fit of Residuals\n1880-2010", sub=paste("http://rhinohide.wordpress.com       ",dstamp), cex.lab=0.75)
sinb <- (A/100) * sin(((Tann[,1]-b)/T)*(2*pi))
lines(Tann[,1],y3+sinb,col="red") # 
points(Tann[,1],Tann[,2],col="green") # 
lines(target[,1],y1 + sina) # 
points(target[,1],target[,2]) # 

res3 <- Tann[,2]-(y3 + sinb)
s3 <- sum(res3^2)
s3
cor(Tann[,2],y3+sinb) # 0.889


# add leading and trailing
CairoPNG(paste("img/tt-e-exp-sine-all-",src[idx-2],"-1750-2100.png",sep=""),col="white")
plot(tyear,rep(0,length(tyear)), ylim=c(ymin,ymax), col="white", xlab="Year", ylab="Anomaly (deg C)") 
title(main="20th Century Exponential Trend with Sine Fit of Residuals\n1750-2100", sub=paste("http://rhinohide.wordpress.com       ",dstamp), cex.lab=0.75)
sinb <- (A/100) * sin(((1750:2100-b)/T)*(2*pi))
lines(Tann[,1],y4+sinb,col="red") # 
points(Tann[,1],Tann[,2],col="green") # 

# clean up
#rm(list=ls(all=TRUE))