#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.

#read SOI index data
enso <- matrix(scan("/Users/apsmith/energy/climate/two_box/SOI.signal.ascii", 0, skip=0, nlines=124), ncol=13, byrow=TRUE)
enso_av = rowSums(enso[,2:13])/12.

# loop over shifts:
shifts = 1:50

vshifted <- list()
length(vshifted) <- 20*124
dim(vshifted) <- c(20, 124)

for (shift in shifts) {
	vshifted[[shift]][1:shift] <- 0
	vshifted[[shift]][shift:124] <- vv[1:(125-shift)]
}

# fit regression
h<-lm(ss ~ enso_av + vshifted[[1]] + vshifted[[2]] + vshifted[[3]] + vshifted[[4]])
a4 = h$coefficients[3:6]

h<-lm(ss ~ enso_av + vshifted[[1]] + vshifted[[2]] + vshifted[[3]] + vshifted[[4]] + vshifted[[5]] + vshifted[[6]] + vshifted[[7]] + vshifted[[8]])
a8 = h$coefficients[3:10]

h<-lm(ss ~ enso_av + vshifted[[1]] + vshifted[[2]] + vshifted[[3]] + vshifted[[4]] + vshifted[[5]] + vshifted[[6]] + vshifted[[7]] + vshifted[[8]] + vshifted[[9]] + vshifted[[10]] + vshifted[[11]] + vshifted[[12]])
a12 = h$coefficients[3:14]

h<-lm(ss ~ enso_av + vshifted[[1]] + vshifted[[2]] + vshifted[[3]] + vshifted[[4]] + vshifted[[5]] + vshifted[[6]] + vshifted[[7]] + vshifted[[8]] + vshifted[[9]] + vshifted[[10]] + vshifted[[11]] + vshifted[[12]] + vshifted[[13]] + vshifted[[14]] + vshifted[[15]] + vshifted[[16]])
a16 = h$coefficients[3:18]

h<-lm(ss ~ enso_av + vshifted[[1]] + vshifted[[2]] + vshifted[[3]] + vshifted[[4]] + vshifted[[5]] + vshifted[[6]] + vshifted[[7]] + vshifted[[8]] + vshifted[[9]] + vshifted[[10]] + vshifted[[11]] + vshifted[[12]] + vshifted[[13]] + vshifted[[14]] + vshifted[[15]] + vshifted[[16]] + vshifted[[17]] + vshifted[[18]] + vshifted[[19]] + vshifted[[20]])
a20 = h$coefficients[3:22]

h<-lm(ss ~ enso_av + vshifted[[1]] + vshifted[[2]] + vshifted[[3]] + vshifted[[4]] + vshifted[[5]] + vshifted[[6]] + vshifted[[7]] + vshifted[[8]] + vshifted[[9]] + vshifted[[10]] + vshifted[[11]] + vshifted[[12]] + vshifted[[13]] + vshifted[[14]] + vshifted[[15]] + vshifted[[16]] + vshifted[[17]] + vshifted[[18]] + vshifted[[19]] + vshifted[[20]] + vshifted[[21]] + vshifted[[22]] + vshifted[[23]] + vshifted[[24]] + vshifted[[25]] + vshifted[[26]] + vshifted[[27]] + vshifted[[28]] + vshifted[[29]] + vshifted[[30]])
a30 = h$coefficients[3:32]

h<-lm(ss ~ enso_av + vshifted[[1]] + vshifted[[2]] + vshifted[[3]] + vshifted[[4]] + vshifted[[5]] + vshifted[[6]] + vshifted[[7]] + vshifted[[8]] + vshifted[[9]] + vshifted[[10]] + vshifted[[11]] + vshifted[[12]] + vshifted[[13]] + vshifted[[14]] + vshifted[[15]] + vshifted[[16]] + vshifted[[17]] + vshifted[[18]] + vshifted[[19]] + vshifted[[20]] + vshifted[[21]] + vshifted[[22]] + vshifted[[23]] + vshifted[[24]] + vshifted[[25]] + vshifted[[26]] + vshifted[[27]] + vshifted[[28]] + vshifted[[29]] + vshifted[[30]] + vshifted[[31]] + vshifted[[32]] + vshifted[[33]] + vshifted[[34]] + vshifted[[35]] + vshifted[[36]] + vshifted[[37]] + vshifted[[38]] + vshifted[[39]] + vshifted[[40]])
a40 = h$coefficients[3:42]

