fun.MEF_WXTmet2011a <- function(PLOT=T,GRAD=T,RAWPLOT=F,N=5,short=F)
	{

## Function designed to get 5 min. averages of the WXT510 met sensors
## at Manitou exp. forest.  Uses the  raw 1 sec data.
## Used to check the accuracy of my averaging program.
######## will then do some plotting of met profiles

### for data after Nov. 1, 2009

## PLOT - shows data plots
## GRAD == T -- we are computing profiles
## RAWPLOT == T, plots raw data, == F, plots filtered data
## FORMAT == 1	data format before 10/14/2008
## N == # of minutes you want to average (usually either 1, 5 or 10)
## Sensor = 1-4, which sensor you want to get averages for.
## short == T/F, if True - only what to calculate for a day's worth of data (for testing) -- 24 datafiles

### Input text file with names of input/output files in it.
### Will need to update "c:/data/Manitou08/Datafiles.txt" to analyze different months.
### will create monthly files (either 1 or 5 min. averages)

### for 2009
#datafile <- scan("c:/data/Manitou/Manitou09/Datafiles09.txt",what="character")
### For 2010
#datafile <- scan("c:/data/Manitou/Manitou010/Datafiles010.txt",what="character")
### For 2011
datafile <- scan("c:/data/Manitou/Manitou011/Datafiles011.txt",what="character")


## Parse "datafile" into the input and output filenames for the WXTs
infile <- datafile[31]		### put in correct number
MONTH <- datafile[30]		### put in correct number

### Loop through the 4 sensors.
Sensor <- 1
## read in the names of data files
datanames <- scan(infile,skip=1,sep=",",list(index=0,starttime="",filename1="",filename2="",filename3="",filename4=""))

while(Sensor <= 4)	{
	### Output data - will depend upon Sensor number:
	if(Sensor == 1)	{
		outfile <- datafile[32]
		WXTfiles <- datanames$filename1
		data.5minfile <- datafile[46]
		z <- 2
			}
	if(Sensor == 2)	{
		outfile <- datafile[33]
		WXTfiles <- datanames$filename2
		data.5minfile <- datafile[47]
		z <- 8
			}
	if(Sensor == 3) {
		outfile <- datafile[34]
		WXTfiles <- datanames$filename3
		data.5minfile <- datafile[48]
		z <- 15
			}
	if(Sensor == 4) {
		outfile <- datafile[35]
		WXTfiles <- datanames$filename4
		data.5minfile <- datafile[49]
		z <- 24
			}

t.starts <- datanames$starttime
t.startA <- strptime(t.starts, "%m/%d/%Y %H:%M")
t.starts <- as.POSIXct(t.startA)
index <- datanames$index
cat("\n Sensor #", Sensor)

cat("\n Correct Input files?")
#browser()

#### Get the 5 min avg. files first!! - then we will average the 1 sec data and substitute in.
## read in the names of data files
data.5min <- matrix(scan(data.5minfile,skip=1,sep=",",list(datetime="",wdir=0,wspd=0,T=0,RH=0,P=0,R1=0,R2=0,H1=0,sUdir=0,sU=0,sT=0,sRH=0,sP=0,z=0)),byrow=T,ncol=15)

date5m <- as.ts(data.5min[[1]])
date5m <- strptime(date5m, "%m/%d/%Y %H:%M")
date5m <- fun.Ltimetag(date5m)
#### now have numeric dates - matrix of yy,jday,hrx100,min,sec
## get decimal day:
dday5m <- date5m[,2] + (date5m[,3]/2400) + (date5m[,4]/1440)
dday5m <- signif(dday5m,d=7)
date5m <- cbind(date5m,dday5m)		## add decimal day to date.
#####################################################################################
### other 5 min. data:
wdir.5m <- as.ts(data.5min[[2]])
wspd.5m <- as.ts(data.5min[[3]])
Temp.5m <- as.ts(data.5min[[4]])
RelH.5m <- as.ts(data.5min[[5]])
Pres.5m <- as.ts(data.5min[[6]])
PPT1.5m <- as.ts(data.5min[[7]])
PPT2.5m <- as.ts(data.5min[[8]])
Hai1.5m <- as.ts(data.5min[[9]])
swdir.5m <- as.ts(data.5min[[10]])
swspd.5m <- as.ts(data.5min[[11]])
sT.5m <- as.ts(data.5min[[12]])
sRH.5m <- as.ts(data.5min[[13]])
sP.5m <- as.ts(data.5min[[14]])
z.5m <- as.ts(data.5min[[15]])
Flag.5m <- rep(5,length=length(z.5m))
######### Put back together as matrix//same format as the final output matrix from teh 1 sec data.// put in wind spd for vector wind spd!
DATA.5min <- cbind(date5m,wdir.5m, wspd.5m,Temp.5m,RelH.5m,Pres.5m,PPT1.5m,PPT2.5m,Hai1.5m,swdir.5m,swspd.5m,sT.5m,sRH.5m,sP.5m,wspd.5m,z.5m,Flag.5m)
cat("\n Correct 5 min avg'd data?")
#browser()


#########################################################################################
##### now start working with the 1 sec data.
i <- 1
if(short == F)	Lfile <- length(WXTfiles)
else		Lfile <- i + 24			## use only 25 datafiles
cat("\n Month=",MONTH," Sensor=",Sensor)
k <- 1		## index for Labview 5 min. avg'd data

#### Main Loop begin:
while(i <= Lfile) {

cat("\n Index:",index[i],"   Loop #:",i)
cat("\n\n")

cat("\n i=",i,"k=",k)
filename.i <- as.numeric(WXTfiles[i])
p <- is.na(filename.i)

#if(Sensor == 2 & i == 575)	browser()

#if(i == 33)	browser()
#if(i == 4)	browser()

if(p != 0)	{

d1 <- matrix(scan(WXTfiles[i],sep="\t",list(date="",msec=0,WD=0,WS=0,T=0,RH=0,P=0,Rain=0,RR=0,Hail=0)),byrow=T,ncol=10)

### parse the data
date1 <- as.ts(d1[[1]])
WD1 <- as.ts(d1[[3]])
WS1 <- as.ts(d1[[4]])
T1 <- as.ts(d1[[5]])
RH1 <- as.ts(d1[[6]])
P1 <- as.ts(d1[[7]])
Rain1 <- as.ts(d1[[8]])
RR1 <- as.ts(d1[[9]])
Hail1 <- as.ts(d1[[10]])
DATA1 <- cbind(WD1,WS1,T1,RH1,P1,Rain1,RR1,Hail1)

date1 <- strptime(date1, "%m/%d/%Y %H:%M:%S")
date1 <- fun.Ltimetag(date1)
#### now have numeric dates - matrix of yy,jday,hrx100,min,sec
## get decimal day:
dday <- date1[,2] + (date1[,3]/2400) + (date1[,4]/1440)
dday <- signif(dday,d=7)
date1 <- cbind(date1,dday)		## add decimal day to date.
##############################################################################

## Loop through data based on "N" (number of minutes/average)
###### Do we want to initialize this to a multiple of 5 mins.????
j <- 1		## index for 1sec calc. data

Flag.1m <- 1	## will let us know what data we have substituted.

#if(i == 401)	browser()

Ldate <- length(date1[,1]) - 2
while(j <= Ldate)	{
	
	
	
	## initially - need to get to the nearest 01 seconds and multiple of 5min.:
	InitMin <- date1[j,4]/5
	Intmin <- as.integer(InitMin)

	#if(Sensor == 2 & k == 5572)	browser()

	if(InitMin == Intmin)	{
		Nsamples <- N*60
		end <- j + Nsamples
		### for screwed up files:
		if(end > Ldate)		end <- Ldate
		if(j >= Ldate)		j <- Ldate - 20
		########################################
		ttag.i <- date1[end,]		## time tag at END of period!
		## get subarray of data:
		
		#cat("\n j=",j,"   end=",end)
		M <- DATA1[j:end,]
		## get simple averages, summs and std deviations:
		Mavg <- apply(M,2,mean, na.rm=T)
		Msum <- apply(M,2,sum, na.rm=T)
		Mstd <- (apply(M, 2,var, na.rm=T))^.5
		
		### now need to get vector average for wind direction and speed.
		Windstuff <- fun.windaverage(M[,1],M[,2])		## uses wind dir and spd.
		MeanWD <- signif(Windstuff$avgWD,d=5)
		StdWD <- signif(Windstuff$stdWD,d=4)
		Vws <- signif(Windstuff$Vws,d=4)		## vector-avg'd wind speed
		avgWS <- signif(Windstuff$U,d=3)
		stdWS <- signif(Windstuff$stdU,d=3)

		############ In Aug. 2010 - realized that WXT #4 was off by 30o
		############ relative to the sonic and cup/vane
		if(Sensor == 4)	{
			MeanWD <- MeanWD - 30
			MeanWD <- ifelse(MeanWD < 0, (MeanWD + 360),MeanWD)
				}		

		## now combine the averages you want and time tag.
		## columns: ttag.i, 
		Mout <- c(ttag.i,MeanWD,avgWS,signif(Mavg[3:5],d=4),Mavg[6:8],StdWD,stdWS,signif(Mstd[3:5],d=2),Vws,z,Flag.1m)
		#cat(ttag.i)
		#cat("\n i=",i,"j=",j)
		#if(Sensor == 1 & i == 332 & k > 3974)	browser()
##########################################################################################
		### Now compare time tags between original 5 min. data and what we have just computed.
		if(k > length(DATA.5min[,1])) 	{
					cat("\n k greater than 5 min file - #1 spot")
					k <- length(DATA.5min[,1])
					#browser()
						}

		## if you want to see the time tags on the screeen:
		#cat("\n Time tags: 1 sec:",ttag.i[6],"5 min.data:",DATA.5min[k,6])
		d.ttag <- DATA.5min[k,6] - ttag.i[6]

		#if(k == 481)		browser()

		#if(j < 100)	browser()

		#if(j > 1000 & j < 1400)	browser()

		## if they are the same - then put computed data in the row.	
		if(abs(d.ttag) < 0.0017)	DATA.5min[k,] <- Mout
		else 	{
			### 2 cases - common case - d.ttag will be negative- so 1 sec data is later than 5 min.
			### Case 1:
			if(d.ttag < 0)	{
				### Special Case - k > the length of the 5 min. file
				if(k >= length(DATA.5min[,1])) 	{
					cat("\n k greater than 5 min file")
					## just append Mout to the 5 min data file
					DATA.5min <- rbind(DATA.5min,Mout)
						}
				else		{
			### need to increment k until they match times again.
				cat("\n Time tag diff: k=",k,dday5m[k])
				#if(i == 401)	browser()
				while(abs(d.ttag) > 0.0015)	{
					k <- k + 1
					d.ttag <- DATA.5min[k,6] - ttag.i[6]
						#cat("\n delta ttag", d.ttag)
						### if ttag "jumps" over the matching time.
						if(d.ttag > 0)	 d.ttag <- 0
						}		## end of while loop
					
					#browser()
					DATA.5min[k,] <- Mout		## once matched again - put computed data in row.
						}			## end of else statement
							}		## end of if(d.ttag < 0.0015)
			### Case 2: 1 sec data is before the 5 min data (mostly at beginning of file!)
			if(d.ttag > 0.0015)	{
				if((k-1) == 0)	{
				### special case - at the beginning -  just add the row at the beginning		
					cat("\n Time tag diff - file beginning: k=", k)
					#browser()
					DATA.5min <- rbind(Mout,DATA.5min)
					
							}
				else		{
				### put this row between rows(k - 1) and k
					cat("\n Time tag diff - Insert row: k=", k)
					#browser()
					LK <- length(DATA.5min[,1])
					DATA.5min <- rbind(DATA.5min[1:(k-1),],Mout,DATA.5min[k:LK,])
						}	
					}	## end of if(d.ttag > 0.0015)

			}		## end of else
		## go to the next sample.
		j <- j + Nsamples
		k <- k + 1
		
		#dput(Nout,outfile)

		}		### end of if(Initmin=Intmin)
	else		{
		## just go around the loop again
		j <- j + 1
		}

		}		## end of sub-while loop (within a file)
		}		## end of if(filename.i !=0)

		### go to the next file
		i <- i + 1

		}		### next 1-hour file (end of main while loop)

### at end of while loop
### output file in ASCI format
fun.asci(DATA.5min,outfile)

if(Sensor == 1)		Data.S1 <- DATA.5min
if(Sensor == 2)		Data.S2 <- DATA.5min
if(Sensor == 3)		Data.S3 <- DATA.5min
if(Sensor == 4)		Data.S4 <- DATA.5min

#################################################
### some plots ##################################

## some parsing for easy editing later
DOY <- DATA.5min[,6]
WDir <- DATA.5min[,7]
WSpd <- DATA.5min[,8]
Temp <- DATA.5min[,9]
RHum <- DATA.5min[,10]
Press <- DATA.5min[,11]
ppt <- DATA.5min[,12]
## 4 plots/page
plot.new()
par(mfrow=c(4,1))

if(Sensor == 1)		{
	plot(DOY,Temp,type="l",ylab="T oC",xlab="DOY")
	plot(DOY,RHum,type="l",ylab="RH %",xlab="DOY")
	plot(DOY,WSpd,type="l",ylab="W Spd, m/s",xlab="DOY")
	plot(DOY,Press,type="l",ylab="Pressure,mbar",xlab="DOY")
			}

if(Sensor == 2)	{
	plot(DOY,Temp,type="l",ylab="T oC",xlab="DOY",col=2)
	plot(DOY,RHum,type="l",ylab="RH %",xlab="DOY",col=2)
	plot(DOY,WSpd,type="l",ylab="W Spd, m/s",xlab="DOY",col=2)
	plot(DOY,Press,type="l",ylab="Pressure,mbar",xlab="DOY",col=2)
			}	
if(Sensor == 3)	{
	plot(DOY,Temp,type="l",ylab="T oC",xlab="DOY",col=3)
	plot(DOY,RHum,type="l",ylab="RH %",xlab="DOY",col=3)
	plot(DOY,WSpd,type="l",ylab="W Spd, m/s",xlab="DOY",col=3)
	plot(DOY,Press,type="l",ylab="Pressure,mbar",xlab="DOY",col=3)
			}	
if(Sensor == 4)	{
	plot(DOY,Temp,type="l",ylab="T oC",xlab="DOY",col=4)
	plot(DOY,RHum,type="l",ylab="RH %",xlab="DOY",col=4)
	plot(DOY,WSpd,type="l",ylab="W Spd, m/s",xlab="DOY",col=4)
	plot(DOY,Press,type="l",ylab="Pressure,mbar",xlab="DOY",col=4)
			}	


#plot.new()
#par(mfrow=c(2,1))
#plot(DOY,WDir,type="l",ylab="Wind Direction",xlab="DOY")
#plot(DOY,ppt,type="l",ylab="Precip",xlab="DOY")
cat("\n end of plots")

Sensor <- Sensor + 1
	}	### end of sensor while loop - move on to the next sensor!


###### Now plot some gradients from all of the sensors:
## some parsing for easy editing later
DOY1 <- Data.S1[,6]
WDir1 <- Data.S1[,7]
WSpd1 <- Data.S1[,8]
Temp1 <- Data.S1[,9]
RHum1 <- Data.S1[,10]
Press1 <- Data.S1[,11]
pp1t <- Data.S1[,12]

DOY2 <- Data.S2[,6]
WDir2 <- Data.S2[,7]
WSpd2 <- Data.S2[,8]
Temp2 <- Data.S2[,9]
RHum2 <- Data.S2[,10]
Press2 <- Data.S2[,11]
pp2t <- Data.S2[,12]

DOY3 <- Data.S3[,6]
WDir3 <- Data.S3[,7]
WSpd3 <- Data.S3[,8]
Temp3 <- Data.S3[,9]
RHum3 <- Data.S3[,10]
Press3 <- Data.S3[,11]
pp3t <- Data.S3[,12]

DOY4 <- Data.S4[,6]
WDir4 <- Data.S4[,7]
WSpd4 <- Data.S4[,8]
Temp4 <- Data.S4[,9]
RHum4 <- Data.S4[,10]
Press4 <- Data.S4[,11]
pp4t <- Data.S4[,12]

cat("\n Ready to do some final plots")
browser()

## 4 plots/page
plot.new()
par(mfrow=c(4,1))

plot(DOY1,Temp1,type="l",ylab="T oC",xlab="DOY")
lines(DOY2,Temp2,type="l",col=2)
lines(DOY3,Temp3,type="l",col=3)
lines(DOY4,Temp4,type="l",col=4)

plot(DOY1,RHum1,type="l",ylab="RH %",xlab="DOY",ylim=c(0,100))
lines(DOY2,RHum2,type="l",col=2)
lines(DOY3,RHum3,type="l",col=3)
lines(DOY4,RHum4,type="l",col=4)

plot(DOY4,WSpd4,col=4,type="l",ylab="W Spd, m/s",xlab="DOY",ylim=c(0,12))
lines(DOY2,WSpd2,type="l",col=2)
lines(DOY3,WSpd3,type="l",col=3)
lines(DOY1,WSpd1,type="l",col=1)

plot(DOY1,Press1,type="l",ylab="Pressure,mbar",xlab="DOY",ylim=c(730,800))
lines(DOY2,Press2,type="l",col=2)
lines(DOY3,Press3,type="l",col=3)
lines(DOY4,Press4,type="l",col=4)



browser()

##### now use M3 to get some diurnal averages.
if(1==0)	{
####### BUT WE HAVE 5 MIN. DATA - NOT 30 MIN. DATA !!!!!!! ######################
### make an hrmin vector
hrmin <- Nout[,3]/100 + Nout[,4]/60		## get hrmin (0.0, 0.5, 1.0,... each 1/2 hr)
hrmin <- round(hrmin,d=1)
hrmin <- ifelse(hrmin > 23.75,0.0,hrmin)
M3 <- cbind(dday,hrmin,WSpd,WDir,Temp,RHum,Press,ppt)
WIND <- T
DielAvg <- fun.diurnalavg(M3,WIND)
par(mfrow=c(2,1))
### plot diurnal avgs.
plot(DielAvg[,1],DielAvg[,8],type="l",ylab="Wind Speed")
plot(DielAvg[,1],DielAvg[,9],type="l",ylab="Wind Direction")
browser()
plot(DielAvg[,1],DielAvg[,10],type="l",ylab="Temp")
plot(DielAvg[,1],DielAvg[,11],type="l",ylab="Humidity")
browser()
plot(DielAvg[,1],DielAvg[,12],type="l",ylab="Pressure")
plot(DielAvg[,1],DielAvg[,13],type="l",ylab="Precip")
cat("\n End of plots - do we want to output diurnal data??")
browser()
		}		## end of if(1==0)
######## Do we want to output this !!!!
################################


### end of function #############
}


