##################
# Usage
##################

# this should be executed from the command line, e.g., 
#  > R --vanilla --slave --args barks means scale_values single_plots out_prefix file.txt < vowel_normalize.R

##################
# Set up R session
##################

# clear memory
rm(list=ls(all=TRUE))

##################
# Define Scaling Function (only used by some normalization methods)
##################

# the scaling function assumes norm.list is a data.frame with F1 values in
#  col 3 and 5 and F2 values in 4 and 6.  I.e., don't send data with F3 values!

scale.results <- function (norm.list) {

      min.f1 <- 10000
      max.f1 <- -1000
      min.f2 <- 10000
      max.f2 <- -1000
      for(i in 1:length(norm.list)) {
        spkr.min.f1 <- min(c(norm.list[[i]][,3], norm.list[[i]][,5]), na.rm=TRUE)
	if (spkr.min.f1 < min.f1) min.f1 <- spkr.min.f1
        spkr.max.f1 <- max(c(norm.list[[i]][,3], norm.list[[i]][,5]), na.rm=TRUE)
	if (spkr.max.f1 > max.f1) max.f1 <- spkr.max.f1
        spkr.min.f2 <- min(c(norm.list[[i]][,4], norm.list[[i]][,6]), na.rm=TRUE)
	if (spkr.min.f2 < min.f2) min.f2 <- spkr.min.f2
        spkr.max.f2 <- max(c(norm.list[[i]][,4], norm.list[[i]][,6]), na.rm=TRUE)
	if (spkr.max.f2 > max.f2) max.f2 <- spkr.max.f2
      }
      for (i in 1:length(norm.list)) {
        norm.list[[i]][,3]<-round((500*(norm.list[[i]][,3] - min.f1)/(max.f1 - min.f1)) + 250, 3)
        norm.list[[i]][,5]<-round((500*(norm.list[[i]][,5] - min.f1)/(max.f1 - min.f1)) + 250, 3)
        norm.list[[i]][,4]<-round((1400*(norm.list[[i]][,4] - min.f2)/(max.f2 - min.f2)) + 850, 3)
        norm.list[[i]][,6]<-round((1400*(norm.list[[i]][,6] - min.f2)/(max.f2 - min.f2)) + 850, 3)
      }
      norm.list
}

##################
# Define Normalizing Functions
##################

organize.data <- function() {
	return.vals<-data.frame(all.vowels, c(1), f1s, f2s, f3s, f1.gls, f2.gls, f3.gls)
	names(return.vals)<-c("Vowel", "N", "F1", "F2", "F3", "F1 gl", "F2 gl", "F3 gl")
	return.vals
}

compute.means <- function() {

	# determine individual vowel means
	f1.means<-vector(length=length(vowel.types))
	f2.means<-vector(length=length(vowel.types))
	f3.means<-vector(length=length(vowel.types))
	f1.gl.means<-vector(length=length(vowel.types))
	f2.gl.means<-vector(length=length(vowel.types))
	f3.gl.means<-vector(length=length(vowel.types))
	num.of.tokens<-vector(length=length(vowel.types))
	for (i in 1:length(vowel.types)) {
		f1.means[i]<-round(mean(f1s[all.vowels==vowel.types[i]]),3)
		f2.means[i]<-round(mean(f2s[all.vowels==vowel.types[i]]),3)
		f3.means[i]<-round(mean(f3s[all.vowels==vowel.types[i]]),3)
	
		f1.gl.means[i]<-round(mean(f1.gls[all.vowels==vowel.types[i]]),3)
		f2.gl.means[i]<-round(mean(f2.gls[all.vowels==vowel.types[i]]),3)
		f3.gl.means[i]<-round(mean(f3.gls[all.vowels==vowel.types[i]]),3)

		num.of.tokens[i]<-length(all.vowels[all.vowels==vowel.types[i]])
	}
	return.vals<-data.frame(vowel.types, num.of.tokens, f1.means, f2.means, f3.means, f1.gl.means, f2.gl.means, f3.gl.means)
	names(return.vals)<-c("Vowel", "N", "F1", "F2", "F3", "F1 gl", "F2 gl", "F3 gl")
	return.vals
}

compute.bark <- function(use.indiv = NA) {
	if (is.na(use.indiv)) vals <- compute.means()
	else vals <- organize.data()
	
	# convert formant mean values to Barks using the formula:
	#  =((26.81/(1+1960/x))-0.53)
	z1.means<-round((26.81/(1+1960/vals[,3]))-0.53,3)
	z2.means<-round((26.81/(1+1960/vals[,4]))-0.53,3)
	z3.means<-round((26.81/(1+1960/vals[,5]))-0.53,3)

	z1.gl.means<-round((26.81/(1+1960/vals[,6]))-0.53,3)
	z2.gl.means<-round((26.81/(1+1960/vals[,7]))-0.53,3)
	z3.gl.means<-round((26.81/(1+1960/vals[,8]))-0.53,3)

	# determine z3-z1, z3-z2, and z2-z1
	z3.z1<-round(z3.means - z1.means, 3)
	z3.z2<-round(z3.means - z2.means, 3)
	z2.z1<-round(z2.means - z1.means, 3)

	z3.gl.z1<-round(z3.gl.means - z1.gl.means, 3)
	z3.gl.z2<-round(z3.gl.means - z2.gl.means, 3)
	z2.gl.z1<-round(z2.gl.means - z1.gl.means, 3)	
	return.vals<-data.frame(vals[,1], vals[,2], z3.z1, z3.z2, z2.z1, z3.gl.z1, z3.gl.z2, z2.gl.z1)
	names(return.vals)<-c("Vowel", "N", "Z3-Z1", "Z3-Z2", "Z2-Z1", "Z3-Z1 gl", "Z3-Z2 gl", "Z2-Z1 gl")
	return.vals
}

compute.labov <- function(all.mean) {
	
	# do Labov's Nearey-like algorithm
	ln.f1s <- log(f1s)
	ln.f2s <- log(f2s)
	ln.f1.gls <- log(f1.gls)
	ln.f2.gls <- log(f2.gls)
	# set all NA values to be the same as onset, basically act like 
	# all monophthongs are diphs with glide same as onset 
	ln.f1.gls[is.na(ln.f1.gls)]<-ln.f1s[is.na(ln.f1.gls)]
	ln.f2.gls[is.na(ln.f2.gls)]<-ln.f2s[is.na(ln.f2.gls)]
	grand.mean <- mean(c(ln.f1s, ln.f1.gls, ln.f2s, ln.f2.gls), na.rm=TRUE)
        factor <- exp(all.mean - grand.mean)
	scaled.f1s <- factor*f1s
	scaled.f2s <- factor*f2s
	scaled.f1.gls <- factor*f1.gls
	scaled.f2.gls <- factor*f2.gls
	
	# determine individual normalized vowel means
	no.f1.gls <- vector(length=length(vowel.types))
	no.f2.gls <- vector(length=length(vowel.types))
	norm.f1.means <- vector(length=length(vowel.types))
	norm.f2.means <- vector(length=length(vowel.types))
	norm.f1.gl.means <- vector(length=length(vowel.types))
	norm.f2.gl.means <- vector(length=length(vowel.types))
	num.of.tokens <- vector(length=length(vowel.types))
	for (i in 1:length(vowel.types)) {
		
		no.f1.gls[i]<-is.na(mean(f1.gls[all.vowels==vowel.types[i]]))
		no.f2.gls[i]<-is.na(mean(f2.gls[all.vowels==vowel.types[i]]))

		# compute normalized means
		norm.f1.means[i] <- round(mean(scaled.f1s[all.vowels==vowel.types[i]]),3)
		norm.f2.means[i] <- round(mean(scaled.f2s[all.vowels==vowel.types[i]]),3)
		norm.f1.gl.means[i] <- round(mean(scaled.f1.gls[all.vowels==vowel.types[i]]),3)
		norm.f2.gl.means[i] <- round(mean(scaled.f2.gls[all.vowels==vowel.types[i]]),3)

		num.of.tokens[i] <- length(all.vowels[all.vowels==vowel.types[i]])
	}

	# replace all onset values with NAs in glide vectors for monophthongs
	norm.f1.gl.means[no.f1.gls]<-NA
	norm.f2.gl.means[no.f2.gls]<-NA
	
	return.vals<-data.frame(vowel.types, num.of.tokens, norm.f1.means, norm.f2.means, norm.f1.gl.means, norm.f2.gl.means)
	names(return.vals)<-c("Vowel", "N", "F1'", "F2'", "F1' gl", "F2' gl")
	return.vals
}