h<-lm(ss ~ enso_av + vshifted[[1]] + vshifted[[2]] + vshifted[[3]] + vshifted[[4]] + vshifted[[5]] + vshifted[[6]] + vshifted[[7]] + vshifted[[8]] + vshifted[[9]] + vshifted[[10]] + vshifted[[11]] + vshifted[[12]] + vshifted[[13]] + vshifted[[14]] + vshifted[[15]] + vshifted[[16]] + vshifted[[17]] + vshifted[[18]] + vshifted[[19]] + vshifted[[20]] + vshifted[[21]] + vshifted[[22]] + vshifted[[23]] + vshifted[[24]] + vshifted[[25]] + vshifted[[26]] + vshifted[[27]] + vshifted[[28]] + vshifted[[29]] + vshifted[[30]] + vshifted[[31]] + vshifted[[32]] + vshifted[[33]] + vshifted[[34]] + vshifted[[35]] + vshifted[[36]] + vshifted[[37]] + vshifted[[38]] + vshifted[[39]] + vshifted[[40]] + vshifted[[41]] + vshifted[[42]] + vshifted[[43]] + vshifted[[44]] + vshifted[[45]] + vshifted[[46]] + vshifted[[47]] + vshifted[[48]] + vshifted[[49]] + vshifted[[50]])
a50 = h$coefficients[3:52]

h<-lm(ss ~ enso_av + vshifted[[1]] + vshifted[[2]] + vshifted[[3]] + vshifted[[4]] + vshifted[[5]] + vshifted[[6]] + vshifted[[7]] + vshifted[[8]] + vshifted[[9]] + vshifted[[10]] + vshifted[[11]] + vshifted[[12]] + vshifted[[13]] + vshifted[[14]] + vshifted[[15]] + vshifted[[16]] + vshifted[[17]] + vshifted[[18]] + vshifted[[19]] + vshifted[[20]] + vshifted[[21]] + vshifted[[22]] + vshifted[[23]] + vshifted[[24]] + vshifted[[25]] + vshifted[[26]] + vshifted[[27]] + vshifted[[28]] + vshifted[[29]] + vshifted[[30]] + vshifted[[31]] + vshifted[[32]] + vshifted[[33]] + vshifted[[34]] + vshifted[[35]] + vshifted[[36]] + vshifted[[37]] + vshifted[[38]] + vshifted[[39]] + vshifted[[40]] + vshifted[[41]] + vshifted[[42]] + vshifted[[43]] + vshifted[[44]] + vshifted[[45]] + vshifted[[46]] + vshifted[[47]] + vshifted[[48]] + vshifted[[49]] + vshifted[[50]])

#plotting
#jpeg(width=500, height=500)
n=1:50
plot(n,a50[n],type="l",xlim=c(1,200),ylim=c(-0.05,0.1),xlab="t (years)",ylab="r(t) (K per W/m^2 per yr)",log="x",bty="n",main=paste("Direct response fit, up to 50 terms"))

lines(n,0+0.0*n,col="black")

n=1:4
lines(n,a4[n],col="red")

n=1:8
lines(n,a8[n],col="brown")

n=1:12
lines(n,a12[n],col="green")

n=1:16
lines(n,a16[n],col="blue")

n=1:20
lines(n,a20[n],col="yellow")

n=1:30
lines(n,a30[n],col="orange")

n=1:40
lines(n,a40[n],col="violet")

n=1:50
lines(n,a50[n],col="black")

legend(x=60,y=0.09,legend=c("4 terms", "8 terms","12 terms", "16 terms", "20 terms", "30 terms", "40 terms", "50 terms"), col=c("red","brown","green", "blue", "yellow", "orange", "violet", "black"),lwd=1)

#dev.off()