########################################################################################
##############3 subroutines ########################

fun.despike <- function(x, name, n)
	{

# Subroutine which takes a time series, x, with a given "name" 
# and filters out any obvious spikes in the data. 
# Can be used for sonic data(name = "u","v",etc), CO2 data (name =
# "co2","h2o","t.ec",etc), Kr Hygrometer data (name = "h2o.kr"), etc.
# The filter used is a essentially n*standard deviation.  Anything
# values that fall outside of n*sigma will be replaced with the
# median value.

# Define the filter - Our cut-off :
		y <- n*(var(x,na.rm=T)^0.5)

# Look at positive spikes first :
	test.1 <- x > (median(x,na.rm=T) + y)
		z2 <- sum(test.1,na.rm=T)
	if(z2 > 0)	{
		cat(" :\n replacing ", z2)
		cat(" (+) spikes in ",name)
		x <- ifelse(x > median(x,na.rm=T)+y, median(x,na.rm=T), x)	}

# Look at negative spikes :
	test.1 <- x < (median(x,na.rm=T) - y)
		z3 <- sum(test.1,na.rm=T)
	if(z3 > 0)	{
		cat(" :\n replacing ", z3)
		cat(" (-) spikes in ",name)
		x <- ifelse(x < median(x,na.rm=T)-y, median(x,na.rm=T), x)	}

# Generate a flag which tells when greater than 5% of the time series
# has spikes in it.
	per.spike <- (z2 + z3)/ length(x)
	if(per.spike>0.04)	{
		SPIKE <- F	}
	else	{
		SPIKE <- T	}

# return list of filtered x and SPIKE flag
list(x=x, SPIKE=SPIKE)

}

