# gaussian smoothing with back padding -0.4, front padding obs then 0.2
gsmoothb<-function(v,d){ # u is the smoothed v with decay exp(-n^2/d^2)
smlength=100
n=-smlength:smlength
e=n/d
e=exp(-e*e)
e=e/sum(e)

vcount = length(v)
avglength = 15
if(vcount < avglength) avglength = vcount
avg1 = -0.4
avg2 = 0.2

obs30yr = c(0.15, 0.05, 0.02, -0.23, -0.35, -0.06, -0.01, -0.18, -0.2, -0.02, -0.14, -0.3, 0.05, -0.1, 0.19, -0.23, 0.23, -0.13, -0.08, 0.07, 0.32, 0, 0.18, 0.03, 0.01, 0.07, 0.11, 0.5, 0.37, 0.5)

v2 = c(rep(avg1, times=smlength), v, obs30yr, rep(avg2, times=smlength-30))

u=filter(v2,e,method="convolution",sides=2)
u[(smlength+1):(vcount+smlength)]
}

# define gaussian smoothing with front and back padding by mean
gsmooth<-function(v,d){ # u is the smoothed v with decay exp(-n^2/d^2)
smlength=100
n=-smlength:smlength
e=n/d
e=exp(-e*e)
e=e/sum(e)

vcount = length(v)
avglength = 15
if(vcount < avglength) avglength = vcount
avg1 = sum(v[1:avglength])/avglength
avg2 = sum(v[(vcount-avglength+1):vcount])/avglength

v2 = c(rep(avg1, times=smlength), v, rep(avg2, times=smlength))

u=filter(v2,e,method="convolution",sides=2)
u[(smlength+1):(vcount+smlength)]
}


#read reconstruction data
s <- matrix(scan("/Users/apsmith/energy/climate/temp_recon/briffa2001_data.txt", 0, skip=9), ncol=8, byrow=TRUE)

#   1: Jones et al. (1998) Holocene
#   2: Mann et al. (1999) Geophys Res Lett
#   3: Briffa et al. (2001) J Geophys Res
#   4: Briffa (2000) Quat Sci Rev
#   5: Overpeck et al. (1997) Science
#   6: Crowley & Lowery (2000) Ambio
#   7: Observed temperatures from Jones et al. (1999) Rev Geophys
years = s[,1]
jones = s[,2]
xjones = years[jones > -10]
yjones = gsmooth(jones[jones > -10], 10.)

mann = s[,3]
xmann = years[mann > -10]
ymann = gsmooth(mann[mann > -10], 10.)

briffa1 = s[,4]
xbriffa1 = years[briffa1 > -10]
ybriffa1 = gsmoothb(briffa1[briffa1 > -10], 10.)

briffa2 = s[,5]
xbriffa2 = years[briffa2 > -10]
ybriffa2 = gsmooth(briffa2[briffa2 > -10], 10.)

overpeck = s[,6]
xoverpeck = years[overpeck > -10]
yoverpeck = gsmooth(overpeck[overpeck > -10], 10.)

crowley = s[,7]
xcrowley = years[crowley > -10]
ycrowley = gsmooth(crowley[crowley > -10], 10.)

obs = s[,8]
xobs = years[obs > -10]
yobs = gsmooth(obs[obs > -10], 10.)

#plotting
t=1000:2000
k=1:1000
#jpeg(width=500, height=500)
plot(xjones,yjones,type="l",xlab="Year",ylab="temp anomaly",asp=300.0,xlim=c(1000,2000),ylim=c(-0.8,0.4),xaxt="n",yaxt="n",col="red",main="Instrumental padding")
lines(xmann,ymann,col="darkblue")
lines(xbriffa1,ybriffa1,col="blue")
lines(xbriffa2,ybriffa2,col="green")
lines(xoverpeck,yoverpeck,col="brown")
lines(xcrowley,ycrowley,col="orange")
lines(xobs,yobs,col="black")
lines(c(1000,2000),c(0.0, 0.0))

q=0:6*0.2-0.8
axis(2,at=q)
axis(4,at=q)
q=0:10*100+1000
axis(1,q)

legend(x=1000,y=-0.35,legend=c("Jones", "Mann", "Briffa 2001", "Briffa 2000", "Overpeck", "Crowley", "Observations"), col=c("red","darkblue","blue","green","brown","orange","black"),lwd=1,cex=0.5)
#dev.off()