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 SINE WAVE FOR 20thC 

s <- cbind(0,0,0,10000)
for (A in 5:20) {
for (b in 25:29) {
for (T in 30:60) {
	sina <- (A/100) * sin(((x-b)/T)*(2*pi))
    res1 <- target[,2] - sina
	ss <- sum(res1^2)
	s <- rbind(s,cbind(A,b,T,ss))
}}}
s[s[,4]==min(s[,4]),]
A <- s[s[,4]==min(s[,4]),1]
b <- s[s[,4]==min(s[,4]),2]
T <- s[s[,4]==min(s[,4]),3]
A
b
T
#A = 16
#b = 28
#T = 50

sina <- (A/100) * sin(((x-b)/T)*(2*pi))
res1 <- target[,2] - sina
s2 <- sum(res1^2)
s2 # 2.7
cor(target[,2],sina) # 0.578

# file name to write plot
CairoPNG(paste("img/tt-c-sine-",src[idx-2],"-",y1,"-",t0,".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 Sine Fit", sub=paste("http://rhinohide.wordpress.com       ",dstamp), cex.lab=0.75)
lines(target[,1],sina) # residuals for annual global T
points(target[,1],target[,2]) # residuals for annual global T
dev.off()

# BEST SINE ON RESIDUAL
rm(s)
#res1
s <- cbind(0,0,0,10000)
for (A1 in seq(1,20,1)) {
for (b1 in seq(110,130,1)) {
for (T1 in seq(130,150,1)) {
	sinb <- (A1/100) * sin(((x-b1)/T1)*(2*pi))
    res2 <- res1 - sinb
	ss <- sum(res2^2)
	s <- rbind(s,cbind(A1,b1,T1,ss))
}}}
s[s[,4]==min(s[,4],na.rm=T),]
A1 <- s[s[,4]==min(s[,4],na.rm=T),1]
b1 <- s[s[,4]==min(s[,4],na.rm=T),2]
T1 <- s[s[,4]==min(s[,4],na.rm=T),3]
A1
b1
T1
#A1 = 16
#b1 = 121
#T1 = 141

sinb <- (A1/100) * sin(((x-b1)/T1)*(2*pi))
res2 <- res1 - sinb
s2 <- sum(res2^2)
s2 # 1.24
cor(res1,sinb) # 0.735

# file name to write plot
CairoPNG(paste("img/tt-c-sine-with-sine-residuals-",src[idx-2],"-",y1,"-",t0,".png",sep=""), bg="white")

# plot
# hidden - it just sets the graph axis
plot(tyear,rep(0,length(tyear)), ylim=c(ymin,ymax), col="white", xlab="Year", ylab="Anomaly (deg C)") 
title(main="20th Century Sine Fit of Residuals", sub=paste("http://rhinohide.wordpress.com       ",dstamp), cex.lab=0.75)

# points in the trend
lines(target[,1],sinb) # data points for annual global T
points(target[,1],res1) # data points for annual global T
dev.off()

# COMBINE LINEAR AND SINE

# file name to write plot
CairoPNG(paste("img/tt-c-sine-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 Sine Fit with Sine Fit of Residuals\n1900-2000", sub=paste("http://rhinohide.wordpress.com       ",dstamp), cex.lab=0.75)
lines(target[,1],sina + sinb) # residuals for annual global T
points(target[,1],target[,2]) # residuals for annual global T
dev.off()

res3 <- target[,2]-(sina + sinb)
s3 <- sum(res3^2)
s3 # 1.11
cor(target[,2],sina+sina) # 0.858

# add leading and trailing
CairoPNG(paste("img/tt-c-sine-sine-all-",src[idx-2],"-1880-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 Sine Fit with Sine Fit of Residuals\n1880-2000", sub=paste("http://rhinohide.wordpress.com       ",dstamp), cex.lab=0.75)
sina1 <- (A/100) * sin(((target2[,1]-b)/T)*(2*pi))
sinb1 <- (A1/100) * sin(((target2[,1]-b1)/T1)*(2*pi))
lines(target2[,1],sina1+sinb1,col="red") # residuals for annual global T
points(target2[,1],target2[,2],col="green") # residuals for annual global T
lines(target[,1],sina + sinb) # residuals for annual global T
points(target[,1],target[,2]) # residuals for annual global T

res3 <- target2[,2]-(sina1 + sinb1)
s3 <- sum(res3^2)
s3 # 1.11
cor(target2[,2],sina1+sinb1) # 0.858

# add leading and trailing
CairoPNG(paste("img/tt-c-sine-sine-all-",src[idx-2],"-1880-2010.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 Sine Fit with Sine Fit of Residuals\n1880-2010", sub=paste("http://rhinohide.wordpress.com       ",dstamp), cex.lab=0.75)
sina1 <- (A/100) * sin(((Tann[,1]-b)/T)*(2*pi))
sinb1 <- (A1/100) * sin(((Tann[,1]-b1)/T1)*(2*pi))
lines(Tann[,1],sina1+sinb1,col="red") # residuals for annual global T
points(Tann[,1],Tann[,2],col="green") # residuals for annual global T
lines(target[,1],sina + sinb) # residuals for annual global T
points(target[,1],target[,2]) # residuals for annual global T

res3 <- Tann[,2]-(sina1 + sinb1)
s3 <- sum(res3^2)
s3 # 1.11
cor(Tann[,2],sina1+sinb1) # 0.858

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