########################################################################

fun.fastplot <- function(t,x, Label)
{

# subroutine to plot profiles - can change xax and yax

## t = time tag matrix - has six columns with six time tags
## x = matrix with 6 levels (columns 1-6, levels 1-6)
## Label = characters for ylabel

p <- is.na(x[,1])
ind <- seq(from=1,to=length(p),by=1)
samples <- ind[p==0]
I.start <- min(samples,na.rm=T)
I.end <- max(samples,na.rm=T)
xax <- c(t[I.start,1],t[I.end,1])
y.min <- min(x,na.rm=T)
y.max <- max(x,na.rm=T)
yax <- c(y.min,y.max)

#browser()

plot(t[,1],x[,1],type="l",xax,ylab=Label,ylim=yax)
lines(t[,2],x[,2],col=2)
lines(t[,3],x[,3],col=3)
lines(t[,4],x[,4],col=4)
lines(t[,5],x[,5],col=5)
lines(t[,6],x[,6],col=6)

#### create legend
leg.names <- c("z=2m","z=10.9m","z=16.7m","z=23.9m","z=30.3m","z=39.7m")
#colors <- c("black","red","green","blue","cyano","magenta")
#legend(locator(1),leg.names,col=c(1,2,3,4,5,6),fill=T)
#legend(5,1,leg.names,col=c(1,2,3,4,5,6),fill=T)
text(locator(1),"black=2m, red=10.9m, grn=16.7m, blue=23.9m, cyano=30.3m, mag.=39.7m")


}

