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,]
target3 <- Tann

# 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 <- cbind(0,0,0,0,0,10000)
for (A in seq(7,12,1)) {
for (b in seq(25,35,1)) {
for (T in seq(60,65,1)) {
for (m in seq(40,50,1)) {
for (a in seq(85,95,1)) {

	sina <- (A/100) * sin(((x-b)/T)*(2*pi))
	lina <- x*m/10000 - a/10
	res1 <- target[,2] - (lina + sina)
	ss <- sum(res1^2)
	s <- rbind(s,cbind(A,b,T,m,a,ss))
}}}}}

s[s[,6]==min(s[,6],na.rm=T),]

#A <- s[s[,6]==min(s[,6]),1]
#b <- s[s[,6]==min(s[,6]),2]
#T <- s[s[,6]==min(s[,6]),3]
#m <- s[s[,6]==min(s[,6]),4]
#a <- s[s[,6]==min(s[,6]),5]
#A 
#b
#T
#m 
#a

quit()
A = 9
b = 33
T = 61
m = 46
a = 90

sina <- (A/100) * sin(((x-b)/T)*(2*pi))
lina <- x*m/10000 - a/10


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

res1 <- target[,2] - (lina + sina)
s2 <- sum(res1^2)
s2 # 0.889
cor(target[,2], lina+sina) # 0.886

# add leading and trailing
CairoPNG(paste("img/tt-d-both-line-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 Line and Sine Fit\n1880-2000", sub=paste("http://rhinohide.wordpress.com       ",dstamp), cex.lab=0.75)
linb <- target2[,1]*m/10000 - a/10
sinb <- (A/100) * sin(((target2[,1]-b)/T)*(2*pi))
lines(target2[,1],linb+sinb,col="red") # residuals for annual global T
points(target2[,1],target2[,2],col="green") # residuals for annual global T
lines(target[,1],lina+sina) # residuals for annual global T
points(target[,1],target[,2]) # residuals for annual global T

res1 <- target2[,2] - (linb + sinb)
s2 <- sum(res1^2)
s2 # 0.889
cor(target2[,2], linb+sinb) # 0.886

# add leading and trailing
CairoPNG(paste("img/tt-d-both-line-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 Line and Sine Fit\n1880-2010", sub=paste("http://rhinohide.wordpress.com       ",dstamp), cex.lab=0.75)
linb <- Tann[,1]*m/10000 - a/10
sinb <- (A/100) * sin(((Tann[,1]-b)/T)*(2*pi))
lines(Tann[,1],linb+sinb,col="red") # residuals for annual global T
points(Tann[,1],Tann[,2],col="green") # residuals for annual global T
lines(target[,1],lina+sina) # residuals for annual global T
points(target[,1],target[,2]) # residuals for annual global T

res1 <- Tann[,2] - (linb + sinb)
s2 <- sum(res1^2)
s2 # 0.889
cor(Tann[,2], linb+sinb) # 0.886

#rm(s)
#s <- cbind(0,0,0,0,0,0,10000)
#for (A1 in seq(12,16,1)) {
#for (b1 in seq(1,5,1)) {
#for (T1 in seq(65,75,1)) {
#for (A2 in seq(20,25,1)) {
#for (b2 in seq(120,130,1)) {
#for (T2 in seq(135,145,1)) {
#
#	res2 <- target[,2] - (sina + sinb)
#	ss <- sum(res2^2)
#	s <- rbind(s,cbind(A1,b1,T1,A2,b2,T2,ss))
#}}}}}}

#s[s[,7]==min(s[,7],na.rm=T),]

#A1 <- s[s[,7]==min(s[,7]),1]
#b1 <- s[s[,7]==min(s[,7]),2]
#T1 <- s[s[,7]==min(s[,7]),3]
#A2 <- s[s[,7]==min(s[,7]),4]
#b2 <- s[s[,7]==min(s[,7]),5]
#T2 <- s[s[,7]==min(s[,7]),6]
A1=15
b1=2
T1=71
A2=23
b2=123
T2=141

sin1 <- (A1/100) * sin(((x-b1)/T1)*(2*pi))
sin2 <- (A2/100) * sin(((x-b2)/T2)*(2*pi))

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

res1 <- target[,2] - (sin1 + sin2)
s2 <- sum(res1^2)
s2 # 0.920
cor(target[,2], sin1+sin2) # 0.885


# add leading and trailing
CairoPNG(paste("img/tt-d-both-sine-",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 Two Sine Fit\n1880-2000", sub=paste("http://rhinohide.wordpress.com       ",dstamp), cex.lab=0.75)
sin1a <- (A1/100) * sin(((target2[,1]-b1)/T1)*(2*pi))
sin2a <- (A2/100) * sin(((target2[,1]-b2)/T2)*(2*pi))
lines(target2[,1],sin1a+sin2a,col="red") # residuals for annual global T
points(target2[,1],target2[,2],col="green") # residuals for annual global T
lines(target[,1],sin1+sin2) # residuals for annual global T
points(target[,1],target[,2]) # residuals for annual global T

res1 <- target2[,2] - (sin1a + sin2a)
s2 <- sum(res1^2)
s2 # 0.920
cor(target2[,2], sin1a+sin2a) # 0.885


# add leading and trailing
CairoPNG(paste("img/tt-d-both-sine-",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 Two Sine Fit\n1880-2010", sub=paste("http://rhinohide.wordpress.com       ",dstamp), cex.lab=0.75)
sin1a <- (A1/100) * sin(((Tann[,1]-b1)/T1)*(2*pi))
sin2a <- (A2/100) * sin(((Tann[,1]-b2)/T2)*(2*pi))
lines(Tann[,1],sin1a+sin2a,col="red") # residuals for annual global T
points(Tann[,1],Tann[,2],col="green") # residuals for annual global T
lines(target[,1],sin1+sin2) # residuals for annual global T
points(target[,1],target[,2]) # residuals for annual global T


res1 <- Tann[,2] - (sin1a + sin2a)
s2 <- sum(res1^2)
s2 # 0.920
cor(Tann[,2], sin1a+sin2a) # 0.885