#LOAD MBH NH

load("c:/climate/data/mann/nhmann.tab")

a<-nhmann[,2]

#LOAD MBH PC1

load("c:/climate/data/mann/rpc.tab")

b<-rpc[,1]

#CALCULATE RECONSTRUCTION OF MBH USING MBH DATA

weights1<-scan("c:/climate/data/mann/weights.txt",skip=1)

load("c:/climate/data/mann/pc.tab")

pc<-pc[1:79,]

load("c:/climate/data/mann/proxy.tab")

MBH<-RPC(proxy)[[1]]

c<-ts(MBH[,1],start=1400,end=1980)

load("c:/climate/data/mann/proxy4.tab") #h

new<-RPC(proxy)[[1]]

d<-ts(new[,1],start=1400,end=1980)

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))

loc.out<-"c:/climate/images/figure7.png"

png(file=loc.out,width=800,height=800,pointsize=120,bg="white");

plot(combine,ann=FALSE,lab=c(7,3,12),las=1,type="l",mai=c(0.5,0.5,0,0))

dev.off();

#

write.table(combine,file="c:/climate/data/articles/figure6.txt",quote=FALSE,sep="\t")

##RECONSTRUCT NH AVERAGE

load("c:/climate/data/mann/eof.tab")

load("c:/climate/data/mann/tpca-eigenvals.tab")

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

load("c:/climate/data/mann/gridpoints.tab")

nhtemp<-(grid2[,4]>0)

load("c:/climate/data/mann/nil.tab")

load("c:/climate/data/mann/mann.scale.tab")

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.2<-ts.union(a,b)

plot.ts(combine.2,ann=FALSE)

loc.out<-"c:/climate/images/figure8.png"

png(file=loc.out,width=800,height=400,pointsize=1000,bg="white");

plot.ts(combine.2,ann=FALSE)

dev.off();

combine<-ts.union(a,b)

write.table(combine,file="c:/climate/data/articles/figure7.txt",quote=FALSE,sep="\t")