".First" <-
function() 
{
Sys.putenv("TCL_LIBRARY"="c:\\Program Files\\Tcl\\lib\\tcl8.3")
library(chron)
tts.gui()
q("no")
}
"read.campbell" <-
function(file,path="C:\\newdata")
{
	# Reads a Campbell data file and creates a data frame
	data <- read.table(paste(path,file,sep="\\"), sep = ",")[, 2:10]
	names(data) <- c("year", "day", "time", "dump", "bottle", "code1",
		"code2", "stg", "turb")
	data$chr <- campbell.date(data$year, data$day, data$time)
	nc <- dim(data)[[2]]
	data[, c(nc, 5:(nc - 1))]
}
"get.thresh" <-
function(tstart, tmax, n, f = sqrt, finv = function(x)
x^2)
{
	# Returns a vector of thresholds 
	# Vector starts with 0 and tstart, and ends with tmax
	dx <- (f(tmax) - f(tstart))/(n - 1)
	x <- seq(from = f(tstart), by = dx, length = n)
	c(0, round(finv(x)))
}

"get.upthresh" <- function(...)  c(get.thresh(...),9999)
"get.dnthresh" <- function(...) get.thresh(...)

"campbell.date" <-
function(year, day, time = 0)
{
	# Returns chron array for given year, day, and time vectors
	# year = 4-digit year
	# day = day number from start of year (Jan.1 = 1)
	# time = military time
	# Sets common origin to 1/1/1970 for axis.time to work properly
	daynum <- numeric(length(year))
	for(yr in unique(year)) {
		tmp <- dates(day[year == yr], origin = c(12, 31, yr - 1))
		daynum[year == yr] <- as.numeric(dates(tmp, origin = c(1, 1,
			1970)))
	}
	chron <- dates(daynum) + mt2msm(time)/1440.
	chron
}

