#NHgamma
	##used in E&E 2005
	##this is to be used with config3
	#series - this is a table listing the number of PC series by 
	sd.rescaling<-rep(1,1082)
	gaspe.method<-"default"
	bristlecone.tag<-FALSE
	NOAMER.tag<-"null"
	series<-series.MBH
	directory<-directory.MBH
	pcsubdirectory<-"MBH98"

   NHgamma<-function(proxy,pcsubdirectory="MBH98",series=series.MBH,directory=directory.MBH,
		gaspe.method="archived",bristlecone.tag=FALSE,NOAMER.tag="null",sd.rescaling=rep(1,1082)) {
	if (gaspe.method=="archived") proxy[1:4,53]<-NA
	if (gaspe.method=="updated") { proxy[,53]<-NA; load("c:/climate/data/erren/gaspe.tab"); proxy[178:581,53]<-gaspe[1:404] }
	proxy0<-proxy

	np<-length(select) #number of steps
	mm<-period[1]-period[length(period)]
	G<-array(c(rep(NA,16*112)),dim=c(16,112))    #matrix of calibration coefficients
	D<- array ( rep(NA,np*112)  , dim=c(112,np)  ) #weighting factors returned as second item in list
	recon<-array(NA,dim=c(581,np)) #reconstruction stepwise
	recon.sparse<-array(NA,dim=c(581,np))  #sparse reconstruction stepwise
	H<-rep(list(NA),11)  #reconstructed RPCs
	recon.spliced<-rep(NA,581)
  
        #small function to do matrix algebra with different masks
  		f<-function(A,mask){ 	if (length(selection)>1 ) E<-diag( lambda[1:16][selection] )  else E<- as.matrix(lambda[1])
		f<- A %*% E  %*% t(eof) [selection,mask] %*% (mann.scale[mask]*sd.rescaling[mask])/sum(mu[mask])		
	
	f} 
      #small function to do matrix algebra with different masks
  		g<-function(selection,mask){ 	if (length(selection)>1 ) E<-diag( lambda[1:16][selection] )  else E<- as.matrix(lambda[1])
		g<- E %*% t(eof) [selection,mask] %*% (mann.scale[mask]*sd.rescaling[mask])/sum(mu[mask])		
		g} 

      #STEPWISE CALCULATIONS
	for (k in 1:11) {
		proxy<-proxy0

		#read in principal component series #tags for special studies on NOAM.CENSORED and PC4 effect
		for (m in 1:6)  {
			if (series[m,k]>0) {if ((m==3) & (k==11) & !(NOAMER.tag =="null")) txtid<-paste(NOAMER.tag,"tab",sep=".") else txtid<-"tab";
				load(file.path(pcstore,pcsubdirectory,paste(regionid[m],directory[m,k],txtid,sep="."))) } #this proides for using CENSORED version
			if (series[m,k]>0) { if (m==3 & bristlecone.tag) proxy[,(mannid[m]-1+(1:(series[m,k]-1)))]<-pca[,c(1:3,5:series[m,k])] else proxy[,(mannid[m]-1+(1:series[m,k]))]<-pca[,(1:series[m,k])] }						
			}
		m1<-apply(proxy[503:581,],2,mean,na.rm=TRUE);## means
		s1<-apply(proxy[503:581,],2,sd,na.rm=TRUE);##standard deviations
		for (j in 1:112) {
			temp<-!is.na(proxy[,j]);
			proxy[temp,j]<-(proxy[temp,j]-m1[j])/s1[j]
			}
		selection<-select[[k]]
		roster<-!is.na(proxy[period[k+1],]) #length 112
		
		U<-as.matrix(pc[1:79,selection])  #not standardized #standardized in older version
		Ut<-t(U)
  		B<-proxy[503:581,roster]
		G[selection,roster]<- solve(Ut %*% U ) %*% Ut %*% B;  #this is array of regression coefficients in calibration
		P<-diag(weights1[roster]^0.5)  #this is used in weighted regression
		y<-proxy [period[k+1]:(period[1]-1),roster] #stepwise proxy rosters

	#calculate weighting factors	
		C <- P %*%  ( t(G[selection,roster] %*% P) %*% solve ((G[selection,roster] %*% P) %*% t((G[selection,roster] %*%P)))   )
			#dimension (#proxies in step x #TPCs in step e.g. 101 11)

	#calculate re-scaling factors for reconstructed temperature principal components
		#this is a very slight non-linearity 
		H[[k]]<-y %*% C 
		scale.est<-apply(as.matrix(H[[k]] [(1902-period.m[k+1]+1):(1980-period.m[k+1]+1),]),2,sd) #adjusted May 2005
		scale.adj<-scale.est/apply(as.matrix(pc[1:79,selection]),2,sd)
		if (length(selection)>1) E<-diag(1/scale.adj) else E<-as.matrix(1/scale.adj)
		H[[k]]<- H[[k]] %*% E #this matches sd as per scaling
				#dimension (#years in step x #TPCs in step e.g. 181 11)

    	#do weighting of RPC
		weights.nh<- f(C%*% E,nhtemp.nil)	#C %*% g(selection,nhtemp.nil)
			### y%*% weights.nh-recontemp[period[k+1]:(period[1]-1)]  #this is exactly equal
		recon[period[k+1]:(period[1]-1),k]<- y %*% weights.nh #edited May 2005 #simplified Jan 2006
		recon.sparse[period[k+1]:(period[1]-1),k]<- y %*% f( C %*% E,mask.logical&nhtemp.nil) 	#y %*% C %*% g(selection,mask.logical&nhtemp.nil)	
		recon.spliced[period[k+1]:(period[k]-1)]<-recon[period[k+1]:(period[k]-1),k]

	#save weighting factor for NH temperature construction
		 D[roster,k]<- weights.nh
		}

    ##COLLECT TERMS
	tags<-list(case,gaspe.method,bristlecone.tag,NOAMER.tag,series,directory,sd.rescaling)
	names(tags)<-c("case","gaspe.method","bristlecone.tag","NOAMER.tag","series","directory","sd.rescaling")
	NHgamma<-list(ts(recon.spliced,start=1400),D,ts(recon[,11],start=1400),ts(recon[,10],start=1400),ts(recon,start=1400),ts(recon.sparse,start=1400),H,proxy,tags)
	names(NHgamma)<-c("spliced","weights","AD1400","AD1450","recon","recon.sparse","RPC","proxy","tags")
	NHgamma
	} 


H1<-H
recon1<-recon
D1<-D


##README
 	#much of the interest has developed from different NOAMER PC series
	#these can be studied by separately creating the PC series in the "UVA" or "WDCP" directories
	#the default tag is NOAMER.1400.tab, etc.
	#the variations are tagged NOAMER.1400.CENSORED.tab; NOAMER.1400.nogaspe.tab, etc...

	# spliced - this is time series spliced from each step
	# weights - this is a matrix of weights for each proxy in each step col 1- 1820,..., col 11- 1400
	#AD1400   - reconstruction using AD1400 proxies
	#AD1450  - reconstruction using AD1450 proxies
	#recon	- 581x11 matrix of "dense" reconstructions for steps 	col 1- 1820,..., col 11- 1400
	# recon.sparse - 581x11 matrix of "sparse" reconstructions for steps 	col 1- 1820,..., col 11- 1400