sp.plot <- function(sp, defree, sp.bootind, h, d, interval, conf = 0.95,
add = F, col = 2, cex = 2)
{
# sp.plot:
# This function wraps everything up and performs some of the 
# "invisible" calculations (e.g. the second derivative estimates and
# changepoints).
# It plots the estimated abundance indices and confidence intervals, and
# marks significant upturns and downturns (i.e. "changepoints") on the
# curve in open and closed circles. 
# 
# "sp" is the name of the species data frame, and must be a character
#    string: e.g. "cb". 
#    It is assumed that the R object "indsp" (e.g. indcb) containing
#    the GAM fits to the original data for the required df has been
#    obtained previously (using function indsp.func), and is available
#    BY THAT NAME on the permanent database. 
#    If it is not there, then form it first at the command line by typing:
#    indsp <- indsp.func(sp, dfvec=defree)
#    (e.g. indcb <- indsp.func(cb, defree) where defree is replaced by
#    the required value). 
# "defree" is the df of the required GAM fit: e.g. if defree=10 then
#    the indices indsp[,"df10"] will be extracted for the main curve. 
# "sp.bootind" is a matrix of bootstrapped indices, e.g.
#    cb.bootind.399.10.  It should be the relevant object for the GAM
#    with specified df (i.e. if df=10 then sp.bootind should be
#    sp.bootind.Nrep.10 for some Nrep).
#   
# See function diffop.func for explanation of "h", "d", and "interval". 
# "conf" gives confidence level for the ABUNDANCE INDICES confidence
#    intervals: *not* the confidence level for the changepoints.
#    The changepoints are automatically calculated at the 5% significance
#    level. This can be changed if required by manually editing the value
#    "deriv.conf" in the function body. 
# "add" determines whether to add the plot to a current plot, or
#    whether to start a new plot.  The default is to start a new plot. 
# "col": if add=T, col determines what colour the new plot (superposed on 
#    the existing plot) should be plotted in.
# "cex" determines the size of the text and changepoints on the plot:
#    increase it for larger text, decrease it for smaller text. 
#
# EXAMPLE: 
# # plot for 10 df:
#      sp.plot("cb", 10, cb.bootind.399.10, 1, 6, 1, conf=0.95)
# # then superpose the plot for 20 df: 
#      sp.plot("cb", 20, cb.bootind.399.20, 1, 6, 1, conf=0.95, add=T)
#
	deriv.conf <- 0.95
	par(font = 3)	#
# housekeeping:
	if(!is.character(sp))
		stop("first argument must be a character string.")
	dfname <- paste("df", defree, sep = "")
	sp.ind.ci <- index.ci.func(sp.bootind, conf = conf)
	Nyears <- length(sp.ind.ci[1,  ])
	yearvec <- seq(1, Nyears)	#
# Might want to change yearvec to the actual years of the survey: e.g.
# yearvec <- seq(1962, 1995) instead of seq(1, 34).  
#
	lower <- min(sp.ind.ci[1,  ])
	upper <- max(sp.ind.ci[2,  ])
	indsp <- eval(parse(text = paste("ind", sp, sep = "")))    
	if(!add) {
		plot(yearvec, indsp[, dfname], type = "l", ylim = c(lower, 
			upper), xlab = "", ylab = "", cex = cex, lwd = 2, bty
			 = "l", font = 3)
		textint <- 0.2/1.7 * (upper - lower)	#
# textint is for beautification and nothing else
		text(min(yearvec) - 1, upper + textint, "index", cex = cex, 
			font = 3)
		text(max(yearvec), lower - textint, "year", cex = cex, font = 3
			)
	}
	if(add) {
		par(col = col)
		lines(yearvec, indsp[, dfname], lwd = 2)
	}
	lines(yearvec, sp.ind.ci[1,  ], lty = 4, lwd = 2)
	lines(yearvec, sp.ind.ci[2,  ], lty = 4, lwd = 2)	#
# Now calculate matrix of bootstrapped second derivatives:
	sp.bootderiv <- indmat.derivmat.func(sp.bootind = sp.bootind, h = h, d
		 = d, interval = interval)
	sp.deriv.ci <- deriv.ci.func(sp.bootderiv, conf = deriv.conf)	#
# Calculate the changepoints:
# the lower changepoints are those where the *upturn* of the curve
# is statistically significant (i.e. the lower confidence limit is 
# greater than 0);
# the upper changepoints are those where the *downturn* of the curve
# is statistically significant (i.e. the upper confidence limit is 
# less than 0);
        changepts.lower <- seq(1, Nyears)[sp.deriv.ci["lower",  ] > 0]
        changepts.upper <- seq(1, Nyears)[sp.deriv.ci["upper",  ] < 0]
        if((length(changepts.lower > 0) & (!is.na(changepts.lower[1])))) {
                points(changepts.lower + min(yearvec) - 1, indsp[   
                        changepts.lower, dfname], pch = 1, cex = cex)
                changepts.lower <- changepts.lower + min(yearvec) - 1
        }
        else changepts.lower <- "none"
        if((length(changepts.upper > 0) & (!is.na(changepts.upper[1])))) {
                points(changepts.upper + min(yearvec) - 1, indsp[
                        changepts.upper, dfname], pch = 16, cex = cex)
                changepts.upper <- changepts.upper + min(yearvec) - 1
        }
        else changepts.upper <- "none"
        par(col = 1)
        list(upturns = changepts.lower, downturns = changepts.upper)
}