"tts.sample" <-
function(data, up, dn, stay = 2, revpct = c(10, 20), repwait = 8, minstg = 0, turblim = 2000, limskip = 2, initcode = 3, revval = 5, basewait = 18)
{
	#  Rev  2.2  JL  051103 added revval and basewait to argument list
        #       Fixed bug where last.rsam was not being set for startup samples
	#       Changed baseflow condition to exclude equality condition (stg=minstg)
	#  Rev  2.1  JL  051028 changed initcode default to 3 (rising status undetermined)
	#       With default initcode, rising status is determined at second interval
	#  Rev  2.0  JL  051027 disabled lowest rising and falling thresholds (normally zero).  
	#       Treatment of zero rising threshold differs from Campbell code which seems to 
	#       be able to sample the zero threshold but rarely does when stg > minstg.
	#  Rev  1.9  JL  041027 added minimum stage as an attribute to the returned data frame
	#  Rev  1.8  JL  040323 moved the "switch to rising" code snippet out of the "if"
	#       block within the falling turbidity "saturation" code
	#  Rev  1.7  JL  040322 made 3 changes in "saturation" logic to match pseudocode
	#     - Comparisons for saturation turbidity include equality (t >= turblim)
	#     - Checks for saturation turbidity during falling condition.  At Rev 1.4,
	#            this block had been mistakenly put in the 'rising' while block
	#     - Suppresses rising-to-falling reversal count when saturated 
	#  Rev  1.6  JL  040130 added initcode option to allow starting in falling condition
	#     And no longer gives error if no samples are produced in the simulation
	#  Rev  1.5  JL  031218 fixed bug: tmax had not been set when coming out of baseflow
	#  Rev  1.4  JL  031212 updated to include overflow sampling
        #  Rev  1.3  JL  030116 updated to TTS Rev 4.0 sampling logic
        #     Now mimics Campbell TTS sampling exactly when not saturated
        #     Returns indexes to samples and list of threshold codes
        #  Rev. 1.2  JL  950807 revroom parameter dropped, 
        #     was no longer needed after Rev 1.1.
        #  Rev. 1.1  JL  950804
        #     This version samples after a reversal if a threshold
        #     has been crossed since the prior peak or trough, 
        #     regardless of revroom.
        #  data is a data frame with stg and turb columns
        #  up: vector of ascending thresholds for rising turbidities
        #  dn: vector of ascending thresholds for falling turbidities
        #  stay: number of intervals required before a threshold    
        #    or reversal condition is deemed to be met
        #  revpct: percent change required for detection of reversal
        #    revpct[1] is for peaks and revpct[2] is for troughs.
	#  revval: absolute change required for detection of reversal
        #    the maximum of the revpct and revval criterion is used
        #  repwait: number of intervals before reusing a threshold
        #  minstg: minimum stage for sampling
        #  turblim: sensor limit
        #  limskip: number of intervals skipped between samples when turblim exceeded
        #  initcode: initial value for tcode
	#      initcode=1 starts simulation in rising condition
	#      initcode=2 starts simulation in falling condition
	#      initcode=3 (default) starts simulation in unknown condition
        #  
        #  i: index to turbidity vector
        #  j: index to next turbidity threshold
        #  thrcount: number of intervals meeting threshold condition
        #  revcount: number of intervals meeting reversal condition
        #  tcode: same as thrcode in TTS, but equals 4 for startup sample
        #  thrcode: vector of tcode values
        sam <- numeric(0)
        thrcode <- numeric(0)
        thrcount <- revcount <- 0
        last.rsam <- last.fsam <-  - repwait
        last.rthr <- last.fthr <- -1
        tmax <- tmin <- data$turb[1]
	tmin <- tmax + 1
	tcode <- initcode
	j <- 1
        if (tcode == 1) 
		j <- 1 + sum(data$turb[1] >= up)
	else if (tcode == 2)
		j <- max(1,sum(data$turb[1] > dn))
        i <- 0
        imax <- length(data$turb)
        while(i < imax) {
                while(data$stg[i + 1] < minstg) {
                        # Process baseflow data
                        tcode <- 0
                        i <- i + 1
                        if(i == imax)
                                break
                }
                if(tcode == 0 && i < imax) {
                        # Emerging from baseflow
                        if(data$turb[i + 1] >= up[2] && (last.rsam == 0 || (
                                i - basewait) > max(last.rsam, last.fsam))) {
                                # Startup sample
                                i <- i + 1
				last.rsam <- i
                                sam <- c(sam, i)
                                thrcode <- c(thrcode, 4)
				j <- 1 + sum(data$turb[i] >= up)
                        }
                        tcode <- 1
			tmax <- data$turb[i]
                }
		else if (tcode == 3) {
			# Rising condition unknown (first or second interval only)		
			if (i == 0)
			   i <- i + 1
			else if (i == 1) {
			    # note: we're at interval 2 (i has not yet been incremented)
			    # compare first two intervals to set turbidity condition			
			   if (data$turb[1] < data$turb[2]) {
				tcode <- 1
				j <- 1 + sum(data$turb[2] >= up)
			   }
			   else {
				tcode <- 2
				j <- max(1,sum(data$turb[2] > dn))
			   }
			}
		}
                while(tcode == 1 && i < imax) {
                        # Process rising turbidities
                        if(data$stg[i + 1] < minstg) break
                        i <- i + 1
                        t <- data$turb[i]
                        if(t >= up[j]) {
                                # Next rising threshold reached
                                thrcount <- thrcount + 1
                                if(thrcount >= stay) {
					if (j > 1) {
					     # Threshold condition satisfied
                                             if(up[j] != last.rthr || i > last.rsam +
                                                repwait) {
                                                # Take a sample
                                                sam <- c(sam, i)
                                                thrcode <- c(thrcode, tcode)
                                                last.rthr <- up[j]
                                                last.rsam <- i
                                             }
					}
                                        j <- 1 + sum(t >= up)
                                        thrcount <- 0
                                }
                        }
			else if (t >= turblim) {
				if (i > max(last.rsam, last.fsam) + limskip) {
					# Collect an 'overflow' sample
					sam <- c(sam, i)
					thrcode <- c(thrcode, 5)
					last.rsam <- i
					last.rthr <- 9999
					revcount <- 0
				}
			}
                        if(t >= tmax) {
                                # Highest turbidity so far
                                tmax <- t
                                revcount <- 0
                        }
                        else if(t < data$turb[i - 1] && t < turblim) {
                                # Turbidity has started to fall
                                revcount <- revcount + 1
                                revthresh <- tmax - max(revval, 0.01 * revpct[1] *
                                        tmax)
                                if(t <= revthresh && revcount >= stay) {
                                        # Turbidity has definitely reversed
                                        tcode <- 2
                                        j <- max(1,sum(t > dn))
                                        tmin <- t
					if (j < length(dn)) {
                                        	fthr <- dn[j + 1]
                                        	if(tmax > fthr && (i > (last.fsam +
                                                	repwait) || fthr != last.fthr)) {
                                                	# Take a sample
                                                	sam <- c(sam, i)
                                                	thrcode <- c(thrcode, tcode)
                                                	last.fthr <- fthr
                                                	last.fsam <- i
						}
                                        }
				}
                        }
                }
                thrcount <- 0
                revcount <- 0
                while(tcode == 2 && (i < imax)) {
                        # Process receding turbidity
                        if(data$stg[i + 1] < minstg) break
                        i <- i + 1
                        t <- data$turb[i]
                        if(t <= dn[j] && j > 1) {
                                # Next falling threshold reached 
                                thrcount <- thrcount + 1
                                if(thrcount >= stay) {
                                        # Threshold condition satisfied
                                        if(dn[j] != last.fthr || i > last.fsam +
                                                repwait) {
                                                # Take a sample
                                                sam <- c(sam, i)
                                                thrcode <- c(thrcode, tcode)
                                                last.fthr <- dn[j]
                                                last.fsam <- i
                                        }
                                        j <- max(1,sum(t > dn))
                                        thrcount <- 0
                                }
                        }
                        if(t <= tmin) {
                                # Lowest turbidity this recession
                                tmin <- t
                                revcount <- 0
                        }
                        else if(t > data$turb[i - 1] && t < turblim) {
                                # Turbidity has started to rise again
                                revcount <- revcount + 1
                                revthresh <- tmin + max(revval, 0.01 * revpct[2] *
                                        tmin)
                                if(t >= revthresh && revcount >= stay) {
                                        # Turbidity has definitely reversed
                                        tcode <- 1
                                        j <- 1 + sum(t >= up)
                                        tmax <- t
                                        rthr <- up[j - 1]
                                        if(tmin < rthr && (i > (last.rsam +
                                                repwait) || rthr != last.rthr)
                                                ) {
                                                # Take a sample
                                                sam <- c(sam, i)
                                                thrcode <- c(thrcode, tcode)
                                                last.rthr <- rthr
                                                last.rsam <- i
                                        }
                                }
                        }
			else if (t >= turblim) {
				# Switch to rising
				tcode <- 1
				tmax <- t
				j <- 1 + sum(t >= up)
                                thrcount <- 0
				if (i > max(last.rsam, last.fsam) + limskip) {
					# Collect an 'overflow' sample
					sam <- c(sam, i)
					thrcode <- c(thrcode, 5)
					last.rsam <- i
					last.rthr <- 9999
					revcount <- 0
				}
			}
                }
                thrcount <- 0
                revcount <- 0
        }
        	result <- data.frame(sam, thrcode)
		attr(result,"minstg") <- minstg
		result
}