##################################################

fun.asci <- function(x,data.file)
{

## function to output an asci file

y <- matrix(NA,nrow=length(x[1,]),ncol=length(x[,1]))
i 	<- 1
while(i <= length(x[1,])) 
{	
	y[i,] 	<- x[,i]
	i 	<- i+1
}

### Output data:
write(y,data.file,ncol=length(x[1,]))


}



######################################################################
fun.Camptimetag <- function(day,hrmin,sec)
#
# Converts Campbell-type time tags into decimal day time tags:
# note - we assume that we know the year!!!!
{

## get hours and minutes separated
hr <- floor(hrmin/100)			## give hour of day
min <- (hrmin/100) - floor(hrmin/100)	## min/100
min <- 100*min

## put everything in days and add it up:
min.sec <- min + sec/60
hr.min <- hr + min.sec/60
ttag <- day + hr.min/24		
ttag <- signif(ttag,d=8)
# output ttag
	ttag

}


###############################################################

	
fun.linearextrap <- function(index,x,xx,yy)
	{
## subroutine to take a time series (x) and and index (0->length)
## and a set of index pts/calibration points (xx and yy) and create a 
## time series of calibration constants.
### does linear extrapolation between different xx points.
### xx == vector of indices
### yy == corresponding signals (for zero's) or gain coeff. (for spans)

cat("\nTest linear extrap")
#browser()

### do line fit between each pt.
Lz <- length(xx) - 1
i <- 1

while(i <= Lz)	{
	
j <- i + 1
x1 <- xx[i:j]
y1 <- yy[i:j]
#browser()
Lfit <- lm.fit <- lm(y1 ~ x1)
Int <- Lfit$coeff[1]
m <- Lfit$coeff[2]
x[xx[i]:xx[j]] <- Int + (m*index[xx[i]:xx[j]])
i <- i + 1

}		###end of loop

#### then use last zero points to fill out the zero time series
Lz <- length(yy)
Last.zero <- yy[Lz]
Li <- length(index)
Lend <- xx[Lz]
x[Lend:Li] <- yy[Lz]



x


}