compute.labov.indiv <- function(all.mean) {
	
	# do Labov's Nearey-like algorithm
	ln.f1s <- log(f1s)
	ln.f2s <- log(f2s)
	ln.f1.gls <- log(f1.gls)
	ln.f2.gls <- log(f2.gls)
	# set all NA values to be the same as onset, basically act like 
	# all monophthongs are diphs with glide same as onset 
	ln.f1.gls[is.na(ln.f1.gls)]<-ln.f1s[is.na(ln.f1.gls)]
	ln.f2.gls[is.na(ln.f2.gls)]<-ln.f2s[is.na(ln.f2.gls)]
	grand.mean <- mean(c(ln.f1s, ln.f1.gls, ln.f2s, ln.f2.gls), na.rm=TRUE)
        factor <- exp(all.mean - grand.mean)
	scaled.f1s <- round(factor*f1s,3)
	scaled.f2s <- round(factor*f2s,3)
	scaled.f1.gls <- round(factor*f1.gls,3)
	scaled.f2.gls <- round(factor*f2.gls,3)
	
	return.vals<-data.frame(all.vowels, c(1), scaled.f1s, scaled.f2s, scaled.f1.gls, scaled.f2.gls)
	names(return.vals)<-c("Vowel", "N", "F1'", "F2'", "F1' gl", "F2' gl")
	return.vals
}

compute.lobanov <- function(f1.all.mean, f2.all.mean) {
	
	if (missing(f1.all.mean)) f1.all.mean <- NA
	if (missing(f2.all.mean)) f2.all.mean <- NA

	# do Lobanov algorithm (subtract mean, divide by stdev, return to normalized range...)
	# for gl, set all NA values to be the same as onset, basically act like 
	# all monophthongs are diphs with glide same as onset 
	alt.f1.gls <- f1.gls
	alt.f2.gls <- f2.gls
	alt.f1.gls[is.na(f1.gls)]<-f1s[is.na(f1.gls)]
	alt.f2.gls[is.na(f2.gls)]<-f2s[is.na(f2.gls)]
	mean.f1 <- mean(c(f1s, alt.f1.gls), na.rm=TRUE)
        if (!is.na(f1.all.mean)) mean.f1<-f1.all.mean
	mean.f2 <- mean(c(f2s, alt.f2.gls), na.rm=TRUE)
        if (!is.na(f2.all.mean)) mean.f2<-f2.all.mean
	stdev.f1 <- sd(c(f1s, alt.f1.gls), na.rm=TRUE)
	stdev.f2 <- sd(c(f2s, alt.f2.gls), na.rm=TRUE)
	norm.f1s <- (f1s-mean.f1)/stdev.f1
	norm.f2s <- (f2s-mean.f2)/stdev.f2
	norm.f1.gls <- (alt.f1.gls-mean.f1)/stdev.f1
	norm.f2.gls <- (alt.f2.gls-mean.f2)/stdev.f2
	sc.f1s <- norm.f1s
	sc.f2s <- norm.f2s
	sc.f1.gls <- norm.f1.gls
	sc.f2.gls <- norm.f2.gls
	#sc.f1s <- (500*(norm.f1s - min(norm.f1s))/(max(norm.f1s) - min(norm.f1s))) + 250
	#sc.f2s <- (1400*(norm.f2s - min(norm.f2s))/(max(norm.f2s) - min(norm.f2s))) + 850
	#sc.f1.gls <- (500*(norm.f1.gls - min(norm.f1.gls))/(max(norm.f1.gls) - min(norm.f1.gls))) + 250
	#sc.f2.gls <- (1400*(norm.f2.gls - min(norm.f2.gls))/(max(norm.f2.gls) - min(norm.f2.gls))) + 850
	
	# determine individual normalized vowel means
	no.f1.gls <- vector(length=length(vowel.types))
	no.f2.gls <- vector(length=length(vowel.types))
	norm.f1.means <- vector(length=length(vowel.types))
	norm.f2.means <- vector(length=length(vowel.types))
	norm.f1.gl.means <- vector(length=length(vowel.types))
	norm.f2.gl.means <- vector(length=length(vowel.types))
	num.of.tokens <- vector(length=length(vowel.types))
	for (i in 1:length(vowel.types)) {
		
		no.f1.gls[i]<-is.na(mean(f1.gls[all.vowels==vowel.types[i]]))
		no.f2.gls[i]<-is.na(mean(f2.gls[all.vowels==vowel.types[i]]))

		# compute normalized means
		norm.f1.means[i] <- round(mean(sc.f1s[all.vowels==vowel.types[i]]),3)
		norm.f2.means[i] <- round(mean(sc.f2s[all.vowels==vowel.types[i]]),3)
		norm.f1.gl.means[i] <- round(mean(sc.f1.gls[all.vowels==vowel.types[i]]),3)
		norm.f2.gl.means[i] <- round(mean(sc.f2.gls[all.vowels==vowel.types[i]]),3)

		num.of.tokens[i] <- length(all.vowels[all.vowels==vowel.types[i]])
	}

	# replace all onset values with NAs in glide vectors for monophthongs
	norm.f1.gl.means[no.f1.gls]<-NA
	norm.f2.gl.means[no.f2.gls]<-NA
	
	return.vals<-data.frame(vowel.types, num.of.tokens, norm.f1.means, norm.f2.means, norm.f1.gl.means, norm.f2.gl.means)
	names(return.vals)<-c("Vowel", "N", "F1<sup>N</sup>", "F2<sup>N</sup>", "F1<sup>N</sup> gl", "F2<sup>N</sup> gl")
	return.vals
}

compute.lobanov.mf <- function(f1.all.mean=NA, f2.all.mean=NA, use.indiv=NA) {

  	# compute means before doing lobanov algorithm
	if (is.na(use.indiv)) vals <- compute.means()
	else vals <- organize.data()
	
	# do Lobanov algorithm on means
	f1.means<-vals[,3]
	f2.means<-vals[,4]
	f1.gl.means<-vals[,6]
	f2.gl.means<-vals[,7]
	f1.gl.means[is.na(f1.gl.means)]<-f1.means[is.na(f1.gl.means)]
	f2.gl.means[is.na(f2.gl.means)]<-f2.means[is.na(f2.gl.means)]
	mean.f1 <- mean(c(f1.means, f1.gl.means))
        if (!is.na(f1.all.mean)) mean.f1<-f1.all.mean
	mean.f2 <- mean(c(f2.means, f2.gl.means))
        if (!is.na(f2.all.mean)) mean.f2<-f2.all.mean
	stdev.f1 <- sd(c(f1.means, f1.gl.means))
	stdev.f2 <- sd(c(f2.means, f2.gl.means))
	norm.f1.means <- (f1.means-mean.f1)/stdev.f1
	norm.f2.means <- (f2.means-mean.f2)/stdev.f2
	norm.f1.gl.means <- (f1.gl.means-mean.f1)/stdev.f1
	norm.f2.gl.means <- (f2.gl.means-mean.f2)/stdev.f2
	sc.f1.means <- round(norm.f1.means, 3)
	sc.f2.means <- round(norm.f2.means, 3)
	sc.f1.gl.means <- round(norm.f1.gl.means, 3)
	sc.f2.gl.means <- round(norm.f2.gl.means, 3)
	#sc.f1.means <- round(((500*(norm.f1.means - min(norm.f1.means))/(max(norm.f1.means) - min(norm.f1.means))) + 250), 3)
	#sc.f2.means <- round(((1400*(norm.f2.means - min(norm.f2.means))/(max(norm.f2.means) - min(norm.f2.means))) + 850), 3)
	#sc.f1.gl.means <- round(((500*(norm.f1.gl.means - min(norm.f1.gl.means))/(max(norm.f1.gl.means) - min(norm.f1.gl.means))) + 250), 3)
	#sc.f2.gl.means <- round(((1400*(norm.f2.gl.means - min(norm.f2.gl.means))/(max(norm.f2.gl.means) - min(norm.f2.gl.means))) + 850), 3)
	sc.f1.gl.means[is.na(vals[,6])]<-NA
	sc.f2.gl.means[is.na(vals[,7])]<-NA
	
	return.vals<-data.frame(vals[,1],vals[,2],sc.f1.means,sc.f2.means,sc.f1.gl.means,sc.f2.gl.means)
	names(return.vals)<-c("Vowel", "N", "F1<sup>N</sup>", "F2<sup>N</sup>", "F1<sup>N</sup> gl", "F2<sup>N</sup> gl")
	return.vals
}