"plot.sim" <-
function(df, simout, sdate, stime, edate, etime, ylim, title = "", numbers=F)
{
#       Revised 051104 J.Lewis   added numbers argument; added legend
#         Previous versions did not plot symbols for sample code = 5
#	Revised 041027 J.Lewis   plots data above minimum stage in red
#		(minimum stage is obtained as an attribute of simout)
#	  Plots a TTS sampling simulation
#	  df = data frame created by read.campbell() or read.flo() function
# 	  simout = object created by tts.sample() simulation function
#	  sdate = optional start date (yymmdd)
#	  stime = optional start time (hhmm)
# 	  edate = optional end date (yymmdd)
#	  etime = optional end time (hhmm)
# 	  ylim = optional ylimits, vector of length 2, e.g. c(0,500)
#	  title = optional title for plot in quotes
#         numbers = T or F; if T then show bottle numbers; but can be VERY slow
        if(missing(sdate)) {
                schron <- min(df$chr)
                sdate <- trunc(schron)
                stime <- msm2mt(1440 * (schron - sdate))
        }
        else {
                if(missing(stime))
                        stime <- 0
                schron <- make.chron(sdate, stime)
                # fractional date
                sdate <- trunc(schron)
        }
        if(missing(edate)) {
                echron <- max(df$chr)
                edate <- trunc(echron)
                etime <- msm2mt(1440 * (echron - edate))
        }
        else {
                if(missing(etime))
                        etime <- 2400
                echron <- make.chron(edate, etime)
                # fractional date
                edate <- trunc(echron)
        }
	row.names(df) <- 1:dim(df)[1]
        df <- df[df$chr >= schron & df$chr <= echron,  ]
        nam <- as.numeric(row.names(df))
        attach(df)
        on.exit(detach(2))
        if (missing(ylim)) ylim _ c(0,max(turb))
        plot(chr, turb, xlab = "", axes = F, ylim=ylim, type = "l")
	turb2 <- turb
	turb2[stg<attr(simout,"minstg")] <- NA
	lines(chr, turb2, col=2)
        axis(2)
        axis.time(chr)
        box()
        t <- simout
	set0 <- which(nam %in% t[, 1])
        set1 <- which(nam %in% t[t[, 2] == 1, 1])
        set2 <- which(nam %in% t[t[, 2] == 2, 1])
        set4 <- which(nam %in% t[t[, 2] == 4, 1])
	set5 <- which(nam %in% t[t[, 2] == 5, 1])
        if(length(set1) > 0)
            points(chr[set1], turb[set1], pch = 2, cex = 0.8)
        if(length(set2) > 0)
            points(chr[set2], turb[set2], pch = 6, cex = 0.8)
        if(length(set4) > 0)
            points(chr[set4], turb[set4], pch = 1, cex = 0.8)
        if(length(set5) > 0)
            points(chr[set5], turb[set5], pch = 16, cex = 0.8)
        if(numbers && length(set0) > 0) {
            dx <- 0.015 * (par()$usr[2.] - par()$usr[1.])
	    botnum <- which(t[, 1] %in% nam)
            text(chr[set0] + dx, turb[set0], botnum, cex = 0.7)
	}
	usr <- par()$usr
	x0 <- .87*usr[2] + .13*usr[1]
	y0 <- .985*usr[4] + .015*usr[3]
	legend(x0,y0,c("rising","falling","startup","fixtime"),pch=c(2,6,1,16),cex=0.75)
        title(title)
}