#####################################################################
fun.skewness <- function(x,OUT=F)
	{
## function designed to take a time series of x and 
## calculate the mean, median, std. dev., skewness and
## kurtosis.

# get rid of NAs for the rest :
xx <- is.na(x)
if(sum(xx,na.rm=T) < length(xx))	{

N <- length(x) - sum(xx,na.rm=T)
xbar <- mean(x,na.rm=T)
xmed <- median(x,na.rm=T)
x.min <- min(x,na.rm=T)
x.max <- max(x,na.rm=T)

sdx <- var(x,na.rm=T)^.5

Skx <- mean((x - mean(x))^3) / (sdx^3)

Krx <- mean((x - mean(x))^4) / (sdx^4)

if(OUT==T)	{
	cat("\n N = ",N)
	cat("\n xbar = ",xbar)
	cat("\n xmed = ",xmed)
	cat("\n sdx = ",sdx)
	cat("\n x.min = ",x.min)
	cat("\n x.max = ",x.max)
		}
	}	## end of if(sum < length)
else	{
	N <- 0
	xbar <- NA
	xmed <- NA
	sdx <- NA
	Skx <- NA
	Krx <- NA
		}
	
list(N=N,xbar=xbar,xmed=xmed,sdx=sdx,Skx=Skx,Krx=Krx,x.min=x.min,x.max=x.max)

#invisible()

# end function
}
##############################################################
######################################################



