## Esper et al 2003
## Emulate Jaemtland Comparison 
#########

source("http://www.climateaudit.org/scripts/utilities.txt") #
f=function(x) filter.combine.pad(x)[,2]

source("http://www.climateaudit.org/scripts/tree/utilities.treering.txt")
source("http://www.climateaudit.org/scripts/yamal/collation.functions.txt")
source("http://www.climateaudit.org/scripts/yamal/analysis.functions.txt")
source("d:/climate/scripts/yamal/analysis.functions.txt")

restate=function(tree,info) { 
	tree$tree=  info$tree[match(tree$id,info$id)] 
 	tree$age=factor(tree$tree)
	levels(tree$age)=tapply(tree$year,tree$tree,min) 
	tree$age=as.numeric(as.character(tree$age))
	tree$age=tree$year-tree$age+1	
	return(tree)
   }


###############
##LOAD DATA
##################
    #the counts match reported counts in Esper 2003
	download.file( "http://www.climateaudit.org/data/esper/jae.rwl.tab","temp.dat",mode="wb");load("temp.dat")
	jae=tree;dim(tree);dim(jae)
	  #  25249     4
	info=info.jae=make.info(tree)
	info$site=factor(substr(info$id,1,3))
	tapply(info$core,info$site,function(x) length(unique(x) ))
	 #545 940 TOR 
	 #  2  10   1 
	  	 #matches Esper 2003: jae: 133= 119 PISYjae + 14 PCAB jae; sk: 24 PISYsk; Tornertrask - 1
	temp= info$site=="940"
	info$tree[temp]=info$id[temp]; 
	info$core[temp]="1"
	nrow(info) # 158 
	jae=restate(jae,info)


##################
#STRATIFY POPULATION
########################
 ##Stratify the population by making Id list

    strata=list()
	strata[[1]]= info$id[info$site=="940"]	#swed023
	strata[[2]]= info$id[info$site=="545"]	#norw002
	strata[[3]]= info$id[info$site=="TOR"]  #tornetrask
	names(strata)=c("swed023","norw003","tornetrask")


###############################
# CALCULATE COUNTS AND RCS CHRONOLOGIES (nls)
###############################

	count.jae=make.count(jae,strata) 
	chron.jae= make.chron(jae,strata[1:2]) #this collects chron for each stratum into a litst
	sapply(strata, function(A) max(tree$age[!is.na(match(tree$id,A))]) );
	 #  swed023    norw003 tornetrask 
       	 #     276        203        201 

	data.frame(t(round(sapply(chron.jae$strata, function(A) A$coefficients) ,4) ) )
	
         #      A        B      C
	#swed023  33.1308 122.7716 0.0340
	#norw003 -74.6444 242.6716 0.0038

#################
## PLOTS
####################

### COUNTS

  # GDD(file="d:/climate/images/2009/esper/jaemtland_count.gif", type="gif", h=360, w=420)
	plot.count(count.jae)
	title ("Esper JAE Core Counts")
    # dev.off()


###  Age Dependence Curve of Subpopulations

#   GDD(file="d:/climate/images/2009/esper/jaemtland_dependence_rev.gif", type="gif", h=360, w=360)
	plot.rcs(chron.jae,ylim0=c(0,250),method="biweight" )

#   dev.off()

#    GDD(file="d:/climate/images/2009/esper/jaemtland_composite.gif", type="gif", h=320, w=420)
	details=Details=list(at=seq(25,75,25),ylim=c(0,150),col=rep(c("grey40","grey80"),4) )
	plot.composite(chron.jae[2:3],count.jae,Details)
	abline(v=c(1120,1283,1340),col=2,lty=3)
	abline(v=c(1789,1828),col="purple",lty=3) # 20%,80% quantiles
	title("Jaemtland Populations")

#	dev.off()


#####
## MISCELLANEOUS 
###########
	download.file( "http://www.climateaudit.org/data/esper/esper.archived.tab","temp.dat",mode="wb");
	load("temp.dat")
	range( time(esper.archived)[!is.na(esper.archived[,"Jae"])] )
		 #1107 1804