"make.chron" <-
function(date, time, out = "day-mon-year", origin = c(1., 1., 1970.))
{
	# Makes a chron object from a date and military time
	# Date and time can be character or numeric vectors
	# out is the format for the date only
	ymd2date(date, out = out, origin = origin) + mt2msm(time)/1440.
}
"msm2mt" <-
function(msm)
{
	hr <- msm %/% 60.
	mn <- msm %% 60.
	round(100. * hr + mn)
}
"ymd2date" <-
function(date, out = "day-mon-year", origin = c(1., 1., 1970.))
{
	# Converts numeric or character date (yymmdd) to "dates" object
	date <- zfill(date, 6.)
	dates(date, format = "ymd", out = out, origin = origin)
}
"up2000" <-
c(0, 20, 77, 170, 300, 467, 670, 910, 1187, 1500, 1850, 9999)
"dn2000" <-
c(0, 30, 62, 105, 159, 225, 302, 391, 491, 602, 724, 858, 1004, 
1160, 1328, 1507, 1698, 1900)
"up1000" <-
c(0, 15, 46, 94, 158, 240, 338, 453, 585, 734, 900)
"dn1000" <-
c(0, 20, 37, 59, 86, 118, 155, 197, 245, 297, 355, 417, 485, 
558, 636, 719, 807, 900)
"read.flo" <-
function(stn, hy, path = ".")
{
	file <- paste(stn, zfill(hy, 2), ".flo", sep = "")
	df <- read.table(paste(path,file,sep="\\"), sep = ",")
	names(df) <- c("year", "mo", "dy", "time", "dump", "bottle", "codes",
		"stg", "stgcode", "q", "turb", "turbcode")
	df$chr <- make.chr(df$year, df$mo, df$dy, df$time)
	df[, c(13, 5:12)]
}