fun.Ltimetag <- function(t1)
	{
# function to convert a POSIXct object to a numeric string of the type:
#	Julian day,  hour of day, min of day,  sec
# output as a vector of numbers:

year <- strftime(t1, "%Y")
year <- (as.numeric(year))
jday <- strftime(t1,"%j")
jday <- as.numeric(jday)
hday <- strftime(t1,"%H")
hday <- (as.numeric(hday))*100
mday <- strftime(t1, "%M")
mday <- as.numeric(mday)
Sday <- strftime(t1, "%S")
Sday <- as.numeric(Sday)

## put in to vector
ttag <- cbind(year,jday,hday,mday,Sday)

}
 
##################################################

########################################################################

fun.windaverage <- function(WD,WS)
	{
## inputs - vectors of wind direction and wind speed that we 
## want to average.

## calculate wind speed and direction 
# (1) calculate scalar wind speed
	U <- mean(WS,na.rm=T)
	stdU <- var(WS, na.rm=T)^0.5

### now 0-360 scale
### Jan. 2001 - this is the correct way to get average wind direction.
	wd <- WD*(pi/180)		# get into radians
	Ux <- (mean(sin(wd),na.rm=T))
	Uy <- (mean(cos(wd),na.rm=T))
	wdir <- (atan(Ux/Uy))*(180/pi)
	## Now get the proper quadrant :
		if(Uy < 0)	wdir <- wdir + 180
		if(Uy > 0)   wdir <- ifelse(wdir < 0, 360+wdir, wdir)
	avgWD <- wdir

	## Vector Wind speed :
	Ue <- mean((WS*sin(wd)),na.rm=T)
	Un <- mean((WS*cos(wd)),na.rm=T)
	Vws <- (Ue^2 + Un^2)^0.5

	# New way of getting std. dev. of wind direction
	## From Campbell 23x book - Yamartino algorithm (pp11-2)
	epsil <- ( 1 - (Ux^2 + Uy^2))^.5
	sdUdir <- asin(epsil)*(1 + (0.1547*epsil^3))
	stdWD <- sdUdir*(180/pi)

list(U=U,avgWD=avgWD,Vws=Vws,stdU=stdU,stdWD=stdWD)
}

