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))

# SOLVE BEST LINE AND SINE WAVE FOR 20thC TOGETHER

#s <- NULL
#for (A in seq(7,11,1)) {
#for (b in seq(32,36,1)) {
#for (T in seq(59,63,1)) {
#for (k in seq(36,40,1)) {
#for (bb in seq(9,13,1)) {
#for (a in seq(143,147,1)) {
#
	#sina <- (A/100) * sin(((target[,1]-b)/T)*(2*pi))
    	#y1= -a/100 + (bb/10) * exp((k/10000)*(target[,1]-1880))
	#res1 <- target[,2] - (y1 + sina)
	#ss <- sum(res1^2)
	#s <- rbind(s,cbind(A,b,T,k,bb,a,ss))
#}}}}}}
#
#s[s[,7]==min(s[,7],na.rm=T),]
#
#A <- s[s[,7]==min(s[,7]),1]
#b <- s[s[,7]==min(s[,7]),2]
#T <- s[s[,7]==min(s[,7]),3]
#k <- s[s[,7]==min(s[,7]),4]
#bb <- s[s[,7]==min(s[,7]),5]
#a <- s[s[,7]==min(s[,7]),6]
#A 
#b
#T
#k
#bb
#a
#quit()
A = 9
b = 34
T = 61
k = 37
bb = 11
a = 145

# file name to write plot
sina <- (A/100) * sin(((target[,1]-b)/T)*(2*pi))
y1= -a/100 + (bb/10) * exp((k/10000)*(target[,1]-1880))

CairoPNG(paste("img/tt-e-both-exp-and-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 Simultaneous Exponential and Sine Fit\n1900-2000", sub=paste("http://rhinohide.wordpress.com       ",dstamp), cex.lab=0.75)
lines(target[,1],y1+sina) # residuals for annual global T
points(target[,1],target[,2]) # residuals for annual global T
dev.off()

res1 <- target[,2] - (y1 + sina)
s2 <- sum(res1^2)
s2 # 0.819
cor(target[,2], y1+sina) # 0.891

# add leading and trailing
sinb <- (A/100) * sin(((target2[,1]-b)/T)*(2*pi))
y2= -a/100 + (bb/10) * exp((k/10000)*(target2[,1]-1880))

CairoPNG(paste("img/tt-e-both-exp-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 Simultaneous Exponential and Sine Fit\n1880-2000", sub=paste("http://rhinohide.wordpress.com       ",dstamp), cex.lab=0.75)
lines(target2[,1],y2+sinb,col="red") # residuals for annual global T
points(target2[,1],target2[,2],col="green") # residuals for annual global T
lines(target[,1],y1+sina) # residuals for annual global T
points(target[,1],target[,2]) # residuals for annual global T

res1 <- target2[,2] - (y2 + sinb)
s2 <- sum(res1^2)
s2 # 0.819
cor(target2[,2], y2+sinb) # 0.891

# add leading and trailing
sinc <- (A/100) * sin(((Tann[,1]-b)/T)*(2*pi))
y3= -a/100 + (bb/10) * exp((k/10000)*(Tann[,1]-1880))

CairoPNG(paste("img/tt-e-both-exp-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 Simultaneous Exponential and Sine Fit\n1880-2010", sub=paste("http://rhinohide.wordpress.com       ",dstamp), cex.lab=0.75)
lines(Tann[,1],y3+sinc,col="red") # residuals for annual global T
points(Tann[,1],Tann[,2],col="green") # residuals for annual global T
lines(target[,1],y1+sina) # residuals for annual global T
points(target[,1],target[,2]) # residuals for annual global T

res1 <- Tann[,2] - (y3 + sinc)
s2 <- sum(res1^2)
s2 # 0.819
cor(Tann[,2], y3+sinc) # 0.891

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