compute.nearey <- function(all.mean = NA, formant.int = NA) {
	
	if (missing(all.mean)) all.mean <- NA
	if (missing(formant.int)) formant.int <- NA
	
	# do Nearey algorithm (convert to natural log, modify, return to normalized range...)
	ln.f1s <- log(f1s)
	ln.f2s <- log(f2s)
	ln.f1.gls <- log(f1.gls)
	ln.f2.gls <- log(f2.gls)
	# set all NA values to be the same as onset, basically act like 
	# all monophthongs are diphs with glide same as onset 
	ln.f1.gls[is.na(ln.f1.gls)]<-ln.f1s[is.na(ln.f1.gls)]
	ln.f2.gls[is.na(ln.f2.gls)]<-ln.f2s[is.na(ln.f2.gls)]
	#mean.ln.f1 <- mean(c(ln.f1s, ln.f1.gls), na.rm=TRUE)
	#mean.ln.f2 <- mean(c(ln.f2s, ln.f2.gls), na.rm=TRUE)
	#exp.ln.f1s <- exp(ln.f1s - mean.ln.f1)
	#exp.ln.f2s <- exp(ln.f2s - mean.ln.f2)
	#exp.ln.f1.gls <- exp(ln.f1.gls - mean.ln.f1)
	#exp.ln.f2.gls <- exp(ln.f2.gls - mean.ln.f2)
	#grand.mean <- mean(c(ln.f1s, ln.f1.gls, ln.f2s, ln.f2.gls), na.rm=TRUE)
	grand.mean.f1 <- mean(c(ln.f1s, ln.f1.gls), na.rm=TRUE)
	grand.mean.f2 <- mean(c(ln.f2s, ln.f2.gls), na.rm=TRUE)
	if (is.na(formant.int)) {
	   grand.mean.f1 <- mean(c(ln.f1s, ln.f1.gls, ln.f2s, ln.f2.gls), na.rm=TRUE)
	   grand.mean.f2 <- mean(c(ln.f1s, ln.f1.gls, ln.f2s, ln.f2.gls), na.rm=TRUE)
	}
        if (!is.na(all.mean)) {
	   grand.mean.f1 <- all.mean
	   grand.mean.f2 <- all.mean
	}
	exp.ln.f1s <- exp(ln.f1s - grand.mean.f1)
	exp.ln.f2s <- exp(ln.f2s - grand.mean.f2)
	exp.ln.f1.gls <- exp(ln.f1.gls - grand.mean.f1)
	exp.ln.f2.gls <- exp(ln.f2.gls - grand.mean.f2)
        scaled.f1s <- exp.ln.f1s
       	scaled.f2s <- exp.ln.f2s
       	scaled.f1.gls <- exp.ln.f1.gls
        scaled.f2.gls <- exp.ln.f2.gls
	#scaled.f1s <- (500*(exp.ln.f1s - min(exp.ln.f1s))/(max(exp.ln.f1s) - min(exp.ln.f1s))) + 250
	#scaled.f2s <- (1400*(exp.ln.f2s - min(exp.ln.f2s))/(max(exp.ln.f2s) - min(exp.ln.f2s))) + 850
	#scaled.f1.gls <- (500*(exp.ln.f1.gls - min(exp.ln.f1.gls))/(max(exp.ln.f1.gls) - min(exp.ln.f1.gls))) + 250
	#scaled.f2.gls <- (1400*(exp.ln.f2.gls - min(exp.ln.f2.gls))/(max(exp.ln.f2.gls) - min(exp.ln.f2.gls))) + 850
	
	# determine individual normalized vowel means
	no.f1.gls <- vector(length=length(vowel.types))
	no.f2.gls <- vector(length=length(vowel.types))
	norm.f1.means <- vector(length=length(vowel.types))
	norm.f2.means <- vector(length=length(vowel.types))
	norm.f1.gl.means <- vector(length=length(vowel.types))
	norm.f2.gl.means <- vector(length=length(vowel.types))
	num.of.tokens <- vector(length=length(vowel.types))
	for (i in 1:length(vowel.types)) {
		
		no.f1.gls[i]<-is.na(mean(f1.gls[all.vowels==vowel.types[i]]))
		no.f2.gls[i]<-is.na(mean(f2.gls[all.vowels==vowel.types[i]]))

		# compute normalized means
		norm.f1.means[i] <- round(mean(scaled.f1s[all.vowels==vowel.types[i]]),3)
		norm.f2.means[i] <- round(mean(scaled.f2s[all.vowels==vowel.types[i]]),3)
		norm.f1.gl.means[i] <- round(mean(scaled.f1.gls[all.vowels==vowel.types[i]]),3)
		norm.f2.gl.means[i] <- round(mean(scaled.f2.gls[all.vowels==vowel.types[i]]),3)

		num.of.tokens[i] <- length(all.vowels[all.vowels==vowel.types[i]])
	}

	# replace all onset values with NAs in glide vectors for monophthongs
	norm.f1.gl.means[no.f1.gls]<-NA
	norm.f2.gl.means[no.f2.gls]<-NA
	
	return.vals<-data.frame(vowel.types, num.of.tokens, norm.f1.means, norm.f2.means, norm.f1.gl.means, norm.f2.gl.means)
	names(return.vals)<-c("Vowel", "N", "F*1", "F*2", "F*1 gl", "F*2 gl")
	return.vals
}

compute.nearey.mf <- function(all.mean = NA, use.indiv = NA, formant.int = NA) {

  	# compute means before doing natural log, etcet...
	if (is.na(use.indiv)) vals <- compute.means()
	else vals <- organize.data()

	# do Nearey algorithm on the means 
	ln.f1.means <- log(vals[,3])
	ln.f2.means <- log(vals[,4])
	ln.f1.gl.means <- log(vals[,6])
	ln.f2.gl.means <- log(vals[,7])
	# set all NA values to be the same as onset, basically act like 
	# all monophthongs are diphs with glide same as onset 
	ln.f1.gl.means[is.na(ln.f1.gl.means)]<-ln.f1.means[is.na(ln.f1.gl.means)]
	ln.f2.gl.means[is.na(ln.f2.gl.means)]<-ln.f2.means[is.na(ln.f2.gl.means)]
	grand.mean.f1 <- mean(c(ln.f1.means, ln.f1.gl.means), na.rm=TRUE)
	grand.mean.f2 <- mean(c(ln.f2.means, ln.f2.gl.means), na.rm=TRUE)
	if (is.na(formant.int)) {
	   grand.mean.f1 <- mean(c(ln.f1.means, ln.f1.gl.means, ln.f2.means, ln.f2.gl.means), na.rm=TRUE)
	   grand.mean.f2 <- mean(c(ln.f1.means, ln.f1.gl.means, ln.f2.means, ln.f2.gl.means), na.rm=TRUE)
	}
        if (!is.na(all.mean)) {
	   grand.mean.f1 <- all.mean
	   grand.mean.f2 <- all.mean
	}
	#grand.mean <- mean(c(ln.f1.means, ln.f1.gl.means, ln.f2.means, ln.f2.gl.means), na.rm=TRUE)
        #if (!is.na(all.mean)) grand.mean<-all.mean
	exp.ln.f1.means <- exp(ln.f1.means - grand.mean.f1)
	exp.ln.f2.means <- exp(ln.f2.means - grand.mean.f2)
	exp.ln.f1.gl.means <- exp(ln.f1.gl.means - grand.mean.f1)
	exp.ln.f2.gl.means <- exp(ln.f2.gl.means - grand.mean.f2)
	scaled.f1.means <- round(exp.ln.f1.means, 3)
	scaled.f2.means <- round(exp.ln.f2.means, 3)
	scaled.f1.gl.means <- round(exp.ln.f1.gl.means, 3)
	scaled.f2.gl.means <- round(exp.ln.f2.gl.means, 3)
	#scaled.f1.means <- round(((500*(exp.ln.f1.means - min(exp.ln.f1.means))/(max(exp.ln.f1.means) - min(exp.ln.f1.means))) + 250), 3)
	#scaled.f2.means <- round(((1400*(exp.ln.f2.means - min(exp.ln.f2.means))/(max(exp.ln.f2.means) - min(exp.ln.f2.means))) + 850), 3)
	#scaled.f1.gl.means <- round(((500*(exp.ln.f1.gl.means - min(exp.ln.f1.gl.means))/(max(exp.ln.f1.gl.means) - min(exp.ln.f1.gl.means))) + 250), 3)
	#scaled.f2.gl.means <- round(((1400*(exp.ln.f2.gl.means - min(exp.ln.f2.gl.means))/(max(exp.ln.f2.gl.means) - min(exp.ln.f2.gl.means))) + 850), 3)

	# reset NA glides to NA
	scaled.f1.gl.means[is.na(vals[,6])]<-NA
	scaled.f2.gl.means[is.na(vals[,7])]<-NA
	
	return.vals<-data.frame(vals[,1],vals[,2],scaled.f1.means,scaled.f2.means,scaled.f1.gl.means,scaled.f2.gl.means)
	names(return.vals)<-c("Vowel", "N", "F*1", "F*2", "F*1 gl", "F*2 gl")
	return.vals
}