###############################################################


##############################################################
fun.diurnalavg <- function(X,WIND)
	{
## Subroutine to get diurnal avg from time series of profiles for a single species
#  X = time series with columns:  doy, hrmin, + however many columns you want to average.

## if(WIND == T)	then 4th column of X is wind direction and must be vector-avg'd

Lcol <- length(X[1,])		## number of input columns
NumCol <- (Lcol*5) + 1		## number of output columns	
i <- 0.0
j <- 1
Y <- matrix(NA,ncol=(Lcol*5),nrow=48)
hrmin <- X[,2]
#browser()

while(i <= 23.5)	{
	#if(i == 7.0)	browser()
	test <- hrmin == i
	### should index all the points at 0.0 hours
	Xtest <- X[test,]
	### get means, medians and std. dev.and # of points (minus the NAs)
	P <- is.na(Xtest[,3:Lcol])
	Lna <- Lcol - 2
	P[,1:Lna] <- ifelse(P[,1:Lna] == 1,0,1)		### switch around na matrix.
	N <- apply(P,2,sum,na.rm=T)		### number of pts.
	xbar <- apply(Xtest[,3:Lcol], 2,mean,na.rm=T)
	xmed <- apply(Xtest[,3:Lcol],2,median,na.rm=T)
	xvar <- apply(Xtest[,3:Lcol],2,var,na.rm=T)
	xstd <- (xvar)^0.5
	xsterr <- xstd/(N^.5)

	if(WIND == T)	{
		WD <- Xtest[,4]			### should be 4th column
		wd <- WD*(pi/180)		# get into radians
		Ux <- (mean(sin(wd),na.rm=T))
		Uy <- (mean(cos(wd),na.rm=T))
		wdir <- (atan(Ux/Uy))*(180/pi)
	## Now get the proper quadrant :
		if(Uy < 0)	wdir <- wdir + 180
		if(Uy > 0)   	wdir <- ifelse(wdir < 0, 360+wdir, wdir)
	Udir <- wdir		## UDIR is avg. wind direction
	# New way of getting std. dev. of wind direction
	## From Campbell 23x book - Yamartino algorithm (pp11-2)
	epsil <- ( 1 - (Ux^2 + Uy^2))^.5
	sdUdir <- asin(epsil)*(1 + (0.1547*epsil^3))
	sdUdir <- sdUdir*(180/pi)
	sderrUdir <- sdUdir/((N[2])^.5)
			}
	Lxbar <- length(xbar)

	### put into a row of output matrix 
	if(WIND == F)	 Y[j,] <- c(i,N,xbar,xmed,xstd,xsterr)
	else		Y[j,] <- c(i,N,xbar[1],Udir,xbar[3:Lxbar],xmed,xstd[1],sdUdir,xstd[3:Lxbar],xsterr[1],sderrUdir,xsterr[3:Lxbar])
	###########################
	## now move to the next 1/2 hour
	i <- i + 0.5
	j <- j + 1
	}		## end of the while loop

## output Y
Y

}

# end function
##########################################################################