###MM03 CALCULATIONS
   #requires 
	#rm(list = ls(all = TRUE)) 

###LOAD EMULATION FUNCTION
	#url<-"C:/Documents and Settings/Stephen McIntyre/My Documents/My Webs/website"
	url<-"http://www.climate2003.com"
	source(file.path(url,"scripts/MM03","rpc.function.txt"))

#LOAD MBH DATA #this collates and may require a minute or so
	source("http://www.climate2003.com/scripts/MM03/read.mann.txt")
		#collects original data and takes about a minute
	#Read 4328 items
	#Read 256 items
	#Read 1162 items
	#Read 1062 items
	#Read 462 items
	#Read 442 items
	#Read 662 items
	#Read 112 items

 #E&E FIGURE 6
	a<-nhmann[1:581,2] #MBH98 reconstruction
	b<-rpc[,1] #rpc 1

#EMULATION OF MBH RPC1 USING MBH DATA

	pc<-pc[1:79,] #temperature pc1
	proxy<-read.table("http://www.climate2003.com/data/MM03/pcproxy.txt") #dim 598 112
	proxy<-ts(proxy,start=1400)
	MBH<-RPC(proxy)[[1]]
	c<-ts(MBH[,1],start=1400,end=1980)
	proxy.MBH<-proxy

#CORRECTION MBH RPC1 USING EMULATION METHOD
	proxy<-read.table("http://www.climate2003.com/data/MM03/proxy4.txt",sep="\t",header=TRUE)
	new<-RPC(proxy)[[1]]
	d<-ts(new[,1],start=1400,end=1980)

#E&E FIGURE 6
	combine<-ts.union(a,b,c,d)
	plot(combine,ann=FALSE,lab=c(7,3,12),las=1,type="l",mai=c(0.5,0.5,0,0))

#######################
## E&E FIGURE 7  
	lambda<-lambda[1:16,2]

   	#RECONSTRUCT w-estimate
	#the MBH periods are 11 EV 1780-1980; 9 from 1760-1779; 8 from 1750-1759; 5 from 1700 to 1749; 4 from 1600 to 1699; 2 from 1450 to 1599; 1 from 1400 to 1449.
	rpc<-new[1:581,]
	count<-apply(!is.na(rpc),1,sum)
	#cbind(1400:1980,count[1:581])

	w<-array  (c(rep(NA,581*1082))  ,  dim=c(581,1082)  )

	period.m<-c(1981,1970,1820,1800,1780,1760,1750,1730,1700,1600,1500,1450,1400)
	period<-period.m-1399
	select<-list(c(1:5,7,9,11,14:16),c(1:5,7,9,11,14:16),c(1:5,7,9,11,14:16),c(1:5,7,9,11,14:16),
		c(1:5,7,9,11,15),c(1:3,5,7,9,11,15),c(1,2,5,11,15),c(1,2,5,11,15),
		c(1,2,11,15),c(1,2),c(1,2),c(1))
	np<-length(select)

	for (iter in 1:(np-1))  {
		w[period[iter+1]:(period[iter]-1),]<- rpc[period[iter+1]:(period[iter]-1),select[[iter]]]  %*% diag( lambda[1:16][select[[iter]] ] )    %*% t(eof) [select[[iter]],]     
		}

	for (iter in np:np)  {
		w[period[iter+1]:(period[iter]-1),]<- as.matrix (rpc[period[iter+1]:(period[iter]-1),select[[iter]]])  %*%   as.matrix(lambda[1])     %*%    t( as.matrix(  t(eof) [select[[iter]],]))
		}
	w<-ts(w,start=1400,end=1980)
	dim(w) #[1]  581 1082

      #EXPANSION AND NH AVERAGE
	nhtemp<-(grid2[,4]>0)
	nil<-unlist(nil);mann.scale<-unlist(mann.scale)
	mann.scale<-mann.scale[!nil]/100
	w<-w[,!nil]
	nhtemp<-nhtemp[!nil]

	w<-scale(w,center=FALSE,scale=1/mann.scale)
	nh.alt<-apply(w[,nhtemp],1,mean)
	#apply(!is.na(w),1,sum)
	#count<-apply(!is.na(w[,nhtemp]),1,sum)
	b<-ts(nh.alt,start=1400,end=1980)

	#SECOND GRAPHIC
	combine.base<-ts.union(a,b)
	plot.ts(combine.base,ann=FALSE)