compute.nearey.wf3 <- function(all.mean = NA, use.indiv = NA, formant.int = NA){

  	# compute means before doing natural log, etcet...
	if (is.na(use.indiv)) vals <- compute.means()
	else vals <- organize.data()

	# do Nearey algorithm on the value
	ln.f1.means <- log(vals[,3])
	ln.f2.means <- log(vals[,4])
	ln.f3.means <- log(vals[,5])
	ln.f1.gl.means <- log(vals[,6])
	ln.f2.gl.means <- log(vals[,7])
	ln.f3.gl.means <- log(vals[,8])
	# set all NA values to be the same as onset, basically act like 
	# all monophthongs are diphs with glide same as onset 
	ln.f1.gl.means[is.na(ln.f1.gl.means)]<-ln.f1.means[is.na(ln.f1.gl.means)]
	ln.f2.gl.means[is.na(ln.f2.gl.means)]<-ln.f2.means[is.na(ln.f2.gl.means)]
	ln.f3.gl.means[is.na(ln.f3.gl.means)]<-ln.f3.means[is.na(ln.f3.gl.means)]
	grand.mean.f1 <- mean(c(ln.f1.means, ln.f1.gl.means), na.rm=TRUE)
	grand.mean.f2 <- mean(c(ln.f2.means, ln.f2.gl.means), na.rm=TRUE)
	grand.mean.f3 <- mean(c(ln.f3.means, ln.f3.gl.means), na.rm=TRUE)
	if (is.na(formant.int)) {
	   grand.mean.f1 <- mean(c(ln.f1.means, ln.f1.gl.means, ln.f2.means, ln.f2.gl.means, ln.f3.means, ln.f3.gl.means), na.rm=TRUE)
	   grand.mean.f2 <- mean(c(ln.f1.means, ln.f1.gl.means, ln.f2.means, ln.f2.gl.means, ln.f3.means, ln.f3.gl.means), na.rm=TRUE)
	   grand.mean.f3 <- mean(c(ln.f1.means, ln.f1.gl.means, ln.f2.means, ln.f2.gl.means, ln.f3.means, ln.f3.gl.means), na.rm=TRUE)
	}
        if (!is.na(all.mean)) {
	   grand.mean.f1 <- all.mean
	   grand.mean.f2 <- all.mean
	   grand.mean.f3 <- all.mean
	}
	#grand.mean <- mean(c(ln.f1.means, ln.f1.gl.means, ln.f2.means, ln.f2.gl.means, ln.f3.means, ln.f3.gl.means), na.rm=TRUE)
        #if (!is.na(all.mean)) grand.mean<-all.mean
	exp.ln.f1.means <- exp(ln.f1.means - grand.mean.f1)
	exp.ln.f2.means <- exp(ln.f2.means - grand.mean.f2)
	exp.ln.f3.means <- exp(ln.f3.means - grand.mean.f3)
	exp.ln.f1.gl.means <- exp(ln.f1.gl.means - grand.mean.f1)
	exp.ln.f2.gl.means <- exp(ln.f2.gl.means - grand.mean.f2)
	exp.ln.f3.gl.means <- exp(ln.f3.gl.means - grand.mean.f3)
	scaled.f1.means <- round(exp.ln.f1.means, 3)
	scaled.f2.means <- round(exp.ln.f2.means, 3)
	scaled.f3.means <- round(exp.ln.f3.means, 3)
	scaled.f1.gl.means <- round(exp.ln.f1.gl.means, 3)
	scaled.f2.gl.means <- round(exp.ln.f2.gl.means, 3)
	scaled.f3.gl.means <- round(exp.ln.f3.gl.means, 3)

	# reset NA glides to NA
	scaled.f1.gl.means[is.na(vals[,6])]<-NA
	scaled.f2.gl.means[is.na(vals[,7])]<-NA
	scaled.f3.gl.means[is.na(vals[,8])]<-NA
	
	return.vals<-data.frame(vals[,1],vals[,2],scaled.f1.means,scaled.f2.means,scaled.f3.means,scaled.f1.gl.means,scaled.f2.gl.means,scaled.f3.gl.means)
	names(return.vals)<-c("Vowel", "N", "F*1", "F*2", "F*3", "F*1 gl", "F*2 gl", "F*3 gl")
	return.vals
}

compute.wattfabricius <- function() {
        
	# do Watt & Fabricius algorithm

	# get idealized min&max F1 and F2, according to Watt & Fabricius
	
	mean.vals <- compute.means()
	
	# unlike W&F (2003), we determine the corners of the vowel triangle
	#  by choosing the vowels automatically...
	#f2.mid<-min(mean.vals[,4])+((max(mean.vals[,4])-min(mean.vals[,4]))/2)
	#fleece.vowel<-as.character(mean.vals[,1][mean.vals[,3]==min(mean.vals[,3][mean.vals[,4]>f2.mid])])
	#trap.vowel<-as.character(mean.vals[,1][mean.vals[,3]==max(mean.vals[,3])])
	#school.vowel<-as.character(mean.vals[,1][mean.vals[,3]==min(mean.vals[,3][mean.vals[,4]<f2.mid])])
	#i.f1 <- mean.vals[,3][mean.vals[,1]==fleece.vowel]
	#i.f2 <- mean.vals[,4][mean.vals[,1]==fleece.vowel]
	#a.f1 <- mean.vals[,3][mean.vals[,1]==trap.vowel]
	#a.f2 <- mean.vals[,4][mean.vals[,1]==trap.vowel]
	#ul.f1 <- mean.vals[,3][mean.vals[,1]==school.vowel]
	#ul.f2 <- mean.vals[,4][mean.vals[,1]==school.vowel]

        i.f1 <- min(c(mean.vals[,3], mean.vals[,6]),na.rm=TRUE)
	i.f2 <- max(c(mean.vals[,4], mean.vals[,7]),na.rm=TRUE)
	a.f1 <- max(c(mean.vals[,3], mean.vals[,6]),na.rm=TRUE)
	a.f2 <- mean.vals[,4][mean.vals[,3]==a.f1]
	ul.f1 <- i.f1
	ul.f2 <- i.f1

	# compute S values
	S.f1 <- (i.f1 + a.f1 + ul.f1)/3
	S.f2 <- (i.f2 + a.f2 + ul.f2)/3

	# divide f1 and f2 means by S values
	norm.f1s <- round(mean.vals[,3]/S.f1, 3)
	norm.f2s <- round(mean.vals[,4]/S.f2, 3)
	norm.f1.gls <- round(mean.vals[,6]/S.f1, 3)
	norm.f2.gls <- round(mean.vals[,7]/S.f2, 3)

	return.val <- data.frame(mean.vals[,1], mean.vals[,2], norm.f1s, norm.f2s, norm.f1.gls, norm.f2.gls)
	names(return.val) <- c("Vowel", "N", "F1/S(F1)", "F2/S(F2)", "F1g/S(F1g)", "F2g/S(F2g)")
	return.val
}

##################
# Define Plotting Functions
##################

# first, set up the plot speaker-colors globally
pl.c<-c("#ff0000", "#0000ff", "#00ff00", "#ff9900", "#cd00cd","#00cdcd","#cdcd00", "#800000", "#008000","#333399","#ffff00","#808080","#cc99ff")
#pl.c<-c("#cd0000", "#0000cd", "#008b00", "#daa520", "#9400d3", "#8b4513", "#00cdcd", "#556b2f", "#999999")
#pl.c<-c("#3300cc", "#cc3300", "#33cc00", "#997300", "#ffd147", "#9900cc")
# and the points for each speaker
pl.p<-c(15, 16, 17, 18, 7, 8, 6, 14)
# and the default point size for plotting
cex.set <- 1.5

setup.plot <- function(vals, hasf3s) {

  col.plus <- 0
  if (hasf3s) col.plus <- 1  

  # determine scale for axes
  xstart <- -1000
  xend <- 100000
  ystart <- -1000
  yend <- 100000
  for (i in 1:length(vals)) {
     if (xstart < max(vals[[i]][,4], vals[[i]][,(6+col.plus)], na.rm=TRUE)) xstart <- max(vals[[i]][,4], vals[[i]][,(6+col.plus)], na.rm=TRUE)
     if (xend > min(vals[[i]][,4], vals[[i]][,(6+col.plus)], na.rm=TRUE)) xend <- min(vals[[i]][,4], vals[[i]][,(6+col.plus)], na.rm=TRUE)
     if (ystart < max(vals[[i]][,3], vals[[i]][,(5+col.plus)], na.rm=TRUE)) ystart <- max(vals[[i]][,3], vals[[i]][,(5+col.plus)], na.rm=TRUE)
     if (yend > min(vals[[i]][,3], vals[[i]][,(5+col.plus)], na.rm=TRUE)) yend <- min(vals[[i]][,3], vals[[i]][,(5+col.plus)], na.rm=TRUE)
  }
  xunit <- (xstart-xend)/12
  xstart <- xstart+xunit
  xend <- xend-xunit
  yunit <- (ystart-yend)/12
  ystart <- ystart+yunit
  yend <- yend-yunit

  c(xstart, xend, ystart, yend, xunit, yunit)
}

