#FUnctions Kaufman
#make ts - take a data.frame and make an annual time series. Can either interpolate values or not.
  # use interp=FALSE if you merely want to get a time series with NA for missing years  

#extend_window - extend a ts to specified start and end:
#decadal - decadal average of a regularized series 
#rescale - utility rescale of time series to period 980-1800 a la Kaufman 

source("http://www.climateaudit.org/scripts/utilities.txt")
	f=function(x) filter.combine.pad(x,truncated.gauss.weights(21))[,2]
	g=function(X) ts(X[,2],start=X[1,1])


#this looks like it makes 
kaufman=function(A,method) {
	A=A[A$year>0,];names(A)[2]="proxy"
	x=round(tapply(A$proxy,factor(floor(A$year/10),levels=0:199),mean,na.rm=T),2);
	year=seq(5,1995,10)
	h=approxfun(year,x)
	test=ts(h(year),start=5,freq=.1)
	y=rescale(test)
	return(y)
}


rescale=function(x,verbose="default") {
	temp= time(x)>=980& time(x)<=1800
	m0=mean(x[temp],na.rm=T); sd0=sd(x[temp],na.rm=T)
	rescale= round((x-m0)/sd0,2)
	if (verbose=="verbose") {
	  h=function(y) (y-m0)/sd0 
	  inv= function(y) sd0*(y-mean(y[temp])) +m0 
	  rescale=list(x=rescale,mean=m0,sd=sd0,h=h,inv=inv)
	  }
	rescale
	}



make.ts = function(A,h=identity,interp =TRUE) {
  names(A)=c("year","proxy")
  time0=min(A$year):max(A$year)
  test=match(time0, A$year)
  temp=!is.na(test) 
  x=rep(NA,length(time0))
  x[temp]= A$proxy[test[temp]]
  ff=approxfun(time0,x) 
  x=ts( approxfun(time0,x)(time0),start=min(time0),end=max(time0) )
  if(!interp) x[!temp]=NA
  return(x)
}


extend_window=function(x,start0=1,end0=2000) {
	dummy= ts(rep(NA,end0-start0+1),start=start0)
	y=window( ts.union(dummy,x)[,2],start=start0,end=end0)	
	return(y)
}


decadal=function(x,na.rm=T) { 
	decadal=ts(  apply ( array(x,dim=c(10,floor(length(x)/10))),2,mean,na.rm=T) ,start=5,freq=.1)
	decadal
}

#
library(gplots)
opar <- par(mar=c(4,4,2,2),las=1)

#plots- require gplots


plotf3=function(x,col0="pink2",col1=2,ylab0="SD Units" ,xlim1=c(400,2000),ylim1=range(x,na.rm=T),
	start.med=950,end.med=1100,start.mod=1856,end.mod=1980,start0=1980,smoothmethod=TRUE) {   
	year=c(time(x));N=length(x)
	par(mar=c(4,4,2,2),bg="grey80")
	plot(year,x,ylim=ylim1,type="n",xlim=xlim1,axes=FALSE,xlab="",ylab=ylab0)
	usr <- par("usr")
	rect(usr[1],usr[3],start.med  ,usr[4],border=NA,col="lavender")  #   - the part used to fit the model
	rect(start.med   ,usr[3],usr[2],end.med,border=NA,col="lemonchiffon") #   - the part used to make the forecast
	rect(end.med   ,usr[3],usr[2],end.med,border=NA,col="lavender") #   - the part used to make the forecast
	rect(start.mod   ,usr[3],usr[2],end.med,border=NA,col="lemonchiffon") #   - the part used to make the forecast
	rect(end.mod   ,usr[3],usr[2],end.med,border=NA,col="lavender") #   - the part used to make the forecast
	abline(h= (-3:3)*2 , col ="gray" , lty =3)
	abline(v= seq(800,2000,200) , col ="gray" , lty =3)

	temp1=!is.na(x);M=sum(temp1);year=c(time(x))[temp1]
	x=x[temp1]
	lines(year,x,col=col0,type="h")
	if(smoothmethod) lines(year,f(x),col=col1,lwd=2)
	title(paste(title0))
	axis(side=1,col=1,font=2,col.axis=1)
	axis(side=2,col=1,font=2,col.axis=1,las=1)
	#smartlegend(x="left",y="bottom",fill=c(2,4),legend=c("Briffa","Grudd"),inset=0,bg="gray80")
	box();abline(h=0,lty=2)
	#mtext(side=2,ylab0,line=2,font=2)
	}