".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)
{
        # With TTS rev 5.0 no longer adds zero threshold
        # Returns a vector of thresholds 
        # Vector starts with tstart, and ends with tmax
	dx <- (f(tmax) - f(tstart))/(n - 1)
	x <- seq(from = f(tstart), by = dx, length = n)
	round(finv(x))
}


"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
}

"tts5.sample" <-
function(data, up, dn, stay = 2, revpct = c(10, 20), repwait = 18, minstg = 0, turblim = 2000, limskip = 2, initcode = 3, revval = 5, startwait = 72)
{
	#  Simulator to partially mimic sampling logic of Campbell TTS code
	#  Omits logic related to data dumps, equipment changes, and temperature conditions
	#  Rev  3.1  JL  061018 fixed bug that could cause a crash coming out of baseflow
	#       if j=0 because last falling turbidity was below all thresholds.
	#       Just had to move the statement that sets j when exiting baseflow
	#  Rev  3.0  JL  060927 incorporates changes corresponding to Campbell TTS rev 5.0
	#       Incorporates lastris and lastfal arrays
	#       Uses lastsam to remember latest sample of any kind
	#       No longer needs a zero or 9999 threshold
	#       "basewait" parameter renamed "startwait" to match Campbell code
	#       New defaults: repwait = 18, startwait = 72
	#       Added a column to ouput: the sampled threshold
	#       Columns are now: (1) interval number (2) threshold code (3) threshold
	#  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)
	thresh <- numeric(0)
        thrcode <- numeric(0)
        thrcount <- revcount <- 0
	lastris <- rep(-999,length(up))
	lastfal <- rep(-999,length(dn))
	lastsam <- -999
        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 <- 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[1]) && (i - startwait > lastsam)) {
                                # Startup sample
                                i <- i + 1
				lastsam <- i
                                sam <- c(sam, i)
                                thrcode <- c(thrcode, 4)
				thresh <- c(thresh, NA)
                        }
			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 <- 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(j <= length(up) & t >= up[j]) {
                                # Next rising threshold reached
                                thrcount <- thrcount + 1
                                if(thrcount >= stay) {
					# Threshold condition satisfied
                                        if(i > lastris[j] + repwait) {
                                        	# Take a sample
                                        	sam <- c(sam, i)
                                        	thrcode <- c(thrcode, tcode)
						thresh <- c(thresh,up[j])
                                        	lastris[j] <- i
						lastsam <- i
                                        }
                                        j <- 1 + sum(t >= up)
                                        thrcount <- 0
                                }
                        }
			else if (t >= turblim) {
				if (i > lastsam + limskip) {
					# Collect an 'overflow' sample
					sam <- c(sam, i)
					thrcode <- c(thrcode, 5)
					thresh <- c(thresh, NA)
					lastsam <- i
					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 <- sum(t > dn)
                                        tmin <- t
					if ((j < length(dn)) && (tmax > dn[j+1]) &&
                                            (i > lastfal[j+1] + repwait)) {
                                                # Take a sample
                                                sam <- c(sam, i)
                                                thrcode <- c(thrcode, tcode)
						thresh <- c(thresh, dn[j+1])
                                                lastfal[j+1] <- i
						lastsam <- 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(j > 0 && t <= dn[j]) {
                                # Next falling threshold reached 
                                thrcount <- thrcount + 1
                                if(thrcount >= stay) {
                                        # Threshold condition satisfied
                                        if(i > lastfal[j] + repwait) {
                                                # Take a sample
                                                sam <- c(sam, i)
                                                thrcode <- c(thrcode, tcode)
						thresh <- c(thresh,dn[j])
                                                lastfal[j] <- i
						lastsam <- i
                                        }
                                        j <- 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
                                        if ((j > 1) && (tmin < up[j-1]) && 
                                            (i > lastris[j-1] + repwait)) {
                                                # Take a sample
                                                sam <- c(sam, i)
                                                thrcode <- c(thrcode, tcode)
						thresh <- c(thresh,up[j-1])
                                                lastris[j-1] <- i
						lastsam <- i
                                        }
                                }
                        }
			else if (t >= turblim) {
				# Switch to rising
				tcode <- 1
				tmax <- t
				j <- 1 + sum(t >= up)
                                thrcount <- 0
				if (i > lastsam + limskip) {
					# Collect an 'overflow' sample
					sam <- c(sam, i)
					thrcode <- c(thrcode, 5)
					thresh <- c(thresh,NA)
					lastsam <- i
					revcount <- 0
				}
			}
                }
                thrcount <- 0
                revcount <- 0
        }
        	result <- data.frame(sam, thrcode, thresh)
		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 = "m/d/y", 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.
}

"ymd2date" <-
function(date, out = "m/d/y", 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)
}

"msm2mt" <-
function(msm)
{
	hr <- msm %/% 60.
	mn <- msm %% 60.
	round(100. * hr + mn)
}

"mt2msm" <-
function (time) 
{
  if(data.class(time)=="character")
    time<-as.numeric(time)
  min<-time%%100.
  hour<-time%/%100.
  60.*hour+min
}