plot.indiv <- function(vals, spkrs) {

  axes <- setup.plot(vals, TRUE)
  cex.set <- 1
  if ((length(spkrs)==1) && (length(unique(vals[[1]][,1])) >= length(pl.c))) {
     pl.c <- rainbow(length(unique(vals[[1]][,1])))
     l.c <- c("black")
  } else if (length(spkrs)==1) {
     l.c <- c("black")  
  } else {
     l.c <- pl.c
  }

  for (i in 1:length(vals)) {
     if (i > 1) par(new=TRUE)
     if (length(spkrs)>1) j=i else j=match(vals[[i]][,1],unique(vals[[i]][,1]))
     plot(vals[[i]][,4], vals[[i]][,3], xlim=c(axes[1],axes[2]), ylim=c(axes[3],axes[4]), xlab="F2", ylab="F1", pch=pl.p[i], cex=cex.set, cex.main=2.0, cex.axis=1.5, cex.lab=1.25, main="Individual vowel formant values", col=pl.c[j])
     legend("bottomleft", legend=spkrs, col=l.c, text.col=l.c, pch=pl.p, inset=.02)
     arrows(vals[[i]][,4], vals[[i]][,3], vals[[i]][,7], vals[[i]][,6], angle=15, length=0.1, lwd=(cex.set+0.3), col=pl.c[j])
     text(vals[[i]][,4], vals[[i]][,3], vals[[i]][,1], adj=c(0,1.5), cex=(2*cex.set/3), col=pl.c[j])
   }
}

plot.means <- function(vals, spkrs) {

  axes <- setup.plot(vals, TRUE)
  if (length(spkrs)==1) pl.c <- c("black")
  if (length(spkrs)>=4) cex.set <- 1

  for (i in 1:length(vals)) {
     if (i > 1) par(new=TRUE)
     plot(vals[[i]][,4], vals[[i]][,3], xlim=c(axes[1],axes[2]), ylim=c(axes[3],axes[4]), xlab="F2", ylab="F1", pch=pl.p[i], cex=cex.set, cex.main=2.0, cex.axis=1.5, cex.lab=1.25, main="Mean formant values", col=pl.c[i])
     legend("bottomleft", legend=spkrs, col=pl.c, text.col=pl.c, pch=pl.p, inset=.02)
     arrows(vals[[i]][,4], vals[[i]][,3], vals[[i]][,7], vals[[i]][,6], angle=15, length=0.1, lwd=(cex.set+0.3), col=pl.c[i])
     text(vals[[i]][,4], vals[[i]][,3], vals[[i]][,1], adj=c(0,1.5), cex=(2*cex.set/3), col=pl.c[i])
   }
}

plot.barks <- function(vals, spkrs) {

  axes <- setup.plot(vals, TRUE)
  if (length(spkrs)==1) pl.c <- c("black")
  if (length(spkrs)>=4) cex.set <- 1
	
  for (i in 1:length(vals)) {
     if (i > 1) par(new=TRUE)
     plot(vals[[i]][,4], vals[[i]][,3], xlim=c(axes[2],axes[1]), ylim=c(axes[4],axes[3]), xlab="Z3-Z2", ylab="Z3-Z1", pch=pl.p[i], cex=cex.set, cex.main=2.0, cex.axis=1.5, cex.lab=1.25, main="Bark difference normalized formant values", col=pl.c[i])
     legend("bottomleft", legend=spkrs, col=pl.c, text.col=pl.c, pch=pl.p, inset=.02)
     arrows(vals[[i]][,4], vals[[i]][,3], vals[[i]][,7], vals[[i]][,6], angle=15, length=0.1, lwd=(cex.set+0.3), col=pl.c[i])
     text(vals[[i]][,4], vals[[i]][,3], vals[[i]][,1], adj=c(0,1.5), cex=(2*cex.set/3), col=pl.c[i])
   }
}

plot.barks.indiv <- function(vals, spkrs) {

  axes <- setup.plot(vals, TRUE)
  cex.set <- 1
  if ((length(spkrs)==1) && (length(unique(vals[[1]][,1])) >= length(pl.c))) {
     pl.c <- rainbow(length(unique(vals[[1]][,1])))
     l.c <- c("black")
  } else if (length(spkrs)==1) {
     l.c <- c("black")  
  } else {
     l.c <- pl.c
  }

  for (i in 1:length(vals)) {
     if (i > 1) par(new=TRUE)
     if (length(spkrs)>1) j=i else j=match(vals[[i]][,1],unique(vals[[i]][,1]))
     plot(vals[[i]][,4], vals[[i]][,3], xlim=c(axes[2],axes[1]), ylim=c(axes[4],axes[3]), xlab="Z3-Z2", ylab="Z3-Z1", pch=pl.p[i], cex=cex.set, cex.main=2.0, cex.axis=1.5, cex.lab=1.25, main="Bark difference normalized formant values", sub="individual vowels", col=pl.c[j])
     legend("bottomleft", legend=spkrs, col=l.c, text.col=l.c, pch=pl.p, inset=.02)
     arrows(vals[[i]][,4], vals[[i]][,3], vals[[i]][,7], vals[[i]][,6], angle=15, length=0.1, lwd=(cex.set+0.3), col=pl.c[j])
     text(vals[[i]][,4], vals[[i]][,3], vals[[i]][,1], adj=c(0,1.5), cex=(2*cex.set/3), col=pl.c[j])
   }
}

plot.labov <- function(vals, spkrs, sub.title) {

  if (missing(sub.title)) sub.title <- NA  

  axes <- setup.plot(vals, FALSE)
  if (length(spkrs)==1) pl.c <- c("black")
  if (length(spkrs)>=4) cex.set <- 1

  for (i in 1:length(vals)) {
     if (i > 1) par(new=TRUE)
     plot(vals[[i]][,4], vals[[i]][,3], xlim=c(axes[1],axes[2]), ylim=c(axes[3],axes[4]), xlab="F2'", ylab="F1'", pch=pl.p[i], cex=cex.set, cex.main=2.0, cex.axis=1.5, cex.lab=1.25, main="Labov ANAE normalized formant values", sub=sub.title, col=pl.c[i])
     legend("bottomleft", legend=spkrs, col=pl.c, text.col=pl.c, pch=pl.p, inset=.02)
     arrows(vals[[i]][,4], vals[[i]][,3], vals[[i]][,6], vals[[i]][,5], angle=15, length=0.1, lwd=(cex.set+0.3), col=pl.c[i])
     text(vals[[i]][,4], vals[[i]][,3], vals[[i]][,1], adj=c(0,1.5), cex=(2*cex.set/3), col=pl.c[i])
   }
}

plot.labov.indiv <- function(vals, spkrs, sub.title) {

  if (missing(sub.title)) sub.title <- NA  

  axes <- setup.plot(vals, FALSE)
  cex.set <- 1
  if ((length(spkrs)==1) && (length(unique(vals[[1]][,1])) >= length(pl.c))) {
     pl.c <- rainbow(length(unique(vals[[1]][,1])))
     l.c <- c("black")
  } else if (length(spkrs)==1) {
     l.c <- c("black")  
  } else {
     l.c <- pl.c
  }

  for (i in 1:length(vals)) {
     if (i > 1) par(new=TRUE)
     if (length(spkrs)>1) j=i else j=match(vals[[i]][,1],unique(vals[[i]][,1]))
     plot(vals[[i]][,4], vals[[i]][,3], xlim=c(axes[1],axes[2]), ylim=c(axes[3],axes[4]), xlab="F2'", ylab="F1'", pch=pl.p[i], cex=cex.set, cex.main=2.0, cex.axis=1.5, cex.lab=1.25, main="Labov ANAE normalized formant values", sub=sub.title, col=pl.c[j])
     legend("bottomleft", legend=spkrs, col=l.c, text.col=l.c, pch=pl.p, inset=.02)
     arrows(vals[[i]][,4], vals[[i]][,3], vals[[i]][,6], vals[[i]][,5], angle=15, length=0.1, lwd=(cex.set+0.3), col=pl.c[j])
     text(vals[[i]][,4], vals[[i]][,3], vals[[i]][,1], adj=c(0,1.5), cex=(2*cex.set/3), col=pl.c[j])
   }
}

plot.lobanov <- function(vals, spkrs, sub.title) {

  axes <- setup.plot(vals, FALSE)
  if (length(spkrs)==1) pl.c <- c("black")
  if (length(spkrs)>=4) cex.set <- 1

  for (i in 1:length(vals)) {
     if (i > 1) par(new=TRUE)
     plot(vals[[i]][,4], vals[[i]][,3], xlim=c(axes[1],axes[2]), ylim=c(axes[3],axes[4]), xlab=expression(F2^N), ylab=expression(F1^N), pch=pl.p[i], cex=cex.set, cex.main=2.0, cex.axis=1.5, cex.lab=1.25, main="Lobanov normalized formant values", sub=sub.title, col=pl.c[i])
     legend("bottomleft", legend=spkrs, col=pl.c, text.col=pl.c, pch=pl.p, inset=.02)
     arrows(vals[[i]][,4], vals[[i]][,3], vals[[i]][,6], vals[[i]][,5], angle=15, length=0.1, lwd=(cex.set+0.3), col=pl.c[i])
     text(vals[[i]][,4]+(i*axes[5]/20), vals[[i]][,3], vals[[i]][,1], adj=c(0,1.5), cex=(2*cex.set/3), col=pl.c[i])
   }
}

