#define exponential smoothing function
expsmooth<-function(v,d){ # u is the smoothed v with decay exp(-1/d)
n=1:length(ss)
e=(n-1)/d
e=exp(-rev(e))
e=e/sum(e)
u=convolve(v,e,type="open")[n]
u
}

#read GISS forcing data
s <- matrix(scan("/Users/apsmith/energy/climate/two_box/gissforc.txt", 0, skip=4), ncol=11, byrow=TRUE)
vv=rowSums(s[,2:11])

#read GISS temp data
v <- matrix(scan("/Users/apsmith/energy/climate/two_box/gistemp.txt", 0, skip=9, nlines=124), ncol=20, byrow=TRUE)
ss=v[,14]/100.

w1=expsmooth(vv,30.0)
w2=expsmooth(vv,1.0)

rplus = 967.0
rminus = -0.000131

wsplus = 0.000842
wsminus = 0.417
woplus = rplus*wsplus
wominus = rminus*wsminus

ts = wsplus*w1 + wsminus*w2
to = woplus*w1 + wominus*w2

#plotting
t=1880:2003
k=1:124
#jpeg(width=500, height=500)
plot(t[k],ss[k],type="l",ylim=c(-0.6,1.0),xlab="Year",ylab="Temperature anomaly", bty="n", main=paste("Tamino (-) fit with gamma=3.17e-11"))
lines(t[k],ts[k],col="red")
lines(t[k],to[k],col="green")
legend(x=1900,y=0.5,legend=c("Tm", "Ts", "To"), col=c("black","red","green"),lwd=1)

#dev.off()