#INPUT HANSEN A,B,C projections
 	hansen88=read.table("http://www.climateaudit.org/data/hansen/hansen_1988.projections.dat",sep="\t",header=TRUE)
		#manual digitization from graph by Willis Eschenbach to 2005; extended by me manually to 2010
	hansen88=ts(hansen88[,2:4],start=fred[1,1])
	   ts.plot(hansen88,col=1:5,type="b")

###TEMPERATURE SPAGHETTI
	url0="http://www.climateaudit.org/scripts/spaghetti"
	source(file.path(url0,"tlt3.glb.txt")) #returns tlt3.glb RSS
	source(file.path(url0,"giss.glb.txt")) #returns giss.glb

	source(file.path("http://www.climateaudit.org/scripts/utilities","ts.annavg.txt")) #returns giss.glb
	spaghetti=ts.union(ts.annavg(tlt3.glb),giss.glb)
	dimnames(spaghetti)[[2]]=c("tlt3.glb","giss.glb")

### re-center satellite to match giss.glb on 1979-1988
	index=(1979:1988)-tsp(spaghetti)[1]+1
	m0=apply(spaghetti[index,],2,mean);m0
 	  #ts.annavg(tlt3.glb)    giss.glb 
          # -0.07725833          0.17000000 
	spaghetti=scale(spaghetti,center=m0-.17,scale=FALSE) # #raise them relative to giss.glb 
	  index=(1979:1988)-tsp(spaghetti)[1]+1
	  apply(spaghetti[index,],2,mean)
		#ts.annavg(tlt3.glb)            giss.glb 
                # 0.17                0.17 

##re-center projections to match GISS.GLB on first 10 years
	id=dimnames(hansen88)[[2]]
	hansen88=ts.union(hansen88,ts(giss.glb[(1958:2007)-tsp(giss.glb)[1]+1],start=1958) )
	dimnames(hansen88)[[2]]=c(id,"giss.glb")	
	index=(1958:1967)-tsp(hansen88)[1]+1
	m1=apply(hansen88[index,],2,mean);m1
   	   #Scenario.A Scenario.B Scenario.C   giss.glb 
           #-0.051     -0.062     -0.062     -0.002 
	hansen88=scale(hansen88,center=m1+.002,scale=FALSE)
	  index=(1958:1967)-tsp(hansen88)[1]+1
	  apply(hansen88[index,],2,mean); #all centered to match

	index=(1951:1980)-tsp(hansen88)[1]+1
	apply(hansen88[index,],2,mean,na.rm=T)
   	   #Scenario.A Scenario.B Scenario.C   giss.glb 
           #-0.051     -0.062     -0.062     -0.002 
	hansen88=scale(hansen88,center=m1+.002,scale=FALSE)
	  index=(1958:1967)-tsp(hansen88)[1]+1
	  apply(hansen88[index,],2,mean); #all centered to match


##COMPARE TO GAVIN PLOT 
	par(mar=c(3,3,2,1))
	plot(c(time(hansen88)),hansen88[,"Scenario.A"],col=2,ylim=c(-.25,1.5),xlim=c(1958,2010),xlab="",ylab="",tcl=.25,axes=FALSE)
	box();axis(side=1,tcl=.25);axis(side=2,las=1,font=2)
	lines(c(time(hansen88)),hansen88[,"Scenario.A"],col=2,lwd=2)
	lines(c(time(hansen88)),hansen88[,"Scenario.B"],col="green4",lwd=2)
	points(c(time(hansen88)),hansen88[,"Scenario.B"],col="green4",pch=1)
	lines(c(time(hansen88)),hansen88[,"Scenario.C"],col=4,lwd=2)
	points(c(time(hansen88)),hansen88[,"Scenario.C"],col=4,pch=1)
	lines(c(time(giss.glb)),giss.glb,col="grey60",lwd=2)
	points(c(time(giss.glb)),giss.glb,col="grey60",pch=19)
	lines(c(time(spaghetti)),spaghetti[,"tlt3.glb"],col=1,lwd=3)
	points(c(time(spaghetti)),spaghetti[,"tlt3.glb"],col=1,pch=19)
	abline(v=1987,lty=3,col="grey80")
	abline(h=seq(0,1.5,.5),col="grey80",lty=2)
	legend(1958,1.6,fill=c(2,3,4,"grey60",1),legend=c("Hansen A","Hansen B","Hansen C","GISS Surf","RSS Sat"))
	title(main="Hansen et al 1988 Projections")


##COMPARE MEANS
	test=ts.union(hansen88[,1:3],ts(giss.glb[(1951:2007)-tsp(giss.glb)[1]+1],start=1951) )
	dimnames(test)[[2]]=c(id,"giss.glb")	
	index=(1951:1980)-tsp(test)[1]+1
	apply(test[index,],2,mean,na.rm=T);
   	   #Scenario.A Scenario.B Scenario.C   giss.glb 
	   #0.10200000 0.06978261 0.06608696 0.00100000 
	index=(1958:1980)-tsp(test)[1]+1
	apply(test[index,],2,mean,na.rm=T);
	  # Scenario.A  Scenario.B  Scenario.C    giss.glb 
	  #0.102000000 0.069782609 0.066086957 0.009565217 
   
	c(mean(giss.glb[(1958:1980)-1879]),mean(giss.glb[(1951:1980)-1879]),mean(giss.glb[(1958:1967)-1879]))
		#   0.009565217  0.001000000 -0.002000000