plot.lobanov.indiv <- function(vals, spkrs, sub.title) {

  axes <- setup.plot(vals, FALSE)
  cex.set <- 1
  if ((length(spkrs)==1) && (length(unique(vals[[1]][,1])) >= length(pl.c))) {
     pl.c <- rainbow(length(unique(vals[[1]][,1])))
     l.c <- c("black")
  } else if (length(spkrs)==1) {
     l.c <- c("black")  
  } else {
     l.c <- pl.c
  }

  for (i in 1:length(vals)) {
     if (i > 1) par(new=TRUE)
     if (length(spkrs)>1) j=i else j=match(vals[[i]][,1],unique(vals[[i]][,1]))
     plot(vals[[i]][,4], vals[[i]][,3], xlim=c(axes[1],axes[2]), ylim=c(axes[3],axes[4]), xlab=expression(F2^N), ylab=expression(F1^N), pch=pl.p[i], cex=cex.set, cex.main=2.0, cex.axis=1.5, cex.lab=1.25, main="Lobanov normalized formant values", sub=sub.title, col=pl.c[j])
     legend("bottomleft", legend=spkrs, col=l.c, text.col=l.c, pch=pl.p, inset=.02)
     arrows(vals[[i]][,4], vals[[i]][,3], vals[[i]][,6], vals[[i]][,5], angle=15, length=0.1, lwd=(cex.set+0.3), col=pl.c[j])
     text(vals[[i]][,4]+(i*axes[5]/20), vals[[i]][,3], vals[[i]][,1], adj=c(0,1.5), cex=(2*cex.set/3), col=pl.c[j])
   }

}

plot.nearey <- function(vals, spkrs, sub.title, version) {

  axes <- setup.plot(vals, FALSE)
  if (length(spkrs)==1) pl.c <- c("black")
  if (length(spkrs)>=4) cex.set <- 1

  main.title <- "Nearey2 normalized formant values"
  if (version==1) main.title <- "Nearey1 normalized formant values"
  for (i in 1:length(vals)) {
     if (i > 1) par(new=TRUE)
     plot(vals[[i]][,4], vals[[i]][,3], xlim=c(axes[1],axes[2]), ylim=c(axes[3],axes[4]), xlab="F*2", ylab="F*1", pch=pl.p[i], cex=cex.set, cex.main=2.0, cex.axis=1.5, cex.lab=1.25, main=main.title, sub=sub.title, col=pl.c[i])
     legend("bottomleft", legend=spkrs, col=pl.c, text.col=pl.c, pch=pl.p, inset=.02)
     arrows(vals[[i]][,4], vals[[i]][,3], vals[[i]][,6], vals[[i]][,5], angle=15, length=0.1, lwd=(cex.set+0.3), col=pl.c[i])
     text(vals[[i]][,4]+(i*axes[5]/20), vals[[i]][,3], vals[[i]][,1], adj=c(0,1.5), cex=(2*cex.set/3), col=pl.c[i])
   }

}

plot.nearey.indiv <- function(vals, spkrs, sub.title, version) {

  axes <- setup.plot(vals, FALSE)
  cex.set <- 1
  if ((length(spkrs)==1) && (length(unique(vals[[1]][,1])) >= length(pl.c))) {
     pl.c <- rainbow(length(unique(vals[[1]][,1])))
     l.c <- c("black")
  } else if (length(spkrs)==1) {
     l.c <- c("black")  
  } else {
     l.c <- pl.c
  }

  main.title <- "Nearey2 normalized formant values"
  if (version==1) main.title <- "Nearey1 normalized formant values"
  for (i in 1:length(vals)) {
     if (i > 1) par(new=TRUE)
     if (length(spkrs)>1) j=i else j=match(vals[[i]][,1],unique(vals[[i]][,1]))
     plot(vals[[i]][,4], vals[[i]][,3], xlim=c(axes[1],axes[2]), ylim=c(axes[3],axes[4]), xlab="F*2", ylab="F*1", pch=pl.p[i], cex=cex.set, cex.main=2.0, cex.axis=1.5, cex.lab=1.25, main=main.title, sub=sub.title, col=pl.c[j])
     legend("bottomleft", legend=spkrs, col=l.c, text.col=l.c, pch=pl.p, inset=.02)
     arrows(vals[[i]][,4], vals[[i]][,3], vals[[i]][,6], vals[[i]][,5], angle=15, length=0.1, lwd=(cex.set+0.3), col=pl.c[j])
     text(vals[[i]][,4]+(i*axes[5]/20), vals[[i]][,3], vals[[i]][,1], adj=c(0,1.5), cex=(2*cex.set/3), col=pl.c[j])
   }

}

plot.nearey.wf3 <- function(vals, spkrs, sub.title, version) {

  axes <- setup.plot(vals, TRUE)
  cex.set <- 1
  if (length(spkrs)==1) {
     pl.c <- rainbow(length(unique(vals[[1]][,1])))
     l.c <- c("black")
  } else {
     l.c <- pl.c
  }

  main.title <- "Nearey2 normalized formant values"
  if (version==1) main.title <- "Nearey1 normalized formant values"
  for (i in 1:length(vals)) {
     if (i > 1) par(new=TRUE)
     if (length(spkrs)>1) j=i else j=match(vals[[i]][,1],unique(vals[[i]][,1]))
     plot(vals[[i]][,4], vals[[i]][,3], xlim=c(axes[1],axes[2]), ylim=c(axes[3],axes[4]), xlab="F*2", ylab="F*1", pch=pl.p[i], cex=cex.set, cex.main=2.0, cex.axis=1.5, cex.lab=1.25, main=main.title, sub=sub.title, col=pl.c[j])
     legend("bottomleft", legend=spkrs, col=l.c, text.col=l.c, pch=pl.p, inset=.02)
     arrows(vals[[i]][,4], vals[[i]][,3], vals[[i]][,7], vals[[i]][,6], angle=15, length=0.1, lwd=(cex.set+0.3), col=pl.c[j])
     text(vals[[i]][,4]+(i*axes[5]/20), vals[[i]][,3], vals[[i]][,1], adj=c(0,1.5), cex=(2*cex.set/3), col=pl.c[j])
   }

}

plot.wattfabricius <- function(vals, spkrs) {

  axes <- setup.plot(vals, FALSE)
  if (length(spkrs)==1) pl.c <- c("black")
  if (length(spkrs)>=4) cex.set <- 1

  for (i in 1:length(vals)) {
     if (i > 1) par(new=TRUE)
     plot(vals[[i]][,4], vals[[i]][,3], xlim=c(axes[1],axes[2]), ylim=c(axes[3],axes[4]), xlab="F2/S(F2)", ylab="F1/S(F1)", pch=pl.p[i], cex=cex.set, cex.main=2.0, cex.axis=1.5, cex.lab=1.25, main="Watt & Fabricius normalized formant values", col=pl.c[i])
     legend("bottomleft", legend=spkrs, col=pl.c, text.col=pl.c, pch=pl.p, inset=.02)
     arrows(vals[[i]][,4], vals[[i]][,3], vals[[i]][,6], vals[[i]][,5], angle=15, length=0.1, lwd=(cex.set+0.3), col=pl.c[i])
     text(vals[[i]][,4], vals[[i]][,3], vals[[i]][,1], adj=c(0,1.5), cex=(2*cex.set/3), col=pl.c[i])
   }
}

##################
# Define script parameters and read in data
##################

# first attempt to get parameters, like file.name, from script args
cmd.args <- commandArgs(trailingOnly = TRUE)
left.method <- cmd.args[1]
right.method <- cmd.args[2]
scale.values <- 0
if (cmd.args[3]==1) scale.values <- 1
single.plots <- 1
if (cmd.args[4]==0) single.plots <- 0
out.prefix <- paste(getwd(),"/",cmd.args[5], sep="")
file.name <- paste(getwd(),"/",cmd.args[6], sep="")

# are methods not set (like if we are running this through the R GUI)?
# if so, load some defaults
if (is.na(left.method)) left.method<-"means"
if (is.na(right.method)) right.method<-"nearey"
# does file exist and is readable? if not, try prompting user
if (!file.exists(file.name)) file.name <- file.choose()
#if (is.na(out.prefix)) { out.prefix<-file.name

vdata<-read.table(file.name, header=TRUE, sep="\t", quote="", comment.char="", blank.lines.skip=TRUE, fill=TRUE)

