#utility function to count by strata and year
make.count= function(tree,strata) {  
      	K=length(strata)
	count=NULL
	for ( i in 1:K) {
	  temp=!is.na(match(tree$id,strata[[i]]));sum(temp) #21362
	  count= ts.union(count,countf(tree[temp,]) )
	}
	dimnames(count)[[2]]=names(strata)
	count[is.na(count)]=0
	return(count)
}


#utility function to collate RCS chronologies by strata
 make.chron=function(tree,strata) {
    	K=length(strata)
	chron=list()
	for ( i in 1:K) {
	  temp=!is.na(match(tree$id,strata[[i]]));sum(temp) #21362
	  chron[[i]]= RCS.chronology(tree[temp,],method="nls")
	}
	names(chron)=names(strata)
	composite= sort( unique(c(unlist(strata) )))
	temp=!is.na(match(tree$id,composite));sum(temp) #21362
	chron.combo = RCS.chronology(tree[temp,],method="nls")
	return(list(strata=chron,combo=chron.combo) )
}

plot.count=function(count) {
	par(mar=c(3,4,2,1))
	year=c(time(count));N=length(year)
	x=apply(count,1,sum,na.rm=T)
	plot(year,x,type="l")
	polygon(xy.coords(x=c(year,rev(year)),y=c(x,rep(0,N))),col="red",border="red")
	x=apply(count[,1:2],1,sum,na.rm=T)
	polygon(xy.coords(x=c(year,rev(year)),y=c(x,rep(0,N))),col="cyan",border="cyan")
	x=count[,1]
	polygon(xy.coords(x=c(year,rev(year)),y=c(x,rep(0,N))),col="green4",border="green4")
	legend("topleft",fill=c("green4","cyan","red"),legend=dimnames(count)[[2]]) 
}  

plot.rcs=function(chron,ylim0=NULL,method="nls" ,plotsmooth=FALSE) {
 #chron=chron.jae # first series is combo
	Coef=data.frame(t(round(sapply(chron$strata, function(A) A$coefficients) ,4) ) )
	 	##N= sapply(strata, function(A) max(tree$age[!is.na(match(tree$id,A))]) )
	N =  sapply(chron$strata, function(A) length(A$mean) )
	if(is.null(ylim0)) ylim0=c(0,1.1* max(sapply(chron, function(A) max(A$mean) )) )	
		K=nrow(Coef)
	xlim0= c(0,ceiling(max(N)/100)*100)  
	col0=c("grey80","salmon","limegreen","lightblue")
	par(mar=c(4,4,2,1))
      if(method=="nls") {
	age=1:N[1]	
	i=1; plot(age,chron$strata[[i]]$mean ,type="l",ylab="Ring Width",xlab="age" ,
	  main="",lwd=1,ylim=ylim0,col="grey60",lty=i,yaxs="i",xlim=xlim0)
	  if(plotsmooth) lines(age,(Coef$A[i]+Coef$B[i]*exp(-Coef$C[i]*age)),lwd=2,col=i )  
	for(i in 2:K) {
	  age=1:N[i]; 
	  lines(age,chron$strata[[i]]$mean,lwd=1,col=col0[i] )  
	  if(plotsmooth) lines(age,(Coef$A[i]+Coef$B[i]*exp(-Coef$C[i]*age)) ,col=i,lty=2,lwd=2)
	  }
        } #nls
      if(method=="biweight") {
	i=1; x= chron$strata[[i]]$biweight.mean ;age=1:length(x)	
	plot(age,x,type="l",ylab="Ring Width",xlab="age" ,
	  main="",lwd=1,ylim=ylim0,col=i,lty=i,yaxs="i",xlim=xlim0)
	  fm <- nls(x ~ A+B*exp(-C*age),data = data.frame(age,x),
		 start = list(  A=mean(x,na.rm=T)/4,B = mean(x,na.rm=T), C= .01 ),
		 trace = TRUE,control=nls.control(maxiter=200, tol=1e-02, minFactor=1e-10)); 
  	  B<-coef(fm); 
  	  if(plotsmooth) lines(age,(B[1]+B[2]*exp(-B[3]*age)),lwd=1,col=i )  
	for(i in 2:K) {
	  x= chron$strata[[i]]$biweight.mean ;age=1:length(x)	
	  lines (age,x,col=i,lty=i)
	  fm <- nls(x ~ A+B*exp(-C*age),data = data.frame(age,x),
		 start = list(  A=mean(x,na.rm=T)/4,B = mean(x,na.rm=T), C= .01 ),
		 trace = TRUE,control=nls.control(maxiter=200, tol=1e-02, minFactor=1e-10)); 
  	  B<-coef(fm); 
  	  if(plotsmooth) lines(age,(B[1]+B[2]*exp(-B[3]*age)),lwd=1,col=i )  
	  }
        } #biweight
       legend("topright",fill=1:K,legend=names(chron$strata),cex=.8)
}

plot.composite=function( chron,count,details=Details) {
	par(mar=c(3,4,2,4))
 	year=c(time(count));N=length(year)
 
	x=apply(count,1,sum,na.rm=T);maxcount=max(x)
	ylim0= c(0,10* ceiling(maxcount/10) *2 )
 	plot(year,x,type="l",ylim=ylim0,ylab="",axes=FALSE)
	axis(side=1); axis(side=4,las=1,at=details$at,labels=details$lat);
	mtext(side=4,"Count",line=1.8,font=2)
 	polygon(xy.coords(x=c(year,rev(year)),y=c(x,rep(0,N))),col=details$col[1],border=details$col[1])
 	x= count[,1]
 	polygon(xy.coords(x=c(year,rev(year)),y=c(x,rep(0,N))),col=details$col[2],border=details$col[2])

 	x= chron$combo$series;
	A= data.frame( x=range(x),y= c(1,2)*10* ceiling(maxcount/10) )
	fm=lm(y~x,data=A)
	h=function(x) fm$coef[1]+fm$coef[2]*x  
	lines(c(time(x)), h(x),col="grey80" ) 
	lines(c(time(x)), h(f(x)),col=1,lwd=2 ) 
	axis(side=2,las=1,at=h(seq(0,2,.5)),labels= seq(0,2,.5) )
	box()
	mtext(side=2,"Chronology",line=2.4,font=2)
 	abline(h=h (seq(0,2.5,.5)),col="grey40",lty=3)
}