##REPLICATE AD200 PC1 in MANN AND JONES 2003

	load("d:/climate/data/tree/northamerica.details.tab") 
	write.table(details.old,file="d:/climate/data/tree/northamerica.details.dat",sep="\t",row.names=FALSE,quote=FALSE)
	details=read.table("d:/climate/data/tree/northamerica.details.dat",sep="\t",header=TRUE,fill=TRUE,quote="");nrow(details)


###IDENTIFY SITES IN PC1 STARTING IN 200
	mbhdata="http://www.climateaudit.org/data/MBH99"
		#url5=file.path(mbhdata,"TREE/ITRDB/NOAMER/BACKTO_1000")
	loc=file.path(mbhdata,"BACKTO_1000","noamer-itrdb-ad1000.txt")
	roster5=read.table(loc);nrow(roster5)##27 series in AD1000 roster
	strip= function(h,n) {
        	h<-unlist(strsplit(h, "\\."))
        	h<-array(h,dim=c(n,length(h)/n))
        	h<-t(h)
        	strip<-h[,1]
        	strip
        	}
	roster5=strip(as.character(roster5[,2]),2) #strips off .txt suffix

	details=read.table("http://www.climateaudit.org/data/tree/northamerica.details.dat",header=TRUE,fill=TRUE,sep="\t",quote="")
		##load("d:/climate/data/tree/northamerica.details.tab") 
		#this is my own collation of WDCP data and you'll have to bear with me on this
		#the short output from this calculation listed below is saved and is readable
	start0=200
	temp=(details$start<= start0)&!is.na(details$start);sum(temp) #31 series identified
		#temp=(details$start<1001)&!is.na(details$start);sum(temp) #31 series identified
	test=roster5[!is.na(match(roster5,details$id[temp]))];test #6 series isolated
		#"ca534" "ca535" "nm572" "nv515" "nv516" "ut509"  #Mann and Jones [2003] mention 6 series
	
	temp=!is.na(match(details$id,test))
	treeinfo=details[temp,1:5];treeinfo
	jones =function(lat,long) {
        	i<-18-floor(lat/5);
        	j<- 36+ceiling(long/5);
        	jones<-72*(i-1) + j
        	jones
        }
	treeinfo$jones=jones(treeinfo[,4],test1[,5]);test1$jones # 733 733 807 733 734 734
		##write.table(details[temp,1:5],file="d:/climate/data/osborn06/bristle.site.txt",sep="\t")
		#this is saved at http://www.climateaudit.org/data/osborn06/bristle.site.txt

	 #       id                     location type   lat    long jones
	#476  ca534               SHEEP MOUNTAIN PILO 37.22 -118.13   733
	#477  ca535              METHUSELAH WALK PILO 37.26 -118.10   733
	#2335 nm572 EL MALPAIS NATIONAL MONUMENT PSME 34.58 -108.06   807
	#2390 nv515                INDIAN GARDEN PILO 39.05 -115.26   733
	#2391 nv516                   HILL 10842 PILO 38.56 -114.14   734
	#2934 ut509                MAMMOTH CREEK PILO 37.39 -112.40   734


###COLLATE 6 SITES
     #collate site chronologies: I've saved R-tables for all WDCP chronologies and this is used here
     #this shows the collation; the collation itself is saved at climateaudit
	#combine=NULL
	#for (i in 1:nrow(test1)) {
	#	load(file.path("d:/climate/data/tree/chronologies/r/northamerica",paste(test1$id[i],"crn.tab",sep=".")))
	#	M=tsp(chron.crn)[1]
	#	chron.crn=ts(chron.crn[(start0-M+1):length(chron.crn)],start=start0)
	#	combine=ts.union(combine,chron.crn)
	#	}
	#combine=extend.persist(combine)
	#combine=cbind(tsp(combine)[1]:tsp(combine)[2],combine)	
	#dimnames(combine)[[2]]=c("year",as.character(test1$id))
	#dim(combine) # 1781    7
	#write.table(combine,file="d:/climate/data/osborn06/treecollation.AD200.txt",sep="\t",row.names=FALSE)
		#this is saved at http://www.climateaudit.org/data/osborn06/tree.collation.AD200.txt

	combine=read.table("http://www.climateaudit.org/data/osborn06/tree.collation.AD200.txt",sep="\t",fill=TRUE,header=TRUE)
	combine=ts(combine[,2:ncol(combine)],start=combine[1,1])
	combine[1,]
		#ca534 ca535 nm572 nv515 nv516 ut509 
		# 1071  1048  1087  1106   791  1257 

###CALCULATE PC1 USING MBH98 "METHOD"
 	options(digits=4)
	mannomatic= function(b,M=78,sdmethod="detrended") {
        	N<-nrow(b)
        	m2<-apply(b[(N-M):N,],2,mean)
        	s2<-apply(b[(N-M):N,],2,sd)
        	if (sdmethod=="unscaled") d<-scale(b,center= m2,scale=FALSE)
        	if (sdmethod=="undetrended")  d<-scale(b,center= m2, scale= apply(b[(N-M):N,],2,sd) )   
        	if (sdmethod=="detrended")  { d<-scale(b,center= m2, scale= apply(b[(N-M):N,],2,sd) );
        	        sdprox<-apply(d[(N-M):N,],2,sd.detrend)
                	d<-scale(d,center= FALSE,scale= sdprox)}
        	mannomatic<-d
        	mannomatic
            }

	combine.mannomatic=mannomatic(combine)
	pca.mannomatic=svd(combine.mannomatic)
	round(pca.mannomatic$v[,1],4)
		## 0.8950 0.1185 0.1297 0.2063 0.2749 0.2235
	pca.mannomatic$v[,1]^2
		## 0.80105 0.01404 0.01682 0.04256 0.07559 0.04994
		#the Sheep Mountain site has over 80% of the eigenvector weighting

	treeinfo$eof_mann=round(pca.mannomatic$v[,1],4); treeinfo
	barplot(treeinfo$eof_mann,names.arg=as.character(treeinfo$id),ylim=c(0,1));box()
		title(main="Mann and Jones 2003 PC1 Weights")

###
#     id        location type   lat   long jones annual decadal eof.mann
#1 ca534  SHEEP MOUNTAIN PILO 37.22 -118.1   733   0.03    0.06   0.8950
#2 ca535 METHUSELAH WALK PILO 37.26 -118.1   733  -0.14   -0.29   0.1185
#3 nm572     EL MALPAIS  PSME 34.58 -108.1   807  -0.28   -0.27   0.1297
#4 nv515   INDIAN GARDEN PILO 39.05 -115.3   733  -0.13   -0.44   0.2063
#5 nv516      HILL 10842 PILO 38.56 -114.1   734  -0.36   -0.48   0.2749
#6 ut509   MAMMOTH CREEK PILO 37.39 -112.4   734  -0.27   -0.75   0.2235