fun.MEFcalibrate011 <- function(o3.zero=T,so2.zero=T,NEW=T)
	{

## Function designed to go through trace gas data from Manitou Forest and :
## (1) find calibrations - either manual or automatic
## (2) Calculate zero's and gains (mainly gain coefficients) for several trace gases.
## (3) 
######## These coefficients will be used in MEFtracegas.R to calculate concentrations.

###### does linear extrapolation between points!!
##  NEW ==T, new cal. scheme with zero offsets
##	==F, old cal. scheme.
####################################################
### Works:  4/28/11
### 
### then updated my calibration scheme and need to use fun.MEFcalibrate010b.R!!!
### only zero the NO analyzer with the pre-rxn chamber.
## 
### Input text file with names of input/output files in it.
### Will need to update "c:/data/Manitou/Manitou09/Datafiles09.txt" to analyze different months.
### will create monthly files (either 1 or 5 min. averages)
#datafile <- scan("c:/data/Manitou/Manitou09/Datafiles09.txt",what="character")
datafile <- scan("c:/data/Manitou/Manitou011/Datafiles011.txt",what="character")

## Parse "datafile" into the input and output filenames
infile <- datafile[53]
### Output data
calfile <- datafile[21]
calflags <- datafile[20]
outflags <- datafile[19]
MONTH <- datafile[30]

#### REAd in data #######
#### works for Oct. 2009 data to present #############
cat("\n Month:",MONTH)
d <- matrix(scan(infile,sep=",",fill=T,skip=1,list(yr=0,jday=0,hr=0,mn=0,sec=0,doy=0,Vso2=0,Vo3=0,Vno=0,Vco2=0,Vh2o=0,Zero=0,Span=0,NOstate=0,Level=0,Vfstd=0,Vftot=0,sdVso2=0,sdVo3=0,sdVno=0,sdVco2=0,sdVh2o=0,sdNOs=0,sdL=0,Lflag=0,Flag1s=0)),byrow=T,ncol=26)
## parse data:
########## Time tag
yr <- as.ts(d[[1]])
jday <- as.ts(d[[2]])
hr <- as.ts(d[[3]])
mn <- as.ts(d[[4]])
sc <- as.ts(d[[5]])
dday <- as.ts(d[[6]])
#### Raw voltages
Vso2 <- as.ts(d[[7]])
Vo3 <- as.ts(d[[8]])
Vno <- as.ts(d[[9]])
Vco2 <- as.ts(d[[10]])
Vh2o <- as.ts(d[[11]])
### STatus ########
Zero <- as.ts(d[[12]])
Span <- as.ts(d[[13]])
NOstate <- as.ts(d[[14]])
Level <- as.ts(d[[15]])
### Flows (voltages) ###3
VFstd <- as.ts(d[[16]])
VFtot <- as.ts(d[[17]])
#### Standard deviations ####
sdSO2 <- as.ts(d[[18]])
sdO3 <- as.ts(d[[19]])
sdNO <- as.ts(d[[20]])
sdCO2 <- as.ts(d[[21]])
sdH2O <- as.ts(d[[22]])
###########################################
d <- matrix(scan(calflags,sep=",",fill=T,skip=1,list(date=0,Zso2=0,Sso2=0,Zo3=0,Zno=0,Sno=0,Zco2=0,Sco2=0,Zh2o=0,Cso2=0,Cno=0,Cco2=0)),byrow=T,ncol=12)
## parse data:
#### Calibration Flags
Fdate <- as.ts(d[[1]])
F.Zso2 <- as.ts(d[[2]])
F.Sso2 <- as.ts(d[[3]])
F.Zo3 <- as.ts(d[[4]])
F.Zno <- as.ts(d[[5]])
F.Sno <- as.ts(d[[6]])
F.Zco2 <- as.ts(d[[7]])
F.Sco2 <- as.ts(d[[8]])
F.Zh2o <- as.ts(d[[9]])
SO2.cal <- as.ts(d[[10]])
NO.cal <- as.ts(d[[11]])
CO2.cal <- as.ts(d[[12]])
F.Zno2 <- F.Zno

### could be a few values of "60" in minutes. replace with "0"
mn <- ifelse(mn > 59.5,0,mn)

### get time tags
#date1 <- strptime(date1, "%m/%d/%Y %H:%M")
#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(yr,jday,hr,mn,sc,dday)		## matrix with 6 columns.
## also create an index - sequential numbers to reference
Ldata <- length(dday)
index <- seq(from=1,to=Ldata,by=1)
cat("\n Check on time tags!")
#browser()
############################################

### quick data filter for flags
p <- is.na(Zero)
Zero <- ifelse(p == T,0,Zero)
p <- is.na(Span)
Span <- ifelse(p==T,0,Span)

par(mfrow=c(5,1))
## Plot the raw voltages one at a time.
plot(index,Vso2,type="l",ylab="SO2")
#browser()
plot(index,Vo3,type="l",ylab="O3")
#browser()
plot(index,Vno,type="l",ylab="NO")
#browser()
plot(index,Vco2,type="l",ylab="CO2")
#browser()
plot(index,Vh2o,type="l",ylab="H2O")
#browser()
cat("\n End of voltage plots")
browser()

### Now go through and find zero's and span's and set some T/F indices #########################
## For Jan.7-April 13, 2011, used the old scheme.
## otherwise the new scheme.

if(NEW==T)	{
## get indices for when a zero occurs (will be same for all species)
Z1 <- Zero > 0.3 & Zero < 1.4			### T/F vector of where the zero's occur, T== zero occurring
		### as of 9/9/10:  this is the zero air method!!!!  All species!!!
Z2 <- Zero > 1.4
		### as of 9/9/10:  this is the pre-rxn chamber method for NO/NOx
## get indices for when a span occurs (will be same for all species)
S1 <- Span > 1.3			### T/F vector of where the span's occur, T== span occurring
	### as of 9/9/10:  this is only the Span for CO2
S2 <- Span > 0.4 & Span < 1.4			### T/F vector of where the span's occur, T== span occurring
S <- S2 == T & Level == 3
	### as of 9/9/10:  this is now the span for NO and SO2.
		}
else		{
## get indices for when a zero occurs (will be same for all species)
Z1 <- Zero > 0.3			### T/F vector of where the zero's occur, T== zero occurring
Z2 <- Z1
## get indices for when a span occurs (will be same for all species)
S <- Span > 1.3			### T/F vector of where the span's occur, T== span occurring
S1 <- S
		}

########## SPECIAL CASES:  Jan. and April ##############3
##### both of these months - the cal. scheme was split somewhere in the middle.
##### treat them separately.
if(MONTH=="jan")	{
	# run with NEW==T
	# break into 2 sections - new code from Jan. 1-Jan. 5., old code after that
	# we know that this break occurred around index==6000.
	## Z1 should be OK!!!  and we really don't use Z2 in the loop later!! (will just compute Z2 incorrectly and NO will not be spanned!)
	##just have to deal with CO2 span differently.  S1 should give correct times for either method (Span=2), but
	## this will not be correct for NO/SO2 spans with old method.
	## S1 is OK for CO2 (always), need to remake "S" for NO and SO2
	 Stest1 <- S					## Level 3 and Span==1
	 Stest2 <- Span > 1.3 & Level == 3		## Level 3 and Span==2
 	 newS <- ifelse(Stest1==T,Stest1,Stest2)
	 S <- newS
		}

if(MONTH=="apr")	{
	# run with NEW==T
	# break into 2 sections - old code from April. 1-Apr. 15., new code after that
	# 
	## Z1 should be OK!!!  and we really don't use Z2 in the loop later!! (will just compute Z2 incorrectly and NO will not be spanned!)
	##just have to deal with CO2 span differently.  S1 should give correct times for either method (Span=2), but
	## this will not be correct for NO/SO2 spans with old method.
	## S1 is OK for CO2 (always), need to remake "S" for NO and SO2
	 Stest1 <- S					## Level 3 and Span==1
	 Stest2 <- Span > 1.3 & Level == 3		## Level 3 and Span==2
 	 newS <- ifelse(Stest1==T,Stest1,Stest2)
	 S <- newS
		}

### also some monthly-based filters for calibration data
	Fil.o3 <- 100
	Fil.co2 <- 150
	Fil.h2o <- 500
	Fil.so2 <- 600
	FilB.so2 <- -250
if(MONTH=="jan")	{
	Fil.o3 <- 45
	Fil.co2 <- 115
	Fil.h2o <- 275
	Fil.so2 <- 400
	FilB.so2 <- 0
			}
if(MONTH=="feb")	{
	Fil.o3 <- 40
	Fil.co2 <- 100
	Fil.h2o <- 250
	Fil.so2 <- 70
	FilB.so2 <- 0
			}


######################################################################################
######## Get some initial values/time series for calibration coeff's. ################
### if no auto zeros - get some default values
if(o3.zero == T)	Z.ozone <- 23		## median of aug 2009
else			{
			Z.ozone <- 25		## Cal. on 4/16/08, zero = 21.88 mV
			##3 then make all O3 zero flags == 0
			F.Zo3 <- rep(0,length=Ldata)
			cat("\n O3 zero =",Z.ozone)
			}

if(so2.zero == T)	Z.sulf <- 0		## median of aug 2009
else			{
			testZ <- sort(Vso2, na.rm=T)
			### use 500 lowest value pts to get the zero.
			Z.sulf <- mean(testZ[1:500],na.rm=T)
			### then make all SO2 zero flags == 0
			F.Zso2 <- rep(0,length=Ldata)
			cat("\n SO2 zero =",Z.sulf)			
			}

Z.o3 <- rep(Z.ozone,Ldata)
Z.so2 <- rep(Z.sulf,Ldata)
Z.no <- rep(446,Ldata)
Z.nox <- rep(446,Ldata)
Z.noy <- rep(446,Ldata)
Z.co2 <- rep(118,Ldata)
Z.h2o <- rep(342,Ldata)
Z.Fstd <- rep(0,Ldata)
Z.Ftot <- rep(0,Ldata)

Z2.no <- rep(446,Ldata)
Z2.nox <- rep(446,Ldata)
Z2.noy <- rep(446,Ldata)


S.co2 <- rep(3.125,Ldata)
#### Some start values: #########################################################
##### NOTE - THESE WILL BE THE VALUES USED IF NO AUTOCALIBRATION!!!! ############
#################################################################################
zO3 <- 25
zSO2 <- 0
zNO <- 446
zNOx <- 446
zNOy <- 446
zCO2 <- 118
zH2O <- 340
zFstd <- 0
zFtot <- 0

## create some vectors with default values in them
S.o3 <- rep(4.958,Ldata)
S.so2 <- rep(96.3,Ldata)
S.no <- rep(65,Ldata)
S.nox <- rep(65,Ldata)
S.noy <- rep(65,Ldata)
S.co2 <- rep(3.125,Ldata)
S.h2o <- rep(83.333,Ldata)

############# Now index these by Zero(Z) and Spans (S)
######### need to think about this with the new cal. scheme!!!!!!!
######################################################################

## Index the offsets:
## Zero air offsets:
index.zero <- index[Z1]
index.diff <- diff(index.zero)		## gives big number on last pt. of Zero
index.diff <- c(500,index.diff)		## this shifts the big number to the first point of the zero.
dday.zero <- dday[Z1]
Vo3.zero <- Vo3[Z1]
Vso2.zero <- Vso2[Z1]
Vh2o.zero <- Vh2o[Z1]
Vco2.zero <- Vco2[Z1]
Vno.zero <- Vno[Z1]
VFstd.zero <- VFstd[Z1]
VFtot.zero <- VFtot[Z1]

Fo3.zero <- F.Zo3[Z1]
Fso2.zero <- F.Zso2[Z1]
Fh2o.zero <- F.Zh2o[Z1]
Fco2.zero <- F.Zco2[Z1]
Fno.zero <- F.Zno[Z1]
Lzeros <- length(index.zero)

### now some indexes for Z2,  prerxn zeros of NO/NOy
Vno.zero2 <- Vno[Z2]
index.Z2 <- index[Z2]
index.diff2 <- diff(index.Z2)		## gives big number on last pt. of Zero
index.diff2 <- c(500,index.diff2)		## this shifts the big number to the first point of the zero.
dday.zero2 <- dday[Z2]
Lzeros2 <- length(index.Z2)

## Span indexes:  "S" for NO,SO2,  "S1" for Co2
### Index the Spans
index.span <- index[S]
index.Sdiff <- diff(index.span)		## gives big number on last pt. of Zero
index.Sdiff <- c(500,index.Sdiff)		## this shifts the big number to the first point of the zero.
dday.span <- dday[S]
Vo3.span <- Vo3[S]
Vso2.span <- Vso2[S]
Vno.span <- Vno[S]
VFstd.span <- VFstd[S]
VFtot.span <- VFtot[S]

Fso2.span <- F.Sso2[S]
Fno.span <- F.Sno[S]
SO2c.span <- SO2.cal[S]
NOc.span <- NO.cal[S]
Lspans <- length(index.span)

index.span1 <- index[S1]
index.Sdiff1 <- diff(index.span1)		## gives big number on last pt. of Zero
index.Sdiff1 <- c(500,index.Sdiff1)		## this shifts the big number to the first point of the zero.
dday.span1 <- dday[S1]
Vh2o.span <- Vh2o[S1]
Vco2.span <- Vco2[S1]
Fco2.span <- F.Sco2[S1]
CO2c.span <- CO2.cal[S1]
Lspans1 <- length(index.span1)

#### Start with the offsets #####################################################
#################################################################################
## Now loop until you reach the end of index.zero
####################################################################################
i <- 1
k <- 0
istart <- 1
lastVso2 <- zSO2
cat("\n Start the zero air offset loop!")
#browser()
### - key - to average the 4 1-min data together.

### Zero's first:
while(i <= Lzeros)	{

	if(index.diff[i] > 1.5)	{
		j <- i + 4
		m <- j + 5
		cat("\n Index numbers:",index.diff[i:j],index.zero[i:j])
		#cat("\n Flags:",Fso2.zero[i],Fo3.zero[i],Fno.zero[i],Fco2.zero[i],Fh2o.zero[i])
		
			if(Fo3.zero[i] == 1)	zO3 <- mean(Vo3.zero[(i+1):j],na.rm=T)
			if(Fso2.zero[i] == 1)	zSO2 <- mean(Vso2.zero[(i+1):j],na.rm=T)
			if(Fno.zero[i] == 1)		{
				Vno.i <- Vno.zero[i:j]
				zNO <- Vno.i[1]
				#if(i == 1)	zNO <- Vno.zero[i]		### in case we start during a zero
				zNOx <- mean(Vno.i[2:3],na.rm=T)
				zNOy <- mean(Vno.i[4:5],na.rm=T)
					}
		if(Fco2.zero[i] == 1)	zCO2 <- mean(Vco2.zero[(i+1):j],na.rm=T)		
		if(Fh2o.zero[i] == 1)	zH2O <- mean(Vh2o.zero[(i+1):j],na.rm=T)
		zFstd <- mean(VFstd.zero[i:(j-1)],na.rm=T)
		zFtot <- mean(VFtot.zero[i:(j-1)],na.rm=T)
		
		#if(index.zero[i] > 23220)	{
		#			cat("not updating zSO2!")
		#			browser()
		#				}	

		#if(MONTH=="apr")	{
		#	if(i > 500 & i < 520)	browser()
		#			}		


		######## now look for second NO/NOx zero values (pre-rxn chamber method)
		if(F.Zno[(j+1)] == 1)		{
				V2no <- index.zero[j] + 1
				Vno.2i <- Vno[V2no:(V2no+4)]
				zNO.2 <- Vno.2i[1]
				zNOx.2 <- mean(Vno.2i[2:3],na.rm=T)
				zNOy.2 <- mean(Vno.2i[4:5],na.rm=T)
					}
	
		## create a matrix of the offsets, etc.
		if(k==0)	{
				M <- c(index.zero[i], dday.zero[i],zSO2,zO3,zNO,zNOx,zNOy,zCO2,zH2O,zFstd,zFtot,zNO.2,zNOx.2,zNOy.2)
				k <- k + 1
				}
		else		{
				N <- c(index.zero[i], dday.zero[i],zSO2,zO3,zNO,zNOx,zNOy,zCO2,zH2O,zFstd,zFtot,zNO.2,zNOx.2,zNOy.2)
				M <- rbind(M,N)
				}	
			i <- i + 4
			#lastVso2 <- zSO2in
			}	## end of if(index.diff[i]==T)	

		#lastVso2 <- zSO2
		i <- i + 1
			
		}	## end of while loop.


## quick filters of calibration data:
M[,3] <- ifelse(M[,3] > Fil.so2,mean(M[,3],na.rm=T),M[,3])		# so2
M[,3] <- ifelse(M[,3] < FilB.so2,mean(M[,3],na.rm=T),M[,3])
M[,4] <- ifelse(M[,4] > Fil.o3,median(M[,4],na.rm=T),M[,4])		# O3
M[,8] <- ifelse(M[,8] > Fil.co2,median(M[,8],na.rm=T),M[,8])		# co2
M[,9] <- ifelse(M[,9] > Fil.h2o,median(M[,9],na.rm=T),M[,9])		# h2o

M[,5] <- ifelse(M[,5] < 375,median(M[,5],na.rm=T),M[,5])		#NO
M[,6] <- ifelse(M[,6] < 375,median(M[,6],na.rm=T),M[,6])		#NOx
M[,7] <- ifelse(M[,7] < 375,median(M[,7],na.rm=T),M[,7])		#NOy

M[,5] <- ifelse(M[,5] > 725,median(M[,5],na.rm=T),M[,5])		#NO
M[,6] <- ifelse(M[,6] > 725,median(M[,6],na.rm=T),M[,6])		#NOx
M[,7] <- ifelse(M[,7] > 725,median(M[,7],na.rm=T),M[,7])		#NOy

M[,12] <- ifelse(M[,12] < 375,median(M[,12],na.rm=T),M[,12])		#NO
M[,12] <- ifelse(M[,12] > 725,median(M[,12],na.rm=T),M[,12])
M[,13] <- ifelse(M[,13] < 375,median(M[,13],na.rm=T),M[,13])		#NO
M[,13] <- ifelse(M[,13] > 725,median(M[,13],na.rm=T),M[,13])
M[,14] <- ifelse(M[,14] < 375,median(M[,14],na.rm=T),M[,14])		#NO
M[,14] <- ifelse(M[,14] > 725,median(M[,14],na.rm=T),M[,14])


####### now have matrix (M) of index numbers,ddays and the calculated offsets at that time.
cat("\n End of zero while loop")
#browser()
#### Quick plots:
par(mfrow=c(2,1))
## Plot #1, zeros
plot(M[,2],M[,4],type="o",ylab="O3")
title("Zero's")
plot(M[,2],M[,3],type="o",ylab="SO2")
browser()

plot(M[,2],M[,8],type="o",,ylab="CO2")
plot(M[,2],M[,9],type="o",,ylab="H2O")
browser()
plot(M[,2],M[,5],type="o",ylab="NO")
	points(M[,2],M[,6],lty=5,col=2)
  	points(M[,2],M[,7],lty=2,col=3)
plot(M[,2],M[,12],type="o",ylab="NO")
	points(M[,2],M[,13],lty=5,col=2)
  	points(M[,2],M[,14],lty=2,col=3)

cat("\n End of zeros")
browser()

##### Now use fun.linear extrap to fill in zeros for the long (1 min.) time series.
######### Note - deal with initial points in subroutine!!! #################
Z.o3 <- fun.linearextrap(index,Z.o3,M[,1],M[,4])
Z.so2 <- fun.linearextrap(index,Z.so2,M[,1],M[,3])
Z.no <- fun.linearextrap(index,Z.no,M[,1],M[,5])
Z.nox <- fun.linearextrap(index,Z.nox,M[,1],M[,6])
Z.noy <- fun.linearextrap(index,Z.noy,M[,1],M[,7])
Z.co2 <- fun.linearextrap(index,Z.co2,M[,1],M[,8])
Z.h2o <- fun.linearextrap(index,Z.h2o,M[,1],M[,9])
Z.Fstd <- fun.linearextrap(index,Z.Fstd,M[,1],M[,10])
Z.Ftot <- fun.linearextrap(index,Z.Ftot,M[,1],M[,11])
Z2.no <- fun.linearextrap(index,Z2.no,M[,1],M[,12])
Z2.nox <- fun.linearextrap(index,Z2.nox,M[,1],M[,13])
Z2.noy <- fun.linearextrap(index,Z2.noy,M[,1],M[,14])
######################################################


####################################################
###### Not sure what to do about missing data???
####################################################
	### check the time tags:
	#diff.time <- dday[j] - dday[i]
	### if time difference is > 0.33 (8hours), then re-set k and istart again
	#if(diff.time > 0.33 & k > 1)	{
	#	cat("\n Break in data",i)
	#	k <- 0
	#	istart <- index[i]
	#			}

############# should have zero's now. ###########################################33
cat("\n check on zeros")
#browser()

###########################################################
### filter out some obviously bad values.
## SO2
F.Zso2 <- ifelse(Z.so2 > 250,F,F.Zso2)
F.Zso2 <- ifelse(Z.so2 < -250,F,F.Zso2)
## O3
F.Zo3 <- ifelse(Z.o3 > 100,F,F.Zo3)
F.Zo3 <- ifelse(Z.o3 < 0,F,F.Zo3)
## NO
F.Zno <- ifelse(Z.no > 650,F,F.Zno)
F.Zno <- ifelse(Z.no < 250,F,F.Zno)
## CO2
F.Zco2 <- ifelse(Z.co2 > 150,F,F.Zco2)
F.Zco2 <- ifelse(Z.co2 < 65,F,F.Zco2)
## H2O
F.Zh2o <- ifelse(Z.h2o > 500,F,F.Zh2o)
F.Zh2o <- ifelse(Z.h2o < 130,F,F.Zh2o)
## NO #2
F.Zno2 <- ifelse(Z2.no > 650,F,F.Zno2)
F.Zno2 <- ifelse(Z2.no < 400,F,F.Zno2)
####### then filter ###################################
#### put in medians for calibrations that are flagged as bad. 
Z.o3 <- ifelse(F.Zo3 == F,median(Z.o3,na.rm=T),Z.o3)
Z.so2 <- ifelse(F.Zso2 == F,median(Z.so2,na.rm=T),Z.so2)
Z.no <- ifelse(F.Zno == F,median(Z.no,na.rm=T),Z.no)
Z.nox <- ifelse(F.Zno == F,median(Z.nox,na.rm=T),Z.nox)
Z.noy <- ifelse(F.Zno == F,median(Z.noy,na.rm=T),Z.noy)
Z.co2 <- ifelse(F.Zco2 == F,median(Z.co2,na.rm=T),Z.co2)
Z.h2o <- ifelse(F.Zh2o == F,median(Z.h2o,na.rm=T),Z.h2o)
Z2.no <- ifelse(F.Zno2 == F,median(Z2.no,na.rm=T),Z2.no)
Z2.nox <- ifelse(F.Zno2 == F,median(Z2.nox,na.rm=T),Z2.nox)
Z2.noy <- ifelse(F.Zno2 == F,median(Z2.noy,na.rm=T),Z2.noy)
cat("\n check on zeros")
#browser()

cat("\n Zeros:")
cat("\n Ozone:",median(Z.o3,na.rm=T),"std=",var(Z.o3,na.rm=T)^.5)
cat("\n SO2:",median(Z.so2,na.rm=T),"std=",var(Z.so2,na.rm=T)^.5)
cat("\n NO:",median(Z.no,na.rm=T),"std=",var(Z.no,na.rm=T)^.5)
cat("\n NOx:",median(Z.nox,na.rm=T),"std=",var(Z.nox,na.rm=T)^.5)
cat("\n NOy:",median(Z.noy,na.rm=T),"std=",var(Z.noy,na.rm=T)^.5)
cat("\n CO2:",median(Z.co2,na.rm=T),"std=",var(Z.co2,na.rm=T)^.5)
cat("\n H2O:",median(Z.h2o,na.rm=T),"std=",var(Z.h2o,na.rm=T)^.5)


par(mfrow=c(4,1))
## Plot #1, zeros
plot(date1[,6],Z.o3,type="l",ylab="O3")
title("Zero's")
plot(date1[,6],Z.so2,type="l",ylab="SO2")
plot(date1[,6],Z.co2,type="l",ylab="CO2")
plot(date1[,6],Z.h2o,type="l",ylab="H2O")
browser()

limtest <- cbind(Z.no,Z.nox,Z.noy)
nolim <- range(limtest,na.rm=T)
plot(date1[,6],Z.noy,type="l",ylab="NO",ylim=nolim)
	points(date1[,6],Z.nox,lty=5,col=2)
  	points(date1[,6],Z.no,lty=2,col=3)
limtest <- cbind(Z2.no,Z2.nox,Z2.noy)
nolim <- range(limtest,na.rm=T)
plot(date1[,6],Z2.noy,type="l",ylab="NO #2",ylim=nolim)
	points(date1[,6],Z2.nox,lty=5,col=2)
  	points(date1[,6],Z2.no,lty=2,col=3)

### get difference:  Z(0air) - Z(pre-rxn)
Zdiff.no <- Z.no - Z2.no
Zdiff.nox <- Z.nox - Z2.nox
Zdiff.noy <- Z.noy - Z2.noy
limtest <- cbind(Zdiff.no,Zdiff.nox,Zdiff.noy)
nolim <- range(limtest,na.rm=T)
plot(date1[,6],Zdiff.noy,type="l",ylab="Zero difference, mV",ylim=nolim)
	points(date1[,6],Zdiff.nox,lty=5,col=2)
  	points(date1[,6],Zdiff.no,lty=2,col=3)
cat("\n Zero Diff, NO:",median(Zdiff.no,na.rm=T),"std=",var(Zdiff.no,na.rm=T)^.5)
cat("\n Zero Diff, NOx:",median(Zdiff.nox,na.rm=T),"std=",var(Zdiff.nox,na.rm=T)^.5)
cat("\n Zero Diff, NOy:",median(Zdiff.noy,na.rm=T),"std=",var(Zdiff.noy,na.rm=T)^.5)
cat("\n End of zeros")
################# End of zero's ################################
browser()

################# Start Spans ##################################
####################################################################################
i <- 1
j <- i + 4
k <- 0
istart <- 1

#### Some start values: #########################################################
##### NOTE - THESE WILL BE THE VALUES USED IF NO AUTOCALIBRATION!!!! ############
#################################################################################
sO3 <- 4.958
sSO2 <- 96.3
sNO <- 55
sNOx <- 55
sNOy <- 55
sCO2 <- 3.125
sH2O <- 83.333

############## working on spans!!!!!!!!!!!!
cat("\n REAdy for Span loop - NO and SO2 spans")
#browser()
### go through the big file - line-by-line and create vectors of zero's and spans
### - key - to average the 4 1-min data together.
while(i <= Lspans)	{

	if(index.Sdiff[i] > 1.5)	{
		j <- i + 4
		cat("\n Index numbers:",index.Sdiff[i:j], index.span[i:j])
		## get the actual index value minus 1 (gets last actual sampled value
		i.data <- index.span[i]		## i.data = index number in 1 min. file.
		lastsample <- i.data - 1
		nextsample <- i.data + 5
		#if(i > 20 & i < 30)	browser()
			#### get mean flows:
			## 0-20 sccm MFC  
			## 2/11/09 calibration and 6/11/09 Voltage cal.
			#sFstd <- (mean(VFstd.span[i:j],na.rm=T) + 65.98)/502.4	## Voltage to REadout
			#sFstd <- (sFstd - 0.264)/0.247				## REAdout to sccm// 2/11/09
			# 9/27/11 calibrations - directly from voltage to flow:
			sFstd <- (mean(VFstd.span[i:j],na.rm=T) - Z.Fstd[i.data])/115.8

			## 0-20000 sccm MFC
			## 2/11/09 calibration and 6/11/09 Voltage cal.					
			sFtot.old <- (mean(VFtot.span[i:j],na.rm=T) - Z.Ftot[i.data])*2	## mV from MFC
			sFtot.old <- (sFtot.old/5000)*20000
			# 9/27/11 calibrations - directly from voltage to flow:
			sFtot <- (mean(VFtot.span[i:j],na.rm=T) - Z.Ftot[i.data])/122.4
			sFtot <- sFtot*1000		### convert to sccm from sLpm
			
			#### Test new calibrations equations:
			#if(i == 1)	browser()

			#### get mean signals:
			if(Fso2.span[i] == 1)	{
				sSO2 <- mean(Vso2.span[(i+1):j],na.rm=T) - Vso2[lastsample]
					## subtract ambient SO2 from last level.
								
				SO2.calc <- SO2c.span[i]*(sFstd/sFtot)
				sSO2 <- sSO2/SO2.calc			## mV/ppb
					}
				
			if(Fno.span[i] == 1)		{
		
				sNO <- Vno.span[i] - Vno[nextsample]		## only use the NO channel
										## subtract NO from next Level
				NO.calc <- NOc.span[i]*(sFstd/sFtot)
				sNO <- sNO/NO.calc			## mV/ppb
				sNOx <- sNO
				sNOy <- mean(Vno.span[(j-1):j],na.rm=T) - Vno[lastsample]		## For NOy - subtract NOy from last level
				sNOy <- sNOy/NO.calc
					}	
			

		

		## create a matrix of the offsets, etc.
		if(k==0)	{
				P <- c(index.span[i], dday.span[i],sSO2,sNO,sNOx,sNOy,sFstd,sFtot,sFtot.old)
				k <- k + 1
				}
		else		{
				N <- c(index.span[i], dday.span[i],sSO2,sNO,sNOx,sNOy,sFstd,sFtot,sFtot.old)
				P <- rbind(P,N)
				}	

		i <- i + 4

			}	## end of if(index.diff[i]==T)			
		i <- i + 1

		}	## end of while loop.

######## Now do a loop for CO2 spans using the S1 index.
i <- 1
j <- i + 4
k <- 0
istart <- 1
cat("\n REAdy for Span loop - CO2 spans")
while(i <= Lspans1)	{

	if(index.Sdiff1[i] > 1.5)	{
		j <- i + 4
		cat("\n Index numbers:",index.Sdiff1[i:j], index.span1[i:j])
		## get the actual index value minus 1 (gets last actual sampled value
		i.data <- index.span1[i]		## i.data = index number in 1 min. file.
		lastsample <- i.data - 1
		nextsample <- i.data + 5
		#if(i > 20 & i < 30)	browser()
			if(Fco2.span[i] == 1)	{
				sCO2 <- mean(Vco2.span[(i+1):j],na.rm=T) - Z.co2[i.data]
				sCO2 <- sCO2/CO2c.span[i]			## mV/ppm
					}
		
		## create a matrix of the offsets, etc.
		if(k==0)	{
				D <- c(index.span1[i], dday.span1[i],sCO2)
				k <- k + 1
				}
		else		{
				E <- c(index.span1[i], dday.span1[i],sCO2)
				D <- rbind(D,E)
				}	

		i <- i + 4

			}	## end of if(index.diff[i]==T)			
		i <- i + 1

		}	## end of while loop.


cat("\n check on spans")

### filter out the bad values before linear extrapolation.
P[,3] <- ifelse(P[,3] > 170 | P[,3] < 60,median(P[,3],na.rm=T),P[,3])
P[,4] <- ifelse(P[,4] > 120 | P[,4] < 20,median(P[,4],na.rm=T),P[,4])
P[,5] <- ifelse(P[,5] > 120 | P[,5] < 20,median(P[,5],na.rm=T),P[,5])
P[,6] <- ifelse(P[,6] > 120 | P[,6] < 20,median(P[,6],na.rm=T),P[,6])

if(MONTH=="aug")	{
	## not many NO spans - must get rid of bad spans first before taking the median.
	test <- P[,4] < 20
	P[,4] <- ifelse(test==T,NA,P[,4])
	P[,5] <- ifelse(test==T,NA,P[,5])
	P[,6] <- ifelse(test==T,NA,P[,6])
	test <- is.na(P[,4])
	P[,4] <- ifelse(test==T,median(P[,4],na.rm=T),P[,4])
	P[,5] <- ifelse(test==T,median(P[,4],na.rm=T),P[,5])
	P[,6] <- ifelse(test==T,median(P[,4],na.rm=T),P[,6])	
		}

D[,3] <- ifelse(D[,3] > 3.5 | D[,3] < 2.9,median(D[,3],na.rm=T),D[,3])

####### now have matrix (M) of index numbers,ddays and the calculated offsets at that time.
#### Quick plots:
par(mfrow=c(3,1))
## Plot #1, Spans
plot(P[,2],P[,3],type="o",ylab="SO2")
title("Span's")
plot(D[,2],D[,3],type="o",,ylab="CO2")
plot(P[,2],P[,4],type="o",ylab="NO")
	points(P[,2],P[,5],lty=5,col=2)
  	points(P[,2],P[,6],lty=2,col=3)
cat("\n End of spans")
browser()
##### Now use fun.linear extrap to fill in zeros for the long (1 min.) time series.
######### Note - deal with initial points in subroutine!!! #################
S.so2 <- fun.linearextrap(index,S.so2,P[,1],P[,3])
S.no <- fun.linearextrap(index,S.no,P[,1],P[,4])
S.nox <- fun.linearextrap(index,S.nox,P[,1],P[,5])
S.noy <- fun.linearextrap(index,S.noy,P[,1],P[,6])
S.co2 <- fun.linearextrap(index,S.co2,D[,1],D[,3])
#S.Fstd <- fun.linearextrap(index,S.Fstd,P[,1],P[,8])
#S.Ftot <- fun.linearextrap(index,S.Ftot,P[,1],P[,9])
####################################################
###### Not sure what to do about missing data???
####################################################
### AS filtering is an interative process - filter out some bad values.
## SO2
F.Sso2 <- ifelse(S.so2 > 200,F,F.Sso2)
F.Sso2 <- ifelse(S.so2 < 40,F,F.Sso2)
## NO
F.Sno <- ifelse(S.no > 150,F,F.Sno)
F.Sno <- ifelse(S.no < 20,F,F.Sno)
## CO2
F.Sco2 <- ifelse(S.co2 > 4,F,F.Sco2)
F.Sco2 <- ifelse(S.co2 < 2.3,F,F.Sco2)

####### then filter again ###################################
#### put in medians for calibrations that are flagged as bad. 
S.so2 <- ifelse(F.Sso2 == F,median(S.so2,na.rm=T),S.so2)
S.no <- ifelse(F.Sno == F,median(S.no,na.rm=T),S.no)
S.nox <- ifelse(F.Sno == F,median(S.nox,na.rm=T),S.nox)
S.noy <- ifelse(F.Sno == F,median(S.noy,na.rm=T),S.noy)
S.co2 <- ifelse(F.Sco2 == F,median(S.co2,na.rm=T),S.co2)

cat("\n Spans:")
cat("\n Ozone:",median(S.o3,na.rm=T))
cat("\n SO2:",median(S.so2,na.rm=T),"std=",var(S.so2,na.rm=T)^.5)
cat("\n NO:",median(S.no,na.rm=T),"std=",var(S.no,na.rm=T)^.5)
cat("\n NOx:",median(S.nox,na.rm=T),"std=",var(S.nox,na.rm=T)^.5)
cat("\n NOy:",median(S.noy,na.rm=T),"std=",var(S.noy,na.rm=T)^.5)
cat("\n CO2:",median(S.co2,na.rm=T),"std=",var(S.co2,na.rm=T)^.5)
cat("\n H2O:",median(S.h2o,na.rm=T))

## Plot #1, Span coeff.
par(mfrow=c(3,1))

title("Span Coeff's")
plot(date1[,6],S.so2,type="l",ylab="SO2")
#browser()
plot(date1[,6],S.no,type="l",ylab="NO")
points(date1[,6],S.noy,col=2,cex=0.5)
plot(date1[,6],S.co2,type="l",ylab="CO2")
#browser()
cat("\n S(O3)=",mean(S.o3,na.rm=T))
cat("\n S(H2O)=",mean(S.h2o,na.rm=T))
#plot(date1[,6],S.h2o,type="l",ylab="H2O")
#plot(index,S.o3,type="l",ylab="O3")

cat("\n End of spans")
browser()


#######################################################################################################
### put all together in a file
# columns:  ttag, index, Z.o3,G.o3, Z.so2, G.so2, Z.no, Z.nox, Z.noy, G.no, Z.co2, G.co2,Z.h2o,G.h2o/// not in order...
### Old output files:
#Cal <- cbind(date1,index,Z.o3,Z.so2,Z.no,Z.nox,Z.noy, Z.co2, Z.h2o, S.o3,S.so2,S.no,S.nox,S.noy,S.co2,S.h2o)
### new method with zero air offsets:
Cal <- cbind(date1,index,Z.o3,Z.so2,Z2.no,Z2.nox,Z2.noy, Z.co2, Z.h2o, S.o3,S.so2,S.no,S.nox,S.noy,S.co2,S.h2o,Z.no,Z.nox,Z.noy)

### already have calibration flag file.  may want to update it later.
Flags <- cbind(date1,index,F.Zso2,F.Sso2,F.Zo3,F.Zno,F.Sno,F.Zco2,F.Sco2,F.Zh2o)


Lp <- length(P[,1])
Lm <- length(M[,1])
## Put together the zero and span calculations:
if(Lm == Lp)	{
	CalCoeffs <- cbind(M,P)
	cat("\n CalCoeffs matrix is OK")
		}
else		{
		cat("\n Different number of Spans/Zeros")
		cat("\n Fix Manually")
		browser()
		}

cat("\n make some plots? - will then output files")
cat("\n Output the Cal coeffiicients matrix??,  output CalCoeffs")


browser()

fun.asci(Cal, calfile)
fun.asci(Flags,outflags)

#browser()

}

##############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])


#browser()

plot(t[,1],x[,1],type="l",xax,ylab=Label)
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=1.m","z=5m","z=9m","z=11m","z=14m","z=23m")
#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=1.5m, red=5m, grn=9m, blue=11m, cyano=14m, mag.=23m")


}

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

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("\n Test linear extrap")
#browser()

### do line fit between each pt.
Lz <- length(xx) - 1
i <- 1
### first - fill in first few points with the first value (no extrapolation)
x[1:xx[i]] <- yy[i]


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.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)

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