spkr.names<-vector()
left.norms<-list()
right.norms<-list()
j<-1
for(spkr in unique(as.vector(vdata[,1]))) {
  # vdata<-all.data[,2:length(names(all.data))][all.data[,1]==spkr]
	
  spkr.names[j]<-spkr
  all.vowels<-as.vector(vdata[,2][vdata[,1]==spkr])
  vowel.types<-unique(all.vowels)

  f1s<-as.vector(vdata[,4][vdata[,1]==spkr])
  f2s<-as.vector(vdata[,5][vdata[,1]==spkr])
  f3s<-as.vector(vdata[,6][vdata[,1]==spkr])

  f1.gls<-as.vector(vdata[,7][vdata[,1]==spkr])
  f2.gls<-as.vector(vdata[,8][vdata[,1]==spkr])
  f3.gls<-as.vector(vdata[,9][vdata[,1]==spkr])

  ##################
 ## Main Processing
  ##################

  normed <- switch(left.method,
		   indiv = organize.data(),
	           means = compute.means(),
		   barks = compute.bark(),
		   barks_indiv = compute.bark(use.indiv=TRUE),
                   labov = compute.labov(mean(log(c(vdata[,4],vdata[,5])))),
		   labov_g=compute.labov(6.896874),
		   labov_indiv=compute.labov.indiv(mean(log(c(vdata[,4],vdata[,5])))),
		   labov_indiv_g=compute.labov.indiv(6.896874),
		   lobanov =  compute.lobanov(),
		   lobanov_gm=compute.lobanov(mean(vdata[,4]), mean(vdata[,5])),
		   lobanov_mf=compute.lobanov.mf(),
		   lobanov_gm_mf=compute.lobanov.mf(f1.all.mean=mean(vdata[,4]), f2.all.mean=mean(vdata[,5])),
		   lobanov_indiv=compute.lobanov.mf(use.indiv=TRUE),
	   	   nearey = compute.nearey(),
		   nearey_adank = compute.nearey(formant.int=TRUE),
		   nearey_adank_indiv = compute.nearey.mf(use.indiv=TRUE, formant.int=TRUE),
		   nearey_adank_wf3 = compute.nearey.wf3(use.indiv=TRUE, formant.int=TRUE),
	   	   nearey_gm=compute.nearey(all.mean=mean(log(c(vdata[,4],vdata[,5])))),
		   nearey_mf=compute.nearey.mf(),
		   nearey_gm_mf=compute.nearey.mf(all.mean=mean(log(c(vdata[,4],vdata[,5])))),
		   nearey_indiv=compute.nearey.mf(use.indiv=TRUE),
		   nearey_wf3=compute.nearey.wf3(use.indiv=TRUE),
		   wattfabricius = compute.wattfabricius()
	   	   )
  left.norms[[j]] <- normed

  normed <- switch(right.method,
		   indiv = organize.data(),
	   	   means = compute.means(),
	  	   barks = compute.bark(),
		   barks_indiv = compute.bark(use.indiv=TRUE),
                   labov = compute.labov(mean(log(c(vdata[,4],vdata[,5])))),
		   labov_g=compute.labov(6.896874),
		   labov_indiv=compute.labov.indiv(mean(log(c(vdata[,4],vdata[,5])))),
		   labov_indiv_g=compute.labov.indiv(6.896874),
		   lobanov =  compute.lobanov(),
		   lobanov_gm=compute.lobanov(mean(vdata[,4]), mean(vdata[,5])),
		   lobanov_mf=compute.lobanov.mf(),
		   lobanov_gm_mf=compute.lobanov.mf(f1.all.mean=mean(vdata[,4]), f2.all.mean=mean(vdata[,5])),
		   lobanov_indiv=compute.lobanov.mf(use.indiv=TRUE),
	   	   nearey = compute.nearey(),
		   nearey_adank = compute.nearey(formant.int=TRUE),
		   nearey_adank_indiv = compute.nearey.mf(use.indiv=TRUE, formant.int=TRUE),
		   nearey_adank_wf3 = compute.nearey.wf3(use.indiv=TRUE, formant.int=TRUE),
	   	   nearey_gm=compute.nearey(mean(log(c(vdata[,4],vdata[,5])))),
		   nearey_mf=compute.nearey.mf(),
		   nearey_gm_mf=compute.nearey.mf(all.mean=mean(log(c(vdata[,4],vdata[,5])))),
		   nearey_indiv=compute.nearey.mf(use.indiv=TRUE),
		   nearey_wf3=compute.nearey.wf3(use.indiv=TRUE),
		   wattfabricius = compute.wattfabricius()
	   	   )   
  right.norms[[j]] <- normed
  j<-j+1
}

if (scale.values) {
   left.norms <- switch(left.method,
			indiv = left.norms,
			means = left.norms,
			barks = left.norms,
			barks_indiv = left.norms,
			labov = left.norms,
			labov_g=left.norms,
			labov_indiv = left.norms,
			labov_indiv_g=left.norms,
			lobanov = scale.results(left.norms),
			lobanov_gm = scale.results(left.norms),
			lobanov_mf = scale.results(left.norms),
			lobanov_gm_mf = scale.results(left.norms),
			lobanov_indiv = scale.results(left.norms),
			nearey = scale.results(left.norms),
		   	nearey_adank = scale.results(left.norms),
		   	nearey_adank_indiv = scale.results(left.norms),
		   	nearey_adank_wf3 = scale.results(left.norms),
			nearey_gm = scale.results(left.norms),
			nearey_mf = scale.results(left.norms),
			nearey_gm_mf = scale.results(left.norms),
			nearey_indiv = scale.results(left.norms),
			nearey_wf3 = scale.results(left.norms),
			wattfabricius = scale.results(left.norms)
			)

   right.norms <- switch(right.method,
			 indiv = right.norms,
			 means = right.norms,
			 barks = right.norms,
  			 barks_indiv = right.norms,
			 labov = right.norms,
			 labov_g=right.norms,
			 labov_indiv = right.norms,
			 labov_indiv_g=right.norms,
			 lobanov = scale.results(right.norms),
			 lobanov_gm = scale.results(right.norms),
			 lobanov_mf = scale.results(right.norms),
			 lobanov_gm_mf = scale.results(right.norms),
			 lobanov_indiv = scale.results(right.norms),
			 nearey = scale.results(right.norms),
		   	 nearey_adank = scale.results(right.norms),
		   	 nearey_adank_indiv = scale.results(right.norms),
		   	 nearey_adank_wf3 = scale.results(right.norms),
			 nearey_gm = scale.results(right.norms),
			 nearey_mf = scale.results(right.norms),
			 nearey_gm_mf = scale.results(right.norms),
			 nearey_indiv = scale.results(right.norms),
		 	 nearey_wf3 = scale.results(right.norms),
			 wattfabricius = scale.results(right.norms)
			 )
}

for (k in 1:length(spkr.names)) {
  write.table(left.norms[[k]], paste(out.prefix,"left_out",k,".txt",sep=""), sep="\t", row.names=FALSE, quote=FALSE)
  write.table(right.norms[[k]], paste(out.prefix,"right_out",k,".txt",sep=""), sep="\t", row.names=FALSE, quote=FALSE)

  if (single.plots) {
    # plot left results
    postscript(file=paste(out.prefix, "left_out",k,".eps",sep=""), family="sans", paper="special", width=8, height=8, horizontal=FALSE)
    switch(left.method,
	   indiv = plot.indiv(list(left.norms[[k]]), spkr.names[k]),
           means = plot.means(list(left.norms[[k]]), spkr.names[k]),
           barks = plot.barks(list(left.norms[[k]]), spkr.names[k]),
	   barks_indiv = plot.barks.indiv(list(left.norms[[k]]), spkr.names[k]),
           labov = plot.labov(list(left.norms[[k]]), spkr.names[k]),
           labov_g = plot.labov(list(left.norms[[k]]), spkr.names[k], "using Telsur grand mean"),
           labov_indiv = plot.labov.indiv(list(left.norms[[k]]), spkr.names[k], "individual vowels"),
           labov_indiv_g = plot.labov.indiv(list(left.norms[[k]]), spkr.names[k], "normalized individual vowels, using Telsur grand mean"),
           lobanov = plot.lobanov(list(left.norms[[k]]), spkr.names[k], "mean normalized values"),
           lobanov_gm=plot.lobanov(list(left.norms[[k]]), spkr.names[k], "mean normalized values, grand means from all speakers"),
	   lobanov_mf=plot.lobanov(list(left.norms[[k]]), spkr.names[k], "normalized mean values"),
	   lobanov_gm_mf=plot.lobanov(list(left.norms[[k]]), spkr.names[k], "normalized mean values, grand means from all speakers"),
	   lobanov_indiv=plot.lobanov.indiv(list(left.norms[[k]]), spkr.names[k], "normalized individual values"),
           nearey = plot.nearey(list(left.norms[[k]]), spkr.names[k], "mean normalized values", 2),
	   nearey_adank = plot.nearey(list(left.norms[[k]]), spkr.names[k], "mean normalized values", 1),
	   nearey_adank_indiv = plot.nearey.indiv(list(left.norms[[k]]), spkr.names[k], "normalized individual vowels", 1),
	   nearey_adank_wf3 = plot.nearey.wf3(list(left.norms[[k]]), spkr.names[k], "normalized individual vowels", 1),
           nearey_gm=plot.nearey(list(left.norms[[k]]), spkr.names[k], "mean normalized values, grand mean from all speakers"),	
           nearey_mf=plot.nearey(list(left.norms[[k]]), spkr.names[k], "normalized mean values", 2),
	   nearey_gm_mf=plot.nearey(list(left.norms[[k]]), spkr.names[k], "normalized mean values, grand mean from all speakers", 2),
	   nearey_indiv=plot.nearey.indiv(list(left.norms[[k]]), spkr.names[k], "normalized individual vowels", 2),
	   nearey_wf3=plot.nearey.wf3(list(left.norms[[k]]), spkr.names[k], "normalized individual vowels", 2),
           wattfabricius = plot.wattfabricius(list(left.norms[[k]]), spkr.names[k])
           )
    dev.off()

    # plot right results
    postscript(file=paste(out.prefix, "right_out",k,".eps",sep=""), family="sans", paper="special", width=8, height=8, horizontal=FALSE)
    switch(right.method,
	   indiv = plot.indiv(list(right.norms[[k]]), spkr.names[k]),
           means = plot.means(list(right.norms[[k]]), spkr.names[k]),
           barks = plot.barks(list(right.norms[[k]]), spkr.names[k]),
	   barks_indiv =plot.barks.indiv(list(right.norms[[k]]), spkr.names[k]),
           labov = plot.labov(list(right.norms[[k]]), spkr.names[k]),
           labov_g = plot.labov(list(right.norms[[k]]), spkr.names[k], "using Telsur grand mean"),
           labov_indiv = plot.labov.indiv(list(right.norms[[k]]), spkr.names[k], "individual vowels"),
           labov_indiv_g = plot.labov.indiv(list(right.norms[[k]]), spkr.names[k], "normalized individual vowels, using Telsur grand mean"),
           lobanov = plot.lobanov(list(right.norms[[k]]), spkr.names[k], "mean normalized values"),
           lobanov_gm=plot.lobanov(list(right.norms[[k]]), spkr.names[k], "mean normalized values, grand mean from all speakers"),
           lobanov_mf=plot.lobanov(list(right.norms[[k]]), spkr.names[k], "normalized mean values"),
	   lobanov_gm_mf=plot.lobanov(list(right.norms[[k]]), spkr.names[k], "normalized mean values, grand means from all speakers"),
	   lobanov_indiv=plot.lobanov.indiv(list(right.norms[[k]]), spkr.names[k], "normalized individual values"),
           nearey = plot.nearey(list(right.norms[[k]]), spkr.names[k], "mean normalized values", 2),
	   nearey_adank = plot.nearey(list(right.norms[[k]]), spkr.names[k], "mean normalized values", 1),
	   nearey_adank_indiv = plot.nearey.indiv(list(right.norms[[k]]), spkr.names[k], "normalized individual vowels", 1),
	   nearey_adank_wf3 = plot.nearey.wf3(list(right.norms[[k]]), spkr.names[k], "normalized individual vowels", 1),	
           nearey_gm=plot.nearey(list(right.norms[[k]]), spkr.names[k], "mean normalized values, grand mean from all speakers", 2),	
           nearey_mf=plot.nearey(list(right.norms[[k]]), spkr.names[k], "normalized mean values", 2),	
           nearey_gm_mf=plot.nearey(list(right.norms[[k]]), spkr.names[k], "normalized mean values, grand mean from all speakers", 2),
	   nearey_indiv=plot.nearey.indiv(list(right.norms[[k]]), spkr.names[k], "normalized individual vowels", 2),
	   nearey_wf3=plot.nearey.wf3(list(right.norms[[k]]), spkr.names[k], "normalized individual vowels", 2),
           wattfabricius = plot.wattfabricius(list(right.norms[[k]]), spkr.names[k])
           )
    dev.off()
  }
}

if (!single.plots) {
  # plot left results
  postscript(file=paste(out.prefix, "left_out_all.eps",sep=""), family="sans", paper="special", width=8, height=8, horizontal=FALSE)
  switch(left.method,
	 indiv = plot.indiv(left.norms, spkr.names),
         means = plot.means(left.norms, spkr.names),
         barks = plot.barks(left.norms, spkr.names),
	 barks_indiv = plot.barks.indiv(left.norms, spkr.names),
         labov = plot.labov(left.norms, spkr.names),
         labov_g =plot.labov(left.norms, spkr.names, "using Telsur grand mean"),
         labov_indiv = plot.labov.indiv(left.norms, spkr.names, "normalized individual vowels"),
         labov_indiv_g =plot.labov.indiv(left.norms, spkr.names, "normalized individual vowels, using Telsur grand mean"),
         lobanov = plot.lobanov(left.norms, spkr.names, "mean normalized values"),
         lobanov_gm=plot.lobanov(left.norms, spkr.names, "mean normalized values, grand means from all speakers"),
         lobanov_mf=plot.lobanov(left.norms, spkr.names, "normalized mean values"),
         lobanov_gm_mf=plot.lobanov(left.norms, spkr.names, "normalized mean values, grand means from all speakers"),
	 lobanov_indiv=plot.lobanov.indiv(left.norms, spkr.names, "normalized individual values"),
         nearey = plot.nearey(left.norms, spkr.names, "mean normalized values", 2),
	 nearey_adank = plot.nearey(left.norms, spkr.names, "mean normalized values", 1),
	 nearey_adank_indiv = plot.nearey.indiv(left.norms, spkr.names, "normalized individual vowels", 1),
	 nearey_adank_wf3 = plot.nearey.wf3(left.norms, spkr.names, "normalized individual vowels", 1),
         nearey_gm=plot.nearey(left.norms, spkr.names, "mean normalized values, grand mean from all speakers", 2),	
         nearey_mf=plot.nearey(left.norms, spkr.names, "normalized mean values", 2),
	 nearey_gm_mf=plot.nearey(left.norms, spkr.names, "normalized mean values, grand mean from all speakers", 2),
	 nearey_indiv=plot.nearey.indiv(left.norms, spkr.names, "normalized individual vowels", 2),
	 nearey_wf3=plot.nearey.wf3(left.norms, spkr.names, "normalized individual vowels", 2),
         wattfabricius = plot.wattfabricius(left.norms, spkr.names)
         )
  dev.off()

  # plot right results
  postscript(file=paste(out.prefix, "right_out_all.eps",sep=""), family="sans", paper="special", width=8, height=8, horizontal=FALSE)
  switch(right.method,
	 indiv = plot.indiv(right.norms, spkr.names),
         means = plot.means(right.norms, spkr.names),
         barks = plot.barks(right.norms, spkr.names),
	 barks_indiv = plot.barks.indiv(right.norms, spkr.names),
         labov = plot.labov(right.norms, spkr.names),
         labov_g =plot.labov(right.norms, spkr.names, "using Telsur grand mean"),
         labov_indiv = plot.labov.indiv(right.norms, spkr.names, "normalized individual vowels"),
         labov_indiv_g =plot.labov.indiv(right.norms, spkr.names, "normalized individual vowels, using Telsur grand mean"),
         lobanov = plot.lobanov(right.norms, spkr.names, "mean normalized values"),
         lobanov_gm=plot.lobanov(right.norms, spkr.names, "mean normalized values, grand mean from all speakers"),
         lobanov_mf=plot.lobanov(right.norms, spkr.names, "normalized mean values"),
         lobanov_gm_mf=plot.lobanov(right.norms, spkr.names, "normalized mean values, grand means from all speakers"),
	 lobanov_indiv=plot.lobanov.indiv(right.norms, spkr.names, "normalized individual values"),
         nearey = plot.nearey(right.norms, spkr.names, "mean normalized values", 2),
	 nearey_adank = plot.nearey(right.norms, spkr.names, "mean normalized values", 1),
	 nearey_adank_indiv = plot.nearey.indiv(right.norms, spkr.names, "normalized individual vowels", 1),
	 nearey_adank_wf3 = plot.nearey.wf3(right.norms, spkr.names, "normalized individual vowels", 1),
         nearey_gm=plot.nearey(right.norms, spkr.names, "mean normalized values, grand mean from all speakers", 2),	
         nearey_mf=plot.nearey(right.norms, spkr.names, "normalized mean values", 2),
	 nearey_gm_mf=plot.nearey(right.norms, spkr.names, "normalized mean values, grand mean from all speakers", 2),
	 nearey_indiv=plot.nearey.indiv(right.norms, spkr.names, "normalized individual vowels", 2),
	 nearey_wf3=plot.nearey.wf3(right.norms, spkr.names, "normalized individual vowels", 2),
         wattfabricius = plot.wattfabricius(right.norms, spkr.names)
         )
  dev.off()
}

# make a txt file describing the results of this script,
#  this makes it easier for php script to know the outcome of this script
write.table(spkr.names, paste(out.prefix,"summary.txt",sep=""), sep="\t", col.names=FALSE, row.names=FALSE, quote=FALSE)