rundcf3<-function(x1,x2,plt=F,zz=2,zzz=2,clim=1,title = "", BP=1,tran1=1.5,tran2=1.5,spa=1/2,bin=30,mlag=150,nrun=2000,corr=0,p=0.95)
 {
   	### Based in part on script written by Williams and Post and based on Kirchner 
	### Does a geospatial equiv. to correlation
	##_________________________________________________________________________________
	##
	## x1 and x2 are original data frames for input
	## zz=2 and zzz=2 specify use of second column of data in the x1 and x2 dataframes as Y Data
	## spa = 1/2 specifies span for Loess smoother if data are transformed
	## bin = 0 tells program to used maximum median differences as choice of bin size
	## BP=0 is flag (used in DCF script) indicating time as years AD vs. BP=1 for years as before present
	## tran=1 standard tran and detrend, tran=0 no tran or detrend,tran=2 is no tran but gls for detrend,tran=-1 is detrend by differencing
	## mlag is the maximum lag desired and mlag/bin determines (+and-) the number of bins that lag calculated for
	## nrun = 1 if no bootstapped CI desired.  Else 1000 is recommended number of bootstraps.  Program bootstraps the y values and notthe original sample intervals
	## corr=0 is default for Pearson Correlation.  corr=1 for Spearman
	## p=0.95 is the default significance level for the bootstrapped CI
	##
	##_________________________________________________________________________________
	##
	##	Choice of transforms (with loess smooth) is based on smallest RSS after fitting the smooth to z(lambda) (see Venables & Ripley, p. 170)
	##_________________________________________________________________________________
	#
	# tya and tyb will start as in original matrix and then be replaced with residuals (if transformed)
	# tyta and tytb and these will be reshuffled if bootstrapping is chosen
	# xname1 and xname2 collect names of variables in x1 and x2 for use in plots
	# Program assumes that time1,val1,time2,val2 in order
	#
	##_________________________________________________________________________________
tempx1_as.matrix(x1)
tempx2_as.matrix(x2)
xname1_name.cols(x1)
xname2_name.cols(x2)
# txa and txb are the time vectors
# tya and tyb are the data vectors
		txa_as.vector(tempx1[,2][!is.na(tempx1[,zz])])
		txb_as.vector(tempx2[,2][!is.na(tempx2[,zzz])])
		tya_as.vector(tempx1[,zz][!is.na(tempx1[,zz])])
		tyb_as.vector(tempx2[,zzz][!is.na(tempx2[,zzz])])
# DCF and Program originally written assuming years BP (1950 is 0 BP and 1870 is 80BP)
# so if dates are AD need to adjust to get lags correct in DCF
#Does two loops to smooth/detrend the tya and tyb data
runs_2
a_0
for (k in 1:runs) {
if(k==1){tran_tran1}
if(k==2) {tran_tran2}
a_a+1
	if (a==1) {
		y <- tya
		year<-txa
	}
	if (a==2) {
		y <-tyb
		year<-txb
	}
if (tran ==1) {
	# Geometric mean
  yold <- y
  const <- min(y)
  if(min(y)<=0)  {y <- y-const+1}
  ydot <- exp(mean(log(y)))
	## No transform
	z <- y - 1
	fit <- loess(z ~ year, span = spa)
	detrend.y <- fit$resid
	detrend.rss <- sum(fit$resid^2)

	## Square root transform
	lambda <- 0.5
	z <- (y^lambda - 1)/(lambda * ydot^(lambda - 1))
	fit.sqrt <- loess(z ~ year, span = spa)
	detrend.sqrt.rss <- sum(fit.sqrt$resid^2)
	fit.sqrt <- loess(z * ydot^(lambda - 1) ~ year, span = spa)
	detrend.sqrt.y <- fit.sqrt$resid

	## Fourth root transform
	lambda <- 0.25
	z <- (y^lambda - 1)/(lambda * ydot^(lambda - 1))
	fit.4thrt <- loess(z ~ year, span = spa)
	detrend.4thrt.rss <- sum(fit.4thrt$resid^2)
	fit.4thrt <- loess(z * ydot^(lambda - 1) ~ year, span = spa)
	detrend.4thrt.y <- fit.4thrt$resid

	## Log transform 
	z <- log(y) * ydot
	fit.log <- loess(z ~ year, span = spa)
	detrend.log.rss <- sum(fit.log$resid^2)

	fit.log <- loess(z/ydot ~ year, span = spa)
	detrend.log.y <- fit.log$resid
	rss.best <- min(detrend.rss, detrend.sqrt.rss, detrend.4thrt.rss,
		detrend.log.rss)

	if(rss.best == detrend.rss) {
		zt <- detrend.y
		transform <- "None"
		fit <- fit
		trend <- fit$fitted.values + 1
	}	
	else if(rss.best == detrend.sqrt.rss) {
		zt <- detrend.sqrt.y
		transform <- "Square root"
		fit <- fit.sqrt
		trend <- (fit$fitted.values/2 + 1)^2
	}
	else if(rss.best == detrend.4thrt.rss) {
		zt <- detrend.4thrt.y
		transform <- "Fourth root"
		fit <- fit.4thrt
		trend <- ((fit$fitted.values)/4 + 1)^4
	}
	else if(rss.best == detrend.log.rss) {
		zt <- detrend.log.y
		transform <- "Natural log"
		fit <- fit.log
		trend <- exp(fit$fitted.values)
	}
## Bracket below closes tran loop to skip transforms
}
## No transformation or detrending
	if(tran==0) {
		zt<-y
		transform<-"No Transformations or Detrending"
		trend_0
	}
## Detrend by differencing.  Note that I divide difference by time step (this may be bogus).
	if(tran==-1) {
		nyear_(length(year)-1)
		zx_diff(year)
		zt_diff(y)
#		zt_zt/zx
		zt_zt[1:nyear]
		year_year[1:(nyear)]
		year_year[1:(nyear)]+(zx[1:(nyear)]/2)
		transform_"Detrend by Differencing"
		trend_zt
	}
## Detrend by linear fit and no transformations	
	if(tran==2) {
		yy_as.vector(y)
		yeary_as.vector(year)
		fit_gls(yy~yeary,na.action=na.omit)
		zt_residuals(fit)
		transform<-"Linear Regression"
		trend<-(fitted.values(fit))
	}
## Detrend by locfit.raw and no transformations	
	if(tran==1.5) {
		fit_locfit.raw(year,y,alpha=.2,deg=c(0,3))
		zt_residuals(fit)
		transform<-"Locfit"
		trend<-(fitted.values(fit))
}		
if (a==1) {tyta_zt}
if (a==1) {transform1_transform}
if (a==1) {
	if(tran1==-1){txaa_year}
	trend1_trend
#	if (min(yold)<=0) {trend1_trend1+const-1}
	}
if (a==2) {tytb_zt}
if (a==2) {transform2_transform}
if (a==2) {
	if(tran2==-1){txbb_year}
	trend2_trend
#	if (min(yold)<=0) {trend2_trend2+const-1}
	}
}
	
### Plot the data
	if(plt==T) {
		graphsheet(pages="every graph") 
		par(mfrow = c(4, 2))
		if(BP==0) lab = "Years AD"  
		if(BP==1) lab = "Years BP"
		plot(txa, tya, xlab = lab, ylab = xname1[zz])
		title(xname1[zz])
		if(tran1>0) {lines(txa, trend1)}
		if(tran1==-1){txa<-txaa}	
		plot(txa, tyta, type = "l", xlab = lab, ylab = "transformed density")
		plot(txb, tyb, xlab = lab, ylab = xname2[zzz])
		title(xname2[zzz])
		if(tran2>0) {lines(txb, trend2)}
	if(tran2==-1){txb<-txbb}
		plot(txb, tytb, type = "l", xlab = lab, ylab = "transformed density")
}
### Replaces original time points with points intermediate to observations that were differenced
	if(tran1==-1){txa<-txaa}		
	if(tran2==-1){txb<-txbb}
### Call the DCF Program to do actual analysis on transformed data
dcf1_as.data.frame(dcf5(BP,txa,tyta,txb,tytb,bin,mlag,nrun,p,corr))

binsize_dcf1[2,1]-dcf1[1,1]
dcf_sort.col(dcf1,c("bin","rbin","lower","upper","numinbin"),"bin")

	if(plt) {
		plot(dcf$bin,dcf$rbin,type="l",ylim=c((min(min(dcf$lower),min(dcf$rbin))),(max(max(dcf$upper),max(dcf$rbin)))),lty=1,xlab="lag(years)",ylab="corr")
		title("Lag plot")
		lines(dcf$bin,dcf$lower,lty=2)
		lines(dcf$bin,dcf$upper,lty=2)
		plot(0:1, 0:1, axes = F, type = "n", xlab = "", ylab = "")
		text(0.5, 0.8, paste("Title: ", title))
		text(0.5, 0.6, paste("Transform1: ", transform1,", Transform2: ", transform2))
		text(0.5, 0.4, paste("Bin Size: ", binsize,", Bootstrap Replicates: ",nrun,", P Value: ",p))
	}

### Calculate the DCF using the detrended series
	list(